基于不完全分解预处理的共轭梯度法光束法平差方法与流程

本发明涉及摄影测量与遥感技术领域，尤其涉及一种基于不完全分解预处理的共轭梯度法光束法平差方法。

背景技术：

基于共线条件方程式的光束法平差解法，是一种把控制点的像点坐标、待定点的像点坐标以及其他内业、外业量测数据的一部分或全部均视作观测值，以整体且同时解求它们的最或是值和待定点空间坐标的解算方法。在航天航空摄影测量、低空摄影测量以及近景摄影测量等区域网平差工作中，一方面由于摄影方式呈现多样化的特点，另一方面也常常需要引入多种类型的附加参数(如定位和定姿设备的安装矩阵、自检校模型参数和附加几何约束条件等)，法方程系数阵的稀疏带状结构遭到了破坏，难以确定带宽且为加边矩阵。这给法方程求解快速算法的设计带来了显著的困难。

针对这一问题，国内外学者提出了多种解决方案。如McGlone C.(1996)考虑各种几何约束条件，构建统一模式的带附加参数的条件平差模型，区域网平差求解时采用整体广义逆解法；Wang X.和Clarke T.A.(2001)基于分组最小二乘估计方式，对像方参数和物方坐标进行相对独立估计的分组平差；Bartoli A.(2002)采用通常的代价函数，将不同类型的模型统一到同一框架下进行准线性优化；冯其强,等(2008)根据工业摄影测量对运算速度的需求，提出逐点修正未知近似值的自检较点松弛迭代解法；Ni K.,等(2007)将全局的非线性代价函数转化为各自坐标系下的子问题，尽可能利用外存空间和实现并行处理等。这些方法都在一定程度上兼顾了解算精度和运算效率的要求，但模型的设计和解算方法相对复杂，且不容易通过程序实现。

共轭梯度法是一种在可预测模式中不需要非零输入的，处理大规模稀疏系统的非常有用的迭代估计方法。该方法由最速下降法演化而来，最早由Hestenes M.R.和Stiefel E.(1952)提出，沿着梯度所构建的相互共轭的方向计算，从而最快速度得到目标函数的最小值。但由于对舍入误差的敏感性问题，算法的应用受到了极大的限制。近些年来，该问题通过与各种预处理算法相结合而得到解决。由于算法简便、存储需求小等优势，共轭梯度法已经成为求解大规模、稀疏线性方程组最常用的一类方法，并且开始引入到测量平差领域。

因而，利用共轭梯度法在正交基下的不变性进行基于共轭梯度法的光束法平差算法设计，并采用基于Cholesky分解的预处理方法减低法方程系数阵的条件数，能够有效地提高区域网平差的稳定性和计算效率，更好地适应大区域影像产品生产和三维重建等工作。

公开内容

(一)要解决的技术问题

本公开提供了一种基于不完全分解预处理的共轭梯度法光束法平差方法，以至少部分解决以上所提出的技术问题。

(二)技术方案

根据本公开的一个方面，提供了一种基于不完全分解预处理的共轭梯度法光束法平差方法，包括：步骤S1，构建光束法平差的法方程的统一形式，并构建采用多类未知数交替趋近解法的多类法方程；步骤S2，采用对角线矩阵块构造与光束法平差法方程系数阵近似的系数阵形式，将其作为预处理矩阵进行不完全Cholesky分解预处理，构造光束法平差新的法方程结构，降低法方程系数阵的条件数；步骤S3，给定共轭梯度法光束法平差的迭代初值，依次沿着相互共轭的方向进行相关目标函数解的搜索，并将目标函数改造为能够直接求得原法方程的解X的值的形式；步骤S4，给定一小值ε作为迭代收敛条件进行迭代计算，直至满足收敛条件迭代终止，输出基于不完全分解预处理的共轭梯度法光束法平差法方程的最优解。

在本公开一些实施例中，所述步骤S1包括：子步骤S101，获取区域内具有重叠度的多景光学影像和对应的成像外方位参数，包括：与成像姿态相关的3个外方位角元素：俯仰角，翻滚角和偏航角，与成像位置相关的3个外方位线元素：像主点相对于摄影测量坐标系在X，Y和Z方向的坐标；子步骤S102，通过多景光学影像和对应的成像参数，利用共线方程和最小二乘平差方法，构建光束法平差法方程的统一形式：

NX＝W

其中，X为法方程的未知数矩阵；N为法方程的系数矩阵，满足n阶对阵正定矩阵；W为法方程的常数项矩阵；

子步骤S103，构造采用多类未知数交替趋近解法的n类未知数的法方程式：

其中，n≥2，所述n类未知数包括：影像外方位元素的改正数、物方点空间坐标的改正数，Xn为第n类法方程的未知数矩阵；Nnn为第n类法方程的系数矩阵，满足n阶对阵正定矩阵；Wn为第n类法方程的常数项矩阵。

在本公开一些实施例中，所述n类未知数包括：影像外方位元素的改正数、物方点空间坐标的改正数，以及各非摄影测量附件参数的改正数。

在本公开一些实施例中，所述步骤S2包括：子步骤S201，采用对角线矩阵块构造与光束法平差法方程系数阵近似的系数阵形式，子步骤S103中多类未知数交替趋近解法，则交替一步的过程用如下方程组作为近似表示：

其中，M为子步骤S103中n类未知数的法方程式的系数矩阵的对角阵；

子步骤S202，将子步骤S201中的对角线矩阵块构造的系数阵作为共轭梯度法光束法平差的预处理矩阵，进行不完全Cholesky分解预处理，得到：

M＝LL^T

其中，L为一个下三角矩阵；

子步骤S203，对系数阵M不完全Cholesky分解后，将子步骤S102中的法方程转换为如下形式：

[L^-1N(L^-1)^T]L^TX＝L^-1W

其中，法方程的未知数X替换为Y＝L^TX；法方程的系数阵N替换为F＝L^-1N(L^-1)^T，仍满足正定对称矩阵条件；法方程的常数项W替换为G＝L^-1W，将上式简化为如下形式：

FY＝G

经过变换，将子步骤S102中共轭梯度法光束法平差原法方程的求解问题转化为对新法方程的求解问题，其中，系数阵F的条件数远低于原法方程系数阵N的条件数。

在本公开一些实施例中，所述步骤S3包括：子步骤S301，给定子步骤S102中法方程的初值为X⁰，将预处理后法方程的初值替换为Y⁰＝L^TX⁰，设定共轭梯度法光束法平差法的迭代第一步为k＝0；子步骤S302：在n维向量空间中，共轭梯度法光束法平差原法方程依次沿着相互共轭的方向X1,X2,…Xn，进行相关目标函数解的搜索；当搜索至第k步时的解为X^k＝α1X1+α2X2+...+αkXk，其中，相关目标函数的解满足Y^k＝L^TX^k；子步骤S303：为使共轭梯度法光束法平差能够直接求得原法方程的解X的值，将目标函数中的X^k替换为(L^-1)^TY^k。

在本公开一些实施例中，所述步骤S302包括：对于新的法方程，引入残差变量R^k：

R^k＝G-FY^k＝L^-1[W-N(L^-1)^TL^TX^k]＝L^-1r^k

其中，Y^k＝L^TX^k，r^k是针对原法方程引入的一个残差变量，满足：r^k＝W-N(L^-1)^TL^TX^k；再计算残差变量R^k的初值：

R⁰＝L^-1[W-N(L^-1)^TL^TX⁰]＝L^-1r⁰。

在本公开一些实施例中，所述步骤S303包括：引入一个变量Z^k：

Z^k＝(L^-1)^TL^-1r^k＝M^-1r^k

解算方程MZ⁰＝R⁰得到变量Z^k的初值Z⁰，并令X1＝Z⁰。

在本公开一些实施例中，所述步骤S4包括：子步骤S401：给定一小值，设定迭代收敛条件对平差是否继续迭代进行判断，若满足收敛条件则转至子步骤S404，若不满足则转至子步骤S402至子步骤S403；子步骤S402：计算迭代至该步时与目标函数取值相关的参数α、方程组的搜索结果X、残差变量R和参数Z的取值；子步骤S403：计算迭代至该步时与目标函数取值相关的参数β以及下一步的迭代方向；子步骤S404：迭代终止，输出目标函数取得最小值时X^k的值，作为共轭梯度法光束法平差法方程的最优解。

在本公开一些实施例中，所述子步骤S401包括：给定一小值ε，进行判断，若||R^k||≥ε则执行子步骤S402至子步骤S403，并进行迭代；否则迭代中止，执行子步骤S404。

在本公开一些实施例中，所述子步骤S402包括：令k＝k+1，在第k步沿着方向Xk进行搜索，Xk同时满足与全部的X1,X2,…,Xk-1均共轭，计算方程组解的搜索结果：

X^k＝X^k-1+αkXk

其中，与目标函数取值相关的参数αk由以下公式计算得出：

再计算变量Z^k：

Z^k＝(L^-1)^TL^-1R^k＝M^-1R^k，

其中，残差变量R^k＝R^k-1-αkNXk；

所述子步骤S403包括：计算第k+1步的搜索方向Xk+1：

Xk+1＝Z^k+βkXk

其中，参数βk满足：

(三)有益效果

从上述技术方案可以看出，本公开基于不完全分解预处理的共轭梯度法光束法平差方法至少具有以下有益效果其中之一：

(1)采用多类未知数交替趋近解法生成与光束法平差法方程系数阵近似的预处理矩阵，该矩阵结构简单且满足正定对称条件，便于对法方程的求解进行优化；

(2)采用不完全Cholesky分解预处理方法对预处理矩阵进行分解，构造新的法方程结构，能够显著降低法方程系数阵的条件数，减弱方程组的病态性，并提高平差算法的稳定性；

(3)相比于传统的光束法平差方法，在保证平差结果精度的前提下，具有算法设计和实现简单，空间和时间复杂度低，并便于实现并行处理等优势，适合解决大规模的区域网平差问题；

(4)相比于其他共轭梯度法光束法平差方法，具有模型参数间的相关性更低，平差迭代的收敛速率更快和算法执行效率更高等优势，能有效解决含各类附加条件的复杂区域网平差问题。

附图说明

图1为本公开第一实施例基于不完全分解预处理的共轭梯度法光束法平差方法的流程图。

具体实施方式

本公开提供了基于不完全分解预处理的共轭梯度法光束法平差方法，首先采用多类未知数交替趋近解法，构造与光束法平差法方程系数阵近似且结构更加简单的预处理矩阵；然后采用不完全Cholesky分解预处理方法对预处理矩阵进行分解，构造新的法方程结构；接下来采用共轭梯度法，构造新的目标函数和n维向量空间中关于法方程系数阵的线性无关的向量组。沿着相互共轭的方向进行方程组解的搜索；最后采用最小二乘平差迭代计算，求得共轭梯度法光束法平差法方程的最优解。

为使本公开的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本公开进一步详细说明。

本公开某些实施例于后方将参照所附附图做更全面性地描述，其中一些但并非全部的实施例将被示出。实际上，本公开的各种实施例可以许多不同形式实现，而不应被解释为限于此数所阐述的实施例；相对地，提供这些实施例使得本公开满足适用的法律要求。另外，虽然本文可提供包含特定值的参数的示范，但应了解，参数无需确切等于相应的值，而是可在可接受的误差容限或设计约束内近似于相应的值。

在本公开的第一个示例性实施例中，提供了一种基于不完全分解预处理的共轭梯度法光束法平差方法。图1为本公开第一实施例基于不完全分解预处理的共轭梯度法光束法平差方法的流程图。如图1所示，本公开基于不完全分解预处理的共轭梯度法光束法平差方法包括：

步骤S1，获取区域内具有重叠度的多景光学影响和对应的成像参数，构建光束法平差的法方程的统一形式，及采用多类未知数交替趋近解法的多类法方程；

步骤S2，采用对角线矩阵块构造与光束法平差法方程系数阵近似的系数阵形式，并将其作为预处理矩阵进行不完全Cholesky分解预处理；再采用分解后的系数阵构造光束法平差新的法方程结构，降低法方程系数阵的条件数；

步骤S3，给定共轭梯度法光束法平差的迭代初值，依次沿着相互共轭的方向进行相关目标函数解的搜索，构造残差变量R的形式并计算其初值；将目标函数改造为可以直接求得原法方程的解X的值的形式，给定新的参数Z并计算其初值；

步骤S4，给定一小值ε作为迭代收敛条件进行迭代计算，直至迭代终止，输出基于不完全分解预处理的共轭梯度法光束法平差法方程的最优解。

以下分别对本实施例的各个步骤进行详细描述。

步骤S1，获取区域内具有重叠度的多景光学影响和对应的成像参数，构建光束法平差的法方程的统一形式，及采用多类未知数交替趋近解法的多类法方程。具体包括：

子步骤S101，获取区域内具有重叠度的多景光学影像和对应的成像外方位参数，包括：与成像姿态相关的3个外方位角元素(俯仰角pitch，翻滚角roll和偏航角yaw),与成像位置相关的3个外方位线元素(像主点相对于摄影测量坐标系在X，Y和Z方向的坐标Xs,Ys和Zs)。

子步骤S102，构建光束法平差法方程的统一形式，包括:

通过多景光学影像和对应的成像参数构建光束法平差的法方程，其形式为：

NX＝W

其中，X为法方程的未知数矩阵；N为法方程的系数矩阵，满足n阶对阵正定矩阵；W为法方程的常数项矩阵。

若给定法方程未知数矩阵X的一组初值为X0，则法方程可转换为如下形式：

N(X-X0)＝W-NX0

为了简便表示，这里仍分别将X-X0和W-NX0记为X和W。则法方程的形式保持不变。

子步骤S103，构造采用多类未知数交替趋近解法的多类法方程式。光束法平差法方程的各类型未知数(包括影像外方位元素的改正数、物方点空间坐标的改正数，以及各非摄影测量附件参数的改正数等)的近似值可以采用多类未知数交替趋近解算得到，例如，当光束法平差法方程的未知数为三类(各景影像的外方位元素，所有物方点的空间坐标和各非摄影测量附加参数)时，即X＝[X1 X2 X3]。

三类未知数交替趋近解法的一次迭代过程包括：

第一步，假定已知所有物方点空间坐标和各非摄影测量附加参数，通过共线方程和最小二乘平差方法求解各景影像的外方位元素，即法方程的未知数矩阵仅包括各景影像的6个外方位元素。

第二步，假定已知各景影像的外方位元素和各非摄影测量附加参数，通过共线方程和最小二乘平差方法求解所有物方点空间坐标，即法方程的未知数矩阵仅包括各物方点的三维空间坐标。

第三步，假定已知各景影像的外方位元素和各物方点的空间坐标，通过共线方程和最小二乘平差方法求解各非摄影测量附加参数，即法方程的未知数矩阵仅包括各非摄影测量附加参数。

三类未知数交替趋近解法相当于依次迭代解算如下三类法方程式：

步骤S2，采用对角线矩阵块构造与光束法平差法方程系数阵近似的系数阵形式，并将其作为预处理矩阵进行不完全Cholesky分解预处理；再采用分解后的系数阵构造光束法平差新的法方程结构，降低法方程系数阵的条件数。具体包括：

子步骤S201，采用对角线矩阵块构造与光束法平差法方程系数阵近似的系数阵形式。利用数学分析中连续函数的性质，如果所给定的各个参数的初值比较接近真值，那么子步骤S103中多类未知数交替趋近解法，交替一步的过程可以用如下方程组作为近似表示：

而此时，子步骤S102中采用光束法区域网平严密解法的法方程结构为：

由“后方交会-前方交会”迭代交替解法的有效性表明其系数阵M是光束法平差中法方程系数阵N的近似值，即M≈N，因此可以采用对角线矩阵块构造与光束法平差法方程系数阵近似的系数阵形式。

子步骤S202，将子步骤S201中的对角线矩阵块构造的系数阵作为共轭梯度法光束法平差的预处理矩阵，进行不完全Cholesky分解预处理。由于子步骤S201中迭代交替解法的系数阵是由子步骤S102中光束法平差法方程系数阵N的对角块矩阵组成的，与系数阵N的值近似，并且系数阵M同样为正定对阵矩阵，且具有更加简洁的矩阵结构，因此可将系数阵M作为共轭梯度法光束法平差的预处理矩阵。

对系数阵M作不完全Cholesky分解预处理可得：

M＝LL^T

其中，L为一个下三角矩阵。

子步骤S203，采用分解后的系数阵构造光束法平差新的法方程结构，降低法方程系数阵的条件数，包括：

对系数阵M不完全Cholesky分解后，将子步骤S102中的法方程转换为如下形式：

[L^-1N(L^-1)^T]L^TX＝L^-1W

其中，法方程的未知数X替换为Y＝L^TX；法方程的系数阵N替换为F＝L^-1N(L^-1)^T，仍满足正定对称矩阵条件；法方程的常数项W替换为G＝L^-1W。上式简化为如下形式：

FY＝G

经过上述变换，子步骤S102中共轭梯度法光束法平差原法方程的求解问题转化为对新法方程的求解问题。由于系数阵F的条件数远低于原法方程系数阵N的条件数，因此，上式可实现基于共轭梯度法光束法平差的法方程求解的快速迭代收敛。

步骤S3，给定共轭梯度法光束法平差的迭代初值，依次沿着相互共轭的方向进行相关目标函数解的搜索，构造残差变量R的形式并计算其初值；将目标函数改造为可以直接求得原法方程的解X的值的形式，给定新的参数Z并计算其初值。具体包括：

子步骤S301，给定子步骤S102中法方程的初值为X⁰，则预处理后法方程的初值替换为Y⁰＝L^TX⁰。设定共轭梯度法光束法平差法的迭代第一步为k＝0。

子步骤S302：在n维向量空间中，共轭梯度法光束法平差原法方程依次沿着相互共轭的方向X1,X2,…Xn，进行相关目标函数解的搜索。当搜索至第k步时的解为X^k＝α1X1+α2X2+...+αkXk。

针对原法方程，假定引入的一个残差变量r^k满足如下关系：

r^k＝W-N(L^-1)^TL^TX^k

而对于新的法方程，其相关目标函数的解满足Y^k＝L^TX^k。引入残差变量R^k满足如下关系：

R^k＝G-FY^k＝L^-1[W-N(L^-1)^TL^TX^k]＝L^-1r^k

计算残差变量R^k的初值为R⁰＝L^-1[W-N(L^-1)^TL^TX⁰]＝L^-1r⁰

子步骤S303：为使共轭梯度法光束法平差能够直接求得原法方程的解X的值，这里将目标函数中的X^k替换为(L^-1)^TY^k。

引入一个变量Z^k，满足如下关系：

Z^k＝(L^-1)^TL^-1r^k＝M^-1r^k

解算方程MZ⁰＝R⁰得到变量Z^k的初值Z⁰，并令X1＝Z⁰。

步骤S4，给定一小值ε作为迭代收敛条件进行迭代计算，直至迭代终止，输出基于不完全分解预处理的共轭梯度法光束法平差法方程的最优解。具体包括：

子步骤S401：给定一小值，设定迭代收敛条件对平差是否继续迭代进行判断，具体包括：

给定一小值ε，进行判断，若||R^k||≥ε则执行子步骤S402至子步骤S403，并进行迭代。否则迭代中止，执行子步骤S404。

子步骤S402：计算迭代至该步时，与目标函数取值相关的参数α，方程组的搜索结果X，残差变量R和参数Z的取值，包括：

令k＝k+1，计算与目标函数取值相关的参数αk，满足：

第k步应该沿着方向Xk进行搜索，Xk同时满足与全部的X1,X2,…,Xk-1均共轭，计算方程组解的搜索结果X^k，满足：

X^k＝X^k-1+αkXk

并且，残差变量R^k满足如下关系：

R^k＝R^k-1-αkNXk

计算变量Z^k＝(L^-1)^TL^-1R^k＝M^-1R^k。

子步骤S403：计算迭代至该步时，与目标函数取值相关的参数β以及下一步的迭代方向，包括：

计算第k+1步的搜索方向Xk+1，满足如下关系：

Xk+1＝Z^k+βkXk

由于Xk和Xk+1共轭，满足条件则其中参数βk满足：

子步骤S404：迭代终止，输出目标函数取得最小值时X^k的值，作为共轭梯度法光束法平差法方程的最优解。

综上所述，本发明基于不完全分解预处理的共轭梯度法光束法平差方法主要采用多类未知数交替趋近解法，构造与光束法平差法方程系数阵近似且结构更加简单的预处理矩阵，并采用不完全Cholesky分解预处理方法对预处理矩阵进行分解，构造新的法方程结构。经过上述优化后，采用共轭梯度法构造新的目标函数和n维向量空间中关于法方程系数阵的线性无关的向量组，沿着相互共轭的方向进行方程组解的搜索，并采用最小二乘平差迭代计算，求得共轭梯度法光束法平差法方程的最优解。

至此，已经结合附图对本实施例进行了详细描述。依据以上描述，本领域技术人员应当对本发明基于不完全分解预处理的共轭梯度法光束法平差方法有了清楚的认识。

至此，本公开第一实施例基于不完全分解预处理的共轭梯度法光束法平差方法介绍完毕。

至此，已经结合附图对本公开实施例进行了详细描述。需要说明的是，在附图或说明书正文中，未绘示或描述的实现方式，均为所属技术领域中普通技术人员所知的形式，并未进行详细说明。此外，上述对各元件和方法的定义并不仅限于实施例中提到的各种具体结构、形状或方式，本领域普通技术人员可对其进行简单地更改或替换。

再者，单词“包含”不排除存在未列在权利要求中的元件或步骤。位于元件之前的单词“一”或“一个”不排除存在多个这样的元件。

此外，除非特别描述或必须依序发生的步骤，上述步骤的顺序并无限制于以上所列，且可根据所需设计而变化或重新安排。并且上述实施例可基于设计及可靠度的考虑，彼此混合搭配使用或与其他实施例混合搭配使用，即不同实施例中的技术特征可以自由组合形成更多的实施例。

在此提供的算法和显示不与任何特定计算机、虚拟系统或者其它设备固有相关。各种通用系统也可以与基于在此的示教一起使用。根据上面的描述，构造这类系统所要求的结构是显而易见的。此外，本公开也不针对任何特定编程语言。应当明白，可以利用各种编程语言实现在此描述的本公开的内容，并且上面对特定语言所做的描述是为了披露本公开的最佳实施方式。

本公开可以借助于包括有若干不同元件的硬件以及借助于适当编程的计算机来实现。本公开的各个部件实施例可以以硬件实现，或者以在一个或者多个处理器上运行的软件模块实现，或者以它们的组合实现。本领域的技术人员应当理解，可以在实践中使用微处理器或者数字信号处理器(DSP)来实现根据本公开实施例的相关设备中的一些或者全部部件的一些或者全部功能。本公开还可以实现为用于执行这里所描述的方法的一部分或者全部的设备或者装置程序(例如，计算机程序和计算机程序产品)。这样的实现本公开的程序可以存储在计算机可读介质上，或者可以具有一个或者多个信号的形式。这样的信号可以从因特网网站上下载得到，或者在载体信号上提供，或者以任何其他形式提供。

本领域那些技术人员可以理解，可以对实施例中的设备中的模块进行自适应性地改变并且把它们设置在与该实施例不同的一个或多个设备中。可以把实施例中的模块或单元或组件组合成一个模块或单元或组件，以及此外可以把它们分成多个子模块或子单元或子组件。除了这样的特征和/或过程或者单元中的至少一些是相互排斥之外，可以采用任何组合对本说明书(包括伴随的权利要求、摘要和附图)中公开的所有特征以及如此公开的任何方法或者设备的所有过程或单元进行组合。除非另外明确陈述，本说明书(包括伴随的权利要求、摘要和附图)中公开的每个特征可以由提供相同、等同或相似目的的替代特征来代替。并且，在列举了若干装置的单元权利要求中，这些装置中的若干个可以是通过同一个硬件项来具体体现。

类似地，应当理解，为了精简本公开并帮助理解各个公开方面中的一个或多个，在上面对本公开的示例性实施例的描述中，本公开的各个特征有时被一起分组到单个实施例、图、或者对其的描述中。然而，并不应将该公开的方法解释成反映如下意图：即所要求保护的本公开要求比在每个权利要求中所明确记载的特征更多的特征。更确切地说，如下面的权利要求书所反映的那样，公开方面在于少于前面公开的单个实施例的所有特征。因此，遵循具体实施方式的权利要求书由此明确地并入该具体实施方式，其中每个权利要求本身都作为本公开的单独实施例。

以上所述的具体实施例，对本公开的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本公开的具体实施例而已，并不用于限制本公开，凡在本公开的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本公开的保护范围之内。