姬宗皓，薄懷志，王東平，高彥生，宋炳忠

(1.濟(jì)寧市采煤塌陷地治理中心，山東 濟(jì)寧 272000；2.自然資源部采煤沉陷區(qū)綜合治理與生態(tài)修復(fù)工程技術(shù)創(chuàng)新中心，山東 濟(jì)寧 272100；3.山東省魯南地質(zhì)工程勘察院，山東 濟(jì)寧 272100；4.山東省地礦局采煤塌陷地綜合治理工程技術(shù)研究中心，山東 濟(jì)寧 272100)

0 引言

濟(jì)寧是魯西煤炭基地的重要組成部分，煤田面積3920km2，主要由兗州煤田、濟(jì)寧煤田、寧汶煤田、滕縣煤田、金鄉(xiāng)煤田、梁山煤田、巨野煤田組成，煤系地層主要為石炭-二疊紀(jì)月門溝群太原組和山西組[1]。經(jīng)過(guò)幾十年高強(qiáng)度開采，濟(jì)寧出現(xiàn)了大面積的采煤沉陷。采煤沉陷造成嚴(yán)重的地質(zhì)環(huán)境問(wèn)題和社會(huì)問(wèn)題，主要表現(xiàn)在：耕地保有量銳減、嚴(yán)重破壞礦區(qū)生態(tài)環(huán)境、引發(fā)諸多復(fù)雜社會(huì)問(wèn)題、嚴(yán)重制約城市建設(shè)發(fā)展等[2-3]。查清濟(jì)寧市采煤沉陷現(xiàn)狀，對(duì)于采煤塌陷地監(jiān)測(cè)監(jiān)管意義重大[4]。與土地利用現(xiàn)狀變更不同，開采沉陷主要表現(xiàn)為垂向空間變化，大范圍監(jiān)測(cè)難度較大。目前對(duì)于單個(gè)開采工作面的沉陷監(jiān)測(cè)多使用傳統(tǒng)的測(cè)量方法，即布設(shè)沉陷監(jiān)測(cè)點(diǎn)，使用水準(zhǔn)儀或GNSS-RTK定期監(jiān)測(cè)沉陷情況[5]。該方法測(cè)量精度較高，但觀測(cè)周期長(zhǎng)、費(fèi)用高、觀測(cè)點(diǎn)離散分布難以反映大范圍全面形變情況和連續(xù)形變規(guī)律。秦曉敏等[6]提出利用多時(shí)相DEM數(shù)據(jù)高程變化、結(jié)合高分辨率遙感影像、礦區(qū)開采資料提取開采沉陷范圍的方法，但精度較低，難以監(jiān)測(cè)沉陷量小于0.5m的塌陷地。

D-InSAR技術(shù)具有全天候、大范圍、精度高等優(yōu)點(diǎn)，近年來(lái)發(fā)展迅速，已經(jīng)在地面沉降監(jiān)測(cè)、地質(zhì)災(zāi)害防治等領(lǐng)域得到了廣泛的應(yīng)用[7-9]。將D-InSAR技術(shù)應(yīng)用于濟(jì)寧市大范圍采煤沉陷區(qū)動(dòng)態(tài)監(jiān)測(cè)，對(duì)開采沉陷監(jiān)管、治理規(guī)劃設(shè)計(jì)、政府決策等方面都具有重要意義。RADARSAT-2是一顆搭載C波段傳感器的高分辨率商用雷達(dá)衛(wèi)星，由加拿大太空署與MDA公司合作，于2007年12月14日在哈薩克斯坦拜科努爾基地發(fā)射升空。該文使用RADARSAT-2的數(shù)據(jù)產(chǎn)品SLC(單視復(fù)型數(shù)據(jù))，研究了D-InSAR技術(shù)在濟(jì)寧市大范圍采煤沉陷區(qū)監(jiān)測(cè)及其時(shí)序變化特征分析中的應(yīng)用。

1 D-InSAR采煤沉陷監(jiān)測(cè)的原理

D-InSAR技術(shù)是由InSAR技術(shù)發(fā)展而來(lái)的，它是以合成孔徑雷達(dá)復(fù)數(shù)圖像的相位信息獲取地表變化信息的技術(shù)。根據(jù)成像時(shí)間分類，InSAR可以分為單次軌道和重復(fù)軌道2種模式。單次軌道干涉是指在同一機(jī)載或星載平臺(tái)上裝載兩幅天線，其中一幅天線發(fā)射信號(hào)，兩幅天線都接受地面回波信號(hào)，并利用獲取的數(shù)據(jù)進(jìn)行干涉處理。重復(fù)軌道干涉是指同一傳感器或相似傳感器按照平行軌道2次對(duì)地成像，分別發(fā)、收信號(hào)，利用得到的數(shù)據(jù)進(jìn)行干涉處理。目前，常用的D-InSAR地表形變監(jiān)測(cè)通常為星載單天線重復(fù)軌道模式。

1.1 InSAR測(cè)量原理

假設(shè)雷達(dá)傳感器S1、S2在不同時(shí)間、兩次相互平行地對(duì)同一區(qū)域重復(fù)觀測(cè)，生成兩幅影像(主影像、從影像)，S1、S2的距離為空間基線，按雷達(dá)信號(hào)入射方向，基線B可投影為平行分量和垂直分量，θ為雷達(dá)側(cè)視角，α為空間基線水平傾角，傳感器S1到地面分辨單元P的斜距為r，與S2到該點(diǎn)的斜距差為Δr，如圖1所示。P在主(從)影像上的灰度、相位值可表示為：

主影像和從影像經(jīng)配準(zhǔn)后共軛相乘生成干涉圖，該分辨單元的干涉相位值為：

ψ=arctan(u1·u2*)=ψ1-ψ2ψ∈[-π,π]

SAR系統(tǒng)觀測(cè)值ψi(i=1,2)與雷達(dá)波至地面單元的距離有關(guān)，同時(shí)受地面分辨單元散射特性的影響。假定2次成像時(shí)的地表散射特性相同，則干涉相位可表示為：

圖1 InSAR幾何示意圖

1.2 D-InSAR沉陷監(jiān)測(cè)基本原理

如圖2所示，兩次成像期間地表像元P由于地下開采發(fā)生了沉陷ΔH,且沉陷量在雷達(dá)視線方向上的形變分量為Δd，為了獲取形變相位，需將地形相位從干涉相位中去除，采用外部DEM差分之前先地理編碼至雷達(dá)坐標(biāo)系，與SAR影像精密配準(zhǔn)，其模擬的地形相位[10-12]：

從干涉相位中去除模擬地形相位后，忽略大氣延遲及噪聲誤差，得到形變相位：

當(dāng)形變相位變化2π時(shí)，對(duì)應(yīng)P的形變?yōu)椋?/p>

d2π為形變模糊度，與SAR成像波長(zhǎng)有關(guān)。獲取P點(diǎn)在雷達(dá)視線方向上的形變分量Δd后，即可根據(jù)ΔH與Δd的幾何關(guān)系計(jì)算得到P點(diǎn)的沉陷量ΔH。

圖2 D-InSAR幾何示意圖

2 采煤沉陷區(qū)范圍的提取

2.1 InSAR數(shù)據(jù)的選擇

根據(jù)開采沉陷的特點(diǎn)，綜合考慮存檔數(shù)據(jù)的覆蓋情況、時(shí)間序列完整性、空間基線等因素，結(jié)合濟(jì)寧煤炭開采范圍，選取15期RADARSAT-2 Wide模式數(shù)據(jù)進(jìn)行開采沉陷監(jiān)測(cè)，詳細(xì)數(shù)據(jù)參數(shù)見(jiàn)表1[13-14]。

表1 數(shù)據(jù)獲取參數(shù)

2.2 數(shù)據(jù)處理流程

根據(jù)以往InSAR數(shù)據(jù)處理的經(jīng)驗(yàn)，采用D-InSAR多基線累積結(jié)果作為濟(jì)寧采煤沉陷的形變監(jiān)測(cè)結(jié)果。多基線干涉相位累積疊加方法以兩軌差分為基礎(chǔ)，選取相鄰影像組合干涉對(duì)，即前一個(gè)干涉對(duì)的輔影像為后一個(gè)干涉對(duì)的主影像，組成首尾相連的干涉對(duì)組合，對(duì)每一個(gè)干涉對(duì)采用兩軌法進(jìn)行處理得到解纏相位，將各個(gè)干涉對(duì)的解纏相位連續(xù)累加并轉(zhuǎn)化為形變量，由此求取的形變即是對(duì)應(yīng)觀測(cè)時(shí)段內(nèi)的累積形變量。該方法利用一系列解纏的差分干涉圖，通過(guò)累積的方法得到相應(yīng)的形變量。處理時(shí)對(duì)獨(dú)立干涉像對(duì)大氣相位進(jìn)行了平均，達(dá)到抑制大氣噪聲的目的。

D-InSAR數(shù)據(jù)處理主要包括主輔影像預(yù)處理、影像配準(zhǔn)及重采樣、干涉圖生成、基線估計(jì)、地形相位差分、差分圖濾波、相位解纏、地理編碼等內(nèi)容[15-16]。

由于采用SAR數(shù)據(jù)為C波段成像，時(shí)間去相干比較嚴(yán)重，為保證干涉的質(zhì)量，對(duì)滿足時(shí)間基線小于200天，空間垂直基線小于300m的34個(gè)干涉對(duì)進(jìn)行干涉、地形相位差分、濾波、解纏等批處理。所選干涉對(duì)的時(shí)間、空間垂直基線見(jiàn)表2。

表2 干涉對(duì)的時(shí)間、空間垂直基線數(shù)據(jù)

為減小SAR成像系統(tǒng)噪聲，在基于強(qiáng)度的SAR影像配準(zhǔn)過(guò)程對(duì)距離向、方位向進(jìn)行1∶5多視，以20140617為主影像進(jìn)行影像配準(zhǔn)，采用強(qiáng)度互相關(guān)的多級(jí)匹配方法，距離向及方位向配準(zhǔn)精度均在0.05像元以內(nèi)，滿足干涉要求。主輔影像重采樣后進(jìn)行干涉處理，得到相干圖。采用分辨率30m的SRTM1 DEM，結(jié)合衛(wèi)星星歷數(shù)據(jù)、成像幾何參數(shù)等與主影像強(qiáng)度圖進(jìn)行精配準(zhǔn)、重采樣，配準(zhǔn)精度優(yōu)于0.1像元(距離向0.06，方位向0.08)，由精配準(zhǔn)并編碼至SAR坐標(biāo)系的高程進(jìn)行地形相位模擬，并從原始干涉圖中去除地形相位，得到差分干涉圖，如圖3所示[17]。

圖3 34個(gè)干涉對(duì)對(duì)應(yīng)的差分干涉圖

經(jīng)解纏處理恢復(fù)真實(shí)的差分相位。解纏參考點(diǎn)位于穩(wěn)定區(qū)域，采用相干性閾值0.3～0.4對(duì)低相干區(qū)域掩膜后使用最小費(fèi)用流法整塊解纏得到初始解纏相位圖。解纏前只對(duì)少數(shù)幾個(gè)差分圖利用條紋頻率法進(jìn)行基線估計(jì)，可以看出大部分解纏相位中含有軌道殘余相位，除此還存在幾處解纏誤差，由局部相干性過(guò)低造成；采用二次曲面擬合的方法去除軌道殘余相位，并對(duì)地表水汽相位采用高程相關(guān)的線性模型進(jìn)行擬合，最后采用最鄰近插值法對(duì)掩膜區(qū)的解纏斷裂進(jìn)行空間插值，得到最終解纏相位，如圖4所示。

圖4 34個(gè)干涉對(duì)對(duì)應(yīng)的最終解纏圖

2.3 結(jié)果分析

采用多基線累積疊加結(jié)果作為該次礦區(qū)沉陷的InSAR監(jiān)測(cè)結(jié)果。根據(jù)多基線累積疊加方法相關(guān)分析，利用未插值累積疊加結(jié)果用于提取超過(guò)InSAR形變監(jiān)測(cè)能力的“空值區(qū)”，“空值區(qū)”包含在5cm范圍線內(nèi)；利用插值累積疊加結(jié)果提取塌陷5cm范圍線。同時(shí)，結(jié)合現(xiàn)場(chǎng)調(diào)查、礦山采掘工程平面圖、概率積分法沉陷預(yù)測(cè)確定塌陷地范圍[18](圖5)。圖中黑色多邊形區(qū)域?yàn)榈V區(qū)范圍線，在礦區(qū)范圍內(nèi)監(jiān)測(cè)到多個(gè)采空區(qū)的塌陷活動(dòng)，采空區(qū)塌陷在空間上表現(xiàn)為漏斗形狀，與礦區(qū)開采引起的地表塌陷空間特征符合的較好。

圖5 基于累積形變提取的塌陷5cm范圍線

2.4 精度分析

為了驗(yàn)證D-InSAR開采沉陷監(jiān)測(cè)的可靠性，選取了某煤礦巖移觀測(cè)數(shù)據(jù)與D-InSAR監(jiān)測(cè)數(shù)據(jù)進(jìn)行對(duì)比。為了真實(shí)反映D-InSAR精度，需保證InSAR數(shù)據(jù)獲取時(shí)間與巖移觀測(cè)時(shí)間基本一致，共選取4期數(shù)據(jù)、150個(gè)巖移觀測(cè)點(diǎn)進(jìn)行對(duì)比，結(jié)果見(jiàn)表3、圖6。

表3 D-InSAR精度檢核

圖6 精度檢核較差范圍圖

3 采煤沉陷區(qū)時(shí)序變化特征

多基線累積疊加在得到礦區(qū)總體形變信息的同時(shí)，也可很好地反映礦區(qū)形變發(fā)生過(guò)程[19]。以花園煤礦為例，分析礦區(qū)地面沉陷時(shí)序變化情況。該處相干性較好，沉降特征明顯，并且沉降中心較為連續(xù)。區(qū)域累積沉陷結(jié)果見(jiàn)圖7。

圖7 區(qū)域累積沉陷結(jié)果圖

分析其沉陷過(guò)程，2013年4月—6月間礦區(qū)形變較為穩(wěn)定；2013年6月—9月出現(xiàn)形變加劇的趨勢(shì)；2013年9月—2014年6月間形變較為平穩(wěn)；2014年6月—9月間形變進(jìn)一步加??；2014年9月—2015年5月間形變平穩(wěn)；2015年5月—10月，形變趨勢(shì)進(jìn)一步加劇，沉陷呈現(xiàn)為非線性逐步增大的過(guò)程。結(jié)合地質(zhì)采礦及工作面布設(shè)情況，經(jīng)分析，A1點(diǎn)周期性間隔沉陷主要是受到臨近工作面開采沉陷影響。A1點(diǎn)沉陷時(shí)序變化曲線見(jiàn)圖8、區(qū)域沉陷時(shí)序變化見(jiàn)圖9。

圖8 A1形變結(jié)果時(shí)序變化曲線圖

圖9 區(qū)域沉陷時(shí)序變化圖

4 結(jié)論

該文闡述分析了使用15期RADARSAT-2 Wide模式數(shù)據(jù)，采用D-InSAR技術(shù)進(jìn)行濟(jì)寧市大范圍采煤沉陷監(jiān)測(cè)的方法和流程，得出如下結(jié)論：

(1)D-InSAR技術(shù)可以精確地提取出采煤沉陷區(qū)范圍(5cm下沉線)，結(jié)合現(xiàn)場(chǎng)調(diào)查和概率積分法沉陷預(yù)測(cè)，得到了濟(jì)寧市采煤沉陷現(xiàn)狀，為采煤塌陷地監(jiān)測(cè)監(jiān)管、治理規(guī)劃提供基礎(chǔ)數(shù)據(jù)。

(2)多基線累積疊加在得到最終沉陷結(jié)果的同時(shí)，也能較好地反映礦區(qū)形變的發(fā)生過(guò)程，全面反映礦區(qū)沉陷形態(tài)，可廣泛應(yīng)用于沉陷監(jiān)測(cè)、概率積分法參數(shù)反演等工作中。

(3)InSAR可以實(shí)現(xiàn)大范圍、高精度的開采沉陷監(jiān)測(cè)，節(jié)省大量人力物力；歷史存檔數(shù)據(jù)豐富，可應(yīng)用于老采空區(qū)殘余變形監(jiān)測(cè)、穩(wěn)定性評(píng)價(jià)等工作中[20]。