基于SBAS 技術(shù)的煤礦采空區(qū)探查與穩(wěn)定性評價

張 豐，刁鑫鵬，譚秀全，李愛軍，楊 雋

（1.山東省魯南地質(zhì)工程勘察院（山東省地質(zhì)礦產(chǎn)勘查開發(fā)局第二地質(zhì)大隊(duì)），山東 濟(jì)寧 272100；2.自然資源部采煤沉陷區(qū)綜合治理工程技術(shù)創(chuàng)新中心，山東 濟(jì)寧 272100；3.中國礦業(yè)大學(xué)環(huán)境與測繪學(xué)院，江蘇 徐州 221116；4.江蘇省資源環(huán)境信息工程重點(diǎn)實(shí)驗(yàn)室，江蘇 徐州 221116）

0 引 言

煤炭作為我國的主體能源，在國民經(jīng)濟(jì)和社會發(fā)展中起著至關(guān)重要的作用。 濟(jì)寧市煤炭資源豐富，開采歷史悠久，轄區(qū)內(nèi)查明的煤炭資源儲量達(dá)140 億t；經(jīng)過長期的開采，地下已形成面積龐大的采空區(qū)。 近年來，城市建設(shè)快速發(fā)展，土地資源日益緊張，煤礦沉陷區(qū)土地成為可持續(xù)發(fā)展的重要資源，故而導(dǎo)致地下采空與地面建設(shè)間的矛盾日益尖銳［1］。 因此，針對該地區(qū)進(jìn)行高精度和大覆蓋面積的采煤沉陷地調(diào)查與采空區(qū)穩(wěn)定性評價工作，有利于緩解兩者間的矛盾，對城市的可持續(xù)發(fā)展亦具有重要的現(xiàn)實(shí)意義。

傳統(tǒng)地表形變監(jiān)測技術(shù)（水準(zhǔn)測量、GPS、全站儀等）的測量精度高，但卻存在監(jiān)測周期長、成果不連續(xù)、空間覆蓋率低的缺陷［2－3］，不易于實(shí)現(xiàn)大范圍采煤沉陷地的調(diào)查。 而當(dāng)前三維激光掃描與無人機(jī)低空攝影測量的形變監(jiān)測精度相對較低，難以滿足采空區(qū)穩(wěn)定性評價的要求［4］。 隨著高分辨率、短重訪周期SAR 衛(wèi)星的連續(xù)升空以及數(shù)據(jù)處理手段的不斷發(fā)展，InSAR 獲取地表形變的精度已達(dá)毫米級［5］，使該技術(shù)在煤礦區(qū)地表變形監(jiān)測領(lǐng)域得到廣泛應(yīng)用［6－13］。

筆者選用短基線集技術(shù)（small baseline subset，SBAS），通過2018 年9 月至2019 年3 月的14 景Sentinel－1 數(shù)據(jù)開展?jié)鷮幊鞘幸?guī)劃區(qū)的地面形變解譯，獲取了研究區(qū)地表形變的時間演變過程和空間分布特征；并結(jié)合相應(yīng)時間段內(nèi)的地表實(shí)測數(shù)據(jù)進(jìn)行了解譯精度的評價。 研究表明，SBAS 能夠?qū)崿F(xiàn)采煤沉陷地毫米級精度的形變監(jiān)測，且能夠顯著區(qū)分不同開采條件下地表殘余變形的特征。 其成果有利于提高煤炭資源采空區(qū)地質(zhì)調(diào)查的有效性，能夠?yàn)椴煽諈^(qū)穩(wěn)定性評價提供數(shù)據(jù)參考。

1 SBAS 基本原理

SBAS 技術(shù)的提出是為削弱時空失相干影響，提高InSAR 地表形變解譯精度。 BERARDINO、LA?NARI 等［14－15］先后探討了該技術(shù)在地表形變監(jiān)測中的應(yīng)用，而后HOOPER 等［16］又對處理過程中高相干點(diǎn)的選取方法進(jìn)行了優(yōu)化。

SBAS 技術(shù)的基本思想是通過設(shè)置時、空基線條件，將同一研究區(qū)的多期SAR 影像分組成若干個短基線集合；利用最小二乘準(zhǔn)則，獲取每個短基線集合的地表形變時間序列；通過奇異值分解方法再將分組后的短基線集進(jìn)行聯(lián)合求解；進(jìn)而得到研究區(qū)整個觀測時間段內(nèi)時序的地表形變信息［17］。 具體可表述［18－19］為：假定同一研究區(qū)成像時間為t0，t1，t2，…，tn－1的N幅SAR 影像，在短基線距的條件下可形成M幅差分干涉圖，且滿足N／2≤M≤N（N－1）／2。 對于干涉圖i，在去除平地相位后，干涉圖中任意位置（x，y）的干涉相位可表示為：

其中，tA、tB為干涉圖i所對應(yīng)SAR 影像的獲取時間；δφdef為tA～tB時間段內(nèi)地表視線向形變相位；δφε為目標(biāo)區(qū)地形相位；δφα、δφn分別為大氣相位和噪聲相位。 其中前3 項(xiàng)可以表示為

式中，λ為雷達(dá)波長；R為斜距；B⊥為垂直基線；Δz為DEM 高程差；θ為入射角。

若不同干涉圖間的形變速率為vk，k＋1，則tA～tB的累積形變量可表達(dá)為

對M幅干涉圖進(jìn)行三維時空相位解纏和地理編碼后，便可求出不同SAR 影像獲取時間的形變速率。

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

2.1 研究區(qū)概況

研究區(qū)域?yàn)闈?jì)寧城市規(guī)劃區(qū)，面積約為950 km2，包括任城區(qū)、兗州區(qū)，鄒城市、曲阜市、嘉祥縣的街道和部分鄉(xiāng)鎮(zhèn)。 研究區(qū)內(nèi)礦井眾多，井田范圍與城市規(guī)劃區(qū)重疊面積約559.36 km2。 研究區(qū)的空間位置、范圍和各礦井分布情況如圖1所示。

圖1 濟(jì)寧城市規(guī)劃區(qū)空間位置、范圍及礦井分布情況Fig.1 Spatial location，scope and coal minesdistribution of Jining urban planning area

據(jù)統(tǒng)計，僅2018 年上述各礦井的煤炭資源產(chǎn)出量即超過3 500 萬t。 煤炭資源的長期大范圍開采，導(dǎo)致地下已形成面積龐大的采空區(qū)。 為緩解城市建設(shè)與地下采空之間的矛盾，開展采煤沉陷地的調(diào)查與老采空區(qū)的穩(wěn)定性評價工作迫在眉睫。

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

用于研究區(qū)形變信息提取的SAR 數(shù)據(jù)為2018年9 月至2019 年3 月的14 期Sentinel－1 影像，C 波段，地面分辨率為5 m×20 m，相關(guān)參數(shù)見表1。 Sen?tinel－1 影像的軌道參數(shù)需根據(jù)影像的成像時間下載，歐空局提供了哨兵影像的Precise Orbit Ephem?erides（POD 精密定軌星歷數(shù)據(jù)），此外，為消除地形相位，外部DEM 采用的是美國宇航局SRTM3 數(shù)據(jù)。

表1 濟(jì)寧Sentinel－1 影像數(shù)據(jù)相關(guān)參數(shù)Table 1 Relevant parameters of sentinel－1 images in Jining

3 數(shù)據(jù)處理與形變解譯過程

SBAS 的數(shù)據(jù)處理主要包括：差分干涉對的生成、點(diǎn)目標(biāo)選取、差分干涉圖的解纏，以及時間形變序列的獲取等幾個步驟，詳細(xì)的處理流程如圖2所示。

圖2 SBAS 地表形變解譯流程Fig.2 Flowchart of SBAS method

本次SBAS 地表形變解譯采用的是SARscape軟件。 數(shù)據(jù)處理過程中參數(shù)設(shè)置如下：①連接圖生成時，為充分利用影像數(shù)據(jù)，時間基線設(shè)置為200 d，空間基線設(shè)置為臨界基線的10%；共連接成91 個像對，其中時間基線最長為180 d，空間基線最長為178.49 m（平均62.32 m）。 ②干涉處理過程中距離向與方位向的視數(shù)比設(shè)置為4 ∶1，影像配準(zhǔn)利用精密軌道數(shù)據(jù)采用強(qiáng)度互相關(guān)算法，濾波方法選用自適應(yīng)濾波法，解纏方法采用最小費(fèi)用流法（相干性閾值設(shè)置為0.3）。 ③軌道精煉和重去平處理過程中，地面控制點(diǎn)選擇相干性高、相位好的點(diǎn)，并且所選控制點(diǎn)位置的形變需為0。④第一次解譯反演選用較為穩(wěn)定的線性模型。 其余參數(shù)均選用SARscape 處理Sentinel－1 數(shù)據(jù)的推薦參數(shù)。

最終解譯的地表形變速率分布和各時間段內(nèi)的地表垂直位移情況分別如圖3 和圖4 所示。

圖3 研究區(qū)地表形變速率分布Fig.3 Distribution of surface deformation rate in study area

4 監(jiān)測結(jié)果分析與精度評價

4.1 監(jiān)測結(jié)果分析

圖3 和圖4 分別展示了沉陷速率大于10 mm／a的區(qū)域和各沉陷位置在2018－09－12—2019－03－11的演化過程。 可以發(fā)現(xiàn)：在SAR 影像監(jiān)測期間濟(jì)寧城市規(guī)劃區(qū)范圍內(nèi)產(chǎn)生了10 余處明顯的沉陷盆地；通過疊加相應(yīng)時間段內(nèi)各礦井地下開采情況（圖3），確定所監(jiān)測的地表變形區(qū)域與開采工作面具有很好的一致性。 驗(yàn)證了SBAS 技術(shù)在監(jiān)測地下煤層開采所引起地表形變的可行性，其監(jiān)測結(jié)果可以有效提高采煤沉陷地調(diào)查的針對性。 通過對圖3 和圖4 進(jìn)一步分析可以得到以下結(jié)論：

1）地表形變速率超過100 mm／a 的位置主要集中在存在地下開采活動的各礦井，如：唐口煤礦、新河煤礦、濟(jì)寧二號煤礦、何崗煤礦、許廠煤礦、岱莊煤礦、運(yùn)河煤礦等。

2）各礦井范圍內(nèi)開采位置以外，仍監(jiān)測到大面積形變速率大于20 mm／a 的區(qū)域（主要為唐口煤礦和濟(jì)寧城市建成區(qū)）；對于唐口煤礦，監(jiān)測所得的地表形變速率約為20 mm／a，其原因主要是該區(qū)域?yàn)樯畈块_采，前期開采影響即老采空區(qū)的殘余變形仍較為顯著；對于濟(jì)寧城市建成區(qū)，監(jiān)測所得的地表形變速率約為15 mm／a，因該區(qū)域不存在地下煤炭資源的開采活動，分析地表變形原因除形變解譯誤差外，主要與城市建設(shè)和地下水開采有關(guān)。

3）研究區(qū)域以外的其他位置，地表沉陷速率大都低于10 mm／a，表明相應(yīng)位置的老采空區(qū)已基本趨于穩(wěn)定。

4.2 精度評價

Sentinel－1 數(shù)據(jù)SBAS 時序處理結(jié)果的精度是通過與相近時間段內(nèi)的水準(zhǔn)測量結(jié)果的比較進(jìn)行評價。 項(xiàng)目組為評價SBAS 解譯結(jié)果的可靠性，在研究區(qū)內(nèi)進(jìn)行了水準(zhǔn)測量工作。 本次精度評定，以唐口煤礦和岱莊煤礦的3 條觀測線（TKL4 測線、TKL9測線、DZL15 測線）的95 個觀測點(diǎn)的外業(yè)水準(zhǔn)結(jié)果作為參考，以10m 范圍內(nèi)的高相干點(diǎn)的平均值作為InSAR 形變監(jiān)測結(jié)果。 其中，唐口煤礦TKL4 測線、TKL9 測線和岱莊煤礦DZL15 測線水準(zhǔn)實(shí)測時間段分別為2019－01－06—2019－03－08、2018－11－25—2019－01－07、2018－12－26—2019－03－01；InSAR 監(jiān)測時間段分別為2019－01－10—2019－03－11、2018－11－23—2019－01－10、2018－12－29—2019－03－11。

圖5 和表2 為各測線的InSAR 解譯與水準(zhǔn)測量的對比結(jié)果與統(tǒng)計情況。

圖5 各時間段地表垂直沉陷圖Fig.5 Surface subsidence curve in different periods

表2 SBAS 地表形變解譯精度統(tǒng)計Table 2 Statistical surface deformation accuracy

從圖5 中各測線SBAS 解算結(jié)果與水準(zhǔn)測量成果的對比情況可以看出：當(dāng)?shù)乇硇巫冊贗nSAR 可探測的量級范圍內(nèi)時，兩者在沉陷的趨勢上基本保持一致，具有較高的相關(guān)性。 表2 中的統(tǒng)計數(shù)據(jù)顯示，TKL4、TKL9 和DZL15 三條測線上兩者最大差值分別為12.1、3.0 和11.6 mm，平均差值小于3.3 mm；監(jiān)測精度高，但所能監(jiān)測的最大沉陷量有限。

綜合以上分析可以得到，SBAS 技術(shù)所監(jiān)測的研究區(qū)沉陷盆地的狀態(tài)、位置和沉陷量，具有良好的準(zhǔn)確性和可靠性，監(jiān)測結(jié)果有利于指導(dǎo)礦山地質(zhì)環(huán)境調(diào)查工作高效、有針對性地開展。

4.3 老采空區(qū)形變監(jiān)測可行性分析

前述SBAS 地表形變的監(jiān)測結(jié)果和精度分析表明，InSAR 時序處理方法可以獲得研究區(qū)地表形變的時空演變過程，且其形變監(jiān)測精度能夠達(dá)到毫米級。 因此，該技術(shù)完全可以應(yīng)用于老采空區(qū)地表殘余形變的監(jiān)測，而且對于變形緩慢的老采空區(qū)地表，SBAS 應(yīng)更具有優(yōu)越性。

另一方面，為驗(yàn)證SBAS 技術(shù)鑒別不同停采時間和不同煤層采厚老采空區(qū)的能力，做了進(jìn)一步分析。 以唐口煤礦為例，其中圖6 為不同停采時間工作面上方地表點(diǎn)的時間演變曲線，圖7 為不同煤層采厚條件下地表點(diǎn)的時間演變曲線。

從圖6 和圖7 中可以得出：①老采空區(qū)上方不同位置的地面點(diǎn)，在SAR 影像監(jiān)測期間，仍然存在著持續(xù)的地表下沉現(xiàn)象。 ②隨著工作面停采時間的增長，地表下沉量逐漸減小；其中，停采時間超過2、5、10 a 的采空區(qū)上方在180 d 內(nèi)發(fā)生的地表殘余形變量分別小于20、15 和10 mm；滿足監(jiān)測點(diǎn)連續(xù)6個月累積下沉值小于30 mm 的指標(biāo)［20］，表明停采時間超過2 年的采空區(qū)上方地表已趨于穩(wěn)定。 ③老采空區(qū)上方地表的殘余變形量，隨煤層采厚的增加而增加，且兩者近似成比例關(guān)系。

上述結(jié)論表明，基于SBAS 的地表形變解譯結(jié)果，不僅能夠區(qū)分不同停采時間的采空區(qū)，對于不同采厚的采空區(qū)亦能有所區(qū)別。 同時，由于SBAS 技術(shù)毫米級的形變監(jiān)測精度，故其解譯結(jié)果可以作為老采空區(qū)穩(wěn)定性評價的依據(jù)（圖6、圖7 中圖注，例如：2008（3.50）含義為2008 年開采，開采厚度為3.50 m）。

圖6 不同停采年限地表點(diǎn)演變曲線Fig.6 Evolution curve of surface point in different stopping years

圖7 不同采厚條件下地表點(diǎn)演變曲線Fig.7 Evolution curve of surface point under different mining thickness conditions

5 結(jié) 論

1）利用14 景Sentinel－1 雷達(dá)影像，基于SBAS時序分析技術(shù)，成功探測到濟(jì)寧城市規(guī)劃區(qū)10 余處采煤塌陷地；監(jiān)測成果與水準(zhǔn)實(shí)測數(shù)據(jù)間的最大中誤差為5.8 mm。 研究區(qū)內(nèi)大部分范圍內(nèi)地表形變速率低于10 mm／a，表明相應(yīng)位置的老采空區(qū)已基本趨于穩(wěn)定。

2）SBAS 時序分析技術(shù)所提取老采空區(qū)地面點(diǎn)的變形曲線符合開采沉陷的基本規(guī)律；形變解譯結(jié)果能夠顯著區(qū)分不同開采條件下地表殘余變形的特征。

3）基于SBAS 時序分析技術(shù)的地表形變解譯成果能夠查明目標(biāo)區(qū)采煤沉陷地的空間分布位置，指導(dǎo)礦山地質(zhì)環(huán)境調(diào)查工作高效、有針對性地開展；可以為老采空區(qū)穩(wěn)定性的初步評價提供基礎(chǔ)依據(jù)。