依據主成分和協整性的大壩變形奇異診斷*

楊光，李姝昱，孫錦

（1.華北水利水電大學水利學院 鄭州，450046）

（2.河海大學水文水資源與水利工程科學國家重點實驗室 南京，210098）

（3.黃河水利科學研究院 鄭州，450003）

（4.華北水利水電大學測繪與地理信息學院 鄭州，450046）

引言

變形是大壩服役性態變化的綜合反映，是衡量結構安全與否的重要標志。科學地分析變形原位監測信息，不僅是現行規范的要求，而且是監控大壩安全的有效手段［1-2］。受多種因素的干擾，原位監測資料中常包含奇異成分，表現為測值在某時刻或某時段的異常突跳，即孤立型或斑點型奇異成分。若未有效診斷出監測資料中的奇異成分，極可能影響大壩安全監控結論的客觀性，造成虛假報警或將危險狀況遺漏。尤其對于運行多年的老壩，監測儀器的性能和穩定性較差，此問題尤為突出。因此，提出科學的大壩變形原位監測數據奇異成分診斷方法，具有重要的理論意義和應用價值。

大壩變形原位監測數據奇異成分診斷包括辨識和估計2個環節。常規的辨識方法有［3-4］：①過程線法，通過繪制變形過程線，直觀辨識奇異成分；②假設檢驗法，如拉依達檢驗、狄克松檢驗。傳統的估計方法包括［3-4］：①忽略法，若無法還原或難度較大，可直接忽略；②似然估計法，如臨近插值、三次Hermite插值、統計模型及自回歸模型。總體來看，常規診斷方法大多以變形時間序列為分析對象，無法有效地捕捉多測點間隱含的變形關系，且主觀性較強，易誤診。

相較于1維時間序列，多維數據序列具有信息量豐富、自由度高及穩定性強的優勢，包含了隱匿、新穎及有潛在價值的信息［5］。PCA［6-9］和CA［10-12］是兩種多維數據特征挖掘方法。PCA可將多維數據序列投影到不同方向上，得到主要規律和無法解釋成分，其中無法解釋成分包含了奇異數據、噪聲等。對于大壩上臨近的監測點，變形監測序列存在關聯性，若未出現超標準洪水、極端天氣等非常規狀況，某監測點變形數據異常，其他臨近點的數據未見異常，則可判定該監測點的變形數據奇異。CA刻畫了多維數據序列的長期關系，若未發生非常規狀況，多測點變形監測序列的協整關系一般不會明顯改變，據此可實現對奇異成分的估計。

綜上所述，筆者利用多維數據的優勢，基于PCA，研究奇異成分辨識準則；依據CA，探究奇異成分似然估計模型；在此基礎上，結合實際工程，檢驗PCA-CA方法的有效性。

1 診斷方法

1.1 奇異成分辨識

依據2維數據格式，建立大壩變形原位監測數據集合，即

其中：m為監測點總數；n為監測時長；xi，xj分別為測點i和j的監測序列；xki，xkj分別為測點i和j在k時刻的監測值。

PCA辨識步驟如下。

1）分析環境量監測資料，若未發生非常規狀況，則執行后續步驟。

2）對矩陣X進行標準化，得到矩陣,為標準化的xi，為 標 準化的xj。相關系數矩陣R為

其 中：rij為和的相關系數，計算公式為rij=(i，j=1，2，…，m)；和為標準化 的xki和xkj。

3）計算相關系數矩陣R的特征向量矩陣P和特征值矩陣λ，計算公式為

其中：0為m階零矩陣；I為m階單位矩陣。

式（3）中的P和λ分別表示為

其 中：λi為矩陣R的第i個特征值，滿足λ1≥λ2≥…≥λm；pi為λi對應的特征向量。

4）建立n×m矩陣T，記為(t1，t2，…，tj，…，tm)

5）計算第j項tj(j=1，2，…，m)對原始序列X的解釋能力ej

由式（7）可知：①前j項對X的累積解釋能力為，所有m項的累積解釋能力為1；②e1≥e2≥…≥em，即t1→tm對X的解釋能力逐漸變弱；③e1越大，X的主要規律越明顯。因此，為保證PCA辨識精度，要求e1盡可能大。

6）建立SPE統計量。由于P為正交矩陣，滿足

則式（6）可變換為

k時刻SPE統計量記為SPEk，計算公式為

7）依據假設檢驗，在顯著性水平α下，建立SPE控制限［13-14］，記為SPEα，即

假設檢驗中α常取為0.05或0.01，即接受原假設時，正確的概率為95%或99%。

8）圖1為SPE統計與控制限的關系。若SPEk未超過SPEα，如圖1（a）所示，則不存在奇異成分；若SPEk在某時段超過SPEα，如圖1（b）所示，則為斑點型奇異；若SPEk在某時刻超過SPEα，如圖1（c）所示，則為孤立型奇異。若奇異，則執行步驟9。

圖1 SPE統計量與控制限的關系Fig.1 Relationships between SPE and its control limit

9）測點j對k時刻SPEk超過SPEα的貢獻度用CSPEkj統計量表征，即

CSPEkj越大，則測點j對k時刻SPEk超限 的貢獻度越大，而CSPEkj最大處，即為存在奇異成分的測點。

1.2 奇異成分似然估計

假設m個變形監測序列中有f個存在奇異成分，記為x′1…x′h…x′f，剩余p=m-f個不存在奇異成分，記為。CA似然估計步驟如下。

1）建立多測點變形原位監測序列的協整模型，即

其中：cqg為的系數；aq為的最高次數，依據變形散點關系確定；ε為余量序列。

2）應用逐步回歸算法，確定式（14）的系數cqg，得到余量序列ε。

3）利用ADF檢驗［15-16］，測試ε的平穩性。若ε平穩，則協整，執行步驟4；若非平穩，則不協整，此時算法不適用。

4）利用式（14），估計x′h的 奇 異成分，修正矩陣X。

5）利用PCA準則，辨識矩陣X中是否仍然存在奇異成分，若存在，重復執行步驟1～4，直至無奇異成分。

2 工程實例分析

2.1 工程實例1

陳村拱壩最大壩高為76.3 m，設計汛限水位為117.0 m，汛后正常高水位為119.0 m，500年校核水位為123.8 m，5 000年校核水位為126.1 m，保壩水位為127.7 m。18正下（1號測點），26正下（2號測點）和29倒（3號測點）測點的位置如圖2所示，徑向變形過程線和置信區間法計算結果如圖3所示，環境量過程線如圖4所示。

圖2 29#，26#和18#壩段垂線監測點Fig.2 Vertical monitoring points of the 29#,26# and 18#dam sections

依據水壓-溫變-時效（hydraulic-seasonal-time，簡稱HST）模式［17］，建立變形監控模型，HST模型精度如表1所示。由表可知，3個模型的復相關系數R均大于0.95，且標準差S較小，擬合精度較高。應用置信區間法［17］，評判變形變化是否正常，2S為正常與基本正常的臨界值，3S為基本正常與異常的臨界值。由圖3（b）～（d）可知：①1號測點未超出3S，在42天超出了2S，未報警；②2號測點未超出3S，在23天超出了2S，未報警；③3號測點在2天超出了3S，且在148天超出了2S，發生報警。

表1 HST模型精度Tab.1 Accuracy of the HST models

圖3 徑向變形過程線和置信區間法計算結果Fig.3 Time series of radial deformation and calculation results of confidence interval approach

由圖4可知，該時段未出現非常規水位和溫度狀況，因此受奇異成分影響，監控模型虛假報警。利用本研究方法，診斷3個測點的奇異監測數據。經計算，e1為95.23%，且余量序列協整，診斷出的奇異監測數據有斑點型和孤立型，如：1號測點的變形監測值在2012-07-11為孤立型奇異；3號測點的測值在2011-05-25至2011-07-19期間為斑點型奇異。

圖4 環境量過程線Fig.4 Time series of environment factors

重新建立監控模型，結果顯示：①1號測點在3天超出了2S，未超過3S；②2號測點在5天超出了2S，未超 過3S；③3號 測點 在14天超 出 了2S，未超過3S。

2.2 工程實例2

斑點型奇異可視為孤立型奇異的集合，因此孤立型奇異辨識和斑點型奇異估計具有分析的代表性。結合錦屏一級拱壩3個垂線點PL11-3，PL11-4和PL13-3的徑向變形監測資料，通過人為構造奇異成分［18］，檢驗PCA-CA方法的性能。測點布置如圖5所示，圖6為徑向變形原位監測過程線。

圖5 11#，13#壩段垂線監測點布置Fig.5 Vertical monitoring points of the 11# and 13# dam sections

圖6 PL11-3，PL11-4和PL13-3的徑向變形監測序列Fig.6 Time series of radial deformation of PL11-3,PL11-4 and PL13-3

2.2.1 PCA辨識性能分析

圖7為PCA方法的診斷結果。依據PCA辨識步驟1～5，得到了t1～t3效應量，如圖7（a）～（c）所示，測點變形時間序列的相關系數矩陣如表2所示，t1～t3的特征值、解釋能力以及累積解釋能力如表3所示。依據步驟6，可得到SPEk過程線，在執行步驟7時，α取為0.01，計算得到控制限SPE0.01=0.11，由圖7（d）可知，SPEk＜SPE0.01，因此3個測點的變形原位監測數據不包含奇異成分。

圖7 PCA方法的診斷結果Fig.7 Diagnosis results of the PCA method

表2 相關系數矩陣Tab.2 Correlation coefficient matrix

表3 特征值、解釋能力及累積解釋能力Tab.3 Eigenvalues,explanatory capabilities and accumulated explanatory capabilities

PL11-4測點在2014-10-03的徑向變形監測值為39.4 mm，將其構造為不同程度的孤立型奇異成分，記作0#～6#情況，分別為39.6，39.9，40.4，40.9，41.4，42.4和43.4 mm。圖8為1#～5#情況下SPEk與SPE0.01的關系圖，圖9為1#～6#情況的CSPE統計量。由圖8可知，在2#～5#情況下，SPEk超過了SPE0.01，且 從 圖9可 以 看 出，在2#～6#情 況 下，PL11-4測點的CSPE統計量明顯高于PL11-3和PL13-3，因此2#～6#可判定為孤立型奇異成分。表4比較了PCA準則、拉依達準則、狄克松準則和t準則的性能，可以看出，4種方法分別可辨識出相對誤差為3.81%，7.61%，7.61%和5.08%的孤立型奇異成分。

表4 4種方法性能對比Tab.4 Performance comparison of four methods

圖8 1#～5#情況下SPEk與SPE0.01的關系Fig.8 Relations between SPEk and SPE0.01 under 1# case～5# case

圖9 1#～6#情況的CSPE統計量Fig.9 CSPE statistics of 1# case～6# case

2.2.2 CA估計性能分析

假設PL11-4測點在2016-07-15至2016-09-15期間的變形監測數據為斑點型奇異，利用2014-01-01至2016-07-14的原位監測資料建立協整模型，對奇異成分進行估計。圖10為3個測點徑向變形散點關系。由圖可知，式（16）的最高次數aq均取1，利用逐步回歸法，可得

圖10 PL11-3，PL11-4和PL13-3測點徑向變形散點關系Fig.10 Scatter relations among radial deformation of PL11-3,PL11-4 and PL13-3

圖11為CA方法的估計結果。PL11-4的監測值、擬合值和余量序列如圖11（a）所示，復相關系數為0.997 1。ADF檢驗結果顯示：t統計量的值為-2.724 4，1%和5%兩個顯著性水平的臨界值分別為-2.400和-2.150，因此t統計量小于臨界值，即原假設均被拒絕，余量序列ε平穩。圖11（b）為自回歸模型、統計模型和CA模型對奇異成分的估計結果，3種模型的復相關系數分別為0.743 2，0.871 5和0.994 5。

圖11 CA方法的估計結果Fig.11 Estimation results of the CA method

3 結論

1）為消除奇異成分對大壩變形安全監控結論客觀性的影響，運用PCA理論，建立了SPE和CSPE，借助假設檢驗，提出了PCA辨識準則，在此基礎上，基于CA原理，應用ADF檢驗和逐步回歸法，提出了CA似然估計模型。利用陳村拱壩和錦屏一級拱壩徑向變形原位監測資料，檢驗了PCA-CA方法的有效性。

2）拉依達準則、狄克松準則和t準則以1維時間序列為分析對象，辨識性能對樣本分布型式的依賴性較大，僅體現了概率意義。PCA準則針對多維數據序列進行分析，不僅包含概率含義，而且考慮了多測點變形的主成分關系，辨識性能最佳。

3）自回歸模型和統計模型著重分析1維時間序列，其中：自回歸模型僅考察了前期監測數據的時間波動規律，外延性較差；統計模型刻畫了水壓、溫變和時效因素，外延性優于自回歸模型，但屬于半經驗模型，估計精度仍不足；CA模型以多維數據序列為分析對象，既刻畫了變形隨時間的波動特征，亦表征了多測點變形間的協整關系，同時也反映出外部作用的綜合效應，外延性較好，取得了最佳的估計精度。

4）PCA-CA方法要求：①t1效應量對原始變形序列X的解釋程度e1盡可能大，通過工程實例1可知，當e1=95.23%時，PCA準則可有效辨識出奇異數據；②余量序列平穩。對于我國已修筑的部分老壩，布設的變形監測儀器較少，不易滿足上述條件，此時，仍需采用常規方法。