李盛 胡久常 周雯 張慧 王慧琳 謝小玲

海南省地震局，?？?570203

0 引言

重力隨時(shí)間的變化可分為潮汐變化和非潮汐變化，前者如地球潮汐與地球自轉(zhuǎn)變化，后者一般是地球物質(zhì)遷移引起的重力隨時(shí)間的變化(王謙身，2003)。連續(xù)重力觀測(cè)是在臺(tái)站固定點(diǎn)上進(jìn)行的連續(xù)觀測(cè)，獲得的臺(tái)站重力變化是臺(tái)站高程以及臺(tái)站和其周圍地殼物性特征隨時(shí)間變化的結(jié)果。連續(xù)重力觀測(cè)能夠精確探測(cè)到地球系統(tǒng)各圈層物質(zhì)遷移引起的重力變化效應(yīng)，包括疊加在一起的長(zhǎng)期(多年的)或短期(季節(jié)性的)的變化(Sun et al，2001、2009、2011； 周江林等，2015； 翟麗娜等，2020)。連續(xù)重力觀測(cè)獲得的數(shù)據(jù)由重力潮汐變化、重力非潮汐變化和誤差構(gòu)成，即

g=gtide+gnon-tide+e(t)

(1)

重力潮汐變化主要是日月的引潮力，而重力非潮汐變化主要由觀測(cè)儀器本身的零漂、自然環(huán)境干擾、人為因素和殘差所構(gòu)成，即

gnon-tide=glinear-drift+genvironment+ghuman+e(t)

(2)

自然環(huán)境干擾引起的重力變化主要是由氣壓負(fù)荷、地下水、極移、海潮等引起，即

genvironment=gair-pressure+gpolar-shift+ggroundwater+gocean-tide+e(t)

(3)

為獲得重力非潮汐變化，需對(duì)重力潮汐變化數(shù)據(jù)進(jìn)行零漂、固體潮、氣壓和海潮等的改正，本文所要研究的就是分別定量分析氣壓和海潮對(duì)瓊中臺(tái)重力非潮汐變化的影響，獲得氣壓和海潮負(fù)荷改正后的瓊中臺(tái)重力變化。海南島位于中國大陸南端，四面環(huán)海，瓊中地震臺(tái)(109.8°E，19.0°N，屬于熱帶地區(qū)，地理位置見圖1)位于海南島中部，2008年3月開始連續(xù)重力觀測(cè)，運(yùn)行儀器為PET重力儀(美國 Micro-g LaCoste，Inc公司生產(chǎn))，該重力儀是一種全自動(dòng)的金屬彈簧(零長(zhǎng)彈簧)連續(xù)重力觀測(cè)儀，是在LCR-G型基礎(chǔ)上改進(jìn)而成，工作原理與LCR-G相同。瓊中臺(tái)距離最近的海岸線約60km，臺(tái)站受熱帶風(fēng)暴、臺(tái)風(fēng)和強(qiáng)降雨氣候影響頻繁，屬于熱帶季風(fēng)氣候，臺(tái)風(fēng)季節(jié)集中在夏末秋初，8至10月是臺(tái)風(fēng)多發(fā)期，氣壓變化具有典型的熱帶特征(李盛等，2016； 李盛等，2020)。基于瓊中臺(tái)這一鮮明的地理環(huán)境和氣候特征，對(duì)瓊中臺(tái)重力非潮汐變化進(jìn)行氣壓與海潮負(fù)荷改正很有必要，可準(zhǔn)確認(rèn)識(shí)瓊中臺(tái)重力變化，充分發(fā)揮瓊中臺(tái)重力基準(zhǔn)作用。

圖1 瓊中地震臺(tái)位置

1 原始資料分析

對(duì)原始數(shù)據(jù)進(jìn)行分析，可掌握儀器的運(yùn)行情況，為后續(xù)的非潮汐分析打下基礎(chǔ)。對(duì)瓊中臺(tái)2009年1月至2019年9月連續(xù)重力原始數(shù)據(jù)進(jìn)行臺(tái)階、突跳等處理，獲得該段時(shí)間瓊中臺(tái)重力潮汐變化數(shù)據(jù)(圖2(a))，由圖可知瓊中臺(tái)重力潮汐變化總體比較穩(wěn)定，固體潮清晰，但漂移存在非線性的情況(數(shù)據(jù)采樣率為分鐘值)。

圖2 瓊中臺(tái)預(yù)處理后的重力潮汐變化數(shù)據(jù)(a)、擬合零漂(b)與零漂改正后殘差(c)時(shí)序圖

2 潮汐分析

1966年保加利亞的Venedikov教授提出了利用48階奇偶濾波器采用濾波方式分別計(jì)算周日、半日和 1/3 日波的潮波分量，并利用上述分量進(jìn)行潮汐參數(shù)的估計(jì)(Venedikov et al，2003)。目前該方法被國家重力臺(tái)網(wǎng)中心(1)http: //www.gncc.ac.cn在結(jié)合Nakai檢驗(yàn)方法和數(shù)字濾波方法(李輝等，1994)基礎(chǔ)上用于進(jìn)行連續(xù)重力觀測(cè)資料的評(píng)價(jià)工作。2003年，Venedikov進(jìn)一步改進(jìn)了程序，引入日本京都大學(xué)Tamura 等(1991) 1200個(gè)分波的潮波表，將氣壓等輔助觀測(cè)序列、分段多項(xiàng)式等方法引入潮汐分析數(shù)學(xué)模型中，使得程序能夠?qū)θ我獬毕^測(cè)數(shù)據(jù)進(jìn)行分析，即為VAV潮汐分析軟件。

運(yùn)用VAV潮汐分析軟件，對(duì)2009年1月至2019年9月瓊中臺(tái)連續(xù)重力觀測(cè)數(shù)據(jù)進(jìn)行潮汐分析，得到瓊中臺(tái)重力潮汐變化觀測(cè)數(shù)據(jù)的各主要潮波參數(shù)(表1)。其中，M2波(月亮的主半日波)的振幅最大，S2波(太陽的主半日波)次之。實(shí)際上，在各潮波中，M2波振幅最大，也最為穩(wěn)定，因此在日常分析中主要跟蹤M2波潮汐因子的變化。

表1 瓊中臺(tái)連續(xù)重力各潮汐波參數(shù)

3 重力非潮汐變化改正

3.1 零漂改正

彈簧重力儀的零漂與彈性材料以及儀器的設(shè)計(jì)方式有著直接的關(guān)系。如果考慮彈性系統(tǒng)的彈性滯后、蠕變變形、彈性后效、彈簧年齡、內(nèi)部傳感器溫度和彈簧的應(yīng)力加載等因素影響，可以模擬彈性系統(tǒng)的零漂。但由于上述因素的差異性以及建模的困難，本文選擇一般多項(xiàng)式擬合方法模擬彈黃重力儀的零漂。

根據(jù)瓊中臺(tái)連續(xù)重力觀測(cè)數(shù)據(jù)的處理情況，將2009年1月至2019年9月重力潮汐變化數(shù)據(jù)分段，然后對(duì)分段數(shù)據(jù)的零漂進(jìn)行一般多項(xiàng)式分段曲線擬合(圖2(b))，即可獲得瓊中臺(tái)重力潮汐變化數(shù)據(jù)的零漂擬合值以及零漂改正值(圖2(c))。

3.2 固體潮改正

采用國際地潮中心推薦的Tsoft固體潮數(shù)據(jù)預(yù)處理軟件作為觀測(cè)數(shù)據(jù)預(yù)處理的分析軟件(van Camp et al，2005)。該軟件可計(jì)算基于WDD模型(非流體靜力平衡的地球模型)的理論固體潮值。運(yùn)用該軟件，計(jì)算可得瓊中臺(tái)理論固體潮，對(duì)零漂改正后瓊中臺(tái)重力殘差(圖2(b))進(jìn)行理論固體潮改正，得到零漂固體潮改正后的重力殘差曲線(圖3)。

圖3 瓊中臺(tái)零漂固體潮改正后的重力殘差時(shí)序曲線

3.3 氣壓改正

表2 不同機(jī)構(gòu)的大氣重力導(dǎo)納值

為了對(duì)瓊中臺(tái)重力非潮汐變化進(jìn)行氣壓改正，運(yùn)用VAV軟件，計(jì)算了瓊中臺(tái)2009年1月至2019年9月每年的大氣重力導(dǎo)納值，計(jì)算結(jié)果包含了各頻段(周期1～7cpd)的大氣重力導(dǎo)納值(表3)。由表3 可知，2009年1月以來，瓊中臺(tái)每年各頻段的氣壓導(dǎo)納值介于(-0.15～-0.49)×10-8m/s2/mbar之間，11年中有5年的氣壓導(dǎo)納值位于(-0.30～-0.40)×10-8m/s2/mbar之間。而若整體計(jì)算2009年1月至2019年9月的氣壓導(dǎo)納值，則該值為-0.34×10-8m/s2/mbar，與表2 中的結(jié)果接近。說明短期來看，瓊中臺(tái)的大氣重力導(dǎo)納值有所差異，但長(zhǎng)期而言，瓊中臺(tái)的大氣重力導(dǎo)納值與其他學(xué)者研究不同區(qū)域的結(jié)果較為一致。

表3 瓊中臺(tái)2009年1月—2019年9月每年各頻段大氣重力導(dǎo)納值(單位：×10-8m/s2/mbar)

基于上述結(jié)果，取2009年1月至2019年9月的大氣重力導(dǎo)納值，即-0.34×10-8m/s2/mbar，對(duì)每年的重力非潮汐變化進(jìn)行改正，最終獲得了氣壓改正后的瓊中臺(tái)重力非潮汐變化(以下簡(jiǎn)稱氣壓改正)。

3.4 海潮改正

海洋潮汐是海水在太陽和月亮的引潮力作用下周期性漲落的自然現(xiàn)象。海洋潮汐是重力固體潮觀測(cè)的重要干擾源。研究表明，在沿海地區(qū)重力固體潮中的海潮負(fù)荷一般可達(dá)幾十微伽量級(jí)。為有效利用重力資料研究地球物理學(xué)和地球動(dòng)力學(xué)問題，對(duì)重力觀測(cè)資料的海潮負(fù)荷信號(hào)進(jìn)行改正就顯得特別重要(孫和平等，2002)。Farrell(1972)給出了計(jì)算海潮負(fù)荷的負(fù)荷格林函數(shù)，并指出海潮負(fù)荷即海潮潮高和負(fù)荷格林函數(shù)的褶積積分

L(θ，λ，t)=R2?ρH(φ，A，t)G(φ)sinφdφdA

(4)

式中，R為地球半徑，ρ為海水密度，θ為計(jì)算點(diǎn)的余緯，λ為計(jì)算點(diǎn)的經(jīng)度，(φ，A)為模板坐標(biāo)，t為時(shí)間，G(φ)為重力負(fù)荷格林函數(shù)。通常，負(fù)荷結(jié)果均以余弦波的振幅Amp和相位Φ給出，即任一時(shí)刻的負(fù)荷值為

L(θ，λ，t)=Amp·cos(Φ+φ)

(5)

因此，總結(jié)海潮負(fù)荷包括三方面內(nèi)容：①海潮模型； ②格林函數(shù)； ③計(jì)算方法。

由于存在多種海潮模型，因此需要對(duì)其進(jìn)行選擇。海潮模型對(duì)重力負(fù)荷計(jì)算的影響具有明顯的地域性特點(diǎn)(孫和平等，2005)，全球區(qū)域沒有絕對(duì)優(yōu)越的海潮模型。研究表明，近海潮汐對(duì)沿海臺(tái)站的影響非常大(Sun，1992； 周江存等，2004)。目前應(yīng)用較為廣泛的有CSR4.0、NAO.99b、FES2004、GOT04、TPXO7.2、DTU10、EOT11a、HAMTIDE11a等8個(gè)海潮模型(Kim et al，2011、2013)。

針對(duì)眾多海潮模型，Yu(2006)認(rèn)為 NAO.99b為最好的潮汐模型，該模型是根據(jù)日本天文臺(tái)應(yīng)用T/P衛(wèi)星數(shù)據(jù)，結(jié)合日本驗(yàn)潮站資料構(gòu)建的區(qū)域模型，精度為0.25°×0.25°。因此，本文基于SPOTL(Agnew，1997)程序，以NAO.99b潮汐模型計(jì)算了瓊中臺(tái)海潮負(fù)荷值，并運(yùn)用該負(fù)荷值對(duì)瓊中臺(tái)重力非潮汐零漂、固體潮、氣壓改正后的數(shù)據(jù)進(jìn)行海潮改正(以下簡(jiǎn)稱海潮改正)。將零漂固體潮改正與海潮改正后的瓊中臺(tái)重力非潮汐變化數(shù)據(jù)進(jìn)行對(duì)比(圖4)，認(rèn)為海潮改正后的重力殘差值變化幅度約為35×10-8m/s2，大于僅進(jìn)行潮汐、零漂改正后的重力殘差值變化幅度。

圖4 瓊中臺(tái)重力潮汐變化零漂固體潮氣壓海潮改正時(shí)序圖

運(yùn)行VAV軟件，分別對(duì)零漂固體潮改正、氣壓改正、海潮改正后的瓊中臺(tái)重力非潮汐變化數(shù)據(jù)進(jìn)行Venedikov調(diào)和分析，獲得以上數(shù)據(jù)的各主要潮汐波參數(shù)， 并將其與所有改正前的瓊中臺(tái)重力潮汐變化數(shù)據(jù)的主要波群參數(shù)進(jìn)行比較，結(jié)果表明經(jīng)海潮改正后的各潮波振幅明顯小于僅進(jìn)行零漂固體潮改正的振幅，其中變化幅度最大的為M2波，其振幅由1.70438×10-8m/s2降至0.44651×10-8m/s2。因此，進(jìn)行海潮負(fù)荷改正可進(jìn)一步消除重力潮汐變化數(shù)據(jù)中的潮汐信號(hào)，提高提取重力非潮汐變化的精度。

表4 瓊中臺(tái)重力潮汐變化、零漂固體潮改正、氣壓改正與海潮改正數(shù)據(jù)各主要潮波振幅(單位：×10-8m/s2)

3.5 小結(jié)

根據(jù)以上計(jì)算的結(jié)果，將瓊中臺(tái)重力非潮汐變化氣壓與海潮負(fù)荷改正值繪制在同一圖中(圖5)，可知?dú)鈮焊恼导s為10×10-8m/s2，海潮改正值約為5×10-8m/s2。

圖5 瓊中臺(tái)重力氣壓負(fù)荷和海潮負(fù)荷改正值

從2012年以來瓊中臺(tái)海潮改正后的重力非潮汐年變化來看(圖6)，2015年以來年變規(guī)律基本呈春冬低、夏秋高的變化趨勢(shì)，其中6—7月為一年中重力最高值，12—1月為重力最低值，2012年以來瓊中臺(tái)海潮改正后重力非潮汐變化未出現(xiàn)趨勢(shì)變化的異常情況。

圖6 瓊中臺(tái)海潮改正后的重力非潮汐年變序列

4 結(jié)論

通過對(duì)2009年1月至2019年9月瓊中臺(tái)重力潮汐變化數(shù)據(jù)的潮汐分析以及非潮汐分析，得到如下認(rèn)識(shí)：

(1)瓊中臺(tái)的大氣重力導(dǎo)納值總體上與不同學(xué)者對(duì)不同區(qū)域研究所得的結(jié)果較為一致，未來對(duì)瓊中臺(tái)重力非潮汐變化數(shù)據(jù)進(jìn)行氣壓改正時(shí)，可直接采用 -0.34×10-8m/s2/mbar作為大氣重力導(dǎo)納值。

(2)楊錦玲等(2016)運(yùn)用8個(gè)海潮模型計(jì)算廈門臺(tái)海潮負(fù)荷，經(jīng)對(duì)比，顯示經(jīng)8個(gè)全球海潮模型改正后，主要潮波的海潮負(fù)荷振幅差異變小。因此，對(duì)瓊中臺(tái)重力非潮汐變化的海潮改正選取NAO.99b潮汐模型是合適的。計(jì)算結(jié)果表明瓊中臺(tái)受到的海潮負(fù)荷變幅約5×10-8m/s2， 小于廈門的海潮負(fù)荷變幅7.3×10-8m/s2，這應(yīng)與瓊中臺(tái)距離海洋稍遠(yuǎn)有關(guān)(廈門臺(tái)距離海岸線不到1km)。

(3)瓊中臺(tái)重力非潮汐變化氣壓改正變幅約為10×10-8m/s2，海潮改正變幅約為 5×10-8m/s2，氣壓改正變幅大于海潮改正變幅。經(jīng)海潮改正后的瓊中臺(tái)重力非潮汐變化數(shù)據(jù)中的潮汐信號(hào)比僅進(jìn)行零漂固體潮改正更加微弱，說明進(jìn)行海潮改正具有一定效果，進(jìn)一步消除了重力殘差中的潮汐信號(hào)，有助于更準(zhǔn)確地提取地球內(nèi)部引起的重力變化，以獲得地球內(nèi)部動(dòng)力學(xué)信息，為地震預(yù)測(cè)服務(wù)。海潮改正后仍然殘留潮汐信號(hào)，可能是因?yàn)樵谶M(jìn)行固體潮改正時(shí)使用的是理論值，且海潮改正基于的是海潮模型。因此計(jì)算的固體潮、海潮負(fù)荷與實(shí)際真實(shí)值仍有誤差，這些均可能導(dǎo)致海潮改正后的重力非潮汐變化數(shù)據(jù)仍殘留潮汐信號(hào)。

需要指出的是，以上重力非潮汐變化改正過程中，零漂改正具有一定的不確定性，主要是因?yàn)槭躊ET/gphone相對(duì)重力儀自身結(jié)構(gòu)特性的限制，很難完全真實(shí)地?cái)M合儀器的零漂。尤其當(dāng)儀器出現(xiàn)諸如停電或儀器故障等原因?qū)е峦y(cè)一段時(shí)間，待儀器恢復(fù)觀測(cè)后，儀器的零漂非線性零漂特征明顯，故真實(shí)擬合儀器零漂更為困難，以致在數(shù)據(jù)處理過程中，可能會(huì)把真實(shí)的重力變化異常當(dāng)做零漂進(jìn)行改正。