1. 引言
近年来,随着新型高分辨率SAR卫星的相继发射以及时序INSAR技术的不断完善和发展,INSAR技术以其难以抵挡的迅猛优势被广泛地应用于地震、地面沉降、山体滑坡、泥石流等地质灾害调查、监测研究 [1] 。随着INSAR技术的广泛应用,时序InSAR技术逐渐凸显其优势,Ferretti等人于2000年提出PS (Permanent Scatterers)以来,以其毫米级的地表形变监测能力 [2] ,成为InSAR技术应用和研究的热点。目前,常用的时序InSAR方法有PS方法、SBAS方法(Small Baseline Subset)、StaMPS方法(Stanford Method for Persistent Scatterers)等。
任何水库大坝在建成投入使用以后,都不可避免地会发生形变 [3] ,形变监测是水库大坝安全评价的基础和重要组成 [4] ,现阶段水库形变监测多是在大坝的重点部位,或坝区的重点区域进行布设监测点,并且监测周期长,监测成本高,只能以点代面的方式来整体评价水库大坝安全运行情况,难免会遗漏一些重大的安全隐患。另外,大坝坝区及周边的地表形变也会对大坝安全运行存在潜在影响。因此,采用时序InSAR技术对水库坝区进行形变监测,能够获取大面积形变监测结果,并有效缩短监测周期,保证了监测精度。2011年,王腾等利用QPS技术(Quasi-PS, QPS)获得了坝区的实施形变结果,并分析了三峡大坝运行时的稳定性问题及其周边区域的形变规律 [3] ;2012年,Voege等利用SBAS方法对挪威斯瓦特湖坝形变监测,探测到坝体的局部形变情况;2013年,裴媛媛等利用改进的PS技术得到了长江口南岸和杭州湾北侧堤坝的平均沉降速率,并与水准数据进行比较,监测形变度达到4.5 mm;2017年,有学者利用PS和SABS两种方法对伊拉克摩苏尔大坝及周边区域的形变监测,并对大坝、滑坡体、排水口等重点区域进行形变分析;2019年,肖儒雅等利用改进的PS方法对广南水库进行形变监测,验证了时序InSAR分析方法对水库库区、大坝、防潮堤等水工建筑物进行形变监测的有效性 [4] 。
本文综合考虑了新疆WLL水库大坝及周边地貌、地物特点,以及新疆冬季积雪深厚并持续时间较长的特殊性,选用SBAS方法对该水库坝区进行形变监测应用研究。SBAS方法可以将较少的影像数据进行合适组合,得到一系列短空间基线差分干涉图来较好地克服空间去相关现象 [5] ,并对高干涉图进行多视处理以提高像元信噪比,使部分原本低相干的像元被有效利用,较适用于自然地表的形变反演 [6] [7] ,特别是在SAR数据较少,时间跨度和时间间隔较大的情况比较适用 [8] 。
2 本文研究方法
SBAS通过简单和高效地合成所有可用的小基线干涉图对,然后再基于形变速率的最小范数准则,通过引用奇异值分解(SVD)方法获取相干目标的形变速率及其时间序列 [7] 。其基本原理如下。
假设获得统一区域按照时间顺序t0,t1,…,tN排列的N + 1幅SAR影像,选取其中一幅影像作为主影像,并将其他SAR影像配准到这幅影像上。N + 1幅SAR影像生成M幅多视差分干涉图。对于从影像tA和tB (tB > tA)时刻获取的第j幅差分干涉图,方位向坐标为x和距离向坐标为r的像素干涉相位可以写为:
(1)
式中:
;λ为电磁波波长;
和
为tB和tA时刻相对于
的雷达视线方向的累积形变量;
为差分干涉图中残余地形相位;
为大气延迟相位;
为去相干噪声。
在不考虑大气延迟相位、残余地形相位和噪声相位的前提下,可以将(1)式简化为:
(2)
将式(2)中相位表示为两个获取时间之间的平均相位速度和时间的乘积:
(3)
第j幅干涉图的相位值可以写为:
(4)
即各时段速度在主、从影响时间间隔上的积分,写成矩阵形式为:
(5)
B为一个M × N的矩阵,v为一个N × 1的矩阵,
为一个M × 1的矩阵。由于SBSA的差分干涉图采用多主影像策略,因此,矩B容易产生秩亏,采用奇异值分解(SVD)可以得到矩阵B的广义逆矩阵,进而得到速度矢量v的最小范数解,最后通过各个时间段内速度的积分可以得到各个时间段的形变量。
3. 研究区概况和获取数据情况
新疆WLL水库建于博格达山北麓和准噶尔盆地南缘的冲积平原上,是典型的平原水库,平均海拔495 m,平均坡度为2%左右,周围地物分布主要以戈壁、农田为主,距水库5 km范围左右分布有工业区和居民区。水库地属大陆性干旱气候,每年10月中旬开始降雪,次年4月份消融,年平均降雨量251.59 mm,年蒸发量1543.84 mm,平均潮湿系数0.163,属于湿度过低带。WLL水库面积约23.8 km2,总库容2.8亿m3,经四面筑坝整个形状呈封闭四边形,属大(2)型水库;大坝建在低液限粉土地基上的碾压式均质土坝,最大坝高28 m,分东、西、南、中四个坝段,四个坝段总长约18 km。
![](//html.hanspub.org/file/9-2840447x21_hanspub.png?20240515083223341)
Figure 1. Schematic diagram of the location of the study 168 area
图1. 研究区域位置示意图
本研究选取了2014年10月~2018年9月23期Sentinel-1A影像数据,去除了每年11月到来年3月份的数据,其他时段保持1个月一期,以2016年8月10日影像为主影像为例,时空基线参数见表1。轨道数据采用欧洲航天局(ESA)发布的精密轨道数据,DSM选用日本宇宙航空研究开发机构免费提供的AW3D30数据,空间分辨率为30 m,高程精度达5 m。图1显示了本研究区域范围,底图为DSM数据,黑色框为数据计算范围,红色框为本研究区域范围。
![](Images/Table_Tmp.jpg)
Table 1. Interferometer image pairs spatiotemporal baseline table
表1. 干涉像对时空基线表
4. 数据处理
本文选用架构于专业遥感处理软件ENVI5.3上的INSAR数据处理模块SARScape进行SBAS方法处理,具体流程如图2所示。
对23幅单视复数影像(Single Look Complex, SLC)建立工作流进行SBAS方法计算,由于研究区域属于干旱性气候,设置时间基线阈值为781天;空间基线阈值设置为临界基线的45%,共有22幅影像作为主影像参与干涉像对的生成,共计生成253对,最大空间基线为207 m,最小空间基线为4 m;选用2016年8月10日影像作为超级影像进行配准,借助外部精密轨道数据和DSM数据配准精度可以达到1/8像素;SBAS方法在去噪上采用前置滤波的多视处理(方位向和距离向视数比1:5)和后置滤波的Goldstein滤波方法;相位解缠相干系数阈值设为0.35;在解缠完成以后,需要人工检查每个干涉像对的相干性图和解缠图,对于区域相干性低,解缠结果较差的干涉像对进行删除处理,以免对形变结果产生较大的误差;轨道精炼和重去平这一步骤需要人为选择地面控制点来进行多项式拟合来消除可能的斜波相位,对卫星轨道和相位偏移进行纠正,地面控制点须选择在地面平坦、相干性高、解缠精度高的稳定区域。在对研究区域整体形变情况未知的前提下,地面控制点选择的精度直接影响最终形变结果精度,本文地面控制点的选取采用自动选取的方法,避免人工选取引入较大的误差,通过相干性、振幅离差指数、形变速率3个阈值提取稳定的PS点作为地面控制点 [9] ,具体参数见表2。第一次形变估计会对轨道精炼和重去平的干涉图进行第二次解缠,大大提高了解缠精度,最后通过大气低通和高通滤波估算大气相位,并对非线性形变进行分离,最终形成研究区形变时间序列。
![](Images/Table_Tmp.jpg)
Table 2. Ground control point parameters table
表2. 地面控制点参数表
5. 数据分析
通过以上处理,获得了研究区的沿雷达视线方向的形变速率和相对于2014年10月20日第一期的形变时间序列。为了验证研究区域SBAS方法形变监测的精度,将雷达视线方向的年变速率时间序列换算到垂直方向,利用同时段大坝迎水坡坝顶处和背水坡坡脚处34个二等水准监测点成果(2014~2018)进行验证,部分结果如表3所示,中误差为±9.3 mm。水准监测点位置分布见图4。
图3截图了坝区周围5 km左右范围的形变速率图,总体上来看,堤坝上存在不均匀沉降,中坝(北面)、南坝中间位置和靠近形变区2的东坝中段沉降量较大。坝区周围存在5个形变区漏斗,距离大坝最近的形变区2是水库水产养殖中心,最大沉降速率(沿雷达视线方向)达到−20.0 mm/a,其余4个形变区漏斗分布在距离大坝5 km左右范围的工业区和居民区,形变区3为某一工业区,其形变影响面积最大,并逐年向水库扩展,该漏斗中心沉降速率到达−76.0 mm/a。
![](//html.hanspub.org/file/9-2840447x23_hanspub.png?20240515083223341)
Figure 3. Diagram of deformation rate in the dam area
图3. 坝区形变速率图
对大坝的不均匀沉降进一步分析,如图5所示,受水库水体外力影响,大坝迎水坡的坝顶沉降年变速率均大于背水坡坡脚的年变速率,平均年变速率分别为−23.0 mm和−12.5 mm;中坝段年变速率明显高于西坝段和东坝段,除了受到大坝自重、水体外力等影响外,与中坝北部区域开荒种植抽取地下水,引起地下水位发生变化有关。为了进一步掌握中坝段和西坝段沉降时间变化情况,在中坝段和西坝段沉降较大区域选择两个水准监测点Z215、Z216与同时段、同位置的SBAS结果进行对比分析,如图6和图7所示,SBAS方法计算两个监测点沉降时间序列变化趋势与水准保持一致,量值上有差别。水库基本上在每年4月份开始上游来水至10月份结束,受水库每年蓄水时段的影响,监测点呈现4月至10月沉降量逐渐增大,11月至来年3月沉降量逐渐减小的周期趋势。西坝受形变区2(水产养殖中心)形变区的影响,2016年至2017年沉降量突变,到2018年沉降趋势开始变缓。自水库蓄水以来,水库运行水位从未触及到南坝坝体,故未在南坝设置形变监测点,但从SBAS监测结果来看,如图3所示,受周边草坪灌溉抽取地下水的影响,南坝段中部位置沉降量较大,沿雷达视线方向沉降年变速率最大达到−18.0 mm/a。
![](//html.hanspub.org/file/9-2840447x24_hanspub.png?20240515083223341)
Figure 4. Distribution map of leveling monitoring points in the dam area
图4. 坝区水准监测点分布图
![](Images/Table_Tmp.jpg)
Table 3. Comparison between dam level and SBAS monitoring results
表3. 大坝水准与SBAS监测结果对比表
![](//html.hanspub.org/file/9-2840447x25_hanspub.png?20240515083223341)
Figure 5. Comparison of annual settlement rate of longitudinal sections of upstream and downstream slopes of the dam
图5. 大坝迎水坡、背水坡纵断面沉降年变速率对比图
![](//html.hanspub.org/file/9-2840447x26_hanspub.png?20240515083223341)
Figure 6. Time series diagram of Z215 level and SBAS settlement
图6. 中坝Z215水准与SBAS沉降时间序列图
![](//html.hanspub.org/file/9-2840447x27_hanspub.png?20240515083223341)
Figure 7. Time series diagram of Z216 level and SBAS settlement
图7. 中坝Z216水准与SBAS沉降时间序列图
对坝区及周边影响区域进一步进行形变分析,距离西坝段不到1 km的水产养殖中心区沉降区,产生的原因主要是由于抽取地下水过度所致,由于西坝段自2016年沉降加速,及时停止了地下水的开采,自2018年形变趋势变缓。受周边工业区、居民区以及农田种植地下水过度开采的影响,距离坝区5 km左右范围的4个沉降区,形变量值和影像面积较大的形变区3是某一工业区,位于水库的西南方向,可能会对水库产生潜在安全隐患,应及时进行人为干预。
6. 结论
本文利用欧洲航天局中等分辨率Sentinel-1A免费影像数据,采用SBAS方法对新疆WLL平原水库坝区进行了形变监测研究,利用大坝二等水准监测结果进行了对比分析和精度验证,并利用SBAS监测结果分别对大坝、坝区、坝区周边区域进行了形变分析。结果表明:1) 坝体变形SBAS监测结果与常规二等水准测量结果基本吻合,后期可持续进行SBAS监测并适当减少二等水准的观测频次,有助于降低常规观测成本;2) 坝区及周边面状区域的SBAS监测对周边农田抽取地下水引起的地表形变“漏斗”的监测效果较好,大坝周边的地表形变也是影响大坝安全的因素之一,SBAS监测较常规测绘手段在监测大坝周边区域方面有明显优势。本次结论也进一步印证了时间序列InSAR方法在水库坝区形变监测上的可靠性和实用性,为整体评价水库大坝安全运行管理提供了可靠的保证。
参考文献
NOTES
*通讯作者。