大壩變形時間序列的奇異譜分析

李世友，王奉偉，沈云中

(同濟大學測繪與地理信息學院，上海 200092)

大壩變形受時效、溫度和水位等確定性因素[1-2]和許多非確定性因素的影響，其變形序列中既有低頻信號也有高頻信號，因此對大壩變形的分析和預測通常十分困難。目前，大壩變形預測的方法主要有灰色模型[3-4]、時間序列分析(ARMA)[5-6]、多元線性回歸[7]和神經網絡模型[8-9]，這些方法純粹從數學或機器學習的角度出發分析變形時間序列。王新洲等[10]用WT-SVM把變形時間序列分解成具有不同頻率特征的分量，分別預測各分量并進行重構，但沒有分析大壩變形的各變形因子。奇異譜分析法(singular spectrum analysis，SSA)是一種與經驗正交函數相聯系的統計技術，具備從含有強背景噪聲的時間序列中有效提取信號的能力，特別適合研究周期振蕩行為的分析方法。本文擬利用SSA分析大壩變形監測時間序列，提取趨勢項和周期項，分析其變形影響因子，并建模預測大壩的變形趨勢。

1 奇異譜分析法

奇異譜分析法由Broomhead等[11]提出并用于處理氣象數據，該方法不需要事先假定信號類型，利用經驗正交函數將時間序列數據分解成以信號為主的子空間和以噪聲為主的子空間，利用信號為主的子空間數據分析序列的時間演變特征。

SSA[12-14]根據一維時間序列X=(x1,x2,…,xN)，構建軌跡矩陣D如下

(1)

式中，M為窗口長度。對于非平穩序列，根據軌跡矩陣計算協方差陣C如下

C=DDT

(2)

對協方差陣C進行特征值分解，確定特征值λk及對應的特征向量Vk(由大到小排列)。第k個特征值對應的主成分為

(3)

(4)

利用所有M個主成分可完全恢復原始序列X=(x1,x2,…,xN)。如果前r個主成分的貢獻統計上屬于信號子空間，則利用這些主成分重構時間序列的信號部分如下

(5)

其余部分屬于噪聲子空間，可用于重構時間序列的噪聲。

2 大壩變形的奇異譜分析

試驗采用我國某大壩2004年5月1日至2010年8月31日單點的徑向形變GPS監測數據(向河流上游為負，下游為正)，同時采集了大壩的水位數據和壩區的溫度數據，采樣間隔為1 d，數據長度N為2314。因為采集的數據序列存在一定的缺失或粗差，試驗前要對數據進行預處理。首先對各項數據進行粗差探測并剔除。對于數據缺失插補問題有很多解決方法，如Shen等[15]在Schoellhamer[16]研究的基礎上提出了一種改進的奇異譜分析插補方案(ISSA)。本文對于缺失的數據序列采用三次樣條插值方法進行插補。

窗口長度M對于提取和分析信號具有很大的影響，選取時通常會考慮周期分量的識別效果，一般取周期的整數倍且滿足N/3≤M≤N/2[17]。通過試驗對比，選取窗口長度M為730(約2年)，以便能較好地提取出各時間序列中的信號成分。用SSA分析大壩變形時間序列(如圖1所示)，前50個特征值曲線如圖2所示，提取的信號分量如圖3所示。由圖2和圖3分析可知，第一個分量是趨勢項，能量貢獻占信號總能量的92.98%；第二和第三個分量為一對周期分量，周期為1.00年，能量貢獻率和為6.29%；第四和第五個分量也是一對周期分量，周期為0.78年，能量貢獻率和為0.11%。趨勢和周期分量的能量和占信號總能量的99.38%，余量占比0.62%即可歸為噪聲部分。為了進一步分析各變形分量與溫度和水位因子的關系，利用SSA分別提取溫度和水位數據中的周期分量，具體如圖4、圖5所示。

圖1 大壩變形時間序列

圖2 特征值曲線(前50個)

圖3 SSA提取趨勢和周期分量及其貢獻率

圖4 溫度數據及其周期分量

圖5 水位數據及其周期分量

從SSA分析結果可以發現：大壩由于長期的外力荷載存在徐變，表現為趨勢分量，主要與時效有關。從大壩變形時間序列中提取出了周期為1.00年的一個分量，同時從溫度和水位數據中也提取出了周期分別為1.00年和0.99年的分量。由此可知，大壩存在約1年的周期彈性形變，且主要與溫度和水位的周期變化相關。比較大壩變形、溫度和水位的周期分量(如圖6所示)可以發現：大壩周期彈性形變與溫度呈明顯的負相關，且變化趨勢幾乎同步；大壩形變與水位變化呈正相關且存在一定的相位差。

圖6 SSA提取形變、溫度、水位的周期分量

大壩在自重和水壓等荷載作用下可能發生徐變，導致大壩明顯地向上游或下游趨勢性變形。溫度和水位因子對大壩變形的影響主要為：溫度升高，大壩向上游的變形位移增大，相反溫度降低，則向下游的位移增大；水位升高，其向下游的變形位移增大，反之則向上游的位移增大[18]。為進一步分析趨勢項和周期項與影響因子的相互關系，求解大壩變形趨勢分量、周期分量(兩個周期分量之和)和余量與影響因子(時效、溫度、水位)的相關系數(見表1)。大壩變形的趨勢分量主要與時效因子有關，相關系數為-0.957 0，表明徐變方向為整體向河流上游(向河流上游為負，下游為正)。大壩變形的周期成分與溫度因子的相關系數為-0.870 4，與考慮相位差的水位因子相關系數為0.489 4。結合圖6可以看出，溫度升高大壩向上游變形，降低則向下游變形；水位升高大壩向下游變形。降低則向上游變形。這與溫度和水位因子影響大壩變形的規律基本一致。

表1 大壩變形與各影響因素相關系數

3 大壩變形預測

受地質、氣候、施工方案等各種條件的影響，大壩在施工或運營過程中總會發生不同程度的變形，及時掌握并準確預測大壩變形狀態，對大壩安全具有重要意義。目前常用的預測方法主要有時間序列模型、灰色理論等，但由于大壩變形具有較強的非平穩性和非線性，這些方法的應用均受到了一定限制。SSA能有效地將時間序列中的信號和噪聲分離，根據圖3可知，前5個分量為趨勢和周期分量，其能量之和占總能量的99.38%，余量可歸為噪聲部分，因此信號的重構階數r=5。

為了避免計算時各因子的數量級相差過大而導致矩陣病態，建模前將時效因子數據歸化為年序列。取樣本數據中前2200個作擬合建模，后114個用于檢核。文中分別對SSA提取的不同分量建模再疊加作為大壩變形的預測結果，流程如圖7所示。

圖7 大壩變形預測流程

采用式(6)擬合變形信號中的趨勢分量

f(t)=P0+P1t+P2t2+P3t3+P4t4

(6)

采用式(7)擬合提取的第二和第三個分量，即第一個周期分量

g(t)=a0+a1cos(ω1t)+a2sin(ω1t)

(7)

采用式(8)擬合提取的第四和第五個分量，即第二個周期分量

h(t)=b0+b1cos(ω2t)+b2sin(ω2t)+b3cos(2ω2t)+b4sin(2ω2t)+b5cos(3ω2t)+b6sin(3ω2t)

(8)

最后根據最小二乘原則解得式(6)—式(8)的各項系數(見表2)。

表2 模型系數

為了驗證上述方法的效果，將文中模型與傳統的多元線性回歸(multiple linear regression，MLR)模型作比較，根據確定性因子時效、溫度、水位建立多元線性回歸模型如下

y=β0+β1t+β2T+β3h+ε

(9)

式中，t、T、h分別代表時效(年)、溫度(℃)和水位(mm)。求得各項系數β0、β1、β2、β3分別為0.802、-0.133、-0.163、-0.070。

從圖8和圖9可知，SSA預測與多元線性回歸預測結果相比，其更逼近真實值，預測殘差明顯小于多元線性回歸模型，預測精度較高。為了定量評價模型的預測性能，利用均方根誤差(RMS)和平均絕對誤差(MAE)兩個指標對模型的預測效果進行評估

(10)

(11)

圖8 兩種模型預測結果與真實數據對比

圖9 兩種模型預測結果殘差對比

SSA方法擬合和預測的RMS值為0.52、0.24 mm，明顯小于多元線形回歸模型的0.81和0.57 mm，具體精度指標見表3和表4。結果表明，SSA方法的擬合及預測精度明顯優于多元線性回歸，能更準確地預測大壩變形。

表3 兩種模型的擬合精度評價指標 mm

表4 兩種模型的預測精度評價指標 mm

4 結 語

本文利用SSA方法對大壩變形時間序列進行分析，提取了趨勢項和周期項。通過分析發現，大壩存在徐變和周期性彈性形變，對周期性彈性形變的分析發現，溫度因素對其影響作用明顯大于水位。最后利用提取的趨勢和周期分量對大壩變形時間序列進行擬合并預測，并與傳統的多元線性回歸進行對比。通過實例分析，SSA與多元線性回歸模型的擬合和預測的均方根誤差分別為0.52、0.24 mm和0.81、0.57 mm；平均絕對誤差分別為0.36、0.20 mm和0.61、0.47 mm。SSA方法的擬合和預測誤差均小于多元線性回歸，表明本文方法能更為準確地預測大壩的變形。但本文僅采用了單一變形點數據，不能從整體上反映大壩變形的內部規律，分析多點的大壩整體變形是下一步研究的重點。