基于滑動(dòng)奇異譜分析方法的重力固體潮提取

王建波 王晉南 孫 潔 齊鑫敏 劉相威 高文宗

1 江蘇海洋大學(xué)海洋技術(shù)與測(cè)繪學(xué)院，江蘇省連云港市蒼梧路59號(hào)，222005 2 菏澤市測(cè)繪研究院，山東省菏澤市珠江路1717號(hào)，274002 3 山東科技大學(xué)測(cè)繪與空間信息學(xué)院，青島市前灣港路579號(hào)，266590

雖然CG-5重力儀自帶軟件可提供固體潮改正值，但理論模型是建立在理想條件和固定參數(shù)的基礎(chǔ)上的，其計(jì)算得到的重力固體潮理論值與實(shí)際值之間存在一定差異。如何快速?gòu)南鄬?duì)重力觀測(cè)數(shù)據(jù)中分離重力固體潮改正，對(duì)相對(duì)重力觀測(cè)具有重要意義。

SSA方法可以對(duì)原先一維時(shí)間序列時(shí)延排序得到的時(shí)滯矩陣進(jìn)行分解和重構(gòu)，提取出原序列中代表不同成分的信號(hào)，如周期信號(hào)、趨勢(shì)信號(hào)、噪聲信號(hào)等。郭金運(yùn)等[1]采用該方法提取重力固體潮，并進(jìn)行周期為10 d的靜態(tài)相對(duì)重力固體潮提取實(shí)驗(yàn)，取得很好的效果。雖然相對(duì)重力儀靜態(tài)零漂線性度很好，但零漂總體呈下降趨勢(shì)[2]。另外，相對(duì)重力儀在長(zhǎng)時(shí)間靜態(tài)重力觀測(cè)時(shí)會(huì)受到觀測(cè)時(shí)間長(zhǎng)、零漂線性度改變、非線性誤差[3]以及復(fù)雜環(huán)境噪聲的影響，這些影響因素會(huì)降低SSA方法提取重力固體潮的精度。為此，本文提出滑動(dòng)SSA方法，通過(guò)設(shè)置奇異譜分析數(shù)據(jù)長(zhǎng)度時(shí)間窗口，解決長(zhǎng)時(shí)間靜態(tài)相對(duì)重力觀測(cè)時(shí)重力固體潮的提取問題。

1 原理與方法

1.1 SSA方法、重力固體潮粗差識(shí)別及修復(fù)

SSA方法的原理、粗差識(shí)別及迭代插值方法可參見文獻(xiàn)[1]。用SSA方法提取重力固體潮時(shí)，一般先利用ω-correlation方法[4]對(duì)重建成分(RCs)進(jìn)行頻譜分析和重構(gòu)。而本文理論重力固體潮值選取與之不同，分為2種：1)對(duì)實(shí)測(cè)重力數(shù)據(jù)采用Eterna調(diào)和分析獲得潮波參數(shù)，然后用TSoft軟件[5]計(jì)算得到理論重力固體潮值；2)用CG-5重力儀自帶軟件進(jìn)行計(jì)算。

粗差識(shí)別時(shí)，利用SSA方法獲得原靜態(tài)相對(duì)重力觀測(cè)時(shí)間序列的主成分，計(jì)算原靜態(tài)相對(duì)重力觀測(cè)時(shí)間序列與主成分之間的殘差，分析殘差值時(shí)間序列，并以3σ閾值原則識(shí)別粗差。具體過(guò)程為：利用SSA方法處理原靜態(tài)相對(duì)重力觀測(cè)時(shí)間序列，利用ω-correlation方法分析各重建成分之間的相關(guān)性，找出趨勢(shì)項(xiàng)、周期項(xiàng)以及第i階次，當(dāng)重建成分階次大于i時(shí)，各個(gè)重建成分之間不能很好地分離；合并前i階項(xiàng)得到主成分，計(jì)算原相對(duì)靜態(tài)重力觀測(cè)值序列和主成分序列之間的殘差并識(shí)別粗差。

本文通過(guò)奇異譜迭代方法插補(bǔ)粗差處的重力值，與文獻(xiàn)[1]不同之處在于：將殘差時(shí)間序列中粗差位置處的殘差值設(shè)置為0，然后通過(guò)SSA方法分解重構(gòu)重建殘差趨勢(shì)項(xiàng)，用趨勢(shì)項(xiàng)的值代替原缺失位置的0值。循環(huán)此過(guò)程，直至缺失位置處前后兩次插補(bǔ)的數(shù)據(jù)殘差RMS小于10-3μGal。這時(shí)，將重建的殘差數(shù)據(jù)序列與之前提取的主成分?jǐn)?shù)據(jù)序列相疊加，得到修正后的靜態(tài)相對(duì)重力時(shí)間觀測(cè)序列。

1.2 零漂和固體潮提取思路

針對(duì)長(zhǎng)時(shí)間靜態(tài)相對(duì)重力觀測(cè)數(shù)據(jù)中零漂和固體潮的變化特點(diǎn)和SSA方法分析有限長(zhǎng)度周期性信號(hào)的特點(diǎn)，本文提出采用滑動(dòng)SSA方法提取長(zhǎng)時(shí)間重力固體潮改正：

1)將原M天的靜態(tài)相對(duì)重力觀測(cè)值時(shí)間序列按L天為一個(gè)時(shí)間段進(jìn)行滑動(dòng)分組，滑動(dòng)窗口長(zhǎng)度為1 d，即任意相鄰兩組數(shù)據(jù)序列的起始時(shí)間相隔1 d，得到每個(gè)時(shí)間段M-L+1組數(shù)據(jù)序列。

2)對(duì)每個(gè)時(shí)間段內(nèi)每組靜態(tài)相對(duì)重力觀測(cè)值序列，利用多項(xiàng)式擬合得到L天內(nèi)的零漂值。從每組觀測(cè)值中減去對(duì)應(yīng)多項(xiàng)式擬合計(jì)算得到的零漂值，對(duì)其結(jié)果采用SSA方法提取固體潮，然后提取每組數(shù)據(jù)序列的固體潮改正。

3)對(duì)每個(gè)時(shí)間段內(nèi)M-L+1組數(shù)據(jù)均通過(guò)多項(xiàng)式擬合提取零漂，將剔除零漂后的重力數(shù)據(jù)觀測(cè)序列通過(guò)SSA方法提取重力固體潮改正值。重力固體潮(或零漂)按同一天取平均值獲得當(dāng)天重力固體潮改正(或零漂)的最終結(jié)果，即第1～M天分別在第1、2、…M-L+1組中找到對(duì)應(yīng)值，并按1 d、2 d、…、M天取平均。為降低SSA中邊界效應(yīng)的影響，在精度對(duì)比統(tǒng)計(jì)中不統(tǒng)計(jì)每組前2 d和最后2 d的數(shù)據(jù)。

2 實(shí)驗(yàn)與分析

2.1 實(shí)驗(yàn)數(shù)據(jù)及數(shù)據(jù)預(yù)處理

相對(duì)靜態(tài)重力觀測(cè)時(shí)間序列中除包含重力固體潮外，還包括海潮負(fù)荷、零漂和噪聲等。本次實(shí)驗(yàn)利用CG-5相對(duì)重力儀(序列號(hào)41221)在固定點(diǎn)(坐標(biāo)120.10°E、36.00°N；大地高25.5 m)進(jìn)行連續(xù)2個(gè)月(2016-03-05 00：00～05-04 23：59)的靜態(tài)相對(duì)重力觀測(cè)。相對(duì)重力儀采樣頻率為6 Hz，儀器前55 s進(jìn)行數(shù)據(jù)采集，后5 s進(jìn)行數(shù)據(jù)處理，因此1 min內(nèi)可得到330個(gè)數(shù)據(jù)。相對(duì)重力儀自帶地震濾波功能，可對(duì)獲取的靜態(tài)相對(duì)重力觀測(cè)數(shù)據(jù)進(jìn)行低頻噪聲過(guò)濾，舍棄高于6倍標(biāo)準(zhǔn)偏差的高頻噪聲。濾波之后對(duì)樣本數(shù)據(jù)每分鐘進(jìn)行一次平均處理并且作為結(jié)果自動(dòng)記錄，61 d共獲得87 840個(gè)數(shù)據(jù)。對(duì)相對(duì)重力觀測(cè)數(shù)據(jù)進(jìn)行預(yù)處理，采用本文提出的方法進(jìn)行重力固體潮改正提取。作為對(duì)比，采用多項(xiàng)式擬合提取零漂，將剔除零漂后的重力數(shù)據(jù)采用SSA提取重力固體潮改正，技術(shù)路線如圖1所示。

圖1 滑動(dòng)SSA提取重力固體潮技術(shù)路線

首先利用SSA方法進(jìn)行粗差探測(cè)，將4 d共5 760個(gè)重力觀測(cè)數(shù)據(jù)作為嵌入維數(shù)，對(duì)87 840個(gè)相對(duì)重力觀測(cè)數(shù)據(jù)時(shí)間序列進(jìn)行SSA分解和重構(gòu)，獲得RCs；通過(guò)ω-correlation方法分析各RC之間的相關(guān)系數(shù)，前30階重構(gòu)成分ω-correlation分析結(jié)果如圖2所示。根據(jù)相鄰兩個(gè)RC之間的相關(guān)系數(shù)大于0.85可以判斷，R2、R3、R4、R5、R7、R8、R9、R10、R16、R17、R20、R21、R26、R27、R28、R29均具有周期性。從重力數(shù)據(jù)序列SSA得到的前30階重構(gòu)部分FFT功率譜分析結(jié)果(圖3)可以看出，從第15階開始，不同頻率混疊較嚴(yán)重，而且振幅非常小，各重構(gòu)成分之間無(wú)法很好地分離，因此合并前14項(xiàng)重構(gòu)成分作為原重力數(shù)據(jù)時(shí)間序列主成分。原靜態(tài)相對(duì)重力觀測(cè)與主成分之間殘差時(shí)間序列如圖4(a)所示，根據(jù)3σ準(zhǔn)則選擇閾值區(qū)間為(-0.009 8 mGal，0.009 8 mGal)，若超出此閾值區(qū)間則判定為粗差，共624個(gè)。根據(jù)奇異譜迭代方法修正殘差處的粗差相對(duì)重力值，經(jīng)過(guò)迭代后得到粗差處相鄰兩次迭代修復(fù)值殘差的RMS為9.57×10-6μGal，修復(fù)后的殘差值如圖4(b)所示?？梢钥闯?，殘差序列較為平穩(wěn)，插補(bǔ)效果良好。修復(fù)后的殘差疊加主成分可得到修復(fù)后的靜態(tài)相對(duì)重力觀測(cè)時(shí)間序列，結(jié)果如圖5(a)所示。海潮和重力固體潮成因相同，因此具有相同的頻率。本文基于NAO.99b海潮模型[6]，結(jié)合TSoft軟件計(jì)算本次實(shí)驗(yàn)同歷元下的重力海潮負(fù)荷值時(shí)間序列，結(jié)果如圖5(b)所示，將海潮對(duì)相對(duì)重力的影響從原時(shí)間序列中去除，結(jié)果如圖5(c)所示。

圖2 靜態(tài)相對(duì)重力數(shù)據(jù)序列SSA前30階重建部分的ω-correlations

圖3 靜態(tài)相對(duì)重力數(shù)據(jù)序列SSA前30階重構(gòu)部分的FFT功率波譜分析

圖4 原靜態(tài)相對(duì)重力觀測(cè)時(shí)間序列與主成分之間的殘差

圖5 靜態(tài)相對(duì)重力觀測(cè)時(shí)間序列

2.2 滑動(dòng)SSA方法提取重力固體潮改正

為了驗(yàn)證滑動(dòng)SSA方法提取重力固體潮改正的效果，對(duì)61 d靜態(tài)相對(duì)重力觀測(cè)時(shí)間序列分別以10 d、15 d、20 d、25 d、30 d、40 d、50 d為時(shí)間段，并以1 d為滑動(dòng)窗口進(jìn)行分組，得到每個(gè)時(shí)間段內(nèi)包含的數(shù)據(jù)序列組數(shù)分別為52、47、42、37、32、22和12。

在提取重力固體潮之前，首先需要剔除每組靜態(tài)相對(duì)重力觀測(cè)數(shù)據(jù)中的零漂項(xiàng)。采用SSA方法識(shí)別CG-5相對(duì)重力觀測(cè)數(shù)據(jù)中的趨勢(shì)，將SSA窗口長(zhǎng)度設(shè)為4 d，對(duì)每組數(shù)據(jù)進(jìn)行分解和重構(gòu)，得到趨勢(shì)項(xiàng)RC1。然后采用多項(xiàng)式擬合提取零點(diǎn)漂移，擬合公式為：

y=ax2+bx+c

(1)

式中，a、b、c為擬合參數(shù)。利用SSA提取零漂時(shí)可剔除噪聲，因此效果會(huì)更好。對(duì)每組數(shù)據(jù)序列剔除擬合后的零漂值，然后選擇嵌入維數(shù)5 760進(jìn)行 SSA，對(duì)重構(gòu)成分進(jìn)行ω-correlation相關(guān)性分析。個(gè)別分組中兩個(gè)相鄰RC之間的相關(guān)系數(shù)大于0.80時(shí)無(wú)數(shù)據(jù)，因此對(duì)比實(shí)驗(yàn)中均選擇相關(guān)系數(shù)大于0.65(當(dāng)相關(guān)系數(shù)為0.75時(shí)結(jié)果類似，相關(guān)系數(shù)為0.65時(shí)提取結(jié)果精度提高更明顯)。從重構(gòu)項(xiàng)i≤15中得到每組中同頻率重構(gòu)序列，對(duì)重構(gòu)序列進(jìn)行合并。對(duì)每組重構(gòu)序列進(jìn)行求和，計(jì)算每組提取的重力固體潮改正。SSA在分解重構(gòu)過(guò)程中存在邊界效應(yīng)，因邊界效應(yīng)不是本文主要研究?jī)?nèi)容，為削弱其影響，將每組前2 d和最后2 d共5 760個(gè)重力固體潮提取結(jié)果舍去，再采用本文方法計(jì)算最終的重力固體潮改正值，得到的結(jié)果分別記為HSSA_10、HSSA_15、HSSA_20、HSSA_25、HSSA_30、HSSA_40和HSSA_50。

為了驗(yàn)證提取精度，將實(shí)測(cè)重力數(shù)據(jù)采用Eterna調(diào)和分析獲得潮波參數(shù)后用TSoft軟件計(jì)算得到重力固體潮值(記為ET_Data)和CG-5重力儀自帶軟件計(jì)算的重力固體潮值(記為CG5_Data)分別作為參考值，提取結(jié)果與同時(shí)刻ET_Data和CGS_Data之間的差值如圖6所示，差值統(tǒng)計(jì)如表1(單位μGal)所示。

圖6 不同方法提取的重力固體潮殘差對(duì)比

表1 滑動(dòng)SSA和SSA提取長(zhǎng)時(shí)間重力固體潮殘差統(tǒng)計(jì)

2.3 多項(xiàng)式擬合方法提取零漂和SSA提取重力固體潮

將61 d靜態(tài)相對(duì)重力觀測(cè)數(shù)據(jù)序列作為一個(gè)時(shí)間段長(zhǎng)度，采用多項(xiàng)式擬合方法提取零漂值，設(shè)嵌入維數(shù)為5 760，利用SSA分解重構(gòu)獲得RC1趨勢(shì)項(xiàng)，并利用式(1)進(jìn)行多項(xiàng)式擬合。將零漂數(shù)據(jù)從靜態(tài)相對(duì)重力序列中剔除，然后選擇嵌入維數(shù)5 760進(jìn)行SSA，對(duì)重構(gòu)成分進(jìn)行ω-correlation相關(guān)性分析，從結(jié)果中選出重構(gòu)項(xiàng)i≤15和相關(guān)系數(shù)大于0.65的兩個(gè)相鄰RC項(xiàng)，得到同頻率重構(gòu)序列，并對(duì)重構(gòu)序列進(jìn)行合并。對(duì)重構(gòu)序列進(jìn)行求和，計(jì)算提取的重力固體潮改正，得到的結(jié)果記為SSA_Data。為削弱邊界效應(yīng)影響，同樣將前2 d和最后2 d共5 760個(gè)重力固體潮提取結(jié)果舍去，其與同時(shí)刻ET_Data和CG5_Data的結(jié)果比較如圖6所示，差值統(tǒng)計(jì)如表1所示。

2.4 不同方法提取重力固體潮殘差分析

由表1可知，無(wú)論是滑動(dòng)SSA還是SSA方法，提取重力潮改正結(jié)果與理論值差值的標(biāo)準(zhǔn)差(STD)均小于11 μGal，說(shuō)明SSA方法提取重力固體潮改正結(jié)果較可靠，精度較高。對(duì)于同一種提取重力固體潮改正方法，選用CG5_Data作為理論參考值，其提取結(jié)果比選用ET_Data作為理論參考值時(shí)大，說(shuō)明SSA或滑動(dòng)SSA方法的提取結(jié)果更接近ET_Data理論值。采用滑動(dòng)SSA提取重力固體潮改正時(shí)，每個(gè)連續(xù)時(shí)間段長(zhǎng)度不同，其提取結(jié)果與理論值殘差的STD和RMS不同。選擇連續(xù)時(shí)間段長(zhǎng)度為30 d和40 d時(shí)，滑動(dòng)SSA提取重力固體潮改正值與理論值殘差的STD和RMS比SSA方法更小，結(jié)果更優(yōu)。選擇滑動(dòng)時(shí)間段長(zhǎng)度為10 d、15 d、20 d、25 d和50 d時(shí)，滑動(dòng)SSA方法提取重力固體潮改正值與理論值殘差的STD和RMS大于SSA方法，表明SSA方法更優(yōu)。選擇滑動(dòng)時(shí)間段長(zhǎng)度為30 d時(shí)，滑動(dòng)SSA方法提取重力固體潮改正值與ET_Data理論值殘差的STD和RMS比SSA方法均小0.43 μGal；滑動(dòng)SSA方法提取值與CG5_Data理論值殘差的STD和RMS比SSA方法分別小0.201 μGal和0.219 μGal。選擇滑動(dòng)時(shí)間段長(zhǎng)度為40 d時(shí)，滑動(dòng)SSA方法提取重力固體潮改正值與ET_Data理論值殘差的STD和RMS比SSA方法均小0.665 μGal；滑動(dòng)SSA方法提取值與CG5_Data理論值殘差的STD和RMS比SSA方法分別小0.614 μGal和0.602 μGal。選擇滑動(dòng)時(shí)間段長(zhǎng)度為10 d時(shí)，采用滑動(dòng)SSA方法提取重力固體潮時(shí)，提取結(jié)果與理論值差值的STD和RMS最大，即提取精度最低。殘差最大值與最小值之間的差值可表明殘差范圍，范圍越小說(shuō)明殘差的離散程度越小。當(dāng)選用CG5_Data作為理論參考值時(shí)，除選擇10 d進(jìn)行滑動(dòng)計(jì)算的結(jié)果外，滑動(dòng)SSA方法提取結(jié)果的離散程度均小于SSA方法，提取結(jié)果殘差離散程度最小的方法為選擇20 d進(jìn)行滑動(dòng)計(jì)算；當(dāng)選用ET_Data作為理論參考值時(shí)，除選擇10 d、50 d滑動(dòng)計(jì)算外，滑動(dòng)SSA方法提取結(jié)果的離散程度均小于SSA方法，提取結(jié)果殘差離散程度最小的方法為選擇20 d進(jìn)行滑動(dòng)計(jì)算。結(jié)果表明，通過(guò)滑動(dòng)SSA方法提取固體潮可降低提取結(jié)果殘差的離散程度。

由圖6可知，殘差序列均在零值附近擺動(dòng)，說(shuō)明提取的零漂與參考信號(hào)擬合程度較高。殘差時(shí)間序列仍然具有明顯的周期性，因?yàn)橄鄬?duì)重力序列中除零漂和固體潮外，還包括其他因素產(chǎn)生的周期信號(hào)，如天氣變化、水文變化以及其他干擾噪聲，這些因素在殘差序列中會(huì)有所反映。由于ET_Data理論值是實(shí)測(cè)重力數(shù)據(jù)采用Eterna調(diào)和分析獲得潮波參數(shù)，再利用TSoft軟件計(jì)算得到的，以ET_Data為參考值時(shí)比以CG5_Data為參考值時(shí)殘差振幅略小。選用不同連續(xù)時(shí)間段進(jìn)行滑動(dòng)計(jì)算的結(jié)果區(qū)別較大，其中選用連續(xù)時(shí)間段30 d、40 d計(jì)算結(jié)果較好。

3 結(jié) 語(yǔ)

本文將時(shí)間長(zhǎng)度為2個(gè)月的靜態(tài)相對(duì)重力觀測(cè)數(shù)據(jù)序列按時(shí)間窗口10 d、15 d、20 d、25 d、30 d、40 d和50 d分別分為52組、47組、42組、37組、32組、22組和12組，采用滑動(dòng)SSA方法進(jìn)行重力固體潮改正提取實(shí)驗(yàn)。結(jié)果表明，SSA方法提取重力固體潮結(jié)果精度較高，采用滑動(dòng)SSA和SSA方法提取的重力固體潮值與理論值(CG5_Data和ET_Data)殘差的標(biāo)準(zhǔn)差均小于11 μGal。與SSA方法相比，選用合適時(shí)間段長(zhǎng)度的滑動(dòng)SSA方法提取固體潮可以降低提取結(jié)果殘差的離散程度。選用CG5_Data值作為參考值、連續(xù)30 d和40 d數(shù)據(jù)作為一組時(shí)，滑動(dòng)SSA方法提取的重力固體潮改正結(jié)果殘差的RMS比未使用滑動(dòng)SSA方法時(shí)分別小0.219 μGal和0.602 μGal；選用ET _Data值作為參考值、連續(xù)30 d和40 d數(shù)據(jù)作為一組時(shí)，滑動(dòng)SSA方法提取的重力固體潮結(jié)果殘差的RMS比未使用滑動(dòng)SSA時(shí)分別小0.430 μGal和0.665 μGal。因此，選擇合適的連續(xù)滑動(dòng)時(shí)間段，滑動(dòng)SSA方法比SSA方法提取重力固體潮改正結(jié)果精度略高。