基于多源衛星數據融合的油氣管道區域三維變形監測

余東亮，周廣，方迎潮，王明波，薛廉，趙雪

（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）顯示的是插值后網格點的分布?？臻g插值后監測數據能夠完整覆蓋研究區域，監測數據稀疏的地區數據更加密集。

圖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）融合后獲得的總體形變態勢能夠較好的與野外實際情況吻合，證明該方法合理且可靠。