摘要：
首先，介紹了InSAR地表形變監測技術的基本原理和技術特點；然后，對當前InSAR礦區地表三維形變監測方法做了系統性分類，并探討了各自的優缺點和適用范圍；之后，綜述了當前InSAR礦區地表三維形變預計方法的研究現狀和進展；最后，歸納梳理了InSAR礦區三維形變監測預計在多源數據融合和沉降機理分析等方面的前沿問題。
關鍵詞： InSAR；開采沉陷；地質災害；綜述
地下礦產資源開采容易打破上覆巖層的原始應力平衡狀態，致使巖層和地表發生移動和變形，從而引發一系列礦山地質環境問題和災害（如含水層破壞、地表塌陷、山體滑坡、建構筑物損毀）。因此，監測和預計礦山地表三維形變對于評估潛在地質災害，分析礦區沉降機理具有至關重要的作用。傳統的礦山地表形變監測通常沿著地下工作面走向和傾向布設兩條觀測線，然后使用水準儀和GPS等傳統大地測量手段定期或不定期觀測兩條剖面線的下沉和水平移動。該方法雖然測量精度較高，但監測范圍小，工作量大，成本高，效率低，測點易破壞，且不能連續監測。因此，傳統方法通常僅監測礦區少數工作面開采導致的地表“線狀”形變，且時空分辨率較低。事實上，傳統“線狀”形變觀測很難完整地反映整個礦區地表三維形變的時空演化規律，從而阻礙了礦山開采沉陷的全面監測及其相關聯的理論研究。此外，僅基于礦山低時空分辨率的“線狀”形變觀測也難以可靠地反演礦區開采沉陷模型參數，從而削弱三維形變的預計精度和可靠性，進而也降低了對潛在地質災害預測的可靠性。
合成孔徑雷達干涉測量是一種新型主動式地表變形監測技術。該技術發展之初主要被用于提取地表高程信息。1989年，首次提出差分InSAR（Differential InSAR，DInSAR）的概念，并通過去除InSAR干涉圖中的地形相位貢獻，實現了地表形變和地形的分離，從此開啟了InSAR地表形變監測的新篇章。相對于傳統大地測量手段，InSAR技術具有全天候、全天時、成本低、覆蓋范圍大、空間分辨率高等優勢。以“哨兵一號”（Sentinel 1）衛星TOPS干涉模式SAR數據(對所有用戶開放)為例，其空間覆蓋范圍約為250 km×160 km，地距向和方位向的空間分辨率約為5 m×20 m，最短重返周期為6天（以Sentinel 1A和1B星座協同）。理想情況下，利用TOPS模式獲取的SAR數據可一次性免費監測約40 000 km2，最短時間間隔為6天，空間分辨率約為20 m的地表形變圖（包括該范圍內的所有礦區）。如此高時空分辨率的地表形變圖是利用傳統大地測量手段幾乎無法獲得的。因此，InSAR技術為礦區三維形變監測和預計提供了全新的手段。
相對于地球物理領域，InSAR技術在礦區地表形變監測中的應用相對較晚。1996年，才首次嘗試利用DInSAR監測地下采礦導致的地表形變。在此之后，如何利用DInSAR或時序InSAR（Multi-Temporal InSAR, MT-InSAR）技術準確獲取礦區地表形變成為了新的研究熱點。然而，由于SAR傳感器斜視成像的特點，InSAR監測的地表形變實際僅為一維視線（line-of-sight，LOS）方向形變，而非地表在垂直、東西和南北方向真實三維形變，這一缺陷影響了InSAR技術在礦區的實際應用(如地質災害風險評估和沉降機理解譯等)。
為了克服以上局限，近年來國內外學者開展了大量基于InSAR一維形變觀測值的礦區地表三維形變重構和預計研究，發展了眾多行之有效的方法。已經對InSAR在礦區地表形變監測和應用進行了較全面描述，但其并未對當前InSAR礦區地表三維形變和預計理論與方法進行詳細地歸納和總結。因此，為了較為系統地梳理該領域的理論進展，本文總結了當前InSAR礦區地表三維形變監測和預計的研究現狀，簡單分析了現有方法的使用范圍，并歸納了存在的問題以及未來可能的發展方向。
InSAR礦區形變監測技術方法概述
1InSAR技術基本原理概述
InSAR可基于兩景覆蓋同一地區且雷達成像幾何具有微小差異的SAR影像的干涉信號，提取地表高程信息或分離地表形變。以目前最為常用的重復軌道InSAR技術為例，兩景SAR影像的干涉相位?int可近似表示為
(1)
式中，?topo和?defo分別為地表地形起伏和地表變形引起的相位；?flat、?orbit、?atm和?noise分別為平地效應、軌道誤差、大氣延遲及噪聲導致的相位，這4項相位貢獻可通過理論公式計算、模型擬合、濾波等方式削弱或消除[16]。至此，式(1)可簡化為
(2)
式中，?′int=?int－?flat－?orbit－?atm－?noise。
1.1 DInSAR技術
從式(2)中可以看出，若兩景SAR影像獲取期間地表變形不明顯（即?defo≈0），則式(2)中的干涉相位?′int可認為主要由地形相位?topo貢獻。這就是利用InSAR技術獲取地表數字高程模型的核心思想。相反，若能估計并去除干涉相位?′int中地形相位貢獻?topo，則可實現地表形變相位?defo的分離。該技術即為DInSAR的核心思想。
DInSAR技術是基于InSAR干涉相位獲取地表變形信息，因此，干涉相位質量直接決定了其獲取地表形變的可靠性。目前，該技術已被廣泛應用于礦區地表形變監測。然而，礦區地表覆蓋變化較快、形變速率較大、沉降范圍較小等因素極易削弱干涉相位質量，甚至導致相位失相關（即干涉相位完全不可靠），從而嚴重降低DInSAR監測的礦區地表形變觀測值精度。因此，在利用DInSAR技術監測礦區地表形變時，應盡量選擇波長較長（如L波段）、分辨率較高的SAR數據，從而盡可能地減少相位失相關對礦區形變監測結果精度的影響。
1.2 MT-InSAR技術
DInSAR技術易受時空失相關、大氣延遲等因素影響，且僅能獲得兩景SAR影像期間的地表差分變形。為了克服該局限，MT-InSAR技術應運而生。該技術根據主影像數量不同，可大致分為兩大類：
(1) 單主影像MT-InSAR技術，以永久散射體（persistent scatterer，PS）InSAR為代表。該技術選取覆蓋同一地區的單軌SAR數據集中的一景為主影像，并將其與其他影像（從影像）差分干涉；然后，根據SAR影像幅度和相位信息選出散射特性相對穩定，且回波信號較強的像素（即PS點）；最后，根據差分干涉相位重構出PS點的時序形變。
(2) 多主影像MT-InSAR技術，以小基線集（small baseline subset, SBAS）InSAR為代表。SBAS-InSAR基于覆蓋同一地區單軌SAR數據集中時空基線均小于一定閾值的DInSAR干涉對（有多個主影像），并根據其干涉相位，重構SAR影像范圍內高相干點的時間序列。
目前，MT-InSAR（包括PS-InSAR和SBAS-InSAR）已被廣泛應用于礦區地表時序形變監測。然而，這兩種MT-InSAR方法在礦區的適應性存在一定的差異。我國大部分礦區位于郊區或農村地區，地表通常被大量的植被或農田覆蓋，因此潛在PS點較少。加之PS-InSAR基于單一主影像，增大了礦區形變梯度和時間失相干，因此，PS-InSAR技術主要用于地表變形較慢的礦區（如每年形變從幾毫米到十幾厘米的區域）。相反，SBAS-InSAR技術由于采用時空基線較短的多主影像模式，因此，相對于PS-InSAR，該技術能夠應用于形變較快的礦區（如幾厘米甚至幾十厘米每年）。
1.3 MAI技術
與DInSAR僅能監測LOS向形變不同，多孔徑雷達干涉（Multiple Aperture Interferometry，MAI）技術通過方位向公共頻譜濾波將兩景SAR影像分別拆分成前視與后視影像。之后，基于拆分的前后視影像分別生成前視和后視干涉圖。差分前后視干涉圖可獲得地表方位向形變。MAI技術目前也被成功引入礦區地表方位向形變監測。2015年，利用MAI技術獲取了撫順西某露天礦的方位向變形。2018年，利用MAI技術估計了山東濟寧東灘煤礦地表方位向形變。然而，由于MAI技術對干涉相位噪聲非常敏感，而礦區地表相干性通常不高，所以該技術目前在礦區（特別是煤礦礦區）地表方位向形變監測中的成功案例不多。
1.4 Offset tracking技術
像素偏移量追蹤（offset tracking，OT）是一種能從兩景SAR影像中同時獲取地表LOS和方位向二維形變的方法。該方法基于兩景SAR影像的強度或干涉條紋信息，利用互相關最優化或相干性追蹤等方式獲取SAR影像在LOS和方位向的像素偏移量。然后，在去除軌道和地形偏移量之后，即可獲得地表沿著LOS和方位向的二維形變。相對于DInSAR和MAI技術，OT可以僅利用SAR影像的幅度信息估計地表二維形變，且不需要相位解纏。鑒于其在地表大變形監測方面的優勢，2013年，首次將OT技術引入礦區地表大形變獲取，并利用ALOS PALSAR數據成功獲得了內蒙古地區兩個煤田開采導致的地表大變形（最大值約為4.5 m）。之后，OT監測礦區地表大形變的潛能在眾多礦區得到證明。
OT技術獲取的形變精度很大程度上取決于所用SAR影像的空間分辨率（約為1/10~1/30個像素）以及影像的紋理特征。因此，在利用該技術監測礦區地表大量級二維形變時應盡量選擇波長較短且分辨率較高的SAR數據（如TerraSAR-X SpotLight模式影像）。需要指出的是，理論上，基于幅度信息的OT技術并未進行相位干涉，因此本不應屬于InSAR技術。但為了便于描述，本文將DInSAR、MT-InSAR、MAI、OT等統稱為InSAR技術。
2InSAR形變觀測值局限性分析
利用DInSAR技術僅能獲取礦區地表一維LOS向形變。即便聯合MAI或使用OT技術也僅能獲得地表沿著LOS和方位向的二維形變。事實上，如圖 1所示，InSAR監測的LOS dLOS和方位向形變值dazi是地表在垂直、東西和南北方向的形變分量W、E和N的投影，即
圖 1 InSAR監測的LOS向形變與地表三維形變的投影
(3)
式中，θ和αh分別為SAR傳感器的入射角和飛行方位角。
由于地表形變發生在三維空間（即垂直、東西和南北三個方向），因此，僅利用InSAR形變觀測值很難反映地表真實形變模式。以LOS方向為例，在極端情況下，僅當地表形變完全沿著SAR衛星的LOS方向移動時，InSAR監測的LOS形變才能完全反映地表真實變形；反之，當地表變形垂直于LOS方向時，InSAR則完全無法監測地表形變。此外，由于礦區地表形變是一個復雜的時空演化過程，因此，僅基于InSAR監測的形變觀測值也難以可靠預計后續開采導致的地表三維形變。鑒于此，研究InSAR礦區地表三維形變監測和預計方法十分必要。
InSAR礦區地表三維形變監測方法
由單軌InSAR一維形變觀測值估計實際的三維形變本質是秩虧問題。因此，增加額外觀測方程是目前InSAR礦區地表三維形變估計最主要思路。根據額外觀測方程構建基礎不同，本文將現有InSAR礦區地表三維形變監測方法分為兩大類，即基于多軌InSAR觀測法和先驗模型+單軌InSAR觀測法。
1基于多軌道InSAR觀測值方法
該方法基于雷達成像幾何具有顯著差異的多軌InSAR干涉對，利用DInSAR、MAI和OT等InSAR技術獲得多個方向的一維或二維形變觀測值，從而增加額外觀測方程，補充三維形變觀測信息。因使用的形變獲取方式不同，該類方法可分為多軌DInSAR觀測法和多軌DInSAR+MAI/OT觀測法兩類。
1.1 多軌DInSAR觀測法
該方法利用DInSAR技術，從3個或以上雷達成像不同平臺/軌道的InSAR干涉對中獲取礦區地表LOS向形變[dLOS1 dLOS2 … dLOSk]，其中k表示InSAR干涉對的個數。然后，將DInSAR監測的多軌LOS形變內插到統一格網，并基于式(3)構建礦區地表三維形變與多軌LOS形變觀測值之間的觀測方程組
(4)
最后，基于式(4)并利用最小二乘或加權最小二乘等方法求解礦區地表三維形變。
基于該思路，2011年，利用升降軌ENVISAT ASAR和ALOS PALSAR-1數據獲取了澳大利亞Appin和Tahmoor煤田開采導致的地表三維形變。2014年，利用同樣的方法獲取了我國某礦區地表三維形變，并與水準測量結果進行了對比。2013年，借鑒SBAS-InSAR思想，將該方法拓展到融合多軌InSAR觀測值的礦區地表三維時序形變監測，并將其應用于德法邊界礦區殘余二維（垂直和東西方向）時序形變監測。2018年，利用概率積分法預計地下開采導致的南北方向的水平移動，并將預計的結果與TerraSAR-X和Radarsat-2數據結合，實現了河北峰峰礦區地表三維形變估計。
1.2 多軌DInSAR+MAI/OT觀測法
2015年，文獻利用DInSAR和MAI技術分別處理了單軌InSAR數據集生成礦區地表沿著LOS和方位向的二維形變觀測值，并利用SBAS-InSAR思想實現了我國撫順西露天礦開采導致的垂直和南北方向二維時序形變。理論上，基于多軌InSAR觀測值并利用的方法可獲取礦區地表高相干地區三維時序形變。2018年，基于升降軌SAR數據，結合DInSAR獲取的LOS向形變以及MAI或OT獲取的方位向形變估計了我國兗州東灘礦地表三維形變。
2先驗模型+單軌InSAR觀測法
與多軌InSAR觀測法的思路不同，基于先驗模型+單軌InSAR觀測法，通過礦區地表三維形變之間的先驗模型構建額外約束，從而實現基于單軌InSAR觀測估計礦區地表三維或三維時序形變。
2.1 先驗模型+單軌DInSAR觀測法
已有研究表明：在水平或近水平煤層開采條件下，礦區地表水平移動與下沉梯度之間呈線性關系。基于該先驗模型，2015年，構建了礦區地表東西形變分量E和南北形變分量N以及下沉梯度之間的理論約束方程
(5)
式中，(i, j)表示形變圖的像素坐標；BE和BN為東西和南北方向下沉梯度ΔWE和ΔWN的比例系數，其具體形式見文獻。融合兩個理論約束方程與InSAR監測的一維LOS形變值，即可實現礦區地表三維形變監測。圖 2為利用該方法獲取的我國安徽錢營孜礦區地表三維形變場。2018年，基于SBAS-InSAR和加權最小二乘平差將該方法擴展到了基于單軌InSAR數據集的礦區地表三維時序形變監測。
圖 2 基于先驗模型+單軌InSAR觀測值方法獲取的錢營孜礦區地表在東西、南北和垂直方向地表三維形變場
2.2 先驗模型+單軌OT觀測法
為了克服文獻提出的方法無法獲取礦區地表大量級（如幾米甚至幾十米）三維形變的局限，2017年，文獻提出基于單個InSAR干涉對的強度信息，并利用OT技術獲取礦區地表LOS向大形變場信息，通過融合構建的兩個理論約束（即式(5)），實現了基于單個InSAR干涉對強度信息的礦區地表大量級三維形變估計。2018年，又提出融合OT獲取的二維形變信息（即LOS和方位向）與式(5)的先驗模型信息，并利用加權最小二乘法提高三維大形變估計的穩健性。同年，利用SBAS-InSAR和穩健估計思想將該方法擴展到了基于單軌SAR強度影像集的礦區地表三維時序大形變監測。
3兩類方法優缺點和使用范圍分析
與基于先驗模型+單軌InSAR觀測法不同，基于多軌InSAR觀測法不需要先驗模型輔助即可直接求解礦區地表三維形變場。理論上，該方法比基于先驗模型+單軌InSAR觀測法應用范圍廣。但在實際應用中，該方法存在許多制約和不足，例如：
①該方法至少需要兩個或以上具有顯著幾何差異的InSAR干涉對，不僅增加監測成本，而且在實際應用中數據難以滿足；
②該方法假設多軌InSAR觀測值獲取期間礦區地表三維形變W，E和N是不變或者線性變化的，對于沉降較快的礦區難以滿足，從而嚴重削弱監測的三維形變精度；
③由于當前SAR衛星均為近南北飛行，導致該方法估計的南北方向形變精度非常差。
與基于多軌InSAR觀測法相比，基于先驗模型+單軌InSAR觀測法存在以下優勢：①減少了對InSAR數據的嚴苛要求，節約了礦區三維形變監測成本(節約一半以上)；
②僅用一個InSAR干涉對，因此不需要假定礦區地表變形為線性，提高了監測精度；
③引入了先驗模型約束，一定程度上能改善南北方向水平移動估計精度。
然而，該方法同樣存在以下局限：
①由于需要利用先驗模型構建額外約束，因此，獲取的礦區地表三維形變精度除了受LOS和方位向形變精度影響之外，還受到先驗模型可靠性的影響，在先驗模型不吻合的地方(如出現較大的臺階或裂縫)，該方法獲取的地表三維形變精度不高，特別是水平方向；
②使用的先驗模型是基于礦區地下水平或近水平煤層開采導致的三維形變推導而來，因此，對于急傾斜煤層開采等導致的地表三維形變區域則無法使用；
③無法應用于露天礦開采導致的地表三維形變。
鑒于以上分析，當礦區地表形變較為緩慢，且獲得了至少3個具有顯著幾何差異的InSAR干涉對時，推薦使用基于多軌InSAR觀測法獲取礦區地表三維形變。當InSAR干涉對無法滿足時，推薦使用基于先驗模型+單軌InSAR觀測法（特指煤礦開采）。
InSAR礦區地表三維形變預計方法
在InSAR礦區地表三維形變預計方法方面，目前研究成果并不豐富。根據文獻檢索發現，2015年，提出了一種基于InSAR和概率積分法（probability integral method，PIM）的礦區地表三維形變預計方法。該方法首先忽略DInSAR監測的LOS形變中的水平移動貢獻，并將LOS形變直接轉換到垂直方向。然后，基于轉換后的垂直形變分量反演概率積分法與下沉有關的模型參數，通過對無法直接估計的模型參數(如水平移動系數)賦予經驗值，從而實現后續開采導致的地表三維形變預計。該方法很大程度上推動了InSAR礦區地表形變預計的研究進展，但仍存在兩個局限：
①地表真實的水平移動被不恰當地轉換到了垂直方向，從而增加了垂直形變觀測值的誤差；
②概率積分法中與水平移動有關的部分參數（如水平移動系數）無法直接估計。
為了克服以上兩個局限，2016年，首先利用概率積分法表示礦區地表沿著垂直、東西和南北方向的三維形變分量，并利用InSAR的空間投影關系，建立了概率積分法全部模型參數與InSAR監測的LOS向形變之間的函數模型；之后，基于大量LOS形變觀測值并利用附加粗差剔除的遺傳算法，估計了概率積分模型的全部參數；最后，利用估計的模型參數預計了礦區后續開采導致的地表三維形變。該方法簡稱InSAR-PIM方法。圖 3為利用此方法預計的我國安徽淮北某礦區的地表三維形變。與的方法略有不同，2016年，提出首先利用概率積分法表示礦區地表垂直方向和衛星地距方向的變形，并基于衛星雷達成像幾何建立了概率積分模型與InSARLOS向形變之間的關系，之后反演模型參數并預計地表三維變形。2018年，基于構建的概率積分模型參數與LOS形變觀測值之間的函數關系，并使用模擬退火法反演模型參數近而預計礦區地表三維變形。
圖 3 InSAR-PIM預計的錢營孜礦區3212工作面全采(圖中藍色矩形)導致的地表沿著垂直、工作面走向和工作面傾向的三維形變
針對概率積分法模型在非充分開采條件下容易過高地預計地表三維形變的問題，2017年，提出利用簡化Boltzmann函數修正概率積分法，并將修正后的模型與InSAR監測的LOS向形變值結合，從而發展了一套適合于不同采動程度（非充分、充分和充分開采）下的礦區地表三維形變預計方法。2018年，將Knothe時間函數修正的概率積分法模型與InSAR監測的LOS向形變值融合，從而發展了一種基于InSAR的礦區地表動態三維形變預計方法和礦區建構筑物動態破壞風險評估體系。然而，這些方法基本都基于單軌InSAR數據獲取的礦區地表形變場，一定程度上也限制了反演的模型參數精度。此外，由于InSAR觀測值數據量大（幾百萬甚至上千萬），因此，模型參數反演耗時較長。
InSAR礦區三維形變監測預計挑戰和展望
1當前的主要挑戰
相位失相干是地表形變InSAR監測的首要挑戰。主要是由于雷達熱噪聲、配準誤差、多普勒質心漂移、垂直基線、地表環境變化、地表形變（或變形梯度）過大等因素造成的。相位相干性是評價InSAR干涉相位可靠性的重要指標，一旦出現失相干，相應區域的InSAR形變監測值是不可靠的。如前所述，我國大部分礦區分布于郊區或者農村地區，地表常常有植被或農作物覆蓋，因此，InSAR相位失相干在我國的大部分礦區比較常見，特別是對于波長較短的SAR數據，如ENVISAT ASAR和TerraSAR-X。此外，礦山地表沉降范圍較小，量級較大，從而加劇了InSAR失相關。在礦區地表形變無法準確獲取的前提下，已有的基于InSAR的礦區地表三維形變預計方法均可能失效。因此，相位失相干是目前InSAR礦區地表三維形變監測和預計的主要挑戰。未來隨著SAR衛星的重返周期的縮短，該現象或許能夠得到一定的改善。
SAR影像分辨率和地表復雜形變是InSAR形變監測的另一大挑戰。對于淺埋深或厚煤層開采礦區，其地表變形量級通常較大（如幾米甚至幾十米）。對于缺乏高分辨率SAR影像的礦區，其地表大量級形變監測精度也將受到極大的限制。
2未來發展
在礦區多源SAR數據融合方面的可能發展。對于基于多軌InSAR觀測估計方法而言，如何在顧及礦區地表時序形變高度非線性的前提下，融合越來越多的衛星數據以提高三維形變精度（特別是南北方向）是一個值得研究的課題。而對于基于先驗模型+單軌InSAR觀測法而言，如何融合多源SAR數據提高三維形變的穩健性、如何將該方法拓展到不同傾角煤層開采（當前僅適用于水平或緩傾斜）導致的地表三維形變監測等均是值得研究的問題。對于InSAR礦區地表三維形變預計而言，當前的方法大都僅基于單軌InSAR觀測值。如何利用不同雷達成像幾何的多源SAR數據提高礦區開采沉陷模型參數反演精度，如何減少模型參數反演的時間消耗等問題均有待解決。
在InSAR與多源異質數據融合方面的可能發展。融合InSAR數據與多源異質數據（如LiDAR、三維激光掃描儀、GNSS、水準等），可以實現多源數據之間的優勢互補。因此，如何對礦區多源觀測數據定權，如何設計合理可行的融合方案，從而實現礦區地表多量級三維形變監測和預計仍有待研究。
InSAR面域形變信息挖掘也存在較大的發展空間。由于InSAR礦區研究主要集中在地表一維或三維形變監測與預計，因此，如何深入挖掘InSAR監測的海量高精度、高時空分辨率的“面域”形變資料在礦山開采沉陷學理論研究及建模、礦區沉降機理解譯、地質災害風險評估等方面的潛在應用仍有待研究。
結論
InSAR技術為解決傳統礦區形變監測和預計的局限帶來了全新的契機。近年來，理論與方法均取得重要進展，尤其是基于多軌InSAR觀測法和先驗模型+單軌InSAR觀測法等眾多切實可行的InSAR三維形變監測方法。與此同時，在InSAR礦區地表三維形變預計方面也提出了多種融合InSAR與概率積分法或修正概率積分法的礦區地表三維或三維時序形變預計方法。這些研究豐富了礦山地表三維形變監測手段，拓展了InSAR技術的應用前景。然而，由于InSAR的固有局限（如相位相關干），工程化和市場化的InSAR礦區地表三維形變監測和預計仍受到一定的限制。未來隨著可用SAR衛星數量的增多及重返周期的縮短，多源SAR數據融合，不依賴于外部形變資料的三維形變精度評定等問題仍有待研究。
參考文獻：略
來源：測繪學報，2019，48(2)：135-144
