一种基于光子点云辅助的卫星摄影测量对地定位方法

1.本发明属于卫星摄影测量与遥感领域，具体涉及一种基于光子点云辅助的卫星摄影测量对地定位方法。


背景技术：

2.随着光学遥感卫星和数字摄影测量技术的发展，卫星摄影测量已成为一项被广泛应用于国防经济建设领域的遥感测绘技术。现今光学遥感卫星如资源三号、高分七号、wordview3/4、pleiades和cartosat-3等，其地面分辨率均已全面迈入分米级水平。各类高分辨率遥感卫星陆续发射运行，为各行业领域提供数量丰富、质量优良的影像数据。同时，国内外学者对利用遥感影像进行地形测绘的理论方法进行了大量研究，形成了较为成熟的卫星摄影测量技术体系。
3.近些年，新兴的卫星激光测高技术上取得重大进展，对地观测定位能力大幅增强，为地形测绘等领域提供了强力的数据支撑。icesat-2卫星采用先进地形激光测高系统(advanced topographic laser altimeter，atlas)，具备强大的激光探测能力，能够采集观测频率更高、定位精度更高的地表激光点云数据。
4.由于光学卫星在成像的过程中，受到姿轨观测误差、平台颤振、大气折光和内部拼接等各种因素影响，其定位结果中含有复杂的系统误差，需要采集一定数量均匀分布的地面控制点，通过区域网平差提高定位精度。但对于如山林、沙漠和高原等区域，获取足够的高精度地面控制点非常困难，作业成本高昂。而卫星激光测高技术可以不受地形地貌的影响，在轨自动采集高精度的测高点数据。因此对icesat-2光子点云和光学遥感影像联合处理，利用激光测高点代替传统人工外业测量点，可大幅提高卫星摄影测量技术提取地表点云的效率。对于星载激光雷达数据与卫星遥感影像的联合处理，已有的研究表明，融合激光测高数据可以显著提高卫星摄影测量技术的定位精度。当前常用的联合处理方案是先从点云中筛选出位于平坦区域的激光点作为控制点，然后在像方坐标中，将控制点作为带权观测值进行区域网平差。该方案可有效提高无外业控制点下的卫星摄影测量定位精度，但面临以下待进一步解决的问题：
5.(1)与传统脉冲式激光雷达不同，当前的激光测高卫星采用新型的光子计数激光雷达，其观测数据为点云形式。由于激光雷达对植被具有一定的穿透性，山林地区光子点云中含有来自树冠反射和地表反射的辐射信息，而光学遥感影像所见主要为树冠目标。
6.(2)建立激光点与影像的准确对应关系是进行两种数据联合处理的前提条件，对此目前通常采用的是像方配准方案，即先选取平坦地区的激光点，利用原始的或经过自由网平差处理后的影像定位模型，计算激光点的像方投影点作为影像对应点，由于定位模型不可避免地存在误差，导致激光点物像对应关系不准确，不适用于对于地势崎岖的观测区域。
7.综上所述，由于光子点云中包含地表反射的辐射信息以及像方配准方案中影像定位模型存在误差，导致现有的联合处理方法的定位精度低。因此对卫星影像数据和光子雷
达数据进行联合处理时，需要从光子点云中分类提取出树冠高程，以保证目标高程的一致性。那么研究光子点云滤波分类算法以提取准确地表高程点，对于联合处理方案具有重要意义。并且鉴于现有联合处理方案的不足，有必要深入研究激光点与遥感影像的物方配准方案，为卫星影像数据和光子雷达数据的联合处理提供必要条件。


技术实现要素：

8.本发明的目的是为解决现有的联合处理方法的定位精度低的问题，而提出的一种基于光子点云辅助的卫星摄影测量对地定位方法。
9.本发明为解决上述技术问题所采取的技术方案是：
10.一种基于光子点云辅助的卫星摄影测量对地定位方法，所述方法具体包括以下步骤：
11.步骤一、获取目标区域的立体影像数据；
12.步骤二、利用立体影像数据生成地物的摄影测量点云；
13.步骤三、获取立体像对区域内的icesat-2光子点云；
14.步骤四、对步骤三中获取的光子点云进行栅格滤波，滤除噪声点后，得到滤波后的光子点云；
15.步骤五、根据滤波后的光子点云等间隔生成高程采样点，再分别计算出每个高程采样点的地面高程值；
16.步骤六、根据各个高程采样点的地面高程值，将高程采样点云和步骤二生成的摄影测量点云进行配准，获得变换矩阵；再根据变换矩阵获得作为控制点的高程采样点；
17.步骤七、根据获得的控制点对摄影测量点云的系统误差进行修正。
18.本发明的有益效果是：
19.1、本发明提出基于栅格连续性约束的光子点云滤波算法，光子点云滤波及高程提取效果优良。对icesat-2/atl03数据的试验结果表明，该算法能高效滤除噪声点并提取地表高程点，滤波成功率高于94％，高程提取精度优于5.4m，对不同地形变化、地物类型和激光能量的光子点云数据均可取得稳定处理效果。
20.2、本发明提出的基于线面状点云配准的联合处理算法，对摄影测量点云定位精度提高显著。试验结果表明引入icesat-2光子点云辅助后，资源三号摄影测点云的x、y、z方向精度分别达到0.46m、2.43m和2.13m，高分七号摄影测点云的x、y、z方向精度分别达到0.36m、1.34m和0.59m，定位精度得到显著提高。
附图说明
21.图1为本发明方法的流程图；
22.图2为icesat-2光子点云滤波及高程提取流程图；
23.图3为栅格连续性判断示意图；
24.图4为点云降趋势处理示意图；
25.图5为高程采样点示意图；
26.图6a为获取高程采样点的邻近点的示意图；
27.图6b为邻域点云高程频数统计示意图；
28.图7为光子点云与摄影测量点云配准流程图；
29.图8为高程相似度计算示意图；
30.图9a为资源三号摄影测量点云与icesat-2点云的示意图；
31.图9b为资源三号摄影测量点云系统误差修正前的示意图；
32.图9c为资源三号摄影测量点云系统误差修正后的示意图；
33.图10a为高分七号摄影测量点云与icesat-2点云的示意图；
34.图10b为高分七号摄影测量点云系统误差修正前的示意图；
35.图10c为高分七号摄影测量点云系统误差修正后的示意图。
具体实施方式
36.具体实施方式一、结合图1说明本实施方式。本实施方式所述的一种基于光子点云辅助的卫星摄影测量对地定位方法，所述方法具体包括以下步骤：
37.步骤一、获取目标区域的立体影像数据；
38.步骤二、利用立体影像数据生成地物的摄影测量点云；
39.步骤三、获取立体像对区域(也即目标区域)内的icesat-2光子点云；
40.步骤四、对步骤三中获取的光子点云进行栅格滤波，滤除噪声点后，得到滤波后的光子点云；
41.步骤五、根据滤波后的光子点云等间隔生成高程采样点，再分别计算出每个高程采样点的地面高程值；
42.步骤六、根据各个高程采样点的地面高程值，将高程采样点云和步骤二生成的摄影测量点云进行配准，获得变换矩阵；再根据变换矩阵获得作为控制点的高程采样点；
43.步骤七、根据获得的控制点对摄影测量点云的系统误差进行修正。
44.具体实施方式二：结合图2说明本实施方式。本实施方式与具体实施方式一不同的是，所述步骤四的具体过程为：
45.步骤四一、栅格划分
46.对整个光子点云剖面范围进行栅格划分，形成大小为w
×
h的格网；
47.步骤四二、分别统计格网中每个栅格内包含的点数量，对于任意的一列栅格，从该列的全部栅格中选取出包含点数量最多的前t个栅格作为该列的候选栅格：
[0048][0049]
其中，gj代表提取出的第j列中的候选栅格，为第j列中第m个栅格内包含的点数量，m为该列中的栅格总个数，代表从第j列中选取出包含点数量最多的t个栅格；
[0050]
同理，分别获得每列的候选栅格；
[0051]
步骤四三、栅格连续性检测
[0052]
步骤四三一、将第j列中的某个候选栅格作为待检测栅格，将第j-1列、第j-2列、
…
、第j-k’列中的候选栅格以及第j+1列、第j+2列、
…
、第j+k’列中的候选栅格作为待检测栅格的邻域范围(即选取待检测栅格左、右各相邻k
′
列栅格中的候选栅格作为邻域范围，且选取的左侧列数与右侧列数可以不相等)，待检测栅格与第j-1列、第j-2列、
…
、第j-k’列
中的候选栅格形成左连通域，待检测栅格与第j+1列、第j+2列、
…
、第j+k’列中的候选栅格形成右连通域；如图3所示。
[0053][0054][0055]
式中，为左连通域内的栅格集合，为右连通域内的栅格集合，(i,j)为待检测栅格在格网中的行列序号，为左连通域内的栅格集合中的第m1个栅格，为右连通域内的栅格集合中的第n1个栅格；
[0056]
步骤四三二、计算待检测栅格的连续性得分continue_score：
[0057][0058]
式中，为左连通域内的栅格数量，为右连通域内的栅格数量；
[0059]
步骤四三三、遍历计算出格网的第j列中每个候选栅格的连续性得分，将连续性得分最高的候选栅格以及连续性得分最高候选栅格的上、下各相邻的t'个栅格作为第j列中的信号栅格；
[0060]
步骤四三四、重复步骤四三一至步骤四三三三的过程，分别获得格网中每列的信号栅格；
[0061]
缩小单位栅格的宽度和高度，再基于缩小后单位栅格的宽度和高度对获得的信号栅格点云剖面范围进行栅格划分；
[0062]
步骤四四、根据步骤四三四获得的栅格划分结果，重复执行步骤四二和步骤四三的过程；直至缩小后单位栅格的宽度和高度均小于给定的最小值时停止，最后一次迭代所获得的信号栅格中包含的点作为滤波后的点云。
[0063]
其它步骤及参数与具体实施方式一相同。
[0064]
具体实施方式三：本实施方式与具体实施方式一或二不同的是，所述t的取值为3。
[0065]
本实施方式中，考虑到信号点在高程方向上分布的连续性及容错性要求，将t的取值设置为3。
[0066]
其它步骤及参数与具体实施方式一或二相同。
[0067]
具体实施方式四：本实施方式与具体实施方式一至三之一不同的是，所述t'的取值为1。
[0068]
本实施方式中，考虑到陡坡地形以及栅格边缘的截断效应，可保留连续性得分最高候选栅格的上、下各相邻的1个栅格。
[0069]
其它步骤及参数与具体实施方式一至三之一相同。
[0070]
具体实施方式五：本实施方式与具体实施方式一至四之一不同的是，所述缩小单位栅格的宽度和高度的具体方法为：
[0071]
width
i’+1
＝width
i’/rwꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)
[0072]
height
i’+1
＝height
i’/rhꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(6)
[0073]
式中，width
i’为第i’次滤波的单位栅格宽度，height
i’为第i’次滤波的单位栅格
高度，rw为单位栅格宽度的缩放比例系数，rh为单位栅格高度的缩放比例系数。
[0074]rw
和rh一般取2即可。
[0075]
其它步骤及参数与具体实施方式一至四之一相同。
[0076]
具体实施方式六：本实施方式与具体实施方式一至五之一不同的是，所述步骤五的具体过程为：
[0077]
步骤五一、点云降趋势处理
[0078]
对滤波后的光子点云进行栅格划分(栅格宽高的取值范围一般在15米以内即可)，形成大小为w0
×
h0的格网；
[0079]
对于格网中的第j列栅格，计算该列中所有点的高程值的中位数，将中位数作为该列的趋势常数项；再分别利用该列中每个点的高程值减去该列的趋势常数项，得到该列中每个点对应的降趋势后的高程值，如图4所示：
[0080][0081][0082]
其中，c_detrendj为第j列的趋势常数项，gridj为w0
×
h0的格网中第j列内包含的点的高程值的集合，代表格网中第j列栅格内的第i”个点的高程值，i”＝1,2，
…
，i，i为第j列栅格内的点的总个数，为第j列栅格内的第i”个点降趋势后的高程值；
[0083]
同理，分别获得格网中每个点降趋势后的高程值；
[0084]
步骤五二、生成等间隔高程采样点
[0085]
在降趋势后点云的沿轨距离方向上，根据降趋势后点云沿轨距离的最小值和最大值，等间隔内插形成高程采样点，如图5所示；
[0086][0087][0088]
式中，为降趋势后点云沿轨距离的最大值，为降趋势后点云沿轨距离的最小值，d
sample
为相邻的两个高程采样点之间的间距，n
sample
为高程采样点的总数量，为第k个高程采样点的沿轨距离；
[0089]
为了使提取的高程曲线尽可能符合实际地形，d
sample
的取值不宜超过10米。
[0090]
步骤五三、获取邻近点
[0091]
如图6a所示，对于第k个高程采样点，按给定的邻域距离，选取第k个高程采样点的邻近点集合：
[0092][0093]
其中，为第k个高程采样点的邻近点集合，表示集合内的第n个点，为邻近点集合内的点的沿轨距离，d
ε
为邻域距离；
[0094]
步骤五四、统计邻近点集合的频数分布
[0095]
根据集合内降趋势后点云高程的最大值、最小值，并设置组距，统计邻近点集合的频数分布直方图；
[0096]
步骤五五、计算点云分类阈值
[0097][0098]
其中，binj′
为频数分布直方图中第j
′
个区间的频数，n
bin
为频数分布直方图中的区间数，ρ为阈值系数，表示第k个高程采样点的点云分类阈值；
[0099]
步骤五六、点云分类
[0100]
在邻近点集合的频数分布直方图中，按高程值由大到小的变化方向，以第一个频数值大于分类阈值的区间作为地物高程区间，将落在地物高程区间内的点云作为地物点云，如图6b所示；
[0101]
步骤五七、计算高程采样点的地面高程
[0102]
计算落在地物高程区间内的点云的降趋势后高程值的平均值，并将计算出的平均值作为高程采样点的地物高程值；
[0103][0104][0105]
式中，为第k个高程采样点的地物高程值，bin
surface
为地物高程区间，n
surface
为地物高程区间内的点云数量，为落在地物高程区间内的第l个点云的降趋势后高程值，c_detrend
l
为第l个点云所对应的趋势常数项；
[0106]
步骤五八、采用步骤五三至步骤五七的方法，分别提取出每个高程采样点的地物高程值。
[0107]
其它步骤及参数与具体实施方式一至五之一相同。
[0108]
具体实施方式七：结合图7说明本实施方式。本实施方式与具体实施方式一至六之一不同的是，所述步骤六的具体过程为：
[0109]
步骤六一、以固定的沿轨距离长度对高程采样点云进行分割，获得各段高程采样点云；
[0110]
步骤六二、对于任意的一段高程采样点云，在点云xoy平面(平面直角坐标系)上，分别设置了代表x方向坐标偏移量的参数dx以及代表y方向坐标偏移量的参数dy，按照x方向搜索范围l
dx
、y方向搜索范围l
dy
和x方向搜索步长α
dx
、y方向搜索步长α
dy
，生成参数dx和dy所有取值的全部组合θ{(dx1,dy1),9dx2,dy2),
…
,(dxj″
,dyj″
),
…
}；
[0111]
[0112][0113]
步骤六三、高程值内插
[0114]
根据坐标偏移量(dx1,dy1)对该段高程采样点云进行整体平移后，以平移后的高程采样点pa的平面坐标为圆心，将落在圆心的半径r范围内的摄影测量点(q1,q2,
…
,qm′
)作为邻域点，再利用邻域点(q1,q2,
…
,qm′
)进行反距离插值，获取内插点的高程，插值的公式如下：
[0115][0116][0117][0118]
式中，为平移后的高程采样点pa的平面坐标，为摄影测量点qb的三维坐标，为pa、qb两点间的平面距离，为邻域点集合，wb为反距离插值的权重，为高程的内插值；
[0119]
步骤六四、高程相似度计算，如图8所示：
[0120][0121][0122][0123]
式中，为平移后的高程采样点pa的高程，n
′
为该段高程采样点云中采样点的个数，ρ为高程相似度；
[0124]
步骤六五、重复步骤六三至步骤六四的过程，使坐标偏移量遍历组合θ中的每组坐标偏移量，取最大的高程相似度ρ
max
对应的坐标偏移量参数组合，作为当前搜索范围和搜索步长下的最优估计值
[0125]
步骤六六、对于该段高程采样点云，不断的按比例系数调整搜索范围和搜索步长，且每次调整后均重复步骤六二至步骤六五的过程，在每组搜索范围和搜索步长下，均获得对应的最优估计值以及对应的最大高程相似度；
[0126]
直至时或达到设置的最大迭代次数时终止迭代，将最后一次迭代所获得的最优估计值作为最终的x方向坐标偏移量和y方向坐标偏移量；
[0127]
其中，为第i次迭代获得的最大高程相似度，为第i+1次迭代获得的最大高程相似度；
[0128]
步骤六七、根据最终的x方向坐标偏移量和y方向坐标偏移量对该段高程采样点云进行平移；
[0129]
步骤六八、同理，采用步骤六二至步骤六七的方法，分别对每段高程采样点云进行
平移；
[0130]
步骤六九、从各段平移后的高程采样点云中随机抽样出若干段高程采样点云，利用icp算法对抽样出的若干段高程采样点云中的高程采样点与摄影测量点云进行配准，得到高程采样点与摄影测量点云间误差距离最小的变换矩阵；
[0131]
利用变换矩阵对全部的高程采样点进行刚性变换后，剔除掉误差距离不满足阈值要求的高程采样点，将剩余的高程采样点作为控制点。
[0132]
其它步骤及参数与具体实施方式一至六之一相同。
[0133]
具体实施方式八：本实施方式与具体实施方式一至七之一不同的是，所述按比例系数调整搜索范围和搜索步长，其具体过程为：
[0134][0135][0136][0137][0138]
式中，为搜索范围的比例系数，一般取值均为2，为搜索步长的比例系数，一般取值均为0.5，为第i次迭代的x方向搜索范围，为第i次迭代的y方向搜索范围，为第i次迭代的x方向搜索步长，为第i次迭代的y方向搜索步长，为第i+1次迭代的x方向搜索范围，为第i+1次迭代的y方向搜索范围，为第i+1次迭代的x方向搜索步长，为第i+1次迭代的y方向搜索步长。
[0139]
其它步骤及参数与具体实施方式一至七之一相同。
[0140]
具体实施方式九：本实施方式与具体实施方式一至八之一不同的是，所述步骤七的具体过程为：
[0141]
步骤七一、建立高程系统误差的拟合模型
[0142]
对于控制点的观测值集合对控制点进行重心化：
[0143][0144][0145][0146]
式中，为重心化控制点集合的观测值，为重心化控制点集合的已知值，为控制点集合的已知值，(xq,yq,zq)为摄影测量点云坐标，为摄影测量点云中第c个摄影测量点的坐标，c为摄影测量点云中包含的摄影测量点个数；
[0147]
以z方向为例，对于高程系统误差的二次多项式拟合模型：
[0148][0149]
其中，为重心化控制点集合中第d个控制点坐标的观测值，为重心化控制点集合中第d个控制点坐标的观测值，和为二次多项式的系数，为重心化控制点集合中第d个控制点坐标的已知值，vd为中间变量；
[0150]
将式(29)记为矩阵形式：
[0151][0152][0153][0154]
其中，d为控制点集合中包含的控制点个数；
[0155]
根据最小二乘原则，式(30)的平差解算结果为：
[0156][0157]
步骤七二、同理，采用步骤七一的方法分别获得x反向和y方向的二次多项式系数；
[0158]
步骤七三、将各个摄影测量点的坐标代入高程系统误差的x反向、y方向和z方向的二次多项式拟合模型以及平面的二次多项式拟合模型中，得到各个摄影测量点在x反向、y方向和z方向的系统误差估计值；
[0159]
再根据系统误差估计值对各个摄影测量点的坐标进行修正。
[0160]
零阶、一次和二次多项式是常用的拟合模型，其具体形式如下：
[0161][0162][0163][0164]
式中，ai为多项式的系数；
[0165]
对于控制点的观测值假设其平面坐标与系统误差(δx,δy,δz函数关系为一般多项式，两者的函数关系具体如下：
[0166][0167][0168][0169]
式中ε为随机误差；
[0170]
再根据估计值修正摄影测量点云坐标(xq+δx,yq+δy,zq+δz)。
[0171]
其它步骤及参数与具体实施方式一至八之一相同。
[0172]
实施例
[0173]
为了验证本发明方法的有效性，本发明采用资源三号影像和高分七号影像进行实验，如表1所示为卫星遥感影像实验数据参数：
[0174]
表1卫星遥感影像实验数据参数
[0175][0176][0177]
如表2所示，为资源三号影像的icesat-2点云滤波及高程提取精度数据：
[0178]
表2
[0179][0180]
如表3所示，为高分七号影像的icesat-2点云滤波及高程提取精度数据：
[0181]
表3
[0182][0183][0184]
由表2、表3可以看出，资源三号影像的icesat-2点云的滤波成功率最高为98.7％，最低为94.8％，高分七号影像的icesat-2点云的滤波成功率最高为99.11％，最低为95.73％，这表明点云滤波算法对不同地形变化、地物类型和激光强度的光子数据均可取得良好的滤波效果。资源三号影像的icesat-2点云提取的地物高程的rmse值最高为5.36m，最低为3.04m，mad值最高为4.20m，最低为2.23m，高分七号影像的icesat-2点云提取的地物高程的rmse值最高为4.06m，最低为2.58m，mad值最高为3.06m，最低为1.64m，这证明光子点云高程提取算法具备有效性，所提取的地物点高程可用作卫星摄影测量的高程控制点。
[0185]
如表4所示为资源三号摄影测量点云定位精度数据：
[0186]
表4
[0187][0188]
如表5所示为高分七号摄影测量点云定位精度数据：
[0189]
表5
[0190][0191]
[0192]
从表4、表5可知，资源三号摄影测量点云x方向的精度由31.49m提升至0.46m，y方向的精度由10.83m提升至2.43m，高程z方向的精度由7.28m提升至2.11m，高分七号摄影测量点云平面x方向的精度由2.95m提升至1.34m，y方向的精度无明显提升，高程z方向的精度由12.34m提升至0.59m。资源三号摄影测量点云与icesat-2点云的示意图如图9a所示，资源三号摄影测量点云系统误差修正前的示意图如图9b所示，资源三号摄影测量点云系统误差修正后的示意图如图9c所示。高分七号摄影测量点云与icesat-2点云的示意图如图10a所示，高分七号摄影测量点云系统误差修正前的示意图如图10b所示，高分七号摄影测量点云系统误差修正后的示意图如图10c所示。由此表明，本发明的一种光子点云辅助卫星摄影测量对地定位技术，对于提高山林地区摄影测量点云定位精度起到显著作用。
[0193]
本发明的上述算例仅为详细地说明本发明的计算模型和计算流程，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动，这里无法对所有的实施方式予以穷举，凡是属于本发明的技术方案所引伸出的显而易见的变化或变动仍处于本发明的保护范围之列。