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

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

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

0 引 言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.2.1 干涉圖生成

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

2.2.2 高相干點選取

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

(1)

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

2.2.3 相位模型建立

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

(2)

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

AΦ=ΔΦ

(3)

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

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

(4)

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

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

(5)

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

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

(6)

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

3 實驗結(jié)果及分析

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

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

圖4 時間基線和空間基線分布

根據(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ū)域，最大沉降速率達到了-40 mm/year。另外在衛(wèi)城鎮(zhèn)也監(jiān)測到比較明顯的沉降。

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

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

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

參考北京地區(qū)地面沉降危險性分級標準[20]，按照地面沉降速率大小將清鎮(zhèn)分為地面沉降輕微區(qū)、一般區(qū)和較為嚴重區(qū)域，對應的沉降速率區(qū)間為：小于10 (mm·a-1),10～30(mm·a-1),大于30(mm·a-1)，統(tǒng)計結(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)計

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)定，整個觀測周期內(nèi)的形變量小于5 mm·a-1。選取了3個點用于形變趨勢分析，位置如圖7(a)，時序形變?nèi)鐖D7(c)。PS點1和點2的形變比較明顯，并且形變趨勢比較一致，觀測周期內(nèi)的形變分辨達到了33 mm和44 mm。從2018年8月至2019年5月形變速率比較緩慢，從2019年5月以后，形變速率明顯加速。通常這種形變的加速分析可能與降雨量有關，特別是在山區(qū)[12]，這些區(qū)域需要注意，特別是在雨季，容易出現(xiàn)地質(zhì)災害隱患。PS3點在整個觀測周期內(nèi)則比較穩(wěn)定。

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

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

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

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

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

圖9(a)是衛(wèi)城鎮(zhèn)的平均沉降速率圖，可以看到大部分區(qū)域都比較穩(wěn)定。同樣選取了PS6和PS7兩個點對該區(qū)域進行分析。其中，PS6點為沉降比較明顯的區(qū)域，PS7點位于穩(wěn)定區(qū)域。可以發(fā)現(xiàn)，在PS6點的累積形變量達到了-20 mm，而PS7點的累積形變量在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]；② 降雨的影響，雨水會進一步壓實和沖刷松散土層，加速地表的形變。在犁倭鎮(zhèn)PS1和PS2兩點觀測到了在雨季形變加速的現(xiàn)象；③ 人類活動的影響。清鎮(zhèn)市有豐富的礦產(chǎn)資源，礦產(chǎn)資源的開采會引起比較明顯的地表形變，如在麥格苗族布依族鄉(xiāng)地區(qū)有豐富的礦產(chǎn)資源，當?shù)氐牡V產(chǎn)開采活動一直在持續(xù)，一定程度上導致了該地區(qū)的沉降。

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

4 結(jié) 語

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

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