基于改進PSInSAR技術的非城區地表形變監測

張紅峰 劉 瀛,2

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

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

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

1 邊緣保持平滑濾波

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

(1)

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

(2)

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

(3)

(4)

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

(5)

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

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

(6)

則濾波時的約束權系數ω(x)為：

(7)

2 PTS目標提取

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

(8)

(9)

(10)

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

(11)

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

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

(12)

式中，p(·)與pR(·)分別為PTS像元與非PTS像元的概率函數，α為調節系數。

3 基于PTS的PSInSAR監測

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

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

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

(13)

(14)

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

4 實驗分析

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

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

4.1 PTS目標點提取實驗

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

4.2 實測數據對比分析

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

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

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

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

表1 實驗區內各算法的監測精度

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

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

5 結 語

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