趙亞博 胡敏章 韋 進(jìn) 江 穎 張曉彤 劉子維

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

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

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

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

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

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

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

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

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

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

表1 重力臺(tái)站觀測(cè)時(shí)長(zhǎng)與數(shù)據(jù)完整度Tab.1 Observation duration and data completeness of gravity stations

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

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

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

3)對(duì)分鐘值數(shù)據(jù)進(jìn)行改正后,再次對(duì)數(shù)據(jù)進(jìn)行降采樣處理,得到小時(shí)值數(shù)據(jù)。最后將經(jīng)過(guò)處理的各月小時(shí)值數(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 觀測(cè)數(shù)據(jù)處理

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

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

gre=gobsf-gth-aast-gotl

(1)

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

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

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

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

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

fM2=δDDW/δM2

(2)

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

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

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

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

(3)

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

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

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

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

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

由圖4可以看出,格值系數(shù)和重力殘差均方根存在二次函數(shù)對(duì)應(yīng)關(guān)系,重力殘差均方根先隨著格值系數(shù)的增加而減小,在到達(dá)一個(gè)最小值后隨著格值系數(shù)的增大而增大。因此,格值系數(shù)存在一個(gè)極值使得重力殘差均方根值最小。對(duì)該二次函數(shù)進(jìn)行求導(dǎo)得出其極值,對(duì)應(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 潮汐因子精度分析

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

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

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

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

對(duì)gPhone重力儀格值系數(shù)標(biāo)定前后的觀測(cè)數(shù)據(jù)分別扣除固體潮、大氣、海潮、儀器漂移項(xiàng)得到重力殘差,然后進(jìn)行振幅頻譜分析,以頻率1.000 CPT作為截止頻率進(jìn)行高通濾波得到殘差振幅譜。圖5為西安站標(biāo)定前后殘差振幅高通濾波頻率譜,圖6為9個(gè)臺(tái)站標(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個(gè)臺(tái)站經(jīng)過(guò)格值系數(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)過(guò)格值系數(shù)改正之后,重力殘差明顯降低,高通濾波值分布在4 μGal左右。

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

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

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

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

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

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

3 結(jié) 語(yǔ)

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

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

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

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