基于成分聚類的高階奇異譜分析及在GNSS監(jiān)測序列分析中的應用

翟長治，岳　順，李小奇

(1.武漢大學 測繪學院,湖北 武漢 430079;2.河海大學 地球科學與工程學院,江蘇 南京 210098；3.上海勘測設計研究院有限公司，上海 200434)

?

基于成分聚類的高階奇異譜分析及在GNSS監(jiān)測序列分析中的應用

翟長治1，岳順2，李小奇3

(1.武漢大學 測繪學院,湖北 武漢 430079;2.河海大學 地球科學與工程學院,江蘇 南京 210098；3.上?？睖y設計研究院有限公司，上海 200434)

摘要：高階奇異譜分析(HSSA)相對于奇異譜分析對不同延時、嵌入維數(shù)變化有較好的魯棒性，在信號處理、工程等領域有著廣泛應用。針對高階奇異譜分析應用中成分序列的選取憑經(jīng)驗過于主觀的問題，在理論推導與分析的基礎上，提出一種基于成分功率譜聚類的方法，并應用于實際GNSS監(jiān)測序列的處理與分析，得到了具有明顯以一天為周期的周期成分、趨勢變化成分以及不規(guī)則變動成分，反映了不同因素造成的大橋索塔變形，為后續(xù)索塔變形的建模預報、安全評估提供了技術支撐，也說明方法的實用性。

關鍵詞：高階奇異譜;成分聚類;功率譜;GNSS監(jiān)測序列

隨著計算機技術、對地觀測技術，以及各種傳感器技術的發(fā)展，自動化的、實時的監(jiān)測技術己經(jīng)成為了可能。其中GNSS系統(tǒng)以精度高、全天候、實時性等優(yōu)點成為這些網(wǎng)絡的重要組成部分。這些GNSS系統(tǒng)在提供其常規(guī)服務的同時，也積累了大量包含區(qū)域或建筑物形變等重要信息的監(jiān)測序列。但由于環(huán)境因素的復雜性以及各種隨機干擾的影響，監(jiān)測序列表現(xiàn)為各種變形與誤差的疊加，因此如何有效分離與提取變形信息，為后續(xù)變形分析、預報、安全評估提供技術支持，具有重要的實用價值[1]。常用的分離與提取變形信息的方法有小波變換、現(xiàn)代譜估計、卡爾曼濾波等方法。近年因奇異譜分析具有穩(wěn)定識別和強化周期信號的優(yōu)良特性，國內外學者也嘗試將其應用于變形序列的分析與變形特征提取[2]。但應用中嵌入維數(shù)、延時以及有效成分的確定存在主觀性，在一定程度上限制了其應用的效果。高階奇異譜分析(HSSA)是一種高階統(tǒng)計方法，與奇異譜分析相比，具有對不同嵌入維數(shù)和延時的魯棒性，并有盲高斯噪聲性和抑制相關噪聲的特性[3]。但在高階奇異譜分析應用中同樣存在重構成分的選取問題。以往研究者主要依據(jù)經(jīng)驗進行分類、選擇，存在較強的主觀性，從而會影響重構結果的客觀性。針對這一問題，本文提出一種根據(jù)成分序列功率譜分析結果聚類的重構成分選擇法，并對實際的GNSS觀測數(shù)據(jù)進行分析處理。

1高階奇異譜分析的原理

高階奇異譜是一種基于高階累積量的非線性信號的處理方法。高階奇異譜是將原狀態(tài)向量的高階累積量進行切片處理，得到二元的高階矩函數(shù)。然后將二元的高階矩函數(shù)進行SVD分解，得到奇異值和左右奇異向量。高階奇異譜重構是將原狀態(tài)向量向奇異向量空間進行投影，分解出與奇異值數(shù)量相同的成分序列，然后從中選擇部分成分序列進行重構，從而獲得感興趣的有用信息。令xN表示某一時段的觀測數(shù)據(jù)，序列長度為N。令m1,m2,m3,m4表示序列的一、二、三、四階矩。

(1)

式中：t1,t2,t3表示時延值;E(·)表示數(shù)學期望。四階累積量為

Cum4(t1,t2,t3)=

(2)

由于四階譜的計算量較大，而時延值t1,t2,t3可以任意給定，為便于進行奇異值分解處理，可以對原高階矩函數(shù)進行切片處理，得到二元高階矩函數(shù)。如令t2=t3，則

Cum4(t1,t2,t3)=

(3)

令Tx(i,j)=Cum4(i,i,j)，便得到高階矩函數(shù)的切片矩陣

Tx=

對Tx進行奇異值分解，求得Tx的奇異值V和奇異向量Ek。將時間序列x1,x2,x3,…,xN中心化處理之后，排列得到共N-M+1個狀態(tài)，N為序列長度，M為嵌入維數(shù)，其中第i個狀態(tài)量為

(5)

矩陣中的元素與原來的時間序列的對應關系為Xji=xj+i。把式(1)看做M個變量N-M+1個樣本。定義原狀態(tài)向量在奇異向量空間的投影

(6)

也稱為時間系數(shù)?？梢杂善娈愊蛄亢蜁r間系數(shù)來重建x的成分，得到M個x的成分序列。

(7)

在實際應用中，可以根據(jù)需要從M個成分序列中選擇部分序列進行重構，得到具有不同性質的重構序列[4]。

(8)

其中每個重構序列反映了原狀態(tài)向量的不同性質。

2基于成分聚類的重構成分選取

應用奇異譜分析進行序列的特征提取最關鍵的是選擇有效的成分序列進行重構。不同成分序列具有不同的特征，代表了原始序列中不同的影響，對于變形監(jiān)測序列來說，就是代表不同因素引起的不同形變特征[5]。因此針對不同的應用，應選擇合適的成分，選擇不當，會使重構序列內容混亂，不能有效提取原序列中不同方面的影響，甚至有可能丟失有用信息，引進不必要的噪聲[6-8]。針對有效成分的選取，以往的研究大都根據(jù)經(jīng)驗確定，具有很大的主觀性。根據(jù)高階奇異譜分析的原理，在分解時奇異值按從大到小排列，因此基于每個奇異值重構的成分序列代表原始序列中不同特征的信息，即具有不同的頻率與振幅。顯然可以根據(jù)成分序列的功率譜來分析成分序列，基于功率譜來選擇重構成分。對于GNSS監(jiān)測序列，由于環(huán)境影響因素以及隨機干擾多而復雜[9-10]，在實際工作中不可能分離每種影響因素引起的變形大小，事實上環(huán)境因素也在不斷變化中，表現(xiàn)在頻率上為一個頻率區(qū)間，而非一個固定的頻率。故此，提出基于成分聚類的重構成分選取，即根據(jù)成分頻率對成分進行聚類分析，然后根據(jù)聚類結果選取成分進行重構，達到分離提取不同特征信號的目的。

首先根據(jù)高階奇異譜分析的原理計算得到每一個成分序列xk,1≤k≤M，采用自相關函數(shù)估計法得到各成分的主頻率fi。若F看作是頻率對于成分序列號的函數(shù)

(10)

由于H(i)并不是嚴格單調遞增，i>I時，f(i)仍有可能小于H1，這時可以結合奇異值對成分序列進行取舍?？紤]到高階奇異譜下降比較迅速，當i>I時，其所對應的奇異值一般都非常小，所以可以將其舍去，即只取前i個作為聚類總體。將成分序列截取得到的聚類總體{fi|0

1)計算樣品兩兩距離，得一距離陣記為D0，開始每個樣品自成一類；

2)找出D0的非對角線最小元素，設為Dpq,則將Dp和Dq合并成一個新類，記為Dr；

3)給出計算新類與其它類的距離公式：Dkr={Dkp,Dkq}，將D(0)中第p,q行及p,q列用上面公式并成一個新行新列，新行新列對應Gr，所得到的矩陣記為D(1)。

4)對D(1)重復上述對D0的2)、3)兩步得D(2)；如此下去，直到所有的元素并成一類為止。

最后根據(jù)聚類結果，把屬于同一類的成分序列進行重構，得到幾個不同的重構序列。

3蘇通大橋GPS觀測數(shù)據(jù)的高階奇異譜分析

本文采用的實際觀測數(shù)據(jù)為蘇通大橋索塔頂部的GPS觀測數(shù)據(jù)，GPS系統(tǒng)使用的是徠卡GRX1200Pro系統(tǒng)，采樣間隔為30 s。系統(tǒng)共5臺GPS接收機，其中4臺為監(jiān)測站，1臺為基準站[11]。本文分析了X方向(橋軸方向)的數(shù)據(jù)，而Y,Z方向的變化性質與X方向相似[12]，可以以此類推。選取其中30 d的觀測數(shù)據(jù)，為了減小數(shù)據(jù)量，采用每10 min平滑得到一個數(shù)據(jù)，最后得到分析數(shù)據(jù)序列長度為4 320，同時為了分析方便，對原數(shù)據(jù)進行了中心化處理。利用基于成分聚類的高階奇異譜方法提取大橋索塔變形特征。

按照上述方法對數(shù)據(jù)進行高階奇異譜分析。序列長度N=4 320，嵌入維數(shù)M=1 500。得到1 500個成分序列。對各成分序列進行功率譜分析，得到各個成分序列的主頻率值。

圖1顯示，第50序列以前，各成分頻率變化平穩(wěn)，主要頻率在0.05 Hz以下，而從第50序列以后，頻率迅速提高，而且變化無規(guī)則。對比奇異值歸一化的結果(見圖2)，第50序列后的奇異值趨于0。因此本文選取前50個成分序列進行聚類分析。

對前50個成分序列進行最短距離聚類分析，得到聚類冰柱狀圖，如圖3所示。

根據(jù)索塔變形特性和有關資料，并結合聚類分析圖，把前50個成分可以分成三類較為合適，具體分類見表1。按分類結果進行重構，分別見圖4～圖6。

表1　最短距離分類

圖1　前100個成分序列的主頻率

圖2　歸一化奇異值圖

圖3　主頻率聚類分析結果

圖4　周期序列

從圖4中可以看出，點位變化在1 d內具有明顯的周期特性，稱為日周期變形，日周期變形幅值為0.08 m。由于日照方向的影響，陽面的混凝土溫度比陰面混凝土溫度高，該溫差將導致陽面的混凝土膨脹，而陰面的混凝土收縮，使塔柱頂部產(chǎn)生無外力作用下的扭轉變形，隨著時間推移，日照的方向在周期性地變化，因此這種扭轉變形在一日24 h內呈周期性變化。另外，由于各天天氣狀況不同，溫度高低有別，使各周期變形的幅值不同。

圖5　趨勢項序列

圖6　不規(guī)則變形序列

圖5顯示了索塔一個月內的趨勢變形，波動范圍為[-0.15,0.05]m，總體具有較好的周期性。

圖6為不規(guī)則變形部分，可見波動范圍較小，且無固定周期性，可能是由于受兩側不對稱索力、風力和橋面荷載所引起的索塔不規(guī)則變形。將三部分合成即反映了一個月內索塔的總體變形，得到圖7。

圖7　原序列與高階奇異譜分析后重建序列

從圖7可見，在一個月內大橋索塔在X方向的趨勢項波動范圍為[-0.2,0.2]m,整體上處于周期變化狀態(tài)，說明索塔具有較好的周期性。原GPS監(jiān)測序列表面粗糙，有奇異點，對監(jiān)測分析帶來了較大的影響，通過高階奇異譜分析得到的曲線連續(xù)光滑，說明該方法具有顯著的濾波效果，對原GPS序列的噪聲進行了有效地去除，可以有效地提取索塔的變形特征，摒棄干擾信息。

4結論

本文針對高階奇異譜分析應用中重構成分的選取問題，提出了一種基于成分功率譜聚類的重構成分選擇方法，克服了以往評經(jīng)驗、主觀性強的缺點。通過對蘇通大橋索塔30 d GNSS實測數(shù)據(jù)處理與分析，提取除了索塔3種變形：周日變形、趨勢變形與不規(guī)則變形。由于日照變化造成的日周期變化、天氣變化造成的總體變形趨勢和風、車輛荷載等因素造成的不規(guī)則變形。分析結果與以往大橋索塔變形分析結果基本一致，并在時域和空間域上更加細致，顯示了該方法的有效性。重構成分聚類在一定程度上基于經(jīng)驗，仍缺乏成熟的理論支持，有待進一步深入研究。

參考文獻：

[1]岳建平,方露,黎昵.變形監(jiān)測理論與技術研究進展[J].測繪通報,2007(7):1-4.

[2]王解先,連麗珍.沈云中.奇異譜分析在GPS 站坐標監(jiān)測序列分析中的應用[J].同濟大學學報(自然科學版),2013,41(2):282-288.

[3]袁堅,肖先賜.非線性時間序列的高階奇異譜分析[J].物理學報,1998,47(6):18-26.

[4]楊向鋒,張效民,李亞安.一種基于高階統(tǒng)計量的相空間重構方法[J].信號處理,2006,22(1):5-8.

[5]岳東杰.水利水電工程變形監(jiān)測中GPS技術與數(shù)據(jù)處理研究[D].南京：河海大學,2006.

[6]黃大偉,戴吾蛟,文斌.相空間重構單通道ICA算法及其在變形監(jiān)測數(shù)據(jù)去噪中的應用[J].測繪工程,2014,23(6):38-42.

[7]李征航,張小紅,朱智勤.利用GPS進行高精度變形監(jiān)測的新模型[J].測繪學報,2002,31(3):206-210.

[8]吳江淮.建筑物靜態(tài)變形監(jiān)測數(shù)據(jù)分析與處理[D].上海：同濟大學,2007.

[9]劉成龍,張德強,黃澤純.大跨懸索橋施工索塔變形成因分析與監(jiān)測[J].西南交通大學學報,2000,35(5):501-504.

[10] 王化光,黃育龍.鄂東長江公路大橋索塔日照變形監(jiān)測及應用[J].鐵道建筑,2010(3):14-17.

[11] 余波,邱洪興,王浩，等.蘇通大橋結構健康監(jiān)測系統(tǒng)設計[J].地震工程與工程振動,2009,29(4):170-177.

[12] 翟長治,岳東杰,岳順.大型橋梁GPS監(jiān)測數(shù)據(jù)的奇異譜分析[J].勘察科學技術,2014(4):20-23.

[責任編輯：劉文霞]

Higher order singular spectrum analysis based on component cluster and its applications of analysis of GNSS monitoring sequence

ZHAI Changzhi1,YUE Shun2,LI Xiaoqi3

(1.School of Geodesy and Geomatics,Wuhan University,Wuhan 430079,China;2.School of Earth Science and Engineering,Hohai University,Nanjing 210098,China;3.Shanghai Investigation,Design & Research Institnte Co.,Ltd.,Shanghai,200434,China)

Abstract：Higher Singular Spectrum Analysis (HSSA) has better robutness than singular spectrum analysis in the different delay,and embedding dimension changes.HSSA has been used in wide range of applications such as signal processing,engineering and other fields.For solving the problem of HSSA component sequences selected too subjectively,a new method is proposed based on component of the power spectral clustering,which is applied to the processing and analysis of the actual sequence of GNSS monitoring.The result has been obvious to one day for periodic components,trends and changes in composition of irregular fluctuation component,reflecting the bridge tower deformation caused by different factors.These works provide technical support for the following tower deformation modeling forecasts,safety assessment and also approve the applicability of this method.

Key words：higher-order singular spectrum;component cluster;power spectrum;GNSS monitoring sequence

中圖分類號：O433.4

文獻標識碼：A

文章編號：1006-7949(2016)04-0046-05

作者簡介：翟長治(1991-)，男，碩士研究生.

收稿日期：2015-01-15