李世友，王奉偉，沈云中

(同濟(jì)大學(xué)測(cè)繪與地理信息學(xué)院，上海 200092)

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

1 奇異譜分析法

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

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

(1)

式中，M為窗口長度。對(duì)于非平穩(wěn)序列，根據(jù)軌跡矩陣計(jì)算協(xié)方差陣C如下

C=DDT

(2)

對(duì)協(xié)方差陣C進(jìn)行特征值分解，確定特征值λk及對(duì)應(yīng)的特征向量Vk(由大到小排列)。第k個(gè)特征值對(duì)應(yīng)的主成分為

(3)

(4)

利用所有M個(gè)主成分可完全恢復(fù)原始序列X=(x1,x2,…,xN)。如果前r個(gè)主成分的貢獻(xiàn)統(tǒng)計(jì)上屬于信號(hào)子空間，則利用這些主成分重構(gòu)時(shí)間序列的信號(hào)部分如下

(5)

其余部分屬于噪聲子空間，可用于重構(gòu)時(shí)間序列的噪聲。

2 大壩變形的奇異譜分析

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

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

圖1 大壩變形時(shí)間序列

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

圖3 SSA提取趨勢(shì)和周期分量及其貢獻(xiàn)率

圖4 溫度數(shù)據(jù)及其周期分量

圖5 水位數(shù)據(jù)及其周期分量

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

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

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

表1 大壩變形與各影響因素相關(guān)系數(shù)

3 大壩變形預(yù)測(cè)

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

為了避免計(jì)算時(shí)各因子的數(shù)量級(jí)相差過大而導(dǎo)致矩陣病態(tài)，建模前將時(shí)效因子數(shù)據(jù)歸化為年序列。取樣本數(shù)據(jù)中前2200個(gè)作擬合建模，后114個(gè)用于檢核。文中分別對(duì)SSA提取的不同分量建模再疊加作為大壩變形的預(yù)測(cè)結(jié)果，流程如圖7所示。

圖7 大壩變形預(yù)測(cè)流程

采用式(6)擬合變形信號(hào)中的趨勢(shì)分量

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

(6)

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

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

(7)

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

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

(8)

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

表2 模型系數(shù)

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

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

(9)

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

從圖8和圖9可知，SSA預(yù)測(cè)與多元線性回歸預(yù)測(cè)結(jié)果相比，其更逼近真實(shí)值，預(yù)測(cè)殘差明顯小于多元線性回歸模型，預(yù)測(cè)精度較高。為了定量評(píng)價(jià)模型的預(yù)測(cè)性能，利用均方根誤差(RMS)和平均絕對(duì)誤差(MAE)兩個(gè)指標(biāo)對(duì)模型的預(yù)測(cè)效果進(jìn)行評(píng)估

(10)

(11)

圖8 兩種模型預(yù)測(cè)結(jié)果與真實(shí)數(shù)據(jù)對(duì)比

圖9 兩種模型預(yù)測(cè)結(jié)果殘差對(duì)比

SSA方法擬合和預(yù)測(cè)的RMS值為0.52、0.24 mm，明顯小于多元線形回歸模型的0.81和0.57 mm，具體精度指標(biāo)見表3和表4。結(jié)果表明，SSA方法的擬合及預(yù)測(cè)精度明顯優(yōu)于多元線性回歸，能更準(zhǔn)確地預(yù)測(cè)大壩變形。

表3 兩種模型的擬合精度評(píng)價(jià)指標(biāo) mm

表4 兩種模型的預(yù)測(cè)精度評(píng)價(jià)指標(biāo) mm

4 結(jié) 語

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