基于BDS及InSAR的龍羊峽水電站庫區地表形變研究

張永滎，張生鵬,2，貢寶扎西，李水平，李亞剛，省天琛

（ 1. 青海省基礎測繪院, 西寧 810001；2. 青海師范大學 地理科學學院, 西寧 810008；3. 合肥工業大學 土木與水利工程學院, 合肥 230009 ）

0 引 言

水電站是通過調控水資源時空分布成為蓄水防旱的重要樞紐工程，其運行直接關系到庫區及下游居民的重大生命財產安全[1]. 中國已建成各類水庫近10萬座，但是多數水庫由于建造年代久遠，受限于當時的條件，病險情況普遍存在，最嚴重時其數量占比達到30%. 形變監測是水庫安全運營評估的基礎，目前水庫形變監測主要采用傳統全站儀、水準測量等手段在重點部位進行監測，傳統大地測量的工作特點決定其監測結果空間采樣率不足及觀測成果具有周期性，導致在大范圍監測結果中無法及時有效地發現新的形變區域，造成重大安全隱患的遺漏[2-3].

全球衛星導航系統(GNSS)和合成孔徑雷達干涉測量(InSAR)技術作為大地測量新技術，為水庫形變監測帶來了新的方案. 王騰等[4]采用InSAR永久散射體(PS)技術主要對三峽水庫大壩進行了穩定性分析，得出利用此種技術可以獲得遠高于傳統手段監測點的結論. 裴媛媛等[5]采用PS技術對上海部分堤壩進行了穩定性分析，得出此方法可以得到毫米級精度，更好地用在人工線狀地物的結論. 張明晶[6]對于橫泉水庫GNSS自動化監測進行研究，得出此類系統可以提高監測的實時性和可靠性的結論. MILILLO等[7]利用InSAR時序方法研究了伊拉克摩蘇爾大壩2004—2016失穩過程，找出了各個時段詳細的形變趨勢. 上述研究都是單獨利用GNSS或InSAR進行監測，其中InSAR基本采用PS技術，聚焦點在水庫人工建筑，GNSS監測未單獨驗證北斗衛星導航系統(BDS)的可用性，也未融合現有連續運行參考站系統(CORS)技術的優勢，其實時定位精度停留在厘米級. 龍羊峽水電站作為黃河上游第一座大型梯級電站，對下游水電站的水量調節起主導作用，由于其特殊的地形、地貌、巖土體類型及氣象因素的誘發，成為大型滑坡的高發育區，并且迄今為止還完全采用傳統測量方式進行形變監測[8]. 本文綜合考慮GNSS的高時間分辨率和InSAR的高空間分辨率，結合青海省連續運行參考站系統(QH-CORS)數據成果采用InSAR小基線集技術對龍羊峽庫區進行了地表形變監測分析，驗證并采用InSAR技術進行大面積形變區排查，基于CORS成果提升BDS實時動態(RTK)差分精度，進行重點區域形變RTK差分監測方法的可行性及精度.

1 研究方法

1.1 BDS RTK差分監測

基于CORS成果進行BDS實時動態差分計算的步驟如下：

步驟一：將基準站的觀測值和實測氣象元素按特定頻率匯集到數據中心，建立新的對流層延遲改正模型進行RTK定位；

步驟二：選擇經驗模型，計算所有站前一段時間每天的全球經驗模型和實測氣象元素的天頂對流層延遲(ZTD)；

步驟三：建立高程、基于全球經驗模型和基于實測氣象元素對流層延遲的模型之間的關系，進行新的修正模型；

步驟四：利用改正后的基準站和監測站ZTD模型得到基準站觀測值的改正量；

步驟五：處理后的基準站觀測值播發給監測站，監測站進行實時差分解算出三維(3D)坐標.

1.2 短基線子集技術

短基線子集(SBAS)技術于2002被提出，是當前較主流的InSAR時間序列分析方法[9]. 該方法依據時空基線的條件，自由組合基線較短的影像對，組合成若干個短基線集合，再利用奇異值分解方法，將組合成的所有短基線集聯合起來求解，從而得到研究區域內整個觀測時間的形變時間序列和平均沉降速率.SBAS的基本步驟如下.

采集同一區域按照時間順序t0， ···，tN排列的N+1幅雷達影像，選取其中一幅影像作為主影像進行配準，則形成M幅差分干涉圖.

對于從影像tA和主影像tB(tB＞tA) 時刻獲取的合成孔徑雷達(SAR)影像生成的第j幅差分干涉圖，在不考慮大氣延遲相位、殘余地形相位和噪聲位的情況下，其中任意像元 (x,y) 的干涉相位為

將式(1)中相位采用兩個獲取時間之間平均相位速度和時間的乘積的方式表示，則有

式中： Φ 為殘余地形相位，得第j幅干涉圖的相位值為

利用矩陣形式表示為

B是M×N矩陣，對B進行奇異值求解，獲取最小范圍下的視線向的相位平均速率，對各時間段的沉降速率進行時間域上的積分，得到視線向上時間序列的形變量，通過三角函數，將它們轉為垂直向上監測時間段內的累計沉降量和形變速率[10-12].

2 研究區域概況

龍羊峽水電站位于青海省共和縣與貴德縣之間的黃河干流上，黃河穿越其中，右岸是險峻的茶納山，左岸是莽原. 屬高原干旱-半干旱大陸性氣候，庫區植被稀少，水位基本每年7月開始上升，11月份達到最高值，然后開始逐漸回落. 水庫右岸發育北東東向活動斷裂和北西向斷裂，巖土體類型為雙層結構巖土體，下部軟弱的層狀半固結巖土體，上部為松散沖洪積一般粘性土[8，13].

研究區及GNSS監測點布設位置如圖1所示，其中紅色框內為龍羊峽水庫庫區，紅色三角點為QHCORS龍羊峽站及GNSS監測點.

圖1 研究區及BDS站點布設圖

3 成果獲取及分析

3.1 GNSS數據處理

在龍羊峽大壩下方布設了5個GNSS監測點，將數據通過無線網絡傳實時輸到QH-CORS核心服務器上. 選取QH-CORS龍羊峽站為基準站，按照1.1節中的研究方法對采集到的BDS數據進行RTK差分定位解算. 其步驟一中特定頻率為1 Hz，步驟二中全球對流層模型為GPT2w，利用薩斯塔莫寧(Saastamoinen)模型計算龍羊峽站每天各時刻的ZTD. 圖2為流程及具體方法解算出2019年BDS數據RTK差分定位形變結果.

圖2 BDS監測點實時動態解算流程圖

利用GNSS定位GAMIT軟件解算收集的與SAR和BDS數據時間相對應的13天的靜態數據[14]，為了最大程度的削弱水域對解算成果的影響，利用商業軟件預算數據得到高精度的先驗坐標，進行潮汐改正.采用目前模擬對流層狀況精度最高的模型進行大氣改正，其主要設置為：

衛星軌道：國際GNSS服務(IGS)快速精密星歷；

衛星截止高度角：15°；

數據采樣間隔：10 s；

坐標約束：CORS水平方向給予5 cm，垂直方向給予10 cm的約束，監測點給予10 m的約束；

潮汐改正：otl_FES2004.grid；

氣象模型：VMF1[15].

基線重復性可以反映出基線解內部精度，是評定數據處理質量的重要指標. 表1顯示基線重復性在各方向的固定誤差均小于2 mm，相對精度都達到10-7水平，滿足形變監測毫米級精度的要求[14].

表1 基線向量重復性統計表

3.2 SAR數據處理

選擇覆蓋研究區域的2017-10-15—2019-10-17獲取的25景Sentinel-1SAR影像，C波段，分辨率5 m×20 m，精密軌道文件，30 m分辨率的數字高程模型數據.

采用SBAS技術解譯獲取的全部影像，構建了一個所有時相配對、連接均勻的SBAS數據集合，如圖3所示，主要參數：

圖3 解纏像對連接圖

極化方式：VV；

臨界基線最大百分比：3%；

最大時間距離：60；

多視視數：1∶5；

解纏方式：3D，Delaunay MCF；

相關系數閾值：0.35；

大氣低通：1 200；

大氣高通：365；

形變投影：垂直.

采用的解纏方法可以很好地處理兩個相干性高的孤立區域，更好地抑制水汽對解譯的影響.

圖4通過解譯獲得了龍羊峽庫區2017—2019年在垂直方向上的形變速率，區域內大面積沉降區主要在入庫區、查西、旅游公司區域. 最嚴重沉降區域位于入庫區，最大沉降速率52.48 mm/a；最大抬升點在查西，最大抬升速率43.60 mm/a. 其全部分布于庫區的南岸，與其他學者研究龍羊峽滑坡集中發育在南岸的論述相符合[16-17].

3.3 成果分析

3.3.1 地表形變原因推斷

入庫區和旅游公司所在區域巖土主要以粘土和湖積粉砂為主. 黃河在龍羊峽內直流而下進入入庫區，在其沖擊作用下將部分泥沙沖入庫區，并在黃河下切影響下，沖刷南岸造成沉降現象. 旅游公司地勢平緩，在庫區水的浸潤、浪蝕作用下形成一個沉降量較小的沉降形變. 查西在黃河下切的影響下形成臨空面，馬蹄形地勢造成沉降和抬升的同時存在水位下降，山體內地下水位下降幅度相對滯后，地下水向庫區滲流，產生動水壓力，增大山體的下滑力[8，13].

3.3.2 成果精度驗證

所得成果的精度采用外符合精度評定. 外符合精度是以外部提供的參考值為比對基準，主要反映觀測值與參考值之間的偏差程度. 外符合精度計算公式如下：

式中：m外是 外符合精度； Δ 是參考值減去觀測值后的殘差；n為監測點參與計算數據期數[15].

將靜態解算成果與BDS RTK監測成果相減得到的殘差值代入式(5)，得到各個點的外符合精度值.然后對各點外符合精度值求平均值得到監測精度. 同理，得到InSAR監測成果的精度.

由表2～3可知，BDS RTK差分監測成果的精度為7.1 mm，InSAR監測成果的精度為4.4 mm. 以GNSS靜態解算成果為參考值，兩種技術的監測精度都達到毫米級.

表3 InSAR監測成果精度計算表 mm

4 結論和討論

通過InSAR提取高空間分辨率的地表形變信息，GNSS獲取高時間分辨率的空間3D信息. 近年來國家對BDS系統的建設投入越來越大，而CORS的迅猛發展推動了其廣闊的應用前景. 本文以龍羊峽水庫為研究區域，通過兩種技術獲取了整個庫區垂直方向的地表形變信息，分析結果表明：在黃河影響下，2017—2019年龍羊峽庫區存在大面積的地表形變區，最嚴重沉降區域位于入庫區，最大沉降速率為52.48 mm/a；最大抬升點在查西，最大抬升速率為43.60 mm/a. 以GNSS靜態解算成果為參考值進行驗證，BDS RTK差分監測成果的精度為7.1 mm，In-SAR監測成果的精度為4.4 mm，都滿足形變監測的精度要求. 采用InSAR進行大面積形變區篩選，再結合CORS利用BDS進行重點區域實時監測的模式對于水庫全方位形變監測是可行的. 由于InSAR對平面信息不敏感，本文只對垂直方向進行了驗證分析，但是庫區全方位形變監測在進行垂直形變監測的同時進行平面形變監測也十分重要[18]，今后結合GNSS空間3D信息及水汽反演成果，對深入研究InSAR水汽干擾消除及3D形變場的建立很有意義.