克拉玛依油田SBAS-DInSAR 长时序形变监测

于 冰 ，佘 杰，赵金洲，张 过，周志伟

1.西南石油大学土木工程与测绘学院，四川 成都 610500；2.油气藏地质及开发工程国家重点实验室·西南石油大学，四川 成都 610500；3.西南石油大学油气空间信息工程研究所，四川 成都 610500；4.武汉大学测绘遥感信息工程国家重点实验室，湖北 武汉 430079；5.中国科学院精密测量科学与技术创新研究院大地测量与地球动力学国家重点实验室，湖北 武汉 430077；6.绵阳市安州区水利发展中心，四川 绵阳 622650

引言

地表形变是指在自然或人为因素影响下地球表面产生形变的一种普遍自然灾害，产生地表形变的因素包括地下资源（地下水、石油、天然气、固体矿藏等）开采、地质运动、工程建设等，其对人类的生命财产安全和生产活动安全有着重要的影响[1-5]。克拉玛依油田作为中国重要的石油开采和生产基地，经过多年开采及注水注气等措施，地下应力发生改变，必然发生显著形变，如何高效地监测克拉玛依油田的复杂形变具有重要的现实意义。

合成孔径雷达干涉测量（Interferometric Synthetic Aperture Rader，InSAR）作为一种高新的遥感技术，被广泛应用于大地测量[6]。随之发展而来的差分合成孔径雷达干涉（Differential InSAR，DIn-SAR）技术通过引入外部数字高程模型（Digital Elevation Model，DEM）进行地形相位差分，可实现地表形变监测，但是DInSAR 容易受时空失相干、地形误差、大气效应、轨道误差等的影响[7-10]。

鉴于此，国内外学者提出了时序差分雷达干涉技术（Time Series DInSAR，TS-DInSAR），也被称为多时相DInSAR 技术（Multi-temporal DIn-SAR，MT-DInSAR），其基本思想是将相同区域的多幅SAR 影像进行联合处理，去除地形、大气、轨道等误差，获取高精度的地表形变信 息[3，11-12]。常见的方法包括永久散射体干涉（Persistent Scatterer DInSAR，PS-DInSAR）、小基线集干涉（Small Baseline Subsets DInSAR，SBASDInSAR）、斯坦福永久散射体干涉（Stanford Method for Persistent Scatterer，StaMPS）、时域相干点差分雷达干涉（Temporarily Coherent Point DInSAR，TCPDInSAR）等方法[13-17]。

克拉玛依油田于1985 年开始进行注水开发[18]，储层应力变化较为复杂，可能引发复杂地表形变。Ji 等[19]利用2010--2016 年的ALOS PALSAR 数据对克拉玛依Hei103 区进行形变监测分析，得到该区域的抬升是由地下流体注入所致；Aimaiti 等[20]利用2007--2009 年的两景ALOS PALSAR 数据进行DInSAR 处理，得到两处明显形变区，再运用20032010 年的21 景ENVISAT ASAR 数据进行PS–DInSAR 和SBAS–DInSAR 处理，比较了3 种方法的相关性并对形变速率进行了分析；Yang 等[21]利用20172018 年的Sentinel–1升降轨数据，采用MSBAS-InSAR 方法研究了形变的变化趋势，运用Moji 和Sill 模型对沉降漏斗进行建模，并且得到该区域地震与油田开采及注水注气无关的结论。

由于以上研究所运用的影像数量有限，无法获取克拉玛依油田精确的长时序地表形变发展规律，本文利用2016–01–122020–05–14（时间跨度1584 d）的110 景哨兵一号（Sentinel–1）影像，通过平均形变速率和长时序形变监测，准确地分析了克拉玛依油田在空间和时间上的形变发展规律，为InSAR 技术用于油田广域长时序非线性形变监测提供了一种较好的参考方案。

1 研究区与数据

1.1 研究区概况

克拉玛依市位于准噶尔盆地西部，由克拉玛依区、独山子区、白碱滩区和乌尔禾区组成，是国家重要的石油石化工业基地。克拉玛依市海拔高度在270～500 m，由西北向东南方向逐渐降低，气候为典型的温带大陆性气候，植被覆盖率低。

克拉玛依油田分布与克拉玛依市区范围基本一致，克拉玛依油田为断裂遮挡的单斜油藏，岩性为砂岩、砾岩、泥岩交互沉积，表现为低孔隙度和低渗透性，具有粗岩性、高非均匀程度的特点[22]，该油田目前已进入高含水开发阶段[18]。研究区概况如图1所示。

图1 研究区位置分布Fig.1 Location of the study area

克拉玛依油田作为中国重要的油气开采及石化产业基地，对其进行形变监测，为油田安全运营提供参考意义重大，本次研究选取克拉玛依区与白碱滩区及周边地区为研究范围，研究区内包含大量的人工建、构筑物，是克拉玛依油田主要的石油开采、炼制生产区，同时也是克拉玛依市主要的人类活动聚集区域，克拉玛依断裂沿东北至西南方向穿过研究区。

1.2 数据集

研究采用干涉宽幅（Interferometric Wide Swath，IW）成像模式下的VV 极化数据，分辨率为5 m×20 m。采用欧洲空间局Sentinel–1 Quality Control提供的精密轨道数据进行精轨校正。差分干涉所需的数字高程模型采用日本航天局提供的AW3D30产品。

2 数据处理原理与方法

2.1 SBAS-DInSAR 原理

常用的TS-DInSAR 技术包括PS-InSAR 及SBAS-DInSAR，考虑到本次研究时间跨度长达1 584 d，并且研究区地面多为裸地，采用SBASDInSAR 能够有效地克服时空失相干影响，识别更多的高相干目标。

SBAS-DInSAR 是由Berardino[14]于2002 年提出的一种TS-DInSAR 处理方法。该方法通过设置时空基线的阈值，将SAR 影像分为若干个小集合，降低了时空失相干的影响，提高了整体的相干性。SBAS-DInSAR 采用奇异值分解法（Singular Value Decomposition，SVD）将小集合之间的形变信息进行联合计算。

采用影像数量，通过设置时空基线阈值，可获得干涉对，其关系如下

式中：

N—影像数量，景；M—干涉对数量，对。

SBAS-DInSAR 生成干涉对后，需经过滤波及相位解缠等处理。根据干涉图及解缠图判定干涉对质量，保留质量较好的干涉对，舍去质量差的干涉对[24]。针对任意t1和t2时刻组成的干涉对，方位向为x、距离向为r的像元处的差分干涉相位可以表示为

式中：

λ—雷达的波长，mm；

Δddef—两景影像成像期间地表发生的形变，mm；

B⊥—垂直基线长度，mm；

Δh—地形误差，mm；

θ—雷达入射角，（°）；

Δdatm—两景影像大气引起的距离差；

Δn—轨道误差、系统误差等引起的噪声相位。

式（2）约等号右边第一项表示形变相位，第二项表示地形误差引起的相位，第三项表示大气延迟相位，第四项为噪声相位。为了获取精确的形变信息，需要将地形误差相位、大气延迟相位以及噪声相位去除。当去除以上相位后，可用矩阵形式表示所有干涉对之间的差分干涉相位

式中：AAA—M×N的系数矩阵。

任一干涉对之间的形变差Δddef可以用该时间内平均形变速率vt2-t1与时间差t2-t1表示如下

基于式（4），可将式（3）变换为

式中：BBB M×N的系数矩阵，vvvT可表示为

当矩阵BBB为满秩时，可以采用最小二乘法计算形变相位，但是这种情况的概率较小，所以SBASDInSAR 在进行求解时引入奇异值分解法进行形变解算，获取最小范数下的最小二乘解。

2.2 SBAS-DInSAR 非线性形变提取流

SBAS-DInSAR 通过设定时空基线阈值，能够有效地解决时空基线限制问题。同时，由于其采用的分段线性模型以及时空滤波提取非线性形变的方法，对于提取长时间跨度且非线性的形变有更好的适用性。

SBAS-DInSAR 的具体处理流程如图2 所示。

图2 SBAS-DInSAR 流程图Fig.2 Flow chart of SBAS-DInSAR

（1）根据时空基线生成连接图：本次研究使用110 景影像，时间间隔为1 584 d。选择2018–03–26影像为主影像进行配准，再根据基线组合共生成440 个干涉对，其中，空间基线最长为91.77 m，时间基线最长为72 d。

（2）干涉数据处理：包括生成干涉图、结合外部DEM 去除地形相位、滤波、相干系数计算、相位解缠。其中，采用Goldstein 方法对干涉图进行滤波，使用最小费用流法（Minimum Cost Flow，MCF）进行相位解缠。

（3）剔除干涉对：受影像质量以及时空基线的影响，部分干涉对效果较差，表现为相干性低以及解缠错误等问题，可通过干涉图、相干系数图以及解缠图进行筛选。本次研究设置最低相干系数阈值为0.3，当相干系数低于0.3 的面积占比超过研究区面积的50%时，剔除该干涉对；同时检查相位解缠图，当解缠存在大面积跳变时剔除该干涉对。共剔除50 对干涉对，剔除比例为11.36%。部分干涉对及其结果如图3 所示。剔除效果不好的干涉对可有效避免误差，提高数据处理的效果。图3 中a 与c表现为相干性较低，并且解缠效果差；图3 中b 干涉条纹较混乱，并且解缠结果中有明显的跳变；图3中d 与e 干涉对效果好。相干性较低的干涉对以及存在跳变的解缠相位将影响结果准确性，应被剔除。基于以上规则，应去除与图3 中a、b 和c 效果相近的干涉对。去除效果较差的干涉对后的时空基线如图4所示。

图3 部分干涉对及相应的干涉处理结果Fig.3 Partial interferometric pairs and the relevant interferometry results

图4 干涉对的时空基线图Fig.4 Spatial-temporal baseline of interferometric pairs

（4）估算形变速率及高程系数：SBAS-DInSAR根据相干性阈值提取高相干目标点，对提取的高相干目标进行模型构建，转换为矩阵形式后不满秩，使用奇异值分解法解算低频形变和高程误差（地形残差）分量。

（5）去除地形残差相位及大气相位：去除残差相位并重新解缠，使用时间域高通滤波以及空间域低通滤波去除大气效应，并获取形变时间序列，进而建模提取形变速率。

（6）地理编码：将SAR 坐标系下的结果进行地理编码转换到地理坐标系中。

3 实验结果与分析

3.1 平均形变速率分析

SBAS-DInSAR 通过对识别的高相干目标进行时序形变建模，能够得到研究区内的平均形变速率，雷达视线（Line of sight，LOS）向和垂直向的平均形变速率能够直观地反映出研究区整体的形变趋势与特点。

3.1.1 LOS 向及垂直向形变分析

利用SBAS-DInSAR 提取研究区LOS 向和垂直向平均形变速率结果如图5a 与图5b 所示。其中，LOS 向形变速率在— 48.1～60.2 mm/a，垂直向形变速率在— 59.3～73.4 mm/a。根据图5 可以得出，形变主要集中于研究区内SW—NE 方向，形变区分布与克拉玛依断裂走势相近，并且与石油开采区域空间分布有高度的一致性。

图5 平均形变速率Fig.5 Average rate of deformation

3.1.2 高形变区提取

高形变区往往对应着快速开采区，也是石油产量较大的区域，因此，单独提取高形变区并进行分析对重点开采区的运营安全监测具有较好的参考价值。标准差作为一种最常使用的数据分布指标，表示数据的稳定性以及数据相对于平均值的偏离程度，通过对平均值和标准差的组合，能够反映数据间的变异性[25]。定义形变超过平均值减两倍标准差的范围为高沉降区

式中：vhsa—高沉降区的形变速率临界值，mm/a；

K—SBAS-DInSAR 识别的所有相干目标点的数量。

定义形变超过平均值加上两倍标准差的范围为高抬升区

式中：vhua—高抬升区的形变速率临界值，mm/a。

提取垂直向平均形变速率图的高沉降区和高抬升区为高形变区，结果如图6 所示。

图6 高形变区Fig.6 High deformation area

由图6 可以看出，研究区内高形变区可以划为6 个主要区域（A—F）。其中，沉降速率最高为区域B 的B2 点，抬升速率最高为区域A 的A3 点。A 区域位于克拉玛依市区西南方向，靠近克拉玛依断裂南侧，形变以抬升为主，表现为小范围剧烈抬升，其形变原因可能与该区域油田注水注气相关；B 区域位于克拉玛依市区东部与金龙镇西北部之间，克拉玛依断裂穿过区域B，B3、B5 表现为抬升，B1、B2、B4、B6 表现为沉降，其中，B4 为龙山滑雪场，其沉降可能是由于建造后形成的地表下沉；区域C 位于金龙镇东北方向与白碱滩区西部之间，克拉玛依断裂穿过区域C，是克拉玛依油田采油三厂所在地，地表密集分布抽油机，仅C5 表现为沉降，其余均表现为抬升，抬升可能与油田注水注气增产相关，沉降可能是由于地下压力挤压所致；区域D 位于白碱滩区西北部，紧邻克拉玛依断裂北侧，表现为抬升，其抬升可能与油田注水注气有关；区域E 位于白碱滩东部，国道217 线南侧与奎北铁路北侧之间，分布与克拉玛依断裂走势一致，该区域为克拉玛依油田的主要开采区，根据雷达强度影像及可见光影像可知该区域采油机分布密集，除E3 外均表现为沉降，并且其沉降区域与采油机分布一致；区域F 位于研究区东南方，奎北铁路南侧，F1—F4 表现为沉降，F5—F10 表现为抬升，其中，F1 沉降速率最高且为人工建筑物，其形变原因需实地确认，F2—F6靠近大拐苇湖，该区域形变可能与盐田生产有关，F7—F10 表现为抬升，由于缺少资料，其形变原因需实地调查。

3.1.3 形变梯度计算

形变梯度通过计算某个方向上的变化速率，能够直观地体现形变在空间上的非均匀性。通过提取研究区垂直向形变速率由西向东和由北到南的梯度，分析研究区的非均匀性形变。由于SBASDInSAR 只能提取研究区内高相干点的形变，为了计算形变梯度，采用反距离权重法对垂直向形变速率进行插值，使用gradient 函数进行梯度计算

式中：GW-E(i,j)—像元（i,j）由西向东的形变梯度，mm/（a·像元）；

GN-S(i,j)—像元（i,j）由北向南的形变梯度，mm（/a·像元）；

v（i,j）— 像元（i,j）的形变速率，mm/a。

插值后的垂直向形变速率以m×n矩阵的形式存储，利用式（9）与式（10）计算由西向东和由北向南的形变梯度，结果如图7 所示。

图7a 表示由西向东的形变梯度，数值在 -1.11～1.01 mm（/a·像元）；图7b 表示由北向南的形变梯度，数值在 -1.23～1.32 mm（/a· 像元），像元大小为11.920 m×11.920 m。形变梯度的绝对值大小能够反映区域内形变空间非均匀性的剧烈程度，绝对值越大，代表该点在该方向上形变速率的变化量越大。研究区由西向东的形变梯度绝对值最大点如图7a红色三角形所示，该区域位于克拉玛依市与金龙镇之间并且表现为抬升；研究区由北向南的形变梯度绝对值最大点如图7b 蓝色三角形所示，并且该区域平均抬升速率也为研究区最大。

图7 垂直向平均形变速率的形变梯度Fig.7 Deformation gradient of vertical average deformation rate

3.2 长时序非线性形变分析

本次研究时间跨度长达1 584 d，且研究区内石油开采及注水注气等活动剧烈，导致研究区内的形变在时间维上表现出显著的非线性，在空间维上表现为非均匀性，为了更直观地分析研究区内长时序非线性形变，本文进一步提取了不同时段及特征点的形变，并分析高形变区的时序形变特征。

3.2.1 等时间间隔累计形变计算

以96 d 为时间间隔，制作16 幅等时间间隔的垂直向累计形变图，并且提取每个时间段内最大抬升点以及最大沉降点位置，结果如图8 所示。可以看出，研究区等时间间隔的形变量是不同的，其中，[384,480）d 时段沉降量最大为 -59.7 mm；在[960,1 056）d 时段抬升量最大为89.4 mm。图8 中红色三角形代表该时间段内最大沉降点，蓝色三角形代表该时间段内最大抬升点。值得注意的是，由图8可以看出，不同时段的最大沉降点和最大抬升点是在变化的。最大沉降点主要分布于区域B 和区域E；最大抬升点主要分布于区域A 和区域D，并且在768 d 之前主要位于D 区，在768 d 之后主要位于A区。通过分析等时间间隔的形变图可以得出，研究区在不同时间段形变发生的主要区域是不同的，且形变速率也是不同的，这与克拉玛依油田开采及注水注气活动有着密切联系。

图8 等时间间隔形变图Fig.8 Ground deformation map with equal time interval

3.2.2 高形变区时序分析

鉴于研究区内形变的复杂性，为了更好地分析研究区内形变总体变化趋势，利用式（7）、式（8）提取累计形变的高形变区，计算其所占研究区面积百分比，制作高形变区百分比时间序列图，其结果如图9 所示。

值得注意的是，在计算高形变区时，为了计算结果的准确性，并未对累计形变结果进行插值处理。研究区内部有30.96%区域的相干性低于所设阈值，图9 中所指百分比为识别出的高形变区占全研究区面积的百分比，不包括由于相干性低于阈值未识别出的高形变区。

图9 高形变区百分比时间序列图Fig.9 Time series diagram of percentage of the high deformation area

如图9 所示，研究区内的高形变区由高沉降区和高抬升区组成。在[0,336）d 内，高沉降区面积变化不大；在[336,432）d 内，高沉降区面积所占比由1.22%迅速增加到2.17%；在[432,1 584）d 内，研究区内高沉降区面积基本表现为降低。在[0,336）d内，高抬升区比例变化不大；在[336,516）d 内，高抬升区所占比例逐渐减小至最小值1.12%；在[516,1 584）d 内，研究区内高抬升区域所占比例逐渐升高，最后达到1.88%。将高沉降区与高抬升区组合定义为高形变区，在[0,660）d 内，高形变区比例高达3.0%以上；在[660,1 164）d 内，高形变区面积表现为降低趋势；在[1 164,1 584]d 内，研究区高形变区表现为高抬升趋势。通过对高形变区所占比例进行分析，能够反映出研究区总体形变趋势仍是变化的，进一步表明了克拉玛依油田形变的时空特征复杂性。

3.2.3 特征点时序分析

分析研究区等时间间隔形变图和高形变区百分比时间序列图可以看出，研究区的形变区域与形变速率是随着时间变化的，仅依靠以上分析难以确定研究区形变具体的时序变化规律，通过分析特征点的垂直向累计形变曲线可以直观地反映研究区内的形变发展特点，对图6 中AF 区域内共42 个特征点进行垂直向累计形变分析，结果如图10 所示。

图10 特征点垂直向累计形变曲线Fig.10 Vertical deformation curves of characteristic points

区域A 中A2、A3、A4 点均为形变中心，在空间上表现为邻接，且形变趋势较一致，在[700,1 200）d剧烈抬升后缓慢沉降；A1 与A5 前期表现为震荡沉降，后期表现为抬升，且抬升速率逐步升高。区域B 中除B5 外形变均接近线性，B5 表现为持续抬升，并且在[800,1 000）d 迅速抬升，该区域为油田开采区域，判断B5 的非线性抬升可能与该区域油田注水有关。区域C 位于克拉玛依采油三厂周边，该区域表现为大范围抬升且伴随着小范围沉降，是研究区中形变最密集的区域。通过分析其特征点的垂直向累计形变曲线可以得出该区域抬升速率较一致，判断其抬升与油田注水注气相关，且C5 在末期有抬升趋势。区域D 位于白碱滩区东北方向，有多处明显的抬升中心，其中，D2 的形变最剧烈；D1、D2在[800,1 000）d 表现为沉降，其余时间段表现为抬升，D3 整体趋势表现为抬升，但是在900 d 之后形变表现为震荡。区域E 为克拉玛依油田的主要开采区域，地面密布抽油机，E1 与E2 前期表现为抬升，后期表现为持续沉降；E3 前期表现为沉降，于600 d 后迅速抬升并在900 d 后表现为相对稳定，其具体形变原因需结合油田开采及注水注气资料进行分析。区域F 中F2、F3、F4 在500 d 前表现为沉降，500 d 后相对稳定。

4 结论

（1）研究区形变主要集中在SW--NE 方向，LOS 向平均形变速率在-48.1～60.2 mm/a，垂直向平均形变速率在 -59.3～73.4 mm/a，形变区分布与克拉玛依断裂走势相近，并且与石油开采区域有高度的一致性。

（2）根据均值和二倍标准差的组合提取了垂直向平均形变速率对应的高形变区，共划分了6 个形变区域，结合地表地物类型分析了高形变区的产生原因多与油田注水注气及油气开采相关。

（3）通过分析研究区内由西向东和由北向南的垂直向平均形变速率的形变梯度可以得出，研究区内的形变在空间上表现出显著的不均匀性。

（4）分析等时间间隔的垂直向累计形变结果，可以得出研究区等时间间隔的形变量在数值上是不同的，同时等时间间隔内的最大沉降点和最大抬升点也在随着时间而发生转移。

（5）根据高形变区百分比时间序列图可以得出，研究区内高沉降区所占比例先上升后下降，高抬升区所占比例先下降后上升，高形变区所占比例先下降后上升，能反映出研究区总体形变趋势仍是变化的。

（6）对研究区内42 个形变特征点进行垂直向累计形变分析，可以得出研究区的形变是非均匀、非线性的，特征点的累计形变量随着时间的推移表现为多样性，与地面地物类型有着相关性。

（7）本研究能够有效证明InSAR 技术在油田长时序非线性形变监测中的可行性，提供了一种高效低成本的油田形变监测方法，可为避免油田形变灾害及保障油田安全运营提供可靠的技术保障。