本发明涉及测绘科学与技术领域，具体涉及一种基于稀疏块状矩阵压缩存储结构的预条件共轭梯度区域网平差方法，主要应用于中大规模测区(500到10000张影像)摄影测量4D产品生产以及三维建模等。
背景技术：
区域网平差是测绘科学与技术摄影测量分支领域的关键步骤，其主要目标是通过区域内影像之间的连接点以及物方控制点信息来恢复摄影时刻相机的位置和姿态，从而实现对地定位，为后端的4D产品生产以及三维建模等应用提供几何定位信息。区域网平差技术经过几十年的发展，其方法和流程已经相对成熟，并且在测绘领域得到了广泛的应用。然而随着科技的加速进步，新的传感器不断涌现，如航天领域的高分辨率卫星，立体测绘卫星，航空领域的倾斜航空摄影系统，无人机、飞艇摄影系统等。同时，全世界范围内三维建模应用需求不断增加，也使得大量地面车载摄影系统，近景摄影系统，普通数码相机，甚至是智能手机，网络图片库等采集的影像都被应用于三维建模。影像数据源越来越丰富的同时，其分辨率也不断提高，以前数十米的卫星影像如今可以达到最高0.35米(WorldView-3)，航空影像的分辨率更是进入了厘米级时代。分辨率的增加必然会带来数据量的增大，摄影时的航线设计也不再满足传统的条带式规则分布，给相应的数据处理方法带来了一定挑战，受法方程大小的限制，传统的区域网平差技术流程已经不能满足中大规模测区数据处理需求。
技术实现要素：
本发明所要解决的技术问题是提供一种基于稀疏块状矩阵压缩存储结构的预条件共轭梯度区域网平差方法，不足解决现有技术的不足。本发明解决上述技术问题的技术方案如下：本发明提供了一种基于稀疏块状矩阵压缩存储结构的预条件共轭梯度区域网平差方法，包括以下步骤：S1、导入区域网平差计算所需要的原始数据，其中，所述原始数据至少包括初始内外方位元素数据以及点位数据；S2、对导入的所述原始数据进行时空基准统一，并统计未知数分组情况以及每一组未知数的个数；S3、根据统计的未知数分组情况以及每一组未知数的个数，构建稀疏块状矩阵压缩存储结构；S4、逐点计算每个像点所生成的子法方程系数矩阵以及法方程常数项子向量，并计算每一个子法方程系数矩阵、法方程常数项子向量的大小以及它们在全局法方程系数矩阵和全局法方程常数项向量中的索引位置；S5、进入区域网平差迭代流程：初始化各类变量，包括全局法方程常数项向量c，未知数向量u，压缩后的全局法方程系数矩阵B以及预条件矩阵M；S6、根据步骤S4中记录的每一个子法方程系数矩阵在全局法方程系数矩阵B中的索引位置，将每一个子法方程系数矩阵累加至压缩后的全局法方程系数矩阵B中，且将位置位于对角线上的子法方程系数矩阵同时累加至预条件矩阵M中，将所有法方程常数项子向量累加至全局法方程常数项向量c中，得到完整的全局法方程常数项向量c，作为S8中残差向量s的初值；S7、对预条件矩阵M中的每一个子块分别求逆并乘以全局法方程常数项向量c中对应的子块，得到矩阵-向量积M-1c的一个分量，并将这个分量累加至矩阵-向量积M-1c中对应的位置，所有分量处理完毕后，得到完整的矩阵-向量积M-1c，作为步骤S8共轭梯度法迭代中方向向量d的初始值；S8、进入预条件共轭梯度法迭代流程：对未知数改正数向量u、残差向量s以及方向向量d进行初始化；S9、将压缩后的全局法方程系数矩阵B中的每一块子法方程系数矩阵与方向向量d中对应子块相乘，得到矩阵-向量积Bd的一个分量，并将这些分量累加至矩阵-向量积Bd中对应的位置，所有分量处理完毕后，得到完整的矩阵-向量积Bd；S10、根据预条件共轭梯度算法，以及本次迭代中的未知数改正数向量u，残差向量s，方向向量d，预条件矩阵M以及全局法方程系数矩阵B，计算新的未知数改正数向量u、新的残差向量s以及新的方向向量d；S11、判断此次预条件共轭梯度迭代是否符合预定收敛条件，若符合，则结束预条件共轭梯度迭代，输出新的未知数改正数向量u，执行步骤S12；否则，返回步骤S8；S12、判断此次区域网平差迭代是否符合预定收敛条件，若符合，则结束区域网平差迭代，执行步骤S13，否则，返回步骤S5；S13、根据步骤S11中输出的新的未知数改正数向量u，更新所有的未知数数值，并输出所有的未知数数值。本发明提供的一种基于稀疏块状矩阵压缩存储结构的预条件共轭梯度区域网平差方法，采用一种大规模稀疏矩阵压缩存储算法，对大数据区域网平差中产生的大规模法方程进行压缩存储和运算，有效减少对内存的需求，同时针对压缩矩阵的求逆问题，引入预条件共轭梯度法迭代求解法方程，避免了对压缩矩阵的直接求逆，克服了传统区域网平差中对于大规模法方程的存储和运算难题，提高了数据处理容量和效率。本发明特别适用于中大型测区(500到10000张影像)数据的区域网平差。附图说明图1为本发明实施例的一种基于稀疏块状矩阵压缩存储结构的预条件共轭梯度区域网平差方法流程图。具体实施方式以下结合附图对本发明的原理和特征进行描述，所举实例只用于解释本发明，并非用于限定本发明的范围。实施例、一种基于稀疏块状矩阵压缩存储结构的预条件共轭梯度区域网平差方法。下面结合图1对本实施例提供的方法进行说明。参见图1，本实施例提供的方法包括：S1、导入区域网平差计算所需要的原始数据，其中，所述原始数据至少包括初始内外方位元素数据以及点位数据；S2、对导入的所述原始数据进行时空基准统一，并统计未知数分组情况以及每一组未知数的个数。具体的，首先导入区域网平差计算需要的原始数据，原始数据主要包括初始内外方位元素数据、连接点数据、控制点数据以及检查点数据(有些时候不存在控制点数据和检查点数据)，其中，连接点数据、控制点数据以及检查点数据统称为点位数据。对原始数据进行预处理，具体包括对导入的原始数据进行时空基准统一，包括时间基准统一和空间基准统一，若原始数据中存在控制点数据，还应该将连接点物方坐标以及初始外方位元素数据的坐标转换至控制点数据所在的坐标系内，以及计算并统计区域网平差计算过程中未知数分组的情况，每一组中未知数的个数以及需要用到的预条件矩阵M的类型和大小等信息。本实施例中的预条件矩阵与法方程系数矩阵大小一致，但构造简单，更易于求逆，使用其逆矩阵左乘法方程系数矩阵后，使得法方程系数矩阵的条件数减少，从而加快共轭梯度法求解法方程的迭代收敛速度，并保持收敛的稳健性，常用的预条件矩阵有Jacobi预条件矩阵。S3、根据统计的未知数分组情况以及每一组未知数的个数，构建稀疏块状矩阵压缩存储结构。具体的，步骤S2中计算并统计了本次区域网平差计算中未知数分组情况，即有多少组未知数，每一组未知数的类型相同，以及每一组中未知数的个数，通常情况下每一组中未知数的个数是相同的，比如，本次区域网平差计算中有m组未知数，每一组中有n个未知数，则需要构建的稀疏块状矩阵为m*m个块组成的对称稀疏块状矩阵，其中，每个块中包括n个元素。S4、逐点计算每个像点所生成的子法方程系数矩阵以及法方程常数项子向量，并计算每一个子法方程系数矩阵、法方程常数项子向量的大小以及它们在全局法方程系数矩阵和全局法方程常数项向量中的索引位置。S5、进入区域网平差迭代流程：初始化各类变量，包括全局法方程常数项向量c，未知数向量u，压缩后的全局法方程系数矩阵B以及预条件矩阵M；S6、根据步骤S3中记录的每一个子法方程系数矩阵在全局法方程系数矩阵B中的索引位置，将每一个子法方程系数矩阵累加至压缩后的全局法方程系数矩阵B中，且将位置位于对角线上的子法方程系数矩阵同时累加至预条件矩阵M中，将法方程常数项子向量累加至全局法方程常数项向量c中，得到完整的全局法方程常数项向量c，作为S8中残差向量s的初值。首先，本实施例涉及的稀疏块状矩阵压缩存储技术如下：稀疏块状矩阵通常包含有很多零元素块，这些零元素块对矩阵的运算没有影响，因此可以放弃存储这些零元素块，仅存储非零元素块，若上述矩阵为对称矩阵时，则只需要存储上三角部分以及对角线部分的非零元素块即可，如表1所示。当稀疏块状矩阵中存在大量的零元素块时，此方法可以达到较高的压缩效率，在摄影测量区域网平差中，其法方程系数矩阵通常就是一个零元素块比例较高的稀疏块状对称矩阵，因此其法方程系数矩阵非常适合于采用本稀疏块状矩阵压缩存储结构进行压缩存储。为了在计算时准确恢复非零元素块在原始矩阵中的位置，还需要存储这些非零元素块的大小和位置，压缩存储结构如表2所示。表1表1中由10*10个块组成的对称稀疏块状矩阵，其中每个块大小一样，且都为6*6，空白区域为零元素块，黑色块和灰色块均为非零元素块，但仅黑色块需要进行存储。表2为了测试本实施例中压缩存储方法的压缩效率，分别统计表1中完整结构和表2中压缩结构在计算机中存储所需的内存空间，其结果如表3所示。压缩矩阵所需内存空间仅是完整矩阵的15.83％，而随着法方程的尺寸进一步增大，矩阵中零元素块占比增大，本方法的压缩效率将会更高。表3表1和表2中的矩阵结构在计算机中存储所需的内存对比本步骤S6具体的实施过程为：以物方点为单位，逐点构建误差方程、法方程以及改化法方程(本实施例中直接针对改化法方程进行预条件共轭梯度法求解，为描述方便，下文中的改化法方程统一简称为法方程)。对每个点计算得到的子法方程(包括子法方程系数矩阵部分和常数项向量部分)，对于其中的子法方程系数矩阵部分，根据步骤S3中记录的每一个子法方程系数矩阵在全局矩阵中的索引位置，将每一个子法方程系数矩阵(包括上三角的子法方程系数矩阵和对角线上的子法方程系数矩阵)累加至压缩后的全局法方程系数矩阵B中，对于位置处于对角线上的子法方程系数矩阵被同时累加至预条件矩阵M中，法方程常数项子向量部分被累加至全局法方程常数项向量c中。S7、对预条件矩阵M中的每一个子块分别求逆并乘以法方程常数项向量c中对应的子块，得到矩阵-向量积M-1c的一个分量，并将这个分量累加至矩阵-向量积M-1c中对应的位置，所有分量处理完毕后，得到完整的矩阵-向量积M-1c，作为步骤S8共轭梯度法迭代中方向向量d的初始值。S8、进入预条件共轭梯度法迭代流程：对未知数改正数向量u、残差向量s以及方向向量d进行初始化；具体的，将未知数改正数向量u的初始值设置为0，将全局法方程常数项向量c作为预条件共轭梯度法迭代中残差向量s的初始值，将矩阵-向量积M-1c作为预条件共轭梯度法迭代中方向向量d的初始值，进入预条件共轭梯度迭代流程。S9、将压缩后的法方程系数矩阵B中的每一块子法方程系数矩阵与方向向量d中对应的子块相乘，得到矩阵-向量积Bd的一个分量，并将这些分量累加至全局矩阵-向量积Bd中对应的位置，所有分量处理完毕后，得到完整的矩阵-向量积Bd。S10、根据预条件共轭梯度算法，以及本次迭代中的未知数改正数向量u，残差向量s，方向向量d，预条件矩阵M以及全局法方程系数矩阵B，计算新的未知数改正数向量u、新的残差向量s以及新的方向向量d。S11、判断此次预条件共轭梯度迭代是否符合预定收敛条件，若符合，则结束预条件共轭梯度迭代，输出新的未知数改正数向量u，执行步骤S12；否则，返回步骤S8；具体的，统计步骤S10中计算得到的新的方向向量d中所有元素的绝对值最大值，若该绝对值最大值小于第一给定阈值或者迭代次数大于第二给定阈值，则结束预条件共轭梯度迭代，输出新的未知数改正数向量u，执行步骤S12，否则，返回步骤S8。本实施例中设置预条件矩阵，减少法方程系数矩阵的条件数，加快迭代收敛速度，共轭梯度法的特征是无需对法方程系数矩阵进行求逆，通过迭代搜索线性方程组的最优解，且每次迭代搜索方向之间时互相共轭的，具有存储量少，计算方便，收敛快等特点。S12、判断此次区域网平差迭代是否符合预定收敛条件，若符合，则结束区域网平差迭代，执行步骤S13，否则，返回步骤S5；具体的，统计步骤S11中输出的新的未知数改正数向量u中所有元素的绝对值最大值，若该绝对值最大值小于第三给定阈值或者迭代次数大于第四给定阈值，则结束区域网平差迭代，执行步骤S13，否则，返回步骤S5。S13、根据步骤S11中输出的新的未知数改正数向量u，更新所有的未知数数值，并输出所有的未知数数值。本方法涉及的利用预条件共轭梯度法迭代的区域网平差计算方法如下所示，具体包括以下几个方面：1)影像几何在经典摄影测量几何中，共线条件方程是区域网平差求解的基本方程，如式(1)-(4)所示，共线条件方程将像点坐标，相机的内部参数，相机的外部位置和姿态以及物方点坐标联系起来，区域网平差过程就是以共线条件方程为基础建立像点观测值方程，然后建立法方程，通过求解法方程得到未知数的估计值。其中，(X,Y,Z)是地面点P的物方坐标，(x,y)是对应的像点坐标，(Xs,Ys,Zs,phi,omega,kappa)为相机的外方位元素，f为相机的焦距，(x0,y0)为相机的主点偏移,(k1,k2)为相机镜头的畸变参数。2)构建误差方程和法方程对每个像点观测值根据上述共线方程列出误差方程：v＝Ax-l(5)其中v为观测值残差向量，A为法方程系数矩阵，由观测值方程对未知数求一阶偏导得到，x为未知数改正数向量,l为误差方程常数项向量，由像点坐标的计算值减去像点坐标观测值得到。根据式(5)可以列出法方程：ATAx＝ATl(6)为了加强法方程求解的稳定性，引入一个阻尼项(Dampingterm)λD，避免法方程的奇异性对求解造成不稳定影响，新的法方程如下所示：(ATA+λD)x＝ATl(7)其中λ为阻尼系数，其取值范围为(0,1)，矩阵D是一个对角矩阵，对角线上的元素与矩阵ATA的对角线上元素相等，根据每次迭代的结果可改变λ的取值，以增强法方程的稳定性。3)改化法方程法方程系数矩阵A可以被分成两个部分，相机参数(包括内参数和外参数)部分和地面点坐标部分，因而矩阵A可以写作A＝[ACAP]，其中AC代表相机参数部分，AP代表地面点坐标部分。同理可以得到，D＝[DCDP]，x＝[xcxp]，此时，方程(7)可以表示为以下形式：令则可得到下式：其中，VC和VP都是块对角矩阵，一般情况下，地面点坐标未知数远远大于相机参数未知数个数，因此可以采用基于块状矩阵的高斯消元法，将地面点坐标未知数消元，得到改化后的法方程：VPxp＝bp-WTuc(11)其中，既是改化后的法方程系数矩阵，此时，相机参数未知数xc可以通过求解改化法方程(10)得到，地面点未知数xp则可以根据回代方程(11)计算得到。由于VP是块对角矩阵，因此其逆矩阵可以通过分块求逆的方法快速计算得到。假设相机个数为m，地面点个数为n，则改化之前的法方程系数矩阵大小为(6m+3n)*(6m+3n)，改化后的法方程系数矩阵的大小为6m*6m，由于地面点个数n通常是一个较大的值，因此改化后的法方程大幅减少了法方程的大小，增加了数据处理容量，提高了求解效率。4)共轭梯度法共轭梯度法(ConjugateGradients)是求解大型线性方程组的有效方法，该方法是由Hestenes和Stiefel在1952提出，其主要优点就是无需存储大规模的法方程，而只需要多次计算矩阵与向量的乘积即可，通过迭代搜索的方法，获取线性方程组的最优解。本发明对改化的法方程采用共轭梯度法求解，将式(10)改写为下式：Bu＝c(12)其中u＝xc(14)共轭梯度法的基本思想是采用迭代搜索的办法求解线性方程组的最优解。对于改化法方程(12)的求解，首先给定初始的未知数向量u0，然后利用线性方程组系数矩阵、常数项矩阵以及上述向量u0，采用共轭梯度算法计算新的未知数向量u1，如此反复循环迭代，当满足迭代收敛条件之后退出，此时得到的未知数向量un即为线性方程组的最优解，n为迭代收敛次数。采用上述方法进行迭代求解时，其理论迭代收敛次数为法方程系数矩阵B的条件数，为了进一步提高迭代收敛次数，可以通过引入预条件矩阵M的方法，降低法方程系数矩阵B的条件数，从而减少收敛所需的迭代次数。5)预条件共轭梯度法预条件共轭梯度法就是在共轭梯度法的基础上，在线性方程组系数矩阵之前，左乘一个预条件矩阵的逆M-1，以达到降低线性方程组系数矩阵的条件数的目的，从而提高收敛速度，对改化的法方程(12)引入预条件矩阵后，法方程变为:M-1Bu＝M-1c(16)此时改化法方程系数矩阵的条件数变为矩阵M-1B的条件数，预条件矩阵M的选取原则是构造简单，易于求逆，且能够有效减少改化法方程系数矩阵的条件数。本发明选取的预条件矩阵为Jacobi预条件矩阵，该矩阵由改化法方程对角线上的块状矩阵构成，出对角线上的块状矩阵外，其余元素均为0，该矩阵构造较为简单，对其进行求逆运算时，采用分块求逆的算法，可大幅提高求逆效率，且该矩阵能够有效减少改化法方程系数矩阵的条件数，因此是一个较为理想的预条件矩阵。预条件共轭梯度法的具体流程和计算公式如下：给定一般线性方程组:Bu＝c；给定预条件矩阵：M；设定初始值:u0；s0＝c-Bu0＝c；d0＝M-1s0＝M-1c；k＝0；While|sk|＜Threshold1:2:xk+1＝xk+αkdk3:sk+1＝sk-αkBdk4:5:dk+1＝M-1sk+1+βkdk6:k＝k+1本发明提供的一种基于稀疏块状矩阵压缩存储结构的预条件共轭梯度区域网平差方法，采用一种大规模稀疏矩阵压缩存储算法，对大数据区域网平差中产生的大规模法方程进行压缩存储和运算，有效减少对内存的需求，同时针对压缩矩阵的求逆问题，引入预条件共轭梯度法迭代求解法方程，避免了对压缩矩阵的直接求逆，克服了传统区域网平差中对于大规模法方程的存储和运算难题，提高了数据处理容量和效率。本发明特别适用于中大型测区(500到10000张影像)数据的区域网平差。在本说明书的描述中，参考术语“实施例一”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体方法、装置或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不必须针对的是相同的实施例或示例。而且，描述的具体特征、方法、装置或者特点可以在任一个或多个实施例或示例中以合适的方式结合。此外，在不相互矛盾的情况下，本领域的技术人员可以将本说明书中描述的不同实施例或示例以及不同实施例或示例的特征进行结合和组合。以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。