地球物理学进展  2015, Vol. 30 Issue (1): 421-424   PDF    
利用新参数和软件改进重力固体潮计算程序
董良, 彭芳苹, 杨涛, 白玉珍, 葛志广, 刘玲    
北京勘察技术工程有限公司, 北京 100085
摘要:重力固体潮的计算无论是在流动重力观测的校正还是在固定台站重力观测值的校正上都是必不可少的, 并且很大程度上决定了最终结果的准确性.固体潮的计算程序比较成熟, 但是目前大多数程序所采用的公式、参数和方法比较陈旧.为了提高计算结果的时效性和准确度, 经过对原公式的详细分析和计算程序的改进, 本文给出一个改进的计算公式和计算方法.特别在计算儒略日和六个天文参数时采用J2000.0系统, 更合时宜.另外, Matlab是非常优秀的数学计算软件, 本文利用Matlab的优势, 重新编写固体潮计算程序, 最后利用改进的计算程序给出一个计算实例, 证明了程序的可靠性和结果的精确性, 得到较好的结果.因此, 本文利用新参数和软件改进重力固体潮计算程序的方法是成功的, 值得在相关计算领域推广.
关键词固体潮     Matlab     儒略日     J2000.0系统    
To improve the earthtide calculating program by new parameters within software
DONG Liang, PENG Fang-ping, YANG Tao, BAI Yu-zhen, GE Zhi-guang, LIU Ling    
Beijing Explo-Tech Engineering Co., Ltd., Beijing 100085, China
Abstract: The earthtide calculating is essential for the observational data correction either in mobile gravitational survey or in stationary gravitational survey. The precision of the earthtide calculating determines the accuracy of the data correction. The programs of earthtide calculating that based on old formulas are well-understood but outdated. A new formula and calculating method are discussed in this paper to update the program and improve the precision of the calculating. J2000.0 system of equations is used in the calculation for Julian day and the specified six astronomical parameters. The earthtide calculating program is reprogramed in Matlab which is brilliant mathematics software. An instance based on the improved program is given at the end of this paper and shows a great result. Therefor the updated method of earthtide calculating with new parameters is effective and worth to be popularized.
Key words: Earthtide     Matlab     Julian day     J2000.0 system    
0 引 言

地球固体潮十分微弱,但是在现代高精度重力勘探和固定重力台站的重力场观测中,相应的固体潮改正却不容忽视(王谦身,2003).因此快速、准确的计算固体潮理论值是十分重要的.固体潮理论发展较早,目前已经比较成熟,但随着计算机技术的发展,原来传统的计算程序已经不合时宜.本文在传统的固体潮计算程序之上,利用优秀的数学计算软件Matlab,参照最新的参数和公式,对固体潮计算程序进行改进.为固体潮理论值的计算提供了更为可靠、便捷、迅速和实用的计算方法.

1 一般固体潮的计算步骤

固体潮的计算公式(中华人民共和国地质矿产部,1997)为

核心的计算步骤(霍志周,2006)大体可分为以下几步:
(1)根据测点所在的纬度φ,求出F(φ),

(2)求出计算时刻的儒略世纪数T,

其中,T为儒略世纪数,T0为儒略日.t为北京时间,tz为所在地时区.

儒略日T0的计算采用自定义函数:

(3)根据计算所得儒略世纪数T计算下列六个天文参数

(4)根据六个天文参数计算月亮 c r和cosZ

其中,

φ为所在地地心纬度,需要注意的是,φ和φ在计算过程中均为角度,计算之后需要转换成弧度.这一点缺乏明确的数学意义.

(5)根据六个天文参数计算太阳cosZs
(6)根据F(φ),月亮cosZ,太阳cosZs和相应时间和经纬度即可以得到固体潮和相应的日、月分潮.

以上具体计算公式和参数意义,在参考文献[3]中有详细阐述,在此不做过多赘述. 2 固体潮计算程序的改进

利用最新的参数,在原有的程序上进行了改进,使之更趋成熟和完善,并利用Matlab软件的优势,使计算程序更为可靠,输出结果更为直观和便捷.

(1)儒略日的计算采用Matlab自带的儒略日计算函数

函数(6)在Matlab中的定义为

利用⑥,避免了人为编写程序产生的错误,直接调用Matlab内置函数,也简化了程序编写过程,提高程序执行效率.笔者曾经接触过不少儒略日计算程序,很多都在计算2000年的时候出现问题.

(2)一般计算所用天文参数标准历元为J1900.0,改进后计算程序采用最新公布的J2000.0系统(郗钦文和候天航,1986),

相应的儒略世纪数计算公式改为

其中T自2000年1月1日12时(JD=2451545.0T DB)算起.

(3)同时,程序中的地心纬度计算采用公式(张少泉,1985):

其中e为地球扁率,等于,公式化为

公式⑤因为缺乏明确的数学意义,并且表现形式也不符合数学计算规则,给初次接触固体潮计算程序的研究人员造成很大的困扰.这应该是在计算机能力十分有限的年代,研究人员为了节约计算机机时而给出的近似公式.公式⑩比起⑤来说虽然较为复杂,但对于目前的计算机能力来说增加的计算量可以忽略,但公式⑩数学意义更明确,使得程序更为严谨和可靠,计算结果也更精确.

(4)最后,利用Matlab的绘图函数,增加了计算程序的绘图功能,使计算程序更直观和快捷的输出结果.而之前,要进行如下文(图 1)中所示图片的绘图工作,是十分困难和复杂的,利用Matlab可以很轻易的快捷的实现.

图 1 改进前后两套程序某地同一天计算结果ET为总固体潮,gs为太阳引起分潮,gm为月亮引起分潮Fig. 1 After the improvement of two sets of procedures in the same day results ET is the total solid tide,GS solar induced tidal component, GM as the moon caused by tidal.
3 结果对比

本文分别利用一般计算方法和改进后的程序编写方法,在Matlab中编写了两套固体潮计算程序,并在程序中增加数据绘图部分,使计算数据能更加直观的展现.对计算结果做了一个比较详细的对比,由此反映出程序改进前后的计算偏差.

图 2可以看出,改进前后,固体潮之差幅值不大,仅在-0.1~0.1微伽之间.但从图 1表 1可以看出,程序改进前后,计算结果总固体潮和太阳、月亮分量潮极值都存在相对明显的相位差,最小差0 s,最大达到11 s.这一偏差对于精确分析固体潮观测值的研究工作来说是值得注意的(Tang et al., 2013).

图 2 程序改进前后固体潮之差Fig. 2 Process improvement and solid tide difference

表 1 改进前后计算结果总固体潮和分量极值时间差(单位:s)Table 1 After the improvement of calculation results of the tide and component extremum time difference.(unit: s)
4 改进后的计算程序在固定重力观测校正中的应用

运用改进后的固体潮计算程序,选取某地重力台站固体潮观测数据(观测数据经过滤波、去零漂等预处理(文武等,2013))在Matlab下进行比对研究,结果如下:

图 3可以看出改进后的固体潮程序所得到的计算值和实际观测值的对应关系很好,计算结果是十分可靠的.

图 3 某地重力台站三天固体潮实测值和改进后计算值对比Fig. 3 A three day tidal gravity station value and improve the calculation value of contrast
5 结 论

本文改进后的固体潮计算程序得到的理论值与实际观测值对照显示结果十分可靠.改进后的程序编写更为严谨,计算更为精确,更具有时效性,并且充分利用了Matlab软件的优点,使得结果输出和展示更为直观和便捷.

致 谢    在本文写作过程中,得到了公司田黔宁、许德树、于国明和何贤明四位教授的指导和建议,以及文武博士无私的帮助和宝贵意见,一并表示衷心的感谢.
参考文献
[1] Chen X D, Sun H P. 2002. New method for pre-processing and analyzing tidal gravity observations[J]. Journal of Geodesy and Geodynamics (in Chinese), 22(3): 83-87.
[2] Fang J. 1984. Earth Tide (in Chinese)[M]. Beijing: Science Press.
[3] He Z Z. 2006. Development of the gravity data processing system (in Chinese)[Thesis]. Xi'an: Department of Geophysics, Chang'an University.
[4] Li G Y, Peng L H, Xu H Z. 1996. Tidal deformation of a Self-rotating, Elliptical, Heterogeneous earth[J]. Chinese J. Geophys. (in Chinese), 39(5): 672 -678.
[5] Li H, Li R H. 1994. Filter preprocessing method for the data of Earth tide[J]. Crustal Deformation and Earthquake (in Chinese), 14(1): 23-31.
[6] Ministry of Geology and Mineral Resources of the People's Republic of China. 1997. Gravity specification of Geology and mineral resources industry standards of the people's Republic of China large scale exploration DZ/T0171-1997 (in Chinese)[S].
[7] Sun H P, Hsu H, Chen W, et al. 2006. Study of Earth's gravity tide and oceanic loading characteristics in Hong Kong area[J]. Chinese J. Geophys. (in Chinese), 49(3): 724-734, doi: 10.3321/j.issn:0001-5733.2006.03.016.
[8] Tang K Y, Hua C C, Wen W, et al. 2013. Observational evidences for the speed of the gravity based on the Earth tide[J]. Chin. Sci. Bull., 58(4-5): 474-477.
[9] Wang Q S. 2003. Gravitology (in Chinese)[M]. Beijing: Seismological Press.
[10] Wen W, Tang K Y, Wang Q S, et al. 2013. To discriminate the gravity anomaly by the subtle gravity observation of the total solar eclipse in 2009[J]. Chinese J. Geophys. (in Chinese), 56(3): 770-782, doi: 10.6038/cjg20130306.
[11] Xi Q W. 1986. Earth tides and generating tide constants[J]. Earthquake Research in China (in Chinese), 2(2): 30-41.
[12] Xu H Z. 2010. Solid Earth Tides (in Chinese)[M]. Wuhan: Science and Technology Press, 111-113.
[13] Xu J Q. 1997. The Theory and Analysis of Gravity Tidal——Superconducting Gravimeter Data Analysis and Processing in Wuchang Station (in Chinese)[Z]. Wuhan: Institute of Geodesy and Geophysics. Chinese Academy of Sciences.
[14] Zhang J, Xi Q W, Yang L Z, et al. 2007. A study on tidal force/stress triggering of strong earthquakes[J]. Chinese J. Geophys. (in Chinese), 50(2): 448-454.
[15] Zhang S Q. 1985. On the conversion of geographical latitude and geocentric latitude[J]. Northwestern Seismological Journal (in Chinese), 7(1): 39-43.
[16] 陈晓东, 孙和平. 2002. 一种新的重力潮汐数据预处理和分析方法[J]. 大地测量与地球动力学, 22(3): 83-87.
[17] 方俊. 1984. 固体潮[M]. 北京: 科学出版社.
[18] 霍志周. 2006. 重力资料预处理系统研制[硕士论文]. 西安: 长安大学地球物理系.
[19] 李国营, 彭龙辉, 许厚泽. 1996. 自转微椭、非均匀地球的潮汐变形[J]. 地球物理学报, 39(5): 672 -678.
[20] 李辉, 李瑞浩. 1994. 固体潮观测数据的滤波预处理方法[J]. 地壳形变与地震, 14(1): 23-31.
[21] 孙和平, 许厚泽, 陈武,等. 2006. 香港地区重力固体潮和海潮负荷特征研究[J]. 地球物理学报, 49(3): 724-734, doi: 10.3321/j.issn:0001-5733.2006.03.016.
[22] 王谦身. 2003. 重力学[M]. 北京: 地震出版社.
[23] 文武, 汤克云, 王谦身,等. 2013. 利用2009年日全食的精细重力观测探寻“引力异常”[J]. 地球物理学报, 56(3): 770-782, doi: 10.6038/cjg20130306.
[24] 郗钦文, 候天航. 1986. 固体潮汐与引潮常数[J]. 中国地震, 2(2): 30-41.
[25] 许厚泽. 2010. 固体地球潮汐[M]. 武汉: 湖北科学技术出版社, 111-113.
[26] 徐建桥. 1997. 重力固体潮汐理论及分析方法——武昌台超导重力仪观测资料的分析处理[Z]. 武汉: 中国科学院测量与地球物理研究所.
[27] 张晶, 郗钦文, 杨林章,等. 2007. 引潮力与潮汐应力对强震触发的研究[J]. 地球物理学报, 50(2): 448-454.
[28] 张少泉. 1985. 有关地理纬度和地心纬度的换算问题[J]. 西北地震学报, 7(1): 39-43.
[29] 中华人民共和国地质矿产部. 1997. 中华人民共和国地质矿产行业标准-大比例尺重力勘探规范DZ/T0171-1997 [S].