1、 引言

激光雷达LiDAR (Light Detection And Ranging)是一种集激光测距仪、全球定位系统GPS (Global Positioning System)和惯性导航系统INS (Inertial Navigation System)技术于一体的用于快速获取地面及地面目标3维信息的主动式观测系统，LiDAR技术与成像光谱、合成孔径雷达一起被列为对地观测系统(EOS，Earth Observing System)计划最核心的信息获取与处理技术(李然等，2007)。激光测高卫星在快速获取地表 3维信息、全球植被高度、两极冰盖高程变化监测等方面具有显著的优势(Schutz等，2005)。采用星载激光测高仪辅助提高卫星立体影像几何定位精度特别是高程精度，已经得到了航天摄影测量界的认可(王任享，2014；Li等，2016)。虽然中国在“嫦娥”探月系列卫星上搭载了激光测高仪，但在对地观测领域的星载激光测高卫星发展方面起步较晚，计划于2017年发射的高分七号卫星上将同时搭载光学立体相机和激光测高仪，将激光测高仪获得的激光足印点作为高程控制点辅助立体影像测图(唐新明和高小明，2012；唐新明等，2016)。此外，中国已经初步建成高精度的平面控制点数据库，基本解决了1:50000的中国DOM快速生产及更新的问题(喻贵银，2014)，而高程控制库的建设相对较弱，控制数据的来源及精度还有待深入的研究(张重阳等，2014)，中国自主研发的覆盖全球的高程控制库更是处于空白状态。

采用激光测高卫星获得全球的高程控制数据是一种比较有效的方式，美国早在1997年和1999年先后采用搭载了激光测高仪的航天飞机SLA-1(Shuttle Laser Altimeter)和SLA-2进行了两次飞行实验(Garvin等，1998)，获得了57°N--57°N范围内的全球高程控制点数据，其激光足印大小为直径100 m，测高精度约1 m。2003年成功发射的ICESat (Ice，Cloud，and land Elevation Satellite)卫星搭载了地球激光测高系统GLAS (Geo-science Laser Altimetry System)，GLAS的足印大小约70 m，地面间隔约170 m，是目前唯一一颗对地观测的激光测高卫星，已于2009年停止工作(Schutz等，2005)。ICESat卫星共采集了全球约20亿个高程数据，在极地冰盖变化监测、全球森林生物量估算等方面得到广泛应用(Wang等，2011；Zwally等，2002；李建成等，2008)，其中大量高精度的陆地高程数据，也被作为高程控制数据用于ASTER-DEM、TerraSAR-DEM等全球DEM数据生产(鄂栋臣等，2009；González等，2010)。由于激光测高卫星发射的激光脉冲在地表对应的光斑较大，一般直径达几十米，即使美国将于2018年和2020年发射的ICESat-2以及LIST (LIdar Surface Topography)卫星的激光足印大小也分别为10 m和5 m，已经显著区别于传统意义上的“点”。此外，激光在往返大气的传输过程中，其测距值易受大气、云层、地表反射率、硬件动态响应范围等因素的影响，导致其测距精度以及最终高程可靠性不高。因此，需要对激光测高卫星获得的数据进行一定规则筛选后，才能作为高程控制数据使用。

ICESat/GLAS数据产品共份3级15类，其中GLA14标准数据产品中包含激光足印3维坐标值、激光回波饱和度改正、增益值等80多个参数，其中部分参数可以作为GLAS激光高程控制点筛选的参照依据(Zwally等，2005)。Huang等人(2012)针对中国伊犁河谷地区的GLAS测高数据，利用地表反射率、增益值和饱和度3种参数进行粗差剔除，并利用筛选后的数据检测研究区内湖泊的水位变化。Wang等人(2012)提出了湖泊激光足印(Lake Water Footprint)的概念，并利用回波波形宽度、波峰数、偏度、峰度以及反射率等参数，并对西藏日喀则地区的佩枯措湖泊2003年-2009年的GLAS数据进行了统计分析。万杰等人(2015)利用GLAS激光测高数据评价青藏高原地区的SRTM数据精度时也利用饱和度、反射率和增益值等参数对GLAS数据进行了处理。黄海兰在极地冰盖研究中采用6种参数，包括增益值，高斯拟合标准差，反射率，饱和度，指向偏差，姿态标志和波峰数，对每个参数的数据保留率进行了统计分析(黄海兰，2011；黄海兰等，2012)。Duong等人(2009)选择高程使用标记、增益值和饱和度改正标记3个参数剔除GLAS激光测高数据误差较大值，并利用高精度的机载激光数据对筛选后的数据进行了精度分析。

本文提出了一种多准则约束的激光高程控制点筛选算法，采用90 m分辨率SRTM-DEM对GLAS数据进行高程粗差剔除，利用云量、饱和度参数、地表反射率、回波波形的波峰数、波形拟合后高斯函数的标准偏差等设定一定阈值对ICESat/GLAS数据进行多重准则筛选研究，将筛选后的激光足印点数据作为高程控制。同时采用天津、河北等地区的ICESat/GLAS激光测高数据以及高精度的DEM成果数据进行了实验验证分析。相应的算法研究及实验结果对于开展国产激光测高数据处理及应用，以及建立全球的高程控制库将会有一定的参考价值。

2、 多准则约束算法

利用SRTM-DEM数据，结合GLAS的GLA14数据产品中的有关测距属性信息以及GLA01中的回波波形特性(Brenner等，2000)，提出了一种多准则约束的星载激光高程控制点筛选算法，实现从GLA01的波形数据以及GLA14的L2级产品中筛选出可靠性好、高程精度高的激光足印点作为高程控制点，用于开展光学影像立体区域网平差或DEM数字高程模型高程精度评价应用。其总体流程如图 1所示，其中，i_ElvuseFlg为高程可用性标记参数，i_satCorrFlg为饱和度改正参数，i_sigmaAtt为姿态质量指标参数，i_reflctUncor为地表反射率参数，i_gval_rcv为大气散射增益参数，i_FRir_qaFlag为云量描述参数。

    (2.1) 基于SRTM-DEM数据的GLAS粗差剔除

美国2000年2月发射的“奋进号”航天飞机历时11天，通过干涉雷达的方式获得南北纬60°间约占全球80%陆地面积的DEM数据，也简称SRTM-DEM (Shuttle Radar Topography Mission)。目前公开版的SRTM-DEM格网间隔约90 m，绝对高程精度是±16 m，绝对平面精度±20 m，高程基准为EGM96大地水准面，平面基准为WGS84(Brown等，2005；Gesch等，2006)。由于SRTM-DEM高程精度具有较高的一致性和可靠性，因此可将其作为GLAS激光足印高程粗差剔除的参考数据。

虽然ICESat/GLAS的高程精度很高，其公布的精度达15 cm，但前提条件是在没有植被、地形平缓、大气透射率良好的情况下，当激光束在传播过程中受云层以及复杂地形影响时，其精度会有一定下降甚至完全错误。考虑到GLAS激光足印大小为直径70 m，平面精度约10 m，而SRTM-DEM的格网大小为90 m，因此可以直接按激光足印的平面坐标在SRTM中内插出高程值，然后与激光足印高程进行比较。对于高程较差绝对值大于16 m的点，即可作为粗差点进行剔除，这样保证那些受云、复杂地形影响严重的高程错误或精度明显不高的激光点被滤掉。

    (2.2) 激光测距属性参数约束的粗粒度筛选

对地观测的激光测高卫星接收的激光脉冲需要往返经过大气层，大气分子及气溶胶会对激光脉冲产生散射、能量衰减、光斑漂移、波形展宽等多重影响；地表的高反射率和云层的前向散射等会引起回波波形过饱和，使得测距精度下降；此外，卫星姿态的调整以及太阳能帆板抖动等也会使激光束的指向精度下降，影响最终激光足印定位精度。为了保证激光测距值的有效性和精确性，需要引入云量、增益值、姿态可靠度、地表反射率、波形饱和度等参数进行描述，同时设定一定的阈值进行筛选，最终保留测距数据有效且精度较高的激光足印点。

本文中主要采用GLA14级产品中的高程可用性标记参数、饱和度改正参数、姿态质量指标参数、地表反射率参数、大气散射增益参数以及云量描述参数等6种参数进行组合对经SRTM-DEM粗差剔除后的激光足印点进行粗粒度的筛选，具体流程见图 1的中间虚线框部分。

其中GLA14产品中每个激光足印点对应一个高程可用性标记参数值，0表示高程数据有效，1表示高程数据无效。与饱和度改正相关的参数有两个，饱和度改正标记(i_satCorrFlg)和饱和度高程改正值(i_satElevCorr)。饱和度改正标记为0或1表示信号受饱和影响可忽略，2表示需要进行饱和改正，其他值表示数据不可用。本文主要保留饱和度改正标记为0或1的激光点。姿态质量指标参数描述的是GLAS测高数据的姿态测量质量，值为0时表示数据良好，值为50表示警告，值为100表示数据质量较差。本文中剔除姿态质量指标50或100的数据。

此外，反射率参数(i_reflctUncor)可以用来剔除异常值，激光足印的反射率由接收激光脉冲的能量除以发射激光脉冲的能量得到，而接收能量通常情况下比发射能量小，因此反射率参数的理论值不可能大于或等于1。本文将反射率大于0.5的GLAS激光足印点剔除。

GLAS激光脉冲受到云层或水蒸气等大气成分引起的散射作用时，返回的信号能量会减小，导致波形变宽，影响测高数据的精度。而增益值参数(i_gval_rcv)和接收激光脉冲的能量成反比，当接收能量较低时，i_gval_rcv的值比较大；当接收能量较高时，i_gval_rcv的值较小。利用增益值参数可以剔除受到大气散射等作用影响的数据，本文保留增益值小于100的数据。云量描述参数(i_FRir_qaFlag)表示在激光往返大气层时天空云量的多少，按GLAS产品的约定，i_FRir_qaFlag=15时表示无云，若小于15则表示有一定云量。为了最大限度减少云层对激光测距值的影响，本文中只保留无云时工作的激光足印点。

    (2.3) 回波波形特征参数约束的精细化筛选

经过2.2节的测距属性参数约束筛选后的激光足印点其测距精度可靠，回波波形也可信，但是由于激光测高卫星的地面激光足印一般为大光斑，对于足印大小约70 m的ICESat/GLAS，当地形有起伏或有建筑、植被等地物时，其单一的高程值无法精确代表足印内地形信息。图 2和图 3都是经初步筛选后保留下来测距值可信的激光足印点，图 2所对应的激光足印正好位于平坦区域，图 3所对应的激光足印则位于森林区域，虽然其距离值没有问题，但显然后者的激光足印点不能作为高程控制点使用，因此还须依赖激光回波波形信息进行精细化筛选。

由于激光发射脉冲在沿发射方向和垂直发射方向均为类高斯分布，而地形地物对回波影响则主要在回波峰值、展宽、波峰对应的位置等，所以对于回波波形仍然可以采用高斯函数进行拟合分析，而对于复杂的地形则需要多个高斯函数组合才能达到最佳拟合效果。

对于实际的回波波形，可以采用式(1)进行拟合(Brenner等，2000)，

$ w(t) = \varepsilon + \sum\limits_{m = 1}^{Np} {{W_m}} ,\;{W_m} = {A_m}{{\rm{e}}^{\frac{{ - {{(t - {t_m})}^2}}}{{2\sigma _m^2}}}} $

(1)

式中，w(t)为t时刻的回波能量值，W_m为第m个高斯函数在t时刻的回波能量值，N_p为高斯函数的个数，ε为噪声，t_m为第m个高斯函数的均值，σ_m为第m个高斯函数的标准差，A_m为第m个高斯函数的峰值。

如果激光回波波形为单峰且展宽在一定范围内，则其对应的地形是平坦的水平面。式(2)对回波波形脉宽与发射波形脉宽、激光发散角、卫星高度等参数的关系进行了定量描述(Gardner，1992)。

$ {\tau _s} = \sqrt {(\tau _h^2 + \tau _f^2) + \frac{{4{h^2}{{\tan }^2}\theta }}{{{c^2}}}({{\tan }^2}S + {{\tan }^2}\theta )} $

(2)

式中，τ_s为回波波形脉宽，τ_f为发射波形脉宽，τ_h为硬件引起的波形展宽，h为卫星高度，θ为激光发散角，c为光速，S为地表斜率。

根据FWHM (Full Width at Half Maxium)脉宽τ的定义：脉冲峰值降低至一半时所对应的两个时刻差，即$w(\tau /2) = {\displaystyle\frac{1}{2}}A$，可得高斯函数的标准差为

$ \sigma = \frac{1}{{\sqrt {2\ln 2} }}\tau $

(3)

按ICESat/GLAS的有关参数进行估算，其中发射脉宽τ_f为6 ns，σ_h约1 ns，θ约110 urad，h约600 km，地面足印大小约70 m，如果其地形起伏不超过0.5 m，则S约为0.007 rad，根据式(2)可求得τ_s约为6.8 ns，由式(3)可得σ=2.89 ns；如果地形起伏约1.0 m，τ_s约为7.5 ns，则σ=3.2 ns。在国标GB/T 12341-2008中对1:50000高程控制点精度要求为1.0 m，为从GLAS数据筛选出高程误差优于1.0 m的作为高程控制点，本文结合GLA01的波形数据以及GLA14中的波形拟合参数，利用回波波峰数为1且拟合的高斯函数标准差σ≤3.2 ns作为约束条件对激光足印进行精细化的筛选。

3、 实验数据

本文选取了天津中东部地区和河北北部的ICESat/GLAS激光测高数据开展实验，参考数据选择该实验区的高精度DEM数据。其中天津地区高程参考数据为1:2000的DEM成果图，该DEM成果平面坐标系统为WGS84，高程基准为WGS84大地高，格网间隔约1.5 m。河北地区高程参考数据由于条件限制，选择的是1:50000的DEM成果图，平面坐标系统为WGS84，高程基准为WGS84大地高。实验数据的总体情况介绍见表 1。

将ICESat/GLAS数据与DEM进行精度对比时，从DEM上获取以GLAS激光足印点的经纬度为中心、35 m为半径范围内的所有高程值并取平均值，将该平均高程值与筛选后的GLAS高程值进行对比。

ICESat/GLAS数据的高程参考基准为TOPEX/Poseidon椭球下的大地高，TOPEX/Poseidon椭球的基本参数为长轴a=6378136.300000 m，偏心率e=1/298.25700000，扁率f=0.081819221456，同一点在该椭球和WGS84椭球上坐标在纬线方向差别较小、经线方向无改正，因此在进行数据转化时可以只考虑对高程值进行改正。可采用式(4)进行高程的转换，转换误差在mm量级(李建成等，2008)。

$ \left\{ {\begin{array}{*{20}{l}} {dh = - {{\cos }^2}Bda - {{\sin }^2}Bdb}\\ {{H_G} = H + dh} \end{array}} \right. $

(4)

式中，B是纬度，da、db表示的是WGS84椭球和TOPEX/Poseidon椭球长轴和短轴之间的差异，dh为高程变换修正值，H、H_G分别表示的是GLAS变换前、后的高程值。为了将GLAS转换后的WGS84大地高和SRTM的高程值进行比较，文章还利用了公开的EGM96大地水准面模型，利用H_g=H-N计算EGM96下的正高，H为大地高，N为大地水准面差距，可从EGM96模型中获得。

采用SRTM-DEM作为参考数据，对实验区的GLAS数据进行粗差剔除，天津和河北实验区的数据剔除率分别为6.11%和5.02%，说明GLAS数据中粗差点存在但数量不是特别多。在进行粗差剔除后，采用激光测距属性参数约束和波形特征参数约束进行进一步的筛选，筛选结果如表 3所示。

对经过激光测距属性参数约束筛选后的激光足印点继续采用回波波形特征参数进行精细化的筛选。为了比较波形拟合后高斯函数不同的标准差阈值筛选条件下，最终GLAS激光点的高程精度，本文采用了波峰数N_P=1和σ₁=8 ns、5 ns、3.2 ns等不同阈值进行对比，最终统计结果如表 4所示。

4、 结论

本文为了剔除ICESat/GLAS激光测高数据中因云层干扰、大气散射、激光回波接收端仪器响应饱和、不可靠的姿态测量、复杂地形地物等因素引起的精度不高的点，保留高精度的激光足印点作为高程控制点使用，提出了一种多准则约束的筛选算法，结合GLA14数据的特点以及ICESat/GLAS的有关硬件参数，对各约束参数的阈值进行了讨论，并利用两个实验区开展了实验验证工作。

由表 3可以看出，利用高程可用性标记、饱和度改正、姿态质量指标、反射率值、增益值以及云量等属性数据对GLAS激光测高数据进行无效值和测距可靠性较低点筛选时，由于增益值较高而被剔除掉的数据最多，其次是饱和度改正标记参数剔除的数据点数量，其他类参数对于数据的剔除率基本相当。例如，河北实验区，因增益值和云量约束剔除了约67.55%的激光点，与DEM相比中误差也从8.859 m减小到5.484 m。由此表明，因云层、大气或地表低反射率等引起的激光高增益对测距值的可靠性影响非常大。为提高激光测距精度和可靠性，在数据处理时应考虑云量的影响，也需从硬件研制以及仪器对不同地表反射率的动态响应范围等方面进行深入探讨。

结合表 3和表 4的统计结果可知，当N_p=1，σ=8 ns时，天津实验区的GLAS数据与1:2000 DEM相比中误差从1.162 m提高到0.803 m，河北实验区与1:50000 DEM数据的中误差由5.484 m减小到3.468 m，说明筛选后GLAS数据的精度有进一步提高。当波形展宽σ进一步设置为3.2 ns时，处于平原区的天津GLAS高程精度提高到0.725 m，处于山区的河北实验区GLAS高程精度提高到3.288 m，考虑到DEM本身的精度，筛选后的天津GLAS高程控制点的精度优于0.5 m，可以达到1:50000以及1:10000部分区域的立体测图高程控制需求(GB/T, 2008, 2012)。而河北实验区由于DEM本身精度的限制，其最终GLAS高程控制点的精度还有待进一步的分析和验证，但文献(Li等，2016)的结果表明可以满足1:50000的高程控制需求。

未来几年，搭载了激光测高仪的国产高分七号以及陆地生态系统碳监测卫星将会陆续发射，两颗卫星上的激光测高仪的主要用途之一就是获取全球的高程控制点。借鉴ICESat/GLAS数据处理经验，结合国产激光测高卫星的特点，开展全球范围高程控制库建设将是下一步的研究方向。