基于DinSAR技術(shù)的寧東煤炭基地區(qū)域性地表沉降監(jiān)測(cè)

程霞 ，王輝 ，王國瑞 ，杜靈通 ，魏采用 ，周峰

（1. 寧夏回族自治區(qū)國土資源調(diào)查監(jiān)測(cè)院，寧夏 銀川 750002；2. 寧夏大學(xué)西北土地退化與生態(tài)恢復(fù)省部共建國家重點(diǎn)實(shí)驗(yàn)室培育基地，寧夏 銀川 750021；3. 寧夏大學(xué)生態(tài)環(huán)境學(xué)院，寧夏 銀川 750021；4. 寧夏回族自治區(qū)遙感調(diào)查院，寧夏 銀川 750021；5. 寧夏回族自治區(qū)礦產(chǎn)地質(zhì)調(diào)查院，寧夏 銀川 750021）

近年來，干旱-半干旱區(qū)煤炭資源大規(guī)模開采所引起的區(qū)域地質(zhì)環(huán)境惡化被廣泛關(guān)注（董英等，2008；杜臻等，2023），寧東煤炭基地開發(fā)對(duì)區(qū)域植被生態(tài)的影響也引起了一些學(xué)者注意，他們利用不同手段從不同方面開展評(píng)估（王景升等，2013；杜靈通等，2018），其中煤炭開采導(dǎo)致的地表沉降塌陷一直備受關(guān)注（謝曉深等，2018）。王國瑞等（2018）利用地面GPS站對(duì)寧東金鳳煤礦的地面變形進(jìn)行監(jiān)測(cè)，結(jié)果表明GPS站可以精準(zhǔn)的監(jiān)測(cè)地面點(diǎn)的沉降過程，特別是能捕捉到突然沉降的時(shí)間點(diǎn)；高冠杰等（2018）利用無人機(jī)光學(xué)遙感技術(shù)獲取的多期地面高程數(shù)據(jù)對(duì)寧東羊場(chǎng)灣煤礦的地表沉降量進(jìn)行監(jiān)測(cè)，結(jié)果表明無人機(jī)光學(xué)遙感技術(shù)不僅能在小區(qū)域范圍內(nèi)對(duì)重點(diǎn)采礦區(qū)的地表沉降進(jìn)行監(jiān)測(cè)，還能解譯提取地表變形引起的裂縫信息；張永庭等（2018a）利用無人機(jī)激光雷達(dá)對(duì)寧東馬蓮臺(tái)煤礦重點(diǎn)沉降區(qū)開展監(jiān)測(cè)，得出無人機(jī)激光雷達(dá)技術(shù)能很好地獲取地表三維沉降數(shù)據(jù)，但受制于無人機(jī)飛行范圍的限制。綜上所述，目前的研究主要從單點(diǎn)開展寧東地表沉降監(jiān)測(cè)，或在小范圍開展沉降監(jiān)測(cè)實(shí)驗(yàn)，而對(duì)于整個(gè)寧東規(guī)劃區(qū)的地表沉降監(jiān)測(cè)工作尚未開展，這導(dǎo)致政府部門無法對(duì)區(qū)域性地質(zhì)環(huán)境變化情況進(jìn)行準(zhǔn)確評(píng)估。

隨著合成孔徑雷達(dá)（Synthetic Aperture Radar，SAR）技術(shù)的發(fā)展，特別是合成孔徑雷達(dá)差分干涉測(cè)量技術(shù)（Differential InSAR，DInSAR）的出現(xiàn)，推動(dòng)了利用衛(wèi)星遙感數(shù)據(jù)開展大范圍地表沉降監(jiān)測(cè)的應(yīng)用（李廣宇等，2018；馮旻譞等，2023），為大型煤炭基地地表沉降監(jiān)測(cè)提供了數(shù)據(jù)和技術(shù)手段。為此，筆者利用歐空局哨兵一號(hào)衛(wèi)星（Sentinel-1A）獲取的覆蓋寧東煤炭基地規(guī)劃區(qū)數(shù)據(jù)，監(jiān)測(cè)2015～2019年該區(qū)域地表的沉降特征，以期為準(zhǔn)確評(píng)估煤炭開發(fā)對(duì)該區(qū)域的地質(zhì)環(huán)境影響提供數(shù)據(jù)支撐和理論依據(jù)。

1 研究區(qū)概況

寧東煤炭基地位于陜、甘、寧、蒙毗鄰地區(qū)，地理坐標(biāo)E 106°21′39″～106°56′34″、N 37°04′48″～38°17′41″，東起鴛鴦湖、馬家灘、萌城礦區(qū)的深部邊界，西至白芨灘東界，南起韋州礦區(qū)和萌城礦區(qū)南端的寧夏與甘肅省界，北至寧夏與內(nèi)蒙古省界，東西寬約為16～41 km，南北長約為127 km，規(guī)劃總面積為3 484 km2，是國家14個(gè)億噸級(jí)煤炭基地之一。行政區(qū)劃覆蓋靈武市、鹽池縣、同心縣、紅寺堡開發(fā)區(qū)等4個(gè)縣市（區(qū)），包括靈武市臨河鎮(zhèn)、寧東鎮(zhèn)、馬家灘鎮(zhèn)、白土崗鄉(xiāng)，鹽池縣惠安堡鎮(zhèn)、馮記溝鄉(xiāng)，同心縣韋州鎮(zhèn)、下馬關(guān)鎮(zhèn)，紅寺堡太陽山鎮(zhèn)等9鄉(xiāng)鎮(zhèn)的部分區(qū)域（圖1）?；氐V區(qū)包括8個(gè)部分，位于北部的橫城和南部的韋州為兩個(gè)古生代石炭紀(jì)—二疊紀(jì)煤礦區(qū)，位于中部的自西向東為石溝驛、碎石井、鴛鴦湖、馬家灘、積家井、萌城6個(gè)中生代侏羅紀(jì)煤礦區(qū)，礦區(qū)井田總面積為1 591.21 km2，占基地總面積的45.67%，其中以馬家灘、積家井礦區(qū)面積最大。規(guī)劃核心區(qū)處于靈武境內(nèi)，核定面積為766 km2，約占規(guī)劃區(qū)總面積的22%。

2 數(shù)據(jù)與方法

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

2.1.1 Sentinel-1A干涉寬幅模式數(shù)據(jù)

2014年4月3日，歐空局發(fā)射了哨兵系列衛(wèi)星之一的Sentinel-1A，該衛(wèi)星使用近極地太陽同步軌道，軌道高度693 km，重訪周期為12天，Sentinel-1A搭載的C波段合成孔徑雷達(dá)，可以提供覆蓋全球任何地方的合成孔徑雷達(dá)（SAR）圖像（李廣宇等，2018）。Sentinel-1A可實(shí)現(xiàn)不同方式的極化和條帶模式、超寬幅模式、干涉寬幅模式及波模式等4種工作模式，其中，以干涉寬幅模式獲得的IW（Interferometric Wide Swath）影像掃描幅寬為250 km，距離向分辨率為5.0 m，方位向分辨率為20.0 m，可實(shí)現(xiàn)HH+HV、VV+VH、HH和VV 4種極化模式，能夠?qū)崿F(xiàn)大范圍覆蓋的干涉測(cè)量測(cè)量（歐陽倫曦等，2017）。筆者選取覆蓋寧東煤炭基地規(guī)劃區(qū)的5景Sentinel-1A IW斜距單視復(fù)數(shù)產(chǎn)品（SLC），下載網(wǎng)址為https://scihub.copernicus.eu/dhus，獲取時(shí)間分別為2015年5月20日、2016年5月14日、2017年5月21日、2018年5月16日和2019年5月23日，成像時(shí)間集中在5月中下旬，植被生長等季相一致，便于地表形變干涉測(cè)量分析。

2.1.2 精密軌道星歷數(shù)據(jù)

軌道信息是InSAR數(shù)據(jù)處理中非常重要的信息，從最初的圖像配準(zhǔn)到最后的形變圖像生成都有著重要的作用。Sentinel-1A衛(wèi)星采用精密軌道控制技術(shù)，其衛(wèi)星質(zhì)量控制子系統(tǒng)網(wǎng)站提供最精確的軌道數(shù)據(jù)，即POD Precise Orbit Ephemerides（POD精密軌道星歷數(shù)據(jù)），該數(shù)據(jù)距離GNSS下行21天之后才可提供（李廣宇等，2018），每天產(chǎn)生一個(gè)文件，每個(gè)文件覆蓋26個(gè)小時(shí)，定位精度優(yōu)于5 cm。本研究為提高形變測(cè)量的精度，使用Sentinel-1A衛(wèi)星精密軌道星歷數(shù)據(jù)對(duì)軌道信息進(jìn)行修正，有效去除因軌道誤差而引起的系統(tǒng)性誤差（Sentinel衛(wèi)星精密軌道星歷數(shù)據(jù)獲取網(wǎng)站https://qc.sentinel1.eo.esa.int/）。

2.1.3 數(shù)字地形高程模型（DEM）

筆者使用了SRTM的DEM數(shù)據(jù)，該數(shù)據(jù)由美國國家航空航天局（NASA）和國防部國家測(cè)繪局（NIMA）聯(lián)合測(cè)量，SRTM覆蓋地球80%以上的陸地表面。SRTM數(shù)據(jù)產(chǎn)品2003年開始公開發(fā)布，經(jīng)歷多次修訂，目前的數(shù)據(jù)修訂版本為V4.1版本。由于覆蓋寧東煤炭基地規(guī)劃區(qū)的Sentinel-lA數(shù)據(jù)跨約經(jīng)度范圍較大，因此下載了跨越2個(gè)經(jīng)度方格的SRTM_57_05和SRTM_58_05數(shù)據(jù)，并進(jìn)行了無縫拼接，數(shù)據(jù)來源于中國科學(xué)院計(jì)算機(jī)網(wǎng)絡(luò)信息中心地理空間數(shù)據(jù)云平臺(tái)（http://www.gscloud.cn）。

2.2 方法

差分干涉測(cè)量（DInSAR）是利用同一地區(qū)不同時(shí)相的SAR影像，通過差分干涉，獲取該地區(qū)地表微小形變信息的技術(shù)手段，常用的有兩軌法、三軌法和四軌法。DInSAR測(cè)量中要雷達(dá)波長、傳感器高度、天線基線距和波束視向等的幾何關(guān)系，開展解算，從而提取研究區(qū)域地面目標(biāo)的微小形變信息，其基本原理如下（鄭美楠等，2017）：

式中：λ為雷達(dá)波長；R1、R2為傳感器到地面目標(biāo)的距離；?f為地球曲面平地效應(yīng)引起的干涉相位；?t為地形引起的干涉相位；?d為地表形變引起的干涉相位；?n為相位噪聲；?a為大氣對(duì)流層及電離層延遲引起的相位。

SARscape是專業(yè)雷達(dá)圖像處理軟件，支持Sentinel-1A數(shù)據(jù)的導(dǎo)入和處理，其中包括DInSAR干涉測(cè)量處理。干涉測(cè)量之前，需先設(shè)置專門針對(duì)Sentinel-1A數(shù)據(jù)做InSAR處理時(shí)的系統(tǒng)參數(shù)，再導(dǎo)入數(shù)據(jù)（包括雷達(dá)數(shù)據(jù)、精密軌道數(shù)據(jù)和DEM數(shù)據(jù)），在計(jì)算干涉相對(duì)的時(shí)空基線的基礎(chǔ)上，開展DInSAR干涉測(cè)量處理（李旺，2018），詳細(xì)步驟流程如圖2所示。

圖2 基于SARscape的DinSAR測(cè)量技術(shù)流程圖Fig. 2 Flowchart of differential interferometry synthetic aperture radar by SARscape

3 結(jié)果與分析

3.1 相位干涉結(jié)果

選取2015年5月20日的影像作為主影像，對(duì)其余副影像相對(duì)于主影像的時(shí)空基線進(jìn)行計(jì)算（表1）。其中，空間基線最長為101.11 m，最短為16.95 m，均小于臨界基線，DInSAR的一個(gè)相位變化周期代表的地形變化為0.028 m。為了監(jiān)測(cè)年際間的地表沉降速率和沉降特征，遙感影像選擇按整年間隔，為2015～2019年每年5月中旬的Sentinel-1A數(shù)據(jù)，做干涉測(cè)量時(shí)，2幅主副影像的時(shí)間基線在360 d或372 d。

表1 Sentinel-1A干涉像對(duì)及其時(shí)空基線表Tab. 1 Temporal and spatial (perpendicular) baselines of Sentinel-1A interferograms

選擇VV極化方式，做差分干涉處理，并生成干涉圖。從2015年5月20日與2016年5月14日的干涉結(jié)果可以看出（圖3），寧東煤炭基地規(guī)劃區(qū)總體干涉信號(hào)較弱，未出現(xiàn)較強(qiáng)的干涉條紋。但去平后的干涉圖開始出現(xiàn)干涉條紋，特別是經(jīng)過濾波處理和重去平處理后，干涉條紋明顯增強(qiáng)，在寧東北部的馬蓮臺(tái)、羊場(chǎng)灣等煤礦和南部的韋州煤礦周圍出現(xiàn)較強(qiáng)的干涉條紋，即在這些區(qū)域出現(xiàn)了因采煤沉降地表形變所引起的雷達(dá)波相位干涉。由此可見，在寧東煤炭基地規(guī)劃區(qū)大范圍的沉降監(jiān)測(cè)中，需要對(duì)雷達(dá)影像進(jìn)行深度的增強(qiáng)處理，以精確提取相干信息。

圖3 不同處理下的干涉條紋圖（2015～2016年）Fig. 3 Interferogram with different processing (from 2015 to 2016)

3.2 地表年際沉降特征

利用2015～2019年的5期Sentinel-1A數(shù)據(jù)做干涉測(cè)量，監(jiān)測(cè)每年寧東煤炭基地規(guī)劃區(qū)范圍內(nèi)的地表沉降情況。從監(jiān)測(cè)結(jié)果可以看出（圖4），寧東煤炭基地規(guī)劃區(qū)內(nèi)，不同年份間的沉降特征存在較大差異。2015～2016年，北部地區(qū)整體沉降量較大，最大沉降量達(dá)到了0.13 m；2016～2017年，整體沉降量最小，其中南部地區(qū)無沉降，北部僅在羊場(chǎng)灣、任家莊、梅花井、清水營等煤礦一帶沉降較為明顯，最大沉降量為0.10 m；2017～2018年，沉降較明顯的地區(qū)也主要圍繞煤礦開采活動(dòng)區(qū)及其周邊，最大沉降量為0.14 m；2018～2019年，寧東煤炭基地區(qū)域性普遍發(fā)生地表沉降，但沉降量較大的依然在煤礦開采活動(dòng)區(qū)，最大沉降量為0.11 m。不同年份間的地表沉降差異，可能與煤礦開采過程有關(guān)，如王國瑞等（2018）監(jiān)測(cè)的金鳳煤礦在2015～2016年有一個(gè)快速的沉降過程，而2016年以后基本進(jìn)入沉降穩(wěn)定期，地表沉降開始變得緩慢。從沉降區(qū)域分布來看，寧東煤炭基地的地表沉降主要與寧東大型礦井開采活動(dòng)區(qū)范圍密切相關(guān)（圖4），特別是煤礦開采集中的北部地區(qū)，也是沉降最為明顯的區(qū)域。

圖4 不同時(shí)間段寧東地表沉降及主要煤礦開采區(qū)Fig. 4 Subsidence map of Ningdong region in different periods and mainly coal mining areas

3.3 地表累積沉降特征

將2015～2019年的沉降量累加，得到寧東煤炭基地區(qū)域近4年的總沉降量（圖5）。從圖中可以看出，總沉降量最大的區(qū)域可達(dá)到0.37 m沉降中心在北部的清水營、靈新、梅花井、羊場(chǎng)灣、石槽村、麥垛山煤礦一帶。為進(jìn)一步分析沉降中心和沉降狀況，選擇了橫穿寧東煤炭基地規(guī)劃區(qū)東西和南北的兩條剖面AB和CD，繪制沉降剖面圖。從AB剖面可以看出，橫穿羊場(chǎng)灣和梅花井礦區(qū)的剖面，在羊場(chǎng)灣和梅花井礦附近出現(xiàn)了多個(gè)超過0.1 m深的沉降中心，從剖面和沉降圖均可看出，這一區(qū)域是寧東煤炭基地地表沉降最嚴(yán)重的區(qū)域。從CD剖面來看，寧東煤炭基地中北部沉降較為嚴(yán)重，且沿著CD剖面出現(xiàn)了4個(gè)沉降中心，到南部以后，沉降量均小于0.05 m，沉降幅度較小。

圖5 2015～2019年累積沉降量及沉降切面圖Fig. 5 Cumulative subsidence and its profile of Ningdong region from 2015 to 2019

利用ARCGIS空間統(tǒng)計(jì)功能，對(duì)寧東煤炭基地規(guī)劃區(qū)的沉降面積和沉降量進(jìn)行統(tǒng)計(jì)計(jì)算，結(jié)果表明，寧東地區(qū)2015～2019年沉降量超過0.1 m的區(qū)域面積為578 km2，為主要的沉降區(qū)；沉降量在0.05～0.1 m的沉降區(qū)有2 484 km2，面積占比最大；沉降量小于0.05 m的微小沉降變形區(qū)面積為313 km2，還有109 km2的區(qū)域未發(fā)生地表沉降或發(fā)生地表抬升。結(jié)合寧東各礦分布和煤炭開采情況，對(duì)比沉降區(qū)分布特征和沉降量，可以得出沉降量超過0.1 m的沉降區(qū)主要為煤炭開采引發(fā)的地表沉降，其他微量地表沉降可能一方面來自雷達(dá)影像的相干失真，另一方面與地表其他風(fēng)沙侵蝕等變形有關(guān)。

3.4 監(jiān)測(cè)精度分析

為評(píng)價(jià)DinSAR監(jiān)測(cè)的寧東煤炭基地區(qū)域地表沉降的精度，本研究利用布設(shè)在寧東金鳳、靈新、梅花井、清水營、任家莊和羊場(chǎng)灣等煤礦的26個(gè)野外GPS沉降監(jiān)測(cè)點(diǎn)的數(shù)據(jù)，對(duì)DinSAR監(jiān)測(cè)沉降量的精度進(jìn)行驗(yàn)證，考慮到與遙感面上數(shù)據(jù)的匹配性，選擇的這些沉降監(jiān)測(cè)點(diǎn)均未在監(jiān)測(cè)期內(nèi)發(fā)生過短期突然沉降事件。首先根據(jù)GPS沉降監(jiān)測(cè)點(diǎn)的經(jīng)緯度坐標(biāo)，在2015～2019年累積沉降量遙感監(jiān)測(cè)圖像中提取對(duì)應(yīng)坐標(biāo)像元的數(shù)值，再將其與野外GPS沉降監(jiān)測(cè)點(diǎn)的對(duì)應(yīng)實(shí)測(cè)沉降數(shù)據(jù)進(jìn)行對(duì)比（圖6）。結(jié)果表明，Din-SAR提取的像元地表沉降與地面實(shí)測(cè)沉降有極顯著的相關(guān)性（R2=0.835 7，p＜0.01），即DinSAR監(jiān)測(cè)出的沉降強(qiáng)度與實(shí)測(cè)一致。但從圖6依然看出，實(shí)測(cè)沉降量要普遍大于DinSAR監(jiān)測(cè)的沉降量，這可能與二者的監(jiān)測(cè)機(jī)理不同有關(guān)，GPS監(jiān)測(cè)的是地表某一點(diǎn)的沉降，而DinSAR獲取的是對(duì)應(yīng)坐標(biāo)點(diǎn)所在區(qū)域幾百平方米區(qū)域內(nèi)的整體雷達(dá)反射，差分干涉測(cè)量出的結(jié)果也是一個(gè)區(qū)域面的沉降，這導(dǎo)致了驗(yàn)證點(diǎn)對(duì)應(yīng)數(shù)據(jù)之間的差異。

圖6 DinSAR監(jiān)測(cè)沉降量的精度驗(yàn)證圖Fig. 6 Precision validation of subsidence by DinSAR

筆者利用DinSAR技術(shù)和Sentinel-1A數(shù)據(jù)開展了寧東煤炭基地規(guī)劃區(qū)2015～2019年的沉降監(jiān)測(cè)。但衛(wèi)星合成孔徑雷達(dá)技術(shù)尚在發(fā)展之中，其中還存在一些技術(shù)難點(diǎn)和精度較低的問題，在地表形變監(jiān)測(cè)中還需進(jìn)一步完善，主要表現(xiàn)在以下2個(gè)方面：①時(shí)間失相干問題。由于兩幅影像成像時(shí)間間隔較長，地面覆蓋物雷達(dá)后向散射特性發(fā)生較大變化。本研究使用了近1 a的時(shí)間基線，一定程度上會(huì)造成兩幅影像的雷達(dá)散射信息不相干，導(dǎo)致干涉條紋不連續(xù)，信噪比降低。由于Sentinel-1A衛(wèi)星使用波長較短的C波段，對(duì)地物的穿透能力較弱，在春季和夏季植被茂密時(shí)期，如農(nóng)田、森林等覆蓋類型區(qū)（白澤朝等，2019），時(shí)間失相干會(huì)比較明顯。但寧東煤炭基地規(guī)劃區(qū)地處西北干旱-半干旱荒漠草原區(qū)，地表植被覆蓋度低，規(guī)劃區(qū)內(nèi)農(nóng)田、森林等土地利用類型的占比非常?。◤堄劳サ?，2018b），故在一定程度上會(huì)減弱時(shí)間失相干問題。②空間失相干。由于兩幅影像成像時(shí)側(cè)視角不同引起的失相干，即干涉對(duì)空間基線較長時(shí)，成像時(shí)刻雷達(dá)側(cè)視角變化大，回波信號(hào)的相干性信號(hào)較低，從而導(dǎo)致兩幅影像的失相干，往往發(fā)生在相對(duì)空間基線長度大于臨界基線情況下。筆者所用的5期影像，計(jì)算出的兩兩相對(duì)空間基線均在臨界基線以內(nèi)，故最大程度的降低了空間失相干現(xiàn)象。

4 結(jié)論

（1）寧東不同區(qū)域的地表沉降存在較大差異，其中北部地區(qū)整體沉降量較大，年最大沉降量達(dá)到了0.13 m，南部整體沉降量較小。

（2）由于煤炭開采階段不同，導(dǎo)致不同年份間的沉降量存在差異，2015～2016年沉降量較大，而2016～2017年沉降量較小。

（3）地表沉降區(qū)與煤礦開采活動(dòng)區(qū)的空間分布一致，2015～2019年沉降量超過0.1 m的區(qū)域面積為578 km2，主要為煤炭開采引發(fā)的地表沉降。

（4）微量地表沉降可能一方面來自雷達(dá)影像的相干失真，另一方面與地表風(fēng)沙侵蝕等變形有關(guān)。