基于時(shí)序InSAR的貴州喀斯特地區(qū)地表形變監(jiān)測(cè)
——以清鎮(zhèn)市為例

王 碩，張正加，范 鵬,王猛猛

(1.自然資源部測(cè)繪發(fā)展研究中心，北京 100830；2.武昌理工學(xué)院 人工智能學(xué)院,湖北 武漢 430223；3.國(guó)網(wǎng)電力科學(xué)研究院武漢南瑞有限責(zé)任公司，湖北 武漢 430074)

0 引 言

地質(zhì)災(zāi)害會(huì)對(duì)人民的生命安全和社會(huì)安全產(chǎn)生嚴(yán)重的威脅，嚴(yán)重阻礙社會(huì)的發(fā)展。我國(guó)自然災(zāi)害頻發(fā)，是世界上地質(zhì)災(zāi)害最嚴(yán)重、受威脅人口最多的國(guó)家之一，其中貴州省是我國(guó)受地質(zhì)災(zāi)害影響最為嚴(yán)重的省份之一[1]。貴州的地質(zhì)災(zāi)害主要以滑坡、泥石流、塌陷等地質(zhì)災(zāi)害為主，對(duì)旅游景區(qū)和人口密集區(qū)危害極大。因此對(duì)貴州地區(qū)進(jìn)行地表形變監(jiān)測(cè)具有十分重要的意義。

雖然傳統(tǒng)的地質(zhì)災(zāi)害監(jiān)測(cè)手段( GPS、水準(zhǔn)測(cè)量等)能夠有效地獲取地面點(diǎn)形變信息，但其通常工作量大，且無法進(jìn)行長(zhǎng)時(shí)間、大范圍觀測(cè)。相對(duì)于傳統(tǒng)地質(zhì)災(zāi)害識(shí)別方法，遙感是一種現(xiàn)代高科技手段，具有客觀、真實(shí)的特點(diǎn)，能為地質(zhì)災(zāi)害調(diào)查監(jiān)測(cè)提供重要的技術(shù)手段，能客觀、真實(shí)、準(zhǔn)確的提取地質(zhì)災(zāi)害數(shù)據(jù),為地質(zhì)災(zāi)害的識(shí)別防治做出巨大貢獻(xiàn)[2]。合成孔徑雷達(dá)干涉測(cè)量(InSAR)技術(shù)作為雷達(dá)遙感的一個(gè)重要分支，具有較強(qiáng)的形變測(cè)量能力，其形變監(jiān)測(cè)精度能夠達(dá)到厘米級(jí)甚至毫米級(jí)[3-4]。InSAR技術(shù)具有監(jiān)測(cè)范圍廣、精度高等優(yōu)點(diǎn)，已經(jīng)逐漸在地質(zhì)災(zāi)害監(jiān)測(cè)應(yīng)用得到廣泛認(rèn)可[1,5]。為了克服傳統(tǒng)InSAR技術(shù)固有的時(shí)間、大氣去相干的影響[6]，研究學(xué)者提出時(shí)序InSAR方法，如永久散射體干涉測(cè)量(PSInSAR)[7]、小基線集(SBAS)[8]，并在城市沉降、地表形變、火山、地震等方面得到廣泛應(yīng)用[9-10]。近些年，時(shí)序InSAR技術(shù)已應(yīng)用于山區(qū)地質(zhì)災(zāi)害監(jiān)測(cè)，特別是在四川、貴州、云南等地區(qū)得到應(yīng)用[1,11-12]。

清鎮(zhèn)市是貴陽市下轄縣級(jí)市，位于貴州省中部，是貴州省地質(zhì)災(zāi)害重災(zāi)區(qū)之一，目前清鎮(zhèn)市地質(zhì)災(zāi)害防治形勢(shì)嚴(yán)峻，開展地質(zhì)災(zāi)害監(jiān)測(cè)是一項(xiàng)重要的工作[13]。因此，本文以貴陽清鎮(zhèn)市為例，利用SBAS技術(shù)，使用26景Sentinel-1A數(shù)據(jù)，對(duì)貴州清鎮(zhèn)市地區(qū)2018～2019年的地表形變監(jiān)測(cè)，并識(shí)別出可能的地質(zhì)災(zāi)害區(qū)域，并對(duì)其中形變嚴(yán)重地區(qū)進(jìn)行詳細(xì)分析。

1 清鎮(zhèn)地理概況

清鎮(zhèn)地處省會(huì)貴陽市以西，是貴陽市下轄的縣級(jí)市，地處苗嶺山脈北坡，烏江干流鴨池河?xùn)|岸，東經(jīng)106°07′～106°32′，北緯26°25′～26°56′之間[14]。東與白云區(qū)、烏當(dāng)區(qū)毗鄰，南與安順市平壩區(qū)接壤，西為東風(fēng)湖與織金縣相連，北為鴨池河、貓?zhí)优c黔西縣、修文縣隔岸相望，總面積1 383 km2。清鎮(zhèn)地區(qū)植被覆蓋茂密，覆蓋面積達(dá)到了48%，如圖1(a)所示。清鎮(zhèn)市海拔變化從782 m至1 738 m，平均海拔1 200 m左右，如圖1(b)所示。清鎮(zhèn)市是貴州省地質(zhì)災(zāi)害重災(zāi)區(qū)之一，高易發(fā)區(qū)主要分布在北部、西北部和中部的新店鎮(zhèn)、流長(zhǎng)鄉(xiāng)、站街鎮(zhèn)等鄉(xiāng)。

圖1 研究區(qū)位置范圍

2 研究數(shù)據(jù)及處理

2.1 數(shù)據(jù)介紹

哨兵1號(hào)(Sentinel-1)衛(wèi)星是歐洲航天局哥白尼計(jì)劃(GMES)中的地球觀測(cè)衛(wèi)星，為雙星系統(tǒng)，搭載C波段合成孔徑雷達(dá)，可全天時(shí)、全天候不受天氣條件影響地提供連續(xù)影像[15]。本文使用數(shù)據(jù)為Sentinel-1A衛(wèi)星干涉寬幅模式的斜距單視復(fù)數(shù)產(chǎn)品(IW SLC)，圖像的方位向和距離向分辨率為20 m和5 m，極化方式為VV。共收集了26景覆蓋貴陽清鎮(zhèn)市的哨兵數(shù)據(jù)，成像時(shí)間從2018年8月1日至2019年7月15日(表1)。清鎮(zhèn)市的雷達(dá)圖像如圖2所示，在雷達(dá)圖像中，水體表現(xiàn)為暗色，城市區(qū)域表現(xiàn)為亮色。本文收集了研究區(qū)的30 m分辨率的SRTM DEM數(shù)據(jù)，用于去除地形相位[16]。另外本文還收集了清鎮(zhèn)市的Landsat-8影像，用于地質(zhì)災(zāi)害識(shí)別與分析。

表1 SAR數(shù)據(jù)列表

圖2 研究區(qū)SAR幅度圖

2.2 SAR數(shù)據(jù)處理

研究區(qū)清鎮(zhèn)市植被覆蓋比較茂密，在干涉處理中容易引起時(shí)間去相干問題，導(dǎo)致反演結(jié)果精度較差。為了減小失相干影響，獲得準(zhǔn)確的地表沉降信息，本文使用小基線集(SBAS)的方法進(jìn)行時(shí)序相干處理。本文的數(shù)據(jù)處理流程如圖3所示，在本節(jié)重點(diǎn)介紹SBAS技術(shù)的關(guān)鍵步驟：干涉圖生成、相干點(diǎn)選取和相位模型建立。

圖3 數(shù)據(jù)處理流程圖

2.2.1 干涉圖生成

獲取同一區(qū)域按照時(shí)間序列(t0,……,tN)排列的N+1幅SAR圖像，選取其中的一幅影像作為主圖像，并將其他SAR影像與主圖像進(jìn)行配準(zhǔn)。根據(jù)干涉條件組合，得到M幅干涉圖和相干圖，以保證相干性的穩(wěn)定。利用引入外部DEM，將干涉條紋圖中的地形相位去掉得到差分干涉條紋圖。為了抑制噪聲的影響，使用自適應(yīng)濾波方法對(duì)干涉相位進(jìn)行濾波[17]。生成干涉條紋圖以后，使用最小費(fèi)用流的方法對(duì)差分干涉相位進(jìn)行解纏[18]。

2.2.2 高相干點(diǎn)選取

選擇適當(dāng)?shù)南喔上禂?shù)閾值，將N幅干涉圖中，平均相干系數(shù)大于該閾值的點(diǎn)認(rèn)為是相位穩(wěn)定的點(diǎn)，即高相干點(diǎn)，如下式：

(1)

式中，ri表示單個(gè)像素在第i幅干涉圖上的相干系數(shù)，thd_r是相干系數(shù)閾值。要注意，高閾值會(huì)舍棄許多有用的相干點(diǎn)，低的閾值又會(huì)帶來許多噪聲。在本文中，選擇0.7作為相干點(diǎn)的選取閾值。

2.2.3 相位模型建立

假設(shè)從tA，tB兩個(gè)時(shí)間獲得的SAR圖像產(chǎn)生第j幅干涉圖，并假設(shè)tB>tA，去除平地和地形相位后，在x處的差分干涉相位[19]：

(2)

式中，d(tB,x)和d(tA,x)是相對(duì)于參考時(shí)間t1的LOS方向累計(jì)形變，d(t1,x)=0，φ(tB,x)和φ(tA,x)分別表示d(tB,x)和d(tA,x)引起的形變相位，λ是雷達(dá)波長(zhǎng)，則利用如下線性模型估計(jì)N幅圖像的形變：

AΦ=ΔΦ

(3)

式中，Φ表示待求點(diǎn)上的N個(gè)時(shí)刻的SAR圖像上的未知形變相位組成的矩陣；ΔΦ為M幅差分干涉圖上相位值組成的矩陣；系數(shù)矩陣A[M×N]每行對(duì)應(yīng)于一幅干涉圖，每列對(duì)應(yīng)于一個(gè)時(shí)間上的SAR圖像，主圖像所在列為+1，輔圖像所在列為-1，其余列為0。如果M>=N，且A的秩是N，則利用最小二乘法可得：

Φ=A#·ΔΦ,A#=(ATA)-1AT

(4)

上式給出了理想情況下的最優(yōu)解。但實(shí)際上對(duì)于多個(gè)小基線集，ATA是一個(gè)奇異矩陣，如假設(shè)有L個(gè)子集，A的秩為N-L+1，方程組就會(huì)有無窮多解(設(shè)N≤M)。為了解決這個(gè)問題，SBAS算法采用矩陣的奇異值分解(SVD)技術(shù)從而求出方程(4)在最小范數(shù)意義上的最小二乘解。

將相位轉(zhuǎn)化到平均相位速度：

(5)

代替式(2)中的相位，于是有：

∑(tk+1-tk)vk=Δφj,j=1,...,M

(6)

上述方程組也可表達(dá)為：Dv=ΔΦ，其中D是一個(gè)M×(N-1)矩陣，表示主輔圖像之間干涉組合的時(shí)間基線，對(duì)其進(jìn)行SVD分解可以得到速度矢量v的最小范數(shù)解。在線性形變估計(jì)的基礎(chǔ)上，繼續(xù)對(duì)殘余相位進(jìn)行合適的時(shí)間和空間濾波處理即能夠分離出大氣相位和非線性形變相位成分。

3 實(shí)驗(yàn)結(jié)果及分析

3.1 清鎮(zhèn)市年平均沉降速率

對(duì)26景Sentinel-1A數(shù)據(jù)進(jìn)行配準(zhǔn)、預(yù)濾波等操作，設(shè)置相關(guān)的時(shí)間基線閾值(小于100 d)和空間基線閾值(小于60 m)，然后對(duì)生成的每個(gè)干涉圖進(jìn)行檢查，剔除相位噪聲較大的干涉對(duì)，最終得到78對(duì)干涉條紋圖(圖4)。然后根據(jù)計(jì)算得到的平均相干系數(shù)圖，設(shè)置相關(guān)的閾值，得到研究區(qū)的高相干點(diǎn)。共選取了577 374個(gè)監(jiān)測(cè)點(diǎn)，監(jiān)測(cè)點(diǎn)密度達(dá)到了386/km2，選取的高相干點(diǎn)主要集中在城鎮(zhèn)及周邊區(qū)域，在山區(qū)選取的高相干點(diǎn)比較少。

圖4 時(shí)間基線和空間基線分布

根據(jù)上述的方法和處理流程得到了清鎮(zhèn)市201808～201907的年沉降速率結(jié)果，如圖5所示。從圖中可以看到，整體上清鎮(zhèn)市的年沉降速率在-40～10 mm/year，大部分沉降速率在-10至10 mm/year區(qū)間，表明整體上比較穩(wěn)定。沉降速率在空間上分布具有不均勻性，在清鎮(zhèn)市的西部(犁倭鎮(zhèn))和中部(麥格苗族布依族鄉(xiāng))發(fā)現(xiàn)了沉降明顯的區(qū)域，最大沉降速率達(dá)到了-40 mm/year。另外在衛(wèi)城鎮(zhèn)也監(jiān)測(cè)到比較明顯的沉降。

圖5 清鎮(zhèn)市201807～201907年平均沉降速率散點(diǎn)圖

圖6是根據(jù)沉降散點(diǎn)差值得到的清鎮(zhèn)市沉降速率差值圖，圖中的線是沉降速率等值線。根據(jù)差值圖，可以很明顯地發(fā)現(xiàn)沉降區(qū)域分布的空間特征。沉降區(qū)域主要集中在清鎮(zhèn)市的西北部和中部，在北部也零星分布了一些沉降區(qū)域，在這些區(qū)域應(yīng)當(dāng)加強(qiáng)地質(zhì)災(zāi)害防治，避免在預(yù)計(jì)發(fā)生山體滑坡等地質(zhì)災(zāi)害。

圖6 清鎮(zhèn)市沉降速率差值圖

參考北京地區(qū)地面沉降危險(xiǎn)性分級(jí)標(biāo)準(zhǔn)[20]，按照地面沉降速率大小將清鎮(zhèn)分為地面沉降輕微區(qū)、一般區(qū)和較為嚴(yán)重區(qū)域，對(duì)應(yīng)的沉降速率區(qū)間為：小于10 (mm·a-1),10～30(mm·a-1),大于30(mm·a-1)，統(tǒng)計(jì)結(jié)果如表2所示。從表2可見，沉降速率超過-30 (mm·a-1)的面積超過了0.42 km2，約占清鎮(zhèn)市面積的0.03%。沉降速率在-30至-10mm/year的面積超過了54 km2，約占清鎮(zhèn)市面積的3%,可以看到有明顯沉降的區(qū)域面積比較大。

表2 沉降區(qū)間面積統(tǒng)計(jì)

3.2 區(qū)域分析

3.2.1 犁倭鎮(zhèn)

本文提取的犁倭鎮(zhèn)的平均沉降速率圖如圖7(a)所示，在犁倭鎮(zhèn)的西南部分，出現(xiàn)了較為明顯的沉降現(xiàn)象，最大的形變速率超多了-30 mm·a-1。在犁倭鎮(zhèn)的東北則比較穩(wěn)定，整個(gè)觀測(cè)周期內(nèi)的形變量小于5 mm·a-1。選取了3個(gè)點(diǎn)用于形變趨勢(shì)分析，位置如圖7(a)，時(shí)序形變?nèi)鐖D7(c)。PS點(diǎn)1和點(diǎn)2的形變比較明顯，并且形變趨勢(shì)比較一致，觀測(cè)周期內(nèi)的形變分辨達(dá)到了33 mm和44 mm。從2018年8月至2019年5月形變速率比較緩慢，從2019年5月以后，形變速率明顯加速。通常這種形變的加速分析可能與降雨量有關(guān)，特別是在山區(qū)[12]，這些區(qū)域需要注意，特別是在雨季，容易出現(xiàn)地質(zhì)災(zāi)害隱患。PS3點(diǎn)在整個(gè)觀測(cè)周期內(nèi)則比較穩(wěn)定。

(a)犁倭鎮(zhèn)平均形變速率圖；(b)犁倭鎮(zhèn)光學(xué)圖像；(c)選取3個(gè)點(diǎn)的累積形變量圖7 犁倭鎮(zhèn)沉降信息

3.2.2 麥格苗族布依族鄉(xiāng)

麥格苗族布依族鄉(xiāng)2018年8月至2019年7月的形變速率如圖8(a)所示，可以看到在鄉(xiāng)鎮(zhèn)區(qū)域監(jiān)測(cè)到了比較明顯的沉降現(xiàn)象。我們選取了在形變區(qū)域選取了2個(gè)點(diǎn)(PS4和PS5)分析其時(shí)序形變規(guī)律。圖8(c)展示了選取2個(gè)點(diǎn)的時(shí)間序列形變量，PS4和PS5兩點(diǎn)的形變速率分別為25.4 mm/year和29.1 mm/year。需要主要的是，PS4和PS5兩點(diǎn)的形變趨勢(shì)比較一致，在觀測(cè)周期內(nèi)呈現(xiàn)線性變化，形變量分別達(dá)到了21 mm和31 mm。圖8(d)是形變區(qū)域比較嚴(yán)重的光學(xué)放大圖，我們可以看到在這些區(qū)域處于開發(fā)狀態(tài)，地表植被遭到破壞，分析地表沉降可能與當(dāng)?shù)氐牡V產(chǎn)資源開采活動(dòng)有關(guān)。

(a)麥格苗族布依族鄉(xiāng)平均形變速率圖；(b)麥格苗族布依族鄉(xiāng)光學(xué)圖像；(c)選取兩個(gè)點(diǎn)的累積形變量；(d)形變區(qū)域的光學(xué)放大圖圖8 麥格苗族布依族鄉(xiāng)沉降信息

3.2.3 衛(wèi)城鎮(zhèn)

圖9(a)是衛(wèi)城鎮(zhèn)的平均沉降速率圖，可以看到大部分區(qū)域都比較穩(wěn)定。同樣選取了PS6和PS7兩個(gè)點(diǎn)對(duì)該區(qū)域進(jìn)行分析。其中，PS6點(diǎn)為沉降比較明顯的區(qū)域，PS7點(diǎn)位于穩(wěn)定區(qū)域。可以發(fā)現(xiàn)，在PS6點(diǎn)的累積形變量達(dá)到了-20 mm，而PS7點(diǎn)的累積形變量在0～10 mm之間起伏。另外，從整體上看，圍城鎮(zhèn)的形變速率明顯小于犁倭鎮(zhèn)和麥格苗族布依族鄉(xiāng)。

3.3 沉降因素分析

根據(jù)上述分析，可以看到清鎮(zhèn)市的地表沉降主要集中在中部和西部。根據(jù)清鎮(zhèn)市的地理環(huán)境，其沉降主要有3方面的因素：① 清鎮(zhèn)市屬于喀斯特巖溶地區(qū)，土地比較脆弱，土地容易石漠化[13-14],土地承載力較差，土層比較松散，在上層作荷載的作用下表面容易發(fā)生形變[12]；② 降雨的影響，雨水會(huì)進(jìn)一步壓實(shí)和沖刷松散土層，加速地表的形變。在犁倭鎮(zhèn)PS1和PS2兩點(diǎn)觀測(cè)到了在雨季形變加速的現(xiàn)象；③ 人類活動(dòng)的影響。清鎮(zhèn)市有豐富的礦產(chǎn)資源，礦產(chǎn)資源的開采會(huì)引起比較明顯的地表形變，如在麥格苗族布依族鄉(xiāng)地區(qū)有豐富的礦產(chǎn)資源，當(dāng)?shù)氐牡V產(chǎn)開采活動(dòng)一直在持續(xù)，一定程度上導(dǎo)致了該地區(qū)的沉降。

(a)衛(wèi)城鎮(zhèn)平均形變速率圖；(b)衛(wèi)城鎮(zhèn)光學(xué)圖像；(c)選取兩個(gè)點(diǎn)的累積形變量圖9 衛(wèi)城鎮(zhèn)沉降信息

4 結(jié) 語

本文利用SBAS-InSAR技術(shù)和26景Sentinel-1A數(shù)據(jù)，提取得到了貴陽清鎮(zhèn)市2018年8月至2019年7月的地表沉降形變場(chǎng)。研究結(jié)果表明，研究區(qū)域內(nèi)大部分地區(qū)表現(xiàn)較為穩(wěn)定的形變，形變速率在-10 mm·a-1以內(nèi)。沉降比較明顯的地方主要集中在犁倭鎮(zhèn)、麥格苗族布依族鄉(xiāng)、衛(wèi)城鎮(zhèn)等3個(gè)鄉(xiāng)鎮(zhèn)地區(qū)，最大形變速率超過了-40 mm·a-1。沉降速率超過-10 mm·a-1的面積得到了54.7 km2,占全市面積的3%左右。在犁倭鎮(zhèn)，發(fā)現(xiàn)了沉降信號(hào)在雨季有增加的趨勢(shì)，分析這種形變?cè)黾于厔?shì)與降雨有關(guān)；另外發(fā)現(xiàn)了地質(zhì)環(huán)境和人類活動(dòng)有關(guān)等因素對(duì)形變的影響。

后期將繼續(xù)利用SBAS技術(shù)對(duì)該地區(qū)的沉降持續(xù)監(jiān)測(cè)，特別是地面沉降嚴(yán)重區(qū)域如衛(wèi)城鎮(zhèn)、犁倭鎮(zhèn)、麥格苗族布依族鄉(xiāng)，為清鎮(zhèn)市的控沉規(guī)劃的制定提供科學(xué)合理的建議。