CN103175516B - 规模化大地控制网平差的分布式计算方法 - Google Patents

Abstract

本发明涉及规模化大地控制网平差的分布式计算方法，可实现规模化大地测量数据的高效处理，使用方便，其技术方案为，控制网平差计算任务分解，利用现代观测技术手段获取的与大地控制网中空间分布的控制点位置信息相关的大量的距离、角度、高差、基线观测数据，每分区单独平差计算，然后把各分区联接成一个整体，各个分区之间的内部参数项彼此相互独立；建立观测方程，将大地控制网分区分解，建立非线性误差方程，对法方程进行解算，建立平差计算服务，平差计算服务分布式的实现，对解算结果进行迭代优化，直到所有未知参数的改正项达到求解精度要求，本发明具有高效的平差计算效率，提高了现有计算资源的利用效率，其经济和社会效益巨大。

Description

规模化大地控制网平差的分布式计算方法

技术领域

本发明涉及“测绘科学与技术”学科中的“大地测量学与测量工程”技术领域，特别是一种规模化大地控制网平差的分布式计算方法。

背景技术

在局域或全球范围内建立由一系列地面控制点构成的大地控制网，通过精密测定控制点的三维位置及其时空变率，研究和测定地球形状、大小及其变化，能够为测制地图、工程设计和地球动力学等科学领域提供控制基础，为人造卫星、远程武器和各种航天器的发射和运行提供精确控制信息，为国家经济与国防建设提供控制保障。

当前全球导航卫星系统(GNSS,Global Navigation Satellite System)逐步成为测定地面点位置的主要技术手段，随着越来越多的连续运行参考站(CORS,Continuously OperatingReference Station)的建立和多个卫星系统的投入使用，大地控制网的数据处理面临日益突出的有限计算资源和不断增长的数据规模的矛盾。大地控制网的数据处理中涉及的平差计算不仅面临海量观测数据的事后精确计算，还要应对大型网络实时数据的实时快速解算。

目前，传统的大地控制网平差的计算方法是将时空分布的测站所采集的观测数据传输到数据处理中心，数据处理中心按照一定的平差模型对观测数据进行平差处理与分析，再将成果发布给用户。然而，集中式的数据处理中心现有计算能力逐渐无法满足急剧增长的规模化大地控制网平差计算的处理需求。如当前国际GNSS服务组织(IGS,International GNSSService)的核心CORS站已有200多个，再加上各种区域CORS网的测站，可用CORS数目已多达2000多个，而IGS数据处理中心目前限于计算能力只能同时处理100个以内的CORS站观测数据。大量的原始观测数据没有得到有效利用，造成相应的平差计算成果精度不稳定或不可靠。

为了实现大地测量控制网平差的高效处理，国内外研究机构和研究人员进行了广泛的关注和研究，发表的文献主要包括：《Journalof Geodesy》的《A new data processing strategyfor huge GNSS global network》，《Journal of Geophysical Research》的《Fixed pointtheorems of GPS carrier phase ambiguity resolution and their application to massivenetwork processing:Ambizap》，《测绘学报》的《超大规模大地网分区平差快速解算方法》。此类研究均是从降低传统集中式平差计算的算法复杂度等角度出发，对大规模大地控制网平差的快速解算带来了一定的效果。但相关的平差模型与算法多是针对单处理器体系架构编写的，没有顾及到当前多机多核硬件平台所带来的高效分布并行计算性能。

近年来，随着计算机硬件平台的高速发展，以分布式计算技术为主的高性能计算成为大规模数据处理的最有效方法。针对规模化大地控制网平差的集中式计算方法的处理瓶颈，本发明利用了分布式计算技术实现了规模化大地控制网平差的分布式计算。本发明的实质是对规模化大地控制网平差计算任务的并行分解与并行计算，易于实现，便于掌握。本发明方法顾及规模化大地控制网观测数据的获取和相关计算资源（CORS站网、联网数据中心等）本身的分布特性，实现规模化大地控制网平差的分布式计算方法，不仅提高对现有计算资源的利用率，而且可以使人们在高性能的计算能力下解决更复杂的大地测量数据处理问题。

发明内容

本发明的目的是克服现有大地控制网平差的集中式计算技术的不足，提供一种规模化大地控制网平差的分布式计算方法，可实现规模化大地测量数据的高效处理，使用方便，其技术方案为：一种规模化大地控制网平差的分布式计算方法，该方法包括以下步骤：

步骤1：控制网平差计算任务分解，利用现代观测技术手段获取的与大地控制网中空间分布的控制点位置信息相关的大量的距离、角度、高差、基线观测数据，具有空间分布和连续采集的特性，根据观测数据的观测条件、观测仪器型号、观测时期、控制点地理位置、控制点点数的原则综合将规模化大地控制网分成若干个区，每分区单独平差计算，然后把各分区联接成一个整体，将联结各区的控制点称为公共点，各分区内的控制点称为内部点，则规模化大地控制网平差计算任务划分为多个子计算任务，采用分区平差方法后，控制网平差计算的未知参数分为公共参数项和内部参数项两类，各个分区之间的内部参数项彼此相互独立；

所述的现代观测技术手段，是指利用电子经纬仪、电磁波测距仪、全站仪、精密水准仪、全球导航卫星系统、激光测卫、甚长基线干涉测量的一种或两种以上的组合的技术手段；

步骤2：建立观测方程，将大地控制网分区分解后共有S个分区，第i个分区的观测数据用L_i表示，由于观测存在误差，其对应的改正项用V_i表示，该分区的内部参数项用表示，i=1,2,…,S，大地控制网的公共参数项用表示，根据大地控制网平差中观测数据L_i、观测数据的改正项V_i和内部参数项公共参数项存在数学函数f，有如下非线性误差方程：

$V_i=f({\hat{X}}_i,\hat{Y})-L_i$ 式(1)

设内部参数项的近似值为X_i0，内部参数项的改正项为公共参数项的近似值为Y₀，公共参数项的改正项为则满足

${\hat{X}}_i=X_{i0}+\mathrm{\δ}{\hat{X}}_i$ $\hat{Y}=Y_0+\mathrm{\δ}\hat{Y}$

则式(1)可写做：

$V_i=f(X_{i0}+\mathrm{\δ}{\hat{X}}_i,Y_0+\mathrm{\δ}\hat{Y}){-L}_i$ 式(2)

按泰勒级数展开式(2)，由于和均很小，舍去二次幂以上项，得：

$V_i=f(X_{i0},Y_0)+{\left(\frac{\∂f}{\∂{\hat{X}}_i}\right)}_0\∂{\hat{X}}_i+{\left(\frac{\∂f}{\∂\hat{Y}}\right)}_0\∂{\hat{Y}-L}_i$ 式(3)

式(3)中，f(X_i0,Y₀)表示取近似值X_i0和Y₀时函数f对应的函数值：

${\left(\frac{\∂f}{\∂{\hat{X}}_i}\right)}_0,$ ${\left(\frac{\∂f}{\∂\hat{Y}}\right)}_0$

分别是函数f在内部参数项和公共参数项取近似值时的偏导数，故皆为常数，令

${\left(\frac{\∂f}{\∂{\hat{X}}_i}\right)}_0=A_i,$ ${\left(\frac{\∂f}{\∂\hat{Y}}\right)}_0=B_i,$ L_i-f(X_i0,Y₀)=l_i

则线性化后的误差方程为：

$V_i=A_i\mathrm{\δ}{\hat{X}}_i+B_i\mathrm{\δ}\hat{Y}-l_i$ 式(4)

在误差理论与平差计算理论中，A_i和B_i被称为系数矩阵，l_i称为常数项，因i=1,2,…,S，则S个分区的误差方程组的矩阵形式如下：

式(5)

其中V_i为第i分区观测数据的改正项，A_i和B_i为第i分区系的系数矩阵，l_i为第i分区的常数项，为第i分区的内部参数项的改正项，为公共参数项的改正项，每个分区的相关常数项l_i的权矩阵为P_i，表示第i分区的观测数据彼此间可信赖程度的相对数值，矩阵中的空白处全为零元素，下同；

步骤3：法方程的解算，按照最小二乘平差原理，得到式(5)的法方程组的矩阵形式如下：

式(6)

为了简化上述法方程组的矩阵形式，设 $A_i^T{P}_i{A}_i=N_{\mathrm{Ai}},$ $A_i^T{P}_i{B}_i=N_i,$ $A_i^T{P}_i{l}_i=d_{\mathrm{Ai}},$ $B_i^T{P}_i{B}_i=N_{\mathrm{Bi}},$ $B_i^T{P}_i{l}_i=d_{\mathrm{Bi}},$ $\underset{i=1}{\overset{S}{\mathrm{\Σ}}}{B}_i^T{P}_i{B}_i=\underset{i=1}{\overset{S}{\mathrm{\Σ}}}{N}_{\mathrm{Bi}},$ $\underset{i=1}{\overset{S}{\mathrm{\Σ}}}{B}_i^T{P}_i{l}_i=\underset{i=1}{\overset{S}{\mathrm{\Σ}}}{d}_{\mathrm{Bi}},$ 式(6)可简写为：

式(7)

将上述矩阵形式的法方程组改写为如下法方程系统：

$\left\{\begin{array}{c}\∀i\∈[1,S]:N_{\mathrm{Ai}}\mathrm{\δ}{\hat{X}}_i+N_i^T\mathrm{\δ}\hat{Y}=d_{\mathrm{Ai}}\\ \underset{i=1}{\overset{S}{\mathrm{\Σ}}}(N_i\mathrm{\δ}{\hat{X}}_i+N_{\mathrm{Bi}}\mathrm{\δ}\hat{Y})=\underset{i=1}{\overset{S}{\mathrm{\Σ}}}\left(d_{\mathrm{Bi}}\right)\end{array}\right.$ 式(8)

由式(8)法方程系统中的第一个方程求解内部参数项的改正项

$\∀i\∈[1,S]:$ $\mathrm{\δ}{\hat{X}}_i=N_{\mathrm{Ai}}^{-1}(d_{\mathrm{Ai}}-N_i^T\mathrm{\δ}\hat{Y})$ 式(9)

代入式(8)法方程系统中的第二个方程得到仅与公共参数项的改正项相关的线性方程组：

$\left(\underset{i=1}{\overset{S}{\mathrm{\Σ}}}(N_{\mathrm{Bi}}-N_i{N}_{\mathrm{Ai}}^{-1}{N}_i^T)\right)\mathrm{\δ}\hat{Y}=\left(\underset{i=1}{\overset{S}{\mathrm{\Σ}}}(d_{\mathrm{Bi}}-N_i{N}_{\mathrm{Ai}}^{-1}{d}_{\mathrm{Ai}})\right)$ 式(10)

设 $(N_{\mathrm{Bi}}-N_i{N}_{\mathrm{Ai}}^{-1}{N}_i^T)=N_i^{\′},$ $(d_{\mathrm{Bi}}-N_i{N}_{\mathrm{Ai}}^{-1}{d}_{\mathrm{Ai}})=U_i^{\′},$ 则式(10)简写为：

$\underset{i=1}{\overset{S}{\mathrm{\Σ}}}\left(N_i^{\′}\right)\mathrm{\δ}\hat{Y}=\underset{i=1}{\overset{S}{\mathrm{\Σ}}}\left(U_i^{\′}\right)$ 式(11)

称N_i′和U_i′分别是第i分区对求解公共参数项的改正项的贡献矩阵和贡献向量，则各个分区的子计算任务可同时消去各个分区的内部参数项的改正项得到各个子计算任务对求解公共参数项的改正项的贡献矩阵N_i′和贡献向量N_i′，累加后求解线性方程组(11)得到改正项和评价其精度的协方差矩阵求解公式如下：

$\hat{Y}={\left(\underset{i=1}{\overset{S}{\mathrm{\Σ}}}\left(N_i^{\′}\right)\right)}^{-1}\left(\underset{i=1}{\overset{S}{\mathrm{\Σ}}}\left(U_i^{\′}\right)\right)$ 式(12)

$Q_{\mathrm{\δ}\hat{Y}}={\left(\underset{i=1}{\overset{S}{\mathrm{\Σ}}}\right)}^{-1}$

通过将求解的公共参数项的改正项带入式(9)中，即可得到各个内部参数项的改正项其对应的协方差矩阵为这里的协方差矩阵和反应出各自改正项的求解精度，矩阵对角线上的元素越小，表明该改正项的精度越高；

步骤4：建立平差计算服务,利用Web Service分布式计算技术编写步骤1到步骤3中涉及的大地控制网平差的子计算服务、求解公共参数项的改正项服务的平差计算服务，并发布在网络中的多台计算机上，其中子计算服务能够在本地计算机上根据观测数据建立线性化误差方程(4)，计算其对求解公共参数项的改正项的贡献矩阵N_i′和贡献向量N_i′，通过Web服务的形式对外提供贡献矩阵N_i′和贡献向量N_i′的结果服务，求解公共参数项的改正项服务提供对各子计算服务得到的对求解公共参数项的改正项的贡献矩阵N_i′和贡献向量N_i′进行累加求和、矩阵求逆、矩阵与向量相乘，对外提供最终的公共参数项的改正项和其协方差矩阵结果的服务；

步骤5：平差计算服务分布式的实现，大地控制网平差用户端或控制端计算机分布并行调用网络中多个子计算服务，得到各个子计算服务的计算结果，即对求解公共参数项的改正项的贡献矩阵N_i′和贡献向量N_i′，通过网络交互，调用求解公共参数项的改正项服务得到公共参数项的改正项和评价其精度的协方差矩阵通过广播公共参数项的改正项到网络中，每个本地计算机利用公共参数项的改正项计算得到内部参数项的改正项和评价其精度的协方差矩阵

步骤6：对解算结果进行迭代优化，将每次得到的大地控制网的公共参数项的改正项和所有内部参数项的改正项作为近似值带入到误差方程(2)中，重复步骤（2）到（5），经过迭代计算，直到所有未知参数的改正项达到求解精度要求。

与现有技术相比，本发明的有益效果：（1）计算高效，对各类规模化大地控制网平差计算任务通过分区分解后，将分布式计算技术带来的高性能计算能力结合分解后的分区平差模型实现其分布处理，极大地缩短了规模化大地控制网平差计算的时间，具有高效的平差计算效率；（2）实现简便，利用分布式计算技术建立平差分布服务体系具有操作简单、跨平台的特点，同时拥有松散耦合、粗粒度、可重用性等优点，可以大幅度降低系统耦合度，提供服务层上的集成应用，将各类平差计算服务发布到网络多台计算机实现分布并行调用，提高了现有计算资源的利用效率，其经济和社会效益巨大。

附图说明

图1为本发明大地控制网分区示意图。

图2为本发明大地控制网平差任务的分布并行计算流程图。

图3为本发明大地控制网平差分布式服务体系结构。

图4为本发明网络中的Web平差服务拓扑图。

具体实施方式

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

1）控制网平差计算任务分解。利用电子经纬仪、电磁波测距仪、全站仪、精密水准仪、全球导航卫星系统、激光测卫、甚长基线干涉测量观测技术手段获取的与大地控制网中空间分布的控制点位置信息相关的大量距离、角度、高差、基线等观测数据具有空间分布和连续采集的特性。大量的观测数据所建立的平差计算任务是一个规模比较大的平差计算任务，普通计算机内存容量有限，难以进行整体平差解算。采用分布式计算技术实现规模化大地控制网平差的解算，首先在于对平差任务所构建的平差模型进行分解，从而实现分布并行计算。

由于大量的观测数据本身是空间分布连续采集的，根据观测数据的观测条件、观测仪器型号、观测时期、控制点地理位置分布、控制点点数的原则综合将大地控制网分成若干个区，每分区单独平差，然后把各分区联接成一个整体。分区后，将联结各区的控制点称为公共点，各分区内的控制点称为内部点。则规模化大地控制网平差任务划分为多个子计算任务。

大地控制网平差计算任务主要是根据观测数据对控制网中的未知参数进行求解，而大地控制网中的未知参数除了包含各个控制点的三维坐标外，还有尺度参数、旋转因子参数等其他未知参数。对规模化大地控制网平差计算任务进行分区分解后，控制点被分为公共点和各区内的内部点两类。参照图1，为某大地控制网分区示意图，以虚线为界，该控制网综合被分为三个分区，编号为Ⅰ、Ⅱ和Ⅲ，其中控制点7、8、9、10、11、12、13点为公共点，控制点1属于分区Ⅰ的内部点，控制点2、3属于分区Ⅱ的内部点，控制点4、5、6属于分区Ⅲ的内部点，三角形符号表示的点A、B、C、D为网中三维坐标已知的控制点，用于提供控制网的基准信息。可见分区后，大地控制网平差计算未知参数分为公共参数项和内部参数项两类，各个分区内之间的内部参数是相互独立。其中公共参数项是由所有公共点的三维坐标参数构成，内部参数项是由各分区内部点的三维坐标参数、尺度参数、旋转因子参数构成。

2）误差方程建立。设大地控制网分区分解后共有S个分区，第i个分区的观测数据用L_i表示，由于观测存在误差，其对应的改正项用V_i表示，该分区的内部参数项用表示，i=1,2,…,S，大地控制网的公共参数项用表示，根据大地控制网平差中观测数据L_i、观测数据的改正项V_i和内部参数项公共参数项存在数学函数f，有如下非线性误差方程组：

$V_i=f({\hat{X}}_i,\hat{Y})-L_i$ 式(1)

设内部参数项的近似值为X_i0，内部参数项的改正项为公共参数项的近似值为Y₀，公共参数项的改正项为则满足

${\hat{X}}_i=X_{i0}+{\mathrm{\δ}\hat{X}}_i$ $\hat{Y}=Y_0+\mathrm{\δ}\hat{Y}$

则式(1)可写做：

$V_i=f(X_{i0}+\mathrm{\δ}{\hat{X}}_i,Y_0+\mathrm{\δ}\hat{Y}){-L}_i$ 式(2)

按泰勒级数展开式(2)，由于和均很小，舍去二次幂以上项，得：

$V_i=f(X_{i0},Y_0)+{\left(\frac{\∂f}{\∂{\hat{X}}_i}\right)}_0\mathrm{\δ}{\hat{X}}_i+{\left(\frac{\∂f}{\∂\hat{Y}}\right)}_0\mathrm{\δ}\hat{Y}-L_i$ 式(3)

式(3)中，f(X_i0,Y₀)表示取近似值X_i0和Y₀时函数f对应的函数值，

${\left(\frac{\∂f}{\∂{\hat{X}}_i}\right)}_0,{\left(\frac{\∂f}{\∂\hat{Y}}\right)}_0$

分别是函数f在内部参数项和公共参数项取近似值时的偏导数，故皆为常数，令

${\left(\frac{\∂f}{\∂{\hat{X}}_i}\right)}_0=A_i,$ ${\left(\frac{\∂f}{\∂\hat{Y}}\right)}_0=B_i,$ L_i-f(X_i0,Y₀)=l_i

则线性化后的误差方程组为：

$V_i=A_i\mathrm{\δ}{\hat{X}}_i+B_i\mathrm{\δ}\hat{Y}-l_i$ 式(4)

在误差理论与平差计算理论中，A_i和B_i被称为系数矩阵，l_i称为常数项，每个分区的相关常数项l_i的权矩阵为P_i，表示各观测数据彼此间可信赖程度的相对数值，因i=1,2,…,S，则S个分区的误差方程组的矩阵形式如下：

式(5)

其中V_i为第i分区观测数据的改正项，A_i和B_i为第i分区系的系数矩阵，l_i为第i分区的常数项，为第i分区的内部参数项的改正项，为公共参数项的改正项。矩阵中的空白处全为零元素，下同。

3）法方程的解算。

由于误差方程个数远远大于未知参数个数，对于不定解方程组（5），为了得到惟一解，引入一限制条件，即最小二乘原理。按照最小二乘原理，得到式(5)的法方程组的矩阵形式如下：

式(6)

为了简化上述法方程组的矩阵形式，设 $A_i^T{P}_i{A}_i=N_{\mathrm{Ai}},$ $A_i^T{P}_i{B}_i=N_i,$ $A_i^T{P}_i{l}_i=d_{\mathrm{Ai}},$ $B_i^T{P}_i{B}_i=N_{\mathrm{Bi}},$ $B_i^T{P}_i{l}_i=d_{\mathrm{Bi}},$ $\underset{i=1}{\overset{S}{\mathrm{\Σ}}}{B}_i^T{P}_i{B}_i=\underset{i=1}{\overset{S}{\mathrm{\Σ}}}{N}_{\mathrm{Bi}},$ $\underset{i=1}{\overset{S}{\mathrm{\Σ}}}{B}_i^T{P}_i{l}_i=\underset{i=1}{\overset{S}{\mathrm{\Σ}}}{d}_{\mathrm{Bi}},$ 式(6)可简写为：

式(7)

将上述矩阵形式的法方程组改写为如下法方程系统：

$\left\{\begin{array}{c}\∀i\∈[1,S]:N_{\mathrm{Ai}}\mathrm{\δ}{\hat{X}}_i+N_i^T\mathrm{\δ}\hat{Y}=d_{\mathrm{Ai}}\\ \underset{i=1}{\overset{S}{\mathrm{\Σ}}}(N_i\mathrm{\δ}{\hat{X}}_i+N_{\mathrm{Bi}}\mathrm{\δ}\hat{Y})=\underset{i=1}{\overset{S}{\mathrm{\Σ}}}\left(d_{\mathrm{Bi}}\right)\end{array}\right.$ 式(8)

由式(8)法方程系统中的第一个方程求解内部参数项的改正项

$\mathrm{\δ}{\hat{X}}_i=N_{\mathrm{Ai}}^{-1}\left(d_{\mathrm{Ai}}{N}_i^T\mathrm{\δ}\hat{Y}\right)$ 式(9)

代入式(8)法方程系统中的第二个方程得到仅与公共参数项的改正项相关的线性方程组：

$\left(\underset{i=1}{\overset{S}{\mathrm{\Σ}}}(N_{\mathrm{Bi}}-N_i{N}_{\mathrm{Ai}}^{-1}{N}_i^T)\right)\mathrm{\δ}\hat{Y}=\left(\underset{i=1}{\overset{S}{\mathrm{\Σ}}}(d_{\mathrm{Bi}}-N_i{N}_{\mathrm{Ai}}^{-1}{d}_{\mathrm{Ai}})\right)$ 式(10)

设 $(N_{\mathrm{Bi}}-N_i{N}_{\mathrm{Ai}}^{-1}{N}_i^T)=N_i^{\′},$ $(d_{\mathrm{Bi}}-N_i{N}_{\mathrm{Ai}}^{-1}{d}_{\mathrm{Ai}})=U_i^{\′},$ 则式(10)简写为：

$\underset{i=1}{\overset{S}{\mathrm{\Σ}}}\left(N_i^{\′}\right)\mathrm{\δ}\hat{Y}=\underset{i=1}{\overset{S}{\mathrm{\Σ}}}\left(U_i^{\′}\right)$ 式(11)

称N_i′是第i分区对求解公共参数项的改正项的贡献矩阵，N_i′是第i分区对求解公共参数项的改正项的贡献向量，则各个分区的子计算任务可同时消去各个分区的内部参数项的改正项得到各个子计算任务对求解公共参数项的改正项的贡献矩阵N_i′和贡献向量U_i′，累加后对矩阵求逆得到改正项和评价其精度的协方差矩阵累加求解公式如下：

$\mathrm{\δ}\hat{Y}={\left(\underset{i=1}{\overset{S}{\mathrm{\Σ}}}\left(N_i^{\′}\right)\right)}^{-1}\left(\underset{i=1}{\overset{S}{\mathrm{\Σ}}}\left(U_i^{\′}\right)\right)$ 式(12)

$Q_{\mathrm{\δ}\hat{Y}}{\left(\underset{i=1}{\overset{S}{\mathrm{\Σ}}}\left(N_i^{\′}\right)\right)}^{-1}$

通过将求解的公共参数项的改正项带入式(9)中，即可得到各个内部参数项项的改正项其对应的协方差矩阵为这里的协方差矩阵和反应出各自改正项的求解精度，矩阵对角线上的元素越小，表明该改正项的精度越高。

4）平差计算服务建立。利用Web Service分布式计算技术的标准化协议开发的应用程序能够在网络环境中自由地进行交互，建立可互操作的分布式应用新平台。将规模化大地控制网平差计算任务进行分解后，利用Web Service分布式计算技术编写步骤1到步骤3中涉及的大地控制网平差子计算服务、求解公共参数项项的改正项服务的平差计算服务，并发布在网络中的多台计算机上，其中子计算服务能够在本地计算机上根据观测数据建立线性化误差方程(4)，计算其对求解公共参数项的改正项的贡献矩阵Ni′和贡献向量N_i′，通过Web服务的形式对外提供贡献矩阵N_i′和贡献向量N_i′的结果服务，该服务可供网络中的平差用户端或控制端调用。求解公共参数项的改正项服务提供对各子计算服务得到的对求解公共参数项的改正项的贡献矩阵N_i′和贡献向量N_i′进行累加求和、矩阵求逆、矩阵与向量相乘，对外提供最终的公共参数项的改正项和其协方差矩阵结果的服务；

按照Web服务规范建立的各类大地控制网平差计算服务体系结构参照图3，在分布式网络环境中，分布的平差计算服务提供者利用Web服务描述语言（WSDL,Web ServicesDescription Language）描述大地控制网平差计算服务的功能并利用通用描述、发现和集成技术（UDDI,Universal Description,Discovery and Integration）发布该服务并注册到平差服务注册中心，用户通过平差注册中心查找相应的平差计算服务，然后进行绑定并利用简单对象访问协议技术（SOAP，Simple Object Access Protocol）调用服务或实现交互；

5）平差计算服务分布式实现，当前分布在各处的大地测量的计算资源通过网络已实现互联互通，而集中式处理使得网络中这些计算资源没有得到有效利用。基于Web Service分布式计算技术的规模化大地控制网平差分布式计算的实现拓扑图参照图4，大地测量观测数据统一放在公共Ftp服务器上或者分散在每个服务节点上，地控制网平差用户端或控制端通过网络绑定计算节点。每台计算机将建立的平差计算服务发布到网络中，即作为服务节点。实际中网络中每台计算机即可作为服务节点也可是客户端。

大地控制网平差用户端或控制端计算机分布并行调用网络中多个子计算服务，得到各个子计算的计算结果，即对求解公共参数项的改正项的贡献矩阵N_i′和贡献向量U_i′，通过网络交互，调用求解公共参数项的改正项服务后得到公共参数项的改正项和评价其精度的协方差矩阵进而广播公共参数项的改正项到网络中，每个本地计算机利用公共参数项的改正项计算得到内部参数项的改正项和评价其精度的协方差矩阵即实现网络中多台计算机分布协同完成规模化大地控制网平差的分布计算。

6）对解算结果进行迭代优化。

由于大地测量数据处理理论缺乏实际经验，需要将每次得到的大地控制网的公共参数项的改正项和所有内部参数项的改正项作为近似值带入到误差方程(2)中，重复步骤2）到5），得到新的一组平差参数改正项，经过迭代计算，直到所有的未知参数的改正项达到求解精度要求。

由于当控制网平差作业计算开销与通信等系统开销的比值越小，即并行粒度越细，通信量就相对增加，系统额外开销就增加，作为网络环境下的计算服务，分布式并行平差服务系统在分配平差作业任务时坚持“粒度粗、通信少”的总原则。用户端的程序根据服务节点个数以及服务能力，以指令的形式将不同运算规模的子计算作业分配到合适的计算节点上，通过多线程技术实现多个计算节点的并行处理。

大地控制网平差计算过程最终可归结为方程组的解算过程，则相应的数值计算均可实现并行计算，如矩阵的并行计算、线性方程组的并行求解等。通过在服务节点发布相应的数值分解计算服务，即可完成各类复杂的规模化控制网平差计算任务的分布式并行处理，从而提高计算效率。

所述的获取的与大地控制网中空间分布的控制点位置信息相关的大量的距离、角度、高差、基线观测数据，对大量的距离、角度、高差和基线观测数据采用分区平差方法进行分解后，控制网平差计算的未知参数项分为公共参数项和内部参数项两类，所述的公共参数项由所有的公共点的三维坐标参数构成，内部参数项由各分区内部点的三维坐标参数、尺度参数、旋转因子参数构成。

所述的同时消去各个分区的内部参数项的改正项即能够在分布式计算环境下并行消去通过并行计算各分区对公共参数项的改正项的贡献矩阵N_i′和贡献向量U_i′，得到仅与公共参数项的改正项相关的线性方程组(11)，然后通过求解线性方程组得到大地控制网的公共参数项的改正项和其协方差矩阵，进而各个分区利用公共参数项的改正项，并行得到各个分区的内部参数项的改正项和其协方差矩阵。

所述的大地控制网平差用户端或控制端计算机分布并行调用网络中多个子计算服务，是指大地控制网平差计算任务分解后，通过分布式计算技术建立控制网平差子计算服务、求解公共参数项的改正项服务的平差计算服务，并将平差计算服务发布到网络中多个计算节点，大地控制网平差用户端或控制端通过网络首先绑定计算节点，进而并行调用多个子计算服务，实现网络中多台计算机分布协同完成规模化大地控制网平差的分布计算。

本发明是一种新的大地控制网平差的分布计算方法，支持对规模化大地控制网平差计算的高效处理。因此本发明与现有技术相比，具有以下突出的有益技术效果：

（1）计算高效

本发明提出的方法具有较广泛的适用性，对各类规模化大地控制网平差计算任务通过分区分解后，将分布式计算技术带来的高性能计算能力结合分解后的分区平差模型实现其分布处理，极大地缩短了规模化大地控制网平差计算的时间，可用于完成规模化大地控制网数据处理涉及的平差计算任务，具有广泛适用性，解决了现有集中式处理面临的计算压力，带来高效的平差计算效率。

（2）实现简便

本发明提出的方法利用分布式计算技术建立平差分布服务体系具有操作简单、跨平台的特点，同时拥有松散耦合、粗粒度、可重用性等优点，可以大幅度降低系统耦合度，提供服务层上的集成应用。本发明方法只需要对现有的控制网平差处理软件进行粗粒度的Web服务开发，将各类控制网平差Web服务发布到网络多台计算机实现分布并行调用，极大地缩短了规模化大地控制网平差计算的时间，提高了现有计算资源的利用效率，其经济和社会效益巨大。

以上所述，仅为本发明较佳的具体实施方式，本发明的保护范围不限于此，任何熟悉本技术领域的技术人员在本发明披露的技术范围内，可显而易见地得到的技术方案的简单变化或等效替换均落入本发明的保护范围内。

Claims (4)

1.规模化大地控制网平差的分布式计算方法，其特征在于，该方法包括以下步骤：

步骤1：控制网平差计算任务分解，利用现代观测技术手段获取的与大地控制网中空间分布的控制点位置信息相关的大量的距离、角度、高差、基线观测数据，具有空间分布和连续采集的特性，根据观测数据的观测条件、观测仪器型号、观测时期、控制点地理位置、控制点点数的原则综合将规模化大地控制网分成若干个区，每分区单独平差计算，然后把各分区联接成一个整体，将联结各区的控制点称为公共点，各分区内的控制点称为内部点，则规模化大地控制网平差计算任务划分为多个子计算任务，采用分区平差方法后，控制网平差计算的未知参数分为公共参数项和内部参数项两类，各个分区之间的内部参数项彼此相互独立；

所述的现代观测技术手段，是指利用电子经纬仪、电磁波测距仪、全站仪、精密水准仪、全球导航卫星系统、激光测卫、甚长基线干涉测量的一种或两种以上的组合的技术手段；

步骤2：建立观测方程，将大地控制网分区分解后共有S个分区，第i个分区的观测数据用L_i表示，由于观测存在误差，其对应的改正项用V_i表示，该分区的内部参数项用表示，i＝1,2,…,S，大地控制网的公共参数项用表示，根据大地控制网平差中观测数据L_i、观测数据的改正项V_i和内部参数项公共参数项存在数学函数f，有如下非线性误差方程：

$V_i=f({\hat{X}}_i,\hat{Y})-L_i$    式(1)

设内部参数项的近似值为X_i0，内部参数项的改正项为公共参数项的近似值为Y₀，公共参数项的改正项为则满足

$\begin{array}{cc}{\hat{X}}_i=X_{i0}+\mathrm{\δ}{\hat{X}}_i& \hat{Y}=Y_0+\mathrm{\δ}\hat{Y}\end{array}$

则式(1)可写做：

$V_i=f(x_{i0}+\mathrm{\δ}{\hat{X}}_i,Y_0+\mathrm{\δ}\hat{Y})-L_i$    式(2)

按泰勒级数展开式(2)，由于和均很小，舍去二次幂以上项，得：

$V_i=f(X_{i0},Y_0)+{\left(\frac{\∂f}{\∂{\hat{X}}_i}\right)}_0\mathrm{\δ}{\hat{X}}_i+{\left(\frac{\∂f}{\∂\hat{Y}}\right)}_0\mathrm{\δ}\hat{Y}-L_i$    式(3)

式(3)中，f(X_i0,Y₀)表示取近似值X_i0和Y₀时函数f对应的函数值：

${\left(\frac{\∂f}{\∂{\hat{X}}_i}\right)}_0{,\left(\frac{\∂f}{\∂\hat{Y}}\right)}_0$

分别是函数f在内部参数项和公共参数项取近似值时的偏导数，故皆为常数，令

${\left(\frac{\∂f}{\∂{\hat{X}}_i}\right)}_0=A_i,{\left(\frac{\∂f}{\∂\hat{Y}}\right)}_0=B_i,L_i-f(X_{i0},Y_0)=l_i$

则线性化后的误差方程为：

$V_i=A_i\mathrm{\δ}{\hat{X}}_i+B_i\mathrm{\δ}\hat{Y}-l_i$    式(4)

在误差理论与平差计算理论中，A_i和B_i被称为系数矩阵，l_i称为常数项，因i＝1,2,…,S，则S个分区的误差方程组的矩阵形式如下：

   式(5)

其中V_i为第i分区观测数据的改正项，A_i和B_i为第i分区系的系数矩阵，l_i为第i分区的常数项，为第i分区的内部参数项的改正项，为公共参数项的改正项，每个分区的相关常数项l_i的权矩阵为P_i，表示第i分区的观测数据彼此间可信赖程度的相对数值，矩阵中的空白处全为零元素，下同；

步骤3：法方程的解算，按照最小二乘平差原理，得到式(5)的法方程组的矩阵形式如下：

   式(6)

为了简化上述法方程组的矩阵形式，设 $A_i^T{P}_i{A}_i=N_{\mathrm{Ai}},$ $A_i^T{P}_i{B}_i=N_i,$ $A_i^T{P}_i{l}_i=d_{\mathrm{Ai}},$

将上述矩阵形式的法方程组改写为如下法方程系统：

$\left\{\begin{array}{c}\∀i\∈[1,S]:N_{\mathrm{Ai}}\mathrm{\δ}{\hat{X}}_{\mathrm{Ai}}\mathrm{\δ}{\hat{X}}_i+N_i^T\mathrm{\δ}\hat{Y}=d_{\mathrm{Ai}}\\ \underset{i=1}{\overset{S}{\mathrm{\Σ}}}(N_i\mathrm{\δ}{\hat{X}}_i+N_{\mathrm{Bi}}\mathrm{\δ}\hat{Y})=\underset{i=1}{\overset{S}{\mathrm{\Σ}}}\left(d_{\mathrm{Bi}}\right)\end{array}\right.$    式(8)

由式(8)法方程系统中的第一个方程求解内部参数项的改正项

$\∀i\∈[1,S]:\mathrm{\δ}{\hat{X}}_i=N_{\mathrm{Ai}}^{-1}(d_{\mathrm{Ai}}-N_i^T\mathrm{\δ}\hat{Y})$    式(9)

代入式(8)法方程系统中的第二个方程得到仅与公共参数项的改正项相关的线性方程组：

$\left(\underset{i=1}{\overset{S}{\mathrm{\Σ}}}(N_{\mathrm{Bi}}-N_i{N}_{\mathrm{Ai}}^{-1}{N}_i^T)\right)\mathrm{\δ}\hat{Y}=\left(\underset{i=1}{\overset{S}{\mathrm{\Σ}}}(d_{\mathrm{Bi}}-N_i{N}_{\mathrm{Ai}}^{-1}{d}_{\mathrm{Ai}})\right)$    式(10)

设 $(N_{\mathrm{Bi}}-N_i{N}_{\mathrm{Ai}}^{-1}{N}_i^T)=N_i^{\′},$ $(d_{\mathrm{Bi}}-N_i{N}_{\mathrm{Ai}}^{-1}{d}_{\mathrm{Ai}})=U_i^{\′},$ 则式(10)简写为：

$\underset{i=1}{\overset{S}{\mathrm{\Σ}}}\left(N_i^{\′}\right)\mathrm{\δ}\hat{Y}=\underset{i=1}{\overset{S}{\mathrm{\Σ}}}\left(U_i^{\′}\right)$    式(11)

称N′_i和U′_i分别是第i分区对求解公共参数项的改正项的贡献矩阵和贡献向量，则各个分区的子计算任务可同时消去各个分区的内部参数项的改正项得到各个子计算任务对求解公共参数项的改正项的贡献矩阵N′_i和贡献向量U′_i，累加后求解线性方程组(11)得到改正项和评价其精度的协方差矩阵求解公式如下：

$\mathrm{\δ}\hat{Y}={\left(\underset{i=1}{\overset{S}{\mathrm{\Σ}}}\left(N_i^{\′}\right)\right)}^{-1}\left(\underset{i=1}{\overset{S}{\mathrm{\Σ}}}\left(U_i^{\′}\right)\right)$    式(12)

$Q_{\mathrm{\δ}\hat{Y}}={\left(\underset{i=1}{\overset{S}{\mathrm{\Σ}}}\left(N_i^{\′}\right)\right)}^{-1}$

通过将求解的公共参数项的改正项代入式(9)中，即可得到各个内部参数项的改正项其对应的协方差矩阵为这里的协方差矩阵和反应出各自改正项的求解精度，矩阵对角线上的元素越小，表明该改正项的精度越高；

步骤4：建立平差计算服务,利用Web Service分布式计算技术编写步骤1到步骤3中涉及的大地控制网平差的子计算服务、求解公共参数项的改正项服务的平差计算服务，并发布在网络中的多台计算机上，其中子计算服务能够在本地计算机上根据观测数据建立线性化误差方程(4)，计算其对求解公共参数项的改正项的贡献矩阵N′_i和贡献向量U′_i，通过Web服务的形式对外提供贡献矩阵N′_i和贡献向量U′_i的结果服务，求解公共参数项的改正项服务提供对各子计算服务得到的对求解公共参数项的改正项的贡献矩阵N′_i和贡献向量U′_i进行累加求和、矩阵求逆、矩阵与向量相乘，对外提供公共参数项的改正项和其协方差矩阵

步骤5：平差计算服务分布式的实现，大地控制网平差用户端或控制端计算机分布并行调用网络中多个子计算服务，得到各个子计算服务的计算结果，即对求解公共参数项的改正项的贡献矩阵N′_i和贡献向量U′_i，通过网络交互，调用求解公共参数项的改正项服务得到公共参数项的改正项和评价其精度的协方差矩阵通过广播公共参数项的改正项到网络中，每个本地计算机利用公共参数项的改正项计算得到内部参数项的改正项和评价其精度的协方差矩阵

步骤6：对解算结果进行迭代优化，将每次得到的大地控制网的公共参数项的改正项和所有内部参数项的改正项作为近似值带入到误差方程(2)中，重复步骤(2)到(5)，经过迭代计算，直到所有未知参数的改正项达到求解精度要求。

2.根据权利要求1所述的规模化大地控制网平差的分布式计算方法，其特征在于，所述的获取的与大地控制网中空间分布的控制点位置信息相关的大量的距离、角度、高差、基线观测数据，对大量的距离、角度、高差和基线观测数据采用分区平差方法进行分解后，控制网平差计算的未知参数项分为公共参数项和内部参数项两类，所述的公共参数项由所有的公共点的三维坐标参数构成，内部参数项由各分区内部点的三维坐标参数、尺度参数、旋转因子参数构成。

3.根据权利要求1所述的规模化大地控制网平差的分布式计算方法，其特征在于，所述的同时消去各个分区的内部参数项的改正项即能够在分布式计算环境下并行消去通过并行计算各分区对公共参数项的改正项的贡献矩阵N′_i和贡献向量U′_i，得到仅与公共参数项的改正项相关的线性方程组(11)，然后通过求解线性方程组得到大地控制网的公共参数项的改正项和其协方差矩阵，进而各个分区利用公共参数项的改正项，并行得到各个分区的内部参数项的改正项和其协方差矩阵。

4.根据权利要求1所述的规模化大地控制网平差的分布式计算方法，其特征在于，所述的大地控制网平差用户端或控制端计算机分布并行调用网络中多个子计算服务，是指大地控制网平差计算任务分解后，通过分布式计算技术建立控制网平差子计算服务、求解公共参数项的改正项服务的平差计算服务，并将平差计算服务发布到网络中多个计算节点，大地控制网平差用户端或控制端通过网络首先绑定计算节点，进而并行调用多个子计算服务，实现网络中多台计算机分布协同完成规模化大地控制网平差的分布计算。