星载激光测高仪采用主动遥感探测方式，能够精确获取地面高程信息，有效地弥补了传统光学遥感系统工作条件苛刻以及三维探测能力弱等的不足。激光作为一种新型的大气和地表遥感手段，未来将广泛地应用于地形测绘、气象、海洋、环境资源等遥感应用领域[1-5]。由于星载激光测高仪轨道、功率和能量等问题，能获取的激光测高数据是采样间隔很大且平面精度较低的大光斑回波波形数据。相对于较高的高程测量精度，目前GLAS等星载激光测高仪的平面定位精度仅为十几米到几十米[6]，不能直接进行高精度定位，且每个单独的光斑回波数据也不能直接与地形实现空间配准，限制了星载激光遥感的应用。目前星载大光斑激光数据还局限在极地冰盖、海洋、植被等空间范围较大且空间分辨率较低的宏观性遥感应用领域内[7-8]。即将发射的ICESat-2，则设计通过增加光束数量，减小光斑直径到10 m以内来提高空间分辨率，减少坡度等对地形观测的影响[9]，但在空间尺度上还是与高分辨率遥感图像有很大差距。在星载激光测高仪检校和测图应用方面[10-13]，可以通过激光测高数据与地形数据匹配，获取控制点，准确找到激光脚点位置，并根据地面控制数据检校激光测高仪的几何定位模型，实现高精度三维定位。另外，多源遥感数据融合应用时，以激光测高数据与地形数据匹配为中介，可以实现激光数据与遥感图像等其他地理数据的空间配准。因此，研究大光斑星载激光测高数据与地形的匹配方法，充分挖掘光斑内所包含地形结构信息的作用，提高星载激光测高仪几何定位能力，实现大光斑激光高精度对地观测具有重大意义。
在地形探测中，星载激光测高仪发射并接收激光光束后向散射信号，可以直接提供观测目标的距离参数和高程信息，再联合由卫星姿轨参数建立的指向矢量，计算得到光斑的位置坐标[7-14]。激光测高数据通过由大量回波观测数据生成的粗格网低分辨率数字表面模型 (digital surface model，DSM) 与已有的高精度DSM配准，来实现几何校正[15-16]。在激光和立体相机联合观测体制下，可利用多源数据联合平差来提高激光观测值平面定位精度[5, 17-18]。但上述方法都不是单独激光光斑回波波形与地形之间的匹配。到目前位置，在摄影测量领域还没有星载大光斑激光回波波形数据与地形匹配的相关研究报道。
本文根据星载大光斑激光测高数据回波信号包含了地形结构信息的特性，在统计信息空间计算激光回波波形及参考DSM的累计剩余分布函数，以交叉累积剩余熵[19-20]为相似性测度，匹配激光回波波形及地形的空间位置。该方法规避了激光回波波形数据与DSM之间的维度差异，在地形参考下直接获取激光光斑空间位置。经过对试验区GLAS01激光测高波形数据与DSM数据的匹配试验，证实本方法能够实现星载大光斑激光测高数据与地形的匹配。
1 星载大光斑激光测高数据与DSM的统计信息特征 1.1 大光斑激光回波信号DSM的统计信息特征星载大光斑激光测高波形数据是由高速模数转换器对探测激光回波进行高分辨率采样和数字量化，通过采样计数统计获得的。如图 1所示，回波信号的横轴是时间 (单位为ms)，纵轴是接收信号数字化采样后量化的结果 (单位为V)。回波信号横轴的正方向表示时间增加，纵轴正方向表示能量强度增大。
|
| 图 1 大光斑激光测高回波信号示意图 Fig. 1 Sketch map of large footprint laser ranging echo waveform |
回波信号波形与时间采样间隔和能量量化方式有关，这两个因素分别作用于波形的横轴和纵轴，决定了波形的形状。回波信号开始接收的部分是背景噪声，一般当信号的强度达到某一阈值时才被视为地面真实观测的有效信号。回波信号是距离的观测值，表示高程值高的信号先被接收到。波形的峰代表光斑内主要的高程值，全波形表示光斑范围内的地形结构信息。回波波形可以通过改变横轴采样间隔，规定化纵轴数值范围，投影到统计空间域，利用统计学手段处理分析。
1.2 DSM高程直方图DSM产品是二维的图像，每个像素的灰度值表示该地理位置的高程值。为了与大光斑激光回波信号维度统一，本文利用高程直方图来表示DSM中的信息。如图 2所示，图 2(a)为DSM图像，图 2(b)为图 2(a)中DSM的高程直方图。高程直方图的横轴为高程值 (单位为m)，纵轴为像素数量。高程直方图的横轴正方向表示高程值增大，纵轴正方向表示该高程值像素的数量增加。
|
| 图 2 DSM高程直方图示意 Fig. 2 Elevation histogram of DSM |
对于相同地区，高程直方图与DSM图像地面采样间隔 (ground sampling distance, GSD) 或地面分辨率和图像像素值量化位数相关。这两个因素分别对应1.1节中时间采样间隔和能量量化方式，作用于高程直方图的横轴和纵轴，决定了直方图的形状。通过计算高程直方图，可以将DSM图像降维到一维信号，再与1.1节中所述一致，通过改变横轴采样间隔，规定化纵轴数值范围，投影到统计特征空间，利用统计学手段处理分析，从而实现激光波形数据与DSM图像信号维度的统一。
1.3 交叉累积剩余熵累计剩余熵 (cumulative residual entropy，CRE)[19-20]仅考虑随机变量的概率分布，统计大于某一阈值的所有信号信息，保持了信号的连续性，有效克服了噪声对局部极值的影响，在连续域和离散域都可以适用。
定义在R上的向量X=(X1, X2…, Xn), n为自然数且n≥1，的累计剩余熵ε(X) 为
(1)
式中，|X|=(|X1|, |X2|, …, |Xn|)，λ=(λ1, λ2，…, λn) 且λi≥0，i=1, 2, …, n。若|X|>λ，则|Xi|>λi。P(|X|>λ) 为|X|>λ的概率，在CRE的计算中称为多元生存函数。
当存在与X定义域相同的向量Y时，则X在条件Y下的条件CRE为
(2)
X与Y的交叉累积剩余熵C(X, Y) 定义为
(3)
计算1.1节中的星载大光斑激光测高回波波形数据和1.2节中的DSM高程直方图两种信号之间的交叉累积剩余熵，首先要在物理意义上统一两种信号的向量方向，即将波形数据横轴反转，按时间从长到短排列，与高程直方图表示的高程高低方向相对应。此时波形数据向量用SL表示，高程直方图向量用SE表示，则SL与SE的交叉累积剩余熵C(SL, SE) 为
(4)
式中, l和e分别为SL与SE中的元素；PL(l) 和PE(e) 为SL与SE的边缘概率密度；P(l, e) 为SL与SE的联合概率密度。C(SL, SE) 值越大，表示SL与SE之间的相似性约大。
2 本文方法匹配流程本文波形数据与地形匹配方法如图 3所示。
|
| 图 3 基于CCRE的大光斑激光波形数据与地形匹配流程 Fig. 3 Big footprint laser waveform and DSM matching method |
基于交叉累积剩余熵的星载激光测高仪大光斑波形数据与地形匹配方法流程如下：
(1) 波形数据提取，减少背景噪声对真实信号的影响。
(2) 根据波形数据的初始地理位置获取DSM匹配搜索区域。
(3) 对去噪后的波形数据按照时间从长到短重新排列，获得与DSM高程直方图一致的高程由低到高的方向矢量。
(4) 设置波形数据和高程直方图采样间隔，规定化每个采样间隔的数值，把波形数据和DSM图像投影到统计特征空间，并确定匹配窗口大小、相似性测度阈值等匹配算法参数。
(5) 遍历DSM搜索区，计算每个匹配窗口的CCRE值，其中CCRE值最大的匹配窗口的中心为匹配点，即为波形数据所表示的光斑重心。
3 试验与分析采用美国ICESat卫星星载激光高度计GLAS (geoscience laser altimeter system) 获取的大光斑激光对地观测回波波形数据GLAS01产品和地面DSM测绘产品匹配，来对本文提出方法进行试验分析。GLAS01数据产品从美国国家冰雪数据中心 (The National Snow and Ice Data Center，NSIDC) 网站 (http://nsidc.ors/data/icesat) 下载得到，光斑直径约为70 m，相邻光斑距离约为170 m。根据GLAS01产品中记录的光斑号等信息，在与之相对应的GLAS14数据产品中获取每一个光斑对应的经纬度坐标，并用GLAS14数据生成的DSM与参考DSM配准，对每个光斑的平面坐标进行修正。为了实现波形数据与地形的匹配，DSM图像的分辨率要尽可能多的高于激光光斑大小，所以本文试验DSM图像空间分辨率为5 m×5 m，为1:50 000标准产品。试验区域为缅甸东北部山区，WGS-84坐标范围为97.250 3°E—97.257 3°E，22.772 1°N—22.785 1°N。
本文激光回波波形试验数据集包括6个光斑，具体信息如表 1所示。
| 序号 | 经度E/(°) | 纬度N/(°) | 采样数 |
| 1 | 97.256 327 | 22.775 236 | 1000 |
| 2 | 97.253 701 | 22.775 656 | 1000 |
| 3 | 97.252 341 | 22.785 008 | 1000 |
| 4 | 97.257 274 | 22.769 096 | 1000 |
| 5 | 97.256 798 | 22.772 161 | 1000 |
| 6 | 97.250 326 | 22.783 437 | 1000 |
GLAS数据的坐标使用Topex/Poseidon椭球，该坐标系与WGS-84椭球差异在平面坐标上仅相差几厘米，与试验数据空间分辨率相比可以忽略。而两个坐标系下高程差异ΔH由纬度B及两个椭球长短半轴 (a1, b1)、(a2, b2) 的函数ΔH=-[(a2-a1) cos2B+(b2-b1) sin2B]来计算，ΔH在一个光斑内是系统性的[21]。对于波形数据来说相当于整个波形在纵轴上平移一个常数值，在映射到统计特征空间时没有影响。
3.1 波形数据提取本文GLAS01试验数据每个光斑回波信号包含1000个采样点，一般前1/3的采样点不包含真实信号。本文选择前300个采样点计算回波信号的背景噪声。根据背景噪声来计算确定真实信号波形起止位置的阈值。阈值通常为背景噪声均值与4倍标准差的和[22]。激光回波波形试验数据集的波形数据提取结果如图 4所示，图中横坐标表示采样数，纵坐标为回波信号强度，单位为伏特 (V)。
|
| 图 4 波形数据提取 Fig. 4 Waveform data detection |
图 4中波形横轴为采样间隔，纵轴为能量，是回波信号电压值。从图 4可以看出，波形提取结果中数据1, 2, 3, 4, 6均较好地保留了波型中的主要信号，数据5相对于其他数据波形较平滑，从物理意义上来说表示光斑中地形结构不明显，所以波形提取结果仅保留了主要波峰，去掉了一些次要的信息。
3.2 统计特征空间变换由1.1节和1.2节可知，激光回波波形数据和DSM图像高程直方图纵坐标物理意义不一致，需要统一到一个量化范围，本文选择8 bit即0~255。匹配窗口应与光斑尺寸相对应，边长为至少70 m/5 m=14(像素)，为了防止直方图疏散陷入局部极值，将高程直方图横轴压缩为8级，则激光回波波形数据横轴也相应压缩为8级。经过3.1节中波形提取后的波形数据转换到统计特征空间后如图 5所示，纵轴表示特征空间变换后的量化的值。
|
| 图 5 激光试验数据统计特征空间分布 Fig. 5 The laser data set in statistics domain |
由图 4和图 5可以看出，激光回波波形数据转换到统计特征空间后保持了原始波形形状。其中激光回波波形数据由1.3节中所述重新排列。而DSM图像映射到统计特征空间为1.2节中的直方图。此时两种数据实现了形式上的统一。
3.3 波形与地形匹配在统计特征空间采用CCRE为相似性测度对激光回波波形数据与DSM图像进行匹配，为了获取更多的地形信息避免地形起伏对光斑大小的影响，匹配窗口大小为85 m×85 m，对应DSM图像为17像素×17像素，把GLAS01波形数据对应的地理坐标，换算得到DSM图像的像素坐标，作为真值评价匹配精度。匹配结果见表 2。
| 试验数据 | X方向精度 (经度)/像素 | Y方向精度 (纬度)/像素 | 平面精度 /像素 |
| 1 | -0.39 | -0.99 | 1.06 |
| 2 | -0.29 | -0.7 | 0.76 |
| 3 | -0.2 | -0.79 | 0.81 |
| 4 | -0.34 | -0.84 | 0.91 |
| 5 | 0.19 | -1.03 | 1.05 |
| 6 | 0.38 | -0.55 | 0.67 |
3.4 匹配结果分析
试验中的匹配结果为最大CCRE值所在的匹配窗口的中心像素，但试验数据GLAS01的地理坐标投影到DSM图像时一般不是整像素，所以表 2中匹配精度为非整像素。从表 2中可以看出，匹配结果沿经度方向优于沿纬度方向，沿经度方向为0.3个像素左右，沿纬度方向具有一定系统性，接近1个像素，整体匹配精度为1个像素左右。
由图 5可以看出，数据1与数据4激光回波波形数据只有一个主峰，即光斑内部地形结构简单，地物起伏较小，投影到统计特征空间分布单一，所以匹配精度较低；而数据5激光回波波形数据波形没有明显的峰也不平滑，光斑内地形结构信息不明显，所以匹配精度较低；数据2与数据3激光回波波形数据波形较平坦，即光斑内部地形高程值连续分布，投影到统计特征空间分布较分散，所以相对匹配精度高于数据1与数据4；数据6激光回波波形数据有几个明显的峰，光斑内地形起伏较大，地形结构明显，所以匹配精度最高。
在本文方法中，基于地理位置设置搜索区域后，大光斑激光回波波形与DSM图像异源数据匹配仅考虑数据的统计信息，所以匹配误差主要来源于匹配方法。交叉累积剩余熵作为相似性测度，属于基于窗口的匹配方法[23]，以匹配窗口内信息作为匹配基元，即使克服了噪声的影响，与其他基于窗口的匹配方法一样，可能陷入局部极值，造成匹配错误。由于激光回波波形数据与DSM图像数据在统计特征空间中，在匹配窗口内作为整体进行匹配，在统计特征空间中向量每一个元素无法与空间位置对应，所以当搜索区出现相同地形统计特征分布区域时会存在误匹配。只有地形结构信息明显时，且根据先验知识获取的搜索区域位置精度较高时，才会取得较好的匹配结果。另外本文试验结果中出现的沿纬度方向的系统性误差可能是DSM图像的系统性误差造成的。
4 结论本文根据星载大光斑激光测高数据回波信号和DSM图像包含地形结构信息的特性，将激光回波波形数据与DSM图像投影到统计特征空间，消除激光回波波形数据与DSM之间的维度差异，以交叉累积剩余熵为相似性测度匹配激光回波波形及地形的空间位置。通过GLAS01激光测高回波波形数据与DSM图像数据的匹配试验，证实本方法能够实现星载大光斑激光测高回波波形数据与地形的匹配，匹配结果可用于星载激光测高仪几何标定以及与其他遥感数据融合应用。
经过试验分析，本文方法激光回波波形数据与DSM图像空间位置的匹配精度在1个像素以内。但是在统计特征空间，激光回波波形数据作为一个整体进行匹配，光斑内地形结构信息不包含绝对地理信息参考，所以当搜索区出现相同统计特征分布时会误匹配。如何在统计特征空间加入地形结构约束，提高匹配精度，同时扩大搜索区域，还需要进一步研究。
| [1] | 李松. 星载激光测高仪发展现状综述[J]. 光学与光电技术, 2004, 2(6): 4–6. LI Song. Recent Progress of Spaceborne Laser Altimeter System[J]. Optics & Optoelectronic Technology, 2004, 2(6): 4–6. |
| [2] | 李然, 王成, 苏国中, 等. 星载激光雷达的发展与应用[J]. 科技导报, 2007, 25(14): 58–63. LI Ran, WANG Cheng, SU Guozhong, et al. Development and Applications of Spaceborne LiDAR[J]. Science & Technology Review, 2007, 25(14): 58–63. |
| [3] | 刘焕玲. 利用ICESat激光测高数据探测南极冰盖的高程变化[D]. 阜新: 辽宁工程技术大学, 2011. LIU Huanling. The Study of Ice Sheet Height Change in Antarctic by Using ICESat Laser Altimeter Data[D]. Fuxin:Liaoning Technical University, 2011. |
| [4] | 岳春宇, 郑永超, 陶宇亮. 星载激光测高仪辅助卫星摄影测量浅析[J]. 航天返回与遥感, 2013, 34(4): 71–76. YUE Chunyu, ZHENG Yongchao, TAO Yuliang. Study on Space-borne Laser Altimeter Supported Satellite Photogrammetry[J]. Spacecraft Recovery & Remote Sensing, 2013, 34(4): 71–76. |
| [5] | WU Bo, HU Han, GUO Jian. Integration of Chang'E-2 Imagery and LRO Laser Altimeter Data with a Combined Block Adjustment for Precision Lunar Topographic Modeling[J]. Earth and Planetary Science Letters, 2014, 391: 1–15. DOI:10.1016/j.epsl.2014.01.023 |
| [6] | BAE S, WEBB C, SCHUTZ B. GLAS PAD Calibration Using Laser Reference Sensor Data[C]//Proceedings of AIAA/AAS Astrodynamics Specialist Conference and Exhibit. Providence, RI:AIAA, 2004:1-10. |
| [7] | ZWALLY H J, SCHUTZ B, ABDALATI W, et al. ICESat'S Laser Measurements of Polar Ice, Atmosphere, Ocean, and Land[J]. Journal of Geodynamics, 2002, 34(3-4): 405–445. DOI:10.1016/S0264-3707(02)00042-X |
| [8] | NELSON R. Model Effects on GLAS-based Regional Estimates of Forest Biomass and Carbon[J]. International Journal of Remote Sensing, 2010, 31(5): 1359–1372. DOI:10.1080/01431160903380557 |
| [9] | ABDALATI W, ZWALLY H J, BINDSCHADLER R, et al. The ICESat-2 Laser Altimetry Mission[J]. Proceedings of IEEE, 2010, 98(5): 735–751. DOI:10.1109/JPROC.2009.2034765 |
| [10] | 唐新明, 李国元, 高晓明, 等. 卫星激光测高严密几何模型构建及精度初步验证[J]. 测绘学报, 2016, 45(10): 1182–1191. TANG Xinming, LI Guoyuan, GAO Xiaoming, et al. The Rigorous Geometric Model of Satellite Laser Altimeter and Preliminarily Accuracy Validation[J]. Acta Geodaetica et Cartographica Sinica, 2016, 45(10): 1182–1191. DOI:10.11947/j.AGCS.2016.20150357 |
| [11] | 王任享, 王建荣. 二线阵CCD卫星影像联合激光测距数据光束法平差技术[J]. 测绘科学技术学报, 2014, 31(1): 1–4. WANG Renxiang, WANG Jianrong. Technology of Bundle Adjustment Using Two-line-array CCD Satellite Image Combined Laser Ranging Data[J]. Journal of Geomatics Science and Technology, 2014, 31(1): 1–4. |
| [12] | 李鑫, 廖鹤, 赵美玲, 等. 激光测绘卫星对不同地表形貌探测能力分析[J]. 测绘学报, 2014, 43(12): 1238–1244. LI Xin, LIAO He, ZHAO Meiling, et al. Research on LiDAR Surveying Satellite Detection Capacity for Different Terrains[J]. Acta Geodaetica et Cartographica Sinica, 2014, 43(12): 1238–1244. DOI:10.13485/j.cnki.11-2089.2014.0188 |
| [13] | 李国元, 唐新明, 王华斌, 等. GLAS激光测高数据辅助的资源三号三线阵区域网平差研究[C]//第三届高分辨率对地观测学术年会分会论文集. 北京: 中国科学院重大科技任务局, 2014. LI Guoyuan, TANG Xinming, WANG Huabin, et al. ZY-3 Three Line Array Images Block Adjustment Supported by GLAS Data[C]//The 3rd Annual Conference on High Resolution Earth Observation. Beijing:Major Science and Technology Bureau of the Chinese Academy of Sciences, 2014. |
| [14] | 李磊, 郑永超, 彭凤超, 等. 地形测绘激光成像雷达技术研究[J]. 红外与激光工程, 2006, 35(S): 294–298. LI Lei, ZHENG Yongchao, PENG Fengchao, et al. Research of Three-dimension Imaging LiDAR on Land Topography[J]. Infrared and Laser Engineering, 2006, 35(S): 294–298. |
| [15] | 范春波, 李建成, 王丹, 等. ICESat/GLAS激光脚点定位及误差分析[J]. 大地测量与地球动力学, 2007, 27(1): 104–106. FAN Chunbo, LI Jiancheng, WANG Dan, et al. ICESAT/GLAS Laser Footprint Geolocation and Error Analysis[J]. Journal of Geodesy and Geodynamics, 2007, 27(1): 104–106. |
| [16] | CARABAJAL C C, HARDING D J, SUCHDEO V P. ICEsat LiDAR and Global Digital Elevation Models:Applications to Desdyni[C]//Proceedings of 2010 IEEE International Geoscience and Remote Sensing Symposium (IGARSS). Honolulu, HI:IEEE, 2010:1907-1910. |
| [17] | 胡文敏, 邸凯昌, 岳宗玉, 等. 嫦娥一号激光高度计数据交叉点分析与平差处理[J]. 测绘学报, 2013, 42(2): 218–224. HU Wenmin, DI Kaichang, YUE Zongyu, et al. Crossover Analysis and Adjustment for Chang'E-1 Laser Altimeter Data[J]. Acta Geodaetica et Cartographica Sinica, 2013, 42(2): 218–224. |
| [18] | 赵双明, 冉晓雅, 付建红, 等. CE-1立体相机与激光高度计数据联合平差[J]. 测绘学报, 2014, 43(12): 1224–1229. ZHAO Shuangming, RAN Xiaoya, FU Jianhong, et al. Combined Adjustment of CE-1 Stereo Camera Image and Laser Altimeter Data[J]. Acta Geodaetica et Cartographica Sinica, 2014, 43(12): 1224–1229. DOI:10.13485/j.cnki.11-2089.2014.0178 |
| [19] | WANG Fei, VEMURI B C. Non-rigid Multi-modal Image Registration Using Cross-cumulative Residual Entropy[J]. International Journal of Computer Vision, 2007, 74(2): 201–215. DOI:10.1007/s11263-006-0011-2 |
| [20] | 江万寿, 彭芳媛, 岳春宇, 等. 利用Ratio梯度和交叉累积剩余熵进行多源遥感影像匹配[J]. 武汉大学学报 (信息科学版), 2009, 34(9): 1047–1050. JIANG Wanshou, PENG Fangyuan, YUE Chunyu, et al. Multi-source Remote-sensing Image Matching Based on Ratio-gradient and Cross-cumulative Residual Entropy[J]. Geomatics and Information Science of Wuhan University, 2009, 34(9): 1047–1050. |
| [21] | 鄂栋臣, 徐莹, 张小红. ICESat卫星及其在南极Dome A地区的应用[J]. 武汉大学学报 (信息科学版), 2007, 32(12): 1139–1142. E Dongchen, XU Ying, ZHANG Xiaohong. ICESat's Performance and Its Application in Dome A Area in Antarctica[J]. Geomatics and Information Science of Wuhan University, 2007, 32(12): 1139–1142. |
| [22] | 宋志英. 利用ICESAT激光测高数据监测南极冰盖高程变化的研究[D]. 阜新: 辽宁工程技术大学, 2009. SONG Zhiying. The Study of Height Change in Antarctic Ice Sheet by Using ICESAT Laser Altimeter Data[D]. Fuxin:Liaoning Technical University, 2009. |
| [23] | 岳春宇. 面向平坦地区的星载SAR与光学图像配准关键技术研究[D]. 武汉: 武汉大学, 2012. YUE Chunyu. Research on Registration of SAR and Optical Satellite Image in Flat Area[D]. Wuhan:Wuhan University, 2012. |