李 柱， 范洪冬， 高彥濤， 許耀宗

(1.中國(guó)礦業(yè)大學(xué)礦山生態(tài)修復(fù)教育部工程研究中心，徐州 221116； 2.河南省地質(zhì)礦產(chǎn)勘查開(kāi)發(fā)局測(cè)繪地理信息院，鄭州 450006)

0 引言

煤火一般是指地下煤層火災(zāi)、近地表煤層火災(zāi)、煤矸石火災(zāi)和煤堆火災(zāi)[1]。大多數(shù)煤火是由于不適當(dāng)?shù)牟傻V造成的，也有一些是因處理不當(dāng)?shù)牡V山廢棄物、閃電和森林火災(zāi)引起的[2]。作為一種全球性災(zāi)害，煤火對(duì)人類(lèi)健康、煤礦安全和基礎(chǔ)設(shè)施造成了嚴(yán)重的負(fù)面影響。煤火的燃燒不僅消耗了大量寶貴的煤炭資源，同時(shí)引發(fā)地面沉降、裂縫、裂隙等地質(zhì)災(zāi)害，釋放的大量有毒氣體(CO，NO2，SO2，H2S)和微量元素(Hg，F(xiàn)，As，Se)還會(huì)造成嚴(yán)重的環(huán)境問(wèn)題[3]。因此，全面、詳細(xì)以及準(zhǔn)確地監(jiān)測(cè)煤火，對(duì)煤火的演變趨勢(shì)、自燃規(guī)律和燃燒狀態(tài)進(jìn)行研究以及開(kāi)展煤火治理活動(dòng)具有重要的意義。

遙感方法具有獲取信息速度快、監(jiān)測(cè)區(qū)域面積大、周期短、全天時(shí)、全天候和經(jīng)濟(jì)效益高的特點(diǎn)，這些獨(dú)特的優(yōu)勢(shì)可以彌補(bǔ)傳統(tǒng)煤火探測(cè)與監(jiān)測(cè)技術(shù)的不足[4]。利用遙感方法進(jìn)行地面沉降監(jiān)測(cè)和溫度反演，已經(jīng)成為煤火探測(cè)與監(jiān)測(cè)領(lǐng)域的一個(gè)重要方向。國(guó)內(nèi)外學(xué)者利用熱紅外遙感技術(shù)對(duì)烏達(dá)煤田煤火區(qū)進(jìn)行了大量研究[5-8]，可以較高精度地探測(cè)出煤火區(qū)及其演變趨勢(shì)，但熱紅外遙感技術(shù)很難實(shí)現(xiàn)對(duì)深部煤火、煤火燃燒階段以及煤火引起的地表形變的監(jiān)測(cè)，這使得熱紅外遙感技術(shù)監(jiān)測(cè)煤火時(shí)存在一定的不足。

差分合成孔徑雷達(dá)干涉測(cè)量(differential interferometric synthetic aperture Radar，D-InSAR)技術(shù)作為一種先進(jìn)的對(duì)地觀測(cè)手段已在各種變形監(jiān)測(cè)領(lǐng)域得到了廣泛的應(yīng)用，但是該技術(shù)易受時(shí)空失相關(guān)及大氣延遲等因素影響，難以在長(zhǎng)時(shí)序上得到高精度的地表形變結(jié)果[9]，而近年來(lái)發(fā)展的分布式目標(biāo)合成孔徑雷達(dá)干涉測(cè)量(distributed scatterer interfero-metric synthetic aperture Radar，DS-InSAR)方法則有效地解決了長(zhǎng)時(shí)序上低相干地區(qū)監(jiān)測(cè)點(diǎn)位密度不足的問(wèn)題[10]。目前，國(guó)內(nèi)外一些學(xué)者利用合成孔徑雷達(dá)干涉測(cè)量(interferometric synthetic aperture Radar，InSAR)技術(shù)對(duì)煤火區(qū)的探測(cè)與監(jiān)測(cè)進(jìn)行了相關(guān)研究，大多數(shù)研究側(cè)重于將InSAR監(jiān)測(cè)的地表形變結(jié)果作為輔助信息與地表熱異常信息相結(jié)合進(jìn)行煤火區(qū)域的探測(cè)[11-14]，而對(duì)煤火引起的地表形變規(guī)律研究較少[15-16]； 另外，對(duì)于烏達(dá)煤田煤火的研究多數(shù)是利用熱紅外遙感技術(shù)進(jìn)行探測(cè)與監(jiān)測(cè)，只有少數(shù)研究采用InSAR技術(shù)[17-19]，并且缺乏近些年的地表形變監(jiān)測(cè)結(jié)果。

針對(duì)上述問(wèn)題以及煤田火區(qū)的地表形變復(fù)雜、失相干嚴(yán)重等現(xiàn)象，本文研究了一種基于DS-InSAR技術(shù)的煤田火區(qū)監(jiān)測(cè)方法，利用該方法對(duì)烏達(dá)煤田火區(qū)2017年3月—2019年4月的63景Sentinel-1A影像進(jìn)行了處理，通過(guò)監(jiān)測(cè)結(jié)果對(duì)烏達(dá)煤田及其煤火區(qū)的地表形變進(jìn)行了分析，驗(yàn)證了該方法的適用性。

1 研究區(qū)域概況及數(shù)據(jù)源

1.1 研究區(qū)域概況

烏達(dá)煤田位于中國(guó)內(nèi)蒙古自治區(qū)西南部的烏海市烏達(dá)區(qū)，如圖1所示。該區(qū)東鄰黃河與鄂爾多斯黃土高原，南部是烏達(dá)工業(yè)園區(qū)，西鄰賀蘭山脈北部邊緣，北鄰烏蘭布和沙漠。煤田呈南北向斜狀，面積約為55 km2，包括五虎山、黃白茨和蘇海圖3個(gè)煤礦，海拔從1 100～1 300 m不等，平均高程為1 200 m。氣候?yàn)榈湫偷臏貛Т箨懶詺夂?，年降水量?68 mm，蒸發(fā)量為3 500 mm左右，常年干旱少雨。

圖1 研究區(qū)域Fig.1 The study area

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

Sentinel-1A衛(wèi)星由歐州航天局于2014年4月3日發(fā)射，采用近極地太陽(yáng)同步軌道，軌道高度為693 km，重訪周期為12 d，搭載了一個(gè)C波段的合成孔徑雷達(dá)(synthetic aperture Radar，SAR)，能提供超寬幅、干涉寬幅、條帶以及波浪4種成像模式。本文選取2017年3月—2019年4月期間的63景干涉寬幅模式成像下的SAR影像為數(shù)據(jù)源，其距離向和方位向的像元尺寸分別為2.3 m和14.0 m。為去除干涉圖中的地形相位，借助空間分辨率為90 m的SRTM DEM數(shù)據(jù)作為外部數(shù)字高程模型(digital elevation model，DEM)進(jìn)行差分干涉處理。

2 DS-InSAR方法原理

Ferretti等[20]于2011年提出了第二代永久散射體技術(shù)SqueeSAR。拉開(kāi)了DS-InSAR方法研究的序幕。與傳統(tǒng)的時(shí)序InSAR方法不同的是，DS-InSAR方法在其基礎(chǔ)上增加了同質(zhì)點(diǎn)選取和相位優(yōu)化2個(gè)步驟，有效地解決了傳統(tǒng)時(shí)序方法相干點(diǎn)選取數(shù)目少和空間分布不均勻的問(wèn)題。

2.1 同質(zhì)點(diǎn)選取

(1)

在勻質(zhì)區(qū)域的單視復(fù)數(shù)SAR影像服從瑞利分布，其變異系數(shù)CV為常數(shù)，即

(2)

假設(shè)像元P后向散射性穩(wěn)定，則式(1)可以轉(zhuǎn)換為：

(3)

2.2 特征值分解相位優(yōu)化

(4)

式中: 向量y=[y1y2,…,yN]為分布式目標(biāo)的同質(zhì)點(diǎn)在N景SAR影像上的復(fù)數(shù)觀測(cè)向量x=[x1,x2,…,xN]經(jīng)過(guò)振幅歸一化處理得到；NSHPs為分布式目標(biāo)同質(zhì)點(diǎn)的數(shù)目； H為Hermitian共軛轉(zhuǎn)置。

(5)

DS-InSAR數(shù)據(jù)處理具體流程如下： ①根據(jù)設(shè)置的時(shí)空基線閾值選擇干涉對(duì)，并生成干涉圖； ②利用FaSHPS方法對(duì)強(qiáng)度影像進(jìn)行同質(zhì)點(diǎn)識(shí)別，并通過(guò)特征值分解法對(duì)干涉圖進(jìn)行相位優(yōu)化； ③利用外部DEM去除干涉圖中的地形相位生成差分干涉圖，并利用優(yōu)化后的相位計(jì)算同質(zhì)點(diǎn)的時(shí)空相干性，設(shè)定相干性閾值選取分布式目標(biāo)(distributed scatterer，DS)點(diǎn)； ④利用差分干涉圖對(duì)DS點(diǎn)上的相位進(jìn)行解纏，并建立解纏相位與地表形變速率、DEM誤差和殘余相位之間的觀測(cè)方程； ⑤利用奇異值分解法估算地表形變相位和DEM殘余相位，并對(duì)殘余相位進(jìn)行時(shí)空濾波分離出非線性形變相位、大氣相位和噪聲相位； ⑥將線性形變相位與非線性形變相位相加解算研究區(qū)域視線向(line of sight，LOS)的地表時(shí)序形變。數(shù)據(jù)處理流程如圖2所示。

圖2 數(shù)據(jù)處理流程Fig.2 Data processing flow chart

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

3.1 監(jiān)測(cè)結(jié)果對(duì)比

分別采用TCP-InSAR與DS-InSAR方法對(duì)覆蓋烏達(dá)煤田的63景Sentinel-1A數(shù)據(jù)進(jìn)行處理，獲得了2017年3月—2019年4月時(shí)間段內(nèi)的地表時(shí)序形變信息，年平均形變速率如圖3所示。從圖3可以看出，二者雖在空間上形變趨勢(shì)具有較好的一致性，但DS-InSAR監(jiān)測(cè)結(jié)果能夠更好地反映出煤田的地表形變情況，并對(duì)研究區(qū)域形變明顯的8個(gè)區(qū)域進(jìn)行A—H編號(hào)。

(a) TCP-InSAR (b) DS-InSAR

TCP-InSAR方法共選出了156 853個(gè)監(jiān)測(cè)點(diǎn)； DS-InSAR方法共選出了194 585個(gè)DS點(diǎn)。從選取的監(jiān)測(cè)點(diǎn)數(shù)量而言，DS-InSAR能夠選出更多的DS點(diǎn)，監(jiān)測(cè)點(diǎn)位密度比TCP-InSAR提高1.24倍，并且在低相干區(qū)域也可以獲得足夠數(shù)量的監(jiān)測(cè)點(diǎn)； 從監(jiān)測(cè)點(diǎn)位的空間分布來(lái)看，DS-InSAR方法選出的DS點(diǎn)在保證點(diǎn)密度的同時(shí)分布更均勻，在一定程度上提高了形變監(jiān)測(cè)結(jié)果的解算精度，為研究區(qū)提供了更為詳細(xì)準(zhǔn)確的形變信息，有效地彌補(bǔ)了在低相干區(qū)域TCP-InSAR等傳統(tǒng)時(shí)序形變監(jiān)測(cè)方法中存在的缺陷。

為驗(yàn)證DS-InSAR監(jiān)測(cè)結(jié)果的可靠性，通過(guò)對(duì)比分析TCP-InSAR和DS-InSAR獲得的形變速率值來(lái)進(jìn)行交叉驗(yàn)證。本文以TCP點(diǎn)為基準(zhǔn)，選取最鄰近的DS點(diǎn)作為同名點(diǎn)，共篩選出95 760對(duì)同名點(diǎn)，通過(guò)同名點(diǎn)對(duì)的形變速率值及其差異的統(tǒng)計(jì)分析，繪制了形變速率相關(guān)性和差異分布直方圖，如圖4所示。通過(guò)形變速率相關(guān)性(圖4(a))可知，DS-InSAR與TCP-InSAR方法獲得的地表形變速率之間具有較強(qiáng)的相關(guān)性，利用Pearson計(jì)算的相關(guān)系數(shù)R為0.84。由形變速率差異分布直方圖(圖4(b))可知，2種方法獲得的形變速率差異較小，兩者之間的標(biāo)準(zhǔn)差為5.83 mm/a。因此，在缺少實(shí)測(cè)水準(zhǔn)數(shù)據(jù)的情況下，與TCP-InSAR方法監(jiān)測(cè)結(jié)果的交叉驗(yàn)證，表明了DS-InSAR地表形變監(jiān)測(cè)結(jié)果是較為可靠的。

3.2 烏達(dá)煤田地表形變分析

為了能夠全面分析烏達(dá)煤田地表形變特點(diǎn)，本文利用DS-InSAR對(duì)烏達(dá)煤田進(jìn)行了長(zhǎng)時(shí)序的地表形變監(jiān)測(cè)，時(shí)序形變累積如圖5所示。圖5中的14幅圖是每隔一段時(shí)間相對(duì)于起始時(shí)間2017年3月17日的累積形變量。從圖5中可以看出，在監(jiān)測(cè)時(shí)間段內(nèi)烏達(dá)煤田存在嚴(yán)重的地表形變現(xiàn)象。在空間上，沉降區(qū)域分布不均勻。在時(shí)間序列上，隨著時(shí)間的推移研究區(qū)域內(nèi)逐漸出現(xiàn)了形變量級(jí)與沉降范圍大小不同的多個(gè)沉降區(qū)域。其中，A，B，C這3個(gè)區(qū)域的形變現(xiàn)象最為嚴(yán)重，沉降范圍最廣，這3個(gè)區(qū)域的地表形變?cè)谡麄€(gè)監(jiān)測(cè)時(shí)間段內(nèi)持續(xù)存在，地表形變量與沉降范圍越來(lái)越大，導(dǎo)致B和C區(qū)域從2018年8月15日之后逐漸連接形成一個(gè)大的沉降區(qū)域； G區(qū)域的地表形變是從2017年5月16日開(kāi)始產(chǎn)生，前期形變變化較小，到后期形變變化加快； D和F區(qū)域從2017年8月20日之后地表開(kāi)始產(chǎn)生形變，并逐漸形成多個(gè)沉降中心； H區(qū)域的地表形變從2017年12月18日開(kāi)始發(fā)生地表形變； E區(qū)域產(chǎn)生地表形變的時(shí)間最晚，但地表形變變化較快，從2018年10月14日—2019年4月12日半年的時(shí)間里產(chǎn)生了與D區(qū)域相同量級(jí)的形變。為更詳細(xì)具體地分析地表形變情況，利用圖3(b)地表形變速率圖對(duì)其進(jìn)行定量分析。從圖5中可知，地表形變最為嚴(yán)重的A，B，C這3個(gè)區(qū)域的最大形變速率分別為-178.4 mm/a，-166.3 mm/a和-215 mm/a，在這3個(gè)沉降區(qū)域的沉降中心，由于形變梯度大，地表形變量在監(jiān)測(cè)時(shí)間段內(nèi)超過(guò)時(shí)序InSAR的監(jiān)測(cè)能力，失相干現(xiàn)象嚴(yán)重，無(wú)法得到足夠密度的DS點(diǎn)； 另外，F(xiàn)和G這2個(gè)區(qū)域同樣存在較為嚴(yán)重的形變，最大形變速率分別為-96.5 mm/a和-100.9 mm/a； 相比之下，D，E和H區(qū)域的形變量相對(duì)較小，最大形變速率分別為-46.5 mm/a，-53.9 mm/a和-48.5 mm/a。

圖5 烏達(dá)煤田地表形變時(shí)序累積圖

3.3 煤田火區(qū)地表形變分析

烏達(dá)煤田嚴(yán)重的地表形變現(xiàn)象主要與煤火燃燒及人工采挖有關(guān)[17,19]，根據(jù)文獻(xiàn)[7]中2018年1月底實(shí)地調(diào)查獲取的部分區(qū)域的明顯煙點(diǎn)和煤火點(diǎn)，這些實(shí)測(cè)的煤火點(diǎn)附近地區(qū)伴隨著地表形變現(xiàn)象。為了研究煤火燃燒引起的地表形變，對(duì)實(shí)測(cè)煤火區(qū)的地表形變進(jìn)行重點(diǎn)分析。

3.3.1 煤田火區(qū)地表形變時(shí)序分析

為分析煤田火區(qū)的時(shí)序地表形變規(guī)律，將實(shí)測(cè)煤火點(diǎn)編號(hào)為H1—H6(圖3(b))，對(duì)應(yīng)D，A和G這3個(gè)沉降區(qū)域。選取3個(gè)區(qū)域中實(shí)測(cè)煤火點(diǎn)附近且累積形變量最大的5個(gè)DS點(diǎn)作為特征點(diǎn)分析地表形變，編號(hào)記為P1—P5(圖3(b))。為了驗(yàn)證特征點(diǎn)時(shí)序形變的可靠性，從TCP-InSAR監(jiān)測(cè)結(jié)果中選取距離特征點(diǎn)最近的TCP點(diǎn)作為驗(yàn)證點(diǎn)對(duì)其進(jìn)行驗(yàn)證，由于TCP-InSAR選點(diǎn)的原因，只在P3和P5點(diǎn)旁選了P3′和P5′點(diǎn)。圖6(a)—(c)分別為D，A和G這3個(gè)沉降區(qū)域的特征點(diǎn)及驗(yàn)證點(diǎn)時(shí)序形變圖，特征點(diǎn)與驗(yàn)證點(diǎn)之間的時(shí)序形變趨勢(shì)基本吻合，具有很好的一致性，2點(diǎn)間的時(shí)序形變值整體存在偏大或偏小的問(wèn)題，但考慮到2點(diǎn)之間具有一定的距離，屬于合理范圍，所以可以利用選取的特征點(diǎn)進(jìn)行地表時(shí)序形變分析。

由圖6可知，3個(gè)沉降區(qū)域的形變量級(jí)與時(shí)序形變趨勢(shì)各不相同，不同時(shí)間段地表形變變化快慢也不相同。P1從2018年5月底才開(kāi)始發(fā)生形變并持續(xù)到監(jiān)測(cè)時(shí)間結(jié)束，這段時(shí)間內(nèi)地表形變變化快，煤火燃燒劇烈。處于同一沉降區(qū)域的P2，P3和P4開(kāi)始發(fā)生形變的時(shí)間各不相同，P2點(diǎn)最晚，從2017年9月開(kāi)始，并且中間一段時(shí)間因煤火處于休眠階段未發(fā)生明顯的地表形變現(xiàn)象； P3點(diǎn)地表形變變化緩慢，但在整個(gè)監(jiān)測(cè)時(shí)間段內(nèi)形變持續(xù)發(fā)生； 相比之下，P4點(diǎn)從2017年6月開(kāi)始出現(xiàn)形變現(xiàn)象，且在之后時(shí)間段內(nèi)出現(xiàn)了3個(gè)未發(fā)生形變的時(shí)期。P5點(diǎn)的時(shí)序形變趨勢(shì)更加復(fù)雜，多個(gè)時(shí)間段內(nèi)出現(xiàn)形變停滯現(xiàn)象，但在發(fā)生形變的時(shí)期內(nèi)形變變化較快。煤火區(qū)的這些復(fù)雜地表時(shí)序形變現(xiàn)象歸因于煤火燃燒階段的復(fù)雜性、燃燒強(qiáng)度的非均質(zhì)性，并且可以從這些特征點(diǎn)時(shí)序上的地表形變變化快慢推斷煤火的燃燒強(qiáng)度及燃燒階段。從整體上分析這5個(gè)特征點(diǎn)時(shí)序地表形變規(guī)律，可以發(fā)現(xiàn)秋冬季節(jié)的地表形變變化相對(duì)較快，這是由于烏達(dá)煤田處于干旱半干旱地區(qū)，秋冬季節(jié)干旱少雨，煤層易發(fā)生火災(zāi)[23]，導(dǎo)致地表形變變化較快。

3.3.2 典型煤火區(qū)地表形變特征分析

為分析煤田火區(qū)的地表形變特征，選取3個(gè)沉降區(qū)域中形變最嚴(yán)重的A區(qū)作為典型煤火區(qū)。為清晰直觀地表征區(qū)域A地表形變的空間形態(tài)，將地表累積形變圖(圖7(a))進(jìn)行克里金插值，并將插值之后的圖轉(zhuǎn)換成三維形變圖(圖7(b))。從圖7(b)可知，區(qū)域A的地表形變?cè)诳臻g上變化較快，整體形變規(guī)律復(fù)雜，邊緣形狀不規(guī)則，且具有多個(gè)發(fā)育程度不同的沉降中心。從三維形變圖的形變等值線中可知，這些沉降中心附近的等值線密集，向四周逐漸變得稀疏，對(duì)應(yīng)的沉降中心凹陷較深，向四周逐漸平緩，邊緣近似長(zhǎng)橢圓形，并且這些沉降中心的形變延伸方向隨著煤火蔓延和燃燒強(qiáng)度而各不相同。

(a) 剖面線示意圖 (b) 三維形變圖

為定量分析典型煤火區(qū)的地表形變特征，在圖7(a)區(qū)域A的沉降中心附近給出了2條剖面線a1和a2，并沿剖面線提取地表累積形變量，結(jié)果如圖8所示。由圖8可知，剖面線a1穿過(guò)了3個(gè)形變量級(jí)與沉降范圍大小不同的沉降中心，最嚴(yán)重的沉降中心累積形變量達(dá)到-311.6 mm； 剖面線a2經(jīng)過(guò)2個(gè)發(fā)育程度相似的沉降中心，最大累積形變量從左到右分別為-369.5 mm和-355.2 mm； 這些沉降中心到邊緣的形變量變化快，形變梯度較大。

圖8 典型煤火區(qū)剖面累積形變圖

4 結(jié)論

本研究利用DS-InSAR對(duì)烏達(dá)煤田火區(qū)2017年3月—2019年4月采集的63幅Sentinel-1A影像進(jìn)行了處理，分析了長(zhǎng)時(shí)序地表形變，得到如下結(jié)論：

1)與TCP-InSAR技術(shù)相比，DS-InSAR監(jiān)測(cè)點(diǎn)密度提高1.24倍，且分布均勻，能夠更好地反映出研究區(qū)域的地表形變情況。研究區(qū)域存在多個(gè)形變量級(jí)與沉降范圍大小不同的沉降區(qū)域，且在空間上分布不均勻，最大地表形變速率約為-215 mm/a。

2)監(jiān)測(cè)時(shí)間段內(nèi)烏達(dá)煤田存在因煤火燃燒導(dǎo)致的地表形變。煤火區(qū)的地表時(shí)序形變規(guī)律復(fù)雜，秋冬季節(jié)地表形變變化相對(duì)較快，并且利用地表時(shí)序形變可以推斷出煤火的燃燒狀態(tài)，形變嚴(yán)重的煤火區(qū)出現(xiàn)多個(gè)形變延伸方向且發(fā)育程度不同的沉降中心。

3)本研究?jī)H利用DS-InSAR獲取的結(jié)果對(duì)烏達(dá)煤田及其煤火區(qū)的地表形變進(jìn)行了分析，將來(lái)還需要與現(xiàn)場(chǎng)實(shí)測(cè)資料以及開(kāi)采沉陷理論知識(shí)相結(jié)合，進(jìn)一步研究煤火燃燒導(dǎo)致的地表形變機(jī)理，最終實(shí)現(xiàn)對(duì)煤田火區(qū)地表形變?nèi)嬖敿?xì)準(zhǔn)確的監(jiān)測(cè)。