龍里某塌陷時(shí)序InSAR變形監(jiān)測(cè)的PS修正

羅雪瑋， 向喜瓊，2， 呂亞東

(1.貴州大學(xué)喀斯特地質(zhì)資源與環(huán)境教育部重點(diǎn)實(shí)驗(yàn)室,貴陽 550025；2.貴州大學(xué)資源與環(huán)境工程學(xué)院，貴陽 550025)

0 引言

為克服差分合成孔徑雷達(dá)干涉測(cè)量技術(shù)(differential interferometric synthetic aperture Radar，D-InSAR)時(shí)空失相干、大氣效應(yīng)及噪聲等因素的影響，一系列的時(shí)序合成孔徑雷達(dá)干涉測(cè)量技術(shù)(interferometric synthetic aperture Radar，InSAR)技術(shù)應(yīng)運(yùn)而生。具有代表性的有Stacking-InSAR、永久散射體(permanent scatter，PS)PS-InSAR、小基線集(small baseline subset，SBAS)SBAS-InSAR以及相干永久目標(biāo)分析(interferometric persistent target analysis)IPTA-InSAR等技術(shù)[1]。Colombo等[2]利用ESA ERS1/2衛(wèi)星采集的合成孔徑雷達(dá)(synthetic aperture Radar，SAR)數(shù)據(jù)，采用D-InSAR和PS-InSAR進(jìn)行了干涉分析，在監(jiān)測(cè)的工業(yè)區(qū)內(nèi)發(fā)現(xiàn)了不同的地面沉降模式； Wu等[3]采用SBAS-InSAR方法對(duì)延安新區(qū)2015—2019年的地表變形進(jìn)行了調(diào)查，發(fā)現(xiàn)區(qū)內(nèi)地基沉降主要來源于填方和城市荷載的作用，而建筑物的應(yīng)力釋放是造成地面隆起的主要因素； 劉一霖等[4]基于SBAS-InSAR技術(shù)，引入偏移量追蹤法、FFT過采樣方法、多窗口迭代自適應(yīng)濾波技術(shù)對(duì)礦區(qū)進(jìn)行了地表時(shí)序監(jiān)測(cè)，分析了研究區(qū)開采工作面內(nèi)地表的時(shí)空變化規(guī)律； 葛大慶等[5]融合了PS-InSAR和SBAS-InSAR兩種方法，以短基線差分干涉生成相位圖后利用點(diǎn)目標(biāo)識(shí)別技術(shù)提取出相干目標(biāo)，分析了德州地區(qū)地面沉降-回彈的動(dòng)態(tài)變化過程； 李金超等[6]將基于灰色模型和支持向量回歸的組合模型(gray model and support vector regression，MG-SVR)與SBAS-InSAR技術(shù)相結(jié)合，對(duì)淮南潘集礦區(qū)發(fā)生沉降的區(qū)域進(jìn)行了監(jiān)測(cè)及預(yù)警。

在使用SBAS方法對(duì)地表形變進(jìn)行監(jiān)測(cè)時(shí)，需要在軌道精煉步驟中引入穩(wěn)定的地面控制點(diǎn)(ground control point，GCP)估算和去除殘余的恒定相位和解纏后還存在的相位坡道，從而使形變監(jiān)測(cè)結(jié)果更為準(zhǔn)確。但在沒有精確GCP的情形下，如果研究區(qū)內(nèi)存在較大高差、或者影像數(shù)據(jù)噪聲過大，無法進(jìn)行消除，人工來選擇GCP可能會(huì)引起較大誤差。

針對(duì)此問題，本文提出了一種基于PS的SBAS-InSAR方法。該方法通過設(shè)置相干系數(shù)的閾值、振幅離差指數(shù)的閾值以及地表形變速率的閾值選出穩(wěn)定的PS點(diǎn)，并將這些點(diǎn)作為SBAS-InSAR軌道精煉中的地面控制點(diǎn)，以期能得到準(zhǔn)確度更高的地表形變監(jiān)測(cè)結(jié)果。最后，以位于貴州省龍里縣洗馬鎮(zhèn)上石坎村的一處地面塌陷為實(shí)例驗(yàn)證了本文方法的有效性。

1 研究方法

1.1 相干系數(shù)法

相干系數(shù)法的基本思路為： 對(duì)于某一像元選定一定大小的移動(dòng)窗口(m,n)，利用該移動(dòng)窗口內(nèi)的像素復(fù)數(shù)信息來估計(jì)相干系數(shù)值。計(jì)算完畢后，移動(dòng)該窗口，重復(fù)計(jì)算[7-8]。相關(guān)系數(shù)γ計(jì)算公式為：

(1)

式中:M為主影像;S為輔影像;*為共軛相乘。

(2)

1.2 振幅離差法

Ferretti等[9]提出可以采用像素中振幅信息的時(shí)間序列來度量干涉相位的穩(wěn)定性。R和I分別為影像的實(shí)部和虛部，若存在標(biāo)準(zhǔn)差為σn的高斯噪聲，則振幅值A(chǔ)服從Rice分布，公式為：

(3)

(4)

式中：v為干涉相位；σv為相位離差指數(shù)；σnI為虛部的標(biāo)準(zhǔn)差；mA和σA為時(shí)序振幅的均值和標(biāo)準(zhǔn)差；DA為振幅離差指數(shù)。

1.3 相干性、振幅離差指數(shù)以及形變速率結(jié)合的選點(diǎn)方法

首先采用相干系數(shù)閾值法初步識(shí)別出具有高相干系數(shù)的像元的PS點(diǎn)目標(biāo)； 然后設(shè)定振幅離差值進(jìn)一步篩選具有穩(wěn)定強(qiáng)散射性的PS點(diǎn)作為目標(biāo)； 最后設(shè)定一個(gè)形變速率區(qū)間來進(jìn)行最終的PS點(diǎn)選取。如圖1所示。

圖1 方法流程Fig.1 Method and procedure

1.4 SBAS-InSAR技術(shù)

SBAS-InSAR技術(shù)可以限制時(shí)空基線失相干和大氣延遲的影響[10]。將N+1景覆蓋同一區(qū)域的SAR影像序列，影像獲取時(shí)間為(t0,…，tN)，根據(jù)空間基線和時(shí)間基線的條件組合干涉對(duì)得到M幅差分干涉圖，M的數(shù)量與影像數(shù)量滿足以下關(guān)系[11]：

(5)

假設(shè)干涉圖j由tA和tB時(shí)刻獲取的影像干涉生成，在去除平地效應(yīng)與地形相位影響后，干涉圖j中距離向?yàn)閞與方位向?yàn)閤的某一像素的干涉相位表示為：

(6)

式中：φ(tB,x,r)和φ(tA,x,r)分別為tB和tA時(shí)刻SAR影像上的相位值；φdef,j(x,r)為tA～tB時(shí)刻內(nèi)衛(wèi)星視線向(los)的形變相位；φtoppo,j(x,r)為地形相位誤差；φatm,j(x,r)為大氣相位誤差；φnoise,j(x,r)為噪聲相位。

對(duì)所有干涉圖構(gòu)建差分干涉相位與影像獲取時(shí)刻形變量之間的線性方程組，其矩陣形式為：

(7)

式中:A為M×N的系數(shù)矩陣;φ為參數(shù)矩陣，由未知像元(x,y)在N個(gè)時(shí)刻對(duì)應(yīng)的未知形變相位值構(gòu)成；δφ為由M個(gè)式(6)所示的δφj組成的觀測(cè)方程組；M為差分干涉圖的數(shù)量。

為求解研究區(qū)域各高相干點(diǎn)的形變速率，可用兩景影像間的平均相位速率代替相位值，則式(7)變?yōu)椋?/p>

(8)

式中：B為M×N的系數(shù)矩陣；v為形變速率向量，vΤ可表示為：

(9)

當(dāng)系數(shù)矩陣B為滿秩(即M≥N)時(shí)，可用最小二乘法求解形變速率； 當(dāng)M

2 實(shí)驗(yàn)數(shù)據(jù)和處理

2.1 數(shù)據(jù)

本文以龍里縣洗馬鎮(zhèn)的一處地面塌陷為實(shí)例(圖2)，選取20景由歐洲航天局提供的C波段Sentinel-1A SAR影像數(shù)據(jù)，時(shí)間范圍為2019年9月1日—2021年4月11日。選用POD精密定軌星歷數(shù)據(jù)去除軌道相位，采用SRTM-DEM去除地形相位并進(jìn)行坐標(biāo)系轉(zhuǎn)換。表1列舉了Sentinel-1A影像數(shù)據(jù)的基本參數(shù)。

圖2 Sentinel-1A影像位置Fig.2 Sentinel-1A image position

表1 Sentinel-1A數(shù)據(jù)參數(shù)Tab.1 Sentinel-1A parameter

2.2 PS點(diǎn)的選取

首先對(duì)20景SAR數(shù)據(jù)進(jìn)行配準(zhǔn)、干涉、去平處理并計(jì)算出數(shù)據(jù)的時(shí)序振幅標(biāo)準(zhǔn)差與平均值。先采用相干系數(shù)法識(shí)別出具有高相干系數(shù)的像元作為PS點(diǎn)選取目標(biāo)，考慮到研究區(qū)為城鎮(zhèn)，建筑物較多，故設(shè)置相干性閾值為0.75； 然后使用振幅離差指數(shù)進(jìn)一步篩選出具有穩(wěn)定相位的PS點(diǎn)，根據(jù)Ferretti等[9]的研究，通常振幅離差指數(shù)只考慮小于0.25的像素點(diǎn)，故設(shè)置振幅離差指數(shù)閾值為0.25； 最后基于形變速率來確定PS點(diǎn)，因?yàn)楸疚倪x取的PS點(diǎn)需要作為SBAS-InSAR軌道精煉時(shí)的GCP點(diǎn)，而GCP點(diǎn)的基本選取原則為遠(yuǎn)離形變區(qū)域，故本文將形變速率區(qū)間設(shè)置為[-0.1，0.1] mm/a。基于以上3個(gè)限制條件，共選出相位穩(wěn)定的PS點(diǎn)36個(gè)，如圖3所示，且將PS點(diǎn)導(dǎo)入地圖對(duì)比查看后發(fā)現(xiàn)大部分PS點(diǎn)都位于穩(wěn)定的建筑物(道路、屋頂?shù)?上。

圖3 基于3個(gè)閾值得到的PS點(diǎn)Fig.3 PS points based on the three thresholds

2.3 SBAS-InSAR處理

將得到的PS點(diǎn)由地理坐標(biāo)轉(zhuǎn)為雷達(dá)斜距坐標(biāo)，然后將其作為GCP點(diǎn)與干涉處理后的像對(duì)進(jìn)行軌道精煉，估算和去除殘余的恒定相位和解纏后還存在的相位坡道； 接著通過2次反演估算形變速率和殘余地形，并去除大氣相位； 最后獲取了研究區(qū)地表形變速率(圖4)和時(shí)序形變結(jié)果(圖5)。可以看出： 區(qū)內(nèi)平均地表形變速率在-6.5～5.7 mm/a之間，年平均最大沉降速率為-6.52 mm，說明該地變形程度輕微，且監(jiān)測(cè)到的形變主要集中在塌陷區(qū)邊界及道路周圍，與實(shí)地地質(zhì)調(diào)查的結(jié)果相符合。除開植被覆蓋區(qū)域出現(xiàn)失相干外，塌陷區(qū)內(nèi)大部分農(nóng)田沒有監(jiān)測(cè)到相關(guān)形變是由于春耕時(shí)節(jié)居民對(duì)農(nóng)田進(jìn)行翻土種地等操作，農(nóng)田里長(zhǎng)出了新作物，所以也導(dǎo)致了該區(qū)域的部分失相干情況。以2019年9月1日為第一時(shí)相，其余各期影像相對(duì)于第一時(shí)相的形變?nèi)鐖D5所示。

圖4 SBAS-InSAR形變速率圖Fig.4 Deformation rate diagram of SBAS-InSAR

圖5 研究區(qū)地表形變時(shí)間序列

3 結(jié)果評(píng)述

為了驗(yàn)證該方法的優(yōu)越性，本文使用相同時(shí)間段的SAR影像數(shù)據(jù)，基于選擇GCP點(diǎn)的重要標(biāo)準(zhǔn)[12]： ①?zèng)]有殘余地形條紋； ②沒有形變條紋，遠(yuǎn)離形變區(qū)域； ③沒有相位躍變； ④至少20～30個(gè)點(diǎn)。通過對(duì)比濾波后的干涉圖和解纏后的干涉圖，在沒有相位躍變和殘余地形相位的區(qū)域，人工手動(dòng)選擇了30個(gè)控制點(diǎn)，得到研究區(qū)的形變速率。并收集了研究區(qū)內(nèi)設(shè)置的3處北斗位移監(jiān)測(cè)站點(diǎn)的監(jiān)測(cè)結(jié)果(選取時(shí)間分別為2020年11月6日，2020年12月12日，2021年1月17日，北斗監(jiān)測(cè)結(jié)果的時(shí)間跨度為2020年10月28日—2021年11月25日，與本文InSAR時(shí)間序列相對(duì)一致)，然后在2種SBAS形變結(jié)果中選取了監(jiān)測(cè)站點(diǎn)周圍不超過50 m的3個(gè)變形點(diǎn)，計(jì)算這3個(gè)變形點(diǎn)在同一位移監(jiān)測(cè)時(shí)間內(nèi)的相對(duì)垂直向形變量，并與位移監(jiān)測(cè)數(shù)據(jù)進(jìn)行比較，比較結(jié)果如表2，圖中點(diǎn)號(hào)為北斗監(jiān)測(cè)的站點(diǎn)名。

表2 2種方法與北斗監(jiān)測(cè)點(diǎn)比較情況Tab.2 Comparison between the two methods and the Beidou monitoring station (mm)

由表2看出，2種方法中，結(jié)合了PS點(diǎn)的SBAS-InSAR方法的監(jiān)測(cè)結(jié)果比人工選取GCP點(diǎn)的傳統(tǒng)SBAS-InSAR方法更接近位移監(jiān)測(cè)數(shù)據(jù)，且2種監(jiān)測(cè)結(jié)果與位移監(jiān)測(cè)數(shù)據(jù)均有一定的誤差。通過查閱文獻(xiàn)再結(jié)合實(shí)地考察，分析造成誤差的原因主要有以下3點(diǎn)： ①由于客觀條件的限制，研究區(qū)內(nèi)北斗位移監(jiān)測(cè)站點(diǎn)的設(shè)置與所選擇的InSAR形變點(diǎn)的位置并不是完全重合的，故選擇了在監(jiān)測(cè)站點(diǎn)50 m以內(nèi)的3個(gè)鄰近InSAR形變點(diǎn)來進(jìn)行對(duì)比分析； ②InSAR形變結(jié)果的分析處理過程是通過像元而不是單個(gè)點(diǎn)來處理的，像元覆蓋了一定的范圍，所以得到的結(jié)果可能會(huì)受到周圍臨近像元及整體區(qū)域上的影響而造成誤差； ③除了InSAR技術(shù)本身存在的誤差源，在時(shí)序InSAR中還存在雷達(dá)影像對(duì)精度的影響、干涉相位對(duì)函數(shù)模型觀測(cè)量與未知量之間函數(shù)關(guān)系的影響[13]。InSAR的相位組成可通過函數(shù)表現(xiàn)為以下形式：

(10)

式中：ε為失相干隨機(jī)噪聲，E{φ}=0。基于時(shí)序InSAR技術(shù)，可以將多幅干涉圖組成其函數(shù)模型觀測(cè)量，增加多余觀測(cè)，從而進(jìn)行誤差的解算與分離。

4 結(jié)論

1)在進(jìn)行SBAS-InSAR處理時(shí)，為了使選擇的GCP點(diǎn)不引起較大誤差從而影響到形變結(jié)果，本文通過設(shè)置相干系數(shù)的閾值、振幅離差指數(shù)的閾值以及地表形變速率的閾值選取了穩(wěn)定的PS點(diǎn)，然后將其作為GCP點(diǎn)進(jìn)行軌道精煉及后續(xù)的反演處理得到地表形變信息。并將本文方法與采用人工選取GCP點(diǎn)的方法所得結(jié)果和北斗位移監(jiān)測(cè)結(jié)果進(jìn)行對(duì)比。由本文方法得到的形變結(jié)果比通過人工選擇GCP的SBAS方法得到的形變結(jié)果更為精準(zhǔn)，能更加有效地分析研究區(qū)地表時(shí)序形變。在實(shí)例研究中，研究區(qū)內(nèi)平均地表形變速率在-6.52～5.66 mm/a之間，年平均最大形變速率為-6.52 mm，該地整體變形程度輕微，變形主要集中在塌陷區(qū)邊界及道路周圍。

2)由于本文方法所選取的PS點(diǎn)大部分都位于建筑物上，當(dāng)建筑物沒有發(fā)生豎直方向上的變形時(shí)，也會(huì)被判定為PS點(diǎn)。那么基于這樣的PS點(diǎn)得到的時(shí)序形變結(jié)果也可能會(huì)出現(xiàn)較大誤差。因此，后續(xù)的研究可以通過將升降軌數(shù)據(jù)進(jìn)行融合得到三維的形變量，再進(jìn)一步解算水平與豎直方向的形變，通過共同判斷x和y方向上的位移，這樣得到的永久散射體就會(huì)更為可靠。