一種顧及永久散射體的SBAS InSAR時序地表沉降提取方法

王舜瑤，盧小平，劉曉幫，付睢寧

(河南理工大學礦山空間信息技術國家測繪地理信息局重點實驗室，河南 焦作 454003)

地表沉降是制約城市規(guī)劃與發(fā)展的重要因素，常規(guī)的地表沉降手段如水準測量、GPS測量等方法，存在監(jiān)測點過于稀疏、效率低等不足。而新興的合成孔徑雷達差分干涉測量技術(differential interferometric synthetic aperture radar,D-InSAR)能夠大范圍、無接觸、全天候、全天時地獲取地表的毫米級形變信息[1-4]，尤其是永久散射體(persist scatterers InSAR,PS InSAR)和小基線集(small BAseline subsets InSAR,SBAS InSAR)技術[5-7]，前者能有效提取非線性形變并消除大氣相位影響，后者可以顯著增加相干影像對的數(shù)目，解決了時間、空間失相關和大氣效應對InSAR的影響；同時，隨著高分辨率SAR衛(wèi)星的成功運行，D-InSAR技術在地表沉降監(jiān)測中發(fā)揮著越來越重要的作用[8]。文獻[9]利用PS InSAR對香港國際機場沉降的提取，監(jiān)測結(jié)果與水準測量一致。文獻[10]借助PS InSAR及高分辨率X波段Cosmo-SkyMed雷達影像，成功獲取了上海市地鐵隧道沿線及部分高速公路的沉降情況。文獻[11]利用SBAS InSAR技術并結(jié)合GPS和地下水水位數(shù)據(jù)，得到了洛杉磯地表時序形變及含水層貯水系數(shù)等物理參數(shù)。文獻[12]利用SBAS InSAR技術監(jiān)測凍土形變，所得到的結(jié)果與物理變化規(guī)律相吻合。PS/SBAS InSAR各有優(yōu)點，這兩種技術的交叉融合也是當前D-InSAR研究的熱點。文獻[13]結(jié)合SBAS和PS的相干點選取方法，提高了相干點的數(shù)量，并提取了某火山的形變信息。文獻[14]融合PS/SBAS InSAR算法，通過構(gòu)建局部Delaunay三角網(wǎng)，從某一參考點開始對形變速率增量作積分，提取了徐州大區(qū)域地表形變。由于在SBAS InSAR的軌道精煉步驟中，引入穩(wěn)定的地面控制點可以有效估計殘留相位并去除平地效應，從而提高結(jié)果的可靠性。然而，在沒有地面控制點信息的情況下，盲目地選擇地面控制點反而可能引入誤差。

針對該問題，本文提出一種顧及永久散射體的SBAS InSAR處理方法，通過相干性、振幅離差指數(shù)、形變速率3閾值提取穩(wěn)定的永久散射體，并將其作為地面控制點引入SBAS InSAR處理流程，最后通過實例驗證了本文方法的有效性。

1 本文方法

1.1 基于振幅離差和相干性雙閾值的永久散射體探測方法

對于復數(shù)SAR影像，若其實部和虛部分別存在標準差為σs的高斯噪聲，則振幅值A服從Rice分布，即

(1)

式中，I0為Bessel函數(shù)；g為目標反射量，是正實數(shù)。當比值g/σs較小時，即該目標為低信噪比目標時，Rice分布演化為Rayleigh分布。而對于高信噪比目標(通常g/σs>4)，該分布接近于高斯分布。當σs?g時，存在以下關系

(2)

式中，σp為相位標準差；σsI為虛部標準差；σA為時序振幅標準差；mA為時序振幅均值；DA為振幅離差指數(shù)。

相干系數(shù)的數(shù)學表達式為

(3)

式中,γ為相干系數(shù)；M、S分別為構(gòu)成干涉對的SAR影像；*為復數(shù)共軛；m、n分別為局部窗口的大?。籭、j則為像素坐標。

設置振幅離差閾值TD對永久散射體目標進行初次探測，而后設置相干性閾值Tγ對永久散射體目標進行二次探測，濾除大部分失相干嚴重的像素點，從而獲得最終探測結(jié)果。

1.2 基于SBAS InSAR的時序地表沉陷提取

在(t0，t1，…，tn)時間內(nèi)，從獲取的同一地區(qū)N景SAR影像中選取1景圖像作為主圖像，將其他SAR圖像與主圖像配準。假設每期SAR圖像都滿足臨界基線要求，則可生成M組干涉像對，并滿足下列條件

N/2≤M≤N(N-1)/2

(4)

假設前后兩景影像獲取時間分別為tA和tB，這兩景影像差分生成第j景干涉條紋圖，在不考慮大氣延遲相位、殘余地形相位和噪聲相位的情況下，干涉相位可表示為

(5)

式中，φ為干涉相位；j為差分干涉圖的景號，j∈(1,…,M)；λ為電磁波波長；d(tA,x,r)和d(tB,x,r)分別為tA和tB時刻相對于d(t0,x,r)=0的雷達視線方向累積形變量。

相位信息是通過相位解纏獲得，已知像素點的形變量對應的未知相位值為

φ=[φ(t0)φ(t1)…φ(tN-1)]T

(6)

解纏后的相位值為

δφdef=[δφdef(Δt0)δφdef(Δt1)…δφ(ΔtM-1)]T

(7)

將形變累積相位與干涉相位序列的線性關系用矩陣表示為

δφdef=Aφ

(8)

式中，A為M×N階系數(shù)矩陣。若M≥N，通過最小二乘法可得

φ=(ATA)-1ATδφdef

(9)

若M

vj=(φj-φj-1)/(tj-tj-1)

(10)

用式(10)替換式(5)中的相位，可得

(11)

各時段速度在主、從圖像時間間隔上的積分用矩陣可表示為

δφdef=Bv

(12)

式中,B為M×N階矩陣。第j行位于主輔影像獲取時間之間的列，B(j,k)=tk-tk-1，其他情況B(j,k)=0。此時，將奇異值分解應用于矩陣B，就可得到速度矢量v的最小范數(shù)解，獲得相位的線性形變信息。對殘余相位的時間和空間進行濾波，可分離出大氣相位、非線性形變相位，將這兩項形變疊加即可求出精確的地表形變量。

2 試驗及處理流程

試驗數(shù)據(jù)為21景Sentinel 1A衛(wèi)星C波段單視復數(shù)影像，成像時間為2016年11月至2017年11月，覆蓋鄭州全市域。本文選擇鄭州市四環(huán)路以內(nèi)作為研究區(qū)域，并對數(shù)據(jù)進行裁剪，試驗區(qū)域如圖1所示。區(qū)域左上角和右下角的地理坐標分別為(34.847 5°N，113.565 1°E)、(34.659 8°N，113.790 8°E)。區(qū)域內(nèi)有3個三等水準點SZ01、SZ02和SZ03，可作為精度驗證的地面控制點。用于去平地處理的DEM數(shù)據(jù)為美國航空航天局(NASA)發(fā)布的30 m分辨率SRTM(shuttle radar topography mission，SRTM)V4版本[15]。精密衛(wèi)星軌道數(shù)據(jù)由Sentinel-1 Quality Control網(wǎng)站提供。試驗流程如圖2所示。

試驗處理軟件采用瑞士sarmap公司研發(fā)的SARScape，該軟件架構(gòu)于專業(yè)的ENVI遙感圖像處理軟件之上，提供用戶友好的圖形化操作界面，具有專業(yè)雷達圖像處理和分析功能[16]，能夠滿足本次試驗的需求。

2.1 地面控制點的提取

首先對原始SLC數(shù)據(jù)進行多視處理，得到強度信息圖，并根據(jù)試驗區(qū)域?qū)姸刃畔D進行裁剪，去掉不需要的區(qū)域，以提高后續(xù)處理的效率。對裁剪后的數(shù)據(jù)進行PS點提取，其中探測PS點的振幅離差閾值TD設置為3.2。由于研究區(qū)為城市區(qū)域，存在大量永久散射體(如大型建筑物、路燈、鐵路、公路等)，相干性閾值設置過低會導致PS點數(shù)目劇增，大大增加數(shù)據(jù)計算的時間，從而增加后續(xù)篩選工作的難度，通過多次試驗最終確定相干性閾值Tγ為0.97，共探測到3580個PS點，如圖3所示。通過空間域低通、時間域高通濾波的方法減小干涉相位中的大氣延遲相位誤差，獲取更加可靠的PS點形變信息，將形變速率閾值Tv設為±0.1 mm/a，再次對這些PS點進行篩選，獲得33個3閾值探測PS點，并將這些點作為地面控制點引入后續(xù)的SBAS InSAR處理流程。由于SBAS InSAR需要的地面控制點為雷達斜距坐標系，因此需要對3次PS探測結(jié)果進行投影變換，將地理坐標系下的PS點轉(zhuǎn)換為雷達斜距坐標備用，如圖4所示。

2.2 改進的SBAS InSAR處理方法

對裁剪后的數(shù)據(jù)進行SBAS InSAR處理，空間基線、時間基線分別如圖5、圖6所示，其中空間基線閾值為5%，時間基線閾值為100 d。

對干涉處理后得到的干涉圖快視圖進行目視解譯，手動去除相干性較差的干涉相對，避免對形變提取結(jié)果產(chǎn)生影響。在軌道精煉和重去平步驟中，導入之前通過3次閾值篩選生成的地面控制點，代替人工目視選取地面控制點，再通過兩次反演去除殘余地形和大氣相位，經(jīng)過地理編碼后獲得最終形變速率和沉降時間序列圖。

3 試驗結(jié)果及說明

試驗獲得了測區(qū)的形變速率圖(如圖7所示)及相對于第一期SAR圖像(2017年10月26號)的沉降量時間序列圖(如圖8所示)。試驗結(jié)果表明，該地區(qū)出現(xiàn)區(qū)域性形變，形變速率從-61～27 mm/a不等，市中心區(qū)域出現(xiàn)毫米級小幅抬升，南部抬升程度較大，最大值為49 mm，而西南部則呈下降趨勢，最大下沉量為-82 mm，西北部和東部也出現(xiàn)不規(guī)則下沉。

為驗證本文方法的優(yōu)越性，使用相同原始數(shù)據(jù)，通過SBAS InSAR及PS InSAR方法分別獲取測區(qū)時序地表沉降，并將結(jié)果與本文方法所獲結(jié)果進行對比，同時引入3個三等水準測量點位的兩次測量結(jié)果進行驗證，具體數(shù)據(jù)情況見表1。假設水準實測數(shù)據(jù)為真值，并將兩次水準測量結(jié)果與時間最接近的SAR影像時間段內(nèi)(2017年3月31日—2017年9月15日)的沉降情況進行對比，結(jié)果見表2，其中正值為抬升，負值為下降。

表1 三等水準點高程測量結(jié)果及形變量統(tǒng)計

表2 本文方法與SBAS InSAR、PS InSAR和水準測量的形變量對比 mm

由表1和表2可知：①PS InSAR和SBAS InSAR提取的時序地表形變精度均達到亞厘米級，在該試驗區(qū)內(nèi)，PS InSAR結(jié)果比SBAS InSAR結(jié)果更接近真值，這是由于城市中的大型建筑物、樓房、高架橋、路燈等可作為永久散射體的情況較多，能充分發(fā)揮PS InSAR的優(yōu)勢；②本文方法提取結(jié)果比SBAS InSAR和PS InSAR提取結(jié)果更接近真值。

但本文方法與水準測量產(chǎn)生了一定的誤差，可能的原因如下：①SAR數(shù)據(jù)分辨率的限制，數(shù)據(jù)的單個像元實地并不是一個點，覆蓋了一定的范圍，無法精確與水準點重合，其沉降情況受到該區(qū)域整體沉降情況的綜合影響；②兩次水準測量的時間與SAR數(shù)據(jù)的獲取時間存在一定差異，無法精準重合；③水準測量本身可能存在誤差，而且時序數(shù)據(jù)量不足，無法充分驗證地表時序沉降量，若有時序GPS高程點測量結(jié)果支撐驗證，則更能體現(xiàn)InSAR方法的時效性。

4 結(jié) 語

為彌補PS InSAR、SBSA InSAR單一處理方法的不足，本文提出了一種顧及永久散射體的SBAS InSAR時序地表形變提取方法，并通過鄭州市區(qū)的實例研究得出以下結(jié)論：①PS InSAR和SBAS InSAR在城市區(qū)域的時序地表形變提取精度均可達到亞厘米級，而且由于城市中存在較多永久散射體，因此PS InSAR比SBAS InSAR更適用；②本文方法提取結(jié)果優(yōu)于SBAS InSAR和PS InSAR，可以更有效地提取地表時序形變。