6S大气传输模型修改源码添加、自定义卫星光谱响应（以HJ-1B CCD为例）

6S大气传输模型修改源码添加、自定义卫星光谱响应（以HJ-1B CCD为例）

最近要做国产卫星的大气校正，打算用6s模型模拟气溶胶的查找表，但是发现6s模型中没有国产卫星的相应光谱响应函数，只能在输入的时候把光谱响应函数整个输进去，觉得麻烦，研究了一下，发现修改源码可以实现添加自定义的卫星光谱响应函数，现在与大家分享一下我的做法：

（p.s. 我是在Linux系统下做的，Windows系统大家自己摸索吧··）

所需文件

6s源码：http://6s.ltdri.org/pages/downloads.html

目前最新的是2015年发布的V2.1版本。

所需卫星的光谱响应函数文件

我用的是HJ-1B卫星的CCD相机的光谱响应函数，可以在

http://www.cresda.com/site1/Downloads/gpxyhs/index.shtml

下载，这个网页还有GF、ZY等卫星的光谱响应函数文件。（吐槽，在网上找GOCI卫星的光谱响应函数怎么也找不到，不知道棒子藏着这些干嘛）。

步骤

1. 修改main.f

解压压缩包之后，找到main.f打开。（建议大家浏览一下这个文件，各种输入参数的说明、注释非常到位，还有彩蛋）

大家可以发现首先出现与传感器有关的代码是在第364行左右：



nsat 这个东西（对Fortran不了解不知道怎么叫它····姑且叫它东西吧···）保存了传感器波段名称字符串，在输出结果的时候用的上。所以我们要在nsat的末尾加上我们自己的传感器：



我加了4个字符串，HJ1/CCD 1-4，代表CCD的4个波段。之所以字符串长度和之前的字符串长度一致（包含空格17个字符），除了看着爽之外，也是有原因的，看到代码292行左右：



这里是nsat的定义，可以看到是一个204*17的character，就是说共有204个XX，每个XX长度是17字符。这个204也是我后来修改的，之前是200，加了4。

接下来，看到代码大约1100行，关于spectral conditions的说明：



不同输入对应不同的spectral condition。当输入大于2时，便是内置的传感器波段了。我们也装模作样的把我们自己的传感器写到这个注释里面：



当然这个编号（200-204）也是有用的。再往下看，大约1370行左右：



这个注释说明了不同传感器输入对应应该goto哪一个语句（尼玛Fortran这么多的goto也是看醉了），我们也装模作样的把自己的传感器加进去。如果输入是200-203，则goto 165编号的行。注意检查一下goto的编号是否已经存在了。

这只是修改了注释，接下来修改代码。就在注释下方编号18的goto，观察一下就是与注释相对应的，我们把自己的传感器加在下方：



然后就是写编号为165的代码行了，在编号164的下方，稍微修改一下函数名就好：



至此，对mian.f的修改完成！

2. 制作传感器文件（HJ1.f）

下面是重头戏，制作我们自己的传感器文件。新建一个文件，我取名为HJ1.f，注意文件名和上一步call的函数名一致。

首先可以打开一个6s自带的传感器文件，就拿MODIS吧，打开看一下，先全部复制粘贴过来吧·····没关系，后面慢慢改。

由于我们只有4个波段，所以定义sr只需要4行（下图标号11的行），我还装模作样的写了一些注释：



下面就是定义每一个波段上的光谱响应了。我先在matlab里面把光谱响应函数的光谱间隔差值到了2.5nm，这也是必须的。下面是CCD第一个波段的例子：



一行一行解释，首先注释行说明了起始波长和终止波长，下面是具体每一个波长上光谱响应值。注意每一行是1501列的，这个不能修改，不然会出错。因为6s中光谱范围是0.25-4um，间隔为0.0025um，所以有(4-0.25)/0.0025+1=1501个波长。14行的数据中第一项69*0说明这一行数据前69个数为0，69=(0.4225-0.25)/0.0025；最后一项1385 * 0说明最后1385个数为0，1385 = (4-0.5375)/0.0025。检验一下，前面有69个0，后面有1385个0，中间有47个光谱响应值，一共是69+1385+47=1501个波段，简直颇费。只要保证起始波长（0.4225）正好对应第一个非0光谱响应值（0.092），而终止波长（0.5735）正好对应最后一个非0光谱响应（0.175），就能按照上面的计算方法保证不多不少正好。不然的话就要绕一些了···总之大家检查一下是不是一共有1501个数据，还有波段上是否每个波长都有光谱响应值就ok。

按照上面的方法继续完成剩余的3个波段，最后需要修改各个波段的上下限：



大功告成！

3. 修改makefile

最后要修改一下makefile，很简单，打开makefile，在那一长串的字符最后加上我们的 HJ1.o就OK了：



4. 编译

下面就是编译了，cd到6s的文件夹下，输入命令 make ，回车，

and good luck！

如果有问题，可以根据错误提示找找哪里出错了，我之前有几次就是在main.f中不小心多打了一些符号或者错删了一些语句造成了失败。可以与原来的main.f比较一下很快就发现问题了。

5. 运行6s

如果编译成功的话，就可以运行6s程序了，

cd到6s的文件夹，输入命令 ./sixsV2.1 回车

此时会等待输入各种参数。



关于6s的参数说明，还是推荐大家看main.f里面的说明，这里贴一张6s User Manual里面的图：



下面是我的输入，注意看第9行，数字201代表了我们自己加入的CCD band2，





然后回车···good luck······



output中可以看到我们自己加入的传感器波段名



以及大气校正的结果。