宜昌臺超導重力觀測中的海潮與大氣負荷效應研究

魏工哲 劉子維,2,3 江 穎,2,3 張曉彤,3 周 浩,3 李 暉

1 中國地震局地震研究所地震大地測量重點實驗室,武漢市洪山側路40號,430071 2 防災科技學院,河北省三河市學院街465號,065201 3 武漢引力與固體潮國家野外科學觀測研究站,武漢市洪山側路40號,430071

地球重力場隨時間的變化可反映地球質量遷移、密度變化等多因素地球物理過程。超導重力儀是一種精密連續重力觀測儀器,在精確觀測重力潮汐信號的同時,也可記錄到大氣、海潮、極移及地震等引起的重力變化[1]。為了對重力觀測數據進行處理與分析,需要識別和分離出各種地球物理信號。重力固體潮作為觀測到的主要地球物理時變現象,其變化幅度可達200～500 μGal。在扣除理論重力固體潮效應后,全球尺度的觀測殘差主要受海潮和大氣潮兩大負荷源影響,而海潮和大氣潮也會受到太陽、月亮等天體引潮力的作用,與固體潮具有相同的頻率特征,無法直接從固體潮觀測中有效分離[2-3]。因此,研究海潮和大氣潮重力負荷效應有助于精確實施重力固體潮觀測改正,為水文地質、地球內核震蕩、地震預警等研究提供更為精確的數據。

本文基于布設在宜昌地震臺的iGrav026超導重力儀在2017-01-01～12-31期間采集到的連續重力觀測數據(圖1(a))和臺站氣壓數據(圖1(b)),結合海潮模型、全球大氣模型等資料,研究海潮和大氣負荷產生的重力效應,從而對超導重力觀測實施精密改正。

圖1 重力潮汐觀測數據預處理

1 數據初步處理和分析

1.1 預處理

使用Tsoft軟件[4]對超導數據進行預處理。采用數字濾波方法將數據采樣率由1 s降為1 min,利用“移去-恢復”技術修正數據。首先采用WDD模型[5]扣除理論固體潮(圖1(c)),然后采用經驗大氣導納值扣除大氣影響,導納值取-3.0 nm·s-2·hPa-1。在殘差數據上進行尖峰、間斷、地震、階躍等干擾的預處理。

1.2 極移改正

由IAGBN絕對重力觀測數據處理標準[6]可知,極移重力影響可表示為:

δg=-1.164ω2a2sinφcosφ

(xpcosλ-ypsinλ)×109

(1)

式中,λ、φ分別表示測點經度、緯度;xp、yp分別表示2個地球自轉參數;ω表示地球自轉速度,取7.292 115×10-5rad·s-1;a表示地球長半軸,取6 378 136 m。IERS(國際地球自轉服務)網站已給出詳盡的地球自轉和極移參數,極移參數文件每天以數據形式給出,因此在計算得到極移重力效應序列后,通過樣條插值得到與原始觀測序列采樣率相當的數據,結果如圖1(d)所示。經計算發現,在2017年內極移效應引起的重力變化最大可達約40 nm·s-2。

1.3 長周期趨勢改正

經過初步的理論重力固體潮、大氣以及極移重力影響改正后,數據中存在一個長期升高趨勢,可能與儀器漂移、水文負荷、地殼垂直位移和地球內部物質遷移有關[7],而這些因素在本文研究中為干擾項。為了盡量削弱該項影響,并避免過擬合造成的信息丟失,本文采用二次多項式方法[8-9]擬合重力長期趨勢:

d(t)=b1+b2t+b3t2

(2)

式中,b1、b2、b3分別表示趨勢多項式的常數項、一次項系數和二次項系數。通過最小二乘擬合方法可以求解各項參數,計算獲得的重力長期趨勢如圖1(e)所示。

對預處理后的重力潮汐數據扣除極移重力影響和長周期趨勢項,得到初步的重力殘差(圖1(f))。由圖可知,數據預處理效果較為理想,數據中殘存較為明顯的周期項,這可能與海潮或大氣負荷有關。基于圖1(f)中殘差數據,恢復先前扣除的理論固體潮和氣壓影響,再將數據降為1 h采樣,作為調和分析的輸入數據。

1.4 調和分析

觀測數據的調和分析是將觀測信號分解為一系列周期性函數之和,求解各個潮波的潮汐特征參數(即振幅因子δ和相位滯后Δφ)并計算其估算精度的過程。振幅因子δ即各個潮波的觀測振幅與理論振幅之比,為3個勒夫數的線性組合;相位滯后Δφ即各個潮波的實測相位與理論相位之差。Eterna[10]是使用最為廣泛的調和分析軟件,其主要原理是對觀測數據使用最小二乘擬合估計潮汐參數。本文采用該軟件對超導重力潮汐觀測數據進行調和分析,結果見表1。

表1 宜昌臺超導重力觀測數據調和分析結果

2 海潮負荷效應

2.1 海潮負荷理論

海潮負荷計算首先從負荷格林函數理論[11]出發,通過將海潮潮高和負荷格林函數在球面進行褶積積分[12-13]得到海潮負荷,具體公式為:

LO(θ,λ,t)=ρw?H(θ′,λ′,t)G(φ,A)ds

(3)

式中,LO(θ,λ,t)為海潮負荷,θ、λ為計算點的余緯、經度;ρw為海水密度;H(θ′,λ′,t)為瞬時海潮潮高,θ′、λ′為海潮負荷點的余緯、經度;G(φ,A)為海潮重力格林函數,φ為負荷點至計算點的角距,A為負荷點至計算點的方向角;s為地球表面,ds為積分面元。

2.2 海潮負荷計算

全球性或區域性海潮模型主要有FES、CSR、GOT、NAO、TPXO、EOT、DTU、HAMTIDE和OSU系列等[14]。其中,NAO.99b是日本國立天文臺在2000年發布的海潮模型,采用同化模型方式構建,主要數據來源為T/P測高衛星以及日本驗潮站資料,格網分辨率為0.5°[14];osu.Chinasea.2010來自于俄勒岡州立大學,吸收了中國南海、臺灣海峽和泰國灣等地區驗潮站資料,在中國近海具有較高精度[15]。研究表明[16],NAO.99b模型在中國地區進行海潮改正是有效的。同時,為進一步提高海潮改正的效果,應結合中國近區海潮模型[13]。因此,本文結合全球NAO.99b模型與近區osu.Chinasea.2010模型,使用SPOTL軟件[17]計算宜昌地震臺海潮重力負荷效應。

表2為采用上述海潮模型計算的周日、半日頻段8個主要潮波重力負荷振幅。由表可知,在8個潮波中,K1波海潮重力負荷影響最大,幅值為4.142 nm·s-2;其次為O1波和M2波,幅值分別為3.735 nm·s-2和3.662 nm·s-2;其他潮波振幅均在2 nm·s-2以下,因此后續分析中重點考察K1、O1、M2三個主要潮波。從不同頻段的影響來看,周日頻段(O1、K1)的海潮負荷振幅大于半日頻段(M2),這是由于宜昌臺站地處中國南方,與區域性的海潮變化性質有關,這與孫和平等[18]的結論相符。

表2 宜昌臺海潮重力負荷計算結果

2.3 海潮負荷改正

圖2(a)為通過上述模型計算的海潮重力負荷時序,可以看出,宜昌臺海潮重力負荷周期性與理論固體潮(圖1(c))基本一致,最大變幅約為25 nm·s-2。圖2(b)為原始觀測數據經過理論固體潮改正和海潮負荷改正的殘差。

圖2 海潮重力負荷效應及改正

為了評價海潮負荷改正的效果,在對原始重力觀測數據進行海潮負荷改正后,再次進行調和分析,改正前后的K1、O1、M2三個主要潮波潮汐因子如表3所示。經海潮負荷改正后,O1、K1、M2潮波振幅因子分別減小1.2%、0.9%、0.6%;O1波相位滯后略有增大,K1、M2波相位滯后分別減小81.8%、24.6%。潮汐因子經海潮負荷改正后可構建實測固體潮模型,用于宜昌臺精確重力固體潮改正。

表3 海潮負荷改正前后O1、K1、M2波潮汐因子對比

3 大氣重力負荷效應

3.1 時域大氣重力導納值

一般情況下,如果可以同時獲取臺站連續重力觀測和同址氣壓觀測數據,便可通過時域導納方法進行氣壓改正。時域導納方法就是將氣壓觀測和重力殘差時間序列作一元線性回歸分析,獲得時域大氣重力導納值[19]。分別對每月重力與氣壓觀測序列進行回歸分析,得到2017-01～12大氣重力導納值(表4)。除1月、7月、8月、9月、10月外,其余7個月份相關系數基本在70%以上,取平均值得到導納值為-3.219 3 nm·s-2·hPa-1,與經驗導納值較為接近。采用該導納值計算氣壓重力負荷效應,結果如圖3(a)所示。由圖可知,該方法計算的大氣重力負荷影響在1 a內峰對峰振幅變化約為120 nm·s-2。

表4 宜昌臺2017年不同月份大氣重力導納值

圖3 大氣重力負荷效應及改正

3.2 全球大氣負荷

為了從重力潮汐觀測中更精確地扣除大氣重力信號,在考慮臺站大氣重力信號的同時,還需顧及全球性的大氣重力效應[20]。本文選取全球大氣模型來計算近區和遠區大氣重力效應。

3.2.1 大氣負荷理論

與海潮負荷計算類似,大氣重力負荷效應計算也從負荷格林函數理論出發,通過對大氣重力格林函數與全球大氣數據進行褶積積分得到大氣重力效應[21],具體公式為:

LP(θ,λ,t)=?P(θ′,λ′,t)G(φ)ds

(4)

式中,LP(θ,λ,t)為大氣重力負荷,θ、λ為計算點的余緯、經度;P(θ′,λ′,t)為大氣負荷點的大氣壓強,θ′、λ′為大氣負荷點的余緯、經度;G(φ)為大氣重力格林函數,φ為負荷點至計算點的角距;s為地球表面,ds為積分面元。

3.2.2 大氣負荷計算

本文選用的全球大氣模型為ERA-Interim,由ECMWF(歐洲中期天氣預報中心)發布(https:∥compreps.ecmwf.int/)。ERA-Interim再分析數據是在ERA40資料基礎上,結合新的四維數據同化處理方法,應用更多衛星和地面觀測資料,空間分辨率為0.75°×0.75°,最大時間分辨率是6 h。ERA-Interim模型包括垂向的37個氣壓位分層,由1 000 hPa上升至1 hPa(約37 km);數據時刻為00:00、06:00、12:00、18:00(UTC時間),其中00:00與12:00為分析數據,06:00和18:00為預報數據[22]。

選擇0.5°為近區與遠區分區閾值,基于ERA-Interim大氣模型,計算得到2017-01～12期間宜昌地震臺全球大氣重力負荷效應,結果如圖3(b)所示。圖中綠線為全球大氣重力效應,包括遠區大氣重力效應(藍線)和近區大氣重力效應(紅線)。計算結果顯示,全球大氣重力效應呈現季節性變化特點。氣壓主要受季節溫度變化影響,春冬高、夏秋低,由于大氣對重力的影響主要為負向吸引,因此重力效應呈現春冬低、夏秋高的特點,其1 a內峰對峰振幅變化約為80 nm·s-2。在全部影響中,近區大氣重力效應占主要成分,其峰對峰幅值變化約為160 nm·s-2,且變化趨勢與時域導納法計算的大氣重力效應基本一致,兩者數值的差異反映出大氣模型在局部地區精度存在欠缺;遠區大氣重力效應相對較小,變化趨勢與近區重力效應呈負相關,在3～7月期間呈下降趨勢,在8～11月期間呈上升趨勢,其1 a內峰對峰幅值變化約為70 nm·s-2。

3.3 大氣負荷改正

原始觀測數據在扣除海潮重力負荷效應和理論重力固體潮基礎上,進一步扣除大氣重力負荷效應。在對由大氣重力導納值計算的重力效應進行改正后,顧及到遠區大氣負荷影響,進一步使用大氣模型計算的遠區大氣重力負荷效應進行改正,得到全球大氣重力負荷改正結果,重力殘差如圖3(c)所示。圖中黑線為采用大氣重力導納值計算的重力效應進行改正的重力殘差,在3～6月殘差數據中存在幅值約40 nm·s-2的下降趨勢;圖中藍線為考慮遠區大氣負荷的全球大氣重力負荷改正殘差。由前文計算結果可知,遠區大氣負荷在3～6月會引起負重力變化,在考慮該項因素后,殘差數據中3～6月下降趨勢得到有效改正。在最終重力殘差結果中,1～2月變化趨勢可能是由超導重力儀安裝后的非線性零漂特征所致;6～10月呈現先上升后下降的趨勢,變幅約為40 nm·s-2,該部分趨勢變化主要反映陸地水[7]、地下水[23]等方面的影響。

4 結 語

本文對宜昌地震臺超導重力觀測中的海潮和大氣重力負荷影響進行計算和分析,主要得到以下結論:

1)宜昌臺地處中國大陸中部,距離海岸線相對較遠,海潮重力影響較小,海潮重力負荷最大變幅為25 nm·s-2,對于高精度的超導重力觀測需要適當考慮海潮重力負荷的影響。在經過海潮負荷改正后,潮汐因子有一定改善。

2)宜昌臺大氣時域導納值為-3.219 3 nm·s-2·hPa-1,與經驗值較為接近。通過理論計算發現,近區大氣重力效應變幅達160 nm·s-2,為大氣重力影響的主要部分;而遠區大氣重力效應變幅達70 nm·s-2,能夠部分解釋觀測數據中的負重力變化。在進行超導重力觀測數據精細化處理時,應考慮全球大氣重力效應。

應當注意的是,超導重力儀在安裝后可能出現非線性零漂,且長期零漂行為與開始安裝時呈現不同的規律,而控制零漂的常用方法是與絕對重力儀進行同址觀測[24]。對宜昌臺超導重力儀進行絕對重力比對觀測,獲取精確格值并研究其長期漂移,是下一步將要開展的重點工作。

致謝：感謝中國科學院精密測量科學與技術創新研究院陳曉東副研究員對本文內容和格式提出的建議和幫助。