基于平差参数快速解算的超大规模区域网平差方法及系统

摘要

本发明提供一种基于平差参数快速解算的超大规模区域网平差方法及系统，包括建立改化法方程，基于共轭梯度快速解算改化法方程；根据改化法方程解算的参数更新附加像方参数，根据像方附加参数和影像RPC参数交会解算连接点物方坐标，根据当前像方附加参数和连接点物方坐标更新改化法方程，直到收敛后根据解算的像方附加参数进行RPC参数的精化，得到平差结果。该技术方案解决了大规模方程组不易解算的难题，同时兼顾计算资源开销，可以满足超大规模区域网平差解算的效率需求。

说明书

技术领域

本发明属于遥感影像处理领域，涉及一种基于平差参数快速解算的超大规模区域网平差 方法及系统。

背景技术

对于包含万景影像以上的超大规模区域网而言，即便先消掉连接点物方坐标这一类未知 数后，各景影像的待平差参数也高达十万个，改化等效法方程系数矩阵的阶数也达到十万阶， 区域网内部影像之间直接或间接的拓扑连接关系错综复杂，使得传统消元法消元过程中，系 数矩阵的稀疏性逐渐降低，稠密性逐渐提高，反而大大增加了解算的复杂度。

在超大规模区域网平差中改化法方程为一个高阶的线性方程组，对其进行求解的效率很 大程度上取决于线性代数方程组的解法。在摄影测量与遥感领域对于误差方程组的解算一般 采用逐点法化，直接解算的方式，直接解算的特点是，对于一个给定的误差方程组，事先可 以按规定的算法步骤计算出它所需的算术运算操作数，直接给出最后的结果，然而直接解算 的方法受到问题规模的限制，所以这种方法对于解算超大规模光学卫星影像区域网平差参数 并不适用。对于大规模的问题多采用迭代解算的方法，迭代解法的特点是，对于一个给定的 线性代数方程组，首先假设一个初始解，然后按一定的算法公式进行迭代。在每次迭代中对 解的误差进行检查，并通过增加迭代次数不断降低解的误差，直至满足解的精度要求，输出 最后的结果。如果能实现平差参数的快速解算，可以提高平差效率。但目前尚未有相关技术 方案出现。

发明内容

本发明所要解决的问题是，针对超大规模光学卫星影像区域网平差问题，提出一种基于 平差参数快速解算的超大规模区域网平差方法及系统。

本发明的技术方案为一种基于平差参数快速解算的超大规模区域网平差方法，包括以下 步骤：

步骤1，建立改化法方程如下，

Ax＝b

其中，A为改化法方程系数矩阵，若区域网中共有m张像片，则A的阶数为6m×6m，x表示 像方附加参数的改正数，为待解求参数，b为一个大小为6m的列向量；

步骤2，基于共轭梯度快速解算步骤1所得改化法方程；

步骤3，判断迭代过程是否收敛，是则进入步骤7，否则进入步骤4；

步骤4，根据改化法方程解算的参数更新附加像方参数如下，

X＝X⁰+dX

其中，X为平差待解求参数，X⁰为平差待解算参数的初值，dX为平差解算的改正数；

步骤5，根据像方附加参数和影像RPC参数交会解算连接点物方坐标；

步骤6，根据当前像方附加参数和连接点物方坐标更新改化法方程，返回步骤2；

步骤7，根据解算的像方附加参数进行RPC参数的精化，得到平差结果。

而且，步骤7中，根据解算的像方附加参数进行RPC参数的精化，实现方式如下，

首先生成地面虚拟控制点，包括在影像像平面上均匀划分规则格网，对每个格网中心像 点p(samp,line)，利用该影像内定向参数，通过前方交会在物方局部一系列高程基准面上交 会得到一系列物方点P_i，此时，像点p与各物方点P_i构成虚拟控制点；

然后以RPC模型为误差模型，RPC参数为待解算参数，在像方加入附加参数，根据原始 RPC列误差方程式，采用最小二乘的方法平差解算未知参数，得到精化后的RPC参数。

本发明相应提供一种基于平差参数快速解算的超大规模区域网平差系统，包括以下模块：

第一模块，用于建立改化法方程如下，

Ax＝b

其中，A为改化法方程系数矩阵，若区域网中共有m张像片，则A的阶数为6m×6m，x表示 像方附加参数的改正数，为待解求参数，b为一个大小为6m的列向量；

第二模块，用于基于共轭梯度快速解算第一模块所得改化法方程；

第三模块，用于判断迭代过程是否收敛，是则命令第七模块工作，否则命令第四模块工 作；

第四模块，用于根据改化法方程解算的参数更新附加像方参数如下，

X＝X⁰+dX

其中，X为平差待解求参数，X⁰为平差待解算参数的初值，dX为平差解算的改正数；

第五模块，用于根据像方附加参数和影像RPC参数交会解算连接点物方坐标；

第六模块，用于根据当前像方附加参数和连接点物方坐标更新改化法方程，命令第二模 块工作；

第七模块，用于根据解算的像方附加参数进行RPC参数的精化，得到平差结果。

而且，第七模块中，根据解算的像方附加参数进行RPC参数的精化，实现方式如下，

首先生成地面虚拟控制点，包括在影像像平面上均匀划分规则格网，对每个格网中心像 点p(samp,line)，利用该影像内定向参数，通过前方交会在物方局部一系列高程基准面上交 会得到一系列物方点P_i，此时，像点p与各物方点P_i构成虚拟控制点；

然后以RPC模型为误差模型，RPC参数为待解算参数，在像方加入附加参数，根据原始 RPC列误差方程式，采用最小二乘的方法平差解算未知参数，得到精化后的RPC参数。

本发明的优点在于：在平差解算时充分利用法方程系数矩阵对称正定的特点，采用基于 共轭梯度的最优解搜索算法迭代解算改化法方程组，并根据解算的参数更新像方附加参数， 再前方交会出连接点的物方坐标，当像方附加参数趋于稳定时，停止计算，该技术方案解决 了大规模方程组不易解算的难题，同时兼顾计算资源开销，可以满足超大规模区域网平差解 算的效率需求。实验表明该技术方案可行、有效，平差解算结果稳定、可靠。

附图说明

图1为本发明实施例的流程示意图。

具体实施方式

以下结合附图和实施例详细说明本发明具体实施方式。

鉴于超大规模光学卫星影像区域网平差中改化法方程系数矩阵的阶数较大且该矩阵为一 正定对称的方阵，本发明结合共轭梯度的迭代解算方法提出一种超大规模光学卫星影像区域 网平差参数的快速解算方法，利用共轭梯度的方法迭代解算改化法方程的参数，再根据解算 的参数根据影像RPC参数前方交会便可得到连接点物方坐标。

参见图1，实施例的流程可以分为七个步骤，每个步骤实施的具体方法、公式以及流程 如下：

1.改化法方程式的建立

对所有连接点像点和控制点像点分别构建观测误差方程，根据最小二乘平差原理进行法 化得到法方程；采用消元改化法方程的策略来进行平差解算，包括先消去连接点坐标，构建 仅包含附加模型参数的改化法方程并简化。

具体实施时，可以预先进行改化法方程式的建立，或者利用已有的改化法方程式。

设非正则化的地面点大地坐标为(Lat,Lon,Height)，像点坐标为(s,l)，根据影像的RPC模 型建立误差方程式，将各景待平差影像的RPC模型像方附加一仿射变换模型，则平差模型公 式为：

$\left(\begin{array}{c}l+\Delta l=F_x(Lat,Lon,Height)\\ s+\Delta s=F_y(Lat,Lon,Height)\end{array}\right)---\left(1\right)$

$\left(\begin{array}{c}\Delta l=a_0+a_{1}l+a_{2}s\\ \Delta s=b_0+b_{1}s+b_{2}s\end{array}\right)---\left(2\right)$

式中，F_x(Lat,Lon,Height)、F_y(Lat,Lon,Height)分别为沿轨方向和垂轨方向的像点坐标函数模型， Δl、Δs沿轨方向和垂轨方向的像方附加参数模型，a_i,b_i(i＝1,2,3)为仿射变换参数。

对于控制点像点而言，由于其对应的物方点坐标精确已知，因此，所构建的误差方程式 中未知参数仅包括该像点所在影像的RPC模型像方附加参数，显然，对于RPC模型像方附 加参数而言，此时式为线性方程而无需进行线性化处理，如式(5)：

$\left(\begin{array}{c}{v}_l=F_x(Lat,Lon,Height)-l-\Delta l\\ v_s=F_y(Lat,Lon,Height)-s-\Delta s\end{array}\right)---\left(3\right)$

式中，v_l、v_s分别为影像上沿轨和垂轨方向的改正数。

对于连接点像点而言，由于其对应的物方点坐标未知，因此，所构建的误差方程式中未 知参数除了包括该像点所在影像的RPC模型像方附加参数外，还包括其对应的连接点物方坐 标(Lat,Lon,Height)。对于连接点物方坐标(Lat,Lon,Height)而言，式(3)为一非线性方程，需要对 其赋予初值(Lat,Lon,Height)⁰并进行线性化处理，如式(4)：

$\left(\begin{array}{c}{v}_l=F_x{\left(Lat,Lon,Height\right)}^0+\frac{\partial F_x}{\partial \left(Lat,Lon,Height\right)}{|}_{{\left(Lat,Lon,Height\right)}^0}d\left(Lat,Lon,Height\right)-l-\Delta l\\ v_s=F_y{\left(Lat,Lon,Height\right)}^0+\frac{\partial F_y}{\partial \left(Lat,Lon,Height\right)}{|}_{{\left(Lat,Lon,Height\right)}^0}d\left(Lat,Lon,Height\right)-s-\Delta s\end{array}\right)---\left(4\right)$

具体实施时，本领域技术人员可自行设定初值(Lat,Lon,Height)⁰，或者设定初值求取规则， 可采用现有技术，例如使用初始定向参数，采用前方交会的方式确定初值。

对所有连接点像点及控制点像点分别构建观测误差方程，并写成矩阵形式：

V＝Mx+Wt-LP(5)

其中，V代表像点坐标观测值残差向量；x＝[X₁…X_i…X_m]^T(i＝1,2…m)代表各 景影像RPC像方附加模型参数向量，X_i＝(a₀,a₁,a₂,b₀,b₁,b₂)_i代表影像Img_i的RPC模型像方 附加参数向量，m代表待平差影像数；t＝[T₁…T_j…T_n]^T(j＝1,2…n)代表各连接点物 方坐标改正值向量，T_j＝d(Lat,Lon,Height)_j代表连接点TP_j的物方坐标改正数，n代表连接 点个数；M、W则分别为对应未知数的偏导数系数矩阵，L和P分别为相应的常向量和权矩 阵。

根据最小二乘平差原理，对观测误差方程进行法化，可得到法方程如式所示：

$\left(\begin{array}{cc}{M}^{T}PM& M^{T}PW\\ W^{T}PM& W^{T}PW\end{array}\right)\left(\begin{array}{c}x\\ t\end{array}\right)=\left(\begin{array}{c}{M}^{T}PL\\ W^{T}PL\end{array}\right)---\left(6\right)$

当进行大规模区域网平差时，由于参与平差的影像以及连接点的数量较大，上式左边的 法方程系数矩阵阶数较高，直接通过对其进行求逆来解算各项未知参数，不论是内存开销还 是解算效率上都无法满足要求。本发明中采用消元改化法方程的策略来进行平差解算，考虑 连接点物方坐标t的维数通常远高于影像附加参数x，可以先消去连接点坐标t，构建仅包含 附加模型参数x的改化法方程，如下式所示：

[M^TM-M^TW(W^TW)^-1W^TM]x＝M^TL-M^TW(W^TW)^-1W^TL(7)

式(7)可简化为：

Ax＝b(8)

其中，A为改化法方程系数矩阵，若区域网中共有m张像片，则A的阶数为6m×6m，x表示 像方附加参数的改正数，为待解求参数，b为一大小为6m的列向量。

2.基于共轭梯度的改化法方程快速解算

对上述改化法方程Ax＝b采用共轭梯度的方法迭代解求其最优解，共轭梯度方法如下：

改化法方程中A∈R对称正定，对于非零向量P₁,P₂若内积(AP₁,P₂)＝0，则称P₁和P₂对 于A是正交性的(或共轭的)。A正交性的向量组是线性无关向量组。

若以线性无关的A的正交性向量组构成R^6m×6m中的一组基(P₁,P₂…P_6m)，则Ax＝b的解x* 可表示为：

$x^{*}=\underset{i=1}{\overset{6m}{\Sigma }}{\alpha }_i{P}_i---\left(9\right)$

其中，α_i为基向量的系数。

代入Ax^*＝b得：

$\underset{i=1}{\overset{6m}{\Sigma }}{\alpha }_i{\mathrm{AP}}_i=b---\left(10\right)$

与P_j做内积得

$(\underset{i=1}{\overset{6m}{\Sigma }}{\alpha }_i{\mathrm{AP}}_i,P_j)=(b,P_i),j=1,2...6m---\left(11\right)$

$\underset{i=1}{\overset{6m}{\Sigma }}{\alpha }_i({\mathrm{AP}}_i,P_j)=(b,P_j)---\left(12\right)$

由P₁，P₂……P_6m，A正交，得：

α_j(AP_j,P_j)＝(b,P_j)(13)

${\alpha }_j=\frac{(b,P_j)}{({\mathrm{AP}}_j,P_j)},j=1,2...6m---\left(14\right)$

代入(9)则可以得到解向量x^*。

若计

$x^{\left(0\right)}=0,x^{\left(k\right)}=\underset{i=1}{\overset{k}{\Sigma }}{\alpha }_i{P}_i,k=1,2...6m---\left(15\right)$

则有

x^(k)＝x^(k-1)+α_kP_k(16)

因为k＝1,2...6m，设n为1,2...6m中某值。当k＝n时，有

$x^{\left(n\right)}=x^{(n-1)}+{\alpha }_n{P}_n=\underset{i=1}{\overset{n-1}{\Sigma }}{\alpha }_i{P}_i+{\alpha }_n{P}_n=\underset{i=1}{\overset{n}{\Sigma }}{\alpha }_i{P}_i=x^{*}---\left(17\right)$

即x⁽ⁿ⁾就是方程组Ax＝b的精确解，这说明按迭代(17)只需n步就得到方程组的精确解。

下面给出构造正交向量组P₁，P₂……P_n的方法：

设r₁，r₂…r_n∈Rⁿ是任意一组线性无关向量，则可利用下面的公式构造出A正交向量组P₁， P₂……P_n：

$P_1=r_1,P_k=r_k-\underset{i=1}{\overset{k-1}{\Sigma }}\frac{({\mathrm{AP}}_i,r_k)}{({\mathrm{AP}}_i,P_i)}{P}_i,k=2,3...n---\left(18\right)$

则有基于共扼梯度法的改化法方程算法过程：

1)令x⁽⁰⁾＝0，计算r₀＝b-Ax⁽⁰⁾，取P₀＝r₀，k＝0，给定误差限ε；

2)若||r_k||≤ε，则x^(k)，为近似解，结束计算。否则计算：

${\alpha }_k=\frac{(b,P_k)}{({\mathrm{AP}}_k,P_k)},x^{\left(k\right)}=x^{(k-1)}+{\alpha }_k{P}_k---\left(19\right)$

3)若k＝n-1，x^(k+1)＝x⁽ⁿ⁾就是Ax＝b的解，结束计算，否则计算：

$r_{k+1}=b-{\mathrm{Ax}}^{(k+1)},P_{k+1}=r_{k+1}-\underset{i=0}{\overset{k}{\Sigma }}\frac{({\mathrm{AP}}_i,r_{k+1})}{({\mathrm{AP}}_i,P_i)}{P}_i---\left(20\right)$

4)转向2)。

3.判断迭代过程是否收敛，是则进入步骤7，否则进入步骤4；

若像方附加参数的变化趋于稳定，可停止迭代。

4.根据改化法方程解算的参数更新附加像方参数

根据上述解算的改化法方程的解加上像方附加参数的当前值便可得到更新后的像方附加 参数的值：

X＝X⁰+dX(21)

其中X为平差待解求参数，X⁰为平差待解算参数的初值，dX为平差解算的改正数。

5.根据像方附加参数和影像RPC交会解算连接点物方坐标

基于附加参数的RPC模型，采用多片前方交会的方式平差解算连接点对应的物方坐标， 并将该坐标作为下次迭代解算的初值。多片前方交会的方式为现有技术，本发明不予赘述。 6.根据当前像方附加参数和连接点物方坐标更新改化法方程组，可视为返回步骤1建立了 新的法方程式，继续重复执行步骤2～5，直到在循环迭代的过程中若像方附加参数的变化趋 于稳定，停止迭代。

7.根据解算的像方附加参数进行RPC参数的精化，相应得到平差解算结果，所述精化具体 实现步骤为：

(1)生成地面虚拟控制点。在影像像平面上按预设间距均匀划分规则格网，对每个格网中心 像点p(samp,line)，利用该影像内定向参数，通过前方交会在物方局部一系列高程基准面(从 -1000m到9000m每1000m一个高程基准面)上交会得到一系列物方点P_i(i＝1,2,3...11)，此 时，像点p与物方点P_i构成11个虚拟控制点。

(2)以RPC模型为误差模型，RPC参数为待解算参数，在像方加入附加参数，根据原始RPC 列误差方程式，采用最小二乘的方法平差解算未知参数，得到精化后的RPC参数。

具体实施时，本发明所提供方法可基于软件技术实现自动运行流程，也可采用模块化方 式实现相应系统。本发明提供一种基于平差参数快速解算的超大规模区域网平差系统，包括 以下模块：

第一模块，用于建立改化法方程如下，

Ax＝b

其中，A为改化法方程系数矩阵，若区域网中共有m张像片，则A的阶数为6m×6m，x表示 像方附加参数的改正数，为待解求参数，b为一个大小为6m的列向量；

第二模块，用于基于共轭梯度快速解算第一模块所得改化法方程；

第三模块，用于判断迭代过程是否收敛，是则命令第七模块工作，否则命令第四模块工 作；

第四模块，用于根据改化法方程解算的参数更新附加像方参数如下，

X＝X⁰+dX

其中，X为平差待解求参数，X⁰为平差待解算参数的初值，dX为平差解算的改正数；

第五模块，用于根据像方附加参数和影像RPC参数交会解算连接点物方坐标；

第六模块，用于根据当前像方附加参数和连接点物方坐标更新改化法方程，命令第二模 块工作；

第七模块，用于根据解算的像方附加参数进行RPC参数的精化，得到平差结果。

各模块具体实现可参见相应步骤，本发明不予赘述。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技 术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不 会偏离本发明的精神或者超越所附权利要求书所定义的范围。