張紅峰 劉 瀛,2

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

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

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

1 邊緣保持平滑濾波

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

(1)

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

(2)

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

(3)

(4)

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

(5)

式中，μ和σ為相應的均值與標準差。由此可得，IMF分界的復合指標為Fo=WR×PR+Wr×Pr，相應的分界閾值為：

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

(6)

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

(7)

2 PTS目標提取

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

(8)

(9)

(10)

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

(11)

進而可得候選目標為可靠PTS的概率為：

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

(12)

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

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

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

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

基于以上分析，以提取的PTS目標對時序εp進行解算，其解算模型為：

(13)

(14)

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

4 實驗分析

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

圖2 實驗區(qū)光學影像Fig.2 Optical image of the test area

4.1 PTS目標點提取實驗

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

4.2 實測數(shù)據(jù)對比分析

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

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

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

圖3 PTS與PS提取實驗對比Fig.3 Comparison of PTS and PS extraction experiments

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

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

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

5 結(jié) 語

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