基于PS-InSAR监测的石油开采地表下沉模型反演

王远坚，姜 岳，MISA R，SROKA A，姜 岩

(1.山东科技大学测绘与空间信息学院，山东 青岛 266590；2.中国矿业大学环境与测绘学院，江苏 徐州 221116；3.波兰科学院岩层力学研究所，克拉科夫 30-059)

石油是重要的基础能源和工业原材料，关系着国民经济命脉和国家安全。石油开采导致的缓慢地表移动变形现象在全球广泛分布，这不仅会破坏油井的生产设施，给石油开采造成巨大的损失，同时会对建筑物和周边生态环境产生损害性影响[1-2]。某油田是一个具有多套生油层系、多种储集类型、多种油气藏的复式油气区，地表移动变形现象在持续发展，因此必须采取相应手段对该油田进行监测。传统的监测手段需要耗费大量人力、物力、财力，而合成孔径雷达干涉测量技术(interferometric synthetic aperture radar，InSAR)是一种空间对地观测技术，具有高空间分辨率、高精度、监测范围广等优势，已经被广泛应用于地表形变监测、灾害监测与预警、地球物理参数反演等工作中[3-4]。针对石油开采导致的地表下沉具有缓慢累积的特点，采用PS-InSAR(permanent scatters InSAR)技术进行监测，该技术能够对长时间以固定重访周期稳定、连续获得的SAR影像进行时间序列处理，利用在长时间和空间基线下能够保持高相干的点目标，克服时间、空间去相干对干涉信号的影响，并利用形变相位和大气效应相位在时间序列上的不同时空特性将两者加以分离，得到地表形变演变参数，这使该技术对于监测缓慢地表形变具有较大的优势[5]。

目前关于油田开采导致的地表形变监测研究较多[6-7]，但对于地表形变的进一步研究则相对较少。在油田地表形变模型研究方面有部分学者将火山形变反演中的Mogi模型、地震形变反演中的Okada

模型引入[8-9]，这些模型能够大致反映油田区域的地表下沉趋势，但对于地表累积下沉量较大区域及下沉边界区域模型拟合结果存在较大误差。因此，根据地表实际下沉情况，确定一个符合目标区域的地表移动和变形模型能够为评估开采引起的地表移动变形的损害风险提供依据。

1 研究区概况

该油田是一个以重油开发为主的复杂断块油田，面积约200 km2，采油厂位置及数据处理范围如图1所示。

本文采用覆盖研究区域的2017年3月至2019年11月的33景C波段升轨Sentinel-1A数据进行下沉监测，数据具体情况见表1。数据的距离向分辨率为5 m，方位向分辨率为20 m，数据处理范围如图1所示，面积约为1 736 km2。

图1 石油开采区域及数据处理区域Fig.1 Oil production area and data processing area

表1 Sentinel-1A数据参数表Table 1 Sentinel-1A data parameters

续表1

2 PS-InSAR技术基本原理及地表移动变形模型

2.1 PS-InSAR技术基本原理

PS-InSAR技术能够克服时间、空间去相关对干涉信号的影响，并利用大气效应相位和形变相位在时间序列上的不同时空特性将两者加以分离，得到地表形变演变参数[10-11]。PS-InSAR技术的基本原理是假设有K+1幅SAR单视复数影像，经配准、辐射定标、PS探测和干涉处理，并借助已知DEM进行差分干涉处理后得到K幅干涉和差分干涉图、H个PS点以及各PS点在各差分干涉图中的差分干涉相位集。在考虑地表形变、高程误差、大气影响及失相关的情况下，每个PS点在每幅差分干涉图上的差分干涉相位组成，见式(1)[12-13]。

φdiff=φtopo+φdef+φorb+φatm+φnoise

(1)

通过设定振幅离差指数、相干系数、旁瓣阈值提取PS候选点，后对这些PS候选点进行噪声计算，去除仅在个别干涉图中相位稳定或受到邻近PS点影响而表现出PS点特征的象元，剩余的即为PS点。估计并剔除PS点相位的空间不相关误差，并更新PS点相位；对更新后的PS点相位进行时空三维解缠；分析各干扰误差相位在时间和空间上的特征差异，对解缠相位进行时、空滤波，分离并剔除各误差相位，最后提取形变相位。PS-InSAR的解算结果是视线向的形变值，为了准确反映研究区域的下沉变化，需将视线向的解算成果转换为垂直向的结果。该转换可通过式(2)来进行[14]。

(2)

式中：d为下沉值；Δr为视线向形变值。

2.2 地表移动变形模型

为了研究开采对地表油井、建筑物及环境的影响，需要根据PS-InSAR监测结果，模拟计算出地表下沉盆地数学模型，为计算地表移动与变形提供基础。关于石油开采地表下沉计算有许多方法[15-17]，但是，部分方法计算比较复杂，本文结合研究区域的地表下沉分布特征，将地表下沉视作椭圆形弹性薄板挠曲，构建基于椭圆薄板模型的地表下沉模型[18-20]。如图2所示，图中q为薄板上方的横向载荷；w为地表下沉量；a为椭圆下沉盆地长轴半长，即椭圆板长半轴；b为椭圆下沉盆地短轴半长，即椭圆板短半轴。

图2 弹性椭圆薄板的坐标和位移示意图Fig.2 Coordinates and displacements ofelastic elliptic thin plates

弹性椭圆薄板的挠度计算模型边界条件表达式见式(3)[21]。

(3)

基于弹性椭圆薄板的挠度计算模型边界条件及下沉区域空间形态，可建立监测区域地表下沉模型表达式见式(4)。

(4)

式中，k为常数，即为需要反演的值。

根据地表移动变形理论[22-23]，由式(4)可以求出，地表倾斜i(x,y,φ),曲率K(x,y,φ),水平移动U(x,y,φ),水平变形ε(x,y,φ)。

沿φ方向的倾斜，见式(5)。

(5)

沿φ方向的曲率，见式(6)。

(6)

沿φ方向的水平移动，见式(7)。

U(x,y,φ)=Bri(x,y,φ)

(7)

沿φ方向的水平变形，见式(8)。

ε(x,y,φ)=BrK(x,y,φ)

(8)

式中：B为水平移动系数；r为主要影响半径；φ为沿x轴正向逆时针计算到指定方向的角值。

3 监测结果及时空分析

以2018年6月17日获取的影像为主影像，其他32景影像为辅影像，运用PS-InSAR技术对获取的影像进行处理，获取了2017年3月至2019年11月的监测成果，在数据处理范围内共提取到401 728个PS点，采油厂监测范围内共提取到17 796个PS点，如图3所示。由于该区域内存在大面积的芦苇荡，造成众多区域失相干，未提取到PS点。从图3中

可以看到，该采油田范围内存在明显的下沉漏斗，下沉现象较严重，下沉区域面积达29.5 km2，最大下沉速率为-210.9 mm/a，最大累积下沉量达-585.6 mm，这会对地表的建构筑物和公路等设施产生影响，甚至造成破坏。

图3 2017年3月至2019年11月下沉速率图Fig.3 Subsidence rate map from March 2017to November 2019

为了更深入准确研究该区域的下沉状况，按季度刻画出该区域下沉时间序列图，如图4所示。由图4可知，从宏观上直观反映出该采油厂在时间和空间上的下沉演变，该区域下沉是一个缓慢发展的过程，随着时间的推移，下沉区域累积下沉量不断增大，下沉范围也在持续扩大。

图4 下沉时间序列图Fig.4 Subsidence time series

根据监测结果可知该下沉区域形状近似椭圆形(图5)，对PS-InSAR监测的矢量成果插值获取下沉区域沿图5中A线和B线的方向的下沉剖面线，分别如图6和图7所示，A线和B线分别视作椭圆形下沉区域的长轴和短轴。根据A线和B两条剖面线数据进一步反映的下沉情况可以看到A线和B线的下沉量关于剖面线上的最大下沉量处对称，可以判断下沉盆地的中心位于A线和B线的交点周围。

研究区域位于辽河盆地，自更新世以来，整个盆地以下沉为主，呈北东向展布，与凹陷区的构造较为一致，在辽阳、鞍山一带下沉速度为-3 mm/a，在临近渤海的盘山，最大下沉速度为-5 mm/a[24]。根据国家重点基础理论研究(攀登)项目研究报告[25]，中国学者对中国东南沿海地区地壳与陆地垂直运动的研究，其中辽宁营口地区辽河平原地壳与陆地平均下沉速度为3 mm/a左右。根据上述文献，研究区域陆地自然下沉量为-3～-5 mm/a，所以3年累计最大下沉为-15 mm，因此在PS-InSAR监测的下沉数据中，包含着石油开采、地下水开采及陆地自然下沉等复杂影响因素。但在约3年的监测时间跨度内，PS-InSAR技术获得的地表最大下沉超过-580 mm，自然下沉量仅占2.6%，自然下沉相对开采引起的下沉量较小，且由图6和图7可以看出，地表下沉是连续的，没有出现非连续台阶状下沉，说明构造对地表下沉的影响没有显现出来，因无法区分开采下沉与自然下沉，因此，把地表下沉监测结果视作开采引起的，忽略了自然下沉等次要因素影响。

图5 下沉盆地区域及剖面线位置示意图Fig.5 Schematic diagram of subsidence basin areaand section line position

图6 下沉盆地长轴剖面线Fig.6 Long axis section line of subsidence basin

4 地面移动变形模型反演及下沉

本文中以下沉量为10 mm的外部边界区域作为下沉盆地边界以求取椭圆下沉盆地长轴和短轴的长度，以A、B两条剖面线的数据为已知量对地表下沉模型进行反演以确定模型参数。根据下沉盆地剖面线的数据可以确定椭圆下沉盆地的长半轴a和短半轴b的长度分别为3 350 m、2 050 m，最终模型反演，见式(9)。

(9)

式中，w(x,y)为累积下沉量，m。

通过利用反演得到的模型对下沉盆地两条剖面线位置的累积下沉量进行正演，比较模型拟合精度，结果如图8和图9所示。

图7 下沉盆地短轴剖面线Fig.7 Short axis section line of subsidence basin

图8 下沉盆地长轴监测结果与正演结果比较Fig.8 Comparison of long axis monitoring resultsand forward results in subsidence basin

图9 下沉盆地短轴监测结果与正演结果比较Fig.9 Comparison of short axis monitoring resultsand forward results in subsidence basin

该模型的相关系数为0.996 5，正演结果的均方根误差为23.4 mm，该模型对该下沉盆地累积下沉量的正演结果具有较高的精度。

基于该下沉模型及地表移动变形相关理论可分别计算得到该下沉区域的倾斜i(x,y,φ)、曲率K(x,y,φ)、水平移动U(x,y,φ)、水平变形ε(x,y,φ)表达式分别为式(10)、式(11)、式(12)和式(13)。

(10)

(11)

结合该区域地质情况和监测数据，取B=0.25，根据主要影响半径定义计算得到r=1 985 m。

U(x,y,φ)=Bri(x,y,φ)=

(12)

ε(x,y,φ)=BrK(x,y,φ)=

(13)

基于上述得到的该采油厂区域地表移动变形模型可计算得到椭圆下沉盆地各点的倾斜、曲率、水平移动和水平变形值。下沉盆地长轴A上各点沿φ=0°方向和短轴B上各点沿φ=90°的倾斜和水平变形分别如图10和图11所示。

图10 A轴倾斜和水平变形Fig.10 A-axis tilt and horizontal deformation

图11 B轴倾斜和水平变形Fig.11 B-axis tilt and horizontal deformation

5 结 论

1) 本文采用PS-InSAR技术对某油田进行了下沉监测，通过对监测结果分析发现：该油田的下沉演变是一个由小到大，不断增加的过程。在监测时间段内，该采油厂监测范围内共提取到17 796个PS点，下沉区域面积达29.5 km2，PS-InSAR技术为对该区域下沉过程监测提供了高效的技术方法。

2) 应用椭圆弹性薄板挠度模型，建立该区域地表移动和变形模型，并结合实测结果反演得到模型参数，模型计算均方根误差为23.4 mm，模型的相关系数为0.996 5，具有较高的拟合精度，所建立的计算模型可以作为地表下沉变形评价的基准。

3) 截至2019年11月，监测区域地表最大下沉585.6 mm，最大下沉速率为-210.9 mm/a。根据反演模型计算，地表最大倾斜0.47 mm/m，最大水平拉伸为0.42 mm/m，最大压缩变形分别为-0.34 mm/m。计算结果为评价该区域开采损害风险评估提供依据。

4) 目前关于石油开采引起的地表下沉规律研究的还不充分，应通过对地表下沉监测获取更多的地表下沉信息与储层变化的相互关系，建立更加精确的地表下沉与变形计算模型，为评价开采对地面设施与环境的影响提供依据。