混凝土壩變形監測缺失數據處理新方法

黃成章，顧沖時，何 菁

(1．河海大學水利水電學院，江蘇 南京 210098； 2.南京市水利建筑工程檢測中心有限公司，江蘇 南京 210036)

混凝土壩運行過程中不僅受到各種荷載作用，還受到各種非荷載因素的影響。混凝土壩的服役安全性態是一個在多因素協同作用下的材料和結構互相影響的非線性動態演化過程，變形作為一個可以直觀反映混凝土壩安全狀態的綜合效應量，是大壩結構性態發生趨勢性變化乃至轉異的一項重要指標[1]，因此，對于混凝土壩變形的監測尤為重要。然而，由于大壩監測儀器發生故障或者大壩自動化監測系統不穩定等原因，可能導致部分監測點存在監測數據缺失的現象，會影響大壩安全監測及評價的效果[2]，因此，研究在監測儀器故障情況下的不完整信息處理方法具有現實意義。目前較為常見的處理方法有加權平均法[3]、貝葉斯估計法[4]、卡爾曼濾波法[5]、遞歸加權最小平方法、極大似然估計法[6]等，但這些方法均以單個測點作為對象進行插值處理，忽略了測點所處位置及各測點之間的相互聯系，影響了處理缺失數據的精度和有效性。

實際上，混凝土壩所有測點的變形序列既包含了時間維度上的變形信息，也包含了橫截面維度上的變形信息，如果僅僅從單一的時間序列或者截面序列來分析大壩變形信息，會難以全面有效地了解大壩的整體變形狀態。隨著時空數據挖掘技術的不斷進步[7-13],與時間或空間相關的模型得到進一步發展[14-16]。例如：顧沖時等[17]通過引入測點位置坐標作為大壩變形的影響因子,建立了大壩時空分布模型；李廣春等[18]通過考慮各測點間的空間相關性,建立了重力壩變形時空自回歸模型。

本文基于聚類和面板數據理論，提出了一種考慮時間和截面兩種維度的混凝土壩變形監測缺失數據處理方法，旨在彌補傳統方法未考慮測點間相關性和非荷載因素影響而帶來的處理缺失數據精度不高的缺陷。

1 測點監測效應量相似性判據及其聚類

1.1 相似性判據

聚類的含義是按照一定的標準和指標將需要處理的數據集分成不同的組，稱之為簇，經聚類后，同一個簇內的數據具有較大的相似性，而不同簇內的數據則具有較大的差異性，通過聚類使得在同一區域內的混凝土壩變形規律相似程度盡可能地“接近”，而不同區域之間的變形規律相似程度盡可能地“遠離”。本文在分析混凝土壩變形監測值時以歐氏距離來[19]刻畫這種相似性，以此來對測點進行分區。在對變形監測數據預處理時，設{δit|i=1,2,…,N;t=1,2,…,T}(N為混凝土壩監測點的總數，T為監測時期總數)為混凝土壩變形監測數據集，即監測時間序列，采用絕對距離和增速距離[20]刻畫各監測點的測值相似性。絕對距離是指不同測點在整個監測時期距離的遠近程度，增速距離是指不同測點的變形增量隨時間變化的趨勢差異，計算公式分別為

(1)

(2)

其中 Δδit=δit-δi(t-1)Δδjt=δjt-δj(t-1)

式中：dij1、dij2分別為測點i與測點j之間的絕對距離和增速距離；δit、δi(t-1)、δjt、δj(t-1)分別為t與t-1時刻測點i、j處的變形量；Δδit、Δδjt分別為t與t-1時刻測點i、j處變形量差值。

從式(2)可以看出，若測點i與測點j的變形值隨時間呈現同向變化的趨勢，且這種變化越協調，那么這兩者也就越相似，dij2也就越小；反之則相似性較差，dij2越大，這也符合相似性度量的基本原則。為了能夠更好地刻畫各測點變形的相似性，構建測點i、j之間的綜合距離dij3：

dij3=w1dij1+w2dij2

(3)

式中w1、w2分別為絕對距離和增速距離的權重，w1+w2=1。

在信息論中，熵是對不確定性的一種度量，信息量越大，不確定性就越小，熵也就越小；信息量越小，不確定性越大，熵也越大。本文基于熵權法來判斷混凝土壩變形變化隨機性及無序程度，進而用于判斷反映變形變化特性的某個指標的離散程度，指標的離散程度越大，該指標對綜合評價的影響越大，據此確定絕對距離和增速距離的權重。假設有m個評價對象，n個評價指標，xab為評價對象b(b=1，2，…，m)在評價指標a(a=1，2，…，n)下的評價值。本文以變形量的絕對距離和增速距離作為評價指標，測點兩兩間的距離作為評價對象。熵權法賦權步驟如下(由于各項指標的計量單位統一，故不需要進行指標的歸一化處理):

步驟1計算評價指標a下評價對象b評價值占該指標評價值的比重pab：

(4)

步驟2計算評價指標a的熵值ea：

(5)

其中c=1/lnm

步驟3計算評價指標a的信息熵冗余度(差異)da：

da=1-ea

(6)

步驟4計算評價指標a的權重wa：

(7)

在n=2時，求得w1與w2代入式(3)即可得到不同測點之間的綜合距離，由此構建所需的相似性判據。

1.2 監測效應量聚類

采用離差平方和法對混凝土壩變形各測點監測效應量進行聚類。假設將混凝土壩N個測點劃分為k個區域，記為G1、G2、…、Gk，記Nl為Gl(l=1,2,…,k)類中測點的個數，Xil為Gl類中測點i(i=1,2,…,Nl)處的變形值，那么混凝土壩Gl區域在T個時期內各測點序列的離差平方和為

(8)

式中：Sl為Nl個測值的總離差平方和；S為k個區域總的離差平方和；Xit為t時刻測點i處變形值；Yit為Gl區域中t時刻測點i處變形相對增量。

通常情況下，要得到分區的最優解，需要選擇使得S取得最小值的方案。在實際工程中，由于混凝土壩結構的復雜性，往往難以求得一個最優分區數。因此，本文假設在聚類過程中，共進行了p次合并，求得第l次與第p次合并的區域間距離之比Rl，即Rl=Dl/Dp，其中Dl、Dp分別為第l次與第p次合并的區域間距離。如果Rl與Rl+1相差較小，但Rl與Rl-1卻相差較大，那么就可以將第l次合并后得到的區域間距離Dl作為變形分區的一個閾值，進而得到混凝土壩變形分區數量[21]。構建混凝土壩變形分區的步驟如下：

步驟1由式(1)、式(2)分別計算混凝土壩變形的絕對距離和增速距離。

步驟2采用熵權法計算兩種距離的權重。

步驟3將計算得到的權重代入式(3)得到不同測點之間的綜合距離以及各區域之間的距離矩陣D0。

步驟4初始將N個測點分為N個區域，區域個數k=N,令D1=D0。

步驟5計算區域間的距離矩陣Di-1，根據離差平方和最小準則，將綜合距離最小的兩個區域合并，組成一個新的區域。

步驟6計算新的區域與其他區域的綜合距離，進而得到新的距離矩陣Di。，不斷重復步驟5和步驟6，直到所有的測點全部被劃分到相應區域為止。

步驟7根據閾值確定變形的最優分區數，得到混凝土壩變形最優分區數K。

步驟8利用譜系聚類樹狀圖畫出混凝土壩變形分區分布圖。

為得到存在缺失數據的混凝土壩變形分區，根據壩體各部位隨時間變化變形趨勢的一致性，取各測點變形監測數據完整時間段的監測效應量進行聚類分區即可。

混凝土壩變形分區的具體流程見圖1。

圖1 混凝土壩變形分區流程

2 混凝土壩變形面板數據隨機系數模型的建模方法

為了對混凝土壩變形監測缺失數據進行處理，需建立表征混凝土壩變形規律的分析模型，本文基于面板數據理論構建大壩變形變化規律分析模型。

2.1 面板數據隨機系數模型的表達形式

面板數據隨機系數模型的一般表達式為

(i=1,2,…,N;t=1,2,…,T)

(9)

式中：βi為隨機系數；Yit為混凝土壩i個測點因變量變形面板數據序列；uit為均值為0、方差為σ2，并且滿足獨立同分布的隨機誤差成分；yit為測點i在t時刻的變形監測值；Xit為混凝土壩i個測點自變量變形面板數據序列。自變量Xit中的元素為xit，xit=(1,Ht,1,Ht,2,Ht,3,Ht,4,T1,t,T2,t,…,Tm,t,θt,lnθt)T，其中，Ht,1,Ht,2,Ht,3,Ht,4為變形監測數據的水壓分量影響因素，T1,t、T2,t、…、Tm,t為溫度分量影響因素，θt和lnθt均為時效分量影響因素。

由于混凝土壩在運行過程中會受到各種隨機因素的影響(如約束條件、荷載作用、材料性質等)，且不同部位的變形受到的影響存在較大差異，故模型中還應考慮這些因素對混凝土壩不同測點處所產生的異質性影響。由于這些影響因素無法明確地在模型自變量中表示，故采用特異效應量αi來代表這些個體特異效應。當βi被看作是隨機系數時，假設βi=β+αi，將參數分解成兩部分，一部分為固定常數，另一部分為混凝土壩不同測點變形的隨機效應量。假定在不同橫截面個體之間的協方差為零，對角線元素的乘積為對應的方差，且隨機系數在不同的橫截面個體之間的協方差也為零，則對于第i個橫截面單元，可以得到Yi=Xiβi+ui。

對于隨機系數回歸模型，采用普通最小二乘法估計將會導致無法估計以及錯誤的統計結果，可采用廣義最小二乘法(GLS)估計來解決上述問題[21]，進而得到βi的漸進正態和有效估計量。

2.2 變形監測缺失數據處理步驟

實際工程中，混凝土壩變形監測中出現的數據缺失現象一般有兩種情況：一是在某時間段存在一個或多個測點測值缺失，一是在某時間段全部測點測值缺失。對于這兩種情況，均可采用前文聚類、建模方法對變形監測缺失數據進行處理。對于第一種情況的處理步驟如下：

步驟1利用監測數據完整的時間段部分(缺失數據之前或之后部分)的監測值，對混凝土壩各測點監測效應量進行聚類分區，使變形規律相似、對荷載具有同質響應的部位被歸為一類，將傳統的點分析方法轉變為區域分析方法，在較大程度上避免僅僅考慮壩體局部變形而造成的偏誤性判斷。

步驟2對于同一區域內測點缺失數據之前部分以及之后部分分別建立面板數據隨機系數模型，引入特異效應量αi來代表個體特異效應，以提高模型的可靠性和有效性。

步驟3基于缺失數據前后兩部分數據建立的模型，分別得到各測點缺失部分數據的擬合值，取兩擬合值的平均值對大壩變形監測缺失數據進行填補。

步驟4對比監測數據完整測點處對應缺失數據點時段的監測值與模型計算值，驗證處理方法的準確性與有效性。

第二種情況處理步驟和第一種情況基本相同，僅步驟4改為：利用位于缺失數據之前或之后部分的連續完整時間段的監測值，人為刪除部分監測值數據，重復步驟2、3，驗證填補值的準確性與有效性。

2.3 評估指標

為衡量混凝土壩變形監測缺失數據填補的有效性，本文將處理后得到的殘差與SL 601—2013《混凝土壩安全監測技術規范》中所規定的拱壩壩體徑向位移量誤差限值(±2 mm)進行比較，若在誤差限值之內，則說明填補數據有效。

3 實例驗證

某混凝土雙曲拱壩壩頂高程1 885.0 m，建基面最低高程1 580.0 m，最大壩高305 m，拱冠梁頂厚16 m、底厚63 m，最大中心角93.12°，頂拱中心線弧長552.23 m，厚高比0.207，弧高比1.811，柔度系數7.99，壩體混凝土方量476萬m3。以布設于大壩11號、13號、16號壩段，高程分別為1 885.00 m、1 829.25 m、1 778.25 m、1 730.25 m和1 664.25 m的15個正垂監測點為研究對象，變形監測序列為2017年12月24日至2018年12月31日的測值，測點編號見表1。

表1 混凝土壩變形測點編號

3.1 某時間段單個測點測值缺失處理及有效性分析

11號壩段PL11-3測點在2018年1月31日至3月22日變形監測數據缺失，采用本文方法對缺失數據進行處理。先對各測點變形監測量進行聚類分區，得到譜系聚類樹狀圖。通過內部有效性指標的多次評估，經過試錯迭代并結合最優分區閾值，確定混凝土壩變形分區的最佳區域數，最終將所有測點劃分為三大類，如圖2所示。經分析區域內的測點變形變化趨勢及規律基本一致，可以描述壩體對應區域的總體變形特征。

圖2 混凝土壩變形分區分布

在PL11-3測點所處的區域(第三類)中，對缺失數據前后兩個時間段各測點變形監測效應量采用面板數據隨機系數模型進行建模分析，分別得到缺失數據的模型計算值，取其平均值作為最終缺失數據處理結果。

為了驗證處理方法的有效性，基于所建立的面板數據隨機系數模型，針對具有完整實測數據的測點PL13-5、PL16-5，計算得到對應測點PL11-3缺失數據時段處PL13-5與PL16-5兩個測點處的變形值，其殘差如圖3所示，可以看出，填補數據殘差絕對值均遠小于2 mm，都在誤差限值之內，可見填補數據具有較高的精度。同時將面板數據隨機系數模型擬合值與三次樣條插值作對比，結果見圖4。由圖4可見，面板數據隨機系數模型擬合值具有更高的可靠性與有效性，面板數據隨機系數模型對于缺失數據處理具有較好的效果，得到的測點PL11-3全時段監測值(含原缺失數據處理后得到的填補值)過程線見圖5。

圖3 測點PL13-5和PL16-5變形量殘差

圖4 測點PL13-5變形量過程線對比

圖5 測點PL11-3全時段變形監測值過程線

3.2 某時間段全部測點測值缺失處理效果分析

為檢驗在某時間段全部測點測值缺失情況下本文方法的有效性，現假設所有測點2018年1月31日至3月22日變形監測數據全部缺失，通過對缺失數據之后時間段的監測數據進行處理來驗證方法的有效性。首先對各測點變形監測效應量進行聚類分區，此處選取圖2中PL11-3所在區域驗證缺失數據處理的有效性；其次選擇缺失數據之后部分(或之前部分)時間段的測點測值進行有效性檢驗。具體處理步驟如下：

步驟1假設2018年8月21日至10月10日變形監測數據缺失，利用對缺失數據前后兩個時間段各個測點監測數據，分別建立面板數據隨機系數模型，據此分別計算得到相應的模型計算值，取兩部分對應缺失數據時段的模型計算值的平均值作為填補的變形監測值。

步驟 2計算得到各時間點的殘差值，判斷是否在誤差限值(±2 mm)以內。

步驟 3對2018年1月31日至3月22日前后時間段的監測值分別進行建模分析，處理得到有效的變形監測數據填補值。

圖6 典型測點變形量殘差

圖6為典型測點處變形量殘差，可以看出填補數據均在誤差限值之內，因此，在某時間段全部測點測值缺失的情況下，本文方法也具有較高的精度。

4 結 語

本文研究了混凝土壩變形各監測效應量相似性判據構建技術，給出了基于熵權法判別大壩變形相似性的指標及各指標的賦權方法，并依據監測數據離差平方和最小原理，提出了混凝土壩變形各測點監測效應量聚類方法。綜合考慮同一聚類分區測點間變形監測效應量的相關性和非荷載因素的影響，提出了基于面板數據理論的混凝土壩變形監測缺失數據的處理方法，彌補了傳統方法的不足，并通過工程實例驗證了方法的有效性，可應用于混凝土壩變形監測缺失數據的處理。