2005, 16 (6): 811-819
由于氧气在空气中有固定容积百分比，用卫星测量氧气0.76 μm 吸收带一些通道地面反射太阳辐射的辐亮度，可遥测反演地面气压场及云顶气压[1] 。业务气象卫星将使用高级甚高分辨率的干涉光谱仪，卫星通道将更窄，取得更多通道的辐射信息。
假设干涉光谱仪卫星遥感通道带宽为1 cm-1 。氧气0.76 μm 吸收带分布在12847.1941 cm-1(0.77949 μm)至13165.2576 cm-1(0.75957 μm)范围内。在12930-13220 cm-1的290 cm-1范围内，选190 个波段，大部分波段在氧气吸收带内，作为对比，有11 个波段在氧气吸收带以外。由气象卫星各通道所测辐亮度物理反演气象要素，必须快速而高精度地计算各通道的透过率。逐线计算(LBL)法计算卫星通道的透过率精度较高，但计算量过大，只能作为其他简捷计算方法计算误差的标准，在卫星资料业务物理反演系统中还无法实用。
为了探索适于业务物理反演计算方法，Clough 等[2] 提出用FASCOD 软件计算通道透过率。NESDIS 业务反演系统使用的是McMinllin 等[3] 建立的透过率统计模式，即计算一些标准温度廓线下透过率与气压的关系，实际大气温度廓线与标准大气温度廓线有差别，卫星通道透过率也和标准大气有差别，用统计方法得出通道透过率差别Δτ与温度差别ΔT 的统计关系，可以快捷地计算透过率。刘全华等[4] 改进了McMinllin等人的统计算法，减少了CO2 通道透过率的计算误差。金心等[5] 用k 分布法计算CO2 通道透过率，透过率的最大误差小于0.0035 。
由于洋面或地面气压的相对变化约10-3量级，预研气压的物理反演，要求计算氧气通道透过率的误差小于万分之一。本文计算透过率的速度和精度都满足反演计算的要求。
1 卫星通道透过率的逐线计算在响应函数没有实测数据前，假设响应函数为抛物线，响应函数在通道中心vc 最大，离中心0.5 cm-1处响应函数是中心的一半，在响应函数大于零的范围内，即离中心波数
|
(1) |
在一定高度以上，需要同时考虑Lorentz 与Doppler 加宽效应，即混合加宽。混合加宽的吸收系数公式为∶
|
(2) |
式(2)中，kv是吸收系数，吸收线强度S (T)是温度T 的函数，取参数x=(v-v0)/αD，y=αL /αD，αL 是Lorentz 半宽度，αD 是Doppler 谱线宽度，有几种计算程序可数值计算式(2)中的函数F(x, y)[6~7] 。
由于加宽效应，每条吸收线在每个波数都会按混合加宽的式(2)对吸收系数产生或多或少的贡献，为了精确计算通道Δv内的平均透过率 τ(Δv)，把所有吸收线对吸收系数产生的贡献累加，就是总的吸收系数，这种计算方法称为逐线计算法。通道Δv内平均透过率为:
|
(3) |
式(3)中，U 是辐射传输路程上吸收气体的数量，y (v)是通道的响应函数。
2 逐线计算数值积分步长的选择用式(3)计算通道内的平均透过率，步长Δv过小，计算量过大；步长过大，平均透过率的精度不足。步长是否合适，要看对应波数一个步长的变化，吸收系数k 的变化是否在适当的范围内，即数值积分时，吸收系数k 的变化是否有足够详细的描写。
对流层吸收线主要是Lorentz 加宽，按Lorentz 加宽公式，在v=v0 ±αL 处，k=S/(2παL)，k′(dk /dv)有极大值k′max=S /(2παL2)，乘步长Δv就是一个步长吸收系数的变化的极大值Δk max, 吸收系数k 在一个步长的相对变化为Δkmax/k=Δv/αL，取ε是一个规定范围的数，若Δkmax /k=Δv/αL <ε，即步长Δv<εαL，就算满足吸收系数k 在一个步长的相对变化Δkmax/k 是在适当的范围ε内。
对流层顶以上要考虑Doppler 加宽，在v=v0 ±αD 处，
离v0 越远，吸收系数越小，所以吸收系数k 的相对变化Δk/k 的极大值所在的波数比v0 ±αL 或v0 ±αD 离吸收线中心v0 更远。以中心在13099 cm-1的通道为例，通道范围为13098.2929-13099.7071 cm-1，有两条氧气吸收线，中心在13098.83735 cm-1的吸收线在温度T=296 K 时，强度S=8.476×10-24(cm/mol)，是氧气吸收带内最强的吸收线之一。若通道均分为30000 段，步长=1.4141/30000 ≈4.7×10-5 cm-1，按混合加宽计算30000 段吸收系数的相对变化Δk /k, 不同高度极大值(Δk/k)max计算结果见表 1。
|
表 1 各高度αL，αD ，(Δk/k)max, Δk 及v的数值 |
表 1 中，Δk 与v(cm-1)分别是(Δk/k)max 处一个步长的吸收系数变化及波数。由表 1 可看出:气压很低处如0.001 hPa, αL 小，吸收系数k 的最大相对变化(Δk/k)max=0.0229，小于2.5 %，当气压p ≥10 hPa 时，(Δk /k)max <1 %。当步长Δv变小，(Δk/k)max大致按比例减小。可以认为:把通道范围1.4141 cm-1等分为30000 段，即取步长Δv≈4.7×10-5 cm-1可以保证足够的精度。
3 计算透过率吸收线上、下限的选择为了取得较好的计算精度，在通道Δv=1.4141 cm-1范围内，等间隔取30000 个波数计算透过率，由于逐线计算法要计算几百条吸收线在30000 个波数吸收系数的和，再按式(3)计算平均透过率。
计算一个通道Δv内吸收系数k 的分布，波段外的吸收线对波段内的吸收系数k 也有影响。吸收线强度S 越小、离该通道越远，影响越小，参与计算的吸收线多些，可提高精确度，但增加计算量。取选择标准为ε，如接近通道Δv的若干条吸收线在波段边缘对吸收系数的贡献为k1，离通道Δv较远的吸收线贡献若小于εk1，则此吸收线可忽略不计。以中心在13099 cm-1的通道为例，通道Δv(13098.2929-13099.7071 cm-1)内有两条吸收线，这两条吸收线在通道边缘13098.2929 cm-1及13099.7071 cm-1处，对吸收系数的贡献分别为k1=5.4481×10-25 cm2 及k2=2.162×10-25 cm2，吸收线中心波数小于13098.2929 cm-1，在13098.2929 cm-1处对吸收系数k 贡献小于εk1，吸收线中心波数大813 6 期尹宏:提高“ k 分布法” 计算遥感通道透过率精度的方法于13099.7071 cm-1，在13099.7071 cm-1处对吸收系数k 贡献小于εk2，取ε分别等于10-8，10-7，10-6，应入选的吸收线数目见表 2 。
由表 2 可知，选择ε=10-7比ε=10-6要多选41 条吸收线，计算量的增加不大。所以在逐线计算前按ε=10-7，先计算各波段入选吸收线的上下限。
|
|
表 2 不同ε入选参加计算吸收系数的波数范围cm-1 |
4 Hitran 86 与Hitran 2000 的比较---FASCOD2 的评价
Hit ran database 是美国空军地球物理实验室(AFGL)公布的气体吸收线参数资料，在1982，1986，1991，1996 年和2000 年各出一个版本，过去常用FASCOD2 程序[8] 计算大气透过率，其说明书注明，FASCOD2 1.0 版公布于1988 年7 月7 日，2.0 版公布于1989 年12 月16 日，FASCOD2 OS/2.1 版公布于1989 年6 月10 日。由此推断FASCOD2 所用的是1986 年以前的Hit ran 版本。最新的Hit ran 2000 与Hit ran 86 相比，以中心在13140-13145 cm-1 范围内吸收线(共7 条)为例，T=296 K 时，线强S(cm/mol)有显著差别。
表 3表明，Hi tran 数据从1986 年到2000 年变化10 % 以上，用Hit ran 86 的数据得出的软件FASCOD2 计算氧气在0.76 μm 透过率将产生显著的误差。
|
|
表 3 Hitran 86 与Hitran 2000 线强S 的比较 |
5 吸收系数随气压、温度变化的模拟
张华等① 计算22 个压力和3 个温度下的吸收系数，用指数多项式模拟吸收系数与温度的关系。22 个压力分别是0.01，0.0158，0.0215，0.0464，0.1，0.158，0.215，0.464 ，1，2.15，4.64，10，21.5，46.4，100，220，460，700，1013.25 hPa ；3 个温度分别是200，260 ，320 K 。其表达式为
|
(4) |
① 张华, 石广玉.两种逐线积分辐射模式大气吸收的比较研究.大气科学, 待发表。
由于外推比内插产生较大的误差，大气下层吸收气体密度较大，吸收系数的误差引起计算透过率的误差也大。海平面气压有时大于1013.25 hPa, 张华等所取的最大气压1013.25 hPa 的数值偏小。只用3 个温度得出的式(4)表示吸收系数随温度的变化，难以保证较优的精度。所以本研究取41 个气压、不同气压(高度)按常见温度变化各取8 个温度作为该气压的温度格点。41 个气压数值为:0.049，0.149，0.348，0.749，1.251，1.751 ，2.503，3.502，4.502，6.005，8.505，12.51，17.50，22.50，27.50，40.03，55.0，77.5，92.5 ，107.5，125.0，142.5，175.0，224.9，274.9，324.9，374.9，415.0，452.4，487.5，534.9 ，594.9，644.9，685.0，739.9，814.9，884.9，935.0，1080 hPa 。各气压取不同的8 个温度，第一个温度T1 比此高度最低温度还低，第八个温度T8 比此高度最高温度还高，如地面的最高和最低温度分别为210 K 及330 K，对流层顶(226 hPa)的最高和最低温度分别为196 K及260 K，平流层顶(1 hPa)的最高和最低温度分别为211 K 及319 K 。在T1，T8之间大致等间隔取6 个温度T2，…，T7 。按上述气压、温度格点事先计算41×8 组吸收系数，使用时按温度、气压双重样条内插，以提高计算不同温度气压条件下吸收系数的精度。
6 吸收线混合加宽计算的简化在主频为1.8 GHz 的Pentium 4 计算机上，计算Lorentz 公式需1.5×10-7s, 计算混合加宽公式(2)需4.4×10-7s, 计算时间为计算Lo rent z 公式的2.93 倍。如计算精度要求不高，当式(2)中αL 比αD 大很多，即y >3.2 时，可以忽略Doppler 加宽效应，用较简单的Lorentz 公式计算吸收系数。
由于用氧气吸收波段遥测气压分布要求通道平均透过率精度在万分之一以内，为考察Lorentz 公式计算的精度，对不同x, y 条件下按Lorentz 公式及混合加宽公式(2)计算两者吸收系数比kL/km, 结果见表 4 。误差大于万分之一的下加横线。
|
|
表 4 不同x, y 条件下，kL/km 的数值 |
由表 4 可知，若要求计算误差小于万分之一，用简化的Lorentz 公式代替计算量较大的混合加宽公式(2)，条件y >3.2 是不够的，需提高到y >80 或x >200，才能保证透过率精度在万分之一以内。目前计算机的速度加快，使大量使用混合加宽公式(2)计算吸收系数成为可行。
由于波数离吸收线中心越远，吸收线对吸收系数的贡献越小，所以一般计算中采取在x 超过一定数值(如1000)时，就忽略此吸收线对吸收系数的影响，以节省计算时间，称为线翼截断。目前计算机的速度加快，Lorentz 公式的每次计算时间为1.5×10-7 s, 不进行或减少线翼截断可提高计算精度，增加的计算量是计算机能承担的。
7 Gauss 数值积分的改进为了较快地计算通道平均透过率，1981 年石广玉[9] 提出，将Gauss 数值积分法用于k分布法，取N 为截点数，为了把0 和1 作为积分上下限，石广玉采取式(5)的变换后，t 的上下限为0 和1 。Δt 代表各分段的权重，其和为1 。
|
(5) |
|
(6) |
|
(7) |
在41×8 个压强和温度组合下，等间隔计算30000 个波数的吸收系数kv值。以响应函数为权重，按大小在通道Δv内从大到小重新排列，使kv成为单调下降的函数，求N=6 及N=14 的Gauss 截点kv值，得出各通道的吸收系数k 分布数组H(N，41，8)，N 为截点数，41 是压强个数，8 是温度个数。
在任意压、温廓线下，已知气压、温度，可由数组H(N，41，8)用双重样条内插计算N个截点的吸收系数k i 值。式(7)中，i 是k 分布各子分段标号，j 是大气分层高度标号，ΔUj 是某层单位截面氧分子数。用式(7)可计算通道平均透过率。在标准大气温度廓线下，将k 分布法与逐线计算法得出各高度的平均透过率相比，求出均方差。
为了缩小k 分布法的误差，变化式(6)的指数r, 不同的指数r, k 分布法的均方差不同。以13022 cm-1为中心的通道为例，将式(6)改为
|
(8) |
不同的r 和N 时的均方差见表 5 。由表 5 知，选用N=14，r=4.3，均方差为6.26×10-6，误差数值较小。对含有S >10-26 (cm/mol)强线的通道，用式(8)变换，以12990 cm-1通道为例，表 6 取N=6 和N=14 均方差都较大。
|
|
表 5 13022 cm-1通道不同变换指数r 及截点数N 时的均方差(×10-8) |
|
|
表 6 12990 cm-1通道不同变换指数r 及截点数N 时的均方差(×10-8) |
石广玉用k 分布法处理的波段较宽(约100 cm-1)，波段内吸收线多达几百条，计算的步长也较大，吸收系数重排的结果，吸收系数最大的区域吸收系数随波数变化最快(图 1 中A 线)。为了使Gauss 数值积分的分段在变化较快处分段更细一些，所以采用式(6)表示的变换。
氧气吸收带的通道波段较窄(1 cm-1)，波段内吸收线较少(多数0-2 条，最多20条)，计算的步长较小，吸收系数按大小重排后，吸收系数的相对变化dk/(k dv)在离吸收线中心Lorentz 或Doppler 宽度以外有极大值，即吸收系数k 不是在k 最大的区域变化最快，而是在离吸收线中心Lorentz 或Doppler 宽度以外变化最快，图 1 中B线是中心波数为13033 cm-1的通道，气压为375 hPa, 温度为235 K，通道均分30000 段，计算所得吸收系数按大小重排的分布实况，B 线左方吸收系数极大处，代表吸收线中心附近，曲线斜率并不是最大。
|
|
|
图 1. 吸收系数按大小重排的图像
A 线为石广玉计算值[9]，B 线为本文计算值 |
|
为了减小k 分布法的误差，在吸收系数相对变化dk/(kdv)极大附近把k 分布分为两段，分别用N=6 和N=14 的式(7)进行变换，使dk/(kdv)极大附近，Gauss 数值积分的分段更细一些，分20 段的k 分布计算平均透过率均方差更小，以12990 cm-1 通道为例，表 6 取N=6 和N=14 均方差都较大，改为三参数，N=20，均方差显著缩小。
取三参数，Q 是百分数分界点，大致在吸收系数相对变化dk/(k d v)极大值附近，r2 是吸收系数较大一侧N=6 的指数，r1 是吸收系数较小一侧N=14 的指数，反复调整Q，r2 及r1 的数值，由表 7 可见，Q=0.945 ，r1=1.92，r2=1.5 这一组，均方差=1.34×10-6最小(小于万分之一)。有的通道内吸收线较多，吸收系数随波数变化较复杂，需把通道分更多段，分段用k 分布模拟，截点数>20，才能使平均透过率的均方差小于万分之一。如13163 cm-1 通道内有16 条吸收线，取N=20 时，均方差为0.00093257，仍大于0.0009，把通道分10 段，分段用k 分布法模拟，共取114 个截点，均方差为0.0000148，缩小63 倍，小于0.0001 。
|
|
表 7 12990 cm-1通道N=20 时不同分界点Q 及变换指数r1 与r2 的均方差(×10-8) |
8 结论
模拟计算了190 个通道，在主频为1.8 GHz 的Pentium4 计算机上，用改进的k 分布法计算41 个高度到大气顶的通道平均透过率:截点数N=6 的有52 个通道，计算机用时0.0095 s；截点数N=14 的有53 个通道，用时0.0175 s ；截点数N=20 的有59 个通道，用时0.023 s ；N=6，14，20 的通道共165 个，占190 个通道的86.8 %，其余25 个通道，截点数N 分别是28，34，42，48，54，60，62，82，88，108，114，120 及136，其中N=48，用时0.048 s ；N=136，用时0.128 s 。与逐线计算法相比，通道透过率的均方差都小于3×10-5 。逐线计算法计算通道透过率所需计算机时间与通道内吸收线的数量有关，所需时间为79 至162 s。k 分布法比逐线计算法快得多。在逐线计算基础上设计的改进的k 分布法，计算透过率的速度和精度都满足反演计算的要求。
| [1] | Barton I J, Scott J C, Remote measurement of surface pressure using the oxygen A-band of absorption. Appl Optics, 1986, 25, (19): 3502–3507. DOI:10.1364/AO.25.003502 |
| [2] | Clough S A, Kneizy F X, Rothman L S, Atmospheric Spectral Transmittance and Radiance:FASCOD1B. Proceedings of SPIE, 1981, 277: 152–166. DOI:10.1117/12.931914 |
| [3] | McMinllin L M, Fleming H E, Atmospheric transimittance of an absorbing gas:A computationally fast and accurate transmittance model for absorbing gases with constant mixing ratio in inhomogeneouce atmospheres. Appl Opt, 1976, 15, (2): 358–363. DOI:10.1364/AO.15.000358 |
| [4] | 刘全华, 董超华, 黎光清. CO2大气透过率的统计算法. 大气科学, 1989, 13, (2): 228–237. |
| [5] | 金心, 尹宏. 计算HIRS/2通道透过率的指数和模式. 应用气象学报, 1996, 7, (1): 96–102. |
| [6] | Diayson S R, Rapid computation of Voigt profile. J Quan Spec Rad Tran, 1976, 16, (7): 611–614. DOI:10.1016/0022-4073(76)90029-7 |
| [7] | Schreier F, The Voigt and complex error function:a comparison of computional methods. J Quan Spec Rad Tran, 1992, 48, (5-6): 743–762. DOI:10.1016/0022-4073(92)90139-U |
| [8] | Ridgway W L,Moose R A,Cogley A C.Atmospheric Transmittance/Radiance:Computer Code FASCOD2.1988,AFGL-TR-82-0392。 |
| [9] | 石广玉. 大气辐射计算的吸收系数分布模式. 大气科学, 1998, 22, (4): 659–676. |