多時相SAR振幅信息的大量級沉陷區識別研究

褚夫玉，胡海峰，廉旭剛，朱利超

(1.太原理工大學礦業工程學院，山西 太原 030024;2.山東省水利勘測設計院，山東 濟南 250000)

煤層開采所引起的地表沉陷是采空區周圍應力場平衡-破壞-再平衡的力學現象，是我國主要的地質災害之一，長期影響著人民生命財產安全和生態環境[1-2]。山西省作為煤礦采出大省，地質災害頻發，開采沉陷區分布廣泛，農田、建筑物損毀、植被退化、地質環境惡化問題時有發生，礦區可持續發展受到嚴重威脅[3]。因此，利用先進的監測手段，研究開采沉陷區的運動特征及規律，識別并監測沉陷區破壞程度，對促進開采沉陷區損害的早期預警、礦區地質災害治理和可持續發展有重要意義。

開采沉陷是煤礦開采工程中常見的地質災害類型。傳統監測礦區地表形變的方法有水準測量、全球定位系統(global positioning system，GPS)等，測點精度高，但在地表溝壑縱橫、地形起伏較大的地區點位布設困難且不能維持監測點的完整性，難以有效開展大范圍、長時間的監測活動[4]。合成孔徑雷達干涉測量(interferometric synthetic aperture radar，InSAR)是監測和評估地質災害的新技術，該技術無需地面控制點，可進行全天候、高精度、大范圍監測，目前已應用于我國煤礦開采沉陷的監測研究中[5]。CARNEC等[6]在1996年首次將差分干涉測量技術(differential interferometry synthetic aperture radar，DInSAR)應用到法國Gardanne煤礦開采沉陷研究中，獲取了單個重訪周期內的最大下沉量，但受大氣延遲、地形誤差和時空失相干的影響，監測結果的有效性和可靠性較低；FERRETTI等[7]在1999年針對常規DInSAR受失相干和大氣影響的缺陷，率先提出永久散射體干涉測量技術(persistent scatterer interferometric synthetic aperture radar，PS-InSAR)，并利用ERS-1/2數據對意大利Camaiore鎮地表沉降進行監測，獲得了較為可靠的形變信息；WALTER等[8]在2004年利用PS技術對德國Prosper-Haniel煤礦進行開采沉陷監測，證明了該技術在礦區沉陷監測的可行性；2005年，CASU等[9]利用短基線集方法(small baseline subset，SBAS)進行地表沉降監測，獲得了與水準及GPS測量一致的結果，但忽略了InSAR形變數據與實測離散點時空尺度上的差異。

目前，國內學者先后利用DInSAR、PS-InSAR、SBAS等多種InSAR技術在我國煤礦開采沉陷區開展了一系列研究[10-16]。其中，PS-InSAR數據量要求較大(一般為20景以上)、PS點的分布呈隨機性且受時空失相干的影響較大，在植被覆蓋密集、地形變化復雜的礦區，難以找到足夠的相位穩定點，導致觀測結果不能客觀反映沉陷盆地整體變化[17]；SBAS技術通過選取短時空基線的干涉對，減小了時空失相干的影響，提高了形變監測點的密度，但在地形變化較大的區域目標的信噪比較低，進而影響形變監測結果的可靠性[18]。DInSAR、PS-InSAR、SBAS等是基于電磁波的相位信息提取微小形變量的技術，在形變劇烈的礦區，相鄰兩像元間的形變量大于雷達波長的1/4時，會產生相位混疊，難以解纏到大量級形變?；赟AR影像振幅信息偏移量追蹤技術的迅速發展，GRAY等[19]在1998年首次采用SAR影像的振幅信息解譯得到研究區地表形變；ZHAO等[20]在2013年采用該技術處理PALSAR影像獲取了內蒙古某煤田4.5 m的垂直位移。 利用SAR影像振幅信息通過多級配準獲取地表形變信息，不需要解纏，不受礦區相干性制約，能夠獲取大量級形變，但對配準精度要求極高，獲取的沉陷盆地邊緣形變精度較低。歐空局(ESA)發射的Sentinel-1A/B SAR系列衛星，具有開放數據源、短重訪周期、超幅寬、多極化等特點[21]，大大縮短了時間基線，在礦區開采沉陷準動態監測中具有很好的應用前景。

本文基于SAR影像振幅信息的動態時序開采沉陷監測方法，利用時序振幅數據提取開采沉陷區地表后向散射系數，分析開采沉陷過程中地表不同區域的后向散射系數與沉降速率的關系，較好地克服了常規SAR技術難以識別大量級沉陷區域的問題，可以提高開采沉陷監測的時空分辨率。以陽煤新元礦某工作面為研究區，利用該方法和16景Sentinel-1A數據，獲取研究區振幅時間序列圖和地表后向散射系數曲線圖，并使用地表后向散射數據分析了沉陷區不同區域時序變化幅度，對地表后向散射系數與沉降速率關系進行研究，發現二者具有極強的相關性。

1 研究區與數據

1.1 研究區概況

陽煤新元礦位于沁水煤田西北部，太原至陽泉的中間地帶，壽陽-陽泉構造堆積盆地地區，海拔1 058.3～1 110.9 m，屬黃土高原溫帶季風氣候區，天氣寒冷干燥，溫差變化大，年均降水量約500 mm，多集中在6～9月。 研究區地表大部分被第三系、第四系松散地層覆蓋，經長期風化剝蝕，溝谷縱橫，地形復雜。礦區南北長4～9 km，東西長15～16 km，地處東經112°58′51″～113°09′33″，北緯37°49′54″～37°52′09″。 工作面傾向長240 m，走向長1 207.08 m，平均采深為579.05 m，煤層厚度為3.40 m，平均傾角為3°，采用走向長壁全陷綜合機械化一次采全高方法開采。研究區上方布設觀測站進行地表巖移監測，點間距為19～38 m，監測日期為2017年5月23日至2017年10月30日。研究區域位置見圖1。

圖1 研究區域地理位置Fig.1 Geographical location of the study area

1.2 影像數據

本文選用歐空局Sentinel-1A衛星自2017年5月13日至2018年2月25日獲取的16景C波段升軌單視復數影像(SLC)，數據重訪周期12 d，空間分辨率為5 m(R)×20 m(A)，極化方式為VV，成像模式為Interferometric Wide Swath (IW)，原始數據寬總幅超過250 km，使用分辨率為90 m的SRTM3 Version 4數據來消除地形相位誤差。 SAR影像信息見表1。

表1 所用Sentinel-1A數據信息Table 1 Sentinel-1A used for data information

2 基于振幅的開采沉陷時序InSAR監測方法

2.1 振幅反演地表形變特征基本原理

煤礦開采引發地表破壞，造成地物目標的物理結構、表面粗糙度或地物目標類型發生變化，向散射能量也會發生相應的變化，這將導致SAR圖像的亮度和色澤發生變化[22]。通過雷達回波信號所攜帶的振幅數據，獲取不同地表散射目標的后向散射系數信息，數字化處理后得到像元亮度值(DN)，利用像元值表示開采沉陷過程中地表目標的后向散射變化來反演地表形變特征[23]。像元DN值與地表的后向散射系數σ0關系如圖2所示。

圖2 像元DN值與地表后向散射系數關系Fig.2 Relation between pixel DN value andsurface backscattering coefficient

像元(i,j)的強度或DN值與相應地面單元的后向散射的能量Pr的關系見式(1)[24]。

(1)

式中，a、b為比例常數。

雷達所接收的能量Pr與地面單元雷達后向散射系數σ0(i,j)成比例關系見式(2)。

Pr(i,j)=k(i,j)σ0(i,j)

(2)

聯合式(1)和式(2)得DN(i,j)與σ0(i,j)的關系式，見式(3)。

(3)

結合分布目標的雷達方程見式(4)。

(4)

由式(4)可得式(5)。

(5)

從式(5)可以看出k(i,j)與發射功率Pt、到地表散射目標的距離R、地面分辨率單元A、波長λ及天線的增益函數G有關。當k(i,j)在絕對或相對尺度上為已知值時，可以生成絕對尺度或相對尺度上的雷達回波圖象。

2.2 雷達影像回波強度的影響因素

1) 傳感器因素。雷達圖像表示地面雷達后向散射的估算值，在高后向散射區域大部分能量反射回雷達系統中表現淺色調，反之，為深色調。從雷達方程可知，雷達回波強度與入射波長直接相關。同時，雷達遙感系統所選擇的波長長短，決定了表面粗糙度的大小以及入射波穿透深度的能力，也直接影響到雷達回波的強弱；不同的極化方式也會導致目標對電磁波有不同的響應，使雷達回波強度不同，需要根據地面不同特點選擇極化方式[25]。除此之外，俯角、照射帶寬度也會影響回波強度。

2) 地面因素。影響回波強度的地面因素主要包括地表特性和負介電常數。開采沉陷過程中地表特性(坡度、坡向等)發生隨機變化，雷達圖像中會出現透視收縮、疊掩、陰影等現象[26]。同時，回波路徑的改變使原有地物的散射機制(面散射、體散射等)發生改變，進而影響回波強度。負介電常數由物質組成及溫度決定，其值越大，回波強度越強，礦區裸露的巖石、建(構)筑物反射能力強，在影像中表現為高亮度區域，植被覆蓋密集區表現灰暗[27]。

3 時序振幅變化及地表沉降速率分析

3.1 研究區時序后向散射系數分析

利用多時相影像處理技術，將研究區影像進行多視、配準、濾波、地理編碼和輻射定標后，獲得2017年5月13日至2018年2月25日期間新元礦某工作面時序后向散射系數圖(圖3)。從圖3可以看出，工作面附近的房屋等建(構)物相干性最高，反射微波能力強，后向散射系數值接近5，在影像中表現較為密集；其次，地表及山坡裸露的巖石相干性較高，反射微波能力較強，數量最多，在影像中分布不均且呈點狀分布，后向散射系數值為3～4。工作面內地物整體后向散射系數值較小，很難選取到密集高的相干點，但仍有部分高相干地物隨機出現。這是因為開采活動的劇烈進行，造成地表破壞，裸露的巖石或高相干地物發生出現-掩埋-再出現的現象，同時受衛星重訪周期和植被季節性變化的影響，無法獲取短暫暴露的高相干地物回波信號。利用開采活動影響范圍內地表后向散射系數的變化，反應整個沉陷盆地內地表沉降速率變化是本文的研究重點。

圖3 后向散射系數時間序列圖Fig.3 Backscattering coefficient time series diagram

3.2 礦區分析模型的建立與地表沉降速率分析

3.2.1 建立分析模型

煤礦開采活動誘發地表發生移動，形成移動盆地，各個區域的移動和變形的性質及大小不盡相同，分為中性區域、壓縮區域和拉伸區域[28]。為分析整個沉陷盆地的沉降特征，依據開采活動過程中工作面及附近地表變化造成其后向散射系數不同，結合地表達到充分采動時采空區的長度和寬度均要達到并超過1.2～1.4H0(H0為平均開采深度)[28]，考慮走向(C線)和傾向(A線)主斷面的監測點實際布設情況，將研究區整體劃分為工作面中心、工作面內邊緣與工作面外邊緣(圖4)。工作面內受開采沉陷影響嚴重，不同區域地表劇烈變化程度不同，為更加準確地研究工作面內地表后向散射系數與沉降速率的關系，利用ArcGIS軟件將工作面劃分為九塊相同大小的分區[29]。其中，工作面外邊緣監測點包含點A1～A20、B1～B14、C1～C19；工作面內邊緣1區包含點B30～B32，2區包含點B22～B29、點C20～C39，3區包含點B15～B21，4區包含點A35～A38，6區包含點A21～A25；工作面中心5區包含點A26～A34。在實際分析過程中，依據開采沉陷盆地的對稱性，選取走向和傾向線部分監測點的像元后向散射系數進行分析，進而反演整個沉陷盆地的地表后向散射系數變化。

圖4 開采沉陷分析模型圖Fig.4 Mining subsidence analysis model

3.2.2 礦區地表沉降速率分析

為了分析研究區的開采活動與后向散射系數的關系，在走向和傾向觀測線上分別以4個監測點為間隔提取出21個監測點內像元的后向散射系數。其中，走向線上的監測點主要集中在工作面外邊緣(C1、C5、C10、C15、C20)與內邊緣(C25、C30、C39)，為保證工作面走向信息的完整性，增加內、外邊緣分界點C19與中心點A30。傾向觀測線上選取工作面外邊緣點(A1、A5、A10、A15、A20)、內邊緣點(A25、A35、A38)與中心點A30,為保證沉降特征區域后向散射系數的特殊性，增加內、外邊緣區分界點A21和內邊緣、中心分界點A26進行研究。 利用簡單的折線繪制出時序后向散射系數圖(圖5和圖6)。

圖5 走向監測點時序后向散射系數Fig.5 Backscattering coefficient of strike monitoring points

由圖5可知，后向散射系數在0.10～0.36之間，整體呈低-高-低趨勢。走向觀測線各點在沉陷開始階段與衰退階段地表后向散射系數較低，2017年5月至2017年6月，地表后向散射系數較小但整體處于上升趨勢，該時間段采礦活動剛開始進行，地表移動變形出現滯后。2017年11月后地表后向散射系數減小并趨于穩定，該時間段采礦活動停止，地表存在微小的殘余變形。2017年7月至2017年11月期間地表后向散射系數較高呈起伏變化，該日期與開采活動日期基本吻合。同時，工作面中心A30、內外邊緣分界點C19、C20、靠近工作面中心點C39，在采礦活動進行時后向散射系數整體較高，且最大后向散射系數依次降低。該研究區域工作面中心地表沉降速率大，地表活動劇烈，實測沉降大于2.7 m；邊界兩側受采礦活動影響，地表活動存在差異。在內外邊緣區的其他各點后向散射系數較低，無明顯規律，同時這些區域受采礦活動影響小，地表移動變形不劇烈。

在傾向時序后向散射系數圖(圖6)中各點與走向各點規律相似，在采礦活動開始階段、活躍階段、衰退階段各點后向散射系數呈現低-高-低變化，沉陷活躍階段后向散射值較高且呈高低起伏變化。2017年5月至2017年6月地表后向散射系數較小但逐漸增加，2017年7月至2017年11月期間地表后向散射系數值較大且變化明顯，2017年11月后地表后向散射系數減小并趨于穩定。傾向線上所有監測點后向散射可分為兩類：①在工作面中心A30、中心與內邊緣分界點A25，后向散射系數整體較高，這兩個區域沉陷量極大，地表移動變形也較為劇烈；②其他各區域點后向散射系數整體較低，中心區域的最外側點A26值較大，最外邊緣點A1值最低。

圖6 傾向監測點時序后向散射系數Fig.6 Backscattering coefficients ofinclination monitoring points

由走向線和傾向線上監測點的后向散射系數變化情況可知，工作面邊緣區域受開采活動影響較小，后向散射系數較低；工作面中心和各分界點處采礦活動影響劇烈，地表后向散射系數值也較大。

進一步分析工作面中心區域(5區)可以發現其后向散射系數分三類(圖7)：①高后向散射系數區(A29-A33)；②中后向散射系數區(A26、A28)；③低后向散射系數區(A27、A34)。由沉陷盆地對稱特征分析可知，工作面中心區域后向散射系數整體由邊緣至中心逐漸增大。

圖7 工作面中心監測點時序后向散射系數圖Fig.7 Backscatter coefficients of monitoring pointsat the working face center

4 開采沉陷區地表沉降速率與后向散射系數相關性分析

由上述分析可知采礦過程中移動盆地不同區域內地物單元后向散射系數不同。為了驗證開采沉陷引起的地表變化速率與后向散射系數的相關性，評估多時相后向散射系數結果的可靠性，結合實測沉降速率，選取走向外邊緣點C1、內外邊緣分界點C20、內邊緣點C39、傾向外邊緣點A1、內邊緣與中心區域分界點A25、最大下沉點A29，分別進行走向和傾向監測點沉降速率與后向散射系數對比分析。

從圖8和圖9可知，只有在部分時段沉降速率與后向散射系數有一定關系，2017-05-13至2017-08-21時段，走向地表沉降速率與后向散射系數變化趨勢不同，傾向上二者變化趨勢一致；2017-08-21至2017-11-29時段，走向地表沉降速率與后向散射系數變化明顯且趨勢相同，傾向上地表沉降速率較小，對應后向散射系數起伏變化明顯，二者相關性較低。由此可知，2017-05-13至2017-08-21時段，走向監測點距推進面較遠，開采活動還未影響該區域地表沉降，地物長勢變化造成后向散射系數變化，造成地表沉降速率與后向散射系數相關性降低；2017-08-21至2017-11-29時段，傾向地表沉降處于衰退期，沉降速率變小，而地物長勢、地物類型變化疊加影響后向散射系數呈伏變化，造成地表沉降速率與后向散射系數相關性降低。由上可得，靠近工作面中心的大量級沉降區域后向散射系數變化主要由地表沉降速率變化引起。

圖8 走向監測點沉降速率與后向散射系數關系Fig.8 Relationship between settlement rate andbackscattering coefficient at strikemonitoring point

圖9 傾向監測點沉降速率與后向散射系數關系Fig.9 Relationship between settlement rate andbackscattering coefficient of inclinationmonitoring point

利用SPSS軟件對工作面特征點后向散射系數與沉降速率作相關性分析見表2。由表2知，走向外邊緣點C1、內外邊緣分界點C20、內邊緣點C39的后向散射系數與沉降速率的相關關系均大于0.850，平均相關性達到0.878，說明存在極強的相關性；傾向方向監測點后向散射系數與沉降速率相關性由外邊緣區至中心區域相關性逐漸增加，最大下沉點A29的相關性為0.960，最外邊緣點A1的相關系數為0.696，相關性較弱，工作面中心與內邊緣分界點A25的相關性達到0.836，相關性較強。傾向主斷面平均相關性為0.831略小于走向，說明傾向主斷面的點相關性整體較走向低，結合地表沉陷規律，發現煤礦開采引發地表移動變形使不同區域地表后向散射系數發生改變，同時地表后向散射系數可以作為描述整個沉陷盆地變化特征的指標。

表2 工作面特征點后向散射系數與沉降速率相關關系Table 2 Correlation between backscattering coefficient ofcharacteristic points of working face and settlement rate

5 結 論

1) 針對傳統InSAR利用相位信息難以解纏礦區大量級沉陷問題，本文提出利用多時相SAR振幅數據提取沉陷區地表后向散射系數分析礦區地表沉降速率規律的方法。利用時序振幅法得到了多時相后向散射系數與地表沉降速率的關系，并發現地表后向散射系數可以描述大量級沉陷區域的地表沉降速率變化。

2) 通過建立研究區后向散射系數分析模型，發現沉陷區地表后向散射系數由工作面邊緣至中心增大，與工作面走向和傾向地表沉降速率的平均相關系數r大于0.8；大量級沉陷區地表沉降速率與后向散射系數變化一致，相關性大，邊緣區地表后向散射系數會受地物類型、地物長勢及臨采影響，兩者相關性較低。