一种基于不精确牛顿解的预条件共轭梯度区域网平差方法

摘要

本发明公开了一种基于不精确牛顿解的预条件共轭梯度区域网平差方法，在传统的区域网平差流程中引入了预条件共轭梯度法求解大规模法方程，避免了对法方程的存储以及直接求逆运算，且本方法不直接存储法方程系数矩阵B，而是逐点计算矩阵‑向量积法方程系数矩阵Bd，并使得区域网平差流程适用于并行化设计方案，同时采用不精确牛顿解替代法方程的精确解，从而减少预条件共轭梯度法迭代次数，减少了数据的存储量,提高整个区域网平差的求解速度。

权利要求

1.一种基于不精确牛顿解的预条件共轭梯度区域网平差方法，其特征在于，包括：S1、导入区域网平差计算需要的原始数据，且将导入的原始数据进行时空基准统一，得到初始数据，所述原始数据至少包括初始内外方位元素数据以及点位数据；

S2、进入区域网平差迭代流程：分别读取每一个像点对应的初始数据，计算每一个像点对应的法方程常数项向量c的一个分量，并将各个法方程常数项向量c的分量累加，得到完整的法方程常数项向量c作为步骤S4中残差向量s的初值；

S3、分别读取每一个像点对应的初始数据，计算每一个像点对应的预条件矩阵M的分-1量，并将该分量求逆并乘以方程常数项向量c得到矩阵-向量积M c的一个分量，由将各个矩阵-向量积M-1c的分量累加，得到完整的矩阵-向量积M-1c作为步骤S4中方向向量d的初值；

S4、进入预条件共轭梯度法迭代流程：对未知数改正数向量u、残差向量s以及方向向量d进行初始化；

S5、分别读取每一个像点对应的初始数据，计算每一个像点对应的法方程系数矩阵B的与方向向量d的矩阵-向量积Bd分量，将各矩阵-向量积Bd分量累加，得到完整的矩阵向量积Bd；

S6、根据预条件共轭梯度法，以及本次迭代中的未知数改正数向量u，残差向量s，方向向量d，预条件矩阵M以及法方程系数矩阵B，计算新的未知数改正数向量u、新的残差向量s以及新的方向向量d；

S7、根据不精确牛顿解法，计算本次迭代的强制序列系数η；

S8、判断此次预条件共轭梯度迭代是否符合预定收敛条件，若符合，则结束预条件共轭梯度迭代，输出新的未知数改正数向量u，执行步骤S9；否则，返回步骤S5；

S9、判断此次区域网平差迭代是否符合预定收敛条件，若符合，则结束区域网平差迭代，执行步骤S10，否则，返回步骤S2；

S10、根据步骤S8中输出的新的未知数改正数向量u，更新所有的未知数数值，并输出所有的未知数数值。

2.如权利要求1所述的基于不精确牛顿解的预条件共轭梯度区域网平差方法，其特征在于，所述步骤S3中的预条件矩阵M为Jacobi预条件矩阵。

3.如权利要求1所述的基于不精确牛顿解的预条件共轭梯度区域网平差方法，其特征在于，所述步骤S8具体包括：根据步骤S7中计算得到的强制序列系数η，若强制序列系数η小于第一给定阈值或者迭代次数大于第二给定阈值，则结束预条件共轭梯度迭代，输出新的未知数改正数向量u，执行步骤S9，否则，返回步骤S5。

4.如权利要求1所述的基于不精确牛顿解的预条件共轭梯度区域网平差方法，其特征在于，所述步骤S9具体包括：统计步骤S8中计算得到的新的未知数改正数向量u中所有元素的绝对值最大值，若该绝对值最大值小于第三给定阈值或者迭代次数大于第四给定阈值，则结束区域网平差迭代，执行步骤S10，否则，返回步骤S2。

5.如权利要求1所述的基于不精确牛顿解的预条件共轭梯度区域网平差方法，其特征在于，所述步骤S7具体为：根据不精确牛顿解法，将步骤S2计算出的法方程常数项向量c的绝对值和步骤S6计算出的新的残差向量s的绝对值进行相除，得到本次迭代的强制序列系数η。

说明书

一种基于不精确牛顿解的预条件共轭梯度区域网平差方法

技术领域

[0001] 本发明涉及测绘科学与技术领域，具体涉及一种基于不精确牛顿解的预条件共轭梯度区域网平差方法，主要应用于超大规模测区摄影测量等领域。

背景技术

[0002] 区域网平差技术经过几十年的发展，其方法和流程已经相对成熟，并且在测绘领域得到了广泛的应用。然而随着科技的加速进步，新的传感器不断涌现，如航天领域的高分辨率卫星，立体测绘卫星，航空领域的规则航空摄影系统，倾斜航空摄影系统，无人机、飞艇摄影系统等。同时，全世界范围内三维建模应用需求不断增加，也使得大量地面车载摄影系统，近景摄影系统，普通数码相机，甚至是智能手机，网络图片库等采集的影像被用于三维建模。影像数据源越来越丰富的同时，其分辨率也不断提高，以前数十米的卫星影像如今可以达到最高0.35米(WorldView-3)，航空影像的分辨率更是进入了厘米级时代。分辨率的增加必然会带来数据量的增大，摄影时的航线设计也不再满足传统的条带式规则分布，给相应的数据处理方法带来了一定挑战，受法方程大小的限制，传统的区域网平差技术流程已经不能满足大规模数据处理需求。特别是当测区数据大小超过1万张影像时，传统的区域网平差方法对内存的需求急剧增大，即使有少部分图形工作站内存容量足够大，但大量的内存占用使得计算效率大幅降低，同时，数据处理的硬件成本也随之增大，上述问题均阻碍了各类新型传感器数据的广泛应用。

发明内容

[0003] 为了解决上述问题，本发明在区域网平差中，引入预条件共轭梯度法求解大规模线性方程组(法方程)，避免存储大规模法方程系数矩阵，使得对大规模测区(1万至10万张影像，下文中统称此类测区为大规模测区)数据的区域网平差成为可能，其流程也更加易于并行化设计。预条件共轭梯度法是一个迭代求解的过程，每一次迭代都需要遍历所有的像点观测值数据，每一次区域网平差迭代内部又包括N次预条件共轭梯度法迭代，N是预条件共轭梯度法迭代收敛次数，因而计算量较传统的摄影测量区域网平差流程要大很多，为了减少预条件共轭梯度法的迭代次数，本发明采用一种不精确牛顿解取代法方程的精确解，使得预条件共轭梯度法可以快速的收敛，从而提高整个区域网平差的求解速度。

[0004] 本发明所要解决的技术问题是提供一种基于不精确牛顿解的预条件共轭梯度区域网平差方法，能够解决现有技术的不足。

[0005] 本发明解决上述技术问题的技术方案如下：

[0006] 一种基于不精确牛顿解的预条件共轭梯度区域网平差方法，包括：

[0007] S1、导入区域网平差计算需要的原始数据，且将导入的原始数据进行时空基准统一，得到初始数据，所述原始数据至少包括初始内外方位元素数据以及点位数据；

[0008] S2、进入区域网平差迭代流程：分别读取每一个像点对应的初始数据，计算每一个像点对应的法方程常数项向量c的一个分量，并将各个法方程常数项向量c的分量累加，得到完整的法方程常数项向量c作为步骤S4中残差向量s的初值；

[0009] S3、分别读取每一个像点对应的初始数据，计算每一个像点对应的预条件矩阵M的分量，并将该分量求逆并乘以方程常数项向量c得到矩阵-向量积M-1c的一个分量，将各个矩阵-向量积M-1c的分量累加，得到完整的矩阵-向量积M-1c作为步骤S4中方向向量d的初值；

[0010] S4、进入预条件共轭梯度法迭代流程：对未知数改正数向量u、残差向量s以及方向向量d进行初始化；

[0011] S5、分别读取每一个像点对应的初始数据，计算每一个像点对应的法方程系数矩阵B的与方向向量d的矩阵-向量积Bd分量，将各矩阵-向量积Bd分量累加，得到完整的矩阵向量积Bd；

[0012] S6、根据预条件共轭梯度算法，以及本次迭代中的未知数改正数向量u，残差向量s，方向向量d，预条件矩阵M以及法方程系数矩阵B，计算新的未知数改正数向量u、新的残差向量s以及新的方向向量d；

[0013] S7、根据不精确牛顿解法，计算本次迭代的强制序列系数η；

[0014] S8、判断此次预条件共轭梯度迭代是否符合预定收敛条件，若符合，则结束预条件共轭梯度迭代，输出新的未知数改正数向量u，执行步骤S9；否则，返回步骤S5；

[0015] S9、判断此次区域网平差迭代是否符合预定收敛条件，若符合，则结束区域网平差迭代，执行步骤S10，否则，返回步骤S2；

[0016] S10、根据步骤S8中输出的新的未知数改正数向量u，更新所有的未知数数值，并输出所有的未知数数值。

[0017] 本发明的有益效果是：本发明提供的基于不精确牛顿解的预条件共轭梯度区域网平差方法，在传统的区域网平差流程中引入了预条件共轭梯度法求解大规模法方程，避免了对法方程的存储以及直接求逆运算，且本方法不直接存储法方程系数矩阵B，而是逐点计算矩阵-向量积法方程系数矩阵Bd，减少了数据的存储量，同时采用不精确牛顿解替代法方程的精确解，从而减少预条件共轭梯度法迭代次数，提高整个区域网平差的求解速度。本发明特别适用于超大规模测区(10万张影像以上，如全省、全国乃至全球作为一整个测区)的区域网平差计算。

附图说明

[0018] 图1为本发明实施例一的基于不精确牛顿解的预条件共轭梯度区域网平差方法流程图。

具体实施方式

[0019] 以下结合附图对本发明的原理和特征进行描述，所举实例只用于解释本发明，并非用于限定本发明的范围。

[0020] 实施例一、一种基于不精确牛顿解的预条件共轭梯度区域网平差方法。下面结合图1对本实施例提供的方法进行说明。

[0021] 参见图1，S1、导入区域网平差计算需要的原始数据，且将导入的原始数据进行时空基准统一，得到初始数据，所述原始数据至少包括初始内外方位元素数据以及点位数据；

[0022] S2、进入区域网平差迭代流程：分别读取每一个像点对应的初始数据，计算每一个像点对应的法方程常数项向量c的一个分量，并将各个法方程常数项向量c的分量累加，得到完整的法方程常数项向量c作为步骤S4中残差向量s的初值；

[0023] 具体的，首先导入区域网平差计算需要的原始数据，原始数据主要包括初始内外方位元素数据、连接点数据、控制点数据以及检查点数据(有些时候不存在控制点数据和检查点数据)。

[0024] 对保存的原始数据进行预处理，具体包括对导入的原始数据进行时空基准统一，包括事件基准统一和空间基准统一，若原始数据中存在控制点数据，还应该将连接点物方坐标以及初始外方位元素数据的坐标转换至控制点数据所在的坐标系内，以及计算并统计区域网平差计算过程中未知数分组的情况，未知数的个数以及需要用到的预条件矩阵M的类型和大小等信息。本实施例中的预条件矩阵与法方程系数矩阵大小一致，但构造简单，更易于求逆，使用其逆矩阵左乘法方程系数矩阵后，使得法方程系数矩阵的条件数减少，从而加快共轭梯度法求解法方程的迭代收敛速度，并保持收敛的稳健性，常用的预条件矩阵有Jacobi预条件矩阵。

[0025] S3、分别读取每一个像点对应的初始数据，计算每一个像点对应的预条件矩阵M的分量，并将该分量求逆并乘以方程常数项向量c得到矩阵-向量积M-1c的一个分量，将各个矩-1 -1阵-向量积M c的分量累加，得到完整的矩阵-向量积M c作为步骤S4中方向向量d的初值；

[0026] 具体的，本实施例中分别读取一个或多个像点对应的初始数据，各线程并行计算得到对应的法方程常数项向量c的一个分量，将各个法方程常数项向量c的分量累加，得到完整的法方程常数项向量c；以及每一个线程分别读取一个或多个像点对应的初始数据，各线程并行计算预条件矩阵M对应的分量，并将该分量求逆乘以法方程常数项向量c得到矩阵-向量积M-1c的一个分量，将各个矩阵-向量积M-1c的分量累加，得到完整的矩阵-向量积M-1c。

[0027] S4、进入预条件共轭梯度法迭代流程：对未知数改正数向量u、残差向量s以及方向向量d进行初始化；

[0028] 具体的，在本实施例中，将未知数改正数向量u的初始值设为0，将法方程常数项向量c作为预条件共轭梯度法迭代流程中的残差向量s的初始值，以及将矩阵-向量积M-1c作为预条件共轭梯度法迭代流程中的残差向量d的初始值，至此进入预条件共轭梯度法迭代流程。本实施例在区域网平差计算过程中引入预条件共轭梯度法，设置预条件矩阵，减少法方程系数矩阵的条件数，加快迭代收敛速度，共轭梯度法的特征是无需对法方程系数矩阵进行求逆，通过迭代搜索线性方程组的最优解，且每次迭代搜索方向之间是相互共轭的，具有存储量少，计算方便，收敛快等特点。

[0029] S5、分别读取每一个像点对应的初始数据，计算每一个像点对应的法方程系数矩阵B的与方向向量d的矩阵-向量积Bd分量，将各矩阵-向量积Bd分量累加，得到完整的矩阵向量积Bd；

[0030] S6、根据预条件共轭梯度算法，以及本次迭代中的未知数改正数向量u，残差向量s，方向向量d，预条件矩阵M以及法方程系数矩阵B，计算新的未知数改正数向量u、新的残差向量s以及新的方向向量d；

[0031] S7、根据不精确牛顿解法，计算本次迭代的强制序列系数η；

[0032] 具体的，是根据不精确牛顿解法，将步骤S2计算出的法方程常数项向量c的绝对值和步骤S6计算出的新的残差向量s的绝对值进行相除，得到本次迭代的强制序列系数η。

[0033] S8、判断此次预条件共轭梯度迭代是否符合预定收敛条件，若符合，则结束预条件共轭梯度迭代，输出新的未知数改正数向量u，执行步骤S9；否则，返回步骤S5；

[0034] 具体的，根据步骤S7中计算得到的强制序列系数η，若强制序列系数η小于第一给定阈值或者迭代次数大于第二给定阈值，则结束预条件共轭梯度迭代，输出新的未知数改正数向量u，执行步骤S9，否则，返回步骤S5。

[0035] S9、判断此次区域网平差迭代是否符合预定收敛条件，若符合，则结束区域网平差迭代，执行步骤S10，否则，返回步骤S2；

[0036] 具体的，统计步骤S8中输出的新的未知数改正数向量u中所有元素的绝对值最大值，若该绝对值最大值小于第三给定阈值或者区域网平差迭代次数大于第四给定阈值，则结束区域网平差迭代，执行步骤S10，否则，返回步骤S2。

[0037] S10、根据步骤S8中输出的新的未知数改正数向量u，更新所有的未知数数值，并输出所有的未知数数值。

[0038] 本方法涉及的利用预条件共轭梯度法迭代的区域网平差计算方法如下所示，具体包括以下几个方面：

[0039] 1)影像几何

[0040] 在经典摄影测量几何中，共线条件方程是区域网平差求解的基本方程，如式(1)-(4)所示，共线条件方程将像点坐标，相机的内部参数，相机的外部位置和姿态以及物方点坐标联系起来，区域网平差过程就是以共线条件方程为基础建立像点观测值方程，然后建立法方程，通过求解法方程得到未知数的估计值。

[0041]

[0042]

[0043]

[0044]

[0045] 其中，(X,Y,Z)是地面点P的物方坐标，(x,y)是对应的像点坐标，(Xs,Ys,Zs,phi,omega,kappa)为相机的外方位元素，f为相机的焦距，(x0,y0)为相机的主点偏移,(k1,k2)为相机镜头的畸变参数。

[0046] 2)构建误差方程和法方程

[0047] 对每个像点观测值根据上述共线方程列出误差方程：

[0048] v＝Ax-l  (5)

[0049] 其中v为观测值残差向量，A为法方程系数矩阵，由观测值方程对未知数求一阶偏导得到，x为未知数改正数向量,l为误差方程常数项向量，由像点坐标的计算值减去像点坐标观测值得到。

[0050] 根据式(5)可以列出法方程：

[0051] ATAx＝ATl  (6)

[0052] 为了加强法方程求解的稳定性，引入一个阻尼项(Damping term)λD，避免法方程的奇异性对求解造成不稳定影响，新的法方程如下所示：

[0053] (ATA+λD)x＝ATl  (7)

[0054] 其中λ为阻尼系数，其取值范围为(0,1)，矩阵D是一个对角矩阵，对角线上的元素T与矩阵AA的对角线上元素相等，根据每次迭代的结果可改变λ的取值，以增强法方程的稳定性。

[0055] 3)改化法方程

[0056] 法方程系数矩阵A可以被分成两个部分，相机参数(包括内参数和外参数)部分和地面点坐标部分，因而矩阵A可以写作A＝[AC AP]，其中AC代表相机参数部分，AP代表地面点坐标部分。同理可以得到，D＝[DC DP]，x＝[xc xp]，此时，方程(7)可以表示为以下形式：

[0057]

[0058] 令 则可得到下式：

[0059]

[0060] 其中，VC和VP都是块对角矩阵，一般情况下，地面点坐标未知数远远大于相机参数未知数个数，因此可以采用基于块状矩阵的高斯消元法，将地面点坐标未知数消元，得到改化后的法方程：

[0061]

[0062] VPxp＝bp-WTuc  (11)

[0063] 其中， 既是改化后的法方程系数矩阵，此时，相机参数未知数xc可以通过求解改化法方程(10)得到，地面点未知数xp则可以根据回代方程(11)计算得到。由于VP是块对角矩阵，因此其逆矩阵可以通过分块求逆的方法快速计算得到。假设相机个数为m，地面点个数为n，则改化之前的法方程系数矩阵大小为(6m+3n)*(6m+3n)，改化后的法方程系数矩阵的大小为6m*6m，由于地面点个数n通常是一个较大的值，因此改化后的法方程大幅减少了法方程的大小，增加了数据处理容量，提高了求解效率。

[0064] 4)共轭梯度法

[0065] 共轭梯度法(Conjugate Gradients)是求解大型线性方程组的有效方法，该方法是由Hestenes和Stiefel在1952提出，其主要优点就是无需存储大规模的法方程，而只需要多次计算矩阵与向量的乘积即可，通过迭代搜索的方法，获取线性方程组的最优解。

[0066] 本发明对改化的法方程采用共轭梯度法求解，将式(10)改写为下式：

[0067] Bu＝c  (12)

[0068] 其中

[0069]

[0070] u＝xc  (14)

[0071]

[0072] 共轭梯度法的基本思想是采用迭代搜索的办法求解线性方程组的最优解。对于改化法方程(12)的求解，首先给定初始的未知数向量u0，然后利用线性方程组系数矩阵、常数项矩阵以及上述向量u0，采用共轭梯度算法计算新的未知数向量u1，如此反复循环迭代，当满足迭代收敛条件之后退出，此时得到的未知数向量un即为线性方程组的最优解，n为迭代收敛次数。

[0073] 采用上述方法进行迭代求解时，其理论迭代收敛次数为法方程系数矩阵B的条件数，为了进一步提高迭代收敛次数，可以通过引入预条件矩阵M的方法，降低法方程系数矩阵B的条件数，从而减少收敛所需的迭代次数。

[0074] 5)预条件共轭梯度法

[0075] 预条件共轭梯度法就是在共轭梯度法的基础上，在线性方程组系数矩阵之前，左乘一个预条件矩阵的逆M-1，以达到降低线性方程组系数矩阵的条件数的目的，从而提高收敛速度，对改化的法方程(12)引入预条件矩阵后，法方程变为:

[0076] M-1Bu＝M-1c  (16)

[0077] 此时改化法方程系数矩阵的条件数变为矩阵M-1B的条件数，预条件矩阵M的选取原则是构造简单，易于求逆，且能够有效减少改化法方程系数矩阵的条件数。本发明选取的预条件矩阵为Jacobi预条件矩阵，该矩阵由改化法方程对角线上的块状矩阵构成，出对角线上的块状矩阵外，其余元素均为0，该矩阵构造较为简单，对其进行求逆运算时，采用分块求逆的算法，可大幅提高求逆效率，且该矩阵能够有效减少改化法方程系数矩阵的条件数，因此是一个较为理想的预条件矩阵。

[0078] 6)不精确牛顿解

[0079] 在采用预条件共轭梯度法求解法方程时，收敛迭代次数取决于法方程系数矩阵的条件数，当测区较大的，未知数数量增大，法方程系数矩阵的条件数也随之增大，虽然采用了预条矩阵来减少其条件数，但其收敛速度仍然不够理想。同时，考虑到区域网平差求解未知数本身也是一个迭代逼近的过程，因此可以提前结束预条件共轭梯度法，得到法方程的一种不精确牛顿解，替代传统的法方程精确解，不影响区域网平差的迭代收敛，从而大幅减少了预条件共轭梯度法的迭代次数，提高了解算速度。

[0080] 为了得到不精确牛顿解，需要在预条件共轭梯度迭代中计算一组强制序列(forcing sequence)η，其计算公式如下：

[0081]

[0082] 其中，ηk是第k次迭代的强制序列系数，c为法方程常数项向量，sk+1为第k+1次迭代的残差向量。

[0083] 结合不精确牛顿解法，预条件共轭梯度法迭代求解法方程的具体流程和计算公式如下：

[0084] 预条件共轭梯度法的具体流程和计算公式如下：

[0085] 给定一般线性方程组:Bu＝c；

[0086] 给定预条件矩阵：M；

[0087] 设定初始值:u0；s0＝c-Bu0＝c；d0＝M-1s0＝M-1c；k＝0；

[0088] While|sk|＜Threshold

[0089] 1:

[0090] 2:xk+1＝xk+αkdk

[0091] 3:sk+1＝sk-αkBdk

[0092] 4:

[0093] 5:dk+1＝M-1sk+1+βkdk

[0094] 6:k＝k+1

[0095] 在本说明书的描述中，参考术语“实施例一”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体方法、装置或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不必须针对的是相同的实施例或示例。而且，描述的具体特征、方法、装置或者特点可以在任一个或多个实施例或示例中以合适的方式结合。此外，在不相互矛盾的情况下，本领域的技术人员可以将本说明书中描述的不同实施例或示例以及不同实施例或示例的特征进行结合和组合。

[0096] 以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。