基于多源卫星数据融合的油气管道区域三维变形监测

余东亮，周广，方迎潮，王明波，薛廉，赵雪

（1.国家管网集团西南管道有限责任公司，成都 610041；2.四川省地质工程勘察院集团有限公司，成都 610032）

滑坡、采空区和不稳定斜坡等地质灾害能够引起地表形变，破坏油气管道周围的环境，威胁油气管道安全，对管道沿线区域进行形变监测为防止地质灾害导致的油气管道损坏，保护区域安全具有重要的意义（Ishwar and Kumar，2017）。

传统地表形变监测通常使用水准仪和GPS 等大地测量方法监测地面下沉和水平移动（何国清，1991），这些方法虽然测量的精度高，但能够监测的范围小，成本高，效率低，工作量较大。而且，传统测量手段都是对离散点进行监测，难以反映大范围的形变趋势，低分辨率的时空监测数据削弱了估计精度和可靠性，降低了地质灾害预测的可靠性。

合成孔径雷达差分干涉测量技术（D-InSAR）能够在各种天气条件下产生高分辨率的远程感知数据，通过数据处理可获得高精度的地表微小形变数据。1989 年，D-InSAR 的概念被首次提出，它通过去除干涉图中的地形相位实现了地表形变和地形数据的分离（Gabriel et al.，1989），从此开启了基于D-InSAR技术的大规模地表形变监测的相关研究（吴立新等，2011），其广泛用于地震、火山运动、滑坡和区域地面沉降等领域（Reeh，et al.，1999；Rignot，2008）。然而，D-InSAR 是基于干涉相位获得地表形变数据，而油气管道区域的地质灾害有可能产生较快的地表形变或较大的形变速率，从而影响了干涉相位数据的相干性，导致数据可靠性下降，影响了D-InSAR 形变观测的精度（Zebker and Villasenor，1992）。

为了克服D-InSAR 技术上的缺陷，PS-InSAR 和SBAS-InSAR 技术应运而生，用于时间序列的地表形变测量（曾琳等，2021）。其中PS-InSAR 技术利用一景SAR 影像为主影像，将主影像与多期从影像差分干涉，基于影像数据幅度信息和相位信息计算散射特征稳定性，对于回波信号强的PS 点，进行差分干涉及相位重构后，获得这些PS 点的时间序列形变。PS-InSAR 技术已被广泛用于滑坡、地震和采矿区等区域的地表形变监测（Abdikan et al.，2014; 朱文峰，2018；薛廉等，2019）。目前油气管道沿线InSAR监测的研究还较少，尤其是基于PS-InSAR 技术对油气管道区域进行三维形变监测的研究。本文重点研究了基于PS-InSAR 技术，获取油气管道区域的PS 点位移，采用空间和时间插值方法计算该区域网格点的位移，并使用多源卫星数据的PS 点位移形变数据，计算获得区域网格点的三维变形，从而监测油气管道地表区域的总体形变，为区域地质灾害预警提供数据基础。

1 BPSTSIM 三维变形监测方法

多源卫星监测数据是指不同平台或同平台不同轨道的InSAR 形变监测数据，例如哨兵升轨数据、哨兵降轨数据以及TerraSAR 数据。根据雷达侧视观测几何可知，InSAR 观测到的是地表的正东、正北和垂直向形变量在雷达视线方向LOS（Line of Sight）的投影和（即矢量和）。单一LOS 向的形变只对沿着这一方向的形变敏感，而对垂直于LOS 方向的形变通常无法监测，单一LOS 向的形变数据不能反映监测区域全方向的形变。因此，为了全面反映地表的三维形变，以及更好的将InSAR 形变信息与管道物联网设备监测到的形变信息进行融合与分析，需要利用多源卫星数据计算获得地表点的三维形变信息。

不同卫星和轨道由于分辨率、入射角等参数的不同，在同一区域内PS 点的分布存在一定差异，由于卫星重访周期、观测时间的差异，获得的PS 点时序形变的时间点也不同，需要对多源卫星数据进行空间对齐和时间对齐。因此，本文提出了BPSTSIM 方法，基于多源卫星PS 点形变数据实现时间和空间数据对齐，然后计算三维形变结果。图1 解释了如何从多源卫星获得的PS 点形变数据，经过逐步计算，最后获得关键点（网格点）三维形变。该方法首先对多源卫星的PS 点形变数据以天为时间间隔进行时域插值，从而能够获得PS 点每天的形变量，同时实现多源形变数据的时间对齐。由于不同数据源处理得到的PS 点不是空间对齐的，因此需要对油气管道监测区域进行网格化。网格化后，需要将所有PS 点的形变转化到临近的网格点形变，然后利用插值方法获得网格点的形变数据，实现多源卫星数据的空间对齐。有了多源卫星网格点的形变数据，就实现了多源卫星数据空间和时间的对齐。最后根据SAR 成像几何特点求解网格点在正东、正北和垂直方向的形变。

图1 BPSTSIM 三维变形监测方法流程图

1.1 PS 点时域插值方法

多源卫星数据是指不同卫星平台或者同一卫星平台不同轨道的SAR 数据，这些卫星数据在同一区域获取的时间并不同步，时间间隔也不同。图2（a）给出了一个PS 点的时间序列形变数据，这些数据并不是以天为单位，而且不同卫星的数据间隔也不同，因此难以实现数据融合，所以首先需要对多源卫星PS 点数据实现时域的插值，获得不同数据来源的所有PS 点以天为单位的形变数据（图2（b））。

图2 PS 点的时域插值

其中t为参变量，为参数，在二维空间它是二维向量。为了确定参数，要求曲线过已知的三个点，并且满足三个条件。即曲线过起始点，当且仅当=0；曲线过终点，当且仅当=1；曲线过，当且仅当=0.5，且切线矢量为−。根据要求的三个条件，显然能够建立三个独立的方程，并求得二次样条曲线为：

则最终得到搭接区间点的值为：

1.2 网格点形变插值方法

时间插值方法能够获得多源卫星PS 点每天的形变值，然而为了能够实现空间点的对齐，需要在监测空间内构建统一的格网，并将各个数据源的PS 点对应的形变值转化到网格点的形变值。整个过程可以分为三个步骤。

第一步：绘制网格。根据卫星监测的PS 点经纬度信息分别查找区域内经纬度的最大值与最小值，然后建立以长宽均为nRange米的格子，并以区域经纬度的最大值与最小值为起点与终点区域绘制网格，所有的PS 点一定位于绘制的网格上或在网格中，同时对PS 点值建树以方便第二步算法中PS 点的查找。建树基本算法如下。

（1）在二维数据集合PS 点中选择具有最大方差的维度，如x 轴方向，然后在x 方向对数据排序，并在该维度上选择其中间值对应的点，以该PS 点为分界点，对PS 点集合进行划分，得到两个子集合和；同时创建一个树结点node，用于存储PS 点数据；

图3 给出了一个建树过程的例子描述。如果有8 个PS 点在二维平面，其分布如图3 所示。依据最大方差维度选择，首先在x方向选择中值，对应的PS点为则为第一个根节点；然后在y 方向以及两个子集合（左右两边形成的两个PS 点集合）找中点，对应的PS 点为和；最后在点形成节点以及叶节点。

图3 PS 点构造树生成过程

第二步：PS 点到网格点转化。在第一步计算的网格基础上，对所有的网格点根据第一步建立的PS树进行网格点周围nRange范围内的PS 点查找。若nRange范围内有PS 点，则计算这些PS 点平均值并赋值给该网格点，同时对该网格点做有效标记（标记为1）；若nRange范围内无PS 点，则该网格点做无效标记（标记为0）。

第三步：插值多源卫星网格点的值。首先在第二步平均值插值的基础上，对所有被有效标记（标记为1）的网格点建树，建树过程同图3，树结构能够优化插值时点的查找速度；然后对于所有的网格点进行第二次插值，若网格点在第二步中被有效标记（标记为1），此时仍然保留原来插入的值不变；若网格点在第二步中标记为无效标记（标记为0），此时采用距离权值法插值，先对所有的无效标记网格点根据先前建立的有效标记网格点树进行N近邻查找，再判断查找到的有效标记网格点是否在有效范围内，若有效范围内存在至少一个查找到的有效标记网格点，对该无效标记网格点做距离权值法插值，计算的值赋给该网格点，并改网格点标记由无效（标记为0）为有效（标记为1）。若有效范围内一个有效标记网格点都没有，该无效标记网格点值仍为0，标记也仍然为无效。

在上述步骤中需要采用距离权值法插值。首先在网格树中查找无效标记网格点的N近邻，计算N近邻数据点与所求网格点的距离，在二维平面空间，N个近邻离散点（,）到预插值网格点（,）的距离为：

其中，Zi为N个近邻离散点上的观测值，Z（x0，y0）为预插值网格点（0,0）上的估算值，N为参与计算的N近邻样本数，Di为预插值点与第个近邻点间的距离，P是距离的幂，这里取2。

在网格点形变插值方法中，除了距离插值算法，还需要用到基于PS 树或网格树的范围搜索和N近邻搜索。这两种搜索算法是类似的，下边阐述了基于网格树的最近邻搜索算法。

（1）从网格树的根节点开始搜索，如果目标点的当前维度小于节点的对应坐标，则向树的左分支搜索，否则到右分支搜索；

（2）重复上述过程，直到搜索的点为叶节点，并存储此网格点为离目标点最近的网格点；

（3）然后在此叶节点向上搜索，并在每个节点重复以下计算过程：①如果当前节点离目标点距离比之前保存的最近距离更近，则将当前节点作为离目标点最近的网格点；②在另一分支搜索是否存在更好的点，首先取以当前最近距离为半径，以目标点为圆心的圆，然后判断是否与该节点对应区域相交，如果相交，那么向下移动到另一边搜索，否则向上到另一级节点。

（4）当向上搜索到根结点，则结束递归，将最后保存的具有离目标点最近距离的网格点作为最近邻点。

基于上述的最近邻算法，则能够实现基于目标点的范围搜索算法和N近邻搜索算法。

1.3 网格点三维形变计算

采用上述方法，能够计算多源卫星监测数据在地表网格点上LOS 方向的位移。如果将LOS 向、正东、正北和垂直向形变用、、、表示，则InSAR 观测的地表形变可以表示为：

式中，θ为InSAR 的LOS 向与垂直方向的夹角；α为卫星航向与正北方向沿顺时针方向的夹角。根据公式，显然用单一的维的LOS 向形变信息是难以求解、、三个方向的形变，但是可以借助多源卫星数据联立求解，则公式（8）可以转化为以下模型：

在公式（9）中，带入三组不同平台和轨道的InSAR 监测数据，即地表形变值，以及它们的飞行方向α及入射角θ，求解方程则能够得到网格点的、、。

因此，基于以上BPSTSIM 算法，能够基于多源卫星的PS 点数据，计算感兴趣区域任意一点的地表三维形变。

2 管道沿线多源数据融合

为了进一步验证BPSTSIM 算法的有效性，本文获取了管道沿线监测区域的多源卫星PS 点数据，包括TerraSAR、哨兵升轨、哨兵降轨的PS-InSAR 监测数据，实施监测区域的网格化和网格点形变插值计算。

图4 给出了部分区域的空间插值结果。图4（a）是PS 点实际位置和分布情况。图4（b）显示的是插值后网格点的分布。空间插值后监测数据能够完整覆盖研究区域，监测数据稀疏的地区数据更加密集。

图4 管道沿线PS 点的空间插值结果（局部）（a.PS 点分布；b.插值后网格点分布）

表1 为InSAR 数据在时域上的部分插值结果，经过时间插值后，原本时间间隔为数十天的InSAR 数据在时间维度上已经细化统一到以天为间隔。

表1 部分时域插值结果

有了多源卫星PS 点及网格点插值的地表形变，则能够基于公式（9）计算所有网格点在、、三个方向的形变值。图5 显示了某一时刻感兴趣的油气管道区域所有网格点的形变值、、，以及它们的矢量和。通过三个方向的形变可视化图，能够清晰地了解在不同方向形变的情况，以及合成以后总体形变的尺度。相比于单一LOS 向的形变信息，三维方向的形变更加立体丰富，也避免了由于InSAR对垂直于LOS 方向的形变不敏感造成的形变漏识别，合成后的总体形变体现了三维方向形变的累积。

图5 多源卫星数据融合结果

我们根据计算出的总体形变量，选取了两处形变量较大、具有一定坡度且距离管道较近的区域进行野外调查（图6），发现A 形变区域存在一处滑坡，B 形变区域存在一处不稳定斜坡（图6）。野外形变区域位置与InSAR 监测形变位置基本吻合，证实了融合方法的合理性和可靠性。

图6 野外调查照片（左：滑坡，右：不稳定斜坡）

3 结论

为了克服油气管道区域地表形变监测困难的问题，本文提出了基于多源卫星数据融合的油气管道区域三维变形监测方法BPSTSIM。它利用多源SAR 数据，获取高相干性的PS 点地表形变位移，然后通过对PS 点的时间插值获取更细时间粒度的PS 点形变值，通过对油气管道区域网格化以及基于PS 点的空间插值，获得多源PS 点数据对应网格点的形变数据，最后计算获得时间维度和空间维度监测油气管道区域的三维形变。实验数据以及可视化结果表明：

（1）利用BPSTSIM 方法，能够基于多源卫星的数据融合，在不同空间粒度和时间粒度监测油气管道的形变状态；

（2）基于BPSTSIM 计算结果，能够分析油气管道区域三维形变及总体形变态势；

（3）融合后获得的总体形变态势能够较好的与野外实际情况吻合，证明该方法合理且可靠。