山西省臨汾市礦區地表形變InSAR大范圍探測與監測

麻學飛， 張雙成， 惠文華， 許 強

(1.長安大學地質工程與測繪學院，西安 710054； 2.地理信息工程國家重點實驗室，西安 710054)

0 引言

煤炭作為我國的支柱能源，關系到國家經濟的命脈，需要不斷地進行開采，這勢必會形成采空區。采空區上部巖層發生斷裂、彎曲等形變，會造成采空區沉陷，破壞當地生態環境，威脅居民的人身財產安全[1-2]。因此對礦區進行大范圍探測和動態監測具有重要的研究意義。水準測量、全球定位系統(global positioning system，GPS)等方式是傳統的礦區地表形變監測手段，但由于這些方法監測成本高、工作量大等原因，難以進行長期監測，且所獲監測結果時空分辨率低，難以反映真實的礦區形變過程[3-5]。

合成孔徑雷達干涉測量(interferometry synthetic aperture Radar，InSAR)作為一種高精度的測量手段，利用相位差獲取研究區地表形變信息，具有監測范圍大、效率高、分辨率高的特點，幾乎不受云霧的影響，成為了礦區地表監測的有效手段[6-7]。 差分合成孔徑雷達干涉測量方法(differential interferometric synthetic aperture Rader,D-InSAR)可獲取到2景合成孔徑雷達(synthetic aperture Radar,SAR)影像覆蓋期間形變體表面的微小形變信息[7]，從而實現對形變體進行監測，處理數據量小，可以快速獲取目標區域的形變信息。Ge等[8]利用D-InSAR技術對Appin，Westcliff和Tower煤礦區地表形變進行監測，所獲監測結果精度可達厘米級； 龍四春等[9]利用D-InSAR與GPS集成實現了對資興唐煤礦區三維形變監測，并對測量結果進行集成與內插，驗證了2種監測結果的一致性。但是D-InSAR技術難以獲取目標區域時序變化過程，且在實際應用過程中常常受到時空去相干、大氣延遲等影響[10-11]。而時序InSAR彌補了D-InSAR技術的不足，在一定程度上解決了時空失相干、數字高程模型(digital elevation model,DEM)誤差、大氣延遲造成的影響，監測精度可達厘米級甚至毫米級[12-13]。永久散射體測量技術(permanent scatters InSAR，PS-InSAR)利用影像上存在的永久散射點來替代整幅影像上的像素點，通過永久散射點的形變來呈現目標區域的形變[14]，但是只適合小形變速率、小范圍的形變，無法很好地適用于沉陷區大范圍探測與監測； 小基線集(small baesline subset,SBAS)-InSAR能以更高的時空密度來提取到相干點目標，其形變結果更加明顯、直觀[15-16]，是目前礦區沉陷監測最主要的方法之一。沙永蓮等[17]使用SBAS-InSAR技術對新疆哈密砂墩子煤田礦區進行監測,發現了一個沉降漏斗，時序結果顯示沉降漏斗先線性下沉后逐漸趨于穩定； 李達等[18]通過對13景Terra SAR-X數據進行SBAS-InSAR技術處理,監測到了陜西榆林某煤礦的地表形變，并提取多個觀測點的時序沉降值進行量化分析，結果表明 SBAS-InSAR技術可為礦區地表形變監測和分析提供新的監測手段。但是時序InSAR技術對影像的數量有一定要求，處理數據量較大，周期長，無法高效率地對大范圍的形變區域進行監測。針對礦區地表形變大范圍探測與監測的問題，可以將D-InSAR技術和SBAS-InSAR技術結合起來使用： 首先利用D-InSAR技術對目標區域進行大范圍探測得到沉陷區位置，再利用時序InSAR技術對目標沉陷區進行時序監測和形變追蹤，高效率地完成對目標區域的形變監測。

本文獲取了覆蓋山西省臨汾市全域不同軌道的C波段Sentinel-1A數據，采用聯合D-InSAR技術和SBAS-InSAR技術開展臨汾市沉陷區的大范圍探測和監測工作，總計探測出沉陷區105處，選取其中的2處礦區進行詳細分析，并結合光學影像進行驗證。

1 研究方法

1.1 D-InSAR基本原理

D-InSAR技術通過采用目標區域形變前后兩景影像共軛相乘獲得干涉圖，與DEM模擬相位差分消除地形相位后獲取到地表形變信息。但是在D-InSAR用于具體的地表形變監測過程中，由于大氣環境和地表形變處于不斷的變化的過程，干涉相位會不可避免地受到影響，進而影響到最終的測量結果[19]。差分相位表達式為：

φ=φflat+φtopo+φdef+φatm+φnoi，

(1)

式中：φ為形變相位；φflat為平地相位，可借助精密軌道數據去除；φtopo為地形相位，可借助DEM數據模擬地形相位去除；φdef為地表形變引起的相位；φatm為大氣延遲相位;φnoi為隨機噪聲相位，可采用濾波的方式進行消除。

1.2 SBAS-InSAR原理

小基線技術利用同一軌道重復觀測獲取到目標區域的SAR影像集，通過設定時空基線將SAR影像集進行分組，形成多個小基線子集，再分別對每個子集內的SAR影像進行差分干涉生成差分干涉圖和相干性圖，利用平均相干性從相干性圖中選取高相干性的像素點，對差分干涉圖依次進行相位解纏得到相位序列，聯合多個小基線子集的相位序列來預估形變信息。若在一定時間內有N+1景 SAR 影像覆蓋目標區域，設定時空基線后可以組合成M幅干涉圖，M需要滿足[20]:

(2)

以t0時刻作為初始時刻，則tk時刻的差分干涉相位φtk為觀察量，其中1

(3)

式中：d(ta,x)和d(tb,x)為ta和tb時刻雷達視線方向上的形變值;φ(ta,x)，φ(tb,x)分別為d(ta,x),d(tb,x)對應得解纏相位值；λ為波長；Ti為主輔影像時間間隔；i為N幅影像獲取時間段內某一時刻；v(x)為像元x處的形變速率。N幅影像可以組成M個方程，其方程組的形式為：

Aφ=δφ，

(4)

式中：A為M×N的系數矩陣，每一行代表一副干涉影像，每一列代表對應時間上的SAR影像。若M>N，且A的秩為N，依據最小二乘法估值為：

(5)

式中ATA為奇異矩陣，存在無數解，可用聯合多個小基線采用奇異值分解法求出上式最小范數下的二乘解。

D-InSAR技術和SBAS-InSAR技術2種方法使用的場景不同，根據需求聯合2種監測技術實現對于臨汾地區礦區的大范圍探測和動態監測，技術流程如圖1所示。

圖1 技術流程圖

2 研究區概況及數據源

山西省臨汾市位于N35°23′～36°57′,E110°22′～112°34′之間，處于太原、鄭州、西安3個省會城市連線的中心點，地理區位優勢明顯，礦產資源豐富； 地處屬溫帶大陸性氣候，四季分明，雨熱同期[21-22]； 地表多為黃土，土質疏松，夏季多暴雨，地表植被容易遭到破壞。臨汾市包含河東煤田北部、霍西煤田中部和沁水煤田西部，煤炭資源豐富，2019年原煤產量6 201.6萬t，占全省原煤產量的6.39%。臨汾市礦區多分布于臨汾斷陷盆地山體之中，地表以植被為主，多年的煤炭開采造成地表沉陷嚴重，容易引發滑坡、地面塌陷等多種地質災害，迫切需要有效的沉陷區監測手段。

實驗選取432景不同軌道觀測獲取的Sentinel-1A影像(圖2)，其中path=11，frame=111有103景，path=11，frame=116有103景，時間跨度為2017年3月12日—2020年11月27日；path=113，frame=111有113景，path=113，frame=116有113景，時間跨度為2017年3月19日—2020年12月28日。影像均為C波段，周期為12 d，模式為升軌，VV極化，距離向分辨率為5 m，方位向分辨率為20 m。DEM數據為AW3D30 DEM，分辨率優于30 m。衛星軌道參數為歐空局提供的POD(AUX_POEORB，POD Precise Orbit Ephemerides)精確軌道參數。干涉雷達處理軟件為瑞士GAMMA軟件。

圖2 影像和DEM覆蓋范圍Fig.2 Image and DEM coverage

3 數據處理結果與分析

3.1 礦區大范圍探測

礦區形變具有連續性，受非人為因素的干擾較小，只要在礦區連續沉陷過程的某一時間段內的形變過程被監測到，則認為該區域為沉陷區。本文選取2017年末、2019年初和2019年末3個時間段共12景不同軌道的Sentinel-1A影像組成覆蓋臨汾市全域的6組干涉對共3個組合(組合a、b、c)進行大范圍探測。較短的時空基線可以在一定程度上減小時空失相干帶來的影響，為了保證干涉對具有較高的相干性，提高觀測精度，選取影像時時間基線設置為12 d，空間基線最大值控制在±100 m以內。

在數據處理過程中，將距離向和方位向的多視比設定為4∶1，地形相位利用外部DEM去除，干涉對中存在的噪聲采用自適應濾波的方法去除，相位使用奇異值分解法進行解纏，大氣相位應用GACOS技術去除，最終獲取到高質量的解纏圖。由于礦區在短時間內形變量級較大，形變相位在解纏圖中特征明顯，可直接依靠相位特征獲取到沉陷區位置，將沉陷區大致范圍用矩形框標識出來。由于沉陷區在地理位置上呈現集群式分布，為了后續的數據方便處理將識別結果劃分成6部分，依次編號為A—F，結果如圖3所示。

圖3 D-InSAR解纏圖識別結果

圖3(a)—(c)中黃色虛線為獲取到的沉降區的大致位置。圖3(d)中紅色矩形框圈畫出來的即為沉陷區，形變相位表現為由青-紫-青變化； 黑色的空洞是因失相干而產生的監測結果缺失。造成失相干的原因主要有2種： ①因沉陷區形變量過大，超出了C波段的監測能力而產生失相干； ②因植被茂密導致地物介電質常數發生變化而導致的失相干。沉陷區的也存在2種狀態： ①形變量級較大的沉陷區在解纏結果上表現為中心呈空洞狀態，周圍有形變相位環繞； ②形變量級較小的沉陷區表現為紫色或者黃色的形變相位環繞。沉陷區主要分布于臨汾斷陷盆地兩側的山體中，沿盆地東西兩側山體走向分布，呈點狀集群分布。分布于盆地中心的附近的沉陷區，地表植被較少，具有較高的相干性，獲得結果完整，易于識別沉陷區； 盆地兩側山體因植被影響造成大面積的失相干，給沉陷區識別帶來了一定的困難。C區域的沉陷區分布區域面積最大，數量最多，后文將著重以C區域作為重點進行表述。

3.2 礦區時序監測

基于D-InSAR識別到的沉陷區位置，依次對6個區域進行SBAS-InSAR時序處理，數據處理詳細參數如表1所示。2個軌道公共主影像時間分別為2019年1月1日和2019年1月8日，將主影像作為統一的基準進行配準并作為時間參考點，設定時間基線不少于90 d，空間基線不超過120 m。 為了提高結果的精度，使用強度互相關算法對主輔影像進行配準，通過設置不同的窗口大小、多項式參數個數等參數提高影像配準精度。通過外部DEM去除地形相位，利用時空濾波方法去除大氣誤差，選取高相干性點位作為基準進行相位解纏，最終獲取到6個區域垂向的形變速率，結果如圖4所示。從圖4中可以看出，識別所得開采沉陷區范圍均在D-InSAR監測結果中，其中灰色線條為手工勾繪的沉陷區邊界。本文此次探測出105處沉陷區，總面積達188.75 km2，其中7處與行政界線相交。識別所得沉陷區數量多，單個沉陷區沉陷面積較小，部分沉陷區形變量級較大，多處開采形變速率超過-100 mm/a，可以反映出臨汾市開采沉陷情況非常嚴重，這足以給當地礦區生態環境造成巨大的破壞。

表1 數據處理詳細參數Tab.1 Detailed data processing parameters

圖4 臨汾市2017年3月—2020年12月SBAS-InSAR沉陷區形變速率圖Fig.4 Deformation rate map of SBAS-InSARsubsidence area in Linfen City fromMarch 2017 to December 2020

選取形變量級最大的C區域進行詳細分析，獲取到的累計形變量和形變速率結果如圖5所示。C區域存在43處沉陷區，2017年3月19日—2020年12月28日該區域最大累計形變量為1 030 mm，最大形變速率為-381 mm/a。區域右側植被覆蓋較少，影像相干性較高，可以獲取到完整沉陷區形變，而左側區域為山區，植被覆蓋率高，白色區域即為影像因植被影響產生失相干所造成的結果缺失。其中C1沉陷區由2個相鄰的礦區構成，左側沉陷區形變量級大于右側，最大累計形變量分別為922 mm和593 mm。C2沉陷區為單沉陷區，形變輪廓呈現為條帶狀，沉陷區的輪廓反映出礦區由西南向東北方向不斷開采，最大累計形變量為942 mm。

從圖6可以看出，C1，C2兩個沉陷區均處于不斷開采過程中。從形變過程來看，4個點位在礦區開采過程中形變量級在不斷增大，形變過程呈現快速沉陷-緩慢沉陷-快速沉陷的周期，以p1點的時序變化表現的最為明顯。在2017年5—8月，由于沒有獲取到有效的干涉對，故時序呈現擬合的直線。p1點整體形變速率最大，形變量總是經歷較長時間快速形變到較短時間緩慢形變的周期循環，形變速率的絕對值也是呈大-小-大的周期變化，但是整體的形變趨勢逐漸緩和，最大沉降量為922 mm。p2點整體形變趨勢較為平穩，周期性形變特征不太明顯，最大沉降量為593 mm。p3點位于C2礦區沉陷中心，p3點的整體形變速率要大于p4點，呈現快速下降-穩定-略微抬升-緩慢下降-快速下降的過程，最大沉降量為942 mm。p4點位于礦區開采沉陷方向的中間位置，與p3點為同一礦區，整體形變趨勢與p3相似，最大沉降量為404 mm。

將C1，C2礦區分別沿圖6中ab，cd剖線截取剖面，結果如圖7所示。C1礦區有2個沉降漏斗，左側的沉降量要大于右側，2個沉陷區中間的位置已經有100 mm的沉降量，已經成為一個雙中心的大型沉陷區。C2礦區沉陷左側形變傾角大于右側，說明礦區開采的方向是沿剖線cd方向進行開采，沉降漏斗中心位于起采區，但是隨著開采工作的進行，沉降漏斗中心會逐漸向開采的方向移動，沉陷區輪廓會呈現為條帶狀。

由于沒有礦區的實測數據，本文利用Google Earth遙感影像數據檢驗監測結果的一致性。首先采用閾值分割法將C1，C2沉陷區邊界提取出來，經過多次嘗試當閾值設置為-10 mm/a時，所獲得沉陷區邊界效果最好，將礦區邊界疊加至光學影像上如圖8所示。結合光學影像進行分析，2個礦區均為地下開采，沉陷范圍內有建筑物存在，地表分布有大量農田，沉陷面積分別為1.54 km2，2.87 km2。C1礦區南側地勢平坦，已被開墾為農田，北側為山地，人工地物較少，中間由一條道路橫穿沉陷區，左側100 m附近有一村落，周圍發現2處疑似采煤作業區。C2礦區地表起伏較大，但地表已被全部開墾為農田，中間由多條道路穿過，東南側與小鎮相接，東北方向50 m有一條國道，北側分布有3個疑似采煤作業區。C1,C2礦區沉陷邊界直接與農田、道路相接，距離住宅不足100 m，產生的地面沉陷會直接威脅到周圍居民的人身財產安全，尤其C2礦區，需要持續進行監測。

4 結語

煤礦開采所造成的沉陷區對地表生態環境和基礎設施帶來了巨大的破壞。本文利用D-InSAR技術對臨汾市沉陷區進行大范圍探測和監測結果工作，將探測結果分成了6個區域并依次使用SBAS-InSAR技術進行時序監測，共發現105處沉陷區； 將其中2處典型沉陷區作為重點區域進一步分析，獲取到了該沉陷區的形變演化情況，為后續礦區沉陷災害防治提供重要依據。通過光學影像發現了沉陷區附近存在疑似采煤作業區，驗證了InSAR技術大范圍探測與監測具有可行性和可靠性。

1)礦區地表形變具有形變速率快、量級大的特點，會對當地的生態環境造成巨大的破壞，難以通過傳統的光學遙感進行探測，InSAR技術可以極易捕獲沉陷區的大量級地表形變，從而實現對研究區沉陷區的大范圍探測。

2)直接利用時序InSAR技術對大范圍區域進行沉陷區探測具有巨大的挑戰性，處理數據量大，工作效率低，結果不穩定。而基于D-InSAR技術的大范圍探測和SBAS-InSAR技術的詳細監測的方法可以很快實現對研究區進行層次化的處理，真正做到對大范圍區域的高效率監測。

3)本文基于InSAR技術使用Sentinel-1A升軌數據對臨汾市礦區地表形變進行了大范圍探測和監測，識別和監測結果對臨汾市礦區形變監測具有較好的參考價值，但是InSAR技術在植被覆蓋率較高的區域適用性不佳，未來還需要借助其他技術來提高大范圍探測和監測結果效率。

志謝：感謝歐洲太空局在科學數據中心網站公布Sentinel-1A衛星雷達影像，感謝日本宇宙航空航天局地球觀測研究中心提供DEM數據。