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

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

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

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

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

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

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

1 研究區(qū)與數(shù)據(jù)

1.1 研究區(qū)概況

陽煤新元礦位于沁水煤田西北部，太原至陽泉的中間地帶，壽陽-陽泉構(gòu)造堆積盆地地區(qū)，海拔1 058.3～1 110.9 m，屬黃土高原溫帶季風(fēng)氣候區(qū)，天氣寒冷干燥，溫差變化大，年均降水量約500 mm，多集中在6～9月。 研究區(qū)地表大部分被第三系、第四系松散地層覆蓋，經(jīng)長期風(fēng)化剝蝕，溝谷縱橫，地形復(fù)雜。礦區(qū)南北長4～9 km，東西長15～16 km，地處東經(jīng)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°，采用走向長壁全陷綜合機械化一次采全高方法開采。研究區(qū)上方布設(shè)觀測站進行地表巖移監(jiān)測，點間距為19～38 m，監(jiān)測日期為2017年5月23日至2017年10月30日。研究區(qū)域位置見圖1。

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

1.2 影像數(shù)據(jù)

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

表1 所用Sentinel-1A數(shù)據(jù)信息Table 1 Sentinel-1A used for data information

2 基于振幅的開采沉陷時序InSAR監(jiān)測方法

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

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

圖2 像元DN值與地表后向散射系數(shù)關(guān)系Fig.2 Relation between pixel DN value andsurface backscattering coefficient

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

(1)

式中，a、b為比例常數(shù)。

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

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

(2)

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

(3)

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

(4)

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

(5)

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

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

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

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

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

3.1 研究區(qū)時序后向散射系數(shù)分析

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

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

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

3.2.1 建立分析模型

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

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

3.2.2 礦區(qū)地表沉降速率分析

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

圖5 走向監(jiān)測點時序后向散射系數(shù)Fig.5 Backscattering coefficient of strike monitoring points

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

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

圖6 傾向監(jiān)測點時序后向散射系數(shù)Fig.6 Backscattering coefficients ofinclination monitoring points

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

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

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

4 開采沉陷區(qū)地表沉降速率與后向散射系數(shù)相關(guān)性分析

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

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

圖8 走向監(jiān)測點沉降速率與后向散射系數(shù)關(guān)系Fig.8 Relationship between settlement rate andbackscattering coefficient at strikemonitoring point

圖9 傾向監(jiān)測點沉降速率與后向散射系數(shù)關(guān)系Fig.9 Relationship between settlement rate andbackscattering coefficient of inclinationmonitoring point

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

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

5 結(jié) 論

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

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