谷幅變形全天時自動化監測的精密氣象改正方法

周 綠，張 凱，劉書明，楊澤發，肖亞子

(1.中國電建集團中南勘測設計研究院有限公司，湖南長沙410014；2.中南大學地球科學與信息物理學院，湖南長沙410083)

近些年來，我國一大批高拱壩投入建設運營，有效促進了社會經濟發展。然而，在壩體運營過程中，尤其在蓄水期庫區兩岸山體受水位上升的影響，應力平衡狀態易受擾動，從而形成大范圍水平變形，簡稱“谷幅變形”[1-2]。數據顯示，我國建成的高拱壩如溪洛渡水電站[3-4]、白鶴灘水電站[5]、錦屏一級水電站[6]等在正常蓄水過程中均監測到不同程度的谷幅變形。因此，準確監測高拱壩谷幅變形對保障水電站安全運營有重要意義。

目前，常用測量機器人等全自動化設備監測谷幅變形，其大都基于光電測距原理。然而，由于庫區內氣溫、氣壓、濕度等因素變化復雜，其會對光電測距精度有顯著影響。目前谷幅變形監測大都采用《中、短程光電測距規范》[7]推薦使用的經典氣象改正模型提高監測精度。然而，該方法需假定觀測時氣象條件(比如氣溫、氣壓和濕度)相對平穩。因此，規范要求光電測距作業時間通常選取空氣溫度垂直變化梯度為零的前后1 h內(比如日出后0.5 h或日落前0.5 h)，而且在正午、雷雨、大霧、大風等復雜氣象條件下氣象模型精度較差。

隨著我國水電站壩體高度和蓄水量的不斷增加，對谷幅變形監測的精度和時間分辨率均提出了更高的要求(比如亞毫米精度、全天時監測等)。然而，受觀測成本、天氣條件等因素影響，人工谷幅變形監測不僅工作強度大，且時間分辨率也不高，難以滿足水電站對于谷幅變形監測的時間分辨率要求。欲實現谷幅變形全天時自動監測，除儀器設備外，適用于全天時的測距氣象改正方法顯得尤為重要。目前，國內外學者在自動化氣象改正方面開展了部分探索，并發展了鄰邊比例法、基線神經網絡法、差分數據處理、似差分改正與模型改正等多種算法[8-10]，取得了一定效果。然而，這些方法大都需要以一條已知且長度不變的精確基線作為參考依據。由于水電站蓄水影響范圍較大，且微氣象學變化顯著，因此，在谷幅觀測線附近選取一條長度不變，且具有代表性的精確基線較為困難。

為此，本文通過充分挖掘谷幅變形自動化監測數據內部的演化趨勢，通過引入現代測量平差理論發展一種谷幅變形全天時自動化監測精密氣象改正方法。該方法首先采用傳統氣象改正模型對谷幅測距數據進行初步改正；隨后，基于Fourier數進一步擬合并去除全天時測距趨勢性偏差，并采用穩健Kalman濾波方法削弱隨機誤差；最后，基于非參數擬合的加權最小二乘方法計算全天時谷幅變形值。

1 研究方法

1.1 基于傳統氣象改正模型的谷幅測距初步改正

根據文獻[7]，經典氣象改正模型可表示為

dD≈(n0-n)D′=(N′g-N″g)×10-6×D′

(1)

式中，D′為待氣象改正的原始測距值，m；N′g為儀器標稱群折射率；N″g為實際作業大氣條件下的群折射率；dD為氣象修正值；n為測量時氣象條件下的實際群折射率；n0為儀器氣象參考點的群折射率。

Barrell和Sears所描述的標準大氣條件(0 ℃，1 013.25 hPa，包含0.03%的二氧化碳的干燥空氣)下可見光和近紅外調制輻射的群折射率Ngo(標準大氣)為

(2)

式中，電磁波波長為微米級時，A=287.604，B=1.628 8和C=0.013 6；λ為真空中光波的有效波長。代入計算得

(3)

波長在560～900 μm間，公式(3)的精度能精確到±0.1 ppm。為了適應測距時的實際大氣條件，Barrell和Sears對公式做如下修改，為

(4)

式中，Ng為一般大氣的群折射率，ppm；Ngo為標準大氣的群折射率，ppm；T為絕對溫度，K；P為大氣壓，hPa；e為水蒸氣壓力，hPa；Q為常數，取0.269 6；V為常數，取11.27。

由于濕度觀測相對容易，因此，大多數谷幅變形自動化監測直接觀測相對濕度h，與水蒸氣壓力e的關系式為

(5)

式中，E為飽和水蒸氣壓，可通過實測濕球溫度tw(單位為℃)估計，即

(6)

1.2 基于Fourier級數的全天時測距趨勢誤差估計

通過分析傳統大氣改正后的距離值，發現其存在周期性系統性誤差。出現該現象的原因可能是現有氣象誤差改正模型要求氣象條件相對平穩，而該要求在全天時觀測時難以準確滿足導致的。由于本文采用加權最小二乘估計全天時谷幅變形，該方法要求觀測數據僅含隨機誤差，因此，需對傳統氣象模型改正后的谷幅距離觀測值進行趨勢誤差去除。考慮到Fourier級數具有較好的趨勢項描述能力[11-12]，因此，本文采用該函數擬合并去除全天時測距趨勢誤差。

Fourier擬合算法可采用一對正交函數通過對Fourier級數的變換，不斷逼近離散數據的平滑處理。通過Fourier級數擬合方法估計全天時谷幅距離值(經過傳統氣象模型校正之后)平均偏差誤差趨勢，發現其可用如下公式表示

(7)

式中，f(x)為周期的連續函數，周期為T；a0、ω為常系數；n為傅里葉展開級數；角頻率ω=2π/T；ai、bi(1≤i≤n)為Fourier系數。

1.3 基于穩健Kalman濾波的全天時谷幅測距隨機誤差削弱

經過趨勢項偏差去除之后，本文采用穩健Kalman濾波方法削弱全天時谷幅測距觀測值的隨機誤差。Kalman濾波[13]算法可根據上一次的估計值與本次的測量值，結合二者的誤差計算出的增益系數，從而計算本次的估計值，即

(8)

(9)

傳統Kalman濾波算法是以可靠的函數模型和隨機模型為前提，當其面對狀態模型不準確或觀測誤差中含有粗差的場景時，其估計結果可能出現錯誤。為了克服該局限，本文采用穩健Kalman濾波算法[14]削弱谷幅測距的隨機誤差。穩健Kalman濾波算法將抗差理論與傳統Kalman濾波方法結合，其重點在于權矩陣的設計。本文采用目前較為常用的IGG-I加權方案，即

(10)

(11)

式中，PXk為預報值的權矩陣，由其協方差矩陣DXk求逆獲得。

1.4 聯合非參數擬合與加權最小二乘的全天時谷幅變形監測

經過穩健Kalman濾波后，谷幅測距觀測值的隨機誤差得到了較好的削弱。然而，此時的測距觀測值精度可能仍難以滿足亞毫米變形監測精度的要求，需對其進行進一步精化。考慮到全天時高時間分辨率測距(理論上可達到近實時)可提供大量的多余觀測，因此，本節將引入加權最小二乘算法，進一步提高谷幅測距精度。

具體地，假定短時間內(比如數小時)谷幅變形值為0，且該時段內有m個距離觀測值，據此可建立該時段內谷幅真實距離S與距離觀測值L之間的函數模型，為

Am×1S1×1=Lm×1+εm×1

(12)

(13)

式中，P為觀測值的權矩陣。

本文采用非參數擬合確定。具體地，首先以全天候測距數據中兩岸氣象條件相對穩定時刻測距觀測值的中值作為參考；然后，將Kalman濾波后的距離觀測值與參考值作差，獲取各時刻觀測值與參考值的差異圖(簡稱“離差”)；之后，采用非參數擬合確定全天時觀測值時序離差，并利用時序離差的倒數為全天時最小二乘權重；最后，根據選取的時段提取相應權值，即利用公式(13)求取該時段內的距離估計值。

在獲得全天候谷幅測距數據后，即可利用相對于首次觀測的距離差作為谷幅變形值。從上述描述中可以發現，本文方法僅依賴于谷幅測距同步的氣溫、氣壓和濕度數據，并通過自動挖掘全天候觀測數據自身蘊含的信息削弱測距誤差，為谷幅變形自動化監測的氣象改正提供了全新的模型。

2 研究區域

本文采用錦屏一級水電站的谷幅變形自動化測線驗證改正模型的可靠性。錦屏一級水電站擁有目前世界最高的蓄水壩體(305 m)，谷幅變形對于壩體安全的影響巨大。為此，水電站安全監測運維及實施相關單位在水電站壩體下游布設了一條自動化觀測線。測距儀采用迪馬斯高精度激光測距，同時同步觀測兩岸的氣溫、氣壓和濕度，采樣頻次為1次/小時。此外，在自動化測線鄰近位置布設了人工高精度觀測線，為本文模型與方法精度驗證提供數據保障。

3 實驗與結果分析

根據前述方法，首先采用了傳統氣象改正模型初步改正全天時自動化測線觀測的谷幅距離值，如圖1所示。圖1展示了該測線2018年12月～2019年11月期間的原始測距值和經過傳統模型校正后的距離值。從圖1中可以看出，由于谷幅實測氣象條件與光電測距參考氣象條件的差異，經過傳統氣象校正模型改正后的距離值與原始觀測值存在顯著差異。此外，初步改正后的距離值仍存在較大的隨機誤差，且有一定的趨勢誤差，需對其進行進一步校正。

圖1 原始測距數據與傳統氣象改正數據對比

為了更清晰地展示全天候谷幅測距值的趨勢性誤差，本文首先選取每天氣象條件相對平穩時刻的測距值作為參考，并計算24 h距離值與參考值的偏差。圖2為錦屏一級水電站在2018年12月～2019年11月期間的小時級測距偏差分布圖。從圖中可以看出，全天測距偏差存在趨勢項誤差，且年際尺度的最大振幅超過1 mm，該誤差可能會嚴重削弱全天時形變的監測精度。因此，削弱全天時谷幅測距的趨勢項誤差至關重要。為此，本文選用了24 h測距離偏差作為樣本，采用Fourier級數擬合并去除谷幅測距趨勢性誤差。圖2中曲線為2018年12月～2019年11月期間的日均值離差傅里葉擬合曲線。

圖2 傳統氣象改正引起的距離日均值偏差

經過趨勢項誤差削弱和穩健Kalman濾波后的測距數據對比圖，如圖3所示。圖3灰色曲線表示經過趨勢項誤差校正后的全天時谷幅距離觀測值。將圖1中的傳統氣象改正與圖3中去趨勢項誤差后的測距數據作對比，如圖4所示。從圖4可以看出，相較于傳統氣象模型校正后的觀測值(圖1中位于上部的曲線)，谷幅距離變化趨勢更清晰，說明本文提出的Fourier級數擬合算法是有效的。然而，從圖3中仍可看出，全天時谷幅距離觀測值仍包含明顯的隨機誤差，且存在個別差異較大的觀測值(稱為粗差)。

圖3 趨勢項誤差削弱和穩健Kalman濾波后的測距數據對比

圖4 傳統氣象改正與去除趨勢線誤差后的測距數據對比

為了更進一步削弱觀測隨機誤差的影響，本文采用了穩健Kalman濾波對去除趨勢誤差后的測距值進行了改正，其結果如圖3紅色曲線所示。從圖中可以看出，除個別觀測點外，大部分觀測值的隨機誤差得到了較好的抑制。與此同時，經過系統性偏差削弱后，谷幅距離變化趨勢也逐漸清晰。

最后，本文采用第1.4節中描述的非參數擬合方法確定全天候24 h距離值的權重，并采用加權最小二乘方法估計了該自動化谷幅測線的距離觀測值，如圖5所示。從圖5可以看出，經過加權最小二乘平差后谷幅變形趨勢信號較為清晰，即2018年12月開始，該測線的水平距離逐漸增大，表明該地區的谷幅呈現擴張趨勢，2019年2月，谷幅迅速收縮，最大收縮幅度達到了6 mm。之后，該谷幅基本保持穩定，到2019年6月有突然的擴張和收縮，隨后平穩直至2019年9月，在此之后，該谷幅再一次快速擴張，直至2019年10月迅速收縮，之后便再一次趨于穩定。

圖5 谷幅變形自動化測線的加權最小二乘解時序

自動化谷幅變形值與人工測線對比，如圖6所示。圖6藍色線為采用加權最小二乘方法估計后該自動化谷幅測線的距離改正值，通過對比和人工谷幅測距結果(圖中DY16紅色測線)發現，本文方法計算的自動化谷幅測距精度約為0.9 mm，基本能夠滿足谷幅變形監測對精度的要求。

圖6 計算的自動化谷幅變形值(藍色實線)與人工測線(紅色菱形)對比

4 結 論

本文提出了一種面向谷幅變形全天時自動化監測的精密氣象改正方法，并在錦屏一級水電站中進行了方法測試。實驗結果表明：

(1)經過傳統氣象模型改正后，全天時谷幅距離觀測值存在趨勢項誤差，且本文提出的Fourier級數擬合算法能夠有效地削弱該趨勢項誤差的影響。

(2)經過趨勢項誤差改正后，全天時谷幅距離觀測值仍存在較為顯著的隨機誤差和粗差，而本文提出的穩健Kalman濾波算法和基于非參數擬合的加權最小二乘算法能較好地抑制隨機誤差和粗差的影響，提高谷幅變形觀測精度。

(3)通過與實測數據對比分析，本文提出的算法精度可達0.9 mm，基本符合谷幅變形監測需求。