基于穩(wěn)健估計的大壩變形監(jiān)測統(tǒng)計模型分析

朱趙輝，劉 健，李 新，尚 層，王國川

（1.中國水利水電科學(xué)研究院，北京 100048；2.北京中水科工程總公司，北京 100048；

3.新疆額爾齊斯河流域開發(fā)工程建設(shè)管理局，新疆 烏魯木齊 830000；4.丹東太平灣發(fā)電廠，遼寧 丹東 118000）

1 研究背景

大壩變形監(jiān)測數(shù)據(jù)中無法避免粗差的存在，粗差的處理常可分兩種途徑［1］：一是將粗差歸為均值漂移模型，在數(shù)據(jù)分析前進(jìn)行定位和剔除，以期得到一簇清潔干凈的數(shù)據(jù)；二是將粗差歸于方差膨脹模型，采用穩(wěn)健估計的方法通過逐次迭代來不斷地改變觀測值的權(quán)，最終使粗差觀測值的權(quán)趨于零，達(dá)到限制、排除粗差對建模分析的不利影響。

粗差定位理論中，傳統(tǒng)“3σ準(zhǔn)則”適用于監(jiān)測量符合正態(tài)分布的情況，無法適應(yīng)壩體變形所面臨的動態(tài)復(fù)雜條件（庫水位、溫度等）的變化；格羅布斯、狄克遜準(zhǔn)則適用于小樣本的粗差剔除［2-3］，在一定顯著性水平下，無法保證粗差被剔除；采用Barrada數(shù)據(jù)探測法由于每次只考慮一個粗差［4］，且未顧及各殘差間的相關(guān)性，檢驗可靠性受到一定的限制?？傮w而言，對于長序列的變形監(jiān)測數(shù)據(jù)粗差的定位與剔除存在理論和應(yīng)用上的諸多困難。

穩(wěn)健估計的目的是挖掘利用有效數(shù)據(jù)，限制使用有用數(shù)據(jù)和清除有害數(shù)據(jù)。自1953年G.E.P.BOX首提穩(wěn)健性（Robustness）的概念以來，國外學(xué)者Huber、Hampel、Rousseeuw等［5-8］對參數(shù)的穩(wěn)健估計進(jìn)行研究，成果頗豐。國內(nèi)學(xué)者周江文、李德仁、歐吉坤、楊元喜等［9-10］也進(jìn)行了大量卓有成效的研究。穩(wěn)健估計的方法最早由Huber引入回歸分析（1981年），目前在水質(zhì)模型［11］、地質(zhì)統(tǒng)計學(xué)［12］等多領(lǐng)域得到了廣泛應(yīng)用。

回歸分析是大壩變形監(jiān)測統(tǒng)計模型的核心，粗差的存在導(dǎo)致回歸參數(shù)不確定，造成統(tǒng)計模型一定程度的失真［13］。在自動化監(jiān)測不斷普及的背景下，如何對有粗差的數(shù)量巨大的數(shù)據(jù)序列進(jìn)行行之有效的建模尤為重要。本文通過對某RCC主壩變形監(jiān)測資料的粗差進(jìn)行處理和建模分析，提出以常規(guī)模型殘差大小進(jìn)行粗差定位的方法，并探討了殘差的分布對于模型參數(shù)的影響、粗差密度與較佳穩(wěn)健權(quán)函數(shù)選取的關(guān)系。

2 大壩安全監(jiān)測統(tǒng)計模型原理

2.1 常規(guī)統(tǒng)計模型原理根據(jù)成因，可將壩體及壩基變形δ分為由水壓分量δH、溫度分量δT和時效分量即：

在壩體已運行多年，溫度場基本趨于穩(wěn)定的條件下，根據(jù)壩工理論和數(shù)學(xué)力學(xué)原理，常規(guī)統(tǒng)計模型可表示為［14］：

式中：Hu、Hu0分別為監(jiān)測日、始測日所對應(yīng)的上游水頭；t為監(jiān)測日到起始監(jiān)測日的累計天數(shù)；t0為建模資料系列第一個監(jiān)測日到始測日的累計天數(shù)；θ等于t除以100，θ0等于t0除以100；其余均為回歸系數(shù)。

基于多元逐步回歸原理，常規(guī)統(tǒng)計模型可表示為：

其中

式中：Y為監(jiān)測值向量；β為待估參數(shù)向量；X為結(jié)構(gòu)矩陣；ε為n維隨機誤差向量，其期望為n維零向量；方差協(xié)方差陣為σ2E（E為n階單位矩陣）。

為了求解式（3）中的參數(shù)β，采用最小二乘法（令結(jié)構(gòu)矩陣X為滿秩矩陣）即可求得β的估計量：

2.2 穩(wěn)健統(tǒng)計模型原理在現(xiàn)階段大壩安全監(jiān)測大力推行自動化系統(tǒng)的過程中，由于無法人工甄別，導(dǎo)致粗差在監(jiān)測序列中普遍存在，因此常使統(tǒng)計模型的參數(shù)求解結(jié)果失真［15-17］。

穩(wěn)健統(tǒng)計模型的建立基于穩(wěn)健估計原理，穩(wěn)健的含義就是在粗差出現(xiàn)（或某一時段存在系統(tǒng)誤差）的情況下使參數(shù)的估計結(jié)果不致失真。與多元回歸不同，穩(wěn)健回歸不追求絕對意義上的最優(yōu)（殘差平方和最?。┓匠蹋塾诳共钜饬x下的最優(yōu)或接近最優(yōu)，即追求回歸參數(shù)的可靠性。二者根本區(qū)別是穩(wěn)健統(tǒng)計模型建立在符合于監(jiān)測數(shù)據(jù)的實際分布模式上，而常規(guī)模型建立在某種理想的分布模式上。

綜上所述，穩(wěn)健統(tǒng)計建模旨在解決常規(guī)統(tǒng)計模型不具有抵抗粗大誤差的缺陷，以適應(yīng)長序列大樣本建模需要，同時也是自動化監(jiān)測系統(tǒng)的進(jìn)行粗差處理及統(tǒng)計建模的解決方案。

穩(wěn)健估計常分三大類型，即M估計、L估計和R估計。其中，M估計又稱極大似然估計，在測量界應(yīng)用成熟廣泛［9-10］。

M估計的準(zhǔn)則為：

ρ稱極值函數(shù)（增長速率小于殘差平方和對未知參數(shù)β求導(dǎo)，并令其為零可得：

考慮到多元的殘差方程是X的第i行向量），有：

令權(quán)函數(shù)則式（7）變?yōu)檫@與最小二乘法估計的方程形式一致，僅是用權(quán)陣代替了觀測矩陣P。值得注意的是與最小二乘所給出的先驗定值權(quán)不同的是，M穩(wěn)健估計的權(quán)函數(shù)是殘差的函數(shù)，計算是未知的，只能通過給其賦予一定的初值，采用迭代的方法估計未知參數(shù)。

由于極值函數(shù)ρ的選取不同，可構(gòu)成的權(quán)函數(shù)也就不同，故稱選權(quán)迭代法。其誤差方程與權(quán)函數(shù)分別為常用的權(quán)函數(shù)迭代方法有Huber法、Andrews法、IGG方案等。

計算流程如圖1所示。

本文穩(wěn)健權(quán)函數(shù)選擇采用Huber法、Andrews法，u表示標(biāo)準(zhǔn)化的殘差指標(biāo)（ui= νi/S ′）。

Huber法：

式中，c取1.0和0.7。

Andrews法：

式中，c取1.0和0.75。

圖1 基于選權(quán)迭代的穩(wěn)健回歸參數(shù)計算步驟

由于常規(guī)統(tǒng)計模型剩余標(biāo)準(zhǔn)差S′受粗差離群幅度及其密度的影響，其值應(yīng)大于無粗差時模型的剩余標(biāo)準(zhǔn)差S。以Huber法為例，當(dāng)滿足條件′時，對殘差數(shù)據(jù)點做降權(quán)處理。按傳統(tǒng)的粗差剔除“2S準(zhǔn)則”或“3S準(zhǔn)則”進(jìn)行粗差剔除，c值應(yīng)取比2小，原因主要有以下兩點：（1）由于S ′＞S ，故2S≤2S′，3S≤3S′；（2）當(dāng)νi≥2S′，作降權(quán)處理，并未剔除。

因此在Huber法中c取1.0和0.7，Andrews法中c取值同理。

3 粗差影響的評價

根據(jù)逐步回歸的理論可知，常規(guī)統(tǒng)計模型參數(shù)的求解基于最小二乘原理（殘差平方和最?。?，而實際上模型的解析力一定程度還與殘差的分布有關(guān)。由于數(shù)據(jù)的采集過程實際上屬于等權(quán)觀測，因此在模型解析力較高的情況下，殘差應(yīng)該大致服從正態(tài)分布。

殘差的正態(tài)性的檢驗方法常用的有：夏皮羅-威爾克（Shapiro-Wilk）法、Pearson法、柯爾莫哥洛夫（Kolmogorov-Smirov）法、Jarque-Bera檢驗及達(dá)戈斯提諾（D′Agostino）法等［18-20］。本文算例均采用Pearson卡方擬合優(yōu)度檢驗法，檢驗方法的步驟為：

（1）將總體X的取值范圍分成k個互不重迭的小區(qū)間，記作A1，A2，…，Ak。

（2）把落入第i個小區(qū)間Ai的樣本值的個數(shù)記作fi，稱為實測頻數(shù)。所有實測頻數(shù)之和f1+f2+…+fk等于樣本容量n。

（3）選定根據(jù)所假設(shè)的理論分布，可以算出總體X的值落入每個Ai的概率pi，于是npi就是落入Ai的樣本值的理論頻數(shù)。

（4）引進(jìn)如下統(tǒng)計量表示經(jīng)驗分布與理論分布之間的差異：

Pearson證明了如下定理，若理論分布F（x）已經(jīng)完全給定，那么當(dāng)n→∞時，式（8）中統(tǒng)計量的分布漸近自由度為k-1的χ2分布。

文中擬合優(yōu)度檢驗顯著性水平α取0.10。

4 工程實例應(yīng)用研究

某工程碾壓混凝土重力壩壩長1489 m，主壩最大壩高121.5 m，電站總裝機容量為140MW。該壩布置有正、倒垂線及引張線等監(jiān)測設(shè)備監(jiān)測壩體變形，壩基傾斜監(jiān)測主要由布設(shè)在壩基橫向廊道的靜力水準(zhǔn)系統(tǒng)完成。

4.1 殘差分布正態(tài)性較差的情況如圖2所示，由于受縱向廊道內(nèi)各類機械施工影響（如混凝土取芯、鉆孔引起的銦鋼絲顫動等）以及觀測間密閉狀況、鋼護(hù)管竄風(fēng)等干擾因素的存在，正垂線上各測點上下游向位移在自動化監(jiān)測系統(tǒng)數(shù)據(jù)采集過程中出現(xiàn)了大量的粗差。

大量粗差的存在給數(shù)據(jù)分析工作帶來了阻礙，常規(guī)的時空分析和特征值分析（極值點常為粗差）已不能準(zhǔn)確反映大壩變形的真實狀況。因此要挖掘出數(shù)據(jù)序列中的有效信息，必須進(jìn)行有效的統(tǒng)計建模。

如圖3為42#壩段EL706.5 m高程PL6-1測點按周頻次取值所建常規(guī)統(tǒng)計模型過程線，可以看出模型整體擬合精度較高（模型復(fù)相關(guān)系數(shù)R=0.9495）。但由于受到粗差的影響，殘差在粗差密度較大的時段明顯增大。圖4為殘差的頻次分布圖，按Pearson卡方擬合優(yōu)度檢驗法，P值為0.0001。因此殘差顯著地不符合正態(tài)分布，對常規(guī)模型擬合曲線存在較大影響。

圖3 PL6-1測點常規(guī)統(tǒng)計模型過程線

圖4 PL6-1測點常規(guī)統(tǒng)計模型殘差頻次圖

表1為該測點常規(guī)統(tǒng)計模型與穩(wěn)健統(tǒng)計模型參數(shù)對比表。其中不含粗差的常規(guī)模型參數(shù)是按“vi＞1.5S′準(zhǔn)則”對粗差點進(jìn)行剔除后再進(jìn)行多元回歸分析得到。由表可知，剔除粗差后，常規(guī)模型各項系數(shù)變化較大，溫度因子系數(shù)b11最大變幅達(dá)23.9%。而不含粗差常規(guī)模型與Huber法（c=1.0）的參數(shù)求解結(jié)果基本一致。表中系數(shù)變幅按下式計算：

表1 PL6-1常規(guī)模型與穩(wěn)健模型參數(shù)

如圖5所示為常規(guī)模型（含粗差、不含粗差）與穩(wěn)健模型Huber（c=1.0）擬合值過程線，可以看出三者擬合過程線基本重合，說明穩(wěn)健統(tǒng)計模型結(jié)果在追求可靠性分析結(jié)果的同時，并未降低擬合效果。而如圖6所示，三者水壓分量過程線存在較為明顯的差異，其中穩(wěn)健模型Huber（c=1.0）水壓分量與粗差剔除的常規(guī)模型水壓分量較為接近。2016年6月6日，二者水壓分量值分別為4.34 mm、4.41 mm，而未剔除粗差常規(guī)模型水壓分量為3.94 mm，其分量評價存在較大偏差。

圖5 PL6-1統(tǒng)計模型擬合過程線

圖6 PL6-1測統(tǒng)計模型水壓分量過程線

此外，還可看到以常規(guī)統(tǒng)計模型的剩余標(biāo)準(zhǔn)差S′作為尺度，并結(jié)合殘差值vi的大小來具體判定粗差點，較之常規(guī)的粗差剔除方法具有以下幾個優(yōu)點：（1）統(tǒng)計模型顧及了壩體變形量與庫水位、溫度、時效等因子的相互影響關(guān)系，采用該法進(jìn)行粗差的定位，不再如“3σ準(zhǔn)則”、格布羅斯及狄克遜準(zhǔn)則那樣完全基于一定的統(tǒng)計分布或根據(jù)一定的顯著性水平孤立地判別變形效應(yīng)量是否為粗差。（2）對于擬合效果較好的統(tǒng)計模型，殘差過程線能夠直觀、準(zhǔn)確地反映出粗差點的位置，并估算粗差密度。（3）可用于大樣本長序列的粗差定位。

4.2 殘差分布正態(tài)性較好的情況如圖7所示，壩基傾斜采用靜力水準(zhǔn)監(jiān)測，儀器抗干擾能力較強，測值粗差數(shù)量相對較小。結(jié)合表2可知，28#壩段（壩下0+007.0至壩下0+058.0）壩基傾斜量常規(guī)統(tǒng)計模型擬合效果良好，模型各分量能夠反映出壩基傾斜量發(fā)展過程的良好規(guī)律。如圖8所示，該壩段傾斜量常規(guī)統(tǒng)計模型殘差基本符合正態(tài)分布，Pearson擬合優(yōu)度檢驗P值達(dá)0.81。

如表3所示，由于殘差分布接近正態(tài)分布，故剔除邊緣粗差尖點后，模型參數(shù)變化相對較小，與穩(wěn)健統(tǒng)計模型Huber（c=1.0）、Andrews（c=1.0）參數(shù)較為接近。而且若采用粗差限制較為嚴(yán)格的權(quán)函數(shù)Huber（c=0.7）和Andrews（c=0.75），將使統(tǒng)計參數(shù)一定程度上失真。

表2 28#壩段壩基傾斜常規(guī)統(tǒng)計模型擬合效果及Pearson檢驗結(jié)果

因此，說明了在進(jìn)行穩(wěn)健統(tǒng)計模型建模前，應(yīng)對常規(guī)統(tǒng)計模型的殘差分布有一定的把握：對于近于正態(tài)分布的殘差序列，常規(guī)統(tǒng)計模型參數(shù)已較為穩(wěn)定，不建議使用穩(wěn)健統(tǒng)計模型；對于與正態(tài)分布相去甚遠(yuǎn)的殘差序列，應(yīng)根據(jù)殘差大小，判定粗差密度，進(jìn)而選擇相應(yīng)的穩(wěn)健統(tǒng)計權(quán)函數(shù)進(jìn)行建模分析。

圖7 28#壩段傾斜量(壩下0+007.0至壩下0+058.0)統(tǒng)計模型過程線

圖8 常規(guī)統(tǒng)計模型殘差頻次分布圖

表3 28#壩段壩基傾斜量常規(guī)模型與穩(wěn)健模型參數(shù)

4.3 權(quán)函數(shù)的選取如前所述，常規(guī)統(tǒng)計模型的剩余標(biāo)準(zhǔn)差S′是有效定位粗差的一個尺度（如按“νi＞1.5S′準(zhǔn)則”或“νi＞2S′準(zhǔn)則”對粗差點進(jìn)行剔除），因此可根據(jù)殘差νi的大小得到粗差的密度（粗差數(shù)/樣本數(shù)）。

而粗差密度可作為權(quán)函數(shù)選取的一個導(dǎo)向，Hampel認(rèn)為，實際測量的數(shù)據(jù)中粗差出現(xiàn)的概率達(dá)10%是屬正常的［18］，如表5所示，根據(jù)某工程正垂線測點統(tǒng)計結(jié)果來看，在干擾條件較強的情況下，粗差密度可達(dá)15.58%。表中粗差判定根據(jù)“νi＞1.5S′準(zhǔn)則”，并由此得到粗差密度。而較佳穩(wěn)健權(quán)函數(shù)的統(tǒng)計模型參數(shù)與剔除粗差后的常規(guī)統(tǒng)計模型參數(shù)最為接近，定義Am為模型系數(shù)平均變幅，按下式計算：

式中：an、bij、cm為穩(wěn)健統(tǒng)計模型回歸參數(shù)；a′n、a′ij、a′m為去粗差后的常規(guī)模型參數(shù)；當(dāng)n取0時，a0與a′0為分別為模型常數(shù)項參數(shù)；N為逐步回歸選入自變量因子個數(shù)；Am最小時對應(yīng)的權(quán)函數(shù)即為較佳穩(wěn)健權(quán)函數(shù)。

此外，表5中Huber及Andrews函數(shù)若未注明c的大小，則為Matlab中默認(rèn)取值，分別為1.345、1.339。

由表5可以看出，粗差密度小于6%時，可選Andrews或Huber作為穩(wěn)健權(quán)函數(shù)；密度介于6%～9%之間時，應(yīng)選取Huber函數(shù)進(jìn)行賦權(quán)迭代計算；粗差密度在10%～11%時，較佳權(quán)函數(shù)為Huber（c=1.0）；密度在12%左右時，宜選取Huber（c=0.7）作權(quán)函數(shù)；密度大于15%時，穩(wěn)健權(quán)函數(shù)應(yīng)選取 Andrews（c=0.75），此時根據(jù)殘差 νi大小給觀測值賦權(quán)有較嚴(yán)格的控制。

表4 PL6-2測點水位相關(guān)系數(shù)矩陣

表中穩(wěn)健模型效果較差兩組數(shù)據(jù)（PL2-2、PL6-2）具有兩個共性的特點：一是常規(guī)統(tǒng)計模型與穩(wěn)健統(tǒng)計模型中的常數(shù)項參數(shù)a0與a′0較大，例如PL2-2測點2016年6月20日位移量為6.03 mm，而其常規(guī)統(tǒng)計模型與穩(wěn)健統(tǒng)計模型常數(shù)項分別為501.52 mm、583.08 mm。由于各分量值遠(yuǎn)遠(yuǎn)小于常數(shù)項，統(tǒng)計模型的解析力弱，造成剔除粗差或穩(wěn)健統(tǒng)計建模后參數(shù)變化較大；二是兩測點均選入了3個水位因子，如表4所示，由于各水位因子間內(nèi)部相關(guān)性較高，存在信息冗余，導(dǎo)致系數(shù)變化相互影響。

表5 某工程正垂測點上下游向水平位移粗差密度與較佳穩(wěn)健權(quán)函數(shù)統(tǒng)計

5 結(jié)論

本文根據(jù)某工程RCC主壩變形監(jiān)測實測資料，對當(dāng)前大壩變形監(jiān)測數(shù)據(jù)序列中出現(xiàn)的粗差處理方法進(jìn)行了探討，并采用穩(wěn)健估計原理進(jìn)行統(tǒng)計建模，得到如下結(jié)論：（1）以常規(guī)統(tǒng)計模型剩余標(biāo)準(zhǔn)差為尺度，根據(jù)殘差的大小來定位粗差的方法可顧及變形量與庫水位、溫度、時效等因子的相互影響關(guān)系，較之傳統(tǒng)的“3σ”、格布羅斯和狄克遜準(zhǔn)則等方法直觀、準(zhǔn)確，適用于大樣本長序列的粗差定位。（2）殘差序列符合正態(tài)分布，常規(guī)統(tǒng)計模型參數(shù)穩(wěn)定，可不采用穩(wěn)健統(tǒng)計模型；若殘差序列與正態(tài)分布相去甚遠(yuǎn)，可根據(jù)粗差密度，選擇相應(yīng)的穩(wěn)健統(tǒng)計權(quán)函數(shù)進(jìn)行建模分析。（3）常數(shù)項參數(shù)過大以及選入多重相關(guān)性較強的因子對統(tǒng)計模型回歸參數(shù)穩(wěn)定不利。（4）實例證明，與常規(guī)統(tǒng)計模型相比，穩(wěn)健統(tǒng)計模型能夠有效地抵抗粗差對于模型參數(shù)求解結(jié)果的影響，并能正確反映和評價水壓、溫度和時效各分量的影響。

參考文獻(xiàn)：

［1］姚宜斌.粗差的定性分析［J］.測繪信息與工程，2002，27（1）：1-3.

［2］楊建潮.測量誤差及粗大誤差的判別與消除［J］.計量與測試技術(shù)，2006，33（11）：4-5.

［3］楊茂興.小樣本容量測量數(shù)據(jù)中粗差的剔除［J］.計量與測試技術(shù)，2005，32（1）：27-28.

［4］陶本藻，姚宜斌.可靠性分析與數(shù)據(jù)探測［J］.武漢大學(xué)學(xué)報（信息科學(xué)版），2002，27（6）：607-610.

［5］HUBER P J.Robust estimation of a location parameter［J］.Annals of Mathematical Statistics，1964，35（1）：73-101.

［6］HUBER P J.Robust Statistics［M］.Wiley-Interscience，2011.

［7］ROUSSEEUW P J，HAMPEL F，RONCHETTI E，et al.Robust Statistics：the Approach Based on Influence Functions［M］.J Wiley，New York，1986.

［8］ROUSSEEUW P J，LEROY A M.Robust Regression and Outlier Detection［M］.Wiley-Interscience，2003.

［9］周江文，黃幼才，楊元喜，等.抗差最小二乘法［M］.武漢：華中理工大學(xué)出版社，1997.

［10］李德仁.測量平差系統(tǒng)的誤差處理和可靠性理論［M］.武漢：武漢測繪科技大學(xué)出版社，1985.

［11］李黎武，施周.隨機噪聲干擾下水質(zhì)模型參數(shù)的魯棒估計方法［J］.水利學(xué)報，2006，37（6）：687-693.

［12］陳亞新，徐英，魏占民，等.基于穩(wěn)健統(tǒng)計學(xué)的水鹽空間變差函數(shù)逼近方法［J］.水利學(xué)報，2004（9）：44-49.

［13］賈超.穩(wěn)健回歸分析方法在變形監(jiān)測中的應(yīng)用［D］.太原：太原理工大學(xué)，2011.

［14］吳中如.水工建筑安全監(jiān)控理論及應(yīng)用［M］.北京：高等教育出版社，2003.

［15］中國水利水電科學(xué)研究院.廣西右江百色水利樞紐工程安全監(jiān)測資料分析評價報告［R］.北京：中國水利水電科學(xué)研究院，2016.

［16］河海大學(xué).李家峽水電站大壩監(jiān)測資料分析［R］.南京：河海大學(xué)，2008.

［17］北京中水科工程總公司.新疆喀臘塑克水利樞紐工程安全監(jiān)測資料分析及評價報告［R］.北京：北京中水科工程總公司，2016.

［18］D′AGOSTINO R B，BELANGER A，D′AGOSTINO JR R B.A suggestion for using powerful and informative tests of normality［J］.The American Statistician，1990，44（4）：316-322.

［19］BERA A K，JARQUE C M.Efficient tests for normality，homoscedasticity and serial independence of regression residuals：Monte Carlo evidence［J］.Economics Letters，1981，7（4）：313-318.

［20］GHASEMI A，ZAHEDIASL S.Normality tests for statistical analysis：a guide for nonstatisticians［J］.Interna?tional Journal of Endocrinology and Metabolism，2012，10（2）：486-489.

［21］HOAGLIN D C，MOSTELLER F，TUKEY J W，et al.探索性數(shù)據(jù)分析［M］.北京：中囯統(tǒng)計出版社，1998.