Sentinel-1A卫星的黄河三角洲近期地表形变分析
Sentinel-1A卫星的黄河三角洲近期地表形变分析
程霞1,2，张永红2，邓敏1，吴宏安2，康永辉2
1.中南大学 地球科学与信息物理学院，长沙410083；
2.中国测绘科学研究院，北京100036
为了实现对黄河三角洲近海沿岸大范围形变监测，该文采用欧空局新一代SAR卫星哨兵一号(Sentinel-1)TOPS模式数据进行监测。该文利用多主影像相干目标小基线干涉技术(MCTSB-InSAR)对覆盖黄河三角洲地区的整景Sentinel-1A数据进行了时序干涉处理，提取了整个黄河三角洲近海沿岸地区2016年1月至2017年5月22景Sentinel-1A的地表沉降信息，并对主要的沉降中心分析了其沉降原因。结果表明，黄河三角洲近海沿岸地区存在10个显著的沉降中心，河口区广河村国星兴盐场内最大年均沉降速率达464.5(mm·a－1)。沉降原因主要是由于抽取地下卤水进行工业制溴、制盐，油田抽采(采油和抽取地下水用于回注)等，并且地下卤水资源的开采对黄河三角洲近海岸区域地面沉降的贡献最大。
本文目录0 引言
1 Sentinel-1A高精度配准及时序分析
2 研究区Sentinel-1A SAR影像及形变结果
3 黄河三角洲近海岸地区地表沉降分析
4 结束语
黄河携带泥沙流入渤海处淤积形成面积约为5 500 km2的现代黄河三角洲，粉质黏土状的沉积物在海岸线附近达15～20 m厚。高度可压的土壤使得三角洲沿岸部分极易受到下沉影响。此外，黄河三角洲地区石油、天然气、地下卤水等资源丰富，此处是我国开发地下卤水成盐最早、规模最大的地方，其盐产量占全国盐产量的一半，并且中国第二大油田胜利油田同样位于该地。除了进行地下卤水资源的开发外，地下卤水还可用于发展养殖业、化工工业，已有研究证明大型三角洲沿岸养殖区每4年可导致1 m的地面沉降[1]。因此，近几十年来黄河三角洲地下卤水的充分开发利用，使得大片盐田、养殖池迅速扩大[2-3]，同时油田的开发利用也在不断进行，在创造了巨大的经济效益的同时也引起了严重的地面沉降。已有研究人员对黄河三角地区特殊的地质构造进行了大量以水准测量和全球定位系统(global positioning system，GPS)为手段的地面沉降监测[4-6] ，但其只是针对单个点，无法提供大区域连续形变监测方案且耗费巨大，同时也有学者利用合成孔径雷达干涉测量(interferometric synthetic aperture radar，InSAR)技术对黄河三角洲地区进行形变监测[7-9]，但这些研究并没有兼顾影像幅宽大且影像获取日期较新的特点。因此无法反映目前黄河三角洲的沉降状况，特别是变化信息极大的沿海地区，其沉降速率几乎是未知的。
2014年欧空局发射的Sentinel-1A卫星随后公开发布免费数据源，可实现大区域地面监测的同时抑制了scalloping效应并增强了成像辐射性能，但其特殊的渐进扫描(terrain observation by progressive scans，TOPS)成像模式使得方位向和距离向同时摆动而引入额外的脉冲重复频率，导致合成孔径雷达(synthetic aperture radar，SAR)影像在方位向上的多普勒中心变化。目前国内外学者已成功将TOPS模式的Sentinel-1A数据进行时序分析应用于监测地表沉降[10]、地震[11]、火山活动[12]等方面，但Sentinel-1A数据公开较晚，干涉处理具有一定的挑战性，目前基于整景宽幅影像的时序分析处理国内外仍较少。本文利用的基于改进的多主影像相干目标小基线干涉(multiple images coherent targets small baselines InSAR，MCTSB-InSAR)技术[13]结合了永久散射体合成孔径雷达干涉测量(persistent scatterer interferometric synthetic aperture radar，PS-InSAR)技术和小基线集(small baseline subset，SBAS)技术的优点，有效提取出城区及郊区的高相干点目标。MCTSB-InSAR技术进行小基线干涉组合后对邻域高相干点进行差分运算，并将相位解缠与形变信息、高程残差信息一起解算[14]，避免了PS-InSAR所需影像数量大(一般大于25景影像)和SBAS对每一幅干涉图都需要解缠的缺点。该方法已成功应用于多个地区的地表沉降监测[15-16]。
1.1 Sentinel-1A高精度配准
影像配准是SAR数据干涉处理的首要任务，配准精度越高，干涉图的相位误差将越小。与条带模式的InSAR数据不同，TOPS模式的InSAR数据由于在方位向上有很大的多普勒中心频率偏移，对配准提出了更高的要求。由多普勒中心频率偏差引入的干涉相位偏差可表示为：
式中：ΔfDC为Burst区域多普勒中心频率变化值；Δβ为相应的方位向上视角差；Δt为主辅影像对应的Burst的配准误差。对于Sentinel-1A，ΔfDC最大为5.2 kHz，Δβ为±0.6°，为使相位坡度小于1/100 rad，方位向的配准精度需优于0.000 9像元[17]，约合1.9μs或1.3 cm。
针对Sentinel-1A数据方位向配准精度需优于0.000 9像元的特点，使用增强谱分集(enhanced spectral diversity，ESD)[18]方法进行方位向精配准。ESD方法利用方位向上相邻Burst之间重叠区域的干涉相位偏差进行方位向精配准，重叠区域任何样本点p的ESD干涉相位差可表示为：
式中：mi、si分别表示主辅影像第i个Burst的重叠区域；mi+1、si+1分别表示主辅影像第i+1个Burst的重叠区域；arg·表示为复数形式的相位。对式(1)进行重写，可求得任何样本点p基于像素级别的方位向配准误差Δy。公式如下：
式中：ΔfovlDC，p表示在重叠区域的多普勒中心频率变化值；faz为影像方位向采样频率。利用二维周期图法[17],并将影像相干性作为权重值来求解逼近真实的方位向配准误差Δy，以下为估计值Δ的简化模型：
由于ESD方法计算得到的解缠差分干涉相位区间为[-π，π),相当于0.05个像素精度，因此在使用ESD方法进行相位解缠前，需要尽可能使用精确的初始配准方法，从而避免相位模糊的问题。本文首先联合几何配准[18-19]和影像相关性配准[20]作为初始配准，在此基础上使用ESD方法计算两幅影像Burst重叠区域之间的相位偏差以估算配准偏差，进一步提高配准精度。需要注意的是，ESD方法的方位向配准精度由两幅SAR影像间的相干性决定，因此时间去相干在一定程度上影响着ESD方法对方位向的配准精度[21]。本文使用基于网络的ESD方法[22]来估计俩相邻影像之间的方位向偏移，解决了两影像之间的时间失相干问题，获得远优于0.000 9像元的方位向配准精度，从而避免了干涉图中出现的相位跳跃现象。在配准的基础上对辅影像进行重采样，最后利用MCTSB_InSAR进行时序处理，具体数据处理流程图如图1所示。
图1 MCTSB-InSAR时序分析技术数据处理流程图
1.2 MTCSB-InSAR时序分析技术
本文仅简述MCTSB-InSAR方法时序处理方面几个关键点，更多细节请参阅参考文献[13-14]。首先根据干涉像对必须连通和干涉相干性最优的原则[13]，对精配准后的SAR单视复数据(single look complex，SLC)设置适当的时空基线阈值来确定小基线对组合；然后用平均相干系数、平均幅度、幅度离差等三参数阈值原则[16]提取稳定的高相干点目标，并采用改进的Delaunay三角网连接高相干点目标[23]；最后利用形变模型将地表形变分解为线性和非线性形变。对于第k个干涉图上两相邻高相干点i、j的干涉相位差可表示为如下公式：
式中：Δψki，j为第k个干涉图中的高相干点的差分相位之差；Δαki，j为相邻高相干点的模型相位，包括两点之间的高程误差之差ΔΖi，j和线性形变速率之差ΔVi，j贡献的相位；Δβki，j为非线性形变相位差；Δφkij，atmo表示大气相位之差；Τki、Βki分别为差分干涉图k的时间基线和空间基线；Δφkij，noise为噪声相位差；λ为雷达波长；θ表示雷达波束的入射角。
式(5)中相邻高相干点之间的相位是缠绕的，因此无法直接利用最小二乘求解，此时，当研究区域中非线性形变的相位变化不是很大的情况下，ΔWi<π这个条件极有可能成立，则可利用周期图法[13]构建时态相关因子求解高程残差之差ΔZ和形变速率之差ΔV，公式如下：
式中：j为单位复数；M为干涉图数目。根据实际经验一般保留满足γij≥0.7的高相干点，选择一个已知数字高程模型(digital elevation model，DEM)误差且稳定不变的点作为基点，对Delaunary三角网中的所有边的速度和高程误差增量进行积分，获得各个点目标的形变速率和高程误差。利用时空滤波分离出大气相位和噪声相位，并将得到的非线性形变相位与线性形变进行叠加,便可得到完整的时序形变信息。最后，将雷达视线向的形变量转换为垂直方向上的地面形变量。
本文选取2016年1月—2017年5月22景覆盖研究区域的Sentinel-1A TOPS模式升轨影像，影像覆盖范围具体位置见图2所示，以Google Earth影像为背景，多边形框为Sentinel-1A影像覆盖范围。数据格式为单视复数据并对其做8×2的多视处理，降低了speckle噪声，提高了相干性，最终距离向和方位向像元大小分别为18.64和27.92 m，选用30 m格网间隔的SRTM DEM去除地形相位。影像的时空基线参数见表1，时空基线阈值分别为90 d，180 m，最终干涉像对时空连接如图3所示。本文选定2016年1月9日的影像作为统一配准基准，并在后续的时序分析处理中，作为时间参考点。通过高精度的影像配准及MCTSB-InSAR时序分析技术后，最终获得整景Sentinel-1A影像的平均形变速率图，如图4所示。
图2 研究区地理位置
表1 Sentinel-1A影像时空基线
图3 干涉像对连接图
图4 2016—2017年研究区平均沉降速率图
通过以上InSAR时序处理，获得2016年1月—2017年5月黄河三角洲平均沉降速率，结果表明在监测时间段内，黄河三角洲仍然存在较为严重的沉降区域，地面沉降分布不均匀，其中在新近形成的陆地、盐场、油田附近的沉降量较大。为了更加详细了解黄河三角洲近年来的地表沉降情况及其影响因素，在研究区选择10个沉降中心点，即图4中的A-J标识作为典型代表进行分析，并于2019年1月14日开始对这10个沉降中心进行了实地调研。经总结，将造成黄河三角洲较大沉降的主要因素分为两类：地下卤水资源开采、油田抽采。下面从地下卤水资源开采和油田抽采两个层面进行定量分析。
长期不合理开采地下卤水将会引起地下卤水水位下降，形成卤水漏斗，引发地面沉降。地下卤水常被化工厂用来制取溴化合物和盐，通常化工厂提取完溴素后，尾水通过排水渠输入盐田晒盐。A、B两沉降中心都处于典型浅层地下卤水分布区(无棣一沾化卤水分布区)，E、F、G沉降中心属于东营浅层卤水分布区。其中，垦利区红光村、六十户村(E、F区)均为大规模抽取地下卤水制溴素随后引入排水沟排放至盐田日晒成盐，其平均沉降速率分别为82.5、118.7(mm·a-1)，最大沉降速率分别为120.6、202.9(mm·a-1)，其排水沟、盐田池、盐堆如图5所示。河口区广河村国星盐场内(C区)存在一个明显的沉降中心，平均沉降速率为160.6(mm·a-1)，最大年均沉降量达464.5(mm·a-1)，沉降最大处位于该厂集中呈队列排布的卤水抽水井处，如图6所示。国星盐场四周均为油田开采区，该盐场极有可能超采地下卤水制溴素晒盐，使得该区域地面沉降量明显比周围油田区大，如图7所示。羊口盐场(G区)于20世纪60年代起就开始对地下卤水进行开发，之后羊口镇开始利用地下卤水发展养殖业、化工工业等。如图8所示，该区最大年均沉降速率为264.3(mm·a-1)，受此沉降中心的影响附近的羊口镇城镇年平均沉降速率37.4(mm·a-1)，化工工厂所在区域最大沉降速率更是达到163.8(mm·a-1)，有些地方出现明显的墙裂缝，如图8(c)所示。由于影像范围的限制，作为典型的因深层地下水开采过度而出现较大沉降漏斗的广饶县[24]，只显示出了该沉降漏斗范围内的广饶县花官乡(H区)，其最大年均沉降速率达80(mm·a-1)。如图9所示，下方的曲线表示A、B、C影像区中实线覆盖区域的每个像素的平均沉降速率，下方的虚直线表示A、B、C影像区中实线覆盖区域的总体平均沉降速率。通过对这些典型浅层地下卤水开采区的总体分析可知，地下卤水资源的开采带来的地面沉降量级大，总体平均沉降速率为50～160(mm·a-1)。
油田区发生地面沉降主要原因在于长期高强度的油气开采导致储油层压力下降，油气区和黏土层压缩，从而发生地面沉降。另外，石油开采过程中一般会抽取附近浅层地下水进行回注，有可能导致注水区域地面抬升而抽水区域地面沉降。80年代末桩西和孤东油田(D区)开始石油开采活动，该区最大沉降速率为449.7(mm·a-1)，年平均沉降速率为286.3(mm·a-1)。如图10所示，图10(c)、图10(d)为沉降最大区域，即图中附带箭头处的小矩形框所指区域的局部放大图，其西北端为桩西油田，东南端为孤东油田，西侧为仙河镇及孤岛油田、孤北水库，该区的沉降原因不仅仅只是石油的抽采，在近海岸区将石油抽采与盐田开发相结合已成为常见形式，具体如图10(f)中位于沉降中心处的实景，抽油机抽采石油的同时，附近会有注水装置(图中设备)将抽取的地下水回注以补偿抽油亏空，抽油机四周同样也存在抽取地下卤水制盐的盐田区域。这种在不同地层同时抽取地下卤水、浅层地下水、油气的现象更是加剧了该区域的地面沉降。石油抽采过程中也会抽取地下水进行回注，回注区域可能会发生地面抬升现象，图10(e)中为孤东油田平均沉降速率图，整个D1区域除西北方向存在部分石油与卤水资源同时抽采的综合区域外，其余区域均只进行石油抽采活动。这里仅讨论石油抽采区域，其最大抬升速率为64(mm·a-1)，最大沉降速率为60.3(mm·a-1)，平均沉降速率为23(mm·a-1)。由图10(a)可以看出，D1区的桩西油田，D3区的孤东油田的沉降量级明显比D2区的盐田、油田综合区域小，为了进一步定量分析抽取地下卤水资源和抽采石油对区域沉降影响程度的不同，对图10(a)中跨过D1、D2、D3区域的黑色实线覆盖区域绘制出线状区域的平均沉降速率曲线，如图11所示，综合区域D2的平均沉降速率量级最大，油田区D1、D3平均沉降速率都比较稳定，且平均沉降速率均值分别为9.2、21.4(mm·a-1)，与之对比，D2区域的沉降速率变化大且平均沉降速率均值达195.0(mm·a-1)。因此，卤水资源的开采比石油抽采活动对D区域造成沉降的影响更大。东营市西城区(I区)建城较早，沉积物基本也已固结压实，区内地质构造较稳定，该区的沉降主要是受石油开采活动影响，最大年均沉降速率达33.4(mm·a-1)，年平均沉降速率为24.1(mm·a-1)。西城区油气比较丰富，处于东章油田、史南油田和六户油田之间，油气密布的西城区在近年沉降量依然比东城区严重，这与文献[4，24]所述一致。桩西油田(D3区)、孤东油田(D1区)、东营市西城区油田(I区)对地面沉降造成的影响年均沉降量为10～25(mm·a－1)，造成D2综合区域年沉降量远大于25(mm·a-1)的主要原因是地下卤水资源的开采。
图5 地下卤水开发利用实景图
图6 盐田区域抽水井
图7 C区平均沉降速率及对应的光学影像
图8 G区平均沉降速率及对应的光学影像
图9 卤水开采区平均沉降速率特征图
图10 D区平均沉降速率及对应的光学影像
对重点分析的典型代表区域C、G、D中，选取最大沉降点c、g、d点(分别在图7、图8、图10中用箭头指出)绘制其累计形变时间序列，如图12所示，从图中可以看出c、d、g 3个最大沉降点在研究时段内沉降量基本呈线性变化，沉降速率趋于稳定。而c、d点的累计沉降量明显大于g点，除了与地下资源的开采量有关外，还可能与土层的固结压实程度有关。
图11 剖线区域平均沉降速率图
图12 累计形变时间序列曲线
语本文阐述了整幅Sentinel-1A TOPS数据进行了干涉处理及时间序列InSAR分析，重点分析了黄河三角洲近海岸地区短周期、近实时沉降监测。研究结果表明：
1)Sentinel-1A数据时间基线短，影像数量多，十分有利于InSAR时序分析，再结合MCTSB-InSAR技术中的干涉组合方法，获得高质量高数量的差分干涉图。
2)在2016年1月—2017年5月监测时间段内，黄河三角洲存在较为严重的沉降，且沉降主要分布在沿海岸的卤水资源开采区、石油抽采区及部分城镇区域，不均匀的沉降可能对重大基础设施造成安全隐患。通过实地调研，本文详细分析了造成黄河三角洲10个显著的沉降漏斗的主要原因，即大量地下卤水开采和广泛的油田抽采活动，并分析出地下卤水开采和石油抽采活动对地面沉降的贡献量的不同，总结出研究时段内造成黄河三角洲大规模大幅度的地面沉降的主要原因是大量地下卤水资源的开发利用。部分区域与已有文献中监测结果比较，沉降分布规律及相对大小基本吻合，与以往文献中的不同之处在于本次研究几乎完整地覆盖了黄河三角的近海岸地区，并对其进行了最新的地面沉降监测分析。
3)文中仅对整景的Sentinel-1A数据的近海岸部分进行了分析，今后可以对处理的其他区域进行分析研究，同时收集区域的地表高精度水准测量数据进一步验证监测结果的正确性和MCTSB-InSAR方法的适用性。
致谢：感谢欧洲太空局提供的Sentinel-1A TOPS影像数据。
参考文献(略)
作者简介：程霞(1992—)，女，湖南长沙人，硕士研究生，主要研究方向为合成孔径雷达遥感。
基金项目：国家重点研发计划项目(2017YFE0107100，2018YFB0505404)；中国测绘科学研究院基本科研业务费项目(7771610，7771624)
通信作者：张永红 研究员
E-mail：yhzhag@casm.ac.cn
引用格式：程霞，张永红，邓敏，等.Sentinel-1A卫星的黄河三角洲近期地表形变分析[J].测绘科学，2020，45(2)：43-51.(CHENG Xia，ZHANG Yonghong，DENG Min,et al.Analysis of recent surface deformation of the Yellow River delta based on Sentinel-1A satellite，2020，45(2)：43-51.)