張 豐，刁鑫鵬，譚秀全，李?lèi)?ài)軍，楊 雋

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

0 引 言

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

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

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

1 SBAS 基本原理

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

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

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

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

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

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

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

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

2.2 數(shù)據(jù)來(lái)源

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

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

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

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

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

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

4 監(jiān)測(cè)結(jié)果分析與精度評(píng)價(jià)

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

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

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

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

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

4.2 精度評(píng)價(jià)

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

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

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

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

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

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

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

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

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

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

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

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

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

5 結(jié) 論

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

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

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