基于改進(jìn)PSInSAR技術(shù)的非城區(qū)地表形變監(jiān)測(cè)

張紅峰 劉 瀛,2

1 山西省第五地質(zhì)工程勘察院，山西省臨汾市廣宣街26號(hào)，041000 2 華北理工大學(xué)礦業(yè)工程學(xué)院，河北省唐山市建設(shè)南路80號(hào)，063210

PSInSAR技術(shù)通過分析SAR影像序列中的時(shí)間穩(wěn)定永久散射體(permanent scatterer,PS)進(jìn)行地表形變監(jiān)測(cè)，廣泛應(yīng)用于城區(qū)、滑坡等形變監(jiān)測(cè)及高程修正[1-2]領(lǐng)域。但非城區(qū)中稀疏的穩(wěn)定散射體導(dǎo)致在地表監(jiān)測(cè)過程中PS點(diǎn)分布密度不足，難以精確估計(jì)速率及殘差等信息，從而使PSInSAR技術(shù)在非城區(qū)地形監(jiān)測(cè)應(yīng)用中受限[3]。為提高非城區(qū)PS點(diǎn)的分布密度，相鄰影像連續(xù)干涉對(duì)[4-5]、小尺寸PS點(diǎn)結(jié)合角反射器優(yōu)化PS點(diǎn)分布[6]、殘差估計(jì)與解纏繞模型優(yōu)化提高空間相干PS點(diǎn)識(shí)別率[7]、多尺度頻率自適應(yīng)相干估計(jì)提高PS點(diǎn)識(shí)別[8]等提高PS點(diǎn)分布密度的改進(jìn)算法被相繼提出。

在已有研究基礎(chǔ)上，為解決非城區(qū)PS點(diǎn)分布不足而導(dǎo)致PSInSAR技術(shù)在非城區(qū)地表形變監(jiān)測(cè)誤差較大的問題，提出基于分時(shí)散射體(partial time scatterer, PTS)的改進(jìn)PSInSAR算法。首先對(duì)SAR影像進(jìn)行基于改進(jìn)經(jīng)驗(yàn)?zāi)B(tài)分解(empirical mode decomposition, EMD)算法的濾波降噪，然后采用雙閾值與可信概率估計(jì)進(jìn)行PTS提取，最后通過參數(shù)差分迭代估計(jì)法對(duì)PTS目標(biāo)進(jìn)行相位分離和形變速率估計(jì)，從而得到非城區(qū)監(jiān)測(cè)區(qū)的地表形變。

1 邊緣保持平滑濾波

在經(jīng)典PSInSAR算法的SAR影像干涉對(duì)(i,k)中，像元P的鄰域像素具有平穩(wěn)隨機(jī)過程特性[9]，因此其相干系數(shù)可表示為：

(1)

SAR影像中通常含有較多的干擾噪聲，因此在PTS目標(biāo)提取前需對(duì)SAR影像進(jìn)行濾波降噪。將傳統(tǒng)EMD算法直接應(yīng)用于SAR影像時(shí)，易造成對(duì)SAR影像邊緣區(qū)域像素值的過度濾波[5]，并且其模式混疊及以經(jīng)驗(yàn)為依據(jù)的分量識(shí)別會(huì)使重建信號(hào)的準(zhǔn)確度較低。為此，在傳統(tǒng)算法行列分解的基礎(chǔ)上，引入45°和135°兩個(gè)分解方向，同時(shí)為優(yōu)化IMF分量分界時(shí)相關(guān)系數(shù)法的復(fù)雜性和單一指標(biāo)的不確定性，本文聯(lián)合SAR影像的細(xì)節(jié)信息與逼近信息，采用基于均方根R和平滑度r估計(jì)的IMF分界K值自適應(yīng)定位方法。假設(shè)干涉圖經(jīng)過改進(jìn)EMD處理后的分量為dj(x)，等權(quán)濾波并重構(gòu)后可得：

(2)

式中，z和ω(x)分別為平滑濾波器的系數(shù)及窗口值。R和r的表達(dá)式為：

(3)

(4)

式中，x(t)為含噪SAR干涉對(duì)影像序列。歸一化R和r，并由變異系數(shù)法進(jìn)行賦權(quán)處理，即

(5)

式中，μ和σ為相應(yīng)的均值與標(biāo)準(zhǔn)差。由此可得，IMF分界的復(fù)合指標(biāo)為Fo=WR×PR+Wr×Pr，相應(yīng)的分界閾值為：

Tk-1=|Fo(k-1)/Fo(k)|

(6)

則濾波時(shí)的約束權(quán)系數(shù)ω(x)為：

(7)

2 PTS目標(biāo)提取

對(duì)SAR影像干涉圖進(jìn)行邊緣保持濾波后，首先由離差指數(shù)和相干系數(shù)對(duì)PTS目標(biāo)點(diǎn)進(jìn)行聯(lián)合初識(shí)，獲得PTS的候選目標(biāo)[9]，以縮小后續(xù)目標(biāo)點(diǎn)的搜索范圍，然后由目標(biāo)點(diǎn)的相位信息推導(dǎo)其可信概率作為PTS可靠篩選的依據(jù)。

(8)

(9)

(10)

式中，Δω為觀測(cè)值與擬合值的相位差，即

(11)

進(jìn)而可得候選目標(biāo)為可靠PTS的概率為：

Pg=1-[(1-α)pR(rg)/p(rg)]

(12)

式中，p(·)與pR(·)分別為PTS像元與非PTS像元的概率函數(shù)，α為調(diào)節(jié)系數(shù)。

3 基于PTS的PSInSAR監(jiān)測(cè)

根據(jù)大氣自相關(guān)特性對(duì)由提取的PTS目標(biāo)構(gòu)建的Delaunay三角網(wǎng)進(jìn)行二次差分，并建立三角網(wǎng)連接邊上像元的觀測(cè)方程[8]。在提取PTS目標(biāo)點(diǎn)時(shí)，由于PTS目標(biāo)的季節(jié)、氣象等的分時(shí)穩(wěn)定性是主要參考因素，傳統(tǒng)的基于共同主影像PSInSAR算法的參數(shù)解算過程在解析PTS貢獻(xiàn)時(shí)較為困難，因此借助連續(xù)干涉對(duì)思想[4]，將干涉圖聚類為N-1個(gè)具有不同時(shí)相特性的組合(圖1)，并以相干度作為組內(nèi)各影像貢獻(xiàn)度的衡量權(quán)重，以削弱低相干影像的貢獻(xiàn)。

圖1 影像干涉對(duì)時(shí)間基線及其分組Fig.1 Image interference pair time baseline and its grouping

基于以上分析，以提取的PTS目標(biāo)對(duì)時(shí)序εp進(jìn)行解算，其解算模型為：

(13)

(14)

采用經(jīng)典PSInSAR算法進(jìn)行初始值解纏，通過Delaunay迭代求解高程修正與形變速率，即可解算出非線性形變相位[2]。

4 實(shí)驗(yàn)分析

為驗(yàn)證本文算法在非城區(qū)形變監(jiān)測(cè)中的有效性，以陜西靖邊縣周邊區(qū)域作為測(cè)試實(shí)驗(yàn)區(qū)(圖2)，實(shí)驗(yàn)數(shù)據(jù)為24景C波段干涉寬模式的SLC影像數(shù)據(jù)，采用ESA精密軌道數(shù)據(jù)處理配準(zhǔn)及相關(guān)相位去除等操作。

圖2 實(shí)驗(yàn)區(qū)光學(xué)影像Fig.2 Optical image of the test area

4.1 PTS目標(biāo)點(diǎn)提取實(shí)驗(yàn)

采用本文算法和傳統(tǒng)PSInSAR算法分別提取PTS點(diǎn)和PS點(diǎn)，實(shí)驗(yàn)結(jié)果如圖3所示。從圖中可以看出，在A區(qū)域，建筑物等可被較好地識(shí)別為PTS點(diǎn)和PS點(diǎn)；在B區(qū)域，植被較為稀少，PS點(diǎn)因失相干性而較少，密度較低，但本文算法仍可識(shí)別出較多的PTS點(diǎn)。對(duì)PTS點(diǎn)和PS點(diǎn)進(jìn)行編號(hào)和經(jīng)緯度重疊分析后發(fā)現(xiàn)，大部分傳統(tǒng)PS點(diǎn)均被標(biāo)記為PTS點(diǎn)，即這些同名的PTS點(diǎn)與PS點(diǎn)具有相同的空間分布特性，說明本文的提取算法具有有效性。

4.2 實(shí)測(cè)數(shù)據(jù)對(duì)比分析

為驗(yàn)證本文算法監(jiān)測(cè)沉降的精度和可靠性，以實(shí)驗(yàn)區(qū)內(nèi)2016年二等水準(zhǔn)資料為基礎(chǔ)，復(fù)測(cè)2016年布設(shè)的水準(zhǔn)路線，以復(fù)測(cè)水準(zhǔn)實(shí)際值作為精度和可靠性實(shí)驗(yàn)數(shù)據(jù)的基準(zhǔn)參數(shù)，并將水準(zhǔn)監(jiān)測(cè)的垂線形變方向投影到SAR衛(wèi)星LOS方向，選取與實(shí)驗(yàn)時(shí)間相同的水準(zhǔn)歷史資料對(duì)復(fù)測(cè)數(shù)據(jù)進(jìn)行精度評(píng)價(jià)。

DSInSAR地形形變監(jiān)測(cè)技術(shù)引入了雷達(dá)最新的分布式識(shí)別和濾波技術(shù)，可將裸土、瀝青路面，特別是新圍墾區(qū)等區(qū)域識(shí)別為DS點(diǎn)，因此將DSInSAR技術(shù)引入到實(shí)驗(yàn)中作為一種性能比較算法。提取監(jiān)測(cè)區(qū)內(nèi)水準(zhǔn)點(diǎn)對(duì)應(yīng)的形變數(shù)據(jù)點(diǎn)，選用平均誤差?和中誤差m作為評(píng)價(jià)指標(biāo)[10]。

表1(單位mm/a)為實(shí)驗(yàn)結(jié)果，從表中可以看出，在A區(qū)域，由于城區(qū)內(nèi)永久相干目標(biāo)較多，3種算法均獲得了較好的誤差精度；在B區(qū)域，PS點(diǎn)較少，形變監(jiān)測(cè)精度的提高需依賴PTS點(diǎn)或DS點(diǎn)，因此本文算法和DSInSAR技術(shù)算法具有明顯的優(yōu)勢(shì)，但本文算法更優(yōu)，原因?yàn)楸疚乃惴ǖ腜TS點(diǎn)可充分利用測(cè)量點(diǎn)在不同時(shí)刻的貢獻(xiàn)優(yōu)勢(shì)；在C區(qū)域，PTS點(diǎn)極少，瀝青路面、新圍墾區(qū)等具有雷達(dá)分布特性的DS點(diǎn)也較少，因此3種算法的監(jiān)測(cè)誤差均較大，而本文算法通過分時(shí)分析，在測(cè)量點(diǎn)較少的情況下充分利用各測(cè)量點(diǎn)在特征最優(yōu)時(shí)刻的監(jiān)測(cè)貢獻(xiàn)度，避免監(jiān)測(cè)點(diǎn)在特征較差時(shí)段對(duì)精度的不利影響，從而有效提高了形變監(jiān)測(cè)的精度。

圖3 PTS與PS提取實(shí)驗(yàn)對(duì)比Fig.3 Comparison of PTS and PS extraction experiments

表1 實(shí)驗(yàn)區(qū)內(nèi)各算法的監(jiān)測(cè)精度

圖4為B區(qū)域水準(zhǔn)測(cè)量結(jié)果與各算法監(jiān)測(cè)結(jié)果的互差直方圖，進(jìn)一步驗(yàn)證了本文算法的監(jiān)測(cè)精度較好，結(jié)果具有可靠性。

圖4 實(shí)測(cè)值與各算法監(jiān)測(cè)結(jié)果互差直方圖Fig.4 Histogram of mutual difference between measurement and monitoring results of each algorithm

5 結(jié) 語

為解決PSInSAR技術(shù)在非城區(qū)監(jiān)測(cè)誤差較大的問題，提出基于PTS的改進(jìn)PSInSAR算法。該算法通過改進(jìn)的EMD對(duì)SAR影像進(jìn)行邊緣保持平滑濾波降噪，通過可信概率聯(lián)合提取PTS目標(biāo)；然后通過迭代估計(jì)PTS相位和形變速率，從而實(shí)現(xiàn)非城區(qū)地表形變的準(zhǔn)確監(jiān)測(cè)。實(shí)驗(yàn)結(jié)果表明，本文算法在非城區(qū)可獲得比PSInSAR技術(shù)和DSInSAR技術(shù)更優(yōu)的監(jiān)測(cè)精度和可靠性，算法具有有效性。