麻學(xué)飛， 張雙成， 惠文華， 許 強(qiáng)

(1.長安大學(xué)地質(zhì)工程與測繪學(xué)院，西安 710054； 2.地理信息工程國家重點(diǎn)實(shí)驗室，西安 710054)

0 引言

煤炭作為我國的支柱能源，關(guān)系到國家經(jīng)濟(jì)的命脈，需要不斷地進(jìn)行開采，這勢必會形成采空區(qū)。采空區(qū)上部巖層發(fā)生斷裂、彎曲等形變，會造成采空區(qū)沉陷，破壞當(dāng)?shù)厣鷳B(tài)環(huán)境，威脅居民的人身財產(chǎn)安全[1-2]。因此對礦區(qū)進(jìn)行大范圍探測和動態(tài)監(jiān)測具有重要的研究意義。水準(zhǔn)測量、全球定位系統(tǒng)(global positioning system，GPS)等方式是傳統(tǒng)的礦區(qū)地表形變監(jiān)測手段，但由于這些方法監(jiān)測成本高、工作量大等原因，難以進(jìn)行長期監(jiān)測，且所獲監(jiān)測結(jié)果時空分辨率低，難以反映真實(shí)的礦區(qū)形變過程[3-5]。

合成孔徑雷達(dá)干涉測量(interferometry synthetic aperture Radar，InSAR)作為一種高精度的測量手段，利用相位差獲取研究區(qū)地表形變信息，具有監(jiān)測范圍大、效率高、分辨率高的特點(diǎn)，幾乎不受云霧的影響，成為了礦區(qū)地表監(jiān)測的有效手段[6-7]。 差分合成孔徑雷達(dá)干涉測量方法(differential interferometric synthetic aperture Rader,D-InSAR)可獲取到2景合成孔徑雷達(dá)(synthetic aperture Radar,SAR)影像覆蓋期間形變體表面的微小形變信息[7]，從而實(shí)現(xiàn)對形變體進(jìn)行監(jiān)測，處理數(shù)據(jù)量小，可以快速獲取目標(biāo)區(qū)域的形變信息。Ge等[8]利用D-InSAR技術(shù)對Appin，Westcliff和Tower煤礦區(qū)地表形變進(jìn)行監(jiān)測，所獲監(jiān)測結(jié)果精度可達(dá)厘米級； 龍四春等[9]利用D-InSAR與GPS集成實(shí)現(xiàn)了對資興唐煤礦區(qū)三維形變監(jiān)測，并對測量結(jié)果進(jìn)行集成與內(nèi)插，驗證了2種監(jiān)測結(jié)果的一致性。但是D-InSAR技術(shù)難以獲取目標(biāo)區(qū)域時序變化過程，且在實(shí)際應(yīng)用過程中常常受到時空去相干、大氣延遲等影響[10-11]。而時序InSAR彌補(bǔ)了D-InSAR技術(shù)的不足，在一定程度上解決了時空失相干、數(shù)字高程模型(digital elevation model,DEM)誤差、大氣延遲造成的影響，監(jiān)測精度可達(dá)厘米級甚至毫米級[12-13]。永久散射體測量技術(shù)(permanent scatters InSAR，PS-InSAR)利用影像上存在的永久散射點(diǎn)來替代整幅影像上的像素點(diǎn)，通過永久散射點(diǎn)的形變來呈現(xiàn)目標(biāo)區(qū)域的形變[14]，但是只適合小形變速率、小范圍的形變，無法很好地適用于沉陷區(qū)大范圍探測與監(jiān)測； 小基線集(small baesline subset,SBAS)-InSAR能以更高的時空密度來提取到相干點(diǎn)目標(biāo)，其形變結(jié)果更加明顯、直觀[15-16]，是目前礦區(qū)沉陷監(jiān)測最主要的方法之一。沙永蓮等[17]使用SBAS-InSAR技術(shù)對新疆哈密砂墩子煤田礦區(qū)進(jìn)行監(jiān)測,發(fā)現(xiàn)了一個沉降漏斗，時序結(jié)果顯示沉降漏斗先線性下沉后逐漸趨于穩(wěn)定； 李達(dá)等[18]通過對13景Terra SAR-X數(shù)據(jù)進(jìn)行SBAS-InSAR技術(shù)處理,監(jiān)測到了陜西榆林某煤礦的地表形變，并提取多個觀測點(diǎn)的時序沉降值進(jìn)行量化分析，結(jié)果表明 SBAS-InSAR技術(shù)可為礦區(qū)地表形變監(jiān)測和分析提供新的監(jiān)測手段。但是時序InSAR技術(shù)對影像的數(shù)量有一定要求，處理數(shù)據(jù)量較大，周期長，無法高效率地對大范圍的形變區(qū)域進(jìn)行監(jiān)測。針對礦區(qū)地表形變大范圍探測與監(jiān)測的問題，可以將D-InSAR技術(shù)和SBAS-InSAR技術(shù)結(jié)合起來使用： 首先利用D-InSAR技術(shù)對目標(biāo)區(qū)域進(jìn)行大范圍探測得到沉陷區(qū)位置，再利用時序InSAR技術(shù)對目標(biāo)沉陷區(qū)進(jìn)行時序監(jiān)測和形變追蹤，高效率地完成對目標(biāo)區(qū)域的形變監(jiān)測。

本文獲取了覆蓋山西省臨汾市全域不同軌道的C波段Sentinel-1A數(shù)據(jù)，采用聯(lián)合D-InSAR技術(shù)和SBAS-InSAR技術(shù)開展臨汾市沉陷區(qū)的大范圍探測和監(jiān)測工作，總計探測出沉陷區(qū)105處，選取其中的2處礦區(qū)進(jìn)行詳細(xì)分析，并結(jié)合光學(xué)影像進(jìn)行驗證。

1 研究方法

1.1 D-InSAR基本原理

D-InSAR技術(shù)通過采用目標(biāo)區(qū)域形變前后兩景影像共軛相乘獲得干涉圖，與DEM模擬相位差分消除地形相位后獲取到地表形變信息。但是在D-InSAR用于具體的地表形變監(jiān)測過程中，由于大氣環(huán)境和地表形變處于不斷的變化的過程，干涉相位會不可避免地受到影響，進(jìn)而影響到最終的測量結(jié)果[19]。差分相位表達(dá)式為：

φ=φflat+φtopo+φdef+φatm+φnoi，

(1)

式中：φ為形變相位；φflat為平地相位，可借助精密軌道數(shù)據(jù)去除；φtopo為地形相位，可借助DEM數(shù)據(jù)模擬地形相位去除；φdef為地表形變引起的相位；φatm為大氣延遲相位;φnoi為隨機(jī)噪聲相位，可采用濾波的方式進(jìn)行消除。

1.2 SBAS-InSAR原理

小基線技術(shù)利用同一軌道重復(fù)觀測獲取到目標(biāo)區(qū)域的SAR影像集，通過設(shè)定時空基線將SAR影像集進(jìn)行分組，形成多個小基線子集，再分別對每個子集內(nèi)的SAR影像進(jìn)行差分干涉生成差分干涉圖和相干性圖，利用平均相干性從相干性圖中選取高相干性的像素點(diǎn)，對差分干涉圖依次進(jìn)行相位解纏得到相位序列，聯(lián)合多個小基線子集的相位序列來預(yù)估形變信息。若在一定時間內(nèi)有N+1景 SAR 影像覆蓋目標(biāo)區(qū)域，設(shè)定時空基線后可以組合成M幅干涉圖，M需要滿足[20]:

(2)

以t0時刻作為初始時刻，則tk時刻的差分干涉相位φtk為觀察量，其中1

(3)

式中：d(ta,x)和d(tb,x)為ta和tb時刻雷達(dá)視線方向上的形變值;φ(ta,x)，φ(tb,x)分別為d(ta,x),d(tb,x)對應(yīng)得解纏相位值；λ為波長；Ti為主輔影像時間間隔；i為N幅影像獲取時間段內(nèi)某一時刻；v(x)為像元x處的形變速率。N幅影像可以組成M個方程，其方程組的形式為：

Aφ=δφ，

(4)

式中：A為M×N的系數(shù)矩陣，每一行代表一副干涉影像，每一列代表對應(yīng)時間上的SAR影像。若M>N，且A的秩為N，依據(jù)最小二乘法估值為：

(5)

式中ATA為奇異矩陣，存在無數(shù)解，可用聯(lián)合多個小基線采用奇異值分解法求出上式最小范數(shù)下的二乘解。

D-InSAR技術(shù)和SBAS-InSAR技術(shù)2種方法使用的場景不同，根據(jù)需求聯(lián)合2種監(jiān)測技術(shù)實(shí)現(xiàn)對于臨汾地區(qū)礦區(qū)的大范圍探測和動態(tài)監(jiān)測，技術(shù)流程如圖1所示。

圖1 技術(shù)流程圖

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

山西省臨汾市位于N35°23′～36°57′,E110°22′～112°34′之間，處于太原、鄭州、西安3個省會城市連線的中心點(diǎn)，地理區(qū)位優(yōu)勢明顯，礦產(chǎn)資源豐富； 地處屬溫帶大陸性氣候，四季分明，雨熱同期[21-22]； 地表多為黃土，土質(zhì)疏松，夏季多暴雨，地表植被容易遭到破壞。臨汾市包含河?xùn)|煤田北部、霍西煤田中部和沁水煤田西部，煤炭資源豐富，2019年原煤產(chǎn)量6 201.6萬t，占全省原煤產(chǎn)量的6.39%。臨汾市礦區(qū)多分布于臨汾斷陷盆地山體之中，地表以植被為主，多年的煤炭開采造成地表沉陷嚴(yán)重，容易引發(fā)滑坡、地面塌陷等多種地質(zhì)災(zāi)害，迫切需要有效的沉陷區(qū)監(jiān)測手段。

實(shí)驗選取432景不同軌道觀測獲取的Sentinel-1A影像(圖2)，其中path=11，frame=111有103景，path=11，frame=116有103景，時間跨度為2017年3月12日—2020年11月27日；path=113，frame=111有113景，path=113，frame=116有113景，時間跨度為2017年3月19日—2020年12月28日。影像均為C波段，周期為12 d，模式為升軌，VV極化，距離向分辨率為5 m，方位向分辨率為20 m。DEM數(shù)據(jù)為AW3D30 DEM，分辨率優(yōu)于30 m。衛(wèi)星軌道參數(shù)為歐空局提供的POD(AUX_POEORB，POD Precise Orbit Ephemerides)精確軌道參數(shù)。干涉雷達(dá)處理軟件為瑞士GAMMA軟件。

圖2 影像和DEM覆蓋范圍Fig.2 Image and DEM coverage

3 數(shù)據(jù)處理結(jié)果與分析

3.1 礦區(qū)大范圍探測

礦區(qū)形變具有連續(xù)性，受非人為因素的干擾較小，只要在礦區(qū)連續(xù)沉陷過程的某一時間段內(nèi)的形變過程被監(jiān)測到，則認(rèn)為該區(qū)域為沉陷區(qū)。本文選取2017年末、2019年初和2019年末3個時間段共12景不同軌道的Sentinel-1A影像組成覆蓋臨汾市全域的6組干涉對共3個組合(組合a、b、c)進(jìn)行大范圍探測。較短的時空基線可以在一定程度上減小時空失相干帶來的影響，為了保證干涉對具有較高的相干性，提高觀測精度，選取影像時時間基線設(shè)置為12 d，空間基線最大值控制在±100 m以內(nèi)。

在數(shù)據(jù)處理過程中，將距離向和方位向的多視比設(shè)定為4∶1，地形相位利用外部DEM去除，干涉對中存在的噪聲采用自適應(yīng)濾波的方法去除，相位使用奇異值分解法進(jìn)行解纏，大氣相位應(yīng)用GACOS技術(shù)去除，最終獲取到高質(zhì)量的解纏圖。由于礦區(qū)在短時間內(nèi)形變量級較大，形變相位在解纏圖中特征明顯，可直接依靠相位特征獲取到沉陷區(qū)位置，將沉陷區(qū)大致范圍用矩形框標(biāo)識出來。由于沉陷區(qū)在地理位置上呈現(xiàn)集群式分布，為了后續(xù)的數(shù)據(jù)方便處理將識別結(jié)果劃分成6部分，依次編號為A—F，結(jié)果如圖3所示。

圖3 D-InSAR解纏圖識別結(jié)果

圖3(a)—(c)中黃色虛線為獲取到的沉降區(qū)的大致位置。圖3(d)中紅色矩形框圈畫出來的即為沉陷區(qū)，形變相位表現(xiàn)為由青-紫-青變化； 黑色的空洞是因失相干而產(chǎn)生的監(jiān)測結(jié)果缺失。造成失相干的原因主要有2種： ①因沉陷區(qū)形變量過大，超出了C波段的監(jiān)測能力而產(chǎn)生失相干； ②因植被茂密導(dǎo)致地物介電質(zhì)常數(shù)發(fā)生變化而導(dǎo)致的失相干。沉陷區(qū)的也存在2種狀態(tài)： ①形變量級較大的沉陷區(qū)在解纏結(jié)果上表現(xiàn)為中心呈空洞狀態(tài)，周圍有形變相位環(huán)繞； ②形變量級較小的沉陷區(qū)表現(xiàn)為紫色或者黃色的形變相位環(huán)繞。沉陷區(qū)主要分布于臨汾斷陷盆地兩側(cè)的山體中，沿盆地東西兩側(cè)山體走向分布，呈點(diǎn)狀集群分布。分布于盆地中心的附近的沉陷區(qū)，地表植被較少，具有較高的相干性，獲得結(jié)果完整，易于識別沉陷區(qū)； 盆地兩側(cè)山體因植被影響造成大面積的失相干，給沉陷區(qū)識別帶來了一定的困難。C區(qū)域的沉陷區(qū)分布區(qū)域面積最大，數(shù)量最多，后文將著重以C區(qū)域作為重點(diǎn)進(jìn)行表述。

3.2 礦區(qū)時序監(jiān)測

基于D-InSAR識別到的沉陷區(qū)位置，依次對6個區(qū)域進(jìn)行SBAS-InSAR時序處理，數(shù)據(jù)處理詳細(xì)參數(shù)如表1所示。2個軌道公共主影像時間分別為2019年1月1日和2019年1月8日，將主影像作為統(tǒng)一的基準(zhǔn)進(jìn)行配準(zhǔn)并作為時間參考點(diǎn)，設(shè)定時間基線不少于90 d，空間基線不超過120 m。 為了提高結(jié)果的精度，使用強(qiáng)度互相關(guān)算法對主輔影像進(jìn)行配準(zhǔn)，通過設(shè)置不同的窗口大小、多項式參數(shù)個數(shù)等參數(shù)提高影像配準(zhǔn)精度。通過外部DEM去除地形相位，利用時空濾波方法去除大氣誤差，選取高相干性點(diǎn)位作為基準(zhǔn)進(jìn)行相位解纏，最終獲取到6個區(qū)域垂向的形變速率，結(jié)果如圖4所示。從圖4中可以看出，識別所得開采沉陷區(qū)范圍均在D-InSAR監(jiān)測結(jié)果中，其中灰色線條為手工勾繪的沉陷區(qū)邊界。本文此次探測出105處沉陷區(qū)，總面積達(dá)188.75 km2，其中7處與行政界線相交。識別所得沉陷區(qū)數(shù)量多，單個沉陷區(qū)沉陷面積較小，部分沉陷區(qū)形變量級較大，多處開采形變速率超過-100 mm/a，可以反映出臨汾市開采沉陷情況非常嚴(yán)重，這足以給當(dāng)?shù)氐V區(qū)生態(tài)環(huán)境造成巨大的破壞。

表1 數(shù)據(jù)處理詳細(xì)參數(shù)Tab.1 Detailed data processing parameters

圖4 臨汾市2017年3月—2020年12月SBAS-InSAR沉陷區(qū)形變速率圖Fig.4 Deformation rate map of SBAS-InSARsubsidence area in Linfen City fromMarch 2017 to December 2020

選取形變量級最大的C區(qū)域進(jìn)行詳細(xì)分析，獲取到的累計形變量和形變速率結(jié)果如圖5所示。C區(qū)域存在43處沉陷區(qū)，2017年3月19日—2020年12月28日該區(qū)域最大累計形變量為1 030 mm，最大形變速率為-381 mm/a。區(qū)域右側(cè)植被覆蓋較少，影像相干性較高，可以獲取到完整沉陷區(qū)形變，而左側(cè)區(qū)域為山區(qū)，植被覆蓋率高，白色區(qū)域即為影像因植被影響產(chǎn)生失相干所造成的結(jié)果缺失。其中C1沉陷區(qū)由2個相鄰的礦區(qū)構(gòu)成，左側(cè)沉陷區(qū)形變量級大于右側(cè)，最大累計形變量分別為922 mm和593 mm。C2沉陷區(qū)為單沉陷區(qū)，形變輪廓呈現(xiàn)為條帶狀，沉陷區(qū)的輪廓反映出礦區(qū)由西南向東北方向不斷開采，最大累計形變量為942 mm。

從圖6可以看出，C1，C2兩個沉陷區(qū)均處于不斷開采過程中。從形變過程來看，4個點(diǎn)位在礦區(qū)開采過程中形變量級在不斷增大，形變過程呈現(xiàn)快速沉陷-緩慢沉陷-快速沉陷的周期，以p1點(diǎn)的時序變化表現(xiàn)的最為明顯。在2017年5—8月，由于沒有獲取到有效的干涉對，故時序呈現(xiàn)擬合的直線。p1點(diǎn)整體形變速率最大，形變量總是經(jīng)歷較長時間快速形變到較短時間緩慢形變的周期循環(huán)，形變速率的絕對值也是呈大-小-大的周期變化，但是整體的形變趨勢逐漸緩和，最大沉降量為922 mm。p2點(diǎn)整體形變趨勢較為平穩(wěn)，周期性形變特征不太明顯，最大沉降量為593 mm。p3點(diǎn)位于C2礦區(qū)沉陷中心，p3點(diǎn)的整體形變速率要大于p4點(diǎn)，呈現(xiàn)快速下降-穩(wěn)定-略微抬升-緩慢下降-快速下降的過程，最大沉降量為942 mm。p4點(diǎn)位于礦區(qū)開采沉陷方向的中間位置，與p3點(diǎn)為同一礦區(qū)，整體形變趨勢與p3相似，最大沉降量為404 mm。

將C1，C2礦區(qū)分別沿圖6中ab，cd剖線截取剖面，結(jié)果如圖7所示。C1礦區(qū)有2個沉降漏斗，左側(cè)的沉降量要大于右側(cè)，2個沉陷區(qū)中間的位置已經(jīng)有100 mm的沉降量，已經(jīng)成為一個雙中心的大型沉陷區(qū)。C2礦區(qū)沉陷左側(cè)形變傾角大于右側(cè)，說明礦區(qū)開采的方向是沿剖線cd方向進(jìn)行開采，沉降漏斗中心位于起采區(qū)，但是隨著開采工作的進(jìn)行，沉降漏斗中心會逐漸向開采的方向移動，沉陷區(qū)輪廓會呈現(xiàn)為條帶狀。

由于沒有礦區(qū)的實(shí)測數(shù)據(jù)，本文利用Google Earth遙感影像數(shù)據(jù)檢驗監(jiān)測結(jié)果的一致性。首先采用閾值分割法將C1，C2沉陷區(qū)邊界提取出來，經(jīng)過多次嘗試當(dāng)閾值設(shè)置為-10 mm/a時，所獲得沉陷區(qū)邊界效果最好，將礦區(qū)邊界疊加至光學(xué)影像上如圖8所示。結(jié)合光學(xué)影像進(jìn)行分析，2個礦區(qū)均為地下開采，沉陷范圍內(nèi)有建筑物存在，地表分布有大量農(nóng)田，沉陷面積分別為1.54 km2，2.87 km2。C1礦區(qū)南側(cè)地勢平坦，已被開墾為農(nóng)田，北側(cè)為山地，人工地物較少，中間由一條道路橫穿沉陷區(qū)，左側(cè)100 m附近有一村落，周圍發(fā)現(xiàn)2處疑似采煤作業(yè)區(qū)。C2礦區(qū)地表起伏較大，但地表已被全部開墾為農(nóng)田，中間由多條道路穿過，東南側(cè)與小鎮(zhèn)相接，東北方向50 m有一條國道，北側(cè)分布有3個疑似采煤作業(yè)區(qū)。C1,C2礦區(qū)沉陷邊界直接與農(nóng)田、道路相接，距離住宅不足100 m，產(chǎn)生的地面沉陷會直接威脅到周圍居民的人身財產(chǎn)安全，尤其C2礦區(qū)，需要持續(xù)進(jìn)行監(jiān)測。

4 結(jié)語

煤礦開采所造成的沉陷區(qū)對地表生態(tài)環(huán)境和基礎(chǔ)設(shè)施帶來了巨大的破壞。本文利用D-InSAR技術(shù)對臨汾市沉陷區(qū)進(jìn)行大范圍探測和監(jiān)測結(jié)果工作，將探測結(jié)果分成了6個區(qū)域并依次使用SBAS-InSAR技術(shù)進(jìn)行時序監(jiān)測，共發(fā)現(xiàn)105處沉陷區(qū)； 將其中2處典型沉陷區(qū)作為重點(diǎn)區(qū)域進(jìn)一步分析，獲取到了該沉陷區(qū)的形變演化情況，為后續(xù)礦區(qū)沉陷災(zāi)害防治提供重要依據(jù)。通過光學(xué)影像發(fā)現(xiàn)了沉陷區(qū)附近存在疑似采煤作業(yè)區(qū)，驗證了InSAR技術(shù)大范圍探測與監(jiān)測具有可行性和可靠性。

1)礦區(qū)地表形變具有形變速率快、量級大的特點(diǎn)，會對當(dāng)?shù)氐纳鷳B(tài)環(huán)境造成巨大的破壞，難以通過傳統(tǒng)的光學(xué)遙感進(jìn)行探測，InSAR技術(shù)可以極易捕獲沉陷區(qū)的大量級地表形變，從而實(shí)現(xiàn)對研究區(qū)沉陷區(qū)的大范圍探測。

2)直接利用時序InSAR技術(shù)對大范圍區(qū)域進(jìn)行沉陷區(qū)探測具有巨大的挑戰(zhàn)性，處理數(shù)據(jù)量大，工作效率低，結(jié)果不穩(wěn)定。而基于D-InSAR技術(shù)的大范圍探測和SBAS-InSAR技術(shù)的詳細(xì)監(jiān)測的方法可以很快實(shí)現(xiàn)對研究區(qū)進(jìn)行層次化的處理，真正做到對大范圍區(qū)域的高效率監(jiān)測。

3)本文基于InSAR技術(shù)使用Sentinel-1A升軌數(shù)據(jù)對臨汾市礦區(qū)地表形變進(jìn)行了大范圍探測和監(jiān)測，識別和監(jiān)測結(jié)果對臨汾市礦區(qū)形變監(jiān)測具有較好的參考價值，但是InSAR技術(shù)在植被覆蓋率較高的區(qū)域適用性不佳，未來還需要借助其他技術(shù)來提高大范圍探測和監(jiān)測結(jié)果效率。

志謝：感謝歐洲太空局在科學(xué)數(shù)據(jù)中心網(wǎng)站公布Sentinel-1A衛(wèi)星雷達(dá)影像，感謝日本宇宙航空航天局地球觀測研究中心提供DEM數(shù)據(jù)。