基于卫星运动物理特性轨道外推的遥感影像几何校正方法与流程

本发明属于卫星遥感影像几何处理领域，具体涉及遥感影像几何校正方法。

背景技术：

基础地理信息是数据量最大、覆盖面最宽、应用面最广的信息资源之一。我国国民经济建设正处于高速发展期，地理环境要素，尤其是交通、城市建设等要素变化很快，对基础地理信息的快速获取和更新提出了新的要求。我国通过大力发展卫星遥感测绘技术，掌握实时获取各种空间信息的能力，为建立和维持国家高精度时空基准、及时更新各种比例尺的基础地理信息、制作现势性强的各类专题地图创造了必要条件，在提高我国管理决策水平、优化空间布局、制定发展战略与规划、加强自然资源管理等诸多方面产生了巨大的推动作用。

利用卫星遥感影像提取空间地理信息，关键技术环节是在实现影像定向即恢复影像获取时所处的空间位置和方位，重构像方坐标和物方坐标之间的几何转换关系基础上，对遥感影像进行几何校正，使校正影像与地理坐标系相符合。影像几何校正精度主要取决于遥感影像定向精度。如何利用适量控制点快速准确地实现高分辨率卫星遥感影像定向，一直是卫星摄影测量的难点。目前主要采用的基于外方位元素多项式的影像定向方法，其基本原理是以某种多项式分别拟合卫星传感器成像过程中线元素和角元素变化，通过一定数量的控制点解算多项式系数。对于最简单的一阶多项式，外方位元素拟合式有12个定向参数，至少要求6个控制点参与求解。而为了取得理想的定向精度，往往会选用二阶以上多项式以及更多控制点，譬如poli的扩展共线方程模型。而德国dlr采用的定向片模型中，使用3张定向片以二阶拉格郎日多项式对一景影像外方位元素拟合，包含18个定向参数，理论上需要9个控制点方可进行最小二乘无偏估计。尽管这种方法理论上可达到较高的定向精度，但实际由于参数多且相互间存在强相关，求解过程不仅要求提供良好的参数初值和分布均匀的控制点，还需另外采取能有效克服参数强相关的方法以保证解算的稳定，导致布设控制点的成本过高、求解过程复杂且最终的影像校正精度较差。

技术实现要素：

本发明是为解决目前卫星遥感影像几何处理中影像定向求解过程不稳定而导致的影像校正精度较差的问题。

基于卫星运动物理特性轨道外推的遥感影像几何校正方法，包括以下步骤：

步骤一、获取遥感卫星的轨道运动参数初值；

步骤二、构建像点坐标观测方程；

步骤三、对像点坐标观测方程进行线性化；

步骤四、利用最小二乘原理迭代计算卫星轨道运动参数估计值；

步骤五、利用轨道运动参数估计值外推生成待校正影像的rpc参数；

步骤六、利用rpc参数进行影像几何校正。

进一步地，步骤二所述构建像点坐标观测方程的过程包括以下步骤：

像点坐标观测方程表示为

其中，x,y为目标物体在卫星影像上的像方坐标，xp、yp、zp为目标物体在地球空间直角坐标系中的物方坐标；i,ω,tp,a,ω,e分别是轨道倾角、升交点赤经、近地点时刻、长半轴、近地点幅角、偏心率的值，α,k是卫星的三轴姿态角；f(i,ω,tp,a,ω,e,α,κ)为以卫星轨道运动参数表示的从物方坐标到像方坐标的转换函数；

在卫星运动轨道经过的控制点区域量测n个控制点目标，从而得到各个控制点的像方坐标观测值和物方坐标观测值，进而通过像方坐标观测值和物方坐标观测值确定像点坐标观测方程。

进一步地，步骤三所述对像点坐标观测方程进行线性化的过程包括以下步骤：

像点坐标观测方程线性化形式为：

v＝ax-l

其中，v＝[vxvy]^t为像方坐标残差向量；a为像方坐标对轨道倾角、升交点赤经、近地点时刻、长半轴、近地点幅角、偏心率和卫星三轴姿态角的偏导系数矩阵；为轨道倾角、升交点赤经、近地点时刻、长半轴、近地点幅角、偏心率和卫星三轴姿态角的改正数向量；l＝[lxly]^t为像方坐标观测值与像方坐标近似值的不符值向量。

进一步地，步骤四所述利用最小二乘原理迭代计算卫星轨道运动参数估计值的具体过程包括以下步骤：

根据最小二乘原理构造法方程组：

a^tax＝a^tl

解法方程，法方程的解即参数改正数为：

x＝(a^ta)^-1a^tl

当卫星三轴姿态角改正数δαδκ的绝对值小于规定的阈值时，即得到卫星轨道运动参数的估计值。

进一步地，步骤五所述利用轨道运动参数估计值外推生成待校正影像的rpc参数的过程包括以下步骤：

对于同一运行轨道上与控制点区域相邻的待校正影像覆盖区域，在该区域内生成m‘个虚拟控制点的物方坐标及其像方坐标；虚拟控制点物方坐标在待校正影像覆盖区域平面和高程方向上均匀分布，其像方坐标由像点坐标观测方程计算得到；

虚拟控制点的物方坐标和像方坐标的关联通过通用有理多项式函数表示为：

其中，(r,c)为正则化的影像坐标，(φ,λ,h)为正则化的物方坐标；影像坐标为(x,y)，物方坐标为(x,y,z)，其正则化公式如下：

其中，xoff、xscale、yoff、yscale为影像坐标的正则化参数；xoff、xscale、yoff、yscale、zoff、zscale为物方坐标的正则化参数；

p1、p2、p3、p4如下：

p1(φ,λ,h)＝a1+a2λ+a3φ+a4h+a5λφ+a6λh+a7φh+a8λ²+a9φ²+a10h²

+a11φλh+a12λ³+a13λφ²+a14λh²+a15λ²φ+a16φ³+a17φh²+a18λ²h

+a19φ²h+a20h³

p2(φ,λ,h)＝b1+b2λ+b3φ+b4h+b5λφ+b6λh+b7φh+b8λ²+b9φ²+b10h²+b11φλh

+b12λ³+b13λφ²+b14λh²+b15λ²φ+b16φ³+b17φh²+b18λ²h+b19φ²h

+b20h³

p3(φ,λ,h)＝c1+c2λ+c3φ+c4h+c5λφ+c6λh+c7φh+c8λ²+c9φ²+c10h²+c11φλh

+c12λ³+c13λφ²+c14λh²+c15λ²φ+c16φ³+c17φh²+c18λ²h+c19φ²h

+c20h³

p4(φ,λ,h)＝d1+d2λ+d3φ+d4h+d5λφ+d6λh+d7φh+d8λ²+d9φ²+d10h²

+d11φλh+d12λ³+d13λφ²+d14λh²+d15λ²φ+d16φ³+d17φh²+d18λ²h

+d19φ²h+d20h³

其中，ai’,bi‘,ci‘,di’即为rpc参数，i‘＝1,2,…,20；b1,d1取值为1。

进一步地，步骤六所述利用rpc参数进行影像几何校正的过程包括以下步骤：

对纠正影像上的每一点，利用rpc参数计算其在原影像上的采样点坐标，采用双线性内插法进行内插计算，然后计算该点的影像灰度值；对纠正影像上的全部点依次采样后，即完成一景影像的几何纠正。

进一步地，步骤六中实施双线性内插时，由采样点p周围4个像素的灰度值参加计算；内插点p的灰度值为

式中，i11，i12，i21，i22为周围4个像素的灰度值，wy1，wy2，wx1，wx2为相应的内插系数。

发明效果：

本发明方法利用遥感卫星运动物理规律，以卫星轨道运动参数代替摄影测量的外方位元素进行影像定向并外推得到待校正影像的rpc参数，再利用rpc参数进行影像几何校正。与传统方法相比，本发明能够极大的提高影像校正精度，经过本发明校正后，实施例中影像2的定位精度从535.95m提高到55.27m，精度提高约90％。效果显著。

如果在精度与传统方法精度相当的情况下，本发明无需校正影像区域的控制点，节省了卫星遥感影像几何处理的外业观测成本，具有良好的实用价值。

附图说明

图1为卫星椭圆轨道图；

图2为卫星轨道运动参数解算流程图；

图3为cbers待校正影像；

图4为cbers校正影像。

具体实施方式

具体实施方式一：结合图2对本实施方式进行说明，

基于卫星运动物理特性轨道外推的遥感影像几何校正方法，包括以下步骤：

步骤一、获取遥感卫星的轨道运动参数初值；

根据开普勒定律，在惯性直角坐标系下卫星运动轨道可通过6个轨道根数描述，分别为：轨道倾角i为卫星轨道平面与地球赤道面之间的夹角；升交点赤经ω为地球赤道平面上升交点与春分点之间的地心夹角；近地点幅角ω为轨道平面上近地点与升交点之间的地心角距；椭圆的长半轴a；椭圆的偏心率e；近地点时刻tp，如图1所示。在所述6个开普勒轨道根数中，轨道倾角i和升交点赤经ω决定轨道平面在空间的位置，近地点幅角ω、长半轴a和偏心率e决定椭圆在轨道平面上的方位、大小和形状。

遥感卫星的轨道运动参数初值通过卫星运营商公布的星历资料获得，但其存在一定系统误差，不能直接用于卫星摄影测量，需要通过后续步骤进行精化。

步骤二、构建像点坐标观测方程；

像点坐标观测方程表示为

其中，x,y为目标物体在卫星影像上的像方坐标，xp、yp、zp为目标物体在地球空间直角坐标系中的物方坐标；i,ω,tp,a,ω,e分别是轨道倾角、升交点赤经、近地点时刻、长半轴、近地点幅角、偏心率的值，α,κ是卫星的三轴姿态角；为以卫星轨道运动参数表示的从物方坐标到像方坐标的转换函数；

为了利用式(5)求解9个卫星轨道运动参数，要在卫星运动轨道经过的控制点区域量测n个控制点目标，n≥5，从而得到各个控制点的像方坐标观测值和物方坐标观测值，进而通过像方坐标观测值和物方坐标观测值确定像点坐标观测方程。

步骤三、对像点坐标观测方程进行线性化；

像点坐标观测方程是非线性函数。为了通过程序计算，将其进行线性化展开。其线性化形式为：

v＝ax-l(6)

其中，v＝[vxvy]^t为像方坐标残差向量；为像方坐标对轨道倾角、升交点赤经、近地点时刻、长半轴、近地点幅角、偏心率和卫星三轴姿态角的偏导系数矩阵；为轨道倾角、升交点赤经、近地点时刻、长半轴、近地点幅角、偏心率和卫星三轴姿态角的改正数向量；l＝[lxly]^t为像方坐标观测值与像方坐标近似值的不符值向量。

步骤四、利用最小二乘原理迭代计算卫星轨道运动参数估计值；

根据最小二乘原理构造法方程组：

a^tax＝a^tl(7)

解法方程，法方程的解即参数改正数为：

x＝(a^ta)^-1a^tl(8)

由于对像点坐标观测方程线性展开时舍弃了二次以上项，且参数初值是近似值，需要进行5-6次迭代计算，当卫星三轴姿态角改正数δαδκ的绝对值小于规定的阈值(取10^-5)时，即得到卫星轨道运动参数的估计值。

步骤五、利用轨道运动参数估计值外推生成待校正影像的rpc参数；

对于同一运行轨道上与控制点区域相邻的待校正影像覆盖区域，在该区域内生成10*10*5＝500个虚拟控制点的物方坐标及其像方坐标；虚拟控制点物方坐标在待校正影像覆盖区域平面和高程方向上均匀分布，其像方坐标由像点坐标观测方程式(5)计算得到；

虚拟控制点的物方坐标和像方坐标的关联通过通用有理多项式函数表示为：

其中，(r,c)为正则化的影像坐标，(φ,λ,h)为正则化的物方坐标；设影像坐标为(x,y)，物方坐标为(x,y,z)，其正则化公式如下：

其中，xoff、xscale、yoff、yscale为影像坐标的正则化参数；xoff、xscale、yoff、yscale、zoff、zscale为物方坐标的正则化参数；经正则化后，像方坐标和物方坐标取值位于(-1.0～+1.0)之间，这样可以减少计算过程中由于坐标数量级差别过大而引起的舍入误差。

p1、p2、p3、p4均为3次多项式，形式如下：

p1(φ,λ,h)＝a1+a2λ+a3φ+a4h+a5λφ+a6λh+a7φh+a8λ²+a9φ²+a10h²

+a11φλh+a12λ³+a13λφ²+a14λh²+a15λ²φ+a16φ³+a17φh²+a18λ²h

+a19φ²h+a20h³

p2(φ,λ,h)＝b1+b2λ+b3φ+b4h+b5λφ+b6λh+b7φh+b8λ²+b9φ²+b10h²+b11φλh

+b12λ³+b13λφ²+b14λh²+b15λ²φ+b16φ³+b17φh²+b18λ²h+b19φ²h

+b20h³

p3(φ,λ,h)＝c1+c2λ+c3φ+c4h+c5λφ+c6λh+c7φh+c8λ²+c9φ²+c10h²+c11φλh

+c12λ³+c13λφ²+c14λh²+c15λ²φ+c16φ³+c17φh²+c18λ²h+c19φ²h

+c20h³

p4(φ,λ,h)＝d1+d2λ+d3φ+d4h+d5λφ+d6λh+d7φh+d8λ²+d9φ²+d10h²

+d11φλh+d12λ³+d13λφ²+d14λh²+d15λ²φ+d16φ³+d17φh²+d18λ²h

+d19φ²h+d20h³

其中，ai’,bi‘,ci‘,di’即为rpc参数，i‘＝1,2,…,20；利用500个虚拟控制点的物方坐标及像方坐标列方程组求解即可，b1,d1取值为1。

步骤六、利用rpc参数进行影像几何校正。

采用间接法进行影像几何校正，对纠正影像上的每一点，利用rpc参数计算其在原影像上的采样点坐标，由于采样点坐标并非整数，采用双线性内插法进行内插计算，然后计算该点的影像灰度值。实施双线性内插时，由采样点p周围4个像素的灰度值参加计算；内插点p的灰度值为

式中，i11，i12，i21，i22为周围4个像素的灰度值，wy1，wy2，wx1，wx2为相应的内插系数。对纠正影像上的全部点依次采样后，即完成一景影像的几何纠正。

实施例：

实验数据使用来自同一运行轨道的两景cbers一级影像，像幅大小为10002pixels×10000pixels。影像1中心成像时间为2004年12月7日11h09min41.657275sec，影像2中心成像时间为2004年12月7日11h11min41.6572750sec。两景影像相差120sec，地面相距900km。在影像1覆盖地区用控制点进行轨道运动参数估计后，对影像2进行轨道外推几何校正，校正前后的定位精度如表1所示。

表1

经外推校正后，影像2的定位精度从535.95m提高到55.27m，精度提高约90％。效果显著。对中巴地球资源卫星cbers影像进行校正前后的结果分别如图3、图4所示。