以潮汐模型為約束的gPhone重力儀格值系數(shù)標(biāo)定方法改進

趙亞博 胡敏章 韋 進 江 穎 張曉彤 劉子維

1 中國地震局地震研究所,武漢市洪山側(cè)路40號,430071

多年的連續(xù)工作使gPhone彈簧相對重力儀發(fā)生老化,導(dǎo)致觀測值和實際重力值之間存在一定誤差,通常用一個無量綱的常數(shù)(格值系數(shù))對觀測值進行改正[1-3]。

gPhone重力儀格值系數(shù)的標(biāo)定方法主要包括長短重力基線法、垂直重力基線法、同址重力比測法、重力潮汐參數(shù)標(biāo)定法以及人工加速度法。重力基線法是指使用儀器在已知段差的重力基線兩端往返觀測,將實測段差與已知段差進行對比得到儀器的格值系數(shù)[4],該方法需要固定不變的重力段差,實際應(yīng)用存在一定誤差。同址重力比測法利用絕對重力儀(AG)或超導(dǎo)重力儀(SG)與gPhone重力儀進行同時、同址觀測,基于最小二乘原理標(biāo)定gPhone重力儀格值系數(shù)[3-6],但我國gPhone臺站較多(約62個)、分布較廣[4],利用AG或SG比測標(biāo)定難度較大、耗時較長、成本較高[4]。目前國內(nèi)也沒有可用于重力儀標(biāo)定的人工加速度重力平臺[4]。因此,以理論潮汐模型為約束,對gPhone重力儀進行格值系數(shù)標(biāo)定是目前最經(jīng)濟、便捷的方法。

通常使用M2潮波標(biāo)定法確定gPhone重力儀格值系數(shù),但該方法只考慮了M2潮波帶來的影響,忽視了其他振幅較大的潮波的影響,如O1、P1、K1、N2、S2潮波,因此存在一定的局限性,改正結(jié)果也不符合海潮負荷改正規(guī)律[3]。

本文使用理論潮汐模型進行約束,對臺站gPhone重力儀進行格值系數(shù)標(biāo)定,并使用潮汐因子中誤差等指標(biāo)評定標(biāo)定結(jié)果。最后,分析我國臺站格值系數(shù)隨經(jīng)緯度的變化規(guī)律。

1 數(shù)據(jù)與格值系數(shù)標(biāo)定

本文選用沿海隨緯度變化以及由沿海至內(nèi)陸隨經(jīng)度變化的2條測線共計9個臺站2018年的觀測資料,使用理論潮汐模型進行約束,對臺站gPhone重力儀進行格值系數(shù)標(biāo)定,并使用潮汐因子相對誤差、重力殘差高通濾波均方根和改正前后的重力殘差矢量等指標(biāo)評定標(biāo)定結(jié)果(圖1)。

圖1 格值系數(shù)標(biāo)定流程Fig.1 The flow of scale factor calibration

1.1 觀測數(shù)據(jù)及其預(yù)處理

為去除不同時段地震、臺風(fēng)等因素造成的影響,使最終分析結(jié)果能夠更好地體現(xiàn)格值系數(shù)的空間分布規(guī)律,選用大連等9個gPhone連續(xù)重力觀測臺站(沿海臺站6個、內(nèi)陸臺站3個)同一時間(2018年)的觀測數(shù)據(jù)。表1給出了各臺站的觀測時長和數(shù)據(jù)完整度。

表1 重力臺站觀測時長與數(shù)據(jù)完整度Tab.1 Observation duration and data completeness of gravity stations

受到地震、降水、環(huán)境噪聲等因素的影響,觀測數(shù)據(jù)中可能出現(xiàn)偏離背景值、粗差較大的情況,并對格值系數(shù)的標(biāo)定結(jié)果產(chǎn)生影響。因此,需要對重力數(shù)據(jù)和氣壓數(shù)據(jù)進行預(yù)處理:

1)使用低通濾波器對原始觀測數(shù)據(jù)進行濾波,去除觀測數(shù)據(jù)中的高頻部分,得到分鐘采樣值。為了保證數(shù)據(jù)修正過程的獨立性和正確性,對觀測數(shù)據(jù)按月份分別進行處理[7]。

2)使用Tsoft軟件[8]計算理論固體潮模型,扣除理論固體潮得到觀測殘差,對觀測數(shù)據(jù)中出現(xiàn)的間斷、突跳、階躍等粗差進行手動改正;然后使用一階多項式擬合對觀測殘差進行零漂改正;最后再加回理論固體潮,得到經(jīng)過預(yù)處理后的分鐘值數(shù)據(jù)。

3)對分鐘值數(shù)據(jù)進行改正后,再次對數(shù)據(jù)進行降采樣處理,得到小時值數(shù)據(jù)。最后將經(jīng)過處理的各月小時值數(shù)據(jù)連接在一起,得到標(biāo)準(zhǔn)格式的數(shù)據(jù)[9]。圖2給出西安站數(shù)據(jù)預(yù)處理流程。

圖2 西安站數(shù)據(jù)預(yù)處理流程Fig.2 The data preprocessing flow of Xi’an station

1.2 觀測數(shù)據(jù)處理

由于gPhone重力儀在出廠前一般都經(jīng)過標(biāo)定,在投入使用后的格值系數(shù)標(biāo)定結(jié)果均在1.009 2±0.004范圍內(nèi)[1],因此本文采用范圍為0.9～1.1、步長為0.001的一系列格值系數(shù)分別對數(shù)據(jù)進行處理。

對經(jīng)過格值系數(shù)改正后的觀測數(shù)據(jù)進行大氣改正以及海潮改正:

gre=gobsf-gth-aast-gotl

(1)

式中,gre為大氣改正與海潮改正后的殘差值,gobs為經(jīng)過預(yù)處理后的原始數(shù)據(jù),f為重力儀的格值系數(shù),a為大氣導(dǎo)納值,ast為儀器測得的大氣壓值,gth為理論固體潮重力潮汐,gotl為計算得到的海潮重力負荷。

采用移去-恢復(fù)法[7]進行改正:1)使用原始重力值乘以格值系數(shù)得到實測重力潮汐值,然后使用Tsoft計算DDW理論固體潮模型重力潮汐模擬信號,在實測重力潮汐變化中扣除理論潮汐信號得到殘差;2)利用Tsoft人機交互界面對殘差進行改正,包括去除尖峰、階躍、掉格和地震等干擾信號的影響,得到改正殘差;3)采用理論大氣導(dǎo)納值(-3.4×10-9m/(s2·hPa))對得到的殘差進行大氣改正,然后再恢復(fù)扣除的DDW理論重力潮汐信號,得到經(jīng)過第1次處理后的重力潮汐信號;4)使用ETERNA軟件[10]對得到的重力潮汐信號進行調(diào)和分析,得到重力潮汐的潮波調(diào)和參數(shù),即重力合成潮;5)利用得到的合成潮潮汐參數(shù)作為模型潮汐參數(shù),重復(fù)步驟1)～2),求得觀測殘差,并與大氣信號作聯(lián)合回歸分析,得到大氣導(dǎo)納值a;6)采用新的大氣導(dǎo)納值進行大氣改正;7)使用海潮計算軟件SPOTL[11],選用全球模型為FES2004、近海模型為osu.chinasea.2010的組合海潮模型計算海潮在臺站引起的海潮重力負荷效應(yīng),隨后扣除海潮負荷引起的重力影響,得到觀測重力殘余;8)將合成重力潮汐模型加回,得到經(jīng)過改正之后的重力觀測數(shù)據(jù)(圖3)。

圖3 西安站數(shù)據(jù)處理流程Fig.3 The data processing flow of Xi’an station

1.3 以理論潮汐模型為約束的格值系數(shù)標(biāo)定方法改進

傳統(tǒng)的以M2潮波為約束的標(biāo)定方法僅考慮振幅最大的M2潮波帶來的影響。其對經(jīng)過預(yù)處理之后的觀測數(shù)據(jù)進行潮汐分析,得到M2潮波的潮汐因子,隨后使用理論潮汐模型(如DDW固體潮模型)M2潮波潮汐因子對觀測數(shù)據(jù)進行改正,將二者的比值作為格值系數(shù)[3]:

fM2=δDDW/δM2

(2)

式中,fM2為格值系數(shù),δDDW為DDW模型M2潮波潮汐因子,δM2為觀測數(shù)據(jù)M2潮波潮汐因子。

M2潮波標(biāo)定法步驟簡單,但由于忽略了其他振幅較大的潮波對大陸的影響,因此得到的格值系數(shù)結(jié)果存在一定誤差。本文在M2潮波標(biāo)定法的基礎(chǔ)上,充分考慮振幅較大的6個潮波(O1、P1、K1、N2、M2、S2)對觀測數(shù)據(jù)產(chǎn)生的影響,并以此為約束對格值系數(shù)進行標(biāo)定。改進方法主要分2步進行:

1)基于潮汐因子均方根(factor RMS,FRMS)確定格值系數(shù)范圍。使用FRMS評價理論潮汐模型DDW和經(jīng)過改正之后的觀測固體潮之間的各潮波累積誤差:

(δmA0sinφm-δ0A0sinφ0)2]}

(3)

式中,δ0、A0和φ0分別為理論固體潮模型的振幅因子、理論振幅和理論相位(0°),δm和φm分別為實測數(shù)據(jù)經(jīng)過大氣和固體潮改正之后的振幅因子和相位,I為潮波的分波數(shù)量。

使用ETERNA軟件對經(jīng)過格值系數(shù)、大氣、海潮改正之后的數(shù)據(jù)進行潮汐調(diào)和分析,得到各重力站不同格值系數(shù)下的潮汐因子;然后利用式(3)計算不同格值系數(shù)各潮波累積均方根;最后篩選FRMS最小區(qū)間的格值系數(shù)作為均方根檢驗的結(jié)果[6]。表2為經(jīng)過均方根檢驗后的格值系數(shù)取值范圍。

表2 經(jīng)過均方根檢驗后的格值系數(shù)取值范圍Tab.2 Value range of scale factors after RMS test

2)基于重力殘差均方根選定最優(yōu)格值系數(shù)。經(jīng)過均方根檢驗之后,進一步對格值系數(shù)標(biāo)定結(jié)果進行檢驗。為分析觀測數(shù)據(jù)的半日波和周日波信號,選用潮汐周日波最高頻率1.000 CPT作為截止頻率進行高通濾波,得到觀測重力殘差高通濾波值,之后計算重力殘差高通濾波值的均方根,即為gPhone重力儀的觀測數(shù)據(jù)精度,最終得到格值系數(shù)與觀測重力殘差均方根一一對應(yīng)的結(jié)果(圖4)。

圖4 西安臺格值系數(shù)與觀測重力殘差均方根Fig.4 RMS of observed gravity residuals and scale factor at Xi’an station

由圖4可以看出,格值系數(shù)和重力殘差均方根存在二次函數(shù)對應(yīng)關(guān)系,重力殘差均方根先隨著格值系數(shù)的增加而減小,在到達一個最小值后隨著格值系數(shù)的增大而增大。因此,格值系數(shù)存在一個極值使得重力殘差均方根值最小。對該二次函數(shù)進行求導(dǎo)得出其極值,對應(yīng)的格值系數(shù)即為固體潮模型約束下的標(biāo)定結(jié)果(表3)。

表3 格值系數(shù)標(biāo)定結(jié)果Tab.3 Results of scale factor calibration

2 標(biāo)定結(jié)果有效性分析

2.1 潮汐因子精度分析

采用時間域內(nèi)改正的方法,使用經(jīng)過檢驗之后的格值系數(shù)標(biāo)定結(jié)果對預(yù)處理數(shù)據(jù)進行改正;之后進行大氣改正與海潮改正,得到改正后的固體潮重力信號;再對這些信號進行調(diào)和分析,得到經(jīng)過改正后的觀測重力數(shù)據(jù)振幅因子和相位。為驗證格值系數(shù)標(biāo)定結(jié)果的準(zhǔn)確性,將改正得到的潮汐因子與理論固體潮模型DDW進行對比,其中DDW模型是由Dehant等[12]提出的使用PREM地球模型[13]計算出的流體非靜力非彈性地球潮汐模型。采用(觀測改正結(jié)果-DDW固體潮模型)/觀測改正結(jié)果的比值來計算相對誤差,并以此評價標(biāo)定結(jié)果的精度。表4給出各臺站觀測模型M2潮波潮汐因子、中誤差以及與DDW模型之間的相對誤差。

表4 各臺站標(biāo)定結(jié)果相對誤差Tab.4 Relative error of calibration results for each station

可以看出,經(jīng)過格值系數(shù)、大氣改正、海潮改正之后,各臺站M2潮波因子潮波中誤差在0.000 4～0.001 2之間,說明本文所使用的臺站gPhone重力儀觀測精度已經(jīng)接近早期超導(dǎo)重力儀的水平[6]。比較經(jīng)過格值系數(shù)標(biāo)定前后的潮汐因子與DDW模型之間的相對誤差可以看出,經(jīng)過格值系數(shù)改正之前的M2潮波平均相對誤差為0.590 1%,而經(jīng)過格值系數(shù)改正之后,M2潮波平均相對誤差為0.215 2%。經(jīng)過格值系數(shù)改正,M2潮波潮汐因子的精度有顯著提高。

2.2 標(biāo)定前后殘差振幅比較

對gPhone重力儀格值系數(shù)標(biāo)定前后的觀測數(shù)據(jù)分別扣除固體潮、大氣、海潮、儀器漂移項得到重力殘差,然后進行振幅頻譜分析,以頻率1.000 CPT作為截止頻率進行高通濾波得到殘差振幅譜。圖5為西安站標(biāo)定前后殘差振幅高通濾波頻率譜,圖6為9個臺站標(biāo)定前后重力殘差均方根。

圖5 西安站格值系數(shù)標(biāo)定前后殘差高通濾波振幅譜Fig.5 Residual high-pass filtered amplitude spectrum before and after scale factor calibration at Xi’an station

圖6 9個臺站經(jīng)過格值系數(shù)改正前后重力殘差均方根Fig.6 The RMS of gravity residuals before and after correction of scale factors at nine stations

由圖5可知,在標(biāo)定之前,西安站的重力殘差高通濾波值分布在7 μGal左右;經(jīng)過格值系數(shù)改正之后,重力殘差明顯降低,高通濾波值分布在4 μGal左右。

從圖6看出,經(jīng)過格值系數(shù)改正之后,各臺站的殘差振幅都明顯減小,且均在2.5 μGal以下,均值為1.690 2 μGal,與韋進等[3]使用FG5絕對重力儀比測方法得到的2 μGal的精度相當(dāng)。改正之后的殘差振幅平均降低15%,其中溧陽站格值系數(shù)改正效果最明顯,殘差降低58.68%。

2.3 標(biāo)定前后重力殘差與殘余矢量比較

將重力臺站按照經(jīng)度大小自西向東排序,計算海潮改正前后的重力殘差矢量和重力殘余矢量,并與M2潮波標(biāo)定法結(jié)果進行對比(表5)。

表5 臺站經(jīng)緯度與重力殘差/殘余矢量Tab.5 Gravity residual vector and latitudes and longitudes of the stations

從表5中可以看出,本文使用以潮汐模型為約束的標(biāo)定方法的結(jié)果不僅符合重力殘差矢量由于海潮影響而出現(xiàn)的隨經(jīng)度自西向東逐漸減小的特征,同時也符合由于海潮負荷改正導(dǎo)致重力殘余矢量小于重力殘差矢量的規(guī)律。M2潮波標(biāo)定法雖然也滿足重力殘差隨經(jīng)度變化的規(guī)律,但是出現(xiàn)了海潮改正后重力殘余矢量大于改正前重力殘差矢量的現(xiàn)象,改正結(jié)果不符合海潮負荷改正規(guī)律。因此,本文使用的以潮汐模型為約束的標(biāo)定法要優(yōu)于M2潮波標(biāo)定法的結(jié)果,且標(biāo)定結(jié)果能用于我國大陸海潮負荷相關(guān)研究。

綜上所述,經(jīng)過格值系數(shù)標(biāo)定之后,gPhone重力儀觀測數(shù)據(jù)M2潮波潮汐因子的中誤差在0.000 4～0.001 2之間,與DDW模型之間的相對誤差也提高了68.82%,達到0.215 2%;同時,觀測重力殘差振幅高通濾波值均方根平均降低15%,數(shù)據(jù)精度和質(zhì)量明顯提高。標(biāo)定結(jié)果符合我國大陸重力殘差矢量自西向東逐漸減小的分布規(guī)律,也符合海潮負荷改正后重力殘余矢量減小的規(guī)律,改正結(jié)果優(yōu)于M2潮波標(biāo)定法。

3 結(jié) 語

本文以潮汐模型為約束,提出一種改進的gPhone重力儀格值系數(shù)標(biāo)定方法,并以9個臺站的gPhone重力儀數(shù)據(jù)為例驗證該方法的精度,主要得到如下結(jié)論:

1)本文使用潮汐模型作為約束對格值系數(shù)進行標(biāo)定,標(biāo)定之后的潮汐因子與DDW模型潮汐因子相對誤差平均為0.215 2%,精度相比標(biāo)定前提高68.82%;重力殘差振幅均方根均在2.5 μGal以下,均值為1.690 2 μGal,相對于改進前平均降低15%。

2)本文改正結(jié)果符合大陸重力殘差分布規(guī)律以及海潮負荷改正規(guī)律,且優(yōu)于M2潮波標(biāo)定法改正結(jié)果。

3)本文使用的以潮汐模型作為約束的格值系數(shù)標(biāo)定方法不僅可以在不依賴AG或SG數(shù)據(jù)比測的情況下對儀器進行標(biāo)定,且標(biāo)定后的潮汐因子中誤差和殘差振幅都符合實際科研工作需要,同時也可以在統(tǒng)一的潮汐基準(zhǔn)下對儀器進行標(biāo)定,為地震監(jiān)測、重力觀測等科學(xué)研究提供支持。