于紅娟 郭金運，2 李九龍 穆大鵬 孔巧麗，2

1 山東科技大學測繪科學與工程學院，青島市前灣港路579號，266590

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

3 中國科學院測量與地球物理研究所，武漢市徐東大街340號，430077

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

1 主成分分析方法

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

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

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

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

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

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

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

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

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

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

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

3 靜態(tài)零漂

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

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

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

將其展開，得：

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

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

4 固體潮

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

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

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

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

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

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

將其展開，得：

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

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

圖4 固體潮于08－15～23的PCA 與SET 時序圖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)計結(jié)果Tab.2 Statistical results of the solid earth tide derived from PCA and SET from 08－15to 08－23

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

5 結(jié) 語

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

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

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

［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］.地球物理學進展，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)，吳建佳.近地面強風不同間隔滑動平均統(tǒng)計特性的對比分析［J］.空氣動力學學報，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）