Progress and prospect of atmospheric scattering correction for laser altimetry satellite

卫星激光测高作为星载激光雷达的一个子类，因其能获得亚米级甚至厘米级高程测量精度，已成为一种新型对地观测技术并受到广泛的关注^[1-2]。2003、2018年美国国家航空航天局（National Aeronautics and Space Administration，NASA）先后发射了冰、云和陆地高程卫星ICESat （Ice、Cloud and land Elevation Satellite）和ICESat-2，在极地冰盖厚度变化、全球树高和生物量反演、全球高程控制点获取、湖泊水位监测等方面取得了很多科研应用成果^[3-4]，其中ICESat-2搭载了“先进地形激光测高系统”（Advanced Topographic Laser Altimeter System，ATLAS），采用了新型光子体制，波束数量和激光点密度均比GLAS有较大提高^[5-6]。NASA还在2018年12月发射了“全球生态系统动力学调查”(Global Ecosystem Dynamics Investigation Lidar，GEDI)^[7]，成功将其布置在国际空间站（International Space Station，ISS）上，主要用于获取热带和温带区域地表高程、森林垂直结构信息。我国先后在2016、2019年发射的资源三号02星、高分七号卫星上配备了激光测高仪^[8-9]，主要用于高程控制点的获取^[10-11]，并计划在2021年发射的陆地生态系统碳监测卫星上搭载多波束的线性体制全波形激光测高仪。

不同于LOLA (Lunar Observer Laser Altimeter)、MOLA (Mars Orbiter Laser Altimeter)对月球及火星等地外星球激光测高探测^[12-13]，对地观测卫星激光测高除受复杂地形影响外，大气对激光测距及最终的测高精度影响不容忽视^[14-15]。目前，对于大气折射延迟分析及其改正已经开展了很多研究，相关算法基本成熟且能工程化应用^[16-17]。但大气散射的影响及改正目前还处于探索阶段，离业务化应用还有较大差距。在大气传输过程中，云、气溶胶等粒子的吸收、散射作用，除导致激光束质心发生一定变形外，还会产生测距偏差^[18]，而目前大气散射改正没有实现业务化，导致因大气散射而发生质量退化的激光点，要么被直接抛弃降低了总的可利用率，要么被应用但降低了精度的可靠性。因此，开展大气散射改正对于有效保证激光点的可靠性、提高总的可利用率具有重要价值。

文中对国内外激光测高载荷的发展现状及技术特点进行了描述，梳理了这些卫星所采用的大气探测算法、大气散射延迟改正方案，在此基础上进行对比分析和总结归纳，为我国后续激光测高卫星大气散射改正技术的发展提出了建议。

在大气传输过程中，卫星对地发射的激光脉冲受到云、气溶胶等粒子散射影响产生测距偏差。对大气散射影响的探索是一个漫长的过程，Platt等人^[19]在1981年首次强调了卷云中多次散射影响，在此基础上，Eloranta^[20]发展了计算激光雷达回波信号中多次散射影响的实际模型。Bissonette^[21]利用辐射传输模型和实验室测量数据，模拟了多次散射的影响。Bruscaglioni等人^[22-23]对多次散射进行了详细的实验论证，表明了大气散射在激光雷达回波信号中的重要作用。Mahesh等人^[24-25]利用蒙特卡洛模型模拟计算了激光测高过程中云对光子多次散射的影响，结果表明，对于薄云和低层云，大气多次散射导致的路径延迟大于20 cm，但是仅涉及仿真结果的模拟及其分析。Brenner等人^[26]通过对比雷达测高系统Envisat和地学激光测高系统GLAS的测量数据发现，大气的多次散射效应是影响激光测距精度的重要因素。渠丽新等人^[27]模拟分析了云对激光脉冲展宽效应的影响，实验结果表明，与云层厚度、单次散射比等因素有关。宋雪平等人^[28]对云雾多次散射条件下所接收的激光数据进行了评估，发现多次散射对系统稳定性有较大的干扰。李颖颖等人^[29-30]通过蒙特卡洛模拟反演多次散射条件下激光雷达消光系数等大气参数，结果表明，大气参数反演必须要考虑散射的影响。陈舒杭等人^[31]基于半解析蒙特卡洛模拟反演特定条件下云多次散射对激光测高数据精度的影响，结果表明，测高偏差与云层高度、云有效粒子半径、云光学厚度、浓度等多种因素有关，必须基于实时的大气观测参数和半解析蒙特卡洛法修正测高偏差，以提高测距值，但实际观测中很难同步获取实时大气相关数据。周辉等人^[32-33]通过蒙特卡洛模拟激光在云中单次散射情况下的回波波形，发现波形会发生拖尾变形而影响测距精度。王婷等人^[34]基于蒙特卡洛模拟分析多个红外波长激光在雾中的衰减效应，结果表明，1 064 nm激光在平流雾中具有很好的传输性，雾对激光的衰减效应与消光系数、散射系数等多个因素密切相关。么嘉棋等人^[35-36]对ICESat的大气范围波形数据进行分析并反演了云光学厚度等参数，为后续开展大气散射延迟改正提供了一定参考。Duda、Palm等人^[37]基于GLAS参数的蒙特卡洛模拟实验表明，云光学厚度（Cloud Optical Depth，COD）为0.2的云会引起20 cm的测高偏差。Yang等人^[38-39]基于ATLAS参数的三维蒙特卡洛模拟实验表明，COD为1的云会引起4~6 cm的测高偏差。

从上述研究可以看出，大气散射延迟影响对于分米级甚至厘米级的高精度激光测量而言，是不容忽视的、必须加以考虑并修正。进入21世纪，国内外先后发射了多种型号的对地观测激光测高载荷，表1列出了相关的技术指标。国产的资源三号02星和高分七号上的激光载荷由于设计及数传等原因，没有获取大气段回波数据，因此不具备大气探测的可能性。CALIOP、CATS获取了大气传输范围内数据，由于主要探测对象并不是地表高程，并没有进行大气散射校正，但采用了各种大气参数进行数据质量控制。GLAS采用动态阈值法实现大气探测，并基于蒙特卡洛模拟结果实现了大气散射误差校正。GEDI通过多参数约束方案，将受大气散射影响的数据准确识别，实现数据质量控制。对于光子计数体制激光雷达而言，多波束、高重频是非常典型的优势，但也为大气散射改正提出新的挑战：如何在兼具技术特点、大气探测范围等因素的基础上实现误差改正^[40]。对此，CATS、ATLAS采用了不同的解决方案。CATS基于阈值法实现大气探测，并基于多参数约束算法获取云相态、气溶胶类型、光学厚度等大气参数。由于高激光重频，数据本身无法实现大气探测，ATLAS采用第三方地表反射率数据集实现了大气散射探测，但是反演结果不包含大气垂直分布、层高度数据，对大气散射改正造成很大的困扰。

目前的对地观测激光载荷在大气散射探测方面，主要有基于后向散射强度的动态阈值探测、多参数约束的直接大气探测以及基于表观反射率的间接大气探测等3种方法。

由于后向散射效应的叠加性，云或气溶胶区域内的整体后向散射系数远远大于天气晴朗区域。基于这种特性，可以通过设置一个合理的阈值进行云或气溶胶探测。以GLAS为例，GLAS大气通道以76.8 m的垂直分辨率记录了从地面到海拔约41 km的后向散射廓线数据。基于后向散射系数经验阈值判断层结构为云或气溶胶，之后定量反演光学厚度，最后通过蒙特卡洛模拟实现大气散射误差校正。GLAS采用不同的阈值算法将数据分别保存成4种产品（40、5、1、0.25 Hz），此节主要介绍了5 Hz产品所采用的大气探测算法，算法流程如图1所示。

Figure 1.  Cloud & aerosol search algorithm flow

如果3个连续的后向散射信号大于阈值，则认为在该高度存在云、气溶胶。计算激光脉冲在当前高度的后向散射系数（β_m）和该海拔高度的散射常数（Ψ）之和，并乘以与当前高度有关的因子（F），通过下式可以计算出层高（Z_t）：

式中：Ψ=1.0×10⁻⁶；F=(0.70+(i_bin /2)/180.0)，i_bin指按大气探测范围每隔76.8 m所划分的节点编号。斜率法通过计算层高度范围内有效消光系数的积分，得到激光光斑范围内柱状界面的光学厚度。斜率法计算流程如下：

斜率法反演大气光学厚度步骤如下：

式中：z为激光脉冲所在海拔高度；S_α(Z)为大气对激光的削弱系数（或称有效消光系数），表示在Z高度散射的全部能量）；β_α(Z)为大气后向散射系数；α_α(Z)为大气消光系数。

式中：τ为光学厚度；Z_top为层顶高度；Z_bottom为穿透层时的高度。至此，实现了基于回波波形数据的光学厚度参数反演，反演流程见图2。

Figure 2.  Optical depth inversion process

GEDI以缓冲区的形式保存大气传输范围内的回波数据，以区别于地物回波并保存在L2产品数据集中。通过结合波形灵敏度、回波能量、土地覆盖类型等多个参数，实现大气探测。

由于大气条件从完全晴朗到多云，GEDI返回信号强度会有很大的变化，穿透遮挡的能力取决于返回信号的强度。当激光穿透云层、比较密集的树冠时，信号强度会迅速下降，甚至衰减至噪声水平。为了识别没有检测到真实地面的波形，GEDI提出了波形灵敏度（Waveform Sensitivity）指标，由公式（4）得到。

式中：WS表示波形敏感度；min_detection_energy表示回波高斯拟合最小积分面积；rx_energy表示激光脉冲的累积能量。为了保证使用用户可以轻松的判断错误数据、较低质量数据，GEDI使用几个单独的质量评价标准组合得到总体质量控制标准，判断数据质量良好（quality_flag=1）的流程如图3所示^[41]。

Figure 3.  Atmospheric detection algorithm of GEDI

图中，Sensitivity_a<n>表示波形敏感度。rx_processing_a<n>是指有效波形数据，其中<n>表示用于检测信号的不同参数。min_detection_energy是指通过有效波形计算的高斯分布下可检测最小积分面积，是Sensitivity_a<n>下的子类别。rx_assess是指接收激光脉冲的特征参数（例如，最小和最大幅度、能量），保存在L2A中。rx_energy是指接收激光脉冲的累积能量，是rx_assess下的子类别。在满足上述条件的基础上，根据激光点所在土地类型，进行进一步判断。如果满足上述条件，那么该数据质量判断为良好。图4是GEDI的一段真实数据，黑色为原始数据中经质量控制后筛掉的点，红色为经筛选后质量良好的点（quality_flag=1），可以看出云散射数据、背景噪声被有效过滤掉。

Figure 4.  Data quality control diagram of GEDI

线性体制激光雷达可以通过接收大气后向廓线信号的变化判断是否受到云散射影响，但是对于ICESat-2/ATLAS光子计数体制而言，激光重频达到10 KHz，会出现激光脉冲混叠现象，无法采用常规的大气探测算法^[42]。ATLAS采用基于表观反射率的大气探测算法^[25]，当激光受到大气散射影响时，表观地表反射率数据与清洁大气时有明显的区别，通过设置合适的阈值可以实现大气探测，进而分析散射影响。

基于表观反射率ASR （Apparent Surface Reflectance，ASR）的大气探测算法流程如下：从25 Hz观测数据中分割出对应激光点每四秒数据的ASR；依据公式（5）计算大气探测阈值；依据公式（6）计算云、气溶胶现象因子，如果P大于40%则存在云、气溶胶，如果P小于等于40%则为清洁大气。

式中：T为大气探测阈值；ACL_Rtru为校正之后的晴空ASR；Φ是用来校正潜在晴空ASR偏差的校正系数，如果地表覆盖类型为陆地则取值为0.9，如果覆盖类型为海洋则取值为1.0。

P为大气现象因子，根据该参数可实现大气探测。在算法流程中还存在其他大气探测方法，如果使用ATLAS数据使用该算法实现大气探测，则将其记录到ATL09参数Use_ATLAS中，0代表未使用该算法，1代表使用该算法实现云检测。在已知表面反射率的前提下，可计算大气柱状截面总体光学厚度，前提是返回信号强度不能被完全衰减。ICESat-2无法穿透光学厚度为3的云或者气溶胶，因此，该技术仅限于大气总光学深度之和小于3的情况^[38]，超过极限阈值，返回信号将太小或完全衰减。

光子的平均路径延迟，与散射能量的分布息息相关，而且光子可能在某一高度发生多次散射，在多种延迟路径下返回接收视场角范围内。结合光子散射、衍射理论以及卫星技术参数，Duda等人^[37]提出大气散射条件下平均延迟距离计算方案，如公式(7)~(9)所示。公式（7）表示光子散射理论框架下能量分布方程，公式（8）表示光子能量高斯分量方程，公式（9）表示光子能量平行分量方程。

式中：S为激光初始能量值；S₁为散射之后能量大小；τ为大气光学厚度；θ_m在高度为Z时视场角所能允许的最大散射角度；θ_s向前衍射峰的宽度；δ_m指在当前轨道高度、接收视场角下理论上最大的延迟距离。f_g、f_i分别表示在垂直方向、同向的能量分量。

综合考虑垂直方向上高斯分量、同向分量，可以将平均延迟距离公式（2）改写为：

Duda等人在上述原理的基础上，经过多次蒙特卡洛模拟实验，分析得到了GLAS数据多次散射误差与多个参数之间的关系，并将OD范围在0.45~1.1的数据归纳为多个查找表。通过数据反演得到的云光学厚度值在查找表中寻找对应的改正值，实现大气散射误差校正。

由大气前向散射引起的激光测距误差改正一般是结合实验室标定的蒙特卡洛模拟结果进行修正，但是蒙特卡洛模拟涉及多个参数的选取，实际生产中无法同步获取这些参数，不能实现误差改正业务化。除了蒙特卡洛模拟之外，在GLAS实际案例中，统计结果表明：在低光学厚度范围内，光学厚度与大气散射引起的测距偏差统计分布规律符合指数函数。基于此，笔者提出了一种基于指数函数的大气散射改正经验模型，初步验证表明该模型能有效提高受大气散射影响的数据精度和数据的可利用率，能为激光测高数据业务化提出一种新的解决方案。

笔者选择青海湖区域的线性体制ICESat/GLAS激光测高数据，初步构建了指数函数大气散射误差改正模型并拟合出相关参数，如公式（11）所示。

式中：OD为光学厚度；P代表测距改正值。

利用兴凯湖以及德国北威州区域数据进行交叉验证，三个区域的云光学厚度与测高偏差的斯皮尔曼相关系数依次是0.925、0.916、0.927，可以看出两个变量之间高度相关。如图5所示，随着光学厚度OD增加，拟合模型的精度逐渐降低。初步统计表明，对于水体区域，该拟合模型可将测高偏差修正到5 cm以内。

Figure 5.  Error correction model and dispersion points of optical depth and ranging deviation distribution in experimental area

经过大气散射校正之后，如果与真实值之间的误差在10 cm (水体区域)/15 cm （陆地区域）以内，将其视为修正后可利用激光点，统计结果如表2所示，实验结果表明，对于线性体制ICESat/GLAS激光测高数据，受大气散射影响的点约占总数的20%，其中OD≤2占比约50%，提出的指数函数模型对OD≤2引起的测高偏差有较好的校正作用，数据可利用率达到90%以上，总体可利用率提高了约10%。

大气散射引起的测距偏差主要是取决于以下几个因素：云高度、云光学厚度、云粒子大小、粒子形状和接收场。无论是蒙特卡洛模拟还是基于指数函数的大气散射改正方法，均需要与激光脉冲同步获取云高、云光学厚度、气溶胶等参数。在这种情况下，仅仅基于单通道、单角度激光脉冲数据反演精度有限，因此，有必要通过辅助设备同步进行大气探测，提高反演大气参数精度，进而提高大气散射误差改正精度。除此之外，受限于没有同步观测数据，反演过程中很多关键参数无法精确确定，比如激光雷达比。比如，采用蒙特卡洛模拟改正法的GLAS将激光雷达比S₁作为定值（S₁=17.8 sr），尽管这一值是在地基、Mie散射雷达数据都广泛使用的值，但是激光雷达比是与激光雷达探测波长、云（或气溶胶）类型、有效粒子半径等因素紧密相关的参数，仅采用固定值对不同类型云、气溶胶反演光学厚度具有较大的误差，进而导致大气散射改正不够精确。

而对于光子计数体制的ATLAS虽然基于第三方数据实现了云检测，但是该算法仅仅判断数据是否受到云影响，而没有实现云垂直分布判断。在没有与数据同步的第三方大气垂直分布数据源的情况下，ATLAS大气散射误差校正研究陷入了瓶颈。因此，同步获取激光测高卫星工作时段的大气相关的有效粒子半径、云相态、垂直分布等参数，进而选择更合理的激光雷达比、精确反演OD值，提高大气散射改正精度显得非常重要。

云、气溶胶等粒子是变化频繁的大气成分，变化特征观测主要围绕着粒子大小、形态、吸收和散射特性等因素进行，偏振测量是现有比较可靠的探测技术，李正强等人^[43]对星载大气探测设备现状与发展进行了梳理，如表3所示。结合前述表1中CALIPSO和CATS的激光偏振观测通道，对于激光测高卫星而言，配置偏振通道既可以提高大气散射改正精度，同时能拓展激光数据在大气方面的应用。

近几年来，对地观测卫星激光测高技术飞速发展，由线性计数体制转变为光子体制，由单波束测量发展为多波束探测，由低重频发展为高重频观测。对于对地观测卫星激光测高而言，大气散射引起的测高误差是不可忽视的，分米级甚至厘米级的精度需求对高精度的大气散射误差改正提出了更高要求。梳理激光测高卫星大气散射延迟改正技术的现状及面临的难点问题，基本可以得出如下结论：

（1）目前大气散射改正算法还处于研究阶段，并没有得到业务化应用，无论是GLAS、GEDI还是ALTAS基本采用对大气散射影响的激光点进行质量控制来确保测量精度，导致点的可利用率受到一定影响。

（2）基于指数函数的大气散射改正算法，在有效反演光学厚度OD值的前提下，可实现规模化应用，是未来业务化应用的一个重要选择。

（3）为了实现激光测高大气散射改正的业务化，同步探测大气参数非常必要且很关键。

大气散射误差改正的前提是具有较高精度的大气层结构数据，GLAS以532 nm单通道获取了75 m分辨率的层高度信息，CALIPSO以532 nm偏振双通道获取了30、60 m分辨率的层高度信息，CATS以532 nm偏振双通道获取了60 m分辨率的层高度信息。而偏振通道能提高退偏比、激光雷达比等大气参数反演精度，这些大气参数与大气散射改正息息相关。因此，后续发展国产激光测高卫星时，建议在常规激光探测通道外，可考虑补充激光大气探测通道或偏振探测模式，如同步配备不超过75 m垂直分辨率的大气探测激光通道和偏振探测通道，提高大气参数反演精度，既提高大气散射改正精度、也进一步拓展激光雷达大气探测应用。

实现大范围三维地形、森林树高、极地冰盖、湖泊水位、浅海水下地形等信息的有效获取，是未来激光测高卫星发展的必然趋势。在以陆海激光测高卫星为代表的国产激光测高卫星观测体系中，大气散射改正是卫星立项论证与后续处理应用时必须考虑的问题，希望文中的有关研究能引起更多学者的关注，最终有效提高国产激光测高卫星的设计和应用水平。