王 成 何 啟

(1.中交第四航務(wù)工程局 第七工程有限公司,廣州 510230;2.河海大學(xué) 水利水電學(xué)院,南京 210098)

混凝土壩可視為一個不斷變化的復(fù)雜系統(tǒng),特殊的結(jié)構(gòu)形式和復(fù)雜的工作環(huán)境,導(dǎo)致大壩的服役性態(tài)呈現(xiàn)出復(fù)雜、多樣、不確定等特點.在諸多反映混凝土壩性態(tài)變化的效應(yīng)量中,壩工界普遍認(rèn)為,大壩的變形能最直觀地反映出大壩的運行性態(tài)[1].在混凝土壩長期服役的過程中,積累了大量的變形監(jiān)測數(shù)據(jù),這些數(shù)據(jù)包含了大量的混凝土壩變形性態(tài)信息,為分析混凝土壩變形性態(tài)提供了有效的資料,也為變形性態(tài)的預(yù)測提供了基礎(chǔ)信息.

在對混凝土壩的變形情況進(jìn)行監(jiān)測時,自動化系統(tǒng)中儀器的數(shù)據(jù)采集一般都是基于一定時間間隔的,如每六小時一測或每天一測等,且通常保持不變；而人工監(jiān)測系統(tǒng)由于人為原因、儀器損壞和數(shù)據(jù)丟失等因素的影響,監(jiān)測信息之間的時間間隔往往并不能保持完全一致,表現(xiàn)為不等間隔的時間序列,這將給建模工作帶來困難.需要特別說明的是,在自動化系統(tǒng)和人工系統(tǒng)中,有時由于儀器長時間損壞,導(dǎo)致監(jiān)測信息會出現(xiàn)長時間的缺失問題,比如某測點兩年的數(shù)據(jù)序列中某一個月或某兩個月的連續(xù)缺失,這種類型的數(shù)據(jù)缺失問題將直接導(dǎo)致無法獲得該時間段內(nèi)該測點的變形情況,當(dāng)該測點位于大壩關(guān)鍵監(jiān)測位置時,很可能錯過大壩異常變形情況的監(jiān)測而給大壩的安全分析工作帶來困難.另外,數(shù)據(jù)的長時間缺失導(dǎo)致信息的中斷不利于變形情況的整體分析,變形預(yù)測工作也不方便展開.因此,需要研究相關(guān)方法實現(xiàn)缺失數(shù)據(jù)的彌補.呂開云[2]指出觀測數(shù)據(jù)插補的方法主要包括內(nèi)在物理聯(lián)系插補法和數(shù)學(xué)插補法兩種,并介紹了線性插補法的原理和過程；李雙平等[3]對比了常用的數(shù)學(xué)插值方法,選擇了充分利用已有數(shù)據(jù)信息、插值曲線光滑的三次Hermite分段插值,實現(xiàn)了數(shù)據(jù)序列的均勻化處理；屠立峰等[4]針對傳統(tǒng)插值函數(shù)在兩端插值區(qū)間易出現(xiàn)“龍格現(xiàn)象”的弊端,充分發(fā)揮了分形插值法能通過物體的部分信息推求整體性態(tài)的優(yōu)勢,將其運用到缺失時間序列的插值計算中；胡添翼[5]運用空間鄰近點的變形信息來回歸目標(biāo)測點的變形值,提出了空間鄰近點回歸插值法和空間反距離加權(quán)插值法.

本文將混凝土壩變形監(jiān)測信息視為時間序列,將信息的不完整問題分為兩類來處理,一類是時間序列的單值缺失問題(不均勻問題),即存在小部分時間間隔與整體不一致的數(shù)據(jù)序列時,對其進(jìn)行均勻化處理；另一類是連續(xù)多個數(shù)據(jù)缺失問題,對其進(jìn)行估計補充.

1 單值缺失處理方法

傳統(tǒng)的混凝土壩單值缺失處理方法常用插值法進(jìn)行估計,主要包括分段線性插值、多項式Lagrange插值、Newton 插值、三次樣條插值和三次Hermite插值等,這些插值方法對單值缺失數(shù)據(jù)進(jìn)行估計時,并不會影響到變形時間序列的整體變化趨勢和規(guī)律,因此當(dāng)缺失信息較少時,可采用此類插值方法進(jìn)行補齊,構(gòu)成等間隔的變形時間序列.但這些傳統(tǒng)插值方法只是基于已知的數(shù)據(jù)本身,并沒有過多考慮實際問題的物理意義,而變形時間序列單值缺失問題的處理是混凝土壩未知時間點變形信息的補充,需要考慮到混凝土壩的實際變形規(guī)律.

為此,本文引入非局部平均(Non-local Means,NLM)思想[6-7],提出了混凝土壩單值缺失處理改進(jìn)的非局部平均插值法,方法思路為：利用變形信息的非局部知識信息,以及變形序列中不同時刻信息規(guī)律的自相似性來對缺失時刻的變形值進(jìn)行估計,在此基礎(chǔ)上,引入完整的與目標(biāo)測點變形趨勢相關(guān)性最強的變形序列作為計算依據(jù).該方法旨在通過綜合考慮變形序列不同時刻變形值之間的自相關(guān)性,以及與目標(biāo)測點位置相當(dāng)?shù)臏y點之間的相關(guān)性來刻畫出缺失信息的特征.

假設(shè)壩體某變形測點A 的測值存在不均勻現(xiàn)象,為了估計出單個缺失數(shù)據(jù),進(jìn)行以下步驟.

首先,從測點的整個變形時間序列來看,找出與測點A 的變形趨勢相關(guān)性最強且序列完整的測點B,可從同一個監(jiān)測垂線上的諸多測點來尋找.本文采用Pearson相關(guān)性檢驗方法來計算測點變形數(shù)據(jù)之間的相關(guān)性,Pearson相關(guān)系數(shù)是一種用來定量衡量變量之間相關(guān)關(guān)系的統(tǒng)計學(xué)參數(shù),其計算公式為：

式中：δAi、δBi分別表示A、B測點同一時刻的變形值；N表示序列的總個數(shù).

從式(1)可以看出,Pearson相關(guān)系數(shù)的值在-1到1之間變化,且相關(guān)系數(shù)的絕對值越大,代表兩變量之間的相關(guān)性越強.當(dāng)相關(guān)系數(shù)越接近于1或-1時,相關(guān)性越強；當(dāng)相關(guān)系數(shù)越接近于0時,相關(guān)性越弱.另外,當(dāng)相關(guān)系數(shù)大于0時,兩變量為正相關(guān)；反之為負(fù)相關(guān).

其次,將測點B 的變形時間序列中與測點A 序列中待求插值點同一時間的變形值稱作假設(shè)插值點,計算測點B的序列中其它點對此假設(shè)插值點變形值的權(quán)重.本文采用歐幾里德距離(又稱歐式距離)的平方(Square of Euclidean Distance,SED)來度量不同時間點對應(yīng)變形值的相似性.歐式距離的平方計算公式為：

式中：δBi、δBj分別表示B 測點i、j時刻對應(yīng)的變形值大小.

通常情況下,不同時刻的變形值和之間的差值越小,說明兩個時刻的變形越相似,計算時賦予的權(quán)重值也越大,權(quán)重采用如下公式計算：

式中：h為控制指數(shù)函數(shù)增減速度的參數(shù),決定著權(quán)重的大小.

最后,將基于測點B完整變形序列求出的各個參考點相對于假設(shè)插值點的權(quán)值賦予測點A 對應(yīng)時刻的變形值,再對其進(jìn)行加權(quán)平均即可求出插值點的數(shù)值,公式如下：

式中：I表示選取的整個時間序列的時刻集.

假設(shè)測點A 和測點B的變形時間序列如圖1所示,其中測點B的序列完整,測點A 的序列中有缺失情況.圖中圓點表示不同時刻對應(yīng)的變形值大小,方點表示測點A 序列中的缺失值,也即待求的插值點.

圖1 單值缺失情況示意圖

2 多值缺失處理方法

當(dāng)變形時間序列中缺失的信息較多時,傳統(tǒng)的插值方法難以進(jìn)行有效的插值計算,而NLM 插值算法雖然可以求出每個缺失點的數(shù)值,但需要逐個計算參考序列中其它點對假設(shè)插值點的權(quán)重值,再計算目標(biāo)序列中各點對插值點的權(quán)重.該方法雖然可行,但計算工作量大.針對上述情況,本文研究混凝土壩變形信息的多值缺失處理方法.

目前常用的多值缺失處理方法主要有非線性回歸分析插值法和空間鄰近點回歸插值法,但兩種方法都有其局限性.當(dāng)回歸模型對變形序列的擬合精度較低或缺失段的環(huán)境量未知時,非線性回歸分析插值法精度較低；空間鄰近點回歸插值法借鑒了統(tǒng)計模型的建模思想,僅對變量的有限個整數(shù)項冪級數(shù)展開進(jìn)行回歸,難以全面刻畫測點變形之間未知的作用關(guān)系.

考慮到空間測點變形之間復(fù)雜未知的作用關(guān)系難以用具體的數(shù)學(xué)表達(dá)式進(jìn)行表征,而BP神經(jīng)網(wǎng)絡(luò)具有強大的非線性映射能力,經(jīng)過對樣本的學(xué)習(xí)訓(xùn)練可以勾畫出數(shù)據(jù)背后復(fù)雜的信息關(guān)系；與此同時,BP神經(jīng)網(wǎng)絡(luò)還具備強大的泛化能力,訓(xùn)練好的網(wǎng)絡(luò)可以實現(xiàn)對新輸入樣本的有效處理,給出合適的輸出結(jié)果.因此,為了提高缺失值估計的精度,盡可能找出最接近缺失時刻變形真實值,本文引入BP神經(jīng)網(wǎng)絡(luò)來處理空間測點變形之間的未知關(guān)系,由此提出相應(yīng)缺失值的估計方法.

假設(shè)某混凝土壩壩體有n個在空間上鄰近且結(jié)構(gòu)上相關(guān)的監(jiān)測點,如混凝土重力壩同一條垂線上的測點或混凝土拱壩同一變形分區(qū)內(nèi)的測點(分區(qū)方法本文不展開說明),當(dāng)?shù)趇個測點的變形信息由于某些原因出現(xiàn)缺失,即可利用其它m=(n-1)個測點的已知信息來估計i點的信息.基于BP神經(jīng)網(wǎng)絡(luò)映射的多值缺失估計方法建立步驟如下：

設(shè)樣本集中含有輸入向量和輸出向量之間的Z個模式對,隨機取一個模式對k,輸入模式向量為A k,期望輸出向量為；中間層單元輸入向量為S=(s1,s2,…,s p)(p為隱含層節(jié)點的數(shù)目,下同),輸出向量為B k=(b1,b2,…,b p)；輸出層單元輸入向量為L k=(l1,l2,…,l p),輸出向量為C=(c)；輸入層與隱含層之間的連接權(quán)為w(w=w ij,i=1,2,…,m；j=1,2,…,p)；隱含層與輸出層之間的連接權(quán)為v(v=v j,j=1,2,…,p)；隱含層各單元輸出閾值為θ(θ=θj,j=1,2,…,p)；輸出層單元的輸出閾值為γ=(γ).

1)網(wǎng)絡(luò)參數(shù)初始化,通過隨機賦值函數(shù)給w、v、θ和γ隨機賦一個(-1,1)之間的較小值；

2)用輸入向量A k、連接權(quán)w和閾值θ計算隱含層的輸入S；用S通過Sigmoid函數(shù)計算隱含層的輸出B k,即：

3)用隱含層的輸出B k、連接權(quán)v和閾值γ計算輸出層單元的輸入L k,再用L k計算輸出層單元的輸出向量C,即：

4)用期望輸出向量Y k、網(wǎng)絡(luò)實際輸出C計算輸出層單元的一般化誤差d k,即：

5)用連接權(quán)v、輸出層的一般化誤差d k和隱含層的輸出Bk計算隱含層各單元的一般化誤差e k,即：

6)用輸出層單元的一般化誤差d k、中間層各單元的輸出B k修正連接權(quán)v和閾值γ,即：

式中：η代表學(xué)習(xí)效率,取η=0.01～0.8；α為動量因子,取α=0.9.

近二三十年來，國內(nèi)學(xué)者通過引進(jìn)和學(xué)習(xí)西方的數(shù)學(xué)教育理論和方法，增強國際間的交流與合作，這是現(xiàn)在與世界進(jìn)行對話的基礎(chǔ)．但是，與數(shù)學(xué)不同的是，數(shù)學(xué)教育有很強的民族性、地域性，如何基于中國的民族性和地域性建立中國數(shù)學(xué)教育理論體系與研究規(guī)范，并以此為基礎(chǔ)建立自己的話語權(quán)，進(jìn)而與世界對話，融入世界學(xué)術(shù)圈？

8)隨機選取訓(xùn)練樣本集中另一個學(xué)習(xí)模式對,重復(fù)步驟3)～6),直至所有的模式對訓(xùn)練完畢；

9)計算網(wǎng)絡(luò)全局誤差函數(shù)E,其計算式為：

若E小于預(yù)先設(shè)定的一個誤差值,則網(wǎng)絡(luò)停止學(xué)習(xí)；否則重復(fù)步驟3)～8),進(jìn)行樣本集的下一輪學(xué)習(xí)訓(xùn)練；

10)將訓(xùn)練好的網(wǎng)絡(luò)保存,輸入新樣本,得出缺失信息估計輸出結(jié)果.

以n=5為例,輸入層即為4個相關(guān)測點的變形序列,輸出層即為目標(biāo)測點的變形序列,網(wǎng)絡(luò)結(jié)構(gòu)如圖2所示.

圖2 BP神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)示意圖

3 工程應(yīng)用

我國某混凝土重力壩壩頂高程384 m,最大壩高162 m；為監(jiān)測大壩的水平變形,在各重要壩段均布置了垂線.本文以泄1壩段正垂線上的各測點監(jiān)測數(shù)據(jù)為例進(jìn)行分析,驗證本文所提出的不完整信息處理方法的有效性,6 個測點的水平變形過程線如圖3 所示.可以發(fā)現(xiàn),這些測點的水平變形過程線具有較強的相關(guān)性.

3.1 單值缺失估計

以圖3中的測點PL5-3為例,首先尋找與該測點變形序列相關(guān)性最強的參考序列,目標(biāo)測點序列與其它各測點序列的相關(guān)性計算結(jié)果見表1.

圖3 泄1壩段正垂線測點水平變形過程線

表1 目標(biāo)測點與各鄰近測點變形序列的相關(guān)性結(jié)果

從上表的計算結(jié)果可以看出,測點PL5-4的變形序列與目標(biāo)序列的相關(guān)性最強,變形規(guī)律最相似,因此,參考測點即為測點PL5-4.

對于2014年8月31日到2014年9月18日的變形序列,共有19個變形數(shù)據(jù),見表2.

表2 目標(biāo)測點與參考測點的變形數(shù)據(jù)匯總表

假設(shè)2014年9月10日的變形數(shù)據(jù)是缺失的,用本文提出的基于非局部平均的插值方法和傳統(tǒng)方法分別對缺失值進(jìn)行估計,結(jié)果見表3.

表3 各插值方法的估計結(jié)果對比 (單位：mm)

從各插值方法的計算結(jié)果可以看出,本文提出的基于非局部平均的單值缺失估計方法的估計結(jié)果接近原監(jiān)測值.同時可以發(fā)現(xiàn),當(dāng)缺失值不在前后兩個數(shù)值區(qū)間范圍內(nèi)時,傳統(tǒng)的插值方法難以有效估計出此類變形值；而NLM 插值方法克服了此局限性,該方法利用變形序列自相似性的同時引入?yún)⒖夹蛄?提高了缺失值估計的精度.

3.2 多值缺失估計

同樣以PL5-3測點的變形序列為例,構(gòu)造長達(dá)一個月的缺失段(2014 年9 月1 日 至2014 年10 月2日).利用BP 神經(jīng)網(wǎng)絡(luò)建立同一垂線上的其它測點與目標(biāo)測點變形序列之間的映射關(guān)系.首先將各測點均已知的變形數(shù)據(jù)組成訓(xùn)練樣本導(dǎo)入BP 神經(jīng)網(wǎng)絡(luò)中進(jìn)行學(xué)習(xí),其次將其它測點缺失時間的變形值組成新樣本導(dǎo)入訓(xùn)練好的網(wǎng)絡(luò)中計算出目標(biāo)測點的缺失數(shù)據(jù).各方法的計算結(jié)果如圖4所示,其估計精度見表4,其中可決系數(shù)

圖4 缺失值估計結(jié)果對比圖

表4 缺失值估計精度對比

由圖4和表4的結(jié)果可以看出,本文提出的基于空間鄰近點BP 映射的多值缺失估計方法的估計精度均高于其它3種方法,可決系數(shù)和均方根誤差均達(dá)到比較理想的效果；空間鄰近點插值方法的估計效果也較好,但該方法僅利用了上下兩個測點的變形信息進(jìn)行回歸分析,沒有充分挖掘出目標(biāo)測點變形的相關(guān)信息；非線性回歸利用了統(tǒng)計模型的思想,在環(huán)境量變化不大的情況下可以采用該方法進(jìn)行估計；線性插值的效果較差,難以適用于多值缺失的估計.

4 結(jié)論

為解決混凝土壩變形信息的缺失問題,本文提出了基于非局部平均思想的單值缺失估計方法和基于空間鄰近點BP映射的多值缺失估計方法；工程應(yīng)用驗證表明,本文提出的兩種方法計算簡單,易于操作,且均具有較高的精確度,可以較好地估計出變形時間序列中的缺失信息.但是,這兩種方法仍有待改進(jìn)的地方：1)基于非局部平均思想的單值缺失估計方法,需要確定時間序列中所有已知時刻變形對缺失時刻變形的權(quán)重,整個時間序列長度的選擇標(biāo)準(zhǔn)值需進(jìn)一步的研究；2)基于空間鄰近點BP映射的多值缺失估計方法,通過BP神經(jīng)網(wǎng)絡(luò)來刻畫目標(biāo)測點變形和鄰近點變形的映射關(guān)系,但難以用數(shù)學(xué)顯式刻畫測點變形之間的關(guān)系,需要進(jìn)一步研究能充分反映混凝土壩各部分變形間相互關(guān)系的直觀表達(dá)式.