一种利用高分七号激光测高数据辅助区域网平差的方法与流程

1.本发明涉及卫星数据处理技术领域，尤其涉及一种利用高分七号激光测高数据辅助区域网平差的方法。


背景技术：

2.提高稀少甚至无地面控制点的区域网平差精度，是实现境外和外业测控困难区域高精度测图的核心问题之一，也是主要的技术难点。高分七号激光测高数据具有极高的高程精度，并与立体影像同步获取，误差特性相近、相对精度较高，为了充分利用上述特点，发挥足印影像激光光斑地面信息的作用，亟待需要一种高分七号卫星激光测高数据辅助的立体影像联合区域网平差方法。


技术实现要素：

3.本发明的目的在于提供一种利用高分七号激光测高数据辅助区域网平差的方法，从而解决现有技术中存在的前述问题。
4.为了实现上述目的，本发明采用的技术方案如下：
5.一种利用高分七号激光测高数据辅助区域网平差的方法，包括以下步骤：
6.s1,采用高分七号激光测高系统获取激光测高原始数据，经过一系列精细改正处理过程后形成可用的激光测高标准产品，包括激光测高点的三维坐标、足印影像以及激光测高点足印影像像点坐标；
7.s2，将获取的立体影像构建区域网，并进行连接点匹配，获取数量足够、分布均匀的连接点，开展基于有理函数模型的自由网平差，实现区域网影像高精度相对定向，得到相对定向后的立体影像；
8.s3，根据步骤s1中获得的激光测高点的三维坐标、足印影像以及激光测高点足印影像像点坐标和步骤s2中获得的相对定向后的立体影像，开展基于激光测高点的高程控制点提取过程，获取符合提取要求的激光高程控制点；
9.s4,利用获取到的激光高程控制点作为高程控制，与相对定向后的立体影像开展联合区域网平差，实现立体影像高程精度提升；
10.s5，更新立体影像有理函数模型，获取精度提升后的立体影像。
11.优选的，步骤s1中所述一系列精细改正处理过程包括但不限于全波形处理、几何定位、大气以及潮汐误差改正处理过程。
12.优选的，步骤s2中所述的自由网平差包括以下步骤：
13.s21，利用卫星影像自带的rpc文件，构建有理函数模型，表达式如公式(1)所示：
14.15.其中，p
i
为一般多项式，表达式为：
16.p
i
＝a0 a1·
y a2·
x a3·
z a4·
y
·
x a5·
y
·
z a6·
x
·
z a7·
y2 a8·
x2 a9·
z2 a
10
·
x
·
y
·
z a
11
·
y3 a
12
·
y
·
x2 a
13
·
y
·
z2 a
14
·
y2·
x a
15
·
x3 a
16
·
x
·
z2 a
17
·
y2·
z a
18
·
x2·
z a
19
·
z317.其中，(x,y,z)为正则化后的地面点三维坐标，(x,y)为对应的正则化像点坐标，a
i
(i＝0,1,
…
,18,19)为有理函数系数,；
18.s22，采用仿射变换模型对公式(1)进行几何定位误差补偿，表达式如公式(2)所示：
[0019][0020]
其中，(δx,δy)为像方改正数，(x,y)为利用公式(1)求得的像点坐标，m
j
,n
j
(j＝0,1,2)为仿射变换模型参数；
[0021]
s23，影像仿射变换参数和连接点物方坐标是要求解的两类未知量，将公式(1)
‑
(2)线性化，可得误差方程，如公式(3)所示：
[0022]
v1＝a1x b1y1‑
l1,p1ꢀꢀ
(13)
[0023]
式中，x＝[δm
0 δm
1 δm
2 δn
0 δn
1 δn2]
t
，是影像的仿射变换参数改正数，y1＝[δlon δlat δh]
t
是连接点物方坐标改正值，a1,b1为对应的系数矩阵，l1为初值计算的常量，p1为权值。
[0024]
优选的，步骤s4中所述的联合区域网平差具体包括：
[0025]
s41，利用卫星影像自带的rpc文件，构建有理函数模型，表达式如公式(1)所示：
[0026][0027]
其中，p
i
为一般多项式，表达式为：
[0028]
p
i
＝a0 a1·
y a2·
x a3·
z a4·
y
·
x a5·
y
·
z a6·
x
·
z a7·
y2 a8·
x2 a9·
z2 a
10
·
x
·
y
·
z a
11
·
y3 a
12
·
y
·
x2 a
13
·
y
·
z2 a
14
·
y2·
x a
15
·
x3 a
16
·
x
·
z2 a
17
·
y2·
z a
18
·
x2·
z a
19
·
z3[0029]
其中，(x,y,z)为正则化后的地面点三维坐标，(x,y)为对应的正则化像点坐标，a
i
(i＝0,1,
…
,18,19)为有理函数系数；
[0030]
s42，采用仿射变换模型对公式(1)进行几何定位误差补偿，表达式如公式(2)所示：
[0031][0032]
其中，(δx,δy)为像方改正数，(x,y)为利用公式(1)求得的像点坐标，m
j
,n
j
(j＝0,1,2)为仿射变换模型参数；
[0033]
s43，影像仿射变换参数和连接点物方坐标是要求解的两类未知量，将公式(1)
‑
(2)线性化，可得误差方程，如公式(3)所示：
[0034]
v1＝a1x b1y1‑
l1,p1ꢀꢀ
(15)
[0035]
式中，x＝[δm
0 δm
1 δm
2 δn
0 δn
1 δn2]
t
，是影像的仿射变换参数改正数，y1＝[δlon δlat δh]
t
是连接点物方坐标改正值，a1,b1为对应的系数矩阵，l1为初值计算的常量，p1为权值；
[0036]
s44，对于激光测高点，其高程作为控制，在平差中仅改正其平面坐标，此时，误差方程变为：
[0037]
v2＝b2y2‑
l2,p2ꢀꢀ
(16)
[0038]
式中，y2＝[δlon δlat]
t
是激光测高点物方坐标改正值，b2是相应的系数矩阵，l2是初值计算常量，p2是相应权值。
[0039]
优选的，步骤s4还包括：当存在其他外业控制点时，开展高程控制点、外业控制点与立体影像的联合区域网平差，实现立体影像高程和平面精度同步提升，此时还需要将外业控制点参与联合区域网平差，所涉及到的误差方程为：
[0040]
v3＝a2x
‑
l3,p3ꢀꢀ
(17)
[0041]
式中，x＝[δm
0 δm
1 δm
2 δn
0 δn
1 δn2]
t
是影像的仿射变换参数改正数，a2是相应的系数矩阵，l3是初值计算常量，p3是相应权值。
[0042]
优选的，步骤s3中在进行高程控制点提取过程包括a方案和b方案，对于每个激光测高点，首先采用a方案进行高程控制点的提取，当a方案高程控制点提取过程失败之后则采用b方案进行高程控制点的提取过程。
[0043]
优选的，所述a方案包括：
[0044]
a31，将激光测高点物方三维坐标代入立体影像有理函数模型，计算得到当前激光测高点在立体影像上的初略位置，设其像点坐标为(p,q)，则该点与步骤s1中的激光测高点足印影像像点坐标(x0,y0)构成共轭点对；
[0045]
a32，以像点坐标(p,q)为中心，以激光光斑半径为最大搜索窗口，按公式(6)逐点计算与以激光测高点足印影像像点坐标(x0,y0)为中心的区域影像相关系数，取最大相关系数点位(x
1'
,y
1'
)作为像素级配准点位：
[0046][0047]
式中，g,g'分别为足印影像和立体影像灰度值，为匹配窗口的平均灰度值，w,h分表表示相关系数搜索窗口的宽和高；
[0048]
a33，以坐标(x0,y0)，(x
1'
,y
1'
)为初值，按公式(7)进行最小二乘匹配，获取子像素级配准点位(x1,y1)，即激光测高点在落点影像上的像点坐标:
[0049][0050]
式中，h0,h1为影像辐射畸变参数，m'
j
,n'
j
为影像几何变换参数,j＝0,1,2；
[0051]
a34，重复步骤a31
‑
s33，获取当前激光测高点在所有落点影像上的像点坐标；
[0052]
a35，重复步骤a31
‑
a34，获取所有激光测高点在立体影像中的像点坐标。
[0053]
优选的，所述b方案包括：
[0054]
b31,根据激光测高点三维坐标，通过影像几何模型正变换，获取激光测高点在立体影像上的像点坐标概略位置；
[0055]
b32,以步骤b31中求得的像点坐标为原点，在激光光斑半径范围内对立体影像与足印影像进行同名点匹配，经过粗差剔除后，得到一定数量的候选点；
[0056]
b33，利用公式(6)计算所有候选点与激光测高点足印影像像点的相关系数，从中选择离激光光斑中心最近且相关系数最大的点对应的像点坐标作为当前激光测高点在落点影像上的像点坐标，；
[0057]
b34，重复步骤b31
‑
b33，获取当前激光测高点在所有落点影像上的像点坐标；
[0058]
b35,重复步骤b31
‑
b34，获取所有激光测高点在立体影像中的像点坐标
[0059]
优选的，步骤s2中的立体影像包括但不限于采用高分七号卫星获取的立体影像。
[0060]
本发明的有益效果是：
[0061]
本发明公开了一种利用高分七号激光测高数据辅助区域网平差的方法，该方法结合高分七号激光测高数据中的足印影像，实现了基于激光测高数据的高程控制点自动提取，建立了激光测高数据高程控制点辅助立体影像的联合区域网平差模型，实现了无控、外业控制以及高程控制条件下的模型统一；将激光测高数据高程控制点作为高程控制应用到平差计算中，实现区域影像几何精度提升。通过本发明中记载的方法能够在覆盖不同地形类型的区域，仅使用激光测高数据作为高程控制，便可实现区域网高程精度的极大提升。
附图说明
[0062]
图1是实施例1中提供的利用高分七号激光测高数据辅助区域网平差的方法流程示意图；
[0063]
图2是实施例1中测区地形与影像、外业控制点和激光测高点分布示意图；
[0064]
图3是实施例1中自由网平差平面与高程残差图；
[0065]
图4是实施例1中激光测高点做高程控制的平面和高程残差图。
具体实施方式
[0066]
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施方式仅仅用以解释本发明，并不用于限定本发明。
[0067]
实施例1
[0068]
本实施例提供了一种利用高分七号激光测高数据辅助区域网平差的方法，如图1所示，包括以下步骤：
[0069]
s1,采用高分七号激光测高系统获取激光测高原始数据，经过一系列精细改正处理过程后形成可用的激光测高标准产品，包括激光测高点的三维坐标、足印影像以及激光测高点足印影像像点坐标；
[0070]
s2，将获取的立体影像构建区域网，并进行连接点匹配，获取数量足够、分布均匀的连接点，开展基于有理函数模型的自由网平差，实现区域网影像高精度相对定向，得到相
对定向后的立体影像；
[0071]
s3，根据步骤s1中获得的激光测高点的三维坐标、足印影像以及激光测高点足印影像像点坐标和步骤s2中获得的相对定向后的立体影像，开展基于激光测高点的高程控制点提取过程，获取符合提取要求的激光高程控制点；
[0072]
s4,利用获取到的激光高程控制点作为高程控制，与相对定向后的立体影像开展联合区域网平差，实现立体影像高程精度提升；
[0073]
s5，更新立体影像有理函数模型，获取精度提升后的立体影像。
[0074]
本实施例中的步骤s1中，获取激光测高数据具体为：在工作状态下，高分七号卫星激光测高仪2个激光束以3hz频率向地面发射激光脉冲，在地面形成沿轨向间隔约2.4km、垂轨向约12.5km的离散激光足印光斑，在幅宽20km的光学影像范围内约有2列共16个激光点。
[0075]
高分七号卫星激光测高仪同时配置了激光足印相机，用于同步获取激光光斑地面影像，经过后期处理，可得出激光光斑中心在足印影像上的位置。
[0076]
本实施例中通过高分七号卫星激光测高系统获得的原始数据在经过全波形处理、几何定位以及大气、潮汐等各项改正处理后形成激光测高产品，并提供各测高点的三维坐标及其足印影像坐标；高分七号卫星激光测高数据平面采用wgs84椭球基准，高程为wgs84椭球高。
[0077]
本实施例中，步骤s2中所述的自由网平差包括以下步骤：
[0078]
s21，利用卫星影像自带的rpc文件，构建有理函数模型，表达式如公式(1)所示：
[0079][0080]
其中，p
i
为一般多项式，表达式为：
[0081]
p
i
＝a0 a1·
y a2·
x a3·
z a4·
y
·
x a5·
y
·
z a6·
x
·
z a7·
y2 a8·
x2 a9·
z2 a
10
·
x
·
y
·
z a
11
·
y3 a
12
·
y
·
x2 a
13
·
y
·
z2 a
14
·
y2·
x a
15
·
x3 a
16
·
x
·
z2 a
17
·
y2·
z a
18
·
x2·
z a
19
·
z3[0082]
其中，(x,y,z)为正则化后的地面点三维坐标，(x,y)为对应的正则化像点坐标，a
i
(i＝0,1,
…
,18,19)为有理函数系数；
[0083]
s22，采用仿射变换模型对公式(1)进行几何定位误差补偿，表达式如公式(2)所示：
[0084][0085]
其中，(δx,δy)为像方改正数，(x,y)为利用公式(1)求得的像点坐标，m
j
,n
j
(j＝0,1,2)为仿射变换模型参数；
[0086]
s23，影像仿射变换参数和连接点物方坐标是要求解的两类未知量，将公式(1)
‑
(2)线性化，可得误差方程，如公式(3)所示：
[0087]
v1＝a1x b1y1‑
l1,p1ꢀꢀ
(22)
[0088]
式中，x＝[δm
0 δm
1 δm
2 δn
0 δn
1 δn2]
t
，是影像的仿射变换参数改正数，y1＝
[δlon δlat δh]
t
是连接点物方坐标改正值，a1,b1为对应的系数矩阵，l1为初值计算的常量，p1为权值。
[0089]
本实施例中的步骤s3中在进行高程控制点提取过程包括a方案和b方案，对于每个激光测高点，首先采用a方案进行高程控制点的提取，当a方案高程控制点提取过程失败之后则采用b方案进行高程控制点的提取过程。
[0090]
所述a方案包括：
[0091]
a31，将激光测高点物方三维坐标代入立体影像有理函数模型，计算得到当前激光测高点在立体影像上的初略位置，设其像点坐标为(p,q),则该点与步骤s1中的激光测高点像点坐标(x0,y0)构成共轭点对；
[0092]
a32，以像点坐标(p,q)为中心，以激光光斑半径为最大搜索窗口，按公式(6)逐点计算与以激光测高点像点坐标(x0,y0)为中心的区域影像相关系数，取最大相关系数点位(x
1'
,y
1'
)作为像素级配准点位：
[0093][0094]
式中，g,g'分别为足印影像和立体影像灰度值，为匹配窗口的平均灰度值，w,h分表表示相关系数搜索窗口的宽和高；
[0095]
a33，以坐标(x0,y0),(x
1'
,y
1'
)为初值，按公式(7)进行最小二乘匹配，获取子像素级配准点位(x1,y1):
[0096][0097]
式中，h0,h1为影像辐射畸变参数，m'
j
,n'
j
为影像几何变换参数,j＝0,1,2；
[0098]
a34，重复步骤a31
‑
s33，获取当前激光测高点在所有落点影像上的像点坐标；
[0099]
a35，重复步骤a31
‑
a34，获取所有激光测高点在立体影像中的像点坐标。
[0100]
所述b方案包括：
[0101]
b31,根据激光测高点三维坐标，通过影像几何模型正变换，获取激光测高点在立体影像上的像点坐标位置；
[0102]
b32,以步骤b31中求得的像点坐标为原点，在激光光斑半径范围内对立体影像与足印影像进行同名点匹配，经过粗差剔除后，得到一定数量的候选点；
[0103]
b33，从候选点中选择离激光光斑中心最近且相关系数最大的点作为最终的激光测高点像点坐标，相关系数由公式(6)计算确定；
[0104]
b34，重复步骤b31
‑
b33，获取当前激光测高点在所有落点影像上的像点坐标；
[0105]
b35,重复步骤b31
‑
b34，获取所有激光测高点在立体影像中的像点坐标。
[0106]
在进行高程控制点的提取过程中，由于a、b方案的混合使用，可能使得一个激光测高点在不同的影像上对应不同的物方地物，也就是一个激光测高点对应多个高程控制点，这就需要对每一个高程控制点进行独立平差操作，如果此时仍作为单个高程控制点进行平差计算，则会产生平差误差，这个误差的计算方法可以采用以下公式进行计算：
[0107][0108]
式中，δh为高程误差，y
e
为影像地面分辨率，δr为像点量测误差，b/h为基高比。
[0109]
对于高分七号卫星来说，前后视沿轨向夹角分别为26度和5度，基高比为0.575，则影像上高程控制点0.5像素的坐标差异便会引起约0.8米的高程误差，这远远大于激光测高点的标称高程精度，因此要特别重视同一激光测高点衍生多个高程控制点的情况。
[0110]
步骤s4中所述的联合区域网平差具体包括：
[0111]
s41，利用卫星影像自带的rpc文件，构建有理函数模型，表达式如公式(1)所示：
[0112][0113]
其中，p
i
为一般多项式，表达式为：
[0114]
p
i
＝a0 a1·
y a2·
x a3·
z a4·
y
·
x a5·
y
·
z a6·
x
·
z a7·
y2 a8·
x2 a9·
z2 a
10
·
x
·
y
·
z a
11
·
y3 a
12
·
y
·
x2 a
13
·
y
·
z2 a
14
·
y2·
x a
15
·
x3 a
16
·
x
·
z2 a
17
·
y2·
z a
18
·
x2·
z a
19
·
z3[0115]
其中，(x,y,z)为正则化后的地面点三维坐标，(x,y)为对应的正则化像点坐标，a
i
(i＝0,1,
…
,18,19)为有理函数系数；
[0116]
s42，采用仿射变换模型对公式(1)进行几何定位误差补偿，表达式如公式(2)所示：
[0117][0118]
其中，(δx,δy)为像方改正数，(x,y)为利用公式(1)求得的像点坐标，m
j
,n
j
为仿射变换模型参数，j＝0,1,2；
[0119]
s43，影像仿射变换参数和连接点物方坐标是要求解的两类未知量，将公式(1)
‑
(2)线性化，可得误差方程，如公式(3)所示：
[0120]
v1＝a1x b1y1‑
l1,p1ꢀꢀ
(27)
[0121]
式中，x＝[δm
0 δm
1 δm
2 δn
0 δn
1 δn2]
t
，是影像的仿射变换参数改正数，y1＝[δlon δlat δh]
t
是连接点物方坐标改正值，a1,b1为对应的系数矩阵，l1为初值计算的常量，p1为权值；
[0122]
s44，对于激光测高点，其高程作为控制，在平差中仅改正其平面坐标，此时，误差方程变为：
[0123]
v2＝b2y2‑
l2,p2ꢀꢀ
(28)
[0124]
式中，y2＝[δlon δlat]
t
是激光测高点物方坐标改正值，b2是相应的系数矩阵，l2是初值计算常量，p2是相应权值。
[0125]
当存在其他外业控制点时，开展高程控制点、外业控制点与立体影像的联合区域网平差，实现立体影像高程和平面精度同步提升,此时还需要将外业控制点参与联合区域
网平差，所涉及到的误差方程为：
[0126]
v3＝a2x
‑
l3,p3ꢀꢀ
(29)
[0127]
式中，x＝[δm
0 δm
1 δm
2 δn
0 δn
1 δn2]
t
是影像的仿射变换参数改正数，a2是相应的系数矩阵，l3是初值计算常量，p3是相应权值。
[0128]
采用上述的方案，选择山东省中部某区域作为实验区域，该区域面积约2.75万平方公里，整体南北长东西窄。区域整体地形起伏从海拔
‑
5米到1000米，其中中部是山地和丘陵，周边以平地为主。共收集8轨共80对高分七号立体影像，拍摄时间集中在2020年4月至2021年1月。共收集同轨获取激光测高点5轨602个，以及异轨获取激光测高点1轨67个。经过高程控制点提取，共获得328个有效高程控制点。为对平差结果进行验证分析，在试验区测量了109个高精度外业控制点(平面和高程精度均优于0.15米)，所有控制点像点坐标均在立体环境下量测，量测精度优于0.3像素。整个测区地形以及控制点、激光测高点分布如图2所示。
[0129]
根据地形坡度以及高差信息，将山东测区划分为平地、丘陵和山地等3个不同地形。将激光测高点作为高程控制点与立体影像进行联合区域网平差，验证激光测高点在不同地形条件下对立体影像几何定位精度的影响效果，其结果如下表所示：
[0130]
山东测区激光测高点对不同地形的平差结果统计表
[0131][0132]
从上表可以看出，在没有任何外部控制条件下，高分七号立体影像整体平面定位精度在9.5米左右，其中山地区域的平面精度好于平地和丘陵，并且最大误差也较小，而平地和丘陵区域精度相当。高程方面，自由网平差后的高程精度在7.9米左右，虽然山地区域高程精度好一些，但是系统误差比较明显，全区域高程系统误差达到了5米。从图3也可以看出，测区平面和高程误差呈现明显的系统误差。经过激光测高点辅助区域网平差后，测区整体平面精度没有较大变化，平面精度依然在8.7米左右，平面最大误差也与自由网平差时相近，同时也没有明显的地形区别，但是高程精度迅速提升，整体达到0.76米，特别是山地区域，高程精度达到0.46米，而高程最大误差也降到了1.72米，相比于自由网平差，高程精度提升非常明显。
[0133]
从图4可以看出，因为没有平面约束，联合平差后平面误差依然显示为系统误差，
整体误差趋势与图3相似，但是高程误差规律不明显，数值降低很多。以上实验结果可以看出本发明提供的方法使用激光测高点对高程精度具有较好的提升作用，特别是山地区域，效果更加明显。
[0134]
通过采用本发明公开的上述技术方案，得到了如下有益的效果：
[0135]
本发明公开了一种利用高分七号激光测高数据辅助区域网平差的方法，该方法结合高分七号激光测高数据中的足印影像，实现了基于激光测高数据的高程控制点自动提取，建立了激光测高数据高程控制点辅助立体影像的联合区域网平差模型，实现了无控、外业控制以及高程控制条件下的模型统一；将激光测高数据高程控制点作为高程控制应用到平差计算中，实现区域影像几何精度提升。通过本发明中记载的方法能够在覆盖不同地形类型的区域，仅使用激光测高数据作为高程控制，便可实现区域网高程精度的极大提升。
[0136]
以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视本发明的保护范围。