天然气作为一种相对清洁的能源，其燃烧时产生的二氧化碳和有害物质少于煤炭、石油等传统能源，但天然气长期开采会引起气田地下压力减小，导致地表沉降甚至诱发地震^[1]。传统的气田地表沉降监测方法包括水准测量和导航卫星测量，这2种技术受到观测条件和成本限制，只能提供离散站点的形变信息，无法全面揭示气田地表沉降状态。InSAR技术具有空间分辨率高、覆盖范围大、可重复观测和无需地面设备等优点，被广泛应用于地下水、油田、矿区、火山等由地下物质变化引起的地表形变研究^[2]。

涩北气田是青海省最大的天然气田，也是我国西部地区重要的天然气基地之一。经过几十年开采，已经引发严重的地表沉降，可能导致采气设备、输气管道和周围建筑物的损坏。因此，需要对涩北气田地表进行定期监测，并根据监测结果采取必要措施，确保设备安全运行。

为监测涩北气田天然气开采引起的地表形变，本文使用时序InSAR技术，基于退役和在轨的多颗SAR卫星数据，恢复涩北气田近22 a开采过程的地表形变特征，并讨论天然气开采导致地表发生形变的机理。

1 研究区概况

柴达木盆地蕴含丰富的油气资源，目前已发现的天然气储量主要集中在盆地东部的三湖地区^[3]。区内目前已探明涩北1号、涩北2号和台南等3个整装大型第四系生物气田及驼峰山、台吉乃尔、伊克雅乌汝等含气构造^[4]，具体分布如图 1所示。

涩北气田开采范围内地层主要为第四系湖泊相沉积，储层岩性主要为泥质粉砂岩，其次为粉砂质泥岩、含粉砂泥岩和泥岩、粉砂岩^[6]。涩北气田主要由涩北1号、涩北2号和台南3大气田组成，其中涩北1号气田发现于1964年，于1995年试采开发，1996年柴达木盆地东部首条天然气长输管道投入运营，标志着涩北天然气田正式开发。由于持续的大规模天然气开采，涩北天然气田已产生显著的地表形变，这对气田开采设备和天然气管道的安全运行形成巨大威胁。

2 InSAR数据处理与形变分析 2.1 SAR数据与InSAR预处理

为研究涩北气田天然气开采引起的地表形变，本文收集5个轨道的ERS-1/2、ENVISAT、ALOS和Sentinel-1A/B等多颗卫星SAR数据，时间跨度为1996~2018年，共168景数据，详细信息见表 1。

由于研究区气候干燥、植被覆盖非常少，2幅相隔数年的SAR影像仍有可能保持相干性。基于此，本文在选取InSAR干涉对时，基于短空间基线和长时间基线原则生成1 030幅干涉图，从中选取597幅高质量干涉图进行后续时间序列分析(图 2)。使用DEOS/DORIS、OPOD等精密轨道数据精化SAR数据的轨道参数，使用30 m分辨率的SRTM DEM数据模拟并去除地形相位。为提高干涉图的相干性，对各轨道的干涉图进行多视处理，使干涉图分辨率与DEM分辨率相同。为进一步提高干涉图的信噪比，对干涉图进行自适应空间域滤波处理，相位解缠使用基于Delaunay三角网的最小费用流算法。

2.2 时序InSAR处理

为恢复涩北气田自正式开发以来地表形变随时间的演化特征，本文采用小基线集技术(small baseline subset, SBAS)求解每个轨道数据的形变时间序列^[7]。SBAS算法通过设置较小的时空基线阈值获取高质量干涉图，以降低时空失相干和DEM误差的影响。此外，通过时间域高通滤波和空间域低通滤波来估计轨道误差和大气延迟等空间相关误差，从而获取高精度的形变速率图(图 3)和形变时间序列(图 4)。

SAR侧视成像仅能获取地表形变在各轨道雷达视线向的投影信息，而联合同一时期的升轨和降轨SAR数据，不仅能提高监测数据的时间分辨率，也可为获取监测区的多维形变信息提供基础数据。考虑到InSAR技术对南北方向形变不敏感，本文采用多维小基线集技术(multidimensional small baseline subset, MSBAS)^[8]对2014年以来收集的升降轨Sentinel-1数据进行东西方向和垂直方向的二维形变分解。MSBAS技术是在SBAS-InSAR技术基础上发展而来，利用在时间和空间上重叠的不同轨道升降轨SAR影像数据，实现多视角影像时间分辨率的融合，解算研究区在垂直向和东西向的二维形变速率及累积形变时间序列。二维时间序列分解公式如下：

$ \left(\begin{array}{c} \boldsymbol{A} \boldsymbol{S} \\ \lambda \boldsymbol{L} \end{array}\right)\left(\begin{array}{l} \boldsymbol{V}_E \\ \boldsymbol{V}_U \end{array}\right)=\left(\begin{array}{l} \boldsymbol{D} \\ \bf{0} \end{array}\right) $

(1)

式中，A为SAR影像的时间间隔矩阵；S={S_E, S_U}={－cosθsinϕ, cosϕ}为投影矩阵，θ为方位向角度，ϕ为视线入射角；λ为正则化参数；V_E和V_U为待求的EW向和UD向分量；D为D-InSAR技术获取的视线向观测值；L为正则化系数矩阵。本文利用MSBAS技术成功获取涩北气田地区二维形变速率(图 5)和形变时间序列(图 6)，可为后续的气田开采时间序列参数反演提供形变数据。

2.3 InSAR形变结果分析

利用1996~2018年5个轨道的SAR数据获取涩北3大气田的形变速率图(图 3)和形变时间序列结果。由图 3可见，天然气开采产生的地表形变主要分布在3大气田开采区内，其空间分布与已探明的涩北气田气藏空间分布(图 1)相似。T362、T269和T490轨道形变结果显示，1996~2010年涩北1号气田是3大气田中形变量和形变范围最大的地区。而T70和T77轨道形变结果显示，2014~2018年台南气田已经成为该地区形变量和形变范围最大的地区，说明从2014年开始，台南气田开采速度明显增加。从5个轨道的形变速率图可以看出，涩北3大气田形变场位置分布存在一定差异，即形变场在不同SAR成像几何中的形变位置不同，说明天然气开采产生的地表形变存在明显的水平分量，而不单是垂直方向形变。

为分析涩北1号气田从正式开发至2018年约22 a的地表形变信息，提取涩北1号气田沉降中心A点的形变时间序列结果(图 4)。InSAR监测结果显示，2003年之前A点基本保持稳定，仅有微弱隆升；2003~2005年A点隆升近6 cm，由于这2 a多时间缺少观测数据，尚不清楚其具体形变过程；2005年A点开始快速下沉，至2018年仍保持下沉趋势，且下沉速率明显增加。

由图 5可见，垂直向形变量约占雷达视线向形变量的80%。由图 6可见，涩北1号气田地表形变中心在垂直方向和东西方向基本表现为线性形变，且垂直方向形变量远大于水平方向形变量。

3 大地测量数据建模

为分析涩北1号气田因天然气开采导致地表产生形变的机理，本文基于均匀弹性半空间位错模型，采用GBIS软件^[9]分析涩北1号气田在2015~2018年间形变源变化过程。该软件采用马尔科夫链蒙特卡洛(MCMC)方法，结合Metropolis-Hasting算法对收敛于贝叶斯后验分布的模型参数进行随机抽样，计算每个参数概率密度函数的估计值，并提供参数在95%置信区间的最优解。

本文选择GBIS软件中Sill模型，利用垂直方向和东西方向形变时间序列结果，计算涩北1号气田在2015~2018年间形变源的长度、宽度、深度、走向、参考点位置(X，Y)和张量等参数。图 7为利用升降轨Sentinel-1卫星垂直方向和东西方向累积形变反演得到的2018-03-11模拟形变和残差形变，图 8为反演得到的涩北1号气田形变源深度和体积随时间变化曲线。

由图 7可见，Sentinel-1卫星垂直方向和东西方向累积形变反演结果的模拟残差均很小，说明该反演方法较为可靠。由图 8可见，形变源的天然气地下体积在2015~2018年逐渐增大，即天然气地下开采体积呈近似线性增加趋势。形变源深度变化则分为两个阶段：2015-01~2016-08为台阶式变化，开采深度增加至2 100 m；2016-08以后开采深度逐步稳定在2 300m左右。

4 讨论 4.1 2015~2017年涩北1号天然气开采量估计

本文通过InSAR形变结果反演得到2015-01~2018-03涩北1号天然气田地下气体体积减少约6 000万m³，2015~2017年地下气体体积分别减少约1 646万m³、1 635万m³和19.63万m³。根据2016~2018年《中国石油青海油田公司年鉴》可知，2015~2017年涩北1号气田实际气体产量分别为13.72亿m³、13.07亿m³和15.98亿m^3[10-12]。在相同时间段内，InSAR反演得到的地下气体体积与实际气体产量之间存在巨大差异，这归因于不同环境条件(如温度和压力)下气体体积的变化。本文估计地下气体体积变化的条件为涩北1号气田2015~2017年平均压力为920 kPa、平均温度为319.26 K^[10-12]；而《中国石油青海油田公司年鉴》统计结果则是基于中国标准状态，即压力为101.325 kPa、温度为293.15 K。为增加评估结果的可靠性，根据理想气体状态方程，将地下条件下的气体体积变化转换到中国标准状态。转换后的气体体积变化结果为13.72亿m³、13.07亿m³和15.98亿m³，由此可见，InSAR反演得到的天然气体积与《中国石油青海油田公司年鉴》统计结果的差异极小，可以定量评估InSAR反演方法在监测地下气体体积变化方面的可靠性。

此外，以下因素也会导致InSAR估计的天然气开采量与《中国石油青海油田公司年鉴》统计结果之间存在差异：1)本文使用均匀反演模型，反演过程中仅采用单一形变源的平均深度和平均张量计算地下气体体积，由于气藏区层数多，层深和层厚不等，因此会产生一定误差；2)天然气开采通常伴随有少量水和砂，会造成一定差异；3)根据孔隙弹性理论，地下气体压力降低可能会抵消一定的地下体积变化^[13]。

4.2 天然气开采过程中的地表形变模式

由于天然气的大量开采，地下物质逐渐减少，导致地表产生明显形变。图 9为天然气开采剖面示意图，从中可看出天然气开采与地表形变之间的关系，红线表示由于地下天然气开采产生的正断层，黄色箭头表示地表水平形变方向。由于地下物质开采导致地面下沉，形变中心周围区域地表也朝着形变中心发生水平位移。从图 5和图 6可以看出，涩北3大气田形变区在快速下沉的同时还伴随有明显的水平形变。

由于SAR为侧视成像，这会导致单一轨道InSAR观测对水平形变的敏感度较低。本文利用同一时期的升降轨SAR数据进行联合解算，获得研究区二维(垂直方向和东西方向)形变时间序列(图 6)，可更全面地展示天然气开采过程中地表形变变化过程。本文选用InSAR二维形变场作为气田形变源反演的观测数据，其反演精度优于基于单一轨道InSAR获得的形变场。

为分析涩北1号气田的开采过程，绘制该气田1997~2013年日产气量曲线图(图 10)。由图可知，2003年以前生产规模较小，2004年开始大规模开采，2007年以后开采速度明显提高。

结合InSAR时间序列形变结果(图 4)可以看出，1996~2003年涩北1号气田地表形变不明显，该时间段内涩北1号气田处于探索试采阶段，开采井数量和产气量较少，因此未产生明显的地表形变。随着涩北1号气田的不断开发，部分老井地层压力逐步下降，涩北1号气田在2003~2006年平均单井日产气量递减明显^[16]。为提高老井天然气产量，通过注气增压方式将底层水排出，这也导致2003~2005年间涩北1号气田地表快速隆升。2005~2007年，涩北气田加快产能建设，引入水平井技术，极大提高了天然气产量。自2005年起，涩北1号气田地表开始出现快速下沉，最大形变区(图 3中A点)在2005~2011年形变速率达到30 mm/a，2014~2018年形变速率更是达到70 mm/a。

5 结语

1) 二维InSAR形变时间序列表明，在天然气开采区域，地表以垂直形变为主、水平形变为辅，而在东台吉乃尔湖地区则表现为单纯的垂直方向地面沉降。

2) 基于弹性位错理论反演的涩北1号气田开采过程显示，开采初期开采量和开采深度成正比。该气田开采深度从2015年的1 500 m逐渐加深至2017年的2 300 m。

3) 根据理想气体状态方程计算InSAR反演的地下天然气储层体积变化在标准状态下的体积，结果显示，InSAR反演结果与实际统计数据差异极小，表明InSAR监测和反演结果具有可靠性。