基于主成分分析的相對重力觀測零漂和固體潮提取

于紅娟 郭金運(yùn)，2 李九龍 穆大鵬 孔巧麗，2

1 山東科技大學(xué)測繪科學(xué)與工程學(xué)院，青島市前灣港路579號(hào)，266590

2 海島（礁）測繪技術(shù)國家測繪地理信息局重點(diǎn)實(shí)驗(yàn)室，青島市前灣港路579號(hào)，266590

3 中國科學(xué)院測量與地球物理研究所，武漢市徐東大街340號(hào)，430077

本文通過實(shí)驗(yàn)，利用主成分分析（principal component analysis，PCA）提取CG－5重力儀［1］靜態(tài)相對觀測的零點(diǎn)漂移和固體潮。將經(jīng)過固體潮改正的9d觀測數(shù)據(jù)作為原統(tǒng)計(jì)量，使用PCA 提取日零漂，并與最小二乘線性擬合值作對比分析，用含有固體潮和零漂影響的觀測數(shù)據(jù)，在進(jìn)行零點(diǎn)漂移改正的基礎(chǔ)上使用PCA 提取固體潮，并與CG－5重力儀內(nèi)部軟件提供的固體潮模型計(jì)算值作對比。通過兩組實(shí)驗(yàn)對比結(jié)果來驗(yàn)證本文方法的可靠性和適用性。

1 主成分分析方法

重力觀測的主成分分析［2－5］是基于協(xié)方差矩陣將重力數(shù)據(jù)作為原統(tǒng)計(jì)變量，通過線性變換得到新變量，新變量之間統(tǒng)計(jì)不相關(guān)。PCA 的主要目的是通過分離不相關(guān)成分提取主要信息并進(jìn)行數(shù)據(jù)降維。文中提取的零點(diǎn)漂移和固體潮為主要信息，即信號(hào)。主成分分析假設(shè)新變量PCAG＝（PCA1，PCA2，…，PCAn）T是由CG－5重力儀觀測數(shù)據(jù)的原統(tǒng)計(jì)變量G＝（G1，G2，…，Gn）T經(jīng)過線性變換得到的，即

式中，E是由重力觀測的統(tǒng)計(jì)變量G的協(xié)方差矩陣計(jì)算得到的線性變換矩陣。假設(shè)協(xié)方差矩陣為ΣG，則它是一個(gè)n階對稱陣。通過求解特征根和特征向量，矩陣ΣG可以分解為：

式中，D是對角陣，對角線上的元素λ1≥λ2…≥λn是矩陣ΣG的特征根，同時(shí)等于對應(yīng)重力觀測數(shù)據(jù)主成分的方差，E是由特征向量構(gòu)成的正交矩陣［3］。因此，可由式（1）得到：

前k個(gè)主成分的綜合能力以及所包含的重力值信息量可用累積貢獻(xiàn)率α表示：

式中，λi表示第i個(gè)特征值，p表示特征值總個(gè)數(shù)，k表示選取的個(gè)數(shù)。

如果選取較大的k個(gè)特征值所對應(yīng)的新變量重構(gòu)原重力觀測值統(tǒng)計(jì)變量，則可以達(dá)到提取主要信息的目的［4－5］：

式中，ET表示由選擇的特征向量構(gòu)成的k階矩陣。

2 數(shù)據(jù)及預(yù)處理

靜態(tài)相對重力數(shù)據(jù)由CG－5重力儀（序列號(hào)140541221）于某一固定觀測站觀測9d（2014－08－15～23）所得，以每天09：21至次日09：20為一個(gè)時(shí)間序列，即作為一個(gè)統(tǒng)計(jì)觀測量。儀器自身具有地震濾波功能，能對低頻噪聲進(jìn)行過濾，并舍棄高于6倍標(biāo)準(zhǔn)偏差的高頻噪聲。采用6 Hz采樣率，可對1min的360個(gè)樣本數(shù)據(jù)進(jìn)行平均得到一個(gè)讀數(shù)，觀測過程中儀器進(jìn)行了溫度補(bǔ)償和傾斜改正。因此，每天可以采集同一個(gè)時(shí)間段（24 h）的1 440個(gè)重力觀測數(shù)據(jù)。

把所關(guān)心的數(shù)據(jù)成分認(rèn)為是信號(hào)，其他的作為噪聲。假設(shè)信號(hào)與噪聲統(tǒng)計(jì)不相關(guān)，嘗試?yán)肞CA 把這9d的時(shí)間序列作為9個(gè)統(tǒng)計(jì)量來分離信號(hào)和噪聲。實(shí)驗(yàn)使用了兩組觀測數(shù)據(jù)，第一組數(shù)據(jù)是進(jìn)行了固體潮改正的靜態(tài)相對觀測數(shù)據(jù)，嘗試運(yùn)用PCA 分離信號(hào)和噪聲來估算CG－5 重力儀的靜態(tài)零漂；第二組數(shù)據(jù)是未進(jìn)行固體潮改正的靜態(tài)觀測數(shù)據(jù)，嘗試使用PCA 分離信號(hào)和噪聲來提取固體潮。

由于觀測站環(huán)境干擾較多，采用滑動(dòng)平均的方法［6］對數(shù)據(jù)進(jìn)行預(yù)處理。通過實(shí)驗(yàn)并根據(jù)重力數(shù)據(jù)的實(shí)際變化情況，選取大小為5的移動(dòng)窗口進(jìn)行平滑處理。

3 靜態(tài)零漂

用第一組數(shù)據(jù)進(jìn)行實(shí)驗(yàn)，把9d觀測數(shù)據(jù)作為原統(tǒng)計(jì)量，使用PCA 將噪聲和信號(hào)進(jìn)行分離。每天的數(shù)據(jù)作為一個(gè)統(tǒng)計(jì)量，每個(gè)統(tǒng)計(jì)量為含有1 440×1個(gè)元素的列向量，這樣9d的觀測值就組成了1 440×9 的矩陣，即，。在得到統(tǒng)計(jì)量G后，通過式（2）計(jì)算其特征向量構(gòu)成的線性變換矩陣E，再由式（1）得到主成分的9個(gè)模式。圖1給出了所得到的9個(gè)特征根的值，它們等于對應(yīng)重力觀測數(shù)據(jù)主成分的方差。特征值相對很小甚至接近于0則認(rèn)為不含有主要的信息，因此，可以從降序排列的特征根中，找到一個(gè)特征根來確定主成分的個(gè)數(shù)，使得剩下的特征根具有接近相同且相對很小的值［4］。從圖1可看出，除第一個(gè)外其他的特征值都非常小，接近于0。并且通過計(jì)算得到第一個(gè)模式的貢獻(xiàn)率為99.8%，即第一模式包含了幾乎所有的數(shù)據(jù)信息。從圖1分析得出，第一模式可以認(rèn)為是信號(hào)部分，而其他8個(gè)模式為噪聲部分。

圖1 零漂提取中9個(gè)模式對應(yīng)的特征值Fig.1 Eigenvalues corresponding to nine modes in the zero drift estimation

因此，選擇第一模式PCA1重構(gòu)統(tǒng)計(jì)量G，分別得到由PCA 估算的9d 的零點(diǎn)漂移。由式（5），選擇第一模式進(jìn)行重構(gòu)得：

將其展開，得：

表1 08－15～23 零點(diǎn)漂移計(jì)算結(jié)果Tab.1 Calculated results of daily zero drifts from 08－15to 08－23

從表1可以看出，PCA 與LSLF 的9d日零漂差值均在10－2μGal／d量級(jí)，具有非常好的一致性，但是PCA 剔除了噪聲，效果更好。因此可以說，PCA 更適用于CG－5靜態(tài)數(shù)據(jù)的日漂移常數(shù)估算。08－15～23 的日漂移常數(shù)有明顯增加，說明CG－5 重力儀的零漂具有隨時(shí)間增加的特點(diǎn)。

4 固體潮

用PCA 得到的日零漂常數(shù)對第二組數(shù)據(jù)進(jìn)行線性改正，得到經(jīng)零漂改正但未進(jìn)行固體潮改正的9d數(shù)據(jù)，將其作為原統(tǒng)計(jì)量并利用主成分分析分離信號(hào)和噪聲來估算固體潮。每個(gè)統(tǒng)計(jì)量為含有1 440×1個(gè)元素的列向量，這樣9d的觀測值就組成了1 440×9的矩陣，即

得到統(tǒng)計(jì)量G之后，通過式（2）計(jì)算其特征向量構(gòu)成的線性變換矩陣E，再由式（1）得到主成分的9個(gè)模式。研究可得，PCA 實(shí)現(xiàn)了信號(hào)和噪聲的分離，圖2給出了所得到的9個(gè)特征根的值，它們分別等于對應(yīng)重力數(shù)據(jù)主成分的方差。

圖2 固體潮提取中9個(gè)模式對應(yīng)的特征值Fig.2 Eigenvalues corresponding to nine modes in the solid earth tide estimation

從圖2可知，前3個(gè)特征值為主要成分［4］。進(jìn)一步計(jì)算得到，第一模式的貢獻(xiàn)率為65.2%，前兩個(gè)模式的累積貢獻(xiàn)率為94.8%，前3個(gè)模式的累積貢獻(xiàn)率為99.9%，幾乎包含了重力數(shù)據(jù)的全部信息。因此，主成分的前3個(gè)模式為分離得到的信號(hào)，選擇該模式來重構(gòu)統(tǒng)計(jì)量G，分別得到PCA 對應(yīng)的每一個(gè)模式的值。圖3 給出了線性改正后的9d數(shù)據(jù)時(shí)序圖。從圖3可見，15和16日的波形變化相近，17～19日的波形變化相近，20～23日的波形變化相近。又因?yàn)槭芄腆w潮影響，重力觀測曲線每天可觀測到兩個(gè)波峰、兩個(gè)波谷［7］。因此，初步推測波形的不同是這9d的固體潮影響不同所致，并推測這3個(gè)模式對應(yīng)固體潮的變化。通過進(jìn)一步計(jì)算PCA 的前3個(gè)模式累積量（見式（9）），即對前3個(gè)模式疊加計(jì)算并進(jìn)行中心化處理后的結(jié)果（PCA）與CG－5重力儀內(nèi)部軟件提供的固體潮模型計(jì)算的固體潮值（solid earth tide，SET）比較發(fā)現(xiàn)，3 個(gè)模式正好對應(yīng)于固體潮的影響量。

圖3 線性漂移改正后的重力數(shù)據(jù)時(shí)序圖Fig.3 Gravimetric datatime series of 9days after the linear correction

于是，由式（5）選擇前3個(gè)模式進(jìn)行重構(gòu)：

將其展開，得：

從表2看出，PCA和SET計(jì)算的固體潮值統(tǒng)計(jì)結(jié)果中，兩者最大差值僅5μGal，具有較好的一致性，說明PCA 可用于靜態(tài)觀測中的固體潮估算。15～23日固體潮的平均值均在0 附近，SET 計(jì)算的均值稍微偏離零軸，但整體上正負(fù)值具有較好的對稱性；15～18日的固體潮高潮值達(dá)到90μGal，低潮值達(dá)－80μGal；19～23日的固體潮高潮值達(dá)124μGal，低潮值達(dá)－78μGal。結(jié)合圖3、圖4 可以看出，固體潮具有一定的規(guī)律性，每天都會(huì)有兩個(gè)高潮和低潮，并且高潮和低潮會(huì)隨日期增加往后推遲一段時(shí)間。PCA 和SET的均方根和標(biāo)準(zhǔn)差較為接近，相差僅在10－1μGal量級(jí)。PCA 和SET 的差值最大均不大于8 μGal，最小都在5μGal內(nèi)；均值在0附近波動(dòng)，波動(dòng)范圍在3μGal內(nèi)；均方根小于5μGal，說明PCA 與SET 計(jì)算結(jié)果具有很好的一致性。

采用以上計(jì)算結(jié)果對觀測數(shù)據(jù)進(jìn)行線性改正和固體潮改正（圖5）。CG－5 相對重力儀的重復(fù)性標(biāo)準(zhǔn)差為5μGal［1］，而從圖5中可以看出，經(jīng)過改正后的15～23日重力數(shù)據(jù)變化很小，均在5 μGal范圍內(nèi)。另外，導(dǎo)致這種變化的原因可能還有海潮、零漂的非線性等影響，最后可對數(shù)據(jù)進(jìn)行平均得到該測站點(diǎn)的重力讀數(shù)值。

圖4 固體潮于08－15～23的PCA 與SET 時(shí)序圖Fig.4 Time series of the solid Earth tide derived from PCA and SET from 08－15to 08－23

表2 固體潮于08－15～23的PCA與SET統(tǒng)計(jì)結(jié)果Tab.2 Statistical results of the solid earth tide derived from PCA and SET from 08－15to 08－23

圖5 固體潮改正后的重力觀測值時(shí)序圖Fig.5 Gravimetric observations time series after zero drift and solid earth tide correction

5 結(jié) 語

1）PCA 提取線性漂移時(shí)，第一個(gè)模式的貢獻(xiàn)率為99.8%，該模式包含了幾乎所有數(shù)據(jù)信息。因此，用該模式估算的線性漂移與LSLF值對比，相差僅在10－2μGal量級(jí)，但PCA 剔除了噪聲，效果更好。說明在估算零點(diǎn)漂移方面，PCA 比LSLF更具有優(yōu)越性。

2）PCA 提取固體潮時(shí)，第一模式的貢獻(xiàn)率為65.2%，前兩個(gè)模式的累積貢獻(xiàn)率為94.8%，前3個(gè)模式的累積貢獻(xiàn)率達(dá)99.9%。可以看出，前3個(gè)模式幾乎包含了重力數(shù)據(jù)的全部信息。前3個(gè)模式的疊加計(jì)算正好對應(yīng)固體潮，并與SET值基本一致。PCA－SET 的統(tǒng)計(jì)值均小于8 μGal，其中均方根都小于5μGal。因此，PCA 用于估算固體潮具有可靠性。

3）根據(jù)以上計(jì)算結(jié)果對觀測數(shù)據(jù)作零漂和固體潮改正后，其變化均在5μGal范圍內(nèi)，導(dǎo)致這種變化的原因可能是儀器的測量精度、海潮和零漂的二次項(xiàng)等微小影響。

［1］Scintrex Limited.CG－5Scintrex Autograv System Operation Manual V5.0［Z］.2009

［2］Wold S，Esbensen K，Geladi P.Principal Component Analysis［J］.Chemometrics and Intelligent Laboratory Systems，1987，2（1）：37－52

［3］穆大鵬.基于主成分分析的GRACE 重力場模型等效水高［J］.地球物理學(xué)進(jìn)展，2014，29（4）：1 512－1 517（Mu Dapeng.Equivalent Water Height from GRACE Gravity Model Based on Principal Component Analysis［J］.Progress in Geophysics，2014，29（4）：1 512－1 517）

［4］Johnson R A，Wichern D W.Applied Multivarariate Statistical Analysis［M］.Englewood Cliffs N J：Prentice Hall，1992

［5］Jolliffe I T.Principal Component Analysis［M］.New York：Springer，2012

［6］史文海，李正農(nóng)，吳建佳.近地面強(qiáng)風(fēng)不同間隔滑動(dòng)平均統(tǒng)計(jì)特性的對比分析［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2013，31（5）：611－615（Shi Wenhai，Li Zhengnong，Wu Jianjia.A Contrastive Study of Different Time Interval Moving Average Statistical Characteristics of Strong Wind Near Ground［J］.Acta Aerodynamica Sinica，2013，31（5）：611－615）

［7］顧兆峰.KSS 31M 海洋重力儀靜態(tài)觀測結(jié)果及分析［J］.海洋測繪，2005，25（2）：66－68（Gu Zhaofeng.The Static Measurement Result of KSS 31M Marine Gravimeter and Its Analysis［J］.Hydrographic Surveying and Charting，2005，25（2）：66－68）