朱傳東 陳 銘 明曉冉 龐柳青 田興宇 蘇建文

1 中國(guó)地震局第一監(jiān)測(cè)中心，天津市耐火路7號(hào)，300180

潮汐是對(duì)日月等天體引潮力的響應(yīng)，包括固體潮、海潮和大氣潮，潮汐改正是重力觀測(cè)數(shù)據(jù)處理的基礎(chǔ)。重力固體潮與海洋潮汐負(fù)荷是影響重力觀測(cè)的主要干擾因素，二者最大變化量分別可達(dá)300 μGal與10 μGal。通常使用引潮位對(duì)重力固體潮進(jìn)行分析計(jì)算，主要計(jì)算方法有天頂距直接法和分波法[1-2]。由于潮汐效應(yīng)具有相同的頻率，常規(guī)的濾波方式無(wú)法分離海洋潮汐負(fù)荷，因此需要利用高精度海潮模型進(jìn)行改正[3-6]。

本文基于MATLAB GUI界面自主編寫(xiě)一套高精度重力測(cè)量的潮汐改正軟件GTIDE，并結(jié)合國(guó)際主流軟件對(duì)其可靠性與適用性進(jìn)行分析，以期為重力觀測(cè)的數(shù)據(jù)處理與分析提供靈活可靠的軟件平臺(tái)。

1 數(shù)據(jù)處理方法

1.1 重力固體潮理論值

根據(jù)潮汐力學(xué)理論可知，重力固體潮理論值是將月亮與太陽(yáng)的引潮位分別取至三階與二階，再由天頂距直接計(jì)算求得(向上為正)[7]。GTIDE軟件計(jì)算給出的重力固體潮理論值為扣除永久潮汐引力直接部分后的零潮汐改正值(單位μGal)，軟件可以輸入由統(tǒng)計(jì)經(jīng)驗(yàn)?zāi)Ｐ瞳@得的平均潮汐因子(默認(rèn)1.16)或者由固體潮觀測(cè)臺(tái)站數(shù)據(jù)計(jì)算得到的實(shí)測(cè)潮汐因子[7-8]。月亮與太陽(yáng)的杜德森常數(shù)以及平均地球半徑都與采用的參考橢球和天文參數(shù)有關(guān)，GTIDE軟件數(shù)據(jù)取自1967年國(guó)際大地測(cè)量協(xié)會(huì)[7]。其他的主要參數(shù)與月球、太陽(yáng)的6個(gè)軌道運(yùn)動(dòng)天文參數(shù)有關(guān)，在J2000.0系統(tǒng)下可表示為儒略世紀(jì)數(shù)(以力學(xué)時(shí)為參改)的線性函數(shù)。此外，GTIDE軟件可以手動(dòng)輸入或者批量導(dǎo)入計(jì)算過(guò)程中需要用到的測(cè)站地理經(jīng)緯度坐標(biāo)以及觀測(cè)時(shí)間等數(shù)據(jù)。

1.2 海洋潮汐負(fù)荷

海洋潮汐負(fù)荷改正可以由海潮模型與重力格林函數(shù)的褶積積分求得(向上為正)[9]：

L(φ,λ,t)=

(1)

式中，φ和λ分別為重力測(cè)站的經(jīng)度和緯度，φ′和λ′分別為負(fù)荷點(diǎn)的經(jīng)度和緯度，ρw為海水密度，r為地球半徑，H(φ′,λ′)為全球或近海區(qū)域的瞬時(shí)潮高，G(ψ)是負(fù)荷點(diǎn)到測(cè)站角距為ψ的重力負(fù)荷格林函數(shù)。

從式(1)可以看出，海潮模型與重力格林函數(shù)是影響海潮負(fù)荷計(jì)算精度的主要因素。GTIDE軟件可以選擇性輸入NAO99b、TPXO、FES、EOT、DTU和CSR等一系列高精度全球海潮模型或NAO99Jb近海海潮模型。選中近海海潮模型后，軟件可將其與全球海潮模型相應(yīng)近海區(qū)域值進(jìn)行替換。GTIDE軟件統(tǒng)一將瞬時(shí)潮高表示為8個(gè)主要調(diào)和諧波分量(M2、S2、N2、K2、K1、O1、P1和Q1)的總和，同時(shí)為保證每個(gè)分潮波均滿足海水質(zhì)量守恒定律，需分別扣除等厚度的水質(zhì)量層以修正全球海潮模型[10]。對(duì)于重力格林函數(shù)，GTIDE軟件可以選擇性輸入特定坐標(biāo)系下的G-B、1066A、PREM、PREMsoft、PREMhard、iasp91和ak135[11]。在實(shí)際計(jì)算過(guò)程中，式(1)中的離散褶積積分區(qū)域共劃分為4個(gè)部分，GTIDE軟件可手動(dòng)輸入每個(gè)積分區(qū)域的邊界范圍及其離散網(wǎng)格大小，也可以在近區(qū)范圍內(nèi)選擇性輸入SRTM DEM數(shù)據(jù)以改善海陸邊界的識(shí)別精度。為改善重力格林函數(shù)直接引力項(xiàng)的計(jì)算精度，采用海水質(zhì)量引力效應(yīng)對(duì)測(cè)站垂向分量進(jìn)行計(jì)算[12]。

基于上述設(shè)置可計(jì)算出各分潮波的負(fù)荷振幅Lp與負(fù)荷相位αp，通過(guò)疊加計(jì)算可獲得海潮負(fù)荷重力效應(yīng)隨時(shí)間變化的特征：

Gp(φ,λ,t)=

(2)

式中，fp與up分別為交點(diǎn)因子和交點(diǎn)訂正角；t為觀測(cè)時(shí)間，可在軟件界面中手動(dòng)輸入或者批量導(dǎo)入，由此得到任意時(shí)刻的負(fù)荷值。各分潮波的負(fù)荷振幅、負(fù)荷相位以及負(fù)荷值時(shí)間序列可在選定的輸出文件中保存。

2 結(jié)果與分析

2.1 重力固體潮理論值

為驗(yàn)證GTIDE軟件計(jì)算得到的理論重力固體潮精度，將其與ETERNA軟件包中PREDICT程序計(jì)算得到的結(jié)果進(jìn)行比較。PREDICT程序采用分波法計(jì)算得到理論重力固體潮，在實(shí)際計(jì)算過(guò)程中可提供7種引潮位展開(kāi)表。本文選用Tamura的引潮位展開(kāi)表，其計(jì)算精度可達(dá)μGal級(jí)[13]。

以分布于不同緯度的4個(gè)重力測(cè)站為例(永興島站、廈門站、佘山站、薊縣站)，對(duì)上述2種軟件計(jì)算得到的理論重力固體潮結(jié)果進(jìn)行比較分析。計(jì)算過(guò)程中，將采樣時(shí)間跨度設(shè)置為2021年全年，采樣時(shí)間間隔設(shè)置為1 min，采樣參考時(shí)間選取世界時(shí)(UTC)，平均潮汐因子設(shè)置為1.16，同時(shí)保留永久潮汐引力直接部分，并從時(shí)間序列中扣除該時(shí)間段內(nèi)的平均值。表1為2種軟件得到的理論重力固體潮結(jié)果，圖1為8月薊縣站2種軟件的理論重力固體潮及差值。從全年采樣數(shù)據(jù)來(lái)看，理論重力固體潮在各個(gè)測(cè)站上的峰對(duì)峰值最大可達(dá)326.9 μGal，說(shuō)明重力固體潮會(huì)對(duì)重力觀測(cè)造成較大影響。2種軟件在各個(gè)測(cè)站上的結(jié)果差值均在0.7 μGal以內(nèi)，相關(guān)系數(shù)均為1.0，呈現(xiàn)出較為一致的變化特征，說(shuō)明本文GTIDE軟件的計(jì)算結(jié)果具有可靠性。2種軟件計(jì)算結(jié)果的差異可能與選用的參考橢球、天文參數(shù)、引潮位截?cái)嘁约耙蔽徽归_(kāi)表有關(guān)[14]。

表1 GTIDE與PREDICT理論重力固體潮計(jì)算結(jié)果

圖1 GTIDE與PREDICT的理論重力固體潮及其差值Fig.1 The theoretical gravity solid tide and its difference between GTIDE and PREDICT

2.2 海洋潮汐負(fù)荷

為驗(yàn)證GTIDE軟件計(jì)算得到的海洋潮汐重力負(fù)荷的精度，將其與Bos-Scherneck網(wǎng)站計(jì)算得到的結(jié)果進(jìn)行比較分析。Bos-Scherneck網(wǎng)站的海洋潮汐重力負(fù)荷值同樣由海潮模型與重力格林函數(shù)的褶積積分計(jì)算得到，在實(shí)際計(jì)算過(guò)程中可提供G-B、STW105重力格林函數(shù)與多種全球海潮模型的褶積計(jì)算結(jié)果。本文仍以上述4個(gè)重力測(cè)站為例(永興島站、廈門站、佘山站、薊縣站距海岸線分別約1 km、25 km、70 km和100 km)，對(duì)2種軟件計(jì)算得到的海洋潮汐重力負(fù)荷進(jìn)行比較分析。實(shí)際計(jì)算過(guò)程中，默認(rèn)重力測(cè)站的高程為0 m，將重力格林函數(shù)設(shè)置為G-B模型，將海潮模型設(shè)置為FES2004模型。表2為2種軟件得到的海洋潮汐重力負(fù)荷結(jié)果。可以看出，各分潮波的重力負(fù)荷最大可達(dá)3.5 μGal，說(shuō)明海洋潮汐會(huì)對(duì)重力觀測(cè)造成一定的影響，且這種影響會(huì)隨臺(tái)站與海洋之間距離的增大而減小。從各分潮波振幅結(jié)果來(lái)看，2種軟件在各個(gè)測(cè)站上的差值均在0.07 μGal以內(nèi)，除K2與N2分潮波外，其余占比均在自身量級(jí)的10%以內(nèi)；從各分潮波的相位結(jié)果來(lái)看，除K2與N2分潮波外，其余差值基本在6.5°以內(nèi)。2種軟件計(jì)算結(jié)果的差異可能與海岸線數(shù)據(jù)以及褶積積分方法有關(guān)[9]。值得注意的是，K2與N2分潮波的相位差值整體較大，可能是因?yàn)槎叩恼穹考?jí)相對(duì)較小，對(duì)上述影響更為敏感，但其對(duì)重力負(fù)荷計(jì)算結(jié)果的影響非常有限。

表2 GTIDE與Bos-Scherneck海洋潮汐重力負(fù)荷計(jì)算結(jié)果

圖2為2021-08四個(gè)測(cè)站2種軟件的海洋潮汐重力負(fù)荷時(shí)間序列，其中采樣時(shí)間間隔設(shè)置為1 min，采樣參考時(shí)間選取世界時(shí)(UTC)。2種軟件在永興島站、廈門站、佘山站和薊縣站的序列差值分別在0.25 μGal、0.19 μGal、0.29 μGal和0.14 μGal以內(nèi)，相關(guān)系數(shù)均為1.0左右，呈現(xiàn)出較為一致的變化特征，說(shuō)明GTIDE軟件的計(jì)算結(jié)果具有可靠性。

圖2 GTIDE與Bos-Scherneck的海洋潮汐重力負(fù)荷時(shí)間序列Fig.2 The time series of ocean tide loading of GTIDE and Bos-Scherneck

基于GTIDE軟件，本文進(jìn)一步分析海洋潮汐模型與重力格林函數(shù)對(duì)海洋潮汐重力負(fù)荷的影響。限于篇幅，本文僅選擇佘山站和廈門站的3種海洋潮汐模型與重力格林函數(shù)模型進(jìn)行比較分析。表3為FES2004、FES2014b和TPXO9_atlas海洋潮汐模型的重力負(fù)荷計(jì)算結(jié)果，重力格林函數(shù)均設(shè)置為G-B模型。其中，F(xiàn)ES2004、FES2014b與TPXO9_atlas模型的空間分辨率分別為1/8°、1/16°和1/30°。從表3可以看出，不同海洋潮汐模型的重力負(fù)荷計(jì)算結(jié)果差異較為明顯。相比于FES2004，F(xiàn)ES2014b與TPXO9_atlas在2個(gè)測(cè)站上的最大振幅差值可達(dá)0.80 μGal，最大相位差值可達(dá)49.4°(對(duì)應(yīng)振幅為0.46 μGal)，上述差異在沿海且海岸線較為復(fù)雜的廈門站更加明顯。綜上可知，海洋潮汐重力負(fù)荷計(jì)算結(jié)果對(duì)海潮模型的精度較為敏感，因此需要采用更高精度的區(qū)域海潮模型、驗(yàn)潮站數(shù)據(jù)以及海岸線數(shù)據(jù)進(jìn)行模型構(gòu)建。

表3 不同海潮模型的海洋潮汐重力負(fù)荷計(jì)算結(jié)果

表4對(duì)不同格林函數(shù)的重力負(fù)荷計(jì)算結(jié)果進(jìn)行比較分析，其中格林函數(shù)的地球模型分別為G-B、PREM和PREMhard，海潮模型均設(shè)置為FES2004。3種地球模型的區(qū)別主要體現(xiàn)在地殼上地幔區(qū)域，不同于G-B、PREM這2種全球平均地球模型，PREMhard是以Crust2.0固結(jié)沉積層為最外層的地球模型，能夠有效改善PREM縱向地殼彈性結(jié)構(gòu)參數(shù)的精度。由表可見(jiàn)，相比于G-B 模型，PREM與PREMhard在2個(gè)測(cè)站上的最大振幅差值可達(dá)0.14 μGal，最大相位差值可達(dá)12.6°(對(duì)應(yīng)振幅為0.25 μGal)。由此可知，雖然不同地球模型得到的重力格林函數(shù)差異明顯(地殼上地幔參數(shù)的改進(jìn)顯著改善了高階負(fù)荷勒夫數(shù)和測(cè)站近區(qū)的格林函數(shù))，但通過(guò)褶積積分計(jì)算可知，其對(duì)海洋潮汐重力負(fù)荷計(jì)算結(jié)果的影響有限。

表4 不同格林函數(shù)的海洋潮汐重力負(fù)荷計(jì)算結(jié)果

3 結(jié) 語(yǔ)

本文基于MATLAB GUI界面，采用天頂距直接法與負(fù)荷格林函數(shù)法自主編寫(xiě)一套高精度重力測(cè)量潮汐改正軟件GTIDE，并將其與國(guó)際主流軟件進(jìn)行比較分析，驗(yàn)證了該軟件的可靠性與適用性。GTIDE軟件易于操作、便于維護(hù)，可為高精度重力觀測(cè)數(shù)據(jù)的處理與分析提供較為可靠的軟件平臺(tái)。

本文分析結(jié)果表明，天頂距直接法得到的結(jié)果可以滿足高精度流動(dòng)重力觀測(cè)精度的要求，而一些沿海地區(qū)以及海島的海潮負(fù)荷影響較大。另外，海潮重力負(fù)荷的改正結(jié)果對(duì)海潮模型的精度較為敏感，需要采用更高精度的區(qū)域海潮模型、驗(yàn)潮站資料以及海岸線數(shù)據(jù)進(jìn)行模型構(gòu)建，以提高重力觀測(cè)中的潮汐改正精度。