Fast Approach to Constructing Normal Equation During the Time of Calculating Height Anomaly Difference by Using Vertical Deflections
-
摘要: 利用垂线偏差等重力格网数据平差计算高程异常差时,施加少量GPS/水准点进行控制,可以确定区域似大地水准面,但是采用传统方法在构造法方程时,需要对系数阵的每个元素逐一进行操作,并全部或者对角存储系数阵,具有计算速度慢、占用内存高等问题。为此提出了在平差解算中对系数阵先进行矩阵分块(操作单元为分块矩阵),再稀疏化处理(仅存非零元素),最后拼接的方法,实现了法方程阵的快速构建及解算。实验表明,相比于传统方法,该方法的计算效率提高了至少两个数量级,并且可快速解算传统方法在一般计算机上难以解算的平差问题,对于解算比较规则的格网数据平差问题具有一定的参考与借鉴意义。Abstract: Vertical deflections can be used tocalculate height anomaly differences, and determine a regional quasi-geoid under the control of GPS leveling points. But the existing typical methods used for structural equation, the elements in coefficient matrix need to be operated one by one, and all of the elements or those on a diagonal must be saved in internal memory for computing. This usually leads to low-speed processing and high-memory occupancy rates. In order to solve these problems, we propose a novel method in the form of a blocked matrix , then processes sparsely (only saving non-zero elements), and finally integrates all the blocked matrixes. This blocked matrix is considered as operational unit in the further steps. Experiments show that: compared to the typical methods, the computation efficiency of our method has been improved at least 2 orders of magnitudes, this makes our method can quickly solve an adjustment problem intractable by traditional methods of computer configuration. Therefore, our method has great reference value when resolving adjustment problems with regular gridded data.
-
利用格网垂线偏差基于天文重力水准的方法计算相邻结点的高程异常差,施加一定数量的GPS/水准点进行控制，求得结点的高程异常，可以确定区域似大地水准面[1-8]。平差计算中涉及到误差方程系数阵A的生成、法方程系数阵N的计算，对于范围4°×6°、分辨率1′×1′的格网垂线偏差数据，平差计算时，系数阵A、N多达十几GB，其规模已超出一般计算机的内存配置。传统方法在构建系数阵时，一般是对系数阵中的每个元素进行操作，多以对角阵的方式存储[9]，随着数据量的增加，生成系数阵的效率迅速降低。从数值分析的角度入手，如Gauss-Seidel迭代法、Gauss-Jordan消去法、Cholesky分解法等，可快速解算法方程，但在计算中仍需要存储大规模数据[10]。Ge提出一种新的参数消去方法[11]，Blewitt提出了Ambizap算法[12]，陈宪东对Ambizap做了验证[13,14]。然而，以上方法主要针对处理大规模的GNSS观测数据[13-17]。宋力杰提出了一种用解算约化代替求逆约化的超大规模大地网分区平差快速解算方法[18]，对系数阵进行压缩存储，但依然是对每个元素进行操作。分布式计算的兴起，有效提高了计算机的利用率[14]，但是对计算平台的性能与数量仍然有一定的要求。由于格网数据排列比较规则，因此平差计算构建的系数阵存在较强的规律性。从这种规律性出发，把操作单元变成分块矩阵并且仅存储矩阵中的非零元素，之后再拼接成系数阵，不仅可以减少对计算机内存的占用率，而且可以显著提高计算速度[19-21]。
1 格网数据平差模型的分析
1.1 格网垂线偏差计算高程异常差的基本原理
本文直接给出利用垂线偏差等数据计算相邻结点高程异常差的方法[4]，如图1所示。
![]()
图 1 利用垂线偏差计算格网高程异常差Figure 1. Using Deflection of the Vertical to Calculate Height Anomaly Difference考虑到重力异常项的影响，计算子午方向和卯酉方向高程异常差的公式分别为：
(1) 
(2) 其中，ζ中、η中分别为相邻两结点子午方向、卯酉方向垂线偏差的平均值；Δg中为两结点重力异常的平均值；Δs为两结点间的距离；为平均正常重力[1]，可以取9.797 644 656 m/s2；Δh为两结点间的海拔高差。
1.2 计算高程异常差的平差模型
利用式(1)与式(2)计算高程异常差近似于水准网中的高差，施加一定的控制点可以采用平差的方式计算格网中每个结点的高程异常[4]。为便于编程实现，一般采用参数平差的方式。
对于m×n格网的每一行，误差方程为：
(3) 其中，Vi,jh为由垂线偏差等数据计算的卯酉方向高程异常差的残差；
对于相邻两行之间的每一列，误差方程为：
(4) 其中，Vi,jl为由垂线偏差等数据计算的子午方向高程异常差的残差；Li,jl为由垂线偏差等数据计算的子午方向相邻两结点的高程异常差。
组成的误差方程为：
(5) 其中，i=1,2，…,m；j=1,2，…,n。
简化为：
(6) 其中，A为误差方程系数矩阵。
相应的法方程为：
(7) 其中，N=ATPA；U=ATPl；P为观测量的权阵[22]。
以上为自由网平差模型，缺少基准，因此需要施加一定数量的GPS/水准点进行控制，才能计算出每个结点的高程异常。为不破坏系数阵的规律性以及易于编程实现，本文将控制点作为“未知点”，但是赋予其较大权值(本文取1030)，即在法方程系数阵N中与控制点相对应的元素加上一个大的数[9]，从而使得在平差计算中控制点的值不会产生变化。
为提高计算速度与精度，在编程时本文求解
(8) 其中,“\\”表示Matlab中的“左除”，用消去的方法解算法方程。
2 平差法方程的构建
2.1 系数阵的特征分析
设误差方程系数阵A、法方程系数阵N的分布形式分别如表1、表2所示。
表 1 3×4格网数据误差方程系数阵ATable 1. Coefficient Matrix A of 3×4 Grid Data Error Equation
表 2 4×4格网数据法方程系数阵NTable 2. Coefficient Matrix N of 4×4 Grid Data Normal Equation
通过对表1、表2系数阵进行分析，不难发现A、N的如下特征(设平差区域格网数据为m行×n列)。
对于误差方程系数阵A,有：(1) 不同颜色分块矩阵的行列数与格网数据的行列数有关。表1横线上半部分，浅色分块矩阵为m个(n－1)×n的矩阵且元素相同。而下半部分，两种颜色的矩阵分块个数是m－1，为n×n的矩阵；(2) 每一行只有两个元素(1和-1)；(3) 整个系数阵、每个分块矩阵都为稀疏矩阵。
对于法方程系数阵N，N=ATPA，取P=1，有：(1) 主对角线上，首尾(浅绿色)的矩阵块完全相同，中间的浅色部分完全相同，浅色部分的个数为m－2，每个分块矩阵大小都是n×n；(2) 主对角线上，首尾矩阵块的首尾元素都为2，中间n－2个元素为3，中间浅色部分矩阵块首尾元素为3，中间n－2个元素为4；(3) 主对角线两侧分块矩阵都是大小为n×n、主对角线元素是-1的对角阵；(4) 整个系数阵和每个分块矩阵都是稀疏矩阵。
2.2 系数阵的构建方法介绍
根据以上特点，本文提出了分块构造,稀疏化处理,最后拼接成系数阵的方法。系数阵的构建过程如下。
2.2.1 误差方程系数阵A的构建过程
(1) 生成一个(n－1)×(n－1)的单位阵
(2) 构造
(3) 生成一个n×n的单位阵
(4) 拼接(2)中生成的分块矩阵，构成系数阵A的上半部分A1：
(9) (5) 拼接(3)中生成的分块矩阵，构成系数阵A的下半部分A2：
(10) (6) 拼接A1和A2，构成误差方程系数阵A：
(11) 2.2.2 法方程系数阵N的构建过程
首先，设对角线上首尾的矩阵块为
(1) 分别生成分块矩阵
(2) 重新拼接
(12) (3) 拼接(2)中分块矩阵构成法方程系数阵N：
(13) 至此，误差方程系数阵A和法方程系数阵N全部构建完毕。
对于格网数据不规则、数据缺失、存在不合理的点(粗差点)的情况，可以添加一定的数据，构造成规则的格网，将法方程系数矩阵N中与参数相对应的元素置为0(即添加的数据权为0，不参与计算)，从而可以利用本文方法完成平差计算。
3 实验计算比较分析
本文选取分辨率相同、范围不同的数据，在处理器为单核E5800 (3.2 GHz)，内存2 GB的联想台式机上采用Matlab R2011a平台[20]进行计算。
采用自由网平差方法分别解算范围为0.5°×0.5°、1°×1.5°、1.5°×1.5°、4°×6°，分辨率为1′×1′的区域，即未知参数的个数分别为900、5 400、8 100、86 400个，采用不同的方法生成误差方程系数阵A、法方程系数阵N，计算参数平差值，耗时统计结果如表3所示(表中“-”表示数据量超出计算机内存，无法继续计算)。
表 3 不同方法计算不同区域耗时/sTable 3. Computing Time of Different Methods/s区域 内容 传统 方法 对A 稀疏化 本文方法 直接 生成N 由A生 成N A 0.045 0 0.017 6 0.5°×0.5° N 0.195 9 0.058 1 0.006 3 0.018 2 0.349 0 0.106 9 0.046 8 0.041 3 A 2.546 9 0.063 0 1°×1.5° N 30.497 8 3.026 4 0.033 7 0.064 6 37.907 0 3.126 3 0.114 1 0.146 5 A 10.580 6 0.125 8 1.5°×1.5° N 104.598 0 11.317 3 0.057 4 0.134 5 - 11.655 9 0.318 3 0.281 1 A - - 4.451 4 4°×6° N - - 2.883 9 4.467 3 - - 8.700 7 5.769 5 由表3可以看出，随着数据量的增大，采用传统方法的耗时迅速增加，当计算范围大到一定的程度时，已无法完成法方程的构建与解算。分析表可以发现：在传统方法的基础上，对误差方程系数阵A做稀疏化处理，可迅速提高构造系数阵的速度。当计算范围扩大到4°×6°时，即解算86 400个未知参数时，该方法已经不能满足计算需要，这是因为初步估计用传统的方法生成的系数阵A将会占用13.86 GB左右的空间，其规模已超出一般计算机内存的容量(8 GB内存，2014)，因此无法完成参数的解算。采用本文提出先矩阵分块，再稀疏化处理，最后拼接的方法，可快速提高解算速度。这是因为采用分块的方法，在生成系数矩阵时，是对小块的矩阵进行操作，减少了循环计算的次数；而对分块矩阵稀疏化处理以后，计算机仅存储分块矩阵中的非零元素，节约了大量的内存空间；而由稀疏化处理后的分块矩阵拼接成完整的系数阵时，计算机操作的是稀疏化处理后的分块矩阵，所生成的完整的系数阵同样是稀疏化处理之后的矩阵，因此即便是生成一个拥有172 200×86 400个元素的误差方程系数阵A，仍然不会超出一般计算机的内存限制。在计算过程中主要是稀疏矩阵之间的运算，因而解算速度得到了很大提高。无论是直接生成法方程系数阵N，还是由A计算生成，都能够快速完成法方程阵的构建和解算。但是由A计算生成更好，这是因为直接生成N，要调用生成N的子函数，同时随着N的规模的增大，生成N所耗费的时间要大于计算AT PA生成N的时间。当考虑观测量权矩阵P时，法方程系数阵N的规律也不会太明显。由误差方程系数阵A计算生成N也更符合法方程构造过程，并且适用范围也会更广。
4 结 语
针对利用垂线偏差等重力格网数据计算高程异常差时法方程规模巨大、解算速度慢等问题，本文提出了先矩阵分块，再稀疏化处理，最后拼接的方法，不仅减小了计算机内存的占用率，而且可快速实现法方程阵的构建与解算等。其具体优势在于，一方面将传统的对系数阵中每个元素进行操作的方式改为对分块的矩阵进行操作，减少计算时的循环次数；另一方面采用了稀疏化处理的方法，仅存矩阵中非零元素，对于一个0元素占主体的矩阵，这无疑缩小了矩阵的存储规模。这两方面的原因使得在一般计算机环境下原本需要耗费大量时间甚至是难以完成解算的高程异常差平差问题可以快速解决。本文提出的方法对于快速解算大规模的重力场扰动场元网数据以及大规模GPS网观测数据具有一定的借鉴和指导意义。
-
![]()
图 1 利用垂线偏差计算格网高程异常差
Figure 1. Using Deflection of the Vertical to Calculate Height Anomaly Difference
表 3 不同方法计算不同区域耗时/s
Table 3 Computing Time of Different Methods/s
区域 内容 传统 方法 对A 稀疏化 本文方法 直接 生成N 由A生 成N A 0.045 0 0.017 6 0.5°×0.5° N 0.195 9 0.058 1 0.006 3 0.018 2 0.349 0 0.106 9 0.046 8 0.041 3 A 2.546 9 0.063 0 1°×1.5° N 30.497 8 3.026 4 0.033 7 0.064 6 37.907 0 3.126 3 0.114 1 0.146 5 A 10.580 6 0.125 8 1.5°×1.5° N 104.598 0 11.317 3 0.057 4 0.134 5 - 11.655 9 0.318 3 0.281 1 A - - 4.451 4 4°×6° N - - 2.883 9 4.467 3 - - 8.700 7 5.769 5 
下载: 导出CSV
-
[1] 李建成,陈俊勇,宁津生,等.地球重力场逼近理论与中国2000似大地水准面的确定[M].武汉:武汉大学出版社,2003 Li Jiancheng, Chen Junyong, Ning Jinsheng, et al. The Earth's Gravity Field Approximation Theory and CQG 2000[M]. Wuhan:Wuhan University Press, 2003
[2] 李建成.我国现代高程测定关键技术若干问题的研究及进展[J].武汉大学学报·信息科学版,2007(11):980-987 Li Jiancheng. Study and Progress in Theories and Crucial Techniques of Modern Height Measurement in China[J]. Geomatics and Information Science of Wuhan University, 2007(11):980-987
[3] 朱雷鸣.基于相对高程异常差平差法的区域似大地水准面精化[D].郑州:信息工程大学,2011 Zhu Leiming. Refining of Quasi-Geoid Based on the Relative Height Anomaly Difference Adjustment[D]. Zhengzhou:Information Engineering University, 2011
[4] 李姗姗,吴晓平,张传定,等.顾及地形与完全球面布格异常梯度项改正的区域似大地水准面精化[J].测绘学报,2012,41(4):510-516 Li Shanshan, Wu Xiaoping, Zhang Chuanding, et al. Regional Quasi-Geoid Refining Considering Corrections of Terrain and Complete Spherical Bouguer Anomaly's Gradient Term[J]. Acta Geodaetica et Cartographica Sinica, 2012, 41(4):510-516
[5] Huang J,Veronneau M. Canadian Gravimetric Geoid Model 2010[J]. J Geod, 2013, 87:771-790
[6] Andrea Gatti, Mirko Reguzzoni, Venuti G, et al. The Height Datum Problem and the Role of Satellite Gravity Models[J]. J Geod, 2013, 87:15-22
[7] Klees R, Prutkin I. The Combination of GNSS-Levelling Data and Gravimetric (Quasi-) Geoid Heights in the Presence of Noise[J]. J Geod, 2010, 84:731-749
[8] 郭东美,许厚泽.应用GPS水准与重力数据联合解算大地水准面[J].武汉大学学报·信息科学版,2011,36(5):621-624 Guo Dongmei, Xu Houze. Determination of Geoid Using GPS Leveling and Gravity Data[J]. Geomatics and Information Science of Wuhan University, 2011, 36(5):621-624
[9] 宋力杰. 测量平差程序设计[M].北京:国防工业出版社,2009 Song Lijie. Program Design of Surveying Adjustment[M]. Beijing:National Defense Industry Press, 2009
[10] 郭飞霄,杨力,刘荣,等.大型稀疏法方程组的代数多重网格解法[J].测绘科学技术学报,2012,29(1):5-8 Guo Feixiao, Yang Li, Liu Rong, et al. Algebraic Multigrid Solution of Large-Scale Sparse Normal Equation[J]. Journal of Geomatics Science and Technology, 2012, 29(1):5-8
[11] Ge M, Gendt G, Dick G. A New Data Processing Strategy for Huge GNSS Global Networks[J]. J Geod, 2006, 80:199-203
[12] Blewitt, Geoffery. Fixed Point Theorems of GPS Carrier Phase Ambiguity Resolution and Their Application to Massive Network Processing:Ambizap[J]. Journal of Geophysical Research, 2008, 113(B12):36-44
[13] 陈宪冬.Ambizap方法在大规模GPS网处理中的应用及结果分析[J].武汉大学学报·信息科学版,2011,36(1):10-13 Chen Xiandong. Application of Ambizap Algorithm in Large GPS Network and Its Test Results[J]. Geomatics and Information Science of Wuhan University, 2011, 36(1):10-13
[14] 崔阳.大规模测量平差分布式计算技术及应[D].郑州:信息工程大学,2013 Cui Yang. Research on the Theory and Application of Distributed Computing for Large-Scale Survey Adjustment[D]. Zhengzhou:Information Engineering University, 2013
[15] 成英燕,程鹏飞,顾旦生,等.天文大地网与GPS2000网联合平差数据处理方法[J].武汉大学学报·信息科学版,2007,32(2):148-151 Cheng Yingyan, Cheng Pengfei, Gu Dansheng, et al. Data Processing Method for Combined Adjustment of National Astro-Geodetic Network and GPS2000 Network[J]. Geomatics and Information Science of Wuhan University, 2007, 32(2):148-151
[16] 姜卫平,赵倩,刘鸿飞,等.子网划分在大规模GNSS基准站网数据处理中的应用[J].武汉大学学报·信息科学版,2011,36(4):389-391 Jiang Weiping, Zhao Qian, Liu Hongfei, et al. Application of the Sub-network Division in Large Scale GNSS Reference Station Network[J]. Geomatics and Information Science of Wuhan University, 2011, 36(4):389-391
[17] 陈俊平,张益泽,谢益炳,等.超大观测网络及多GNSS系统的快速数据处理[J].武汉大学学报·信息科学版,2014,39(3):253-257 Chen Junping, Zhang Yize, Xie Yibing, et al. Rapid Data Processing of Huge Networks and Multi-GNSS Constellation[J]. Geomatics and Information Science of Wuhan University, 2014, 39(3):253-257
[18] 宋力杰,欧阳桂崇.超大规模大地网分区平差快速解算方法[J].测绘学报,2003,32(3):204-207 Song Lijie, Ouyang Guichong. A Fast Method of Solving Partitioned Adjustment for Super Large-Scale Geodetic Network[J]. Acta Geodaetica et Cartographica Sinica, 2003, 32(3):204-207
[19] 侯秋果.矩阵分块的应用[J].科技信息,2008,13(2):190-191 Hou Qiuguo. Application of Blocked Matrix[J]. Science & Technology Information, 2008, 13(2):190-191
[20] 张志涌.精通MATLAB R2011a[M].北京:北京航空航天大学出版社,2011 Zhang Zhiyong. Master MATLAB R2011a[M]. Beijing:Beihang University Press, 2011
[21] 巩学美,魏代勇.大型稀疏正定法方程非零动态存储三角分解法[J].测绘学报,2000,29(3):209-215 Gong Xuemei, Wei Daiyong. Non-zero Element Dynamic Memory Trigonometric Decomposition Method of Large-Scale Sparse Positive Definite Normal Equation[J]. Acta Geodaetica et Cartographica Sinica, 2000, 29(3):209-215
[22] 隋立芬,宋力杰,柴洪洲.误差理论与测量平差基础[M].北京:测绘出版社,2010 Sui Lifen, Song Lijie, Chai Hongzhou. Error Theory and Foundation of Surveying Adjustment[M]. Beijing:Surveying and Mapping Press, 2010
-
期刊类型引用(1)
1. 吴静，傅优杰，程朋根. 基于粗糙集的局部同位模式挖掘算法. 测绘通报. 2022(10): 80-85+104 . 
百度学术
其他类型引用(3)
计量
- 文章访问数: 1608
- HTML全文浏览量: 63
- PDF下载量: 339
- 被引次数: 4
