基于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形变场的建立很有意义.