瓊中臺重力非潮汐變化氣壓與海潮負荷改正

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

海南省地震局，?？?570203

0 引言

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

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

(1)

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

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

(2)

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

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

(3)

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

圖1 瓊中地震臺位置

1 原始資料分析

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

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

2 潮汐分析

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

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

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

3 重力非潮汐變化改正

3.1 零漂改正

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

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

3.2 固體潮改正

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

圖3 瓊中臺零漂固體潮改正后的重力殘差時序曲線

3.3 氣壓改正

表2 不同機構(gòu)的大氣重力導納值

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

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

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

3.4 海潮改正

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

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

(4)

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

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

(5)

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

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

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

圖4 瓊中臺重力潮汐變化零漂固體潮氣壓海潮改正時序圖

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

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

3.5 小結(jié)

根據(jù)以上計算的結(jié)果，將瓊中臺重力非潮汐變化氣壓與海潮負荷改正值繪制在同一圖中(圖5)，可知氣壓改正值約為10×10-8m/s2，海潮改正值約為5×10-8m/s2。

圖5 瓊中臺重力氣壓負荷和海潮負荷改正值

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

圖6 瓊中臺海潮改正后的重力非潮汐年變序列

4 結(jié)論

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

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

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

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

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