一种地表起伏区域激光测高点作为影像高程控制点方法

摘要

本发明公开了一种地表起伏区域激光测高点作为影像高程控制点方法，首先对地表起伏区域激光点回波信号开展全波形处理，获得激光点地面足印内的地表最大和最小高程值。随后利用激光足印影像与立体影像的高精度配准，获取立体影像上激光点地面足印范围。接着在立体影像上激光点地面足印范围内开展影像密集匹配并生成地面三维点云，获取高程值最大和最小的地面三维点在立体影像上对应的像点坐标。最后将激光点地面足印内的地表最大和最小高程值分别赋予它们，亦即通过一个激光点获取了立体影像上两个高程控制点。本发明突破了以往地表起伏区域激光点因无法获取足印中心点的精确高程而无法用做高程控制的局限，大幅提升了地表起伏区域激光点的利用率。

权利要求

1.一种地表起伏区域激光测高点作为影像高程控制点方法，其特征在于，所述方法包括以下步骤：

步骤1，对位于地表起伏区域激光测高点的原始回波信号开展全波形处理和分析，获取该激光测高点地面足印范围内多个不同高度地表的高程值，并统计地表的最大高程值Hmax和最小高程值Hmin；

步骤2，将区域内的立体影像构建区域网，在立体影像内部和相邻景立体影像之间布置连接点，开展自由网平差，实现立体影像的高精度相对定向；

步骤3，将激光测高点足印影像与立体影像进行匹配和配准，获取激光测高点地面足印在立体影像上的位置和范围；

步骤4在立体影像上激光测高点地面足印范围内，采用半全局密集匹配算法开展立体影像密集匹配，生成地面三维点云；找出高程值最大和最小的地面三维点，并分别获取高程值最大和最小的地面三维点在立体影像上对应的像点坐标Pmax和Pmin；

步骤5，将激光测高点地面足印范围内的地表最大高程值Hmax，作为利用立体影像密集匹配生成的激光测高点地面足印范围内的高程值最大的地面三维点在立体影像上对应像点坐标Pmax的高程值，亦即获取了1个立体影像上的高程控制点；将激光测高点地面足印范围内的地表最小高程值Hmin，作为利用立体影像密集匹配生成的激光点地面足印范围内的高程值最小的地面三维点在立体影像上对应像点坐标Pmin的高程值，亦即又获得了1个立体影像上的高程控制点；

所述步骤1中，在地表起伏区域，激光测高仪发射的高斯波形的激光脉冲会被地面足印范围内的不同高度的地物多次反射并被星上激光接收单元接收，从而形成多个回波信号；

所述步骤1具体包括：

步骤1.1，对回波波形进行去噪和平滑预处理，并采用多个高斯函数组合的最小二乘拟合分解方法对多个回波信号进行拟合；

步骤1.2，分别计算各个回波波形的拟合高斯形状的波峰位置与激光发射脉冲高斯波峰位置的时间差值Δt，采用公式(1)计算出导致激光脉冲多次回波的激光足印范围内各个地表点的激光测距值ρ：

式中，c为光速；

步骤1.3，根据激光指向角度、卫星姿态和轨道测量数据，可将步骤1.2获取的多个激光测距值ρ，换算成在地面大地坐标系下的多个高程值；统计最大高程值Hmax和最小高程值Hmin；

所述步骤2包括：

步骤2.1将区域内连续覆盖的卫星立体影像构建一个平差区域网；采用影像自动匹配为主，人工判读/量测为辅的方法布设密集的连接点；设置一个区域网平差模型；

步骤2.2逐点构建连接点的误差方程；

步骤2.3对误差方程进行法化，得到法方程；对法方程利用谱修正法迭代估计求解每张影像RFM的仿射变换参数、连接点对应地面点三维坐标改正数；

步骤2.4利用上步获取的改正数更新影像RFM仿射变换参数和连接点对应地面点三维坐标，重复步骤2.2—步骤2.3，通过平差迭代过程不断更新RFM仿射变换参数和连接点对应地面点三维坐标，直至影像RFM放射变换参数中的平移参数a0，b0小于阈值时，平差迭代结束；

步骤2.5平差迭代完成时，统计所有影像连接点的地面坐标通过RFM模型投影至影像上的投影坐标与其原始影像坐标的差值，即为像点精度，也就是平差后所能达到的像方精度；

所述步骤3包括：

步骤3.1，计算足印影像上激光足印中心点位置处像点的尺度不变特征转换特征SIFT，生成SIFT特征描述符；

步骤3.2，从立体影像中选取一张影像，根据激光测高点的地面经纬度坐标和高程值，利用影像的成像几何模型，计算激光测高点地面足印中心点在选定影像上的大致位置，并以此为中心，向四周外扩一定像素范围作为影像匹配的搜索范围；计算选定影像搜索范围内每个像素的SIFT特征，并生成其SIFT特征描述符；

步骤3.3针对选定影像搜索范围内每个像素的SIFT特征描述符，使用快速最近邻逼近搜索函数库搜索与足印影像上激光测高点足印中心点像素的SIFT特征的最邻近值，则该像素为激光测高点地面足印中心位置在选定影像上的像点坐标；

步骤3.4根据激光测高点足印在地面的直径，以及选定影像的空间分辨率，获得激光测高点地面足印在选定影像上的位置和范围信息；

步骤3.5依次从立体影像中选取其他影像，重复执行步骤3.2—步骤3.4，获取激光测高点地面足印在这些影像上的位置和范围信息；

所述步骤4具体包括：

步骤4.1，在获取的立体影像上激光测高点地面足印范围内，采用半全局密集匹配算法生成视差图，得到密集同名点，再利用立体影像成像几何模型进行空间前方交会获取每个同名点的地面三维坐标，生成地面三维点云；

步骤4.2，统计获取高程值最大的地面三维点，由于地面三维点自带在立体影像上的像素坐标信息，亦即获取了激光测高点地面足印范围内，地表最大高程值位置在立体影像上的像点坐标；

步骤4.3，统计获取高程值最小的地面三维点，同样亦即获取了激光测高点地面足印范围内，地表最小高程值位置在立体影像上的像点坐标。

2.如权利要求1所述的地表起伏区域激光测高点作为影像高程控制点方法，其特征在于，所述步骤3.1具体包括：

步骤3.1.1，将足印影像与使用不同“高斯核”的高斯函数进行卷积运算，得到足印影像的高斯多尺度空间空间序列：

L(x,y,σ)＝G(x,y,σ)*I(x,y)                            (2)式中，L(x,y,σ)表示影像的高斯尺度空间；(x,y)为原图像上的像素坐标；σ为尺度空间因子，它是高斯正态分布的标准差，其值越大图像越模糊，对应的尺度也就越大；I(x,y)为原图像上(x,y)位置处的像素值，G(x,y,σ)为可变参数的高斯函数，其形式如下：步骤3.1.2，采用公式(4)将相邻的两个高斯空间的图像相减得到了高斯差分的响应图像，并寻找高斯差分响应图像局部极大或极小值，作为稳定特征点：式中，k为相邻两个高斯尺度空间的比例因子；

步骤3.1.3，针对步骤3.1.2获得的特征点，利用图像的局部特征为特征点分配一个基准方向，使描述符对图像旋转具有不变性；计算以特征点为中心、以3×1.5σ为半径的区域图像的幅角和幅值，每个点L(x,y)的梯度的模m(x,y)以及方向θ(x,y)通过下式求得：再利用图像的梯度直方图求取关键点局部结构的稳定方向；

步骤3.1.4，针对步骤3.1.3获取的SIFT特征点，将坐标轴旋转为特征点的主方向，求取以特征点为中心的16×16的窗口的像素的梯度幅值和梯度方向，将窗口内的像素分成16块，每块上绘制8个方向的梯度直方图并计算每个梯度方向的累加值，可形成128维的特征向量，亦即SIFT特征描述符。

3.如权利要求1所述的地表起伏区域激光测高点作为影像高程控制点方法，其特征在于，所述步骤4.1具体包括：

步骤4.1.1，利用公式(6)的Census变换算法将图像像素的灰度值编码成二进制码流，以此获取邻域像素灰度值相对于中心像素灰度值的大小关系，亦即影像各像素的匹配代价；

式中，p为某一像点位置；q为以像点p为中心的邻域窗口中的其他像素，Np表示以像素p为中心的邻域窗口； 为比特位的逐位连接运算；I(*)表示像点*处的灰度值；ξ运算定义如下：

通过计算左右影像对应的两个像素的Census变换值的Hamming距离，获取该对像素的匹配代价：

C(p,d)＝Hamming(Csl(p),Csr(d))              (8)式中，d为其他影像上与像素p对应的像素；

步骤4.1.2，采用动态规划的方法，合并至少8个方向的一维动态规划结果实现二维的代价积聚；

利用公式(9)计算像素p沿某一路径r的匹配代价聚合值：式中，d为像素p在立体影像上的左右视差值；Lr(p,d)为像素p路径r上在视差d处的代价积聚值；C(p,d)为像点p在视差为d时的匹配代价；右边等式第2项表示像素p在r方向上的前一个像素代价累计的最小值；右边等式第3项是为了防止数值过大，对最优路径的选择没有影响；P1和P2分别是当前像素p和路径r方向上前一个像素的视差差异为1的惩罚系数和视差差异大于1的惩罚系数，且P2>P1；

将像素p所有方向的匹配代价聚合值相加得到像素的最终的匹配代价，亦即聚合代价立方体：

步骤4.1.3，对于每一个像素，从上式(10)选择聚合代价最小的视差作为最优视差，从而生成最终的视差图D：

D＝mindS(p,d)                                  (11)步骤4.1.4，为了获得更加精细平滑的视差图，对计算得到的视差图开展进一步的优化，包括利用左右一致性检查剔除粗差，使用二次曲线拟合相邻匹配代价进行亚像素插值，基于快速双边算子进行平滑处理；

步骤4.1.5，根据视差图，获取立体影像同名点，通过立体影像的空间前方交会得到密集的地面三维点云数据，并进行必要的粗差剔除。

说明书

一种地表起伏区域激光测高点作为影像高程控制点方法

技术领域

[0001] 本发明属于遥感与摄影测量技术领域，尤其涉及一种地表起伏区域激光测高点作为立体影像高程控制点方法。

背景技术

[0002] 现有激光测高仪作为一种高精度测距仪器，因具有方向性好、测距精度高等特点，可以获取高精度地表高程信息。目前可用的星载激光测高数据主要包括两类，一类是专门
的激光测高卫星获取的激光测高点数据，如xx国于2003年发射的ICESat卫星上搭载的
Geoscience Laser Altimeter System获取的激光点测高点数据；xx国于2018年发射的
ICESat‑02卫星上搭载的Advanced Topographic Laser Altimeter System获取的激光测
高点数据。另一类是立体测图卫星同时搭载了光学立体相机和激光测高仪，可同步获取光
学立体影像和激光测高点数据，如我国于2019年11月发射的高分七号(GF‑7)卫星上同时搭
载了两线阵光学立体相机和一台2波束的激光测高仪，可同步获取亚米分辨率的两线阵光
学立体影像和高精度的激光测高点。我国于2020年7月发射的资源三号03星同时搭载了三
线阵光学立体相机和一台单波束的激光测高仪，可同步获取2米分辨率的三线阵立体影像
和高精度的激光测高点。

[0003] 目前国内外很多学者都将激光测高点作为立体影像的高程控制点，用于提升卫星立体影像的高程精度。将激光测高点作为立体影像的高程控制点的一个重要前提是需要获
取激光测高点在立体影像上的准确像方坐标。此前的所有研究中，一般都是将激光测高点
地面足印中心点作为立体影像的高程控制点，然而由于仅分布于地表平滑(即地形平坦且
地物单一)区域的激光测高点才可以获取地面足印中心点位置的高精度高程值，因此以往
只能挑选地面足印范围内地表平滑的激光测高点作为立体影像控制，而舍弃了足印范围内
地表起伏(即地形起伏或存在不同高程地物)的激光测高点的使用，在激光测高点的使用上
存在较大局限性，大量地面足印范围内地表起伏的激光测高点被迫舍弃了，造成了巨大的
资源浪费。而地形起伏区域往往也是其他类型高精度地面控制点相对匮乏区域，这也降低
了激光测高点作为立体影像高程控制点的适用性和实用性。综上所述，目前技术存在的问
题是：利用激光测高点作为高程控制的常规方法是将激光测高点地面足印中心点作为立体
影像的高程控制点，导致分布于地表起伏(即地形起伏、地物复杂)区域的激光测高点由于
无法获取地面足印中心点位置的高精度高程值，而无法应用于立体影像高程控制，造成大
量激光测高点闲置浪费，导致激光测高点在地表复杂区域的利用效率极为低下。

[0004] 因此，设计一种能够充分利用地表起伏区域激光测高点作为立体影像高程控制点的有效方法，是破解目前地表起伏区域激光测高点无法有效应用于立体影像高程控制这一
难题，提升全区域激光测高点利用率的迫切需求，也必将创造较大的社会和经济效益。

发明内容

[0005] 为解决上述技术问题，本发明充分利用激光测高点获取的地面足印范围内不同高度地物反射的多个回波信号，提出了一种全新的利用激光测高点作为立体影像高程控制的
方法，有效提升激光测高点在地表起伏区域的利用率。

[0006] 本发明的目的通过以下的技术方案来实现：

[0007] 一种地表起伏区域激光测高点作为影像高程控制点方法，包括以下步骤：

[0008] 步骤1，对位于地表起伏区域激光测高点的原始回波信号开展全波形处理和分析，获取该激光测高点地面足印范围内多个不同高度地表的高程值，并统计地表的最大高程值
Hmax和最小高程值Hmin；

[0009] 步骤2，将区域内的立体影像构建区域网，在立体影像内部和相邻景立体影像之间布置连接点，开展自由网平差，实现立体影像的高精度相对定向；

[0010] 步骤3，将激光测高点足印影像与立体影像进行匹配和配准，获取激光测高点地面足印在立体影像上的位置和范围；

[0011] 步骤4在立体影像上激光测高点地面足印范围内，采用半全局密集匹配算法开展立体影像密集匹配，生成地面三维点云；找出高程值最大和最小的地面三维点，并分别获取
它们在立体影像上对应的像点坐标Pmax和Pmin；

[0012] 步骤5，将激光测高点地面足印范围内的地表最大高程值Hmax，作为利用立体影像密集匹配生成的激光测高点地面足印范围内的高程值最大的地面三维点在立体影像上对
应像点坐标Pmax的高程值，亦即获取了1个立体影像上的高程控制点；将激光测高点地面足
印范围内的地表最小高程值Hmin，作为利用立体影像密集匹配生成的激光点地面足印范围
内的高程值最小的地面三维点在立体影像上对应像点坐标Pmin的高程值，亦即又获得了1个
立体影像上的高程控制点。

[0013] 与现有技术相比，本发明的一个或多个实施例可以具有如下优点：

[0014] (1)突破了以往激光地面足印范围内地表起伏的激光测高点无法用于立体影像高程控制的瓶颈，解决了以往山地、高山地地形区域可用激光测高点数量极少的窘境，大幅提
升了激光测高点在全地形(尤其是山地和高山地)区域的可用性和有效性。

[0015] (2)充分利用激光测高点可以较为准确地获取地面足印内地表最大和最小高程值的特性，成功利用一个激光测高点获取了两个用于立体影像的高程控制点，相比传统的一
个激光点仅能获取一个高程控制点的常规方法，增加了控制点数量，将激光测高点的利用
效率提升了一倍。

[0016] (3)提取的高程控制点在立体影像上像素位置较为准确，加之激光测高点高程值精度极高，能有效提升立体影像的高程精度，可媲美传统的采用高精度地面控制点的影像
处理精度；尤其适用于无法获取其它高精度地面控制资料区域(如境外区域、境内山地等控
制获取困难区域等)的卫星影像高精度测图应用，将产生了极大的经济效益和社会效益。

附图说明

[0017] 图1是地表起伏区域激光测高点作为影像高程控制点方法流程图。

具体实施方式

[0018] 为使本发明的目的、技术方案和优点更加清楚，下面将结合实施例及附图对本发明作进一步详细的描述。

[0019] 如图1所示，为本发明提供的地表起伏区域激光测高点作为影像高程控制点方法，包括以下步骤：

[0020] 步骤1，对位于地表起伏区域激光测高点的原始回波信号开展全波形处理和分析，获取该激光测高点地面足印范围内多个不同高度地表的高程值，并统计地表最大高程值
Hmax和最小高程值Hmin；

[0021] 在地表起伏区域，激光测高仪发射的高斯波形的激光脉冲会被地面足印范围内的不同高度的地物多次反射并被星上激光接收单元接收，从而形成多个回波信号。本步骤进
一步包括：

[0022] 步骤1.1，对回波波形进行去噪和平滑等预处理，并采用多个高斯函数组合的最小二乘拟合分解方法对多个回波信进行拟合；

[0023] 步骤1.2，分别计算各个回波波形的拟合高斯形状的波峰位置与激光器发射脉冲高斯波峰位置的时间差值Δt，采用公式(1)计算出导致激光脉冲多次回波的激光足印范围
内各个地表点的激光测距值ρ：

[0024]

[0025] 式中，c为光速。

[0026] 步骤1.3，根据激光指向角度、卫星姿态和轨道测量数据，可将步骤1.2获取的多个激光测距值ρ，换算成在地面大地坐标系下的多个高程值；统计最大高程值Hmax和最小高程
值Hmin。

[0027] 步骤2，将区域内的立体影像构建区域网，在立体影像内部和相邻景立体影像之间布设密集的连接点，开展自由网平差，实现立体影像的高精度相对定向，为后续其他操作提
供基础。

[0028] 本步骤进一步包括：

[0029] 步骤2.1，将区域内连续覆盖的卫星立体影像构建一个平差区域网；采用影像自动匹配为主，人工判读/量测为辅的方法布设密集的连接点；设置一个区域网平差模型。

[0030] 本实例采用像方仿射变换补偿的基于有理函数模型(RFM)的区域网平差模型，亦即构建一个影像像方坐标点的仿射变换模型，对影像RFM计算获得的像点坐标进行误差改
正，其基本方程定义如下：

[0031]

[0032] 式中，为避免计算过程中参数数值量级差别过大而引入舍入误差，需要将影像像点坐标和对应的地面点坐标均归一化到‑1～1之间，以增强参数求解的稳定性。(X,Y,Z)和
(r,c)分别是归一化后的地面点坐标和影像坐标。Pi(i＝1,2,3,4)分别表示一般多项式，式
中各变量Xn,Yn,Zn的幂均不超过3次，所有变量的幂之和也不超过3次。(Δr,Δc)为像点坐
标(r,c)的系统误差像方补偿值，其值为：

[0033]

[0034] 式中，(a0,a1,a2,b0,b1,b2)表示补偿影像RFM系统误差的仿射变换参数。

[0035] 步骤2.2，逐点构建连接点的误差方程如下：

[0036] V＝At+Bx‑L P                        (4)

[0037] 式中，将影像RFM的仿射变换参数和连接点对应的地面三维坐标作为未知数；V＝T
[vR vC] 为连接点在影像上行、列坐标观测值残差向量；t＝[Δa0 Δa1 Δa2 Δb0 Δb1 Δ
T
b2]为影像RFM的像方坐标系统误差仿射变换补偿参数的改正数向量；x＝[ΔXtie ΔYtie Δ
T
Ztie]为连接点对应地面大地坐标的改正数向量； 是未知数t对应
的偏导数系数矩阵； 是未知数x对应的偏导数系数矩阵； 为
0 0
利用初值代入观测方程计算得到的常数项，(r ,c)为利用未知数的近似值通过RFM计算出
的连接点像点行、列坐标；P为各观测值的权矩阵，通常各类观测值的权值由他们的先验信
息确定，一般可取观测值10倍先验标准差，并在每次平差迭代计算后重新计算各类观测值
的权。

[0038] 步骤2.3，对误差方程进行法化，得到法方程；对法方程利用谱修正法迭代估计求解每张影像RFM的仿射变换参数、连接点对应地面点三维坐标改正数。

[0039] 步骤2.4：利用上步获取的改正数更新影像RFM仿射变换参数和连接点对应地面点三维坐标，重复步骤2.2—步骤2.3，通过平差迭代过程不断更新RFM仿射变换参数和连接点
对应地面点三维坐标，直至影像RFM放射变换参数中的平移参数小于阈值时(本实施例为
0.1个像元pixel，但是不限于此)时，平差迭代结束；如果迭代次数达到预设迭代次数(本实
施例设定的迭代次数为20次，但是不限于此)，仍然不能收敛，那么平差失败退出，此时的平
差精度可能有损失。

[0040] 步骤2.5：平差迭代完成时，统计所有影像连接点的地面坐标通过RFM模型投影至影像上的投影坐标与其原始影像坐标的差值，即为像点精度，也就是平差后所能达到的像
方精度。当像方精度满足要求后(本实例中像方精度优于1个像素)，利用计算得到的仿射变
换参数对影像RFM进行像方仿射变换补偿，生成得到新的RFM文件，此时立体影像已经实现
高精度相对定向。

[0041] 步骤3，将激光测高点足印影像与立体影像开展高精度匹配和配准，获取激光测高点地面足印在立体影像上的位置和范围。

[0042] 本步骤进一步包括：

[0043] 步骤3.1，计算足印影像上激光足印中心点位置处像点的尺度不变特征转换(Scale‑invariant feature transform,SIFT)特征，生成SIFT特征描述符。具体包括：

[0044] 步骤3.1.1，将足印影像与使用不同“高斯核”的高斯函数进行卷积运算，得到足印影像的高斯多尺度空间空间序列：

[0045] L(x,y,σ)＝G(x,y,σ)*I(x,y)                            (5)

[0046] 式中，L(x,y,σ)表示影像的高斯尺度空间；(x,y)为原图像上的像素坐标；σ为尺度空间因子，它是高斯正态分布的标准差，其值越大图像越模糊，对应的尺度也就越大；I(x,
y)为原图像上(x,y)位置处的像素值，G(x,y,σ)为可变参数的高斯函数，其形式如下：

[0047]

[0048] 步骤3.1.2，采用公式(7)将相邻的两个高斯空间的图像相减得到了高斯差分的响应图像，并寻找高斯差分响应图像局部极大或极小值，作为稳定特征点：

[0049]

[0050] 式中，k为相邻两个高斯尺度空间的比例因子。

[0051] 步骤3.1.3，针对步骤3.1.2获得的特征点，利用图像的局部特征为特征点分配一个基准方向，使描述符对图像旋转具有不变性。计算以特征点为中心、以3×1.5σ为半径的
区域图像的幅角和幅值，每个点L(x,y)的梯度的模m(x,y)以及方向θ(x,y)通过下式求得：

[0052]

[0053] 再利用图像的梯度直方图求取关键点局部结构的稳定方向。

[0054] 步骤3.1.4，针对步骤3.1.3获取的SIFT特征点，将坐标轴旋转为特征点的主方向，求取以特征点为中心的16×16的窗口的像素的梯度幅值和梯度方向，将窗口内的像素分成
16块，每块上绘制8个方向的梯度直方图并计算每个梯度方向的累加值，可形成128维的特
征向量，亦即SIFT特征描述符。

[0055] 步骤3.2，从立体影像中选取一张影像，根据激光测高点的地面经纬度坐标和高程值，利用影像的成像几何模型，计算激光测高点地面足印中心点在选定影像上的大致位置，
并以此为中心，向四周外扩一定像素范围作为影像匹配的搜索范围；计算选定影像搜索范
围内每个像素的SIFT特征，并生成其SIFT特征描述符。

[0056] 步骤3.3，针对选定影像搜索范围内每个像素的SIFT特征描述符，使用快速最近邻逼近搜索函数库搜索与足印影像上激光测高点足印中心点像素的SIFT特征的最邻近值，则
该像素为激光测高点地面足印中心位置在选定影像上的像点坐标。

[0057] 步骤3.4，根据激光测高点足印在地面的直径，以及选定影像的空间分辨率，获得激光测高点地面足印在选定影像上的位置和范围等信息。

[0058] 步骤3.5，依次从立体影像中选取其他影像，复步骤3.2—步骤3.4，获取激光测高点地面足印在这些影像上的位置和范围信息。

[0059] 步骤4，在立体影像上激光测高点地面足印范围内，采用半全局密集匹配算法开展立体影像密集匹配，生成地面三维点云；找出高程值最大和最小的地面三维点，并分别获取
它们在立体影像上对应的像点坐标Pmax和Pmin。

[0060] 本步骤进一步包括：

[0061] 步骤4.1，在获取的立体影像上激光测高点地面足印范围内，采用半全局密集匹配算法生成视差图，得到密集同名点，再利用立体影像成像几何模型进行空间前方交会获取
每个同名点的地面三维坐标，生成地面三维点云。具体包括：

[0062] 步骤4.1.1，利用公式(9)Census变换算法将图像像素的灰度值编码成二进制码流，以此获取邻域像素灰度值相对于中心像素灰度值的大小关系，亦即影像各像素的匹配
代价；

[0063]

[0064] 式中，p为某一像点位置；q为以像点p为中心的邻域窗口中的其他像素，Np表示以像素p为中心的邻域窗口； 为比特位的逐位连接运算；I(*)表示像点*处的灰度值；ξ运算
定义如下：

[0065]

[0066] 通过计算左右影像对应的两个像素的Census变换值的汉明(Hamming)距离，获取该对像素的匹配代价：

[0067] C(p,d)＝Hamming(Csl(p),Csr(d))              (11)

[0068] 式中，d为其他影像上与像素p对应的像素。

[0069] 步骤4.1.2，采用动态规划的方法，合并至少8个方向的一维动态规划结果实现二维的代价积聚。

[0070] 利用公式(12)计算像素p沿某一路径r的匹配代价聚合值：

[0071]

[0072] 式中，d为像素p在立体影像上的左右视差值；Lr(p,d)为像素p路径r上在视差d处的代价积聚值；C(p,d)为像点p在视差为d时的匹配代价；右边等式第2项表示像素p在r方向
上的前一个像素代价累计的最小值；右边等式第3项是为了防止数值过大，对最优路径的选
择没有影响；P1和P2分别是当前像素p和路径r方向上前一个像素的视差差异为1的惩罚系数
和视差差异大于1的惩罚系数，且P2>P1。

[0073] 将像素p所有方向的匹配代价聚合值相加得到像素的最终的匹配代价，亦即聚合代价立方体：

[0074]

[0075] 步骤4.1.3，对于每一个像素，从上式(13)选择聚合代价最小的视差作为最优视差，从而生成最终的视差图D：

[0076] D＝mindS(p,d)                                  (14)

[0077] 步骤4.1.4，为了获得更加精细平滑的视差图，对计算得到的视差图开展进一步的优化，包括利用左右一致性检查剔除粗差，使用二次曲线拟合相邻匹配代价进行亚像素插
值，基于快速双边算子进行平滑处理等。

[0078] 步骤4.1.5，根据视差图，获取立体影像同名点，通过立体影像的空间前方交会得到密集的地面三维点云数据，并进行必要的粗差剔除。

[0079] 步骤4.2，统计获取高程值最大的地面三维点，由于地面三维点自带在立体影像上的像素坐标信息，亦即获取了激光测高点地面足印范围内，地表最大高程值位置在立体影
像上的像点坐标。

[0080] 步骤4.3，统计获取高程值最小的地面三维点，同样亦即获取了激光测高点地面足印范围内，地表最小高程值位置在立体影像上的像点坐标。

[0081] 步骤5，将激光测高点地面足印范围内的地表最大高程值Hmax，作为利用立体影像密集匹配生成的激光测高点地面足印范围内的高程值最大的地面三维点在立体影像上对
应像点坐标Pmax的高程值，亦即获取了1个立体影像上的高程控制点；将激光测高点地面足
印范围内的地表最小高程值Hmin，作为利用立体影像密集匹配生成的激光点地面足印范围
内的高程值最小额地面三维点在立体影像上对应像点坐标Pmin的高程值，亦即又获得了1个
立体影像上的高程控制点。

[0082] 虽然本发明所揭露的实施方式如上，但所述的内容只是为了便于理解本发明而采用的实施方式，并非用以限定本发明。任何本发明所属技术领域内的技术人员，在不脱离本
发明所揭露的精神和范围的前提下，可以在实施的形式上及细节上作任何的修改与变化，
但本发明的专利保护范围，仍须以所附的权利要求书所界定的范围为准。