2. 西安测绘研究所,陕西 西安 710054
2. Xi'an Research Institute of Surveying and Mapping, Xi'an 710054, China
空间后方交会是卫星摄影测量几何处理中的重要环节,是卫星遥感影像几何定位、区域网平差和自检校光束法平差等计算的关键步骤,其基本思想是利用合适的外方位元素模型描述成像传感器的位置与姿态,建立线阵影像的成像几何模型,对共线条件方程进行线性化后迭代求解外方位元素。当前,大多数卫星采用线阵CCD传感器成像,传感器在一定时间内的外方位元素通常被描述为某一成像行关于时间的多项式函数形式[1, 2, 3],或者若干离散成像时刻外方位元素的插值形式[4, 5, 6]。对于运行比较平稳的卫星,这两类模型均可较好地描述成像传感器的成像位置与姿态,参数之间的相关性成为影响后方交会精度的主要因素,其成因是线阵传感器的焦距较大、CCD线阵较短,使得成像时刻传感器的瞬时视场角较小,共线条件方程经线性化后产生了线元素和角元素的相关。
为了解决此类问题,学者们试图从方程解算方法和外方位元素建模这两个方面来提高外方位元素的求解精度与稳定性。在方程解算方面,文献[7]使用广义岭估计法求解外方位元素;文献[8] 采用岭-压缩组合估计的方法进行了参数解算;文献[9,10]对虚拟观测方程中权值的确定进行了改进。以上文献中均利用欧拉角描述角元素,因而这些算法仅能提高参数解算的精度,却无法消除或者降低参数之间的相关性。在外方位元素建模方面,文献[11]使用姿态四元数代替了欧拉角,建立了基于四元数球面线性(spherical linear interpolation,SLERP)的线阵影像外方位元素模型,实现了外方位元素的高精度解算;文献[12]利用四元数微分理论,提出了一种基于四元数微分方程(quaternion differential equation,QDE)的外方位元素模型,使用Tikhonov正则化法进一步提高了外方位元素求解的精度;文献[13]利用信噪比定量比较了各种四元数外方位元素模型在参数求解中的差异。虽然利用四元数可以有效降低参数之间的相关性,但仍是将平移与旋转分开考虑,位置与姿态之间没有严格的约束条件。
近些年,部分学者引入对偶四元数对外方位元素进行描述,利用其同时描述平移与旋转的特性,已在遥感影像几何处理中取得了初步的成果:文献[14]利用对偶四元数实现了面阵影像的空间后方交会;文献[15]验证了对偶四元数在区域网平差求解中的优势;文献[16]利用对偶四元数进行了面阵影像的直接定向;文献[17]利用线性蒙皮混合算法进行外方位元素建模,提高了GeoEye-1卫星遥感影像的立体定位结果,但是该卫星隐藏了定轨和定姿数据,因而该文献未研究定轨、定姿数据在后方交会中的应用,难以保证迭代在真值附近进行。
实际上,线阵卫星遥感影像的后方交会计算依赖于未知数的初值[18],应当在外方位元素准确建模的基础上,充分利用卫星定轨、定姿设备的观测结果,准确给定各类观测方程的权值[19],才能确保外方位元素的正确求解。本文试图利用对偶四元数同时描述成像传感器的位置和姿态,引入虚拟观测方程,进一步提高外方位元素的求解精度。
1 对偶四元数对偶四元数由单位四元数和对偶数理论发展起来,其定义为
式中,r=[r0 rx ry rz]T和s=[s0 sx sy sz]T分别称为对偶四元数的实部和对偶部;ε为对偶运算符,满足关系:ε≠0且ε2=0。r和s虽有8个参数,但是二者满足关系 因此对偶四元数中仅有6个自由度,与平移和旋转未知数的数量相同。利用对偶四元数的性质[20],旋转矩阵R与平移矢量t可利用对偶四元数的实部与对偶部分别表示为[16] 式中,rV=[rx ry rz]T;E为单位阵。将式(3)右侧进行矩阵乘法运算,可得到对偶四元数组成的旋转矩阵R的具体表达式为式(5)与单位四元数[11]所构成的旋转矩阵是一样的。同理可知,平移矢量t的具体表达式为
式中,(XS,YS,ZS)T为平移矢量的3个分量,在传统摄影测量中表示线元素。当线元素和姿态欧拉角已知时,对偶四元数实部r可按欧拉角计算单位四元数的方法确定[21],对偶部s可由利用矢量t和r计算为由式(5)和式(6)可知,两个坐标系之间的旋转矩阵R可以用由对偶四元数的实部r计算得到,坐标系原点之间的平移矢量t可由实部r与对偶s共同计算得到,因而利用对偶四元数建立线阵影像的外方位元素模型是可行的。
2 基于对偶四元数的外方位元素求解 2.1 基于对偶四元数的成像几何模型设q1=r1+εs1和q2=r2+εs2分别为影像首行和末行外方位元素所对应的对偶四元数,rk=[r0krxk ryk rzk]T与sk=[s0k sxk syk szk]T( =1,2)分别表示影像首行和末行外方位元素的实部和对偶部。在卫星成像过程中,卫星的位置和姿态变化比较平稳,可利用低阶插值计算第i个成像行的外方位元素为
式中,pt=i/n,n为首末两行之间成像行的数量;t1与t2分别表示首末两行的外方位线元素;C1和C2为四元数的插值系数。为了实现四元数在插值时角度匀速变化[6, 11],通常采用四元数球面线性插值,即 显然球面插值会导致C1和C2中包含三角函数与反三角函数的复合运算,使得线性化求解外方位元素的表达式过于复杂。由于卫星平台运行较为平稳,首末两行的姿态四元数差异很小,因而θ→0+,sinθ与θ是等价无穷小,此时插值系数可化简为式(9)表明,当卫星运行比较平稳,两成像行的成像间隔不太长时,可利用线性插值代替球面线性插值计算任意成像行的姿态。本文将利用式(9)所示的线性表达式构建的外方位元素模型称为线性插值的对偶四元数模型(linear interpolation based dual quaternion model,LDQ),将利用球面线性插值构建的外方位元素模型称为球面线性插值的对偶四元数模型(SLERP based dual quaternion model,SDQ)。
本文以LDQ模型为例推导求解外方位元素的误差方程式。不考虑主点偏移与焦距变化时,第个成像行的瞬时构像方程可以写成
式中,f为相机焦距;(X,Y)为像点坐标,可利用传感器沿轨道方向的指向角ψX和垂直于轨道方向的指向角ψY表示为(-ftan ψY,f tan ψX);系数aji、bij、cji(j=1,2,3)为该行对应的旋转矩阵中的元素,可内插计算ri后,由式(5)计算旋转矩阵得到;ti=[Xsi Ysi Zsi]T是成像行摄影中心在物方坐标系内的坐标,可由式(6)和式(8)联立计算得到。式(10)即为利用LDQ模型构建的线阵影像成像几何方程。为求解外方位元素,将式(10)在初值处展开至一次项,可得误差方程为
式中,lx=x-xc;ly=y-yc;xc=-fX/Z;yc=-fY/Z;cij(i=1,2;j=1,2,…,16)为线性化系数,可对m=[x y]T求ri和ti的偏导数,并依据复合函数微分法则计算s1、s2、r1和r2的偏导数为 则中不包含三角函数运算。式中,əm/əti(i=1,2)为线元素的偏导数,与航空摄影测量中的结果一致;əti/əsi可由对式(6)求偏导数计算得到;əm/əri的具体表达式在文献[14]和文献[15] 中有详细描述,在此不再赘述。误差方程(11)的矩阵形式为
式中,X=[ds1T ds2T drT1 dr2T]T为未知数改正数向量;V=[vx vy]T;l=[lx ly]T。此外,影像首行和末行对应的对偶四元数均满足式(2),对该式进行线性化并写成矩阵形式可得
式中误差方程(13)与条件方程(14)构成了外方位元素计算的基本平差模型。
2.2 对偶四元数虚拟观测方程的建立虚拟观测方程在摄影测量中有很高的实际应用价值,在线阵影像的后方交会和光束法平差中被广泛采用。由于该方法可以将外方位元素的改正数限制在一个合理的范围内,因而于有偏估计这种方法更为严密[18]。当前大多数卫星可以提供高精度的定轨数据,而定姿数据通常含有系统误差,且误差的构成较为复杂[22]。本文仅对线元素设置虚拟观测方程,适当放宽姿态改正数的限制。
对式(6)求微分可得
式中,Mr和Ms的表达式分别为对线元素增设的虚拟观测方程可表示为
式中,PM为虚拟观测方程的权矩阵。本文在迭代过程中动态调整权值[18],即先利用验前权值进行预平差,并进行方差分量估计值,然后计算第次迭代中参数di在虚拟观测方程中的权值为 式中,p2为像点观测值的方差估计值;k为正常数,通常取k=1;di2为未知数di第n-1次迭代时计算得到的方差估计值,可由X各分量的方差估计值按误差传播定律得到。将误差方程(13)、条件方程(14)以及虚拟观测方程(16)联立构成总误差方程式利用nG个地面控制点,按照式(18)总共可以列出2nG+6个误差方程和4个条件方程,应有nG≥3。当nG<6时,式(18)实质上是将定轨数据作为控制条件参与方程的求解,仅在精确给定初值时正确收敛,通常情况下应有nG≥6。利用约束条件的参数平差方法可求得改正数X,通过迭代即可求得利用对偶四元数描述的外方位元素值。
3 试验与分析为验证提出方法的正确性与有效性,本文将使用两组高分辨率卫星遥感影像进行对比试验验证。第1组试验数据(影像Ⅰ)为河南某地SPOT-5 HRS立体影像,共有13 376条扫描行,每个扫描行有12 000个像元,影像内包含25个平面高程控制点(图 1(a))。第2组试验数据(影像Ⅱ)为中国西部某地丘陵地区的SPOT-5 HRS立体影像,影像内共有12 000条扫描行,共包含10个平高控制点(图 1(b))。
3.1 试验1——影像I的外方位元素求解本试验用于比较对偶四元数外方位元素模型、欧拉角外方位模型和单位四元数外方位元素模型在外方位元素求解中的差异。为了验证式(9)的合理性,本试验还将比较SDQ模型和LDQ模型在后方交会计算中的区别。试验评价指标有两种:一是利用求解出的外方位元素反算地面点的像点坐标,并计算像点重投影误差的中误差;二是利用两张影像的外方位元素进行空间前方交会,计算三维坐标残差的中误差。试验中分别采用6、7、14和20个控制点求解外方位元素。表 1和表 2为不同控制点数量时检查点的像点坐标重投影的中误差及单位权方差统计结果;表 3和表 4为求解外方位元素值按空间前方交会计算地面点的三维坐标,计算检查点坐标残差的误差统计结果。
影像 | 控制 点数 | 检查 点数 | x/像素 | y /像素 | ||||||||
Euler | SLERP | QDE | SDQ | LDQ | Euler | SLERP | QDE | SDQ | LDQ | |||
前视影像 | 6 | 19 | 2.712 | 1.799 | 1.766 | 1.735 | 1.722 | 1.529 | 1.526 | 1.447 | 1.530 | 1.531 |
7 | 18 | 1.734 | 1.856 | 1.847 | 1.793 | 1.776 | 1.418 | 1.317 | 1.481 | 1.322 | 1.326 | |
14 | 11 | 2.202 | 1.667 | 1.598 | 1.518 | 1.469 | 1.256 | 1.257 | 1.452 | 1.261 | 1.262 | |
20 | 5 | 1.196 | 1.303 | 1.276 | 1.207 | 1.195 | 0.789 | 0.754 | 0.767 | 0.753 | 0.752 | |
后视影像 | 6 | 19 | 不收敛 | 1.582 | 1.414 | 1.554 | 1.538 | 不收敛 | 1.261 | 1.439 | 1.261 | 1.261 |
7 | 18 | 1.576 | 1.517 | 1.536 | 1.495 | 1.481 | 1.768 | 1.238 | 1.335 | 1.238 | 1.238 | |
14 | 11 | 1.676 | 1.656 | 1.485 | 1.614 | 1.585 | 1.554 | 1.444 | 1.487 | 1.444 | 1.444 | |
20 | 5 | 0.929 | 1.210 | 0.997 | 1.142 | 1.092 | 1.436 | 1.193 | 1.183 | 1.191 | 1.191 | |
注:Euler为利用欧拉角求解外方位元素的算法,SLERP为文献[11]中利用四元数球面线性插值求解外方位元素算法,QDE为文献[12]中的四元数微分方程算法,下同。 |
影像 | 控制 点数 | 检查 点数 | 平面/像素 | 0/像素 | ||||||||
Euler | SLERP | QDE | SDQ | LDQ | Euler | SLERP | QDE | SDQ | LDQ | |||
前视影像 | 6 | 19 | 3.114 | 2.359 | 2.283 | 2.313 | 2.304 | 0.391 | 2.573 | 2.007 | 0.994 | 0.965 |
7 | 18 | 2.240 | 2.276 | 2.370 | 2.227 | 2.216 | 0.476 | 1.889 | 1.566 | 1.155 | 1.131 | |
14 | 11 | 2.536 | 2.087 | 2.159 | 1.973 | 1.937 | 0.677 | 1.635 | 1.480 | 1.304 | 1.299 | |
20 | 5 | 1.433 | 1.505 | 1.489 | 1.423 | 1.412 | 0.877 | 1.762 | 1.565 | 1.374 | 1.358 | |
后视影像 | 6 | 19 | 不收敛 | 2.023 | 2.018 | 2.002 | 1.989 | 不收敛 | 0.916 | 0.839 | 0.517 | 0.509 |
7 | 18 | 2.368 | 1.958 | 2.035 | 1.941 | 1.931 | 0.224 | 0.998 | 0.668 | 0.739 | 0.736 | |
14 | 11 | 2.285 | 2.197 | 2.101 | 2.165 | 2.144 | 0.441 | 1.051 | 0.920 | 0.871 | 0.860 | |
20 | 5 | 1.710 | 1.699 | 1.548 | 1.651 | 1.616 | 0.436 | 1.672 | 1.279 | 1.182 | 1.163 |