基于PS-InSAR的曹妃甸填海区域形变监测

李 星，陈星彤

(华北理工大学，河北 唐山 063000)

从2008年起，曹妃甸填海区为了响应‘一港双城’的政策，开始了大规模的加速发展。尤其是华北理工大学等一系列高校的搬入组成了曹妃甸大学城，京唐首钢等各种钢厂、铁厂的搬入所引起的大规模的人流扰动、基础设施的建设，势必会造成地表的扰动。为了更好地建设曹妃甸填海区以及保护人民财产安全，对曹妃甸填海区进行地表形变监测是极其有必要的。

以往的水准测量等监测地表形变受环境天气影响，而合成孔径雷达(SAR)获取数据有全天时、全天候的特点。因此学者们提出了用合成孔径雷达干涉测量(InSAR)技术、合成孔径雷达差分干涉测量技术(D-InSAR)监测地表形变。但该技术受时空失相干和大气延迟问题的影响较为严重，学者们又提出了永久散射体差分干涉技术(PS-InSAR)。Kiseleva对Caucasus的黑海海岸进行地表沉降监测，得到了很好的效果[1]。Gokhan对Istanbul地区采用PS-InSAR技术进行地面运动监测，并对沉降原因进行分析[2]。裴媛媛等、刘鹏等、Peng H等采用永久散射体时序分析(PS-InSAR)技术分别对上海市临港新城主城区、黄河三角洲东营地区、杭州湾南部沿海地区这三个地区进行地表形变监测，反演了各个地表形变过程，分析其沉降原因[3-5]。

本文详细介绍了 PS-InSAR技术用于形变监测的技术原理，在此基础上对20幅覆盖唐山市曹妃甸区的Sentinel_1A影像进行了沉降监测实验，通过取得的沉降规律，加以分析其沉降原因。

1 PS-InSAR的原理与关键技术

1.1 PS-InSAR的原理及技术流程

PS-InSAR技术是从地区一致时间不同的多幅SAR影像中提取出散射特性稳定的永久散射体(persistent scatterers，PS)点，如人工建筑体、裸露的岩石等硬物体，用提取出的PS点构建网络，并提取其相位信息，从而建立关于相位的函数模型，采用合理的算法，精确地反演地表形变。

在K幅SAR影像中，首先对主影像进行选取，将其与剩余影像配准及干涉处理，生成K-1幅干涉条纹图；去除用外部参考DEM数据模拟的高程相位，生成去平后的干涉图；在差分干涉处理后得到的平均强度图像上提取PS点，对PS点进行线性形变相位和DEM误差的估计，从差分干涉相位中将其去除；通过加入在时间域和空间域的高通、低通滤波去除获得的残余相位信息，包括非线性形变相位、大气相位、噪声相位；最后解缠形变相位，对解缠结果时序分析，得到可靠的地表形变信息。

1.2 关键技术

(1)公共主影像的选取

影响干涉相位的精度以及准确性的因素有很多，时间基线、空间基线和Doppler质心频率基线是最直接的，而这三者会随着主影像的不同而变化。因此，主影像的选取是PS-InSAR技术中的一个关键问题。多个影响因素综合考虑，让其成为最佳组合，即综合相关系数γm最大：

(1)

(2)PS点的选取

选取永久散射体(PS)点也是PS-InSAR技术中的关键步骤。该技术通过分析处理这些高相干点目标上的干涉相位，从而得到时间序列上的某地区形变信息。在处理过程中，建模的条件是影像上要有足够数量和密度的PS点。选取PS点的方法主要有相位离差阈值法、相干系数阈值法和振幅离差阈值法的单种阈值方法以及两者结合或三者结合的复合阈值方法。PS点选取的方法中最为经典的是振幅离差阈值法。该方法的基本原理是，信噪比高的点像元的振幅离差指数可以和其相位的噪声水平等价衡量。对于标准差为δ的高斯噪声的SAR影像而言，其振幅信息和 Rice 分布相似，若点具有高斯信噪比的特征，则该点的分布与高斯分布相近，即振幅离差指数DA：

(2)

式中，δA为平均强度图像的振幅标准差，mA为平均强度图像的振幅平均值[7]。

(3)干涉相位模型及形变相位反演

PS-InSAR技术中最主要、最关键的一步就是相位模型的建立以及形变相位的反演。每个提取的PS点的解缠干涉相位φunw：

φunw=φtopo+φdef+φatm+φnoise

(3)

式中，φtopo为DEM高程误差引起的地形相位；φdef为沿LOS方向的地表形变相位；φatm为大气影响相位；φnoise为噪声相位。

在时间域与空间域上DEM高程误差引起的地形相位φtopo、大气影响相位φatm、噪声相位φnoise有不同的表现形式，可通过采用高通滤波和低通滤波的方式来去除，最后，形变相位φdef：

(4)

式中，φlin是线性的形变相位；φnon-lin是非线性的形变相位，V是线性的形变速率；T是时间间隔；Dnon-lin是非线性形变。

由此LOS方向线性形变速率得以获取[8]。

2 曹妃甸填海区形变监测与分析

2.1 研究区概况及数据选取

唐山市曹妃甸区位于唐山南部、环渤海中心地带，地理坐标为北纬39°07′～39°27′，东经118°12′～118°43′。曹妃甸区的地质构造为填海造陆，土质松软，近几年曹妃甸区进行了大规模建设，引进工业、商业和学校。而地表形变会对各种基础设施的安全造成威胁，这是曹妃甸填海区地表变形监测的原因。

本文选取20景2018年8月～2019年4月，覆盖曹妃甸地区的Sentinel_1A影像，其分辨率为5 m×20 m。对曹妃甸填海区进行裁剪，裁剪区域如图1黑色线框所示，外部的DEM数据采用的是分辨率为30 m的STRM数据。

图1 研究区域

2.2 形变信息提取与分析

通过公共主影像的选取方法，最终确定2018年11月30日的图像为公共主影像，并与其他辅影像进行配准、干涉，而干涉数据的时间基线与空间基线的连接图如图2所示。从图中可以看出，空间基线的范围是在-76.246 1 ～ 96.068 9 m，时间基线范围是在-120 ～ 131 d。最后干涉处理完成后得到了19对干涉条纹图，如图3所示。图中大部分干涉条纹图的条纹不是很明显，对后面PS点的提取有极大的影响。去除干涉图中利用外部DEM数据模拟的高程相位后，得到的去平后的差分干涉图效果更好，如图4所示。

图2 数据的连接图

在生成差分干涉图之后，开始选取PS点，由于其地形条件，将振幅离差指数阈值DA设置为0.8，PS点的密度为25个/km2，总共选取了1 107 106个PS点，结果如图5所示。PS点大都落在人工建筑，如房屋、道路等硬面地物上，而不在草地、养殖池等散射特性较弱的地方。

图3 干涉条纹图

图4 去平后的差分干涉图

图5 永久散射体(PS)分布

PS点选取后，对其进行线性形变相位和DEM误差的估计，从差分干涉相位中去除得到残余相位。由于大气引起的相位延迟和噪声属于高频信号，而非线性形变相位属于低频信号，则可以分别采用时间域的高通滤波和空间域的低通滤波去除。2019年2月22日的影像能更直观的表达滤波前后的效果，如图6所示，可以看出滤波后的影像相较于滤波前降低了相干噪声，有效地去除了一些由大气引起的噪声，效果明显。

图6 差分干涉图(a为滤波前，b为滤波后)

由此，对滤波后的影像进行线性形变相位解缠，得到真实的相位信息以及各PS点在LOS方向的形变速率信息，如图7所示。

图7 曹妃甸填海区的年平均沉降速率图

由图7可以看出，整体的平均沉降速率约为-26.56～22.45 mm/a，可见从2018年8月起该地区仍然存在着较为严重的沉降。红圈范围内的A区和B区是沉降严重的区域，C区为抬升区。A区位于曹妃甸区填海区，紧靠海堤，该区域的沉降速率达到了-23 mm/a以上，沉降尤为显著。B区位于曹妃甸区的柳赞镇附近，沉降速率最大为-26.532 5 mm/a，是整个研究区域沉降最为严重的地区。C区为抬升区最明显的地区，最大抬升速率为22.440 8 mm/a。

为了更好的探讨其沉降原因，首先对大部分沉降文献进行了研究，发现除了地震等大型自然灾害的直接因素外，地表沉降的主要原因是地下水的开采和地上人工建筑物建造这两个间接因素。其次，查找该区域的水文地质资料，该区域的含水层形成原因是河流冲积和海积，属第四系松散岩类，其空隙水一般分为Ⅰ、Ⅱ、Ⅲ组，第Ⅰ、第Ⅱ含水组属潜水及微承压水，一般不作为开采层。第Ⅲ含水组为深层承压水，底界深度一般为160～500 m，水位埋深约为l5 m。20世纪80年代由于超采深层地下水形成以南堡开发区和唐海县城为中心的两个沉降漏斗，受这两个沉降漏斗的影响，该区地下水位呈下降趋势。这也是该区地面整体下降的主要原因。就A区而言，其位于曹妃甸工业区，我们对该区域进行了实地考察，发现其附近有原油码头和钢厂(首钢京唐等)。2018年首钢京唐共生产815万t，每天至少需要近10万t淡水(首钢京唐相关报道)。由此可见，该区域地表沉降的原因主要是地下水的开采。B区位于柳赞镇—一个以养殖海产品为主要经济产业的地区，其附近有许多养殖池，如虾池、蟹池。根据资料显示，2019年3月柳赞镇南部5万多亩露天海水养殖区域已全数建成并进入海水养殖准备阶段。由此可见，B区沉降严重的主要原因为养殖池的建造与开挖。C区位于唐山湾生态城，曹妃甸东站附近。生态城复合功能用地规划用地面积387.90 ha，占城市建设用地总面积的20.28%左右[9]，曹妃甸东站于2017年开始建设，2019年运营。唐山湾生态城即将被打造为临海宜居城市，建设面积在不断地增加，可见养殖池的平填与地上建设用地的建造为C区地表提升的主要原因。

为了更加直观的了解PS-InSAR技术所处理的沉降结果，分别在沉降尤为严重的A区、B区选定了两个PS点——Dot1(38°56′10.88′′N、118°31′44.12′′E)、Dot2(39°13′39.32′′N、118°37′25.82′′E)，其位置在图7中标记显示。这两个点的沉降结果如图8所示，在时间区域内，两者的沉降过程以非线性的方式沉降，且其沉降坡度逐渐增大，由此可见，A、B两区域的沉降极为严重，也间接印证了图7沉降速率的准确性。Dot1的最大累积沉降量为21.856 mm，Dot2的最大累积沉降量为15.625 mm。从这张折线图可以看出，在沉降过程中，虽然总体是呈现沉降趋势，但是某一时间区域内会有些许抬升，这主要是与该区域的地质条件、地下水开采及地上建设有极大关系。

图8 A、B区域中Dot1、Dot2的累积沉降量

3 结 语

本文采用PS-InSAR技术对曹妃甸填海区进行地表形变监测，得到了很好的监测效果。从结果看出该地区的确有着极大的地表变形，最大年平均沉降率达到了-26.558 76 mm/a，沉降区域最为严重的地方在钢厂附近，年平均沉降率达到了-23 mm/a。该区域的沉降与本身填海造陆且松软的地质条件有一定的关联，同时，工业区建设运营、长期的地下水开采、道路的开通运营会对这个本来就是填海造陆的地区产生了一定的负担。本文也证实了采用PS-InSAR技术进行地表沉降监测的有效性，为曹妃甸的开发提供了参考。