逐像元大气校正，常预先计算查找表（LUT，LookUp Tabel），6S大气辐射传输模式也可以用来计算LUT。但6S源程序输出信息多，且浮点数输出精度低，不利于提取关键信息生成LUT，本文描述了怎样修改6S源码以生成LUT。
首先确定LUT内容要素。
阅读MODIS M?D04 气溶胶产品生成算法文档（NASA相关网页），归纳了以下查找表要素：
- 1) 太阳天顶角观测天顶角太阳方位角观测方位角之差（相对方位角）散射角
- 2) 大气模式
- 3) 气溶胶模式
- 4) 550nm气溶胶光学厚度
- 5) 波段号
- 6) 大气透过率
- 7) atm. intrin. ref.（个人理解，这是大气程辐射换算后的反射率，用于校正计算）
- 8) total sca. （总散射，包括rayleigh散射和气溶胶散射，用于校正计算）
- 9) 表观反射率
- 10) 校正后的地表反射率
- 11) xa xb xc参数（物理意义不明，用于校正计算）
其次，阅读源码，明确LUT各要素在6S源码中的变量名。
6S大气校正计算源码（Excel验证中采用此公式）
rog=rapp/tgasm
rog=(rog-ainr(1,1)/tgasm)/sutott/sdtott
rog=rog/(1.+rog*sast)
xa=pi*sb/xmus/seb/tgasm/sutott/sdtott
xb=srotot/sutott/sdtott/tgasm
xc=sast
由计算源码确定需输出的变量。下面是输出LUT要素的Fortran77代码示例，建议放在源码main.f中stop一句之前。
write(*,*) asol,phi0,avis,phiv,adif,phi,idatm,iaer,v,taer55,
1 iwave,tgasm,ainr(1,1),sutott*sdtott,rapp,rog,xa,xb,xc
- 其中，asol是太阳天顶角，
- phi0是太阳方位角，
- avis是观测天顶角，
- phiv是观测方位角，
- adif是散射角，
- phi是相对方位角，
- idatm是大气模式号，
- iaer是气溶胶模式号，
- v是水平能见度，
- taer55是550nm气溶胶光学厚度，
- iwave表示波段号，
- tgasm表示大气透过率，
- ainr（1,1）是大气本身的反射率（姑且这么理解），
- sutott*sdtott表示总散射，
- rapp是表观反射率，
- rog是校正后的地表反射率，
- xa，xb，xc是校正计算参数。
Fortran77每行长度不能超过80个字符，续前行需在特定位置指明(示例中的iwave前的1即表示续行）。
该示例源码没有指定任何输出样式。浮点数会按默认样式输出，小数点后的位数比较多，精度较好。
挑选一个查找表文件用Excel验证后表明，excel公式计算值与6S输出值之间最大误差小于1E-7。说明方法是可行的。
再次，编译源码，编写Shell脚本：
编译环境：OpenSUSE操作系统 g95编译器，版本未知。
编译命令：g95 *.f -o sixs（在BRDF相关代码处可能有几个warning，本文不涉及BRDF，故暂不修改调试。在Windows下用f77编译，无warning，编译通过）
生成LUT的bash脚本getLUT.sh：
- 最后调用该脚本
1: function LUTCalc(){
2: #{42,44,48,25,27,30}
3: IBand=$1
4: echo "Band ${IBand} is running"
5: for Iref in {0.1,0.5}
6: do
7: fstat=${IBand}"_"${Iref}.res
8: echo "asol,phi0,avis,phiv,adif,phi,idatm,iaer,v,taer55, iwave,tgasm,ainr,tott,rapp,rog,xa,xb,xc" >> ${fstat}
9: for SolarZ in {0,15,30,45,75}
10: do
11: for SolarAz in {0,45,90,135}
12: do
13: for ViewZ in {0,15,30,45,75}
14: do
15: for Iaer in {1,2,3,5,6,7}
16: do
17: for Idatm in {1,2,3,4,5,6}
18: do
19: for IAod in {0.1,0.2,0.5,1.0};
20: do
21: # Run sixs and output
22: ./sixs <<EOF | tail -n 1 >> ${fstat}
23: 0
24: $SolarZ $SolarAz $ViewZ 0 3 15
25: $Idatm
26: $Iaer
27: 0
28: $IAod
29: -0
30: -1000
31: ${IBand}
32: 0
33: 0
34: 0
35: 0.0
36: -$Iref
37: EOF
38: done
39: done
40: done
41: done
42: done
43: done
44: done
45: }
46: function getLUT()
47: {
48: echo "LUT is Calcing"
49: for iwave in {42,44,48,25,27,30}
50: {
51: LUTCalc ${iwave}
52: } &
53: }
.csharpcode, .csharpcode pre
{
font-size: small;
color: black;
font-family: consolas, "Courier New", courier, monospace;
background-color: #ffffff;
/*white-space: pre;*/
}
.csharpcode pre { margin: 0em; }
.csharpcode .rem { color: #008000; }
.csharpcode .kwrd { color: #0000ff; }
.csharpcode .str { color: #006080; }
.csharpcode .op { color: #0000c0; }
.csharpcode .preproc { color: #cc6633; }
.csharpcode .asp { background-color: #ffff00; }
.csharpcode .html { color: #800000; }
.csharpcode .attr { color: #ff0000; }
.csharpcode .alt
{
background-color: #f4f4f4;
width: 100%;
margin: 0em;
}
.csharpcode .lnum { color: #606060; }
>source getLUT.sh
>getLUT
最好晚上计算早晨看结果，如果CPU给力的话，几个小时后就可以得到结果。
下面是生成的LUT示例:
- asol,phi0,avis,phiv,adif,phi,idatm,iaer,v,taer55, iwave,tgasm,ainr,tott,rapp,rog,xa,xb,xc
- 0.0000000 0.0000000 0.0000000 0.0000000 180.00000 0.0000000 1 1 63.664398 0.10000000 25 0.98984975 6.89081103E-02 0.80637234 0.10000000 3.87301818E-02 1.98730151E-03 8.71319696E-02 0.14777245
- 0.0000000 0.0000000 0.0000000 0.0000000 180.00000 0.0000000 1 1 26.739149 0.20000000 25 0.98984975 7.75793791E-02 0.76532620 0.10000000 2.94530466E-02 2.09388486E-03 0.10336593 0.16389242
- 0.0000000 0.0000000 0.0000000 0.0000000 180.00000 0.0000000 1 1 8.4940033 0.50000000 25 0.98984975 0.10173188 0.64870018 0.10000000 -2.69861287E-03 2.47033220E-03 0.15994170 0.20088956
- 0.0000000 0.0000000 0.0000000 0.0000000 180.00000 0.0000000 1 1 3.5674956 1.0000000 25 0.98984975 0.13688390 0.48083964 0.10000000 -7.89646432E-02 3.33272247E-03 0.29038188 0.24035060
- ……