王道順，李偉*，孫建，鄭國棟

(1.山東省魯南地質(zhì)工程勘察院(山東省地質(zhì)礦產(chǎn)勘查開發(fā)局第二地質(zhì)大隊)，山東 兗州 272100；2.高分辨率對地觀測系統(tǒng)濟(jì)寧數(shù)據(jù)與應(yīng)用中心，山東 兗州 272100；3.自然資源部采煤沉陷區(qū)綜合治理工程技術(shù)創(chuàng)新中心，山東 兗州 272100)

0 引言

煤炭資源作為我國的重要能源之一，極大地推動了國家的經(jīng)濟(jì)和社會發(fā)展。隨著開采強(qiáng)度的日益增加，由此引起的生態(tài)環(huán)境問題更加凸顯[1-3]。煤炭開采造成的地表沉陷屬于典型地質(zhì)災(zāi)害之一，往往導(dǎo)致礦區(qū)及一定范圍內(nèi)的地表產(chǎn)生不同程度的沉降，對地面房屋設(shè)施、道路、農(nóng)田水利基礎(chǔ)設(shè)施等建(構(gòu))筑物造成不同程度的損毀，為經(jīng)濟(jì)發(fā)展和社會穩(wěn)定帶來重大隱患，因此對煤炭開采引起的地表沉降進(jìn)行監(jiān)測具有重要意義。

以往的礦區(qū)地表沉降監(jiān)測多采用布設(shè)巖移觀測站的方式，利用水準(zhǔn)測量、GNSS測量等傳統(tǒng)方法進(jìn)行監(jiān)測，此類方法監(jiān)測精度可觀，但作業(yè)成本較高，獲取的監(jiān)測數(shù)據(jù)量有限，不能完全反映礦區(qū)地表沉降特征，難以實現(xiàn)大范圍的地表沉降監(jiān)測。合成孔徑雷達(dá)差分干涉測量(D-InSAR)通過微波遙感主動測量的方式，可以克服傳統(tǒng)監(jiān)測方法的不足，且具有全天時、全天候作業(yè)的優(yōu)點，可獲取大范圍、高精度的空間延續(xù)性地表形變信息。已有學(xué)者將其應(yīng)用于礦區(qū)的地表形變監(jiān)測研究[4-9]，但該方法使用數(shù)據(jù)量較少，在長時間序列的監(jiān)測過程中容易受到時空失相干、大氣延遲等因素的影響，降低其監(jiān)測精度。在此基礎(chǔ)上發(fā)展而來的時序InSAR技術(shù)彌補(bǔ)了上述不足，極大地提高了形變監(jiān)測精度，已廣泛應(yīng)用于礦區(qū)形變監(jiān)測、城市沉降監(jiān)測等方面，其中以永久散射體技術(shù)(Permanent Scatterer InSAR，PS-InSAR)和短基線集方法(Small Baseline Subset InSAR，SBAS-InSAR)應(yīng)用最為普遍[10-15]。

PS-InSAR技術(shù)已成功應(yīng)用于礦區(qū)的地表沉降監(jiān)測，為了獲得在時間序列上散射特性穩(wěn)定的PS點，PS-InSAR往往需要20景以上的數(shù)據(jù)量參與計算。該方法基于大量PS點的迭代回歸計算，導(dǎo)致運算效率偏低[16]。礦區(qū)多位于郊區(qū)、農(nóng)田等低相干地區(qū)，若所用數(shù)據(jù)量偏少、時間間隔較長，利用PS-InSAR技術(shù)難以獲取足夠的PS點，導(dǎo)致無法完整揭示礦區(qū)沉降規(guī)律。本文在分析SBAS-InSAR技術(shù)方法的基礎(chǔ)上，以濟(jì)寧北部礦區(qū)為研究區(qū)域(圖1)，采用C波段Sentinel-1A雷達(dá)數(shù)據(jù)和SBAS-InSAR時序方法開展礦區(qū)地表沉降監(jiān)測和時序分析，獲取監(jiān)測范圍內(nèi)的地面沉降時空分布特征，可為采煤沉陷區(qū)的調(diào)查和治理提供技術(shù)支撐。

圖1 研究區(qū)域影像圖

1 SBAS-InSAR原理

SBAS-InSAR由Berardino等人[17]于2002年提出，不同于PS-InSAR的單主影像，該方法是一種基于多主影像的InSAR時間序列方法。它通過短基線原則，將大量SAR數(shù)據(jù)組合為具有多個主影像的干涉子集，每個子集內(nèi)的干涉對基線長度均低于臨界基線值，時間基線也盡可能短，集合間的SAR影像基線距大，通過這種方式克服了時間和空間上的失相關(guān)，因此SBAS-InSAR可以通過較少的數(shù)據(jù)量來獲取較可靠的監(jiān)測結(jié)果，且其監(jiān)測對象為分布式散射體，該方法更適合于礦區(qū)的地表形變監(jiān)測。

基于短基線集原則，SBAS-InSAR具有多個主影像，但仍需選取其中一景影像作為公共主影像進(jìn)行配準(zhǔn)；組成各干涉子集后，利用外部參考數(shù)字高程模型(Digital Elevation Model，DEM)數(shù)據(jù)模擬并去除每個干涉相對的地形相位，然后生成時間序列差分干涉圖集；相位解纏后得到每個相干目標(biāo)的相位信息，包括形變相位、大氣延遲相位、軌道誤差相位等信息，各誤差相位可在時間序列上采用濾波方法或多項式模型予以去除；由于SBAS-InSAR具有多個主影像，各干涉子集在聯(lián)合求解時容易出現(xiàn)方程秩虧現(xiàn)象，因此引入奇異值分解方法，利用最小二乘原理得到地表時間序列形變信息。

2 研究區(qū)域與數(shù)據(jù)源

2.1 研究區(qū)域

濟(jì)寧市位于山東省魯西南地區(qū)，礦產(chǎn)資源豐富，是魯西煤田基地的重要組成部分，是全國重要的煤炭基地之一[18-20]。本次實驗研究區(qū)域位于濟(jì)寧市北部地區(qū)，覆蓋范圍大致為E116°28′01″～116°44′34″、N35°29′40″～35°42′35″，研究區(qū)域內(nèi)礦區(qū)集中分布，主要有義橋煤礦、唐陽煤礦、義能煤礦、魯西煤礦、新驛煤礦等礦區(qū)(圖1)。

2.2 數(shù)據(jù)源

本文使用覆蓋研究區(qū)域2019年9月19日—2020年3月29日間的17景Sentinel-1A C波段單視復(fù)數(shù)(SLC)數(shù)據(jù)，該數(shù)據(jù)為干涉寬幅(IW)成像模式，極化方式VV，空間分辨率5m×20m(方位向×距離向)，具體信息見表1。Sentinel-1A衛(wèi)星在軌運行高度693km，重訪周期12d，若與Sentinel-1B衛(wèi)星組網(wǎng)觀測，可將重訪周期縮短至6d，可有效提高監(jiān)測區(qū)域的時間相干性，此外由于Sentinel-1A衛(wèi)星采用先進(jìn)的軌道控制技術(shù)，其干涉對空間基線最大為162m，多數(shù)都在100m范圍內(nèi)[21]，遠(yuǎn)遠(yuǎn)小于其空間臨界基線值，因此可以忽略其干涉對空間基線的差異。外部參考DEM選擇由日本太空發(fā)展署(JAXA)免費分發(fā)的ALOS World 3D數(shù)據(jù)，空間分辨率30m，用于模擬去除地形相位和地理編碼。

表1 研究區(qū)Sentinel-1A影像參數(shù)

3 礦區(qū)形變監(jiān)測與分析

3.1 短基線集干涉組合

實驗選擇2019年11月30日的影像作為公共主影像，其余影像作為從影像與主影像進(jìn)行配準(zhǔn)并采樣到主影像像素空間。設(shè)置空間基線閾值100m，時間基線閾值60d，共得到55組干涉相對，其中最大空間基線87.47m，最小空間基線4.11m，時空基線連接情況見圖2。

圖2 時空基線圖

3.2 干涉圖生成與相位解纏

將得到的55組像對進(jìn)行差分干涉處理，得到時間序列干涉圖集，同時生成對應(yīng)的相干系數(shù)圖，由于設(shè)置的閾值相對嚴(yán)格，多數(shù)差分干涉圖干涉效果較好，可以看到明顯的形變干涉條紋。圖3是2019年11月30日—2019年12月24日內(nèi)生成的差分干涉結(jié)果(SAR坐標(biāo)系)，圖3a為相干系數(shù)圖，顏色越亮代表相干性越高；圖3b為差分干涉圖，可以看出，監(jiān)測區(qū)域內(nèi)有多處明顯的形變干涉條紋。相位解纏選擇最小費用流(MCF)解纏算法，相干性閾值設(shè)為0.2，低于相干閾值的像元進(jìn)行掩膜處理，不參與計算，然后逐個干涉圖進(jìn)行相位解纏，圖3c為對應(yīng)的相位解纏圖，黑色代表相干性過低而沒有解纏結(jié)果的地方。

3.3 軌道精煉和重去平

對得到的55組干涉像對結(jié)果進(jìn)行目視檢查，剔除相干性過低或解纏錯誤的像對。選取一組干涉效果良好的像對作為參考影像，選擇遠(yuǎn)離形變區(qū)域且相干性較高的若干穩(wěn)定地面點作為控制點，使用多項式模型進(jìn)行軌道精煉，去除殘余的恒定相位和相位坡道。

3.4 獲取監(jiān)測結(jié)果

上述步驟完成后反演區(qū)域形變速率，同時對干涉圖進(jìn)行二次解纏優(yōu)化。在第一次形變速率反演的基礎(chǔ)上，利用大氣相位在時間上的高通濾波和空間上的低通濾波進(jìn)行去除，獲取時間序列上的位移結(jié)果(圖4)。

a—相干系數(shù)圖；b—差分干涉圖；c—解纏相位圖圖3 干涉對相干系數(shù)圖、差分干涉圖、解纏相位圖

圖4 監(jiān)測區(qū)域年均沉降速率圖

3.5 結(jié)果分析

從監(jiān)測結(jié)果可以看出，監(jiān)測區(qū)域整體基本處于穩(wěn)定狀態(tài)，濟(jì)寧北部礦區(qū)內(nèi)明顯分布有6處不同沉降程度的區(qū)域(A—F)，分別位于義橋煤礦、唐陽煤礦、義能煤礦、新驛煤礦、魯西煤礦、何崗煤礦，其中C區(qū)域沉降最大，年平均沉降速率達(dá)到-242mm/y，呈現(xiàn)出雙沉降漏斗特征；A、B兩處沉降區(qū)域相距僅640m，已呈現(xiàn)接連成片沉降趨勢，沉降范圍最大；D區(qū)域最大沉降速率達(dá)到-85.7mm/y；E，F(xiàn)兩處區(qū)域沉降范圍最小，E區(qū)域最大沉降速率為-70.1mm/y，F(xiàn)區(qū)域最大沉降速率達(dá)到-91.6mm/y。部分沉降區(qū)域中心不同程度的出現(xiàn)了“監(jiān)測空值”，主要是礦區(qū)沉降過大，超出了InSAR可探測的最大形變梯度，致使沉降中心出現(xiàn)失相干，但其沉降邊緣范圍探測不受影響。

4 結(jié)論

(1)基于2019年9月19日—2020年3月29日的17景Sentinel-1A數(shù)據(jù)，運用SBAS-InSAR技術(shù)獲取了濟(jì)寧北部礦區(qū)的年均沉降速率和沉降漏斗分布。

(2)濟(jì)寧北部礦區(qū)內(nèi)分布有6處明顯沉降區(qū)域，分布于義橋煤礦、唐陽煤礦、義能煤礦、新驛煤礦、魯西煤礦和何崗煤礦，形變范圍與礦區(qū)實際空間分布位置一致，其中義能煤礦沉降特征最為明顯，年均沉降速率達(dá)到-242mm/y。

(3)受形變梯度影響，部分沉降區(qū)域中心出現(xiàn)“監(jiān)測空值”，但沉降范圍仍可以完整探測。相比較于C波段數(shù)據(jù)，L波段數(shù)據(jù)波長更長，在低相干區(qū)域更容易穿透地表植被覆蓋，得到更高質(zhì)量的干涉結(jié)果和更大中心沉降量。下一步工作將收集該區(qū)域L波段數(shù)據(jù)，開展更高質(zhì)量的InSAR礦區(qū)地表形變監(jiān)測，為礦區(qū)的地面沉降防治提供決策依據(jù)。