刘厚通 毛敏娟
1) (安徽工业大学数理科学与工程学院, 马鞍山 243002)
2) (浙江省气象科学研究所, 杭州 310017)
如何对地基米散射激光雷达回波信号进行准确定标一直是激光雷达数据反演中的一个重要研究课题[1]. 对于探测高度较高的米散射激光雷达,标定高度一般选取近乎不含气溶胶粒子的清洁大气层所在的高度来确定, 这个高度一般选在对流层顶附近[2−4]. 但对于探测最高高度在6 km左右的米散射激光雷达, 可以选取4–6 km高度范围内X(z ")/(z ")的最小值所在的高度为标定高度[5−7],X(z "),(z ")的物理意义见参考文献[7], 由于4–6 km高度范围内X(z ")/(z ")的最小值所对应的气溶胶消光系数或后向散射系数是变化的, 所以这种定标法得到的标定值存在一定的误差.
利用激光雷达对雾霾探测数据进行探测, 如果雾霾层上无云且激光能够透过雾霾层到达较高的高度, 可以采用“大气清洁层”对激光雷达信号定标[2−4], 或者采用斜率法进行定标[8,9], 利用斜率法对水平探测的激光雷达信号进行定标时, 能够获得比较准确的标定值, 但对于垂直探测的激光雷达数据, 利用斜率法得到的标定值存在一定的误差, 有时误差还比较大; 如果雾霾上无云但激光束穿透雾霾后只能到达6 km左右的高度, 此时仍然可以利用4–6 km的某个高度区域进行定标[5−7], 但是由于4–6 km高度范围内的气溶胶受到地面雾霾颗粒的“污染”, 定标误差较大; 如果雾霾上有低云,利用激光雷达对低云和雾霾同时出现的大气进行探测时, 由于雾霾和低云都具有很强的消光特性,有时激光很难穿透云层, 也就是说很难利用云上的大气清洁层对激光雷达探测信号进行定标, 而处于云层和雾霾之间的气溶胶层, 由于受云层和雾霾的“污染”, 其消光系数变化较大而且很难找到所谓的“均匀层”, 一般也不能利用斜率法进行定标, 这时利用传统的定标法将无法对激光雷达信号进行准确定标, 从而影响雾霾消光系数垂直分布数据的准确获取.
本文提出了一种基于Fernald前向积分方程的气溶胶消光系数反演方法, 为了表达方便, 该法在本文中简称为气溶胶消光系数迭代算法. 利用气溶胶消光系数迭代算法, 无需对激光雷达信号定标就能准确求出气溶胶消光系数, 这种方法不仅适用于激光雷达对低云下雾霾探测数据的消光系数反演, 也可应用于探测高度较低(小于6 km)的米散射激光雷达信号(例如微脉冲激光雷达)的气溶胶消光系数求解.
当米散射激光雷达垂直向大气中发射532 nm波长的激光时, 对每一发激光脉冲, 望远镜接收到高度处大气后向散射回波功率可以用米散射激光雷达方程表示为[10−13]
根据激光脉冲能量和实验测得的激光雷达几何重叠因子, 可以获得归一化激光雷达距离校正信号为
由(3)式可得到大气后向散射系数表示式为
其中的S1为气溶胶激光雷达比, 在本文中S1= 50 sr.从地面到高度处大气透过率可表示为[14]:
利用米散射激光雷达数据反演气溶胶消光系数, 通常采用Klett算法和Fernald算法, 其中根据Fernald提供的前向积分算法对激光雷达信号进行处理得到的气溶胶消光系数廓线表达式为[15,16]:
Fernald后向积分算法的表达式为[15,16]
激光雷达不但能够探测雾霾的强度, 还能够探测雾霾的垂直分布[17−19], 但是对于无法利用“清洁层”进行定标的雾霾探测数据的普适定标方法研究, 至今没见相关文献报道. 图1所示的是2017年3月30日凌晨1点在安徽工业大学东区, 利用安徽工业大学拉曼−米散射激光雷达探测雾霾时得到的532 nm通道激光雷达归一化距离校正回波信号, 该回波信号已经过激光雷达几何重叠因子订正. 从图1中可以看出, 由于雾霾上方存在水云层,激光束不能同时穿透雾霾和水云层, 难以利用“大气清洁层”对激光雷达信号进行定标. 激光雷达在工作时有一段信号盲区, 在信号盲区内几何重叠因子几乎为 0, 难以获得准确的激光雷达回波信号.本文参考美国标准大气模式, 依据文献[20]中的方法得到盲区的激光雷达回波信号数据. 安徽工业大学拉曼−米散射激光雷达由中国科学院安徽光学精密机械研究所研制, 共有四个探测通道: 532 nm偏振平行通道、532 nm偏振垂直通道、607 nm通道及660 nm通道, 其接收望远镜直径为400 mm.532 nm通道的激光脉冲能量为210 mJ, 对晴天气溶胶进行探测时激光雷达回波信号的有效高度一般在15 km以上.
图1 雾霾的激光雷达距离校正回波信号Fig. 1. The range corrected lidar signal about fog and haze.
下面以这组信号的气溶胶消光系数反演为例说明气溶胶消光系数迭代算法的主要步骤.
1)根据激光功率计监测得到的激光脉冲能量,对图1中的激光雷达信号廓线进行归一化处理, 得到归一化信号.
2)第一次迭代. 以图1中A点为迭代起点, 该点选在激光雷达激光发射口处, 以某高度处的B点的消光系数值为迭代结束项, B点的位置应在激光雷达几何重叠因子范围之外, 当然在具体计算中应该固定一个高度, 例如本文中的B点选在1.02 km高度处(因为激光雷达的垂直探测分辨率为0.03 km), 这样有利于计算程序的编写. A点选在激光发射口处, 也即中的0高度处, 这样可以直接利用(4)式得到B点的气溶胶后向散射系数, 进而得到B高度处气溶胶消光系数.
假定A, B两点之间大气透过率为T1(本文中T1= 0.8), 根据(4)式和(5)式算出B点的气溶胶消光系数. 令地面A点的气溶胶消光系数的取值范围为 0–2 km–1,的取值从 0 km–1开始, 以 0.01 km–1的步长递增. 对应每一个A点的气溶胶消光系数值, 根据(7)式和激光雷达归一化信号, 就会得到一条气溶胶消光系数廓线, 同时也得到B点的气溶胶消光系数. 迭代结束条件为
3)第二次迭代. 根据第一次迭代得到的气溶胶消光系数廓线, 利用(6)式再次求出0–1.02 km高度范围内的大气透过率T2, 然后重复第一次迭代的步骤, 再次得到气溶胶消光系数廓线……..
4)最终反演结果确定. 假设根据第k–1次迭代得到的气溶胶消光系数廓线和第k次迭代得到的消光系数廓线在1.02 km高度处的气溶胶消光系数之差小于某一设定值c =0.00001 km–1, 即
在上面的反演中, 由于大气透过率初始值T1是任意设置的, 可能出现因为T1的初设值和真实值相差太大而无法迭代反演出准确的气溶胶消光系数的情况.
对于图1中的激光雷达数据, 当0–1.02 km高度范围内的大气透过率分别预设为50%, 55%,60%, 65%, 70%和75%时, 第一次迭代获得的对应高度范围内的大气透过率分别为57.43%,61.33%, 64.88%, 68.09%, 71.10%和73.69%. 可以看出, 当大气透过率初设值在50%–70%范围内时, 第一次迭代得到的0–1.02 km大气透过率均比初设值高. 初步反演表明: 利用(7)式和3.1节的迭代步骤进行消光系数反演时, 随着迭代次数的增加, 得到的0–1.02 km范围内的大气透过率有趋于真值的趋势, 由此可知0–1.02 km的大气透过率应该高于70%. 而当大气透过率初设值为75%时, 第一次迭代反演得到的0–1.02 km的大气透过率为73.69%, 小于初设值75%, 也就是说0–1.02 km高度范围内的大气透过率小于75%.从上面的分析可以得到0–1.02 km范围的大气透过率在70%–75%之间. 大气透过率初设值和第一次迭代结果对应的大气透过率值之间的关系如图2所示. 图2中的纵坐标表示大气透过率的初设值, 横坐标表示大气透过率初设值的赋值序号.
图2 大气透过率初始值与第一次迭代值之间的关系Fig. 2. The relationship between the initial values and the first iterative values of atmospheric transmittance.
根究上面的分析, 取大气透过率初始值为0.7,根据3.1节的迭代步骤, 经过7次迭代, 反演得到雾霾的消光系数如图3所示, 图3中第7次迭代得到的气溶胶消光系数即为所求. 由于经过3.2节的估算, 0–1.02 km大气透过率的初设值0.7和真实值相差不大, 所以在图3中的7次迭代反演得到的气溶胶消光系数之间差别较小.
图3 利用气溶胶消光系数迭代算法反演得到的气溶胶消光系数廓线Fig. 3. The aerosol extinction coefficient profiles retrieved by the iterative algorithm of aerosol extinction coefficient.
利用气溶胶消光系数迭代算法反演得到的雾霾的消光系数垂直分布是否正确?由于利用现有的定标方法不能对低云下雾霾的激光雷达探测信号进行准确定标, 使得上述反演结果无法得到验证. 但是, 对于探测高度较高的激光雷达信号, 利用气溶胶消光系数迭代算法进行消光系数反演得到的气溶胶消光系数廓线, 和利用大气清洁层定标得到的气溶胶消光系数廓线比较, 如果两者符合得很好, 则证明利用气溶胶消光系数迭代算法能够准确反演得到气溶胶消光系数廓线.
图4是在2005年2月22日晚8点, 利用中国科学院安徽光学精密机械研究所研制的偏振–米散射激光雷达对合肥上空的气溶胶进行探测时, 根据激光脉冲能量及(1)式–(3)式得到的一组归一化激光雷达距离校正回波信号. 偏振–米散射激光雷达探测波长为532 nm, 单发激光脉冲能量可达180 mJ,接收望远镜孔径为254 mm, 设有两种可切换的探测模式(米散射探测模式和米散射–偏振探测模式)能够对气溶胶消光系数、后向散射系数及偏振特性进行测量.
图4 米散射激光雷达距离校正信号Fig. 4. The range corrected lidar signal of the Mie scatter−ing lidar.
利用3.1节中的气溶胶消光系数迭代算法对图4中的激光雷达信号进行处理, 经过4次迭代,得到气溶胶消光系数廓线如图5中曲线A所示.由于图4中的激光雷达回波信号探测高度较高, 可以利用8–14 km的“大气清洁层”进行定标, 并在标定点上、下分别利用(7)式和(8)式进行气溶胶消光系数反演, 反演结果如图5中的B线所示, 由于利用8–14 km高度范围内“大气清洁层”定标法和Fernald方程(7)式和(8)式获得的气溶胶消光系数比较可靠, 本文把图5中的B曲线作为气溶胶消光系数标准廓线.
图5 利用气溶胶消光系数迭代算法与其他定标方法获得的气溶胶消光系数比较Fig. 5. Comparison of aerosol extinction coefficients ob−tained by the iterative algorithm of aerosol extinction coef−ficient and other calibration methods.
从图5中可以看出, 基于Fernald积分方程,利用“大气清洁层”定标法获得的气溶胶消光系数廓线和利用气溶胶消光系数迭代算法获得气溶胶消光系数廓线基本符合, 说明利用气溶胶消光系数迭代算法反演得到的气溶胶消光系数是准确的.
对于雾霾上有云且激光没有穿透云层的激光雷达探测数据, 斜率法是一种常用的定标方法. 利用斜率法对图4中的激光雷达信号进行定标(拟合高度范围为3–4.6 km, 在此范围内的大气符合均匀层的条件), 在标定点上、下分别采用Fernald前向积分方程和Fernald后向积分方程进行反演, 得到的气溶胶消光系数廓线如图5中曲线C所示.从图中可以看出, 利用斜率法定标和气溶胶消光系数迭代算法获得的气溶胶消光系数廓线之间存在一定的差别, 如果以图5中曲线B的反演结果作为标准, 则利用斜率法定标得到的气溶胶消光系数廓线的精度低于利用气溶胶消光系数迭代算法获得的气溶胶消光系数廓线的精度.
需要说明的是, Fernald前向积分公式在地基激光雷达气溶胶消光系数反演中很少用到, 根本原因是该公式具有“误差发散特性”, 容易得到不稳定的解. 但是根据气溶胶消光系数迭代算法, 对(7)式按照本文3.1节的反演方法设置A, B两点, 尽管开始时A, B两点之间的大气透过率初始值和大气透过率真实值之间相差较大, 但是每经过一次迭代, 反演得到的A, B两高度之间的大气透过率反演值和真实值之间的差值就会相应减小, 经过几次迭代, A, B两高度之间的大气透过率反演值和真实值之间的差值已经很小, 这时得到的B高度处的气溶胶消光系数接近B点气溶胶消光系数的真实值. 尽管气溶胶消光系数迭代算法是以Fernald前向积分方程为基础的, 但是只要参照本文3.1节的设置, 利用气溶胶消光系数迭代算法求解A,B之间的大气透过率时能够获得“稳定解”, 进而得到准确的气溶胶消光系数垂直分布.
由于用于迭代的两点之间的大气透过率所在的高度范围包含整个几何重叠因子区, 激光雷达几何重叠因子的误差会直接导致定标高度范围内的大气透过率的误差, 从而对整条气溶胶消光系数廓线的反演精度产生影响. 对地基米散射激光雷达来说, 如果在晴朗、微风的傍晚进行激光雷达几何重叠因子的测量, 并把多次测量的平均值作为激光雷达几何重叠因子的标准值, 可知每次测量的激光雷达几何重叠因子的误差都不大, 其相对误差一般不会超过10%, 且实验测得的激光雷达几何重叠因子的误差一般在500 m以下.
图6是在晴空、微风的傍晚, 利用偏振–米散射激光雷达水平对大气探测时得到的激光雷达几何重叠因子, 该激光雷达几何重叠因子的距离范围为0–0.72 km, 当激光雷达几何重叠因子在0–0.72 km距离范围内的每个值不存在误差、存在5%及10%的误差时(见图6), 对应的归一化激光雷达距离校正回波信号有明显的差别(图7), 由于在几何重叠因子以上激光雷达回波信号相同, 图7只显示0–3 km高度范围内的激光雷达回波信号.
图6 激光雷达几何重叠校正因子的相对误差Fig. 6. The relative errors of lidar overlap function.
图7 不同的激光雷达几何重叠因子误差所对应的激光雷达距离校正信号Fig. 7. The range corrected lidar signals corresponding to different errors of lidar geometric overlap function.
图8 几何重叠因子误差对气溶胶消光系数反演值的影响Fig. 8. The influence of lidar overlap function error on the retrievals of aerosol extinction coefficients.
根据激光雷达几何重叠因子不存在误差、存在5%及10%的误差时得到的激光雷达距离校正信号, 利用气溶胶消光系数迭代算法反演得到的气溶胶消光系数廓线如图8所示, 可以看出, 尽管几何重叠因子存在一定的误差, 利用气溶胶消光系数迭代算法反演得到的气溶胶消光系数廓线和几何重叠因子不存在误差时的气溶胶消光系数廓线几乎重合. 对激光雷达几何重叠因子的多次实际测量表明, 单次测得的几何重叠因子的误差一般不会大于5%, 所以几何重叠因子的误差对气溶胶消光系数迭代算法的反演结果影响不大.
到目前为止, 难以对低空大气探测激光雷达数据和低云下雾霾的激光雷达探测数据进行准确定标. 本文根据Fernald前向积分方程的特点, 提出了一种新的大气探测激光雷达数据反演方法––气溶胶消光系数迭代算法. 这种算法可以先预设某一高度范围内大气透过率一系列的值, 然后根据气溶胶消光系数迭代算法的第一次迭代结果找到大气透过率初始值的准确范围, 进而利用相对精确的大气透过率初始值及气溶胶消光系数迭代算法得到准确的气溶胶消光系数.
初步反演结果表明: 该算法不用对米散射激光雷达数据进行定标就能准确求出气溶胶消光系数,这对不能凭借“大气清洁层”定标的激光雷达数据的气溶胶消光系数反演具有重要的意义. 激光雷达几何重叠因子对气溶胶消光系数迭代算法的反演结果影响不大, 当然几何重叠因子越准确, 利用气溶胶消光系数迭代算法反演得到的气溶胶消光系数误差越小.