一种基于稀疏矩阵的大区域网平差快速并行解算方法与流程

本发明涉及一种快速并行解算方法，尤其涉及一种基于稀疏矩阵的大区域网平差快速并行解算方法。

背景技术

光束法区域网平差在摄影测量占重要的位置，目前随着卫星发射数量的增加，卫星影像也更容易获取，因此，数据量的急速增加和卫星影像的实时处理对大区域网平差的运算内存占用和运算效率都提出了较高的要求。光束法区域网平差需要对影像的定向参数和连接点物方坐标同时求解，实际作业时，构建的区域网中影像数量较多，待求的连接点数目可达到10⁵，构建的误差方程的系数矩阵规模将达到至少6*10⁵×3*10⁵阶，法方程的系数矩阵运算量为3*10⁵×6*10⁵×3*10⁵阶，通常使用lu分解等方法对法方程进行答解，其计算量大约为另一方面，申请6*10⁵×3*10⁵阶的矩阵，大约需要(6*10⁵×3*10⁵×8byte)1440g内存，这样大的数据存储量与计算量，对一般的计算机是难以实现的。目前区域网平差的算法多对法方程直接运算，多借助lu分解、qr分解等矩阵分解算法进行求解，导致构建法方程系数矩阵时内存占用过大，运算效率也较低。

技术实现要素：

为了解决上述技术所存在的不足之处，本发明提供了一种基于稀疏矩阵的大区域网平差快速并行解算方法。

为了解决以上技术问题，本发明采用的技术方案是：一种基于稀疏矩阵的大区域网平差快速并行解算方法，其整体步骤为：

一、借助稀疏矩阵的性质，对光束法区域网平差进行改化与并行运算；

二、采用矩阵运算后直接填充的方式，仅存储改化后的法方程的上三角矩阵；

三、针对改化后的法方程为实对称矩阵的性质，对法方程矩阵进行分块并行求逆；

四、对连接点物方坐标进行并行解算。

进一步地，步骤一的具体过程为：

ⅰ、摄影测量中，根据卫星影像的构像模型，将光束法区域网平差的误差方程写成公式1：

v＝at+bx-l公式1

其中：v为像点坐标改正数向量；t＝[a0a1a2b0b1b2]^t，为模型定向参数改正向量，t为向量转置；a为定向参数改正系数矩阵；x＝[δxδyδz]^t，为连接点物方改正向量，t为向量转置；b为连接点物方坐标改正数系数矩阵；l为像点坐标观测值残差向量；

对公式1法化得到相应的法方程，如公式2：

令n11＝a^tpa，n12＝a^tpb，n21＝b^tpa，n22＝b^tpb，n1＝a^tpl，n2＝b^tpl，则将公式2变形为公式3：

对公式3消去一类未知数，得公式4：

ⅱ、对法方程进行改化：

将法方程公式3展开，得公式5：

式中，i为连接点点号，j为像片序号，m为总像片数目，n为总连接点数目；

根据公式2到公式3的变形规则，将公式5变形为：

通过公式6得出以下规律：

a、法方程中n11和n22均为块对角矩阵，n12和n21互为转置矩阵；

b、n11中对角线上元素为每景影像上的连接点对应的系数矩阵a自法化累加而得；

c、n22中对角线上元素为单个连接点对应的所有同名像点对应的系数矩阵b自法化累加而得；

d、n12每个元素为当前像点在当前影像上的系数矩阵a和b互法化而得，若像点不在影像上则n12对应的元素为0。

进一步地，步骤二的具体过程为：

根据公式6得到的规律，对定向元素解算的矩阵运算公式进行分析，得公式7：

t＝(n11-n12n22^-1n12^t)^-1(n1-n12n22^-1n2)公式7

根据公式7，进行以下处理：

a、n22为块对角矩阵，所以将n22^-1分解为每个块的矩阵求逆；

b、n11-n12n22^-1n12^t矩阵的存储方式仅与影像编号有关：对该矩阵按照矩阵的大小进行分块，各分块矩阵的元素按照影像编号由该影像上的所有连接点计算后累加存储到相应位置；

c、n11-n12n22^-1n12^t为实对称矩阵，该矩阵的存储采用上三角矩阵存储；

d、同理，向量也根据相应规律进行填充；

e、矩阵填充过程中，由于各连接点的矩阵运算是各自独立的，利用多个线程计算矩阵对应的相乘项，然后进行求和运算。

进一步地，步骤三的具体过程为：

法方程运算过程中，由于n11-n12n22^-1n12^t为实对称矩阵，用分块迭代法实现其并行运算：

令实对称矩阵wt＝n11-n12n22^-1n12^t，则wt的分块形式为：

式中，t为迭代次数，wt-1为wt的前t-1阶方阵，rt^t为矩阵wt的第t行的前t-1个元素，rt为rt^t的转置矩阵，pt为矩阵的第t行第t列元素；

设为wt的逆矩阵，则有：

其中，it-1为t-1阶单位矩阵，根据公式9可推导出：

令：则wt的求逆公式为：

利用公式11，wt的逆矩阵由递推进行并行求解。

进一步地，步骤四的具体过程为：

求得各张像片的外方位元素改正数之后，由改化法方程式求解结果反代回法方程式，求解连接点物方坐标的改正：

式中，tj为第j张影像的定向参数改正数；i为连接点点号，j为像片序号，m为总像片数目，n为总连接点数目。

本发明通过对区域网平差中大型稀疏矩阵运算的分析，针对区域网平差法方程是实对称矩阵的性质，对改化后的法方程进行分块解算，同时对区域网平差算法进行并行化改进，进一步提高运算效率。

附图说明

图1为常规的光束法区域网平差解算流程图。

图2为实对称矩阵并行分块求逆解算的流程图。

图3为本发明的整体流程示意图。

具体实施方式

下面结合附图和具体实施方式对本发明作进一步详细的说明。

实施例一、

一种基于稀疏矩阵的大区域网平差快速并行解算方法，其整体步骤为：

一、借助稀疏矩阵的性质，对光束法区域网平差进行改化与并行运算,具体过程为：

ⅰ、摄影测量中，根据卫星影像的构像模型，将光束法区域网平差的误差方程写成公式1：

v＝at+bx-l公式1

其中：v为像点坐标改正数向量；t＝[a0a1a2b0b1b2]^t，为模型定向参数改正向量，t为向量转置；a为定向参数改正系数矩阵；x＝[δxδyδz]^t，为连接点物方改正向量，t为向量转置；b为连接点物方坐标改正数系数矩阵；l为像点坐标观测值残差向量；

对公式1法化得到相应的法方程，如公式2：

令n11＝a^tpa，n12＝a^tpb，n21＝b^tpa，n22＝b^tpb，n1＝a^tpl，n2＝b^tpl，则将公式2变形为公式3：

对公式3消去一类未知数，得公式4：

依据公式4进行传统光束法区域网平差的整体流程如图1所示。

ⅱ、对法方程进行改化：

光束法区域网平差是一个迭代计算的算法，每次迭代之间存在着强烈的依赖关系，因此，整个算法并不能使用并行运算，但在每个运算步骤中间都存在着大量的矩阵和向量间的相乘与相加运算，且矩阵之间的运算相互独立，因此本方法进一步通过对每个运算步骤中的矩阵运算进行并行优化设计，可有效提升运算速度。

为了方便说明上述算法的填充规律及各部分的意义，将法方程公式3展开，得公式5：

式中，i为连接点点号，j为像片序号，m为总像片数目，n为总连接点数目；

根据公式2到公式3的变形规则，将公式5变形为：

通过公式6得出以下规律：

a、法方程中n11和n22均为块对角矩阵，n12和n21互为转置矩阵；

b、n11中对角线上元素为每景影像上的连接点对应的系数矩阵a自法化累加而得；

c、n22中对角线上元素为单个连接点对应的所有同名像点对应的系数矩阵b自法化累加而得；

d、n12每个元素为当前像点在当前影像上的系数矩阵a和b互法化而得，若像点不在影像上则n12对应的元素为0。

二、采用矩阵运算后直接填充的方式，仅存储改化后的法方程的上三角矩阵,具体过程为：

根据公式6得到的规律，对定向元素解算的矩阵运算公式进行分析，得公式7：

t＝(n11-n12n22^-1n12^t)^-1(n1-n12n22^-1n2)公式7

根据公式7，进行以下处理：

a、n22为块对角矩阵，所以将n22^-1分解为每个块的矩阵求逆；

b、n11-n12n22^-1n12^t矩阵的存储方式仅与影像编号有关：对该矩阵按照矩阵的大小进行分块，各分块矩阵的元素按照影像编号由该影像上的所有连接点计算后累加存储到相应位置；

c、n11-n12n22^-1n12^t为实对称矩阵，该矩阵的存储采用上三角矩阵存储；

d、同理，向量也根据相应规律进行填充；

e、矩阵填充过程中，由于各连接点的矩阵运算是各自独立的，利用多个线程计算矩阵对应的相乘项，然后进行求和运算。

三、针对改化后的法方程为实对称矩阵的性质，对法方程矩阵进行分块并行求逆,具体过程为：

法方程运算过程中，由于n11-n12n22^-1n12^t为实对称矩阵，因此可用分块迭代法实现其并行运算：

令实对称矩阵wt＝n11-n12n22^-1n12^t，则wt的分块形式为：

式中，t为迭代次数，wt-1为wt的前t-1阶方阵，rt^t为矩阵wt的第t行的前t-1个元素，rt为rt^t的转置矩阵，pt为矩阵的第t行第t列元素；

设为wt的逆矩阵，则有：

其中，it-1为t-1阶单位矩阵，qt^t、qt、αt只是一种表现形式，方便后面公式推导，没有具体含义；

根据公式9可推导出：

令：则wt的求逆公式为：

bt、βt同样只是一种表现形式，方便后面公式推导，没有具体含义；bt^t仅为了简化矩阵书写大小，没有具体含义；

利用公式11，wt的逆矩阵由递推进行并行求解。法方程系数的分块迭代求逆属于循环迭代的计算方法，但每个循环都包含大量的矩阵相乘、相加运算，且相互独立，因此对其进行并行化设计可行，该步骤的具体流程如图2所示。

四、对连接点物方坐标进行并行解算,具体过程为：

求得各张像片的外方位元素改正数之后，由改化法方程式求解结果反代回法方程式，求解连接点物方坐标的改正：

式中，tj为第j张影像的定向参数改正数；i为连接点点号，j为像片序号，m为总像片数目，n为总连接点数目。

由于空间点之间相互独立，因此可对连接点物方坐标的解算进行并行优化。综上所述，本发明的整体算法流程如图3所示。

选用四组试验数据对实施例一中的方法进行试验：

试验使用的电脑配置为：dell7040，内存8g。

试验利用传统方法直接申请系数矩阵进行解算内存占用过大，计算失败。利用改化后和并行优化后的算法对4组数据进行区域网平差，并行优化后对内存的占用很小，单次迭代的运行时间统计结果如表1所示。

表1测试数据统计结果

本发明的关键创新点及优势为：

1.借助稀疏矩阵的性质，对光束法区域网平差进行改化与并行运算，提高了运算效率；

2.采用矩阵运算后直接填充的方式，运算过程中只需存储改化后的法方程的上三角矩阵，可进一步降低对内存的要求；

3.针对改化后的法方程为实对称矩阵的性质，对法方程矩阵进行分块并行求逆，进一步提高了运算速度。

上述实施方式并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的技术方案范围内所做出的变化、改型、添加或替换，也均属于本发明的保护范围。