文章快速检索     高级检索
  大地测量与地球动力学  2023, Vol. 43 Issue (9): 974-979  DOI: 10.14075/j.jgg.2023.09.017

引用本文  

易天阳, 郑兵, 李亚伦, 等. 关于相对重力测量数据处理的探讨[J]. 大地测量与地球动力学, 2023, 43(9): 974-979.
YI Tianyang, ZHENG Bing, LI Yalun, et al. Discussion on Data Processing of Relative Gravity Survey[J]. Journal of Geodesy and Geodynamics, 2023, 43(9): 974-979.

第一作者简介

易天阳,高级工程师,主要从事地震监测研究,E-mail:sunsky-yi@163.com

About the first author

YI Tianyang, senior engineer, majors in earthquake monitoring, E-mail: sunsky-yi@163.com.

文章历史

收稿日期:2022-10-30
关于相对重力测量数据处理的探讨
易天阳1     郑兵1     李亚伦1     温军军1     易松泉1     张稔1     
1. 四川省地震局地壳形变观测中心,四川省雅安市上坝路139号,625000
摘要:按照《地震重力测量规范》相关理论和细则,对目前所采用的“重力测量电子记簿系统”数据处理软件存在、实际工作中急需解决的几个问题,如仪器高改正、段差计算、环闭合差计算等进行讨论。
关键词格值表气压改正高度改正潮汐改正零漂改正

相对重力测量是指测站间的重力差测量,可采用动力法和静力法,现在普遍采用基于静力法的弹簧重力仪测定重力差值。重力测量包括野外数据采集和内业平差计算两个方面,中国地震局各野外测量队外业数据采集所使用的“重力测量电子记簿系统”由中国地震局地震研究所研制,可用于地面重力仪各种观测数据的记录,包括数据采集、限差校核、实时预处理;而外业电子记簿系统导出的交换文件为内业平差计算所需的数据格式,为最终的数据处理作准备[1]

对重力网进行计算之前,必须对重力仪读数进行归化和改正。对与固体潮、极移和大气压力变化有关的重力瞬时变化,可借助模型计算。

就数据采集而言,“重力测量电子记簿系统”所记录的关键信息(表 1)非常完善,但是在数据处理上存在值得讨论的地方。

表 1 数据记录信息 Tab. 1 Data record information
1 手簿工程文件到平差交换文件

根据“重力测量电子记簿系统”平差交换文件(表 2)输出的数据内容,结合参与重力计算的相关观测数据并与表 1比较发现,最终会参与到归化改正的关键数据中缺失仪器高、经纬度及海拔高度。虽然经纬度和海拔高度可根据点号联合点之记文件计算,但仪器高不可或缺,而与归化改正无关的环境温度值却被要求记录。

表 2 平差交换文件 Tab. 2 Documents of adjustment exchange

在一个观测点进行多次(多周期)观测仪器安置时,仪器顶面到观测台的顶面高度不一定完全相同,而且多年后可能由于观测台倾斜或顶面不平整需要添加托盘支架才能水平安置仪器,造成观测周期间同一测点的仪器安置高度不一定相同。另外,同一小组、同一观测周期使用多台仪器进行观察,仪器间的高差也可能不一样。为了进行零漂率计算和不同周期间数据的比较分析,将每一次观测值改算到观测台顶面非常必要,因此量测仪器高度是必不可少的。

2 数据的实时预处理

流动重力野外数据采集的仪器读数并不是重力值,需要利用仪器厂家给定的格值表进行转化。在格值转换(读数值到重力值的计算)过程中,按照《地震重力测量规范》[2]的要求,格值转换后加入一次项系数的计算再进行潮汐改正,零漂改正在潮汐改正之后。而“重力测量电子记簿系统”中测线段差在“改前段差”到“改后段差”这一过程加入一次项系数(表 3,单位mGal)的计算,即是在“改后段差”中将潮汐改正、气压改正、高度改正、零漂改正这几项均纳入了一次项系数的计算。这种做法是不合理的,因为仪器的一次项系数仅与观测仪器有关,而与潮汐、气压、高度、零漂改正无关。例如已知如下数据:1)仪器号:G1149;零漂值:0.000 00; 零漂率:-0.006 917; 2)一次项系数:1.000 038;3)格值表(表 4, 单位mGal),可结合表 1观测数据完成相关计算。

表 3 观测手簿完成(往返同天)的重力联测计算 Tab. 3 Calculation of gravity joint survey completed in one day

表 4 格值 Tab. 4 Grid value table
2.1 潮汐改正

计算重力固体潮汐理论值时,应使用测点经、纬度值并结合当前日期、时间,具体可采用直接法或分波叠加法。采用直接法较为简单,从未作谐波展开的引潮位公式出发,计算重力固体潮理论值。计算公式为[3]

$ \begin{aligned} g= & -165.7 F(\varphi)\left(\frac{C}{R}\right)^3\left(\cos ^2 Z-\frac{1}{3}\right)- \\ & 1.37 F^2(\varphi)\left(\frac{C}{R}\right)^4\left(5 \cos ^2 Z-3\right)- \\ & 76.08 F(\varphi)\left(\frac{C_{\mathrm{s}}}{R_{\mathrm{s}}}\right)^3\left(\cos ^2 Z_{\mathrm{s}}-\frac{1}{3}\right) \end{aligned} $ (1)

式中,F(φ)=0.998 327+0.001 676 cos2φφ为计算点的地理纬度;Z为计算时刻月亮地心天顶距;Zs为计算时刻太阳地心天顶距;$\frac{C}{R}$为月地平均距离和瞬时距离之比;$\frac{C_{\mathrm{s}}}{R_{\mathrm{s}}}$为日地平均距离和瞬时距离之比。为计算g,必须使用5个有关的天文参数:月亮平黄经S,太阳平黄经H,月亮近地点经度p,月亮升交点经度N,近日点经度Ps,这5个参数分别有各自的涉及儒略世纪数T的计算公式,具体见文献[3]。

表 1中51052503测点12:01时刻的潮汐改正值依据式(1)计算为0.130 8 mGal(表 5,单位mGal),与表 3潮汐改正数0.129 6 mGal略有差异,可能是由于“重力测量电子记簿系统”所用算法不是直接算法而造成的。

表 5 重力观测(往返同天)数据改正计算 Tab. 5 Correction calculation of gravity observation data completed in one day
2.2 仪器读数转重力值

依据仪器制造厂家给定的仪器格值表(单位10 μGal)将仪器读数转换为相应的重力值。重力值换算按下式进行[3]

$ g=\left[F_1+\left(R-R_1\right) \times F_2\right] \times F-\mathrm{TD} $ (2)

式中,g为相应仪器读数的重力值,F1为仪器读数每100个分划间隔的整数所对应的重力值,R为仪器读数,R1为仪器读数以100分划为间隔的最大整数,F2为仪器格值间隔因子,F为仪器格值一次项系数,TD为潮汐改正值。

表 1中51052503测点12:01时刻读数转换重力值按式(2)的计算为:

$ \begin{gathered} (2248.788+(2291.4170-2200.00) \times \\ 1.02258) \times 1.000038+0.1308= \\ 2342.4890(\mathrm{mGal}) \end{gathered} $

表 3中数据转换重力值计算为:

$ \begin{gathered} (2248.788+(2291.4170-2200.00) \times \\ \text { 1. } 02258)=2342.2692(\mathrm{mGal}) \end{gathered} $
2.3 气压改正

大气压力变化的直接和间接效应会引起重力变化。考虑到测点上实际大气压p(气压计测量的精度为±1 hPa)和正常大气压pn之间的差异,应采用下式进行归一化[3]

$ \delta g_p=3\left(p-p_{\mathrm{n}}\right)_{(\mathrm{hPa})} \mathrm{nms}^{-2} $ (3)

式中,pn为根据测站高度H计算的标准大气压值:

$ p_{\mathrm{n}}=1013.25\left(1-\frac{0.0065 H_{(m)}}{288.15}\right)^{5.2559} \mathrm{hPa} $ (4)

表 1中51052503测点12:01时刻按式(4)进行气压改正:

$ \begin{gathered} \delta g_p=3(848-1013.25 \times \\ \left.\left(1-\frac{0.0065 \times 1561.60}{288.15}\right)^{5.2559}\right) \mathrm{nms}^{-2}= \\ 26.43 \mathrm{nms}^{-2}=0.0026 \mathrm{mGal} \end{gathered} $

表 3表 5的气压改算结果一致。

2.4 高度改正

正常重力和高程的关系通常使用垂向导数描述,展开到f阶并归化到φ=45°为[3]

$ \begin{gathered} \left(\frac{\partial \gamma}{\partial h}\right)_0=-3086 \mathrm{nms}^{-2} / \mathrm{m}= \\ -3.086 \mathrm{mGal} / \mathrm{m} \end{gathered} $ (5)

表 1中51052503测点12:01时刻按式(5)进行高度改正:

$ \begin{gathered} -3086 \mathrm{nms}^{-2} \times \frac{210 \mathrm{~mm}}{1000 \mathrm{~mm}}= \\ -648 \mathrm{mms}^{-2}=-0.0648 \mathrm{mGal} \end{gathered} $

表 4表 5的高度改算结果一致。

2.5 零漂改正

要完成零漂改正,首先需要完成重力仪漂移率的计算,而重力漂移率的计算方法可采用全网最或然漂移率,亦可求取测线拟合漂移率。目前较为常用的是对测线进行漂移率计算,具体公式为[3]

$ K=\frac{\sum \mathrm{d} R_i \times \mathrm{d} t_i}{\sum \mathrm{d} t_i \times \mathrm{d} t_i}=\frac{\sum\left(g_i^{\prime}-g_i\right)\left(t_i^{\prime}-t_i\right)}{\sum\left(t_i^{\prime}-t_i\right)^2} $ (6)

式中,K为漂移率,dRi为第i点经固体潮改正的往返观测重力值之差,dti为第i点往返观测时间之差。

3 计算实例

按照规范要求,经过固体潮改正后进行漂移率的计算,但往返测各测点的重力值不应该是仅进行了固体潮和零漂改正之后的重力值。因为流动周期观测任务中,同一测点在不同时间进行观测,其气压值和仪器高度不一定完全相同,各测点的重力值还应归算到测点同一高度和经过气压改正之后的重力值,这样往返测量结果的较差才有意义。

表 5是对一个测段且往返测在同一天完成的简单重力测量线路进行零漂计算,计算的段差结果与表 3一致。虽然这种做法不会对测点间的段差大小产生较大影响,但会明显影响测点的重力值计算(表 6,单位mGal)。如果测线涉及多个测段在多天完成,或同一天同一测点二次设站,或线路中较常规计划线路有支线或新增测点,则线路拟合零漂率的计算较为复杂。

表 6 段差(往返同天)计算比较 Tab. 6 Comparison of segment difference calculations completed in one day

以往返测不同天完成的测线为例,表 7(单位mGal)中51026400测点结束于11日下午16:01,12日上午8:23接着该点出发,2次读数肯定不同,经过格值转换的重力值也不一定完全相等,即便经过潮汐改正、气压改正、高度改正,重力值也不一定完全相等。但是该点这2 d的重力值理论上应该相等,因此将第2天的开始时间归算到前1天的结束时间,将第2天的开始重力值归算到前一天的结束重力值,在当日起点归算重力值的基础上加上当日测线上测点初测重力值与当天第1个测点初测重力值之差获得后续测点的归一重力值,在当日起点时间基础上加上测点当日时间与当天第1个测点起测时间之差获得归一时间,具体见表 8

表 7 观测手簿完成(往返不同天)的重力联测计算 Tab. 7 Calculation of gravity joint survey completed in different days

表 8 重力观测(往返不同天)数据改正计算表模板 Tab. 8 Template for correction calculation of gravity observation data completed in different days

表 8为同一测点Ni的2次(同天或不同天观测)设站情况。一般而言,2次设站读数值经过潮汐改正后的重力值不一定相等,而理论上同一个测点的重力值应该相等,因此为方便计算漂移率,需要将2次设站时间和重力值统一到上一次设站,即⑤=④,Gi″=Gi′,而后续测点的时间(注意将时分转换为时小数计算):⑥=⑤+③-②,重力值:Gi+1=Gi′+G3+(7)+(8)+(9)-G2-(4)-(5)-(6)。

依据表 8所描述的计算方法,表 9(单位mGal)中51034001测点返测相关计算如下:

表 9 重力观测(往返不同天)数据改正计算 Tab. 9 Correction calculation Table. of gravity observation data completed in different days

1) 51034001测点返测观测的归一时顺时间为:16.016 7+9+13/60-8-23/60=16.850 0(h)。

2) 51034001测点返测观测的归一重力为:2 060.464 5+2 066.387 7-0.034 3+0.003 6+0.066 3-2 060.443 3+0.062 3-0.004 5-0.066 3=2 066.436 0(mGal)。

依据式(6),表 9中的漂移率K值计算为:K=(2.483 4×(-0.011 1)+4.133 3×(-0.000 5))/(2.483 42+4.133 32)=-0.001 27(mGal/) h。将K值再乘以各个观测点与第1个观测点时间差,得到各观测点的零漂改正,最后得到各观测点的改正后重力值。

先期纳入一次项系数与后期纳入一次项系数的测点重力值、测点间段差的结果有没有差别是不同方法应该比较的问题。将表 9表 7的计算结果进行对比,结果见表 10(单位mGal)。通过表 10可以看出,不同计算方法的段差较差基本保持在同一水平,但测点重力值较差随着作业的先后在逐渐增大,表明“重力测量电子记簿系统”所采用的后期纳入一次项系数计算会影响测点间的段差值及测点的重力值。究其原因,是测前段差包含了潮汐改正、气压改正和仪器高改正,也是“重力测量电子记簿系统”将潮汐改正、气压改正和仪器高改正加入一次项系数计算,而仪器一次项系数仅与仪器有关,将潮汐改正、气压改正和仪器高改正纳入一次项系数计算是不合适的。

表 10 段差(往返不同天)计算结果比较 Tab. 10 Comparison of segment difference calculations completed in different days
4 闭合环综合评价

《相对重力联测资料质量评比评分细则》[4]在环闭合差的评价方面有相关得分要求,而现有“重力测量电子记簿系统”尚不能完成闭合环的搜寻和环闭合差统计,需要增加上述功能并在野外数据成环测段进行环闭合差均值的计算。环闭合差均值的计算公式[4]为:

$ A_6=\sum\limits_{n=1}^N\left|B_n\right| \cdot \sqrt{C_n / 5} / N $ (7)

式中,N为闭合环个数,Bn为第n个闭合环的闭合差,Cn为第n个闭合环测段数(或节点数)。

假定有不同测段数的4个闭合环,每个环闭合差及相关计算结果见表 11(单位mGal)。

表 11 环闭合差均值计算比较 Tab. 11 Comparison of mean of loop closure error

表 11可以看出,本文最后得出的环闭合差均值为0.117 8 mGal,其值远大于文献[4]要求的最大值0.036 mGal,然而通过Bn/Cn的比值发现,每一个环内测段的平均闭合差相差并不大,其根本原因在于环内的测段数多,造成乘积的结果偏差大。

假设某个测区只有1个如表 11中的3号闭合环,11个测段的环闭合差0.030 8 mGal,按式(7)计算出的环闭合差均值为0.045 7 mGal;另外假设3号环只有6个测段,其环闭合差也是0.030 8 mGal,计算出的环闭合差均值为0.033 7 mGal。通过11个测段与6个测段的闭合环比较可以看出,在环闭合差相同的情况下,测段数多的闭合环精度较测段数少的闭合环精度差,这显然是不合适的。由此可见,文献[4]中的算法会造成结果值偏大,尤其对于测段数较多的闭合环来说其评价结果是非常不利的。因此,|Bn|后的乘号更改为除号更合适,这样表 11中4个环的环闭合差均值为0.024 5 mGal。

5 结语

1)“重力测量电子记簿系统”中工程文件到平差交换文件缺少仪器高数据,造成野外量测仪器高似乎是多余的,然而仪器架设在每一测点或同一测点往返测的高度并不相同,为将重力值归算到观测墩台面,仪器高改正不可缺少,因此工程文件到平差交换文件需要完善仪器高的输出。

2) 一次项系数仅与仪器有关,而与气压改正、仪器高改正、潮汐改正、零漂改正等无关,因此一次项系数参与计算应该在格值转换后完成,而目前所采用的“重力测量电子记簿系统”后期纳入一次项系数计算会影响测点间的段差值和测点的重力值成果。

3) 简单线路的各项改正计算较为明确和简单,但涉及测点间歇或线路上增设测点后漂移率的计算较为复杂。

4) 对文献[4]中成环的环平均闭合差计算公式存疑,建议将公式中|Bn|后的乘号修正为除号。

参考文献
[1]
中国地震局地震研究所. CG-5重力仪的地震重力测量应用手册[Z]. 2007 (Institute of Seismology, CEA. Application Manual of Seismic Gravity Measurement of CG-5 Gravimeter[Z]. 2007) (0)
[2]
国家地震局. 地震重力测量规范[S]. 北京: 地震出版社, 1979 (China Earthquake Administration. Specification for Seismic Gravity Survey[S]. Beijing: Seismological Press, 1979) (0)
[3]
Torge W. 重力测量学[M]. 北京: 地震出版社, 1993 (Torge W. Gravimetry[M]. Beijing: Seismological Press, 1993) (0)
[4]
中国地震局地震研究所. 相对重力联测资料质量评比评分细则[Z]. 2022 (Institute of Seismology, CEA. Detailed Rules for Quality Evaluation and Grading of Relative Gravity Joint Survey Data[Z]. 2022) (0)
Discussion on Data Processing of Relative Gravity Survey
YI Tianyang1     ZHENG Bing1     LI Yalun1     WEN Junjun1     YI Songquan1     ZHANG Ren1     
1. Crustal Deformation Monitoring Center of Sichuan Earthquake Agency, 139 Shangba Road,Ya’an 625000, China
Abstract: In accordance with the relevant theories and specifications of the Standard for Earthquake Gravity Survey, we discuss several problems urgently need to be solved in practical work, and existing in the currently used data processing software of the gravity survey electronic recording system, such as instrument height correction, section difference calculation, ring closure difference calculation, etc.
Key words: grid value table; barometric correction; height correction; tidal correction; zero drift correction