2. 云南省地震局,昆明市知春街249号,650093;
3. 河南测绘职业学院,郑州市工贸路,451464
流动重力数据的采集一般利用2台重力仪在同一测点进行同步、往返观测,获得2组野外观测数据(CH/T 2003-1999, GB/T 20256-2006),由于观测仪器、观测人员及观测环境等的影响,2组观测数据的精度存在不同。在进行平差处理时,以绝对点上观测的重力值作为基准,控制整个区域的相对点,通过等权、不等权等方法对2组观测数据进行处理,获得测点的相对重力数据。然后对整个测区的重力变化情况进行对比分析,得到重力异常的大小及区域变化情况。不同定权方法的计算结果不同,根据2组观测数据及其权值的数学关系,将观测值及对应测点的绝对重力值作为已知量,各自的权值作为未知量,并组成方程组,超过2组观测量组成超定方程组。最后利用最小二乘方法解算超定方程组[1],计算最优权值,再利用最优权值进行平差处理,得到所有测点最接近绝对重力值的结果。该方法可为地球内部物质的研究及地震监测预报分析提供更为真实可靠的观测数据。
1 测点情况及计算公式云南省目前有相对重力观测点240多个,绝对重力观测点8个,绝对观测点分布均匀,不存在平差计算中短边控制长边或部分测区无绝对点控制的情况。相对观测点与绝对观测点联测,并利用绝对观测点对整个区域进行平差控制,保证重力数据的正确性。云南省重力相对观测点及绝对观测点的分布情况见图 1(联测四川省的乡城测点)。
流动重力观测严格遵守野外测量规范,使用2台相对重力仪进行同步、同精度观测,并按一定的周期重复观测测点之间的重力段差值,对同一测段往返观测,最后对野外观测数据进行固体潮、零漂和格值等改正,与绝对观测点进行整体平差,得到每个测点的重力值。具体计算公式为:
$ {a_1}{p_1} + {a_2}{p_2} = b $ | (1) |
式中,a1、a2为测点A经过固体潮、零漂和格值等改正后2台仪器各自的观测值,P1、P2为2台仪器各自的权值,b为测点A的重力值。
2 不等权和等权计算结果云南地区的流动重力观测工作一般在每年的3~6月及7~10月进行,并将绝对观测点联测到整个观测网中,与相对观测点进行同精度、同步观测。2014年起,使用2台CG-5型相对重力仪进行观测,编号分别为C1169和C1170。多期观测结果表明,C1170在野外观测时稳定性较好且数据精度较高[2]。云南地区的绝对重力观测工作一般在每年的5月和6月进行,由中国地震局地震研究所统一进行观测,观测仪器为A10重力仪,精度为1~2 μGal。本文选取2018年下半年及2019年下半年的观测数据,分别进行等权及不等权计算[3-4](根据2台仪器的观测精度定权),并对比2种定权方式计算结果与绝对值的差值,结果见表 1(单位μGal),2种定权方式计算结果的精度及各仪器方差见表 2(单位μGal)。
从表 1和2可以看出,在定权方式上,不等权的计算结果在绝对观测值及精度方面均优于等权的计算结果,在数据计算方面也较好。
3 最小二乘计算结果 3.1 最小二乘原理对m×n(m>n)的线性方程组有:
$ \left\{ \begin{array}{l} {a_{11}}{P_1} + {a_{12}}{P_2} + \ldots + {a_{1n}}{P_n} = {b_1}\\ {a_{21}}{P_1} + {a_{22}}{P_2} + \ldots + {a_{2n}}{P_n} = {b_2}\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \vdots \\ {a_{m1}}{P_1} + {a_{m2}}{P_2} + \ldots + {a_{mn}}{P_n} = {b_m} \end{array} \right. $ | (2) |
即
$ \mathit{\boldsymbol{AP}} = \mathit{\boldsymbol{b}} $ | (3) |
当方程的个数m大于未知量个数n时,方程组即为超定方程组,一般来说是没有解的,即对于任意一组未知数(P1,P2,…,Pn),有:
$ {\delta _i} = \sum\limits_{j = 1}^n {{a_{ij}} - {b_i}} , i = 1, 2, \ldots , m $ | (4) |
其结果不会全为0。假设求一组未知数P*=(P1* P2* … Pn*),使得
$ \begin{array}{l} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \varphi ({P_1}, {P_2}, \ldots , {P_n}) = \\ \sum\limits_{i = 1}^m {\delta _i^2} {\rm{ = }}\sum\limits_{i = 1}^m {{{\left( {\sum\limits_{j = 1}^n {{a_{ij}}{P_j} - {b_i}} } \right)}^2}} \end{array} $ | (5) |
取得极小值,利用多元函数求极值的方法可得:
$ \frac{{\partial \varphi }}{{\partial {x_k}}} = 2\sum\limits_{i = 1}^m {\sum\limits_{j = 1}^n {({a_{ij}}{P_j} - {b_i}){a_{jk}}} } = 0 $ | (6) |
即
$ \sum\limits_{i = 1}^m {\left( {\sum\limits_{j = 1}^n {{a_{ij}}{a_{jk}}} } \right){P_j}} = \sum\limits_{j = 1}^n {{a_{ij}}{b_i}} $ | (7) |
写为矩阵形式:
$ {\mathit{\boldsymbol{A}}^{\text{T}}}\mathit{\boldsymbol{AP}} = {\mathit{\boldsymbol{A}}^{\text{T}}}\mathit{\boldsymbol{b}} $ | (8) |
可以证明,如果A是满秩的,则该方程存在唯一解,且取得最小值,该解即为超定方程组的最小二乘解,也是2组观测仪器的最优权值。
3.2 最小二乘计算结果对于C1169和C1170的流动重力观测而言,n=2,m始终大于n,线性方程组为超定方程组。采用最小二乘方法求解2018年下半年及2019年下半年2个超定方程组,分别得到2台仪器的权值,并将结果代入式(1)用于平差计算,得到绝对观测值的差值,结果见表 3(单位μGal)。
由表 3可知,最小二乘定权方式计算得到的值大部分与不定权结果接近,约67%的结果优于等权的结果,说明最小二乘定权方法在理论上是可取的。
3.3 重力异常与地震的关系分别采用3种不同的定权方式计算2018~2019年4期数据,并分析期间云南省重力异常与MS4.0及以上地震的对应关系,结果见图 2。
由图 2可知,不等权计算得到的重力异常值区间为-47.0~27.6 μGal,而等权计算得到的重力异常值区间为-69.8~70.2 μGal,最小二乘计算得到的重力异常值区间在-62.1~53.5 μGal。图 2(a)中,根据2台仪器的观测精度取不等权的计算方法,重力异常结果的取值范围缩小,排除了正常干扰的同时也去掉了一些因地壳活动造成的重力异常信息,因此利用该方法计算得到的重力异常等值线在反映MS5.0及以上地震的情况时,与等权及最小二乘定权方式基本相同,但对MS4.0~5.0地震,其重力异常等值线较为稀疏,且地震的发生位置没有等权及最小二乘定权方式计算结果中的零值线、四象限等特征那么明显; 图 2(b)和2(c)中的图形基本一致,但对于靠近绝对观测点处的地震,最小二乘定权方式的计算结果中,重力异常等值线的高梯度带现象更明显。
4 结语综合3种定权方式的计算结果及重力异常等值线与地震的对应关系,得出以下结论:
1) 仪器的权值确定侧面反映了仪器本身的性能,与仪器各自的观测精度一致;
2) 不等权的定权方式虽然是根据仪器各自观测精度来定权,但在对其中一台仪器取较大权值时削弱了另一台仪器的观测作用,虽然得到的观测值精度更高,但也削弱了反映MS4.0~5.0地震的有用的重力异常信息;
3) 采用最小二乘定权方式的计算结果与绝对观测点的计算结果一致,从地震监测预报上来说是一种比较好的方法,相比不等权的计算结果,重力异常与MS4.0以上地震的发生有很好的对应关系,相比等权的计算结果,其精度更好,与绝对观测点得到的重力值更接近。
本文最小二乘定权方式的计算结果在精度方面优势不明显,初步估计是由于各个绝对观测点的相对观测时间与绝对观测时间不同步,地球内部物质流动随时都在变化,受各种因素的影响,同一观测点在不同时间观测的结果都不一样。
[1] |
肖云, 王云鹏, 刘晓刚, 等. 空域最小二乘法用于重力卫星误差分析[J]. 武汉大学学报:信息科学版, 2019, 44(3): 340-346 (Xiao Yun, Wang Yunpeng, Liu Xiaogang, et al. Application of Space-Wise Least Square Method to Error Analysis for Satellite Gravimetry[J]. Geomatics and Information Science of Wuhan University, 2019, 44(3): 340-346)
(0) |
[2] |
汪健, 孙少安, 邢乐林, 等. CG-5重力仪的漂移特征[J]. 大地测量与地球动力学, 2016, 36(6): 556-560 (Wang Jian, Sun Shao'an, Xing Lelin, et al. Drift Characteristics of CG-5 Gravimeter[J]. Journal of Geodesy and Geodynamics, 2016, 36(6): 556-560)
(0) |
[3] |
蒋萍, 于宏亮, 王司光. 不等精度测量中权的确定方法[J]. 济南大学学报:自然科学版, 2017, 30(3): 229-232 (Jiang Ping, Yu Hongliang, Wang Siguang. Weights Determining Method of Unequal Precision Measurement[J]. Journal of University of Jinan: Science and Technology, 2017, 30(3): 229-232)
(0) |
[4] |
高冰. 不等精度测量数据处理中的加权原则[J]. 数学理论与应用, 2002, 22(1): 119-122 (Gao Bing. The Principle of Adding Weight in Unequal Precision Measured Data Processing[J]. Mathematical Theory and Application, 2002, 22(1): 119-122)
(0) |
2. Yunnan Earthquake Agency, 249 Zhichun Street, Kunming 650093, China;
3. Henan College of Surveying and Mapping, Gongmao Road, Zhengzhou 451464, China