熊文秀，馮光財(cái)，李志偉，杜亞男，李 寧

中南大學(xué)地球科學(xué)與信息物理學(xué)院，湖南 長(zhǎng)沙410083

1 引 言

合成孔徑雷達(dá)遙感差分干涉測(cè)量技術(shù)（differential interferometric synthetic aperture radar，DInSAR）作為一種新興的對(duì)地觀測(cè)技術(shù)手段，可以在不需要地面控制點(diǎn)的情況下監(jiān)測(cè)厘米級(jí)甚至毫米級(jí)精度的地表形變［1］，相較于傳統(tǒng)的形變監(jiān)測(cè)方法（如精密水準(zhǔn)、GPS監(jiān)測(cè)等）具有全天時(shí)、全天候、覆蓋面廣、省時(shí)省力的特點(diǎn)［2］。然而由于大氣相位延遲，地表的時(shí)空去相干以及DEM誤差等因素的影響，D－InSAR技術(shù)的應(yīng)用和推廣也受到一定的限制［3］。因此，一些學(xué)者相繼提出了時(shí)序InSAR技術(shù)，如永久散射體技術(shù)［4］（permanent scatterer technology，PS）、小基線集技術(shù)［5］（small baseline subset technology，SBAS），StaMPS技術(shù)［6］（Stanford method of permanent scatterer technology）等，利用SAR影像上時(shí)序穩(wěn)定相位的連續(xù)觀測(cè)，通過(guò)不同相位解纏模型和參數(shù)估計(jì)算法獲得時(shí)序上的形變序列。高質(zhì)量時(shí)序穩(wěn)定相位的選取克服了傳統(tǒng)D－InSAR中時(shí)空失相干現(xiàn)象并減弱了大氣延遲的影響。其中，SBAS技術(shù)相較于永久散射體技術(shù)而言，具有降低噪聲水平、監(jiān)測(cè)大范圍形變的特點(diǎn)，已經(jīng)被廣泛應(yīng)用于國(guó)內(nèi)外地表形變監(jiān)測(cè)［7－10］。SBAS技術(shù)利用時(shí)空基線限制條件將干涉對(duì)分成若干個(gè)小基線集，再利用奇異值分解將各個(gè)小基線集聯(lián)系以來(lái)，在選取的高質(zhì)量點(diǎn)上解算出時(shí)序形變［5］，為了區(qū)別于傳統(tǒng)的高相干點(diǎn)，本文將時(shí)序相位穩(wěn)定的點(diǎn)稱為高質(zhì)量點(diǎn)。

高質(zhì)量點(diǎn)的選取對(duì)于SBAS技術(shù)最終呈現(xiàn)的形變信息至關(guān)重要，而目前已有的SBAS高質(zhì)量點(diǎn)選取算法仍然存在一定的缺陷。最初，文獻(xiàn)［5］利用相干系數(shù)法來(lái)選取高相干點(diǎn)。該方法雖然計(jì)算簡(jiǎn)便，但是在相干系數(shù)計(jì)算的過(guò)程中，窗口太大時(shí)容易產(chǎn)生旁瓣效應(yīng)［11］，窗口太小時(shí)，不滿足零均值的條件，常常造成相干性的高估［12］，并且干涉圖上存在的形變條紋也會(huì)造成相干性的低估，因此相干系數(shù)作為衡量像元相位時(shí)序穩(wěn)定性的準(zhǔn)則是有偏的［12］。文獻(xiàn)［13］在利用條件函數(shù)最優(yōu)解算局部形變時(shí)，同時(shí)也利用條件函數(shù)的值對(duì)通過(guò)相干系數(shù)法預(yù)選出的點(diǎn)進(jìn)行了二次篩選，這種方法不僅保留了相干系數(shù)法的缺點(diǎn)，而且在條件函數(shù)的計(jì)算中沒(méi)有考慮到非線性形變。文獻(xiàn)［14］將小基線集的思想應(yīng)用到Sta MPS中，雖然利用多基線的干涉相位信息進(jìn)行選點(diǎn)，但始終只是對(duì)未解纏的相位進(jìn)行處理，并且對(duì)于軌道誤差較大的影像容易造成誤估。文獻(xiàn)［11］利用多時(shí)相SAR影像對(duì)相干系數(shù)計(jì)算中的偏差進(jìn)行了改正，得到了較好的效果，但是算法過(guò)于復(fù)雜。更重要的是，目前大部分相干系數(shù)法都是利用時(shí)序上相干系數(shù)的均值作為衡量像素相干性的標(biāo)準(zhǔn)［8－10］，沒(méi)有考慮時(shí)序上的統(tǒng)計(jì)信息，因此對(duì)于在一些時(shí)段保持高相干而在其他時(shí)段相干性較低的點(diǎn)很容易產(chǎn)生誤選。

鑒于此，本文探索一種反映相位質(zhì)量的準(zhǔn)則和算法來(lái)進(jìn)行SBAS技術(shù)高質(zhì)量點(diǎn)的選取。根據(jù)干涉相位各個(gè)組成成分的時(shí)空特征，分離得到噪聲相位，最終利用時(shí)序噪聲相位判斷相位的質(zhì)量。

2 顧及時(shí)空特性的SBAS高質(zhì)量點(diǎn)選取算法

SBAS技術(shù)應(yīng)用于解纏后的干涉相位，而干涉相位的各個(gè)組成成分在時(shí)序和空間上表現(xiàn)出的特性有所不同，因此本文從干涉相位在時(shí)空上的特性出發(fā)，探索衡量相位時(shí)序穩(wěn)定性的準(zhǔn)則和空間相關(guān)相位成分估計(jì)算法。

2．1 時(shí)序穩(wěn)定相位判斷準(zhǔn)則

假設(shè)對(duì)于解纏后的差分干涉圖上任一像元x有

式中，i為干涉對(duì)號(hào)；φdiff，i，x為解纏后的差分干涉相位；φdef，i，x為 形 變 相 位；φatm，i，x為 大 氣 相 位；φorbit，i，x為 軌 道 相 位；φtop＿res，i，x為 高 程 殘 余 相 位；φnoi，i，x包含噪聲相位以及數(shù)據(jù)處理誤差。在這些相位成分中，形變相位、大氣相位以及軌道相位在空間上存在較大相關(guān)性，可以利用空間上的低通濾波進(jìn)行估計(jì)，用表示。另外，高程殘余相位與垂直基線具有式（2）所示線性關(guān)系，可以利用最小二乘的方法估計(jì)出［6］

式中，B⊥，i為干涉對(duì)的垂直基線；λ為雷達(dá)波長(zhǎng)；Ri為天線到地面的距離；θi為雷達(dá)波的入射角；Δzx為像元的高程殘余。構(gòu)造如式（3）所示的時(shí)間相干系數(shù)作為判斷像元質(zhì)量好壞的準(zhǔn)則，M為干涉對(duì)數(shù)量；＾φtop＿res，i，x為利用最小二乘準(zhǔn)則估計(jì)的高程殘余相位；γx取值范圍在0到1之間，表示像元質(zhì)量從低到高

空間上相關(guān)的相位可以通過(guò)在一定范圍內(nèi)進(jìn)行空間上的低通濾波獲得，該范圍通常設(shè)定為大氣在空間上的相關(guān)范圍（800 m）［6］。如果空間相關(guān)相位成分估計(jì)不準(zhǔn)確，將引起φnoi，i，x項(xiàng)包含更多的成分，繼而導(dǎo)致時(shí)間相干系數(shù)γx作為判斷像元質(zhì)量好壞的準(zhǔn)則有偏，因此空間相關(guān)相位成分的準(zhǔn)確估計(jì)對(duì)高質(zhì)量點(diǎn)的選取具有重要意義。

2．2 空間相關(guān)相位成分估計(jì)

常用的空間相關(guān)成分的估計(jì)方法為均值濾波，但是均值濾波得到的結(jié)果較模糊，對(duì)邊緣的保持效果較差。文獻(xiàn)［6］利用自適應(yīng)的巴特沃斯濾波器進(jìn)行頻率域的低通濾波，考慮了截止頻率處的光滑響應(yīng)，以及不同頻譜強(qiáng)度處的自適應(yīng)定權(quán)。本文采用考慮同質(zhì)點(diǎn)的Non－Local濾波來(lái)獲得空間相關(guān)相位成分。先利用時(shí)序SAR影像的強(qiáng)度信息進(jìn)行兩個(gè)樣本的KS檢驗(yàn)，確定參考點(diǎn)一定范圍內(nèi)的同質(zhì)點(diǎn)，再利用Non－Local算法對(duì)同質(zhì)點(diǎn)定權(quán)，進(jìn)行空間平均獲得參考點(diǎn)上的空間相關(guān)相位。這樣不僅排除了非同質(zhì)點(diǎn)對(duì)參考點(diǎn)的影響，而且實(shí)現(xiàn)了低通濾波的自適應(yīng)定權(quán)。文獻(xiàn)［15］證明Non－Local算法能夠更好地保存影像中各個(gè)成分的邊緣信息。

2．2．1 兩個(gè)樣本的KS檢驗(yàn)

KS檢驗(yàn)全稱為柯?tīng)柲缏宸颍鼓β宸颍↘olmogorov－Smirnov）檢驗(yàn)，用來(lái)檢驗(yàn)數(shù)據(jù)的分布是否符合某個(gè)已知的分布。兩個(gè)樣本的KS檢驗(yàn)用來(lái)區(qū)分兩個(gè)數(shù)據(jù)是否服從同一分布，易于實(shí)施，沒(méi)有參數(shù)設(shè)置，并且當(dāng)影像數(shù)量不少于8景時(shí)就能得到可靠的檢驗(yàn)結(jié)果［16］。

在雷達(dá)遙感影像上，穩(wěn)定的像元通常對(duì)應(yīng)穩(wěn)定的后向散射特性，而對(duì)于不穩(wěn)定的像元其后向散射特性通常表現(xiàn)得更為隨機(jī)。因此，利用表征雷達(dá)回波后向散射特性的強(qiáng)度值作為KS檢驗(yàn)的觀測(cè)值［16］。

假設(shè)有N張配準(zhǔn)到同一影像空間的SAR影像，對(duì)其中任一像元x有如下從小到大排列的強(qiáng)度觀測(cè)值

則有累積頻率分布函數(shù)

對(duì)于影像上任意兩點(diǎn)x1和x2，為了檢驗(yàn)它們是否是同質(zhì)點(diǎn)，可以建立如下檢驗(yàn)統(tǒng)計(jì)量

不論觀測(cè)值服從什么分布，檢驗(yàn)統(tǒng)計(jì)量DN均近似為KS分布，其累積概率密度函數(shù)為［16，17］

根據(jù)顯著性水平α設(shè)定閾值c，當(dāng)DN≤c時(shí)以置信水平α認(rèn)為點(diǎn)x1與x2是同質(zhì)點(diǎn)。

2．2．2 Non－Local算法

在SBAS技術(shù)中，差分干涉相位包含多個(gè)空間相關(guān)的成分：形變相位、大氣相位、軌道誤差相位，但是它們的相關(guān)范圍有所不同，尤其對(duì)于形變相位，其空間相關(guān)的范圍可大可小。為了獲得準(zhǔn)確的空間相關(guān)相位，除了要考慮參與低通濾波的像元與參考像元是否同質(zhì)外，還要考慮不同空間相關(guān)相位成分的相關(guān)范圍，即對(duì)邊緣的清晰反映。而Non－Local算法，根據(jù)相位空間分布的差異進(jìn)行定權(quán)，可以很好地保留不同成分的空間分布［15］。

假設(shè)像元i為待計(jì)算像元，j為搜索范圍I內(nèi)的任一除參考像元外的其他像元，則對(duì)像元i可以得到經(jīng)Non－Local算法計(jì)算得到的空間相關(guān)相位為［15，18］

式中，φ（j）為像元j的原始相位值；w（i，j）為像元i與像元j之間的權(quán)重，根據(jù)像元之間空間上的相位差異得到

式中，δ為像元i和像元j各自的鄰域；Ga是與鄰域δ大小相同的高斯卷積核（鄰域δ取7×7）。

3 空間相關(guān)相位成分提取方法比較

如第2．1節(jié)中介紹，準(zhǔn)確提取空間相關(guān)相位成分對(duì)于噪聲相位的分離和高質(zhì)量點(diǎn)的選取至關(guān)重要，在本節(jié)中系統(tǒng)比較了不同濾波器對(duì)影像上空間相關(guān)成分的提取能力，分別選擇反距離加權(quán)均值濾波、文獻(xiàn)［6］中的自適應(yīng)巴特沃斯濾波、Non－Local濾波以及考慮同質(zhì)點(diǎn)的Non－Local濾波進(jìn)行模擬試驗(yàn)和真實(shí)數(shù)據(jù)試驗(yàn)。

3．1 模擬試驗(yàn)

圖1（a）模擬了多個(gè)不同范圍的形變，以此來(lái)比較各個(gè)濾波器對(duì)不同空間相關(guān)范圍成分的提取能力，形變漏斗的半徑依次為10、15、20、25、30、40像素，數(shù)值范圍為0～2 cm，并生成20景添加了高斯白噪聲的影像，高斯白噪聲的標(biāo)準(zhǔn)差為0．28 cm。根據(jù)科爾莫格羅夫湍流理論［19］，模擬0均值大氣噪聲，并從20景影像中分別隨機(jī)選取5景加入隨機(jī)的大氣噪聲和線性軌道項(xiàng)，從而使得模擬數(shù)據(jù)更加接近真實(shí)的干涉相位數(shù)據(jù)。試驗(yàn)過(guò)程中第15景影像同時(shí)加入了大氣項(xiàng)和軌道項(xiàng)，其中大氣項(xiàng)如圖1（b）所示，線性軌道項(xiàng)如圖1（c）所示，圖1（d）為加入大氣項(xiàng)、軌道項(xiàng)和噪聲項(xiàng)的模擬影像，分別利用反距離加權(quán)均值濾波、自適應(yīng)的巴特沃斯濾波、Non－Local濾波以及考慮同質(zhì)點(diǎn)的Non－Local濾波對(duì)該影像進(jìn)行低通濾波，選取濾波窗口的半徑為15，即認(rèn)為空間相關(guān)的范圍為15個(gè)像素，得到的結(jié)果如圖1（e）—（h）所示。

從濾波的結(jié)果可以看出，圖1（e）反距離加權(quán)均值濾波的結(jié)果由于受到濾波窗口的影響，影像基本失真，在形變的范圍與數(shù)值上均發(fā)生較大變化；圖1（f）結(jié)果中形變的范圍依稀可見(jiàn)，對(duì)于范圍大的形變基本保留了原有的形變信息，但是對(duì)于小范圍的形變?yōu)V波后的結(jié)果受截止頻率的限制變得模糊；圖1（g）中Non－Local濾波與反距離加權(quán)均值濾波的不同在于，Non－Local濾波利用數(shù)值差異定權(quán)，對(duì)于小范圍的形變依然造成了很大程度上的失真，但是對(duì)大范圍形變邊緣的保存相較于反距離加權(quán)均值濾波更為清晰；圖1（h）中考慮同質(zhì)點(diǎn)的Non－Local濾波不僅能夠較好地保存范圍大的形變，而且可以很好地保存小范圍的形變，基本上不受濾波窗口大小的影響；圖1（i）—（l）顯示了各個(gè)方法濾波結(jié)果與模擬數(shù)據(jù)真值的差值，考慮同質(zhì)點(diǎn)的Non－Local濾波不僅對(duì)形變區(qū)域信息的保存表現(xiàn)出了明顯的優(yōu)勢(shì)，而且對(duì)于大氣項(xiàng)以及軌道項(xiàng)也能夠很好地提取出來(lái)。因此對(duì)于差分干涉圖上各個(gè)空間相關(guān)相位相關(guān)范圍不一致的情況，利用考慮同質(zhì)點(diǎn)的Non－Local濾波能夠更加有效地提取出空間相關(guān)的相位成分。

3．2 真實(shí)數(shù)據(jù)試驗(yàn)

選取1992—1998年覆蓋上海市的24景JERS－1影像中的小部分區(qū)域作為研究對(duì)象（圖2（a）），根據(jù)文獻(xiàn)［20—21］可知，在該區(qū)域內(nèi)沉降比較嚴(yán)重。利用時(shí)序SAR影像的強(qiáng)度圖判斷各個(gè)點(diǎn)的同質(zhì)性（圖2（b）），并選取1994年12月3日與1997年7月28日的兩景影像生成的解纏差分干涉圖進(jìn)行濾波處理，提取空間相關(guān)的相位成分，根據(jù)大氣在空間上的相關(guān)距離可以確定空間相關(guān)范圍為800 m，結(jié)果如圖2（c）—（f）所示。

均值濾波與自適應(yīng)的巴特沃斯濾波得到的結(jié)果均更為平滑，但是形變的形狀發(fā)生了較大的改變；Non－Local濾波與考慮同質(zhì)點(diǎn)的 Non－Local濾波在形變的范圍上較為接近，但是前者在形變邊緣的細(xì)節(jié)保存上沒(méi)有后者的效果好?？紤]同質(zhì)點(diǎn)的Non－Local濾波也會(huì)因?yàn)樵谒阉鞔翱趦?nèi)同質(zhì)點(diǎn)的個(gè)數(shù)較少而不能準(zhǔn)確地平滑掉噪聲信息。

從上述模擬試驗(yàn)以及真實(shí)數(shù)據(jù)試驗(yàn)均可以看出，本文考慮同質(zhì)點(diǎn)的Non－Local濾波相較于其他的低通濾波器而言能夠更好地分離出影像上空間相關(guān)的成分，并且較好地保存了空間相關(guān)成分的邊緣和細(xì)節(jié)信息。

4 試驗(yàn)分析

由上述結(jié)論可知，考慮同質(zhì)點(diǎn)的Non－Local濾波對(duì)空間相關(guān)相位成分的提取具有一定的優(yōu)勢(shì)，因此本文利用該方法以及最小二乘準(zhǔn)則分別分離出干涉相位中的空間相關(guān)相位和高程殘余相位，最終得到噪聲相位計(jì)算時(shí)間相干系數(shù)作為選取高質(zhì)量點(diǎn)的準(zhǔn)則，具體技術(shù)路線如圖3所示。為了比較本文方法與傳統(tǒng)的SBAS技術(shù)選點(diǎn)方法的差別，分別利用本文的方法以及SBAS技術(shù)中常用的相干系數(shù)法對(duì)試驗(yàn)區(qū)進(jìn)行選點(diǎn)，并在所選點(diǎn)的基礎(chǔ)上解算出平均形變速率。對(duì)兩種方法的優(yōu)缺點(diǎn)進(jìn)行了對(duì)比分析，并對(duì)兩種方法所選點(diǎn)的可靠性進(jìn)行了進(jìn)一步的檢核與驗(yàn)證。

4．1 研究區(qū)域介紹

研究選取1992—1998年的JERS－1數(shù)據(jù)，數(shù)據(jù)主要覆蓋上海市的主城區(qū)、浦東新區(qū)以及閔行區(qū)、南匯區(qū)和寶山區(qū)的部分地區(qū)，面積約為1440 km2，如圖4中紅框所示。主城區(qū)、閔行區(qū)和寶山區(qū)的地物類(lèi)型主要表現(xiàn)為密集的住宅區(qū)與高樓，并遍布相對(duì)密集的道路交通網(wǎng)。2000年以前，浦東新區(qū)、南匯區(qū)主要表現(xiàn)為村落與農(nóng)田相交錯(cuò)的地表特征。選擇地物類(lèi)型相對(duì)復(fù)雜的該區(qū)域作為研究區(qū)域可以反映本文算法與傳統(tǒng)算法在不同地物類(lèi)型覆蓋地區(qū)選點(diǎn)上的差別，并且在該區(qū)域發(fā)育有多個(gè)沉降漏斗可以更好地反映本文算法對(duì)于形變區(qū)域的選點(diǎn)能力。

4．2 研究數(shù)據(jù)介紹

JERS－1數(shù)據(jù)采用3：9多視（分辨率約為40 m×40 m），設(shè)置時(shí)間基線閾值為3 a，空間基線閾值為2000 m，最終由24景SAR影像得到62個(gè)干涉對(duì)，組建情況如圖4（b）所示。JERS－1數(shù)據(jù)相干性較好，但是軌道誤差較大，利用多項(xiàng)式擬合的方法加以去除［22］。研究采用90 m分辨率的SRTM數(shù)據(jù)去除地形相位，上海區(qū)域地形平坦，DEM的精度比較有限。上海市屬于海濱城市，受對(duì)流層大氣影響較大，本文采用掩膜的方法處理受大氣影響嚴(yán)重的區(qū)域。

圖1 4種濾波方法模擬數(shù)據(jù)試驗(yàn)比較Fig．1 Comparison of the four filtering methods using the simulated data

圖2 4種濾波方法真實(shí)數(shù)據(jù)試驗(yàn)比較Fig．2 Comparison of the four filtering methods using the real SAR data

圖3 本文方法流程Fi g．3 The flow chart of the proposed algorithm

圖4 研究區(qū)域（a）和干涉對(duì)組建（b）Fig．4 The study area（a）and interferogram network in this study（b）

4．3 試驗(yàn)結(jié)果與分析

圖5（a）用相干系數(shù)閾值法選出高相干點(diǎn)334 072個(gè)（相干系數(shù)閾值為0．7），并得到平均形變速率圖。圖5（b）中先利用相干系數(shù)閾值法剔除水域等完全失相干的地區(qū)（閾值為0．2），然后利用本文的方法進(jìn)行選點(diǎn)（時(shí)間相干系數(shù)閾值為0．7），最終選取318 587個(gè)高相干點(diǎn)。

從圖5的結(jié)果可以看出，兩種方法在主城區(qū)均能選出足夠密度的點(diǎn)，但是在遠(yuǎn)離城區(qū)的郊區(qū)，相干系數(shù)法的塊狀聚集現(xiàn)象更為明顯，而本文的方法選出的點(diǎn)相對(duì)獨(dú)立，不受旁瓣效應(yīng)的影響，分布較為均勻。圖6（a）、圖6（b）對(duì)應(yīng)圖5中的a、b區(qū)域，表現(xiàn)出明顯的沉降特征。利用相干系數(shù)法（左圖）只能提取出零星的信息，本文方法（右圖）在這些區(qū)域選出的點(diǎn)相對(duì)較多并且相對(duì)均勻，可以更詳細(xì)地反映沉降的情況。圖6（c）對(duì)應(yīng)了圖5中的c區(qū)域，由點(diǎn)的分布可以看出相干系數(shù)法（左圖）選出了較大村莊處的高相干點(diǎn)，而對(duì)于范圍小的村莊則受旁瓣效應(yīng)的影響被淹沒(méi)掉，而本文方法（右圖）可以有效地減少周?chē)匚飳?duì)高質(zhì)量點(diǎn)的影響，從而有效地提取出小村莊處的高質(zhì)量點(diǎn)。圖6（d）對(duì)應(yīng)了圖5中的d區(qū)域，具體表現(xiàn)為村落中間夾雜著農(nóng)田。受附近房屋的影響，相干系數(shù)法（左圖）在農(nóng)田處誤選出了高相干點(diǎn)，本文方法（右圖）有效地避免了誤選農(nóng)田處的點(diǎn)。因此本文方法相較于相干系數(shù)法減弱了旁瓣效應(yīng)的影響，提高了非城區(qū)范圍選點(diǎn)的密度和可靠性。然而本文方法相較于相干系數(shù)法的缺陷主要在于利用Non－Local算法定權(quán)進(jìn)行低通濾波時(shí)效率相對(duì)較低。

4．4 精度評(píng)定

為了評(píng)價(jià)本文選點(diǎn)算法的可靠性，分別對(duì)本文的方法進(jìn)行定量以及定性的評(píng)估。根據(jù)SBAS技術(shù)解算得到的形變相位、高程殘余相位、大氣和軌道相位［5］，可以得到各個(gè)點(diǎn)上的噪聲相位，進(jìn)而利用式（3）再次計(jì)算各點(diǎn)的時(shí)間相干系數(shù)。相干系數(shù)法與本文的方法所選點(diǎn)上時(shí)間相干系數(shù)統(tǒng)計(jì)結(jié)果如圖7所示。

從對(duì)時(shí)間相干系數(shù)的統(tǒng)計(jì)可以看出，相干系數(shù)法以及本文方法的結(jié)果均集中于0．7至0．9。通過(guò)計(jì)算可知，相干系數(shù)法選出的點(diǎn)中時(shí)間相干系數(shù)大于0．7的占總數(shù)的63．96%，本文方法選出的點(diǎn)中時(shí)間相干系數(shù)大于0．7的占總數(shù)的70．74%，由此可見(jiàn)，本文利用時(shí)間相干系數(shù)準(zhǔn)則選出的點(diǎn)總體上優(yōu)于相干系數(shù)法選出的點(diǎn)。

為了進(jìn)一步說(shuō)明本文算法的可靠性，從本文的研究區(qū)域中截取一小片區(qū)域（圖8（a））進(jìn)行人工判讀，該區(qū)域位于浦東新區(qū)，以交錯(cuò)的農(nóng)田與村落為地物特征（2000年前）。值得注意的是，圖8（a）中底圖為時(shí)間最近的Google Earth地圖，因此與SAR影像時(shí)間段內(nèi)的地物有所差別，而所選擇進(jìn)行人工判讀的Google Earth影像時(shí)間為2000年11月22日，相較于本文的研究時(shí)段這片區(qū)域地物沒(méi)有太大的變化。隨機(jī)地從該區(qū)域兩種方法已選點(diǎn)中抽取50個(gè)點(diǎn)進(jìn)行人工判讀，找出與所選50個(gè)點(diǎn)對(duì)應(yīng)的地物類(lèi)型，得到的判讀結(jié)果統(tǒng)計(jì)如圖8（b）所示。

根據(jù)SAR影像以及Google Earth影像分辨率情況，將地物類(lèi)型主要分為房屋、道路、農(nóng)田、水面、綠化帶五大類(lèi)，其中房屋、道路等人工地物視為穩(wěn)定的散射體，可以被認(rèn)為是時(shí)序高質(zhì)量點(diǎn)，而農(nóng)田、水面以及綠化帶由于其失相干特性不能被認(rèn)為是時(shí)序高質(zhì)量點(diǎn)。兩種方法中隨機(jī)選出的50個(gè)點(diǎn)沒(méi)有位于水面上的點(diǎn)。從圖8（b）中的統(tǒng)計(jì)結(jié)果中可以看出兩種方法選點(diǎn)的可靠性均為84%，但是由圖5、圖6可以看出相干系數(shù)法選出的點(diǎn)呈現(xiàn)塊狀聚集分布，并主要分布在城鎮(zhèn)以及范圍大的村莊上，而本文方法除了在這些區(qū)域分布有足夠密度的點(diǎn)，在農(nóng)田中的小村莊上依然能夠選出高質(zhì)量點(diǎn)，證明了本文方法對(duì)于農(nóng)田與村落交錯(cuò)分布地區(qū)高質(zhì)量點(diǎn)的選取能力。

5 討論與總結(jié)

本文顧及時(shí)空特性的SBAS高質(zhì)量點(diǎn)選取方法相較于已有相干系數(shù)法具有如下幾點(diǎn)優(yōu)勢(shì)：

（1）相干系數(shù)法與本文的方法所選點(diǎn)可以被歸類(lèi)為3種情況：①同時(shí)被兩種方法選出的點(diǎn)有268 724個(gè)；②被相干系數(shù)法選出，本文方法未選出的點(diǎn)有65 348個(gè)；③被本文方法選出，相干系數(shù)法未選出的點(diǎn)有49 863個(gè)。圖9中的1號(hào)點(diǎn)即為第2種情況，2號(hào)點(diǎn)為第3種情況。相干系數(shù)法根據(jù)各個(gè)干涉對(duì)相干系數(shù)的均值來(lái)選點(diǎn)，對(duì)于在部分干涉對(duì)上低于閾值部分干涉對(duì)上高于閾值但是均值高于閾值的點(diǎn)則會(huì)被選出，1號(hào)點(diǎn)即為這種情況，相干系數(shù)均值為0．702 6，高于閾值0．7，但是1號(hào)點(diǎn)位于靠近岸邊的水面上，不適合作為高質(zhì)量點(diǎn)。本文的方法在時(shí)序上計(jì)算時(shí)間相干系數(shù)作為衡量像元在各個(gè)干涉對(duì)上綜合質(zhì)量的標(biāo)準(zhǔn)，計(jì)算得到1號(hào)點(diǎn)的時(shí)間相干系數(shù)為0．417，有效地將1號(hào)點(diǎn)從選取的高質(zhì)量點(diǎn)集中排除。

（2）從點(diǎn)的分布上來(lái)看，相干系數(shù)法受到相干系數(shù)估計(jì)窗口的影響旁瓣效應(yīng)明顯，選出的高相干點(diǎn)呈現(xiàn)塊狀分布的趨勢(shì)。圖9中的1號(hào)點(diǎn)則是由于受到估計(jì)窗口內(nèi)岸上像元的影響，造成相干性的高估，從而被誤選為高相干點(diǎn)。在農(nóng)田與村落交錯(cuò)的郊區(qū)，受周?chē)r(nóng)田的影響，面積較小的村莊常常會(huì)被淹沒(méi)（如圖9中2號(hào)點(diǎn)相干系數(shù)為0．615，低于閾值）。本文方法能夠很好地緩解這種缺點(diǎn)，可以有效地選出位于小面積村莊上的高質(zhì)量點(diǎn)（時(shí)間相干系數(shù)為0．744 9，高于閾值）。

圖5 本文方法（b）與相干系數(shù)法（a）形變結(jié)果比較Fig．5 Comparison of the deformation results estimated by the proposed approach（b）and coherence approach（a）

圖6 本文方法與相干系數(shù)法局部區(qū)域選點(diǎn)比較Fig．6 Comparison of the extracted points of the proposed approach and coherence approach in local areas

圖7 本文方法與相干系數(shù)法所選點(diǎn)的時(shí)間相干系數(shù)統(tǒng)計(jì)結(jié)果Fig．7 The statistical results of the temporal coherence of the points generated by the proposed approach and coherence approach

圖8 隨機(jī)采樣點(diǎn)分布圖（a）及其人工判讀結(jié)果（b）Fig．8 Distribution of points selected randomly for manual interpretation（left）and the result of manual interpretation（right）

圖9 兩種方法選點(diǎn)情況比較Fig．9 Comparison of the points selected by two methods

（3）對(duì)于沉降區(qū)域，形變引起的條紋會(huì)造成相干性的低估，因此當(dāng)形變條紋對(duì)相干系數(shù)的影響達(dá)到使其低于閾值時(shí)，相干系數(shù)法在這些區(qū)域則不能有效選取出高相干點(diǎn)。本文方法在準(zhǔn)確地分離出形變相位、大氣相位等空間相關(guān)相位的基礎(chǔ)上，能夠有效地利用噪聲相位提取出高質(zhì)量點(diǎn)，從而增加了形變區(qū)域的信息量（如圖6中的a、b區(qū)域）。

SBAS技術(shù)相較于傳統(tǒng)的形變監(jiān)測(cè)技術(shù)而言具有省時(shí)、省力、大面積監(jiān)測(cè)等優(yōu)點(diǎn)，已經(jīng)成為當(dāng)下地面沉降監(jiān)測(cè)的重要手段。如何得到可靠的地表形變數(shù)據(jù)是形變監(jiān)測(cè)中永恒的話題，而高質(zhì)量點(diǎn)的選取則是SBAS技術(shù)中研究這一問(wèn)題的重要一環(huán)。本文利用兩個(gè)樣本的KS檢驗(yàn)方法找出空間上的同質(zhì)點(diǎn)，然后利用Non－Local算法進(jìn)行定權(quán)，從而從解纏相位中去除空間相關(guān)相位，不僅顧及了參與估計(jì)空間相關(guān)相位像元的同質(zhì)性，而且考慮到了各個(gè)空間相關(guān)相位邊緣的保持。繼而從剩余的相位中提取出高程殘余相位，最終剩下噪聲相位以及可能存在的數(shù)據(jù)處理誤差相位，并以該相位作為衡量解纏相位質(zhì)量的基礎(chǔ)。文章采用覆蓋上海地區(qū)的1992—1998年的JERS－1數(shù)據(jù)進(jìn)行試驗(yàn)，并與相干系數(shù)法的結(jié)果比較，最終顯示本文的方法不僅在城區(qū)能夠選取較多的點(diǎn)，而且在遠(yuǎn)離城區(qū)的郊區(qū)依然能夠更大程度上的提取出高質(zhì)量點(diǎn)，人工判讀的結(jié)果表明該方法選取出的高質(zhì)量點(diǎn)具有很高的可靠性。

致謝：感謝JAXA項(xiàng)目（P1390002；P1229002）提供JERS－1數(shù)據(jù)。

［1］ GABRIEL A K，GOLDSTEIN R M，ZEBKER H A．Mapping Small Elevation Changes over Large Areas：Differential Radar Interferometry［J］．Journal of Geophysical Research：Solid Earth （1978—2012），1989，94（B7）：9183－9191．

［2］ LIAO Mingsheng，LIN Hui．Synthetic Aperture Radar Interferometry：Principle and Signal Processing［M］．Beijing：Surveying and Mapping Press，2003：1－26．（廖明生，林琿．雷達(dá)干涉測(cè)量學(xué)：原理與信號(hào)處理基礎(chǔ)［M］．北京：測(cè)繪出版社，2003：1－26．）

［3］ LIU Guoxiang，CH EN Qiang，LUO Xiaojun，et al．Principle and Method of Permanent Scatterer Synthetic Aperture Radar Interferometry［M］．Beijing：Science Press，2012：10－12．（劉國(guó)祥，陳強(qiáng)，羅小軍，等．永久散射體雷達(dá)干涉理論與方法［M］．北京：科學(xué)出版社，2012：10－12．）

［4］ FERRRTTI A，PRA TI C，ROCCA F．Permanent Scatterers in SAR Interferometry［J］．IEEE Transactions on Geoscience and Remote Sensing，2001，39（1）：8－20．

［5］ BERARDINO P，F(xiàn)ORNARO G，LANARI R，et al．A New Algorithm for Surface Deformation Monitoring Based on Small Baseline Differential SAR Interferograms［J］．IEEE Transactions on Geoscience and Remote Sensing，2002，40（11）：2375－2383．

［6］ HOOPER A J．Persistent Scatter Radar Interferometry for Crustal Deformation Studies and Modeling of Volcanic Deformation［D］．California：Stanford University，2006：46－54．

［7］ ZH ANG Yonghong，WU Hong’an，SUN Guangtong．Deformation Model of Time Series Interferometric SAR Technique［J］．Acta Geodaetica et Cartographica Sinica，2012，41（6）：864－869．（張永紅，吳宏安，孫廣通．時(shí)間序列InSAR技術(shù)中的形變模型研究［J］．測(cè)繪學(xué)報(bào)，2012，41（6）：864－869．）

［8］ LI Shanshan，LI Zhiwei，Hu Jun，et al．Investigation of the Seasonal Oscillation of the Permafrost over Qinghai－Tibet Plateau with SBAS－InSAR Algorithm［J］．Chinese Journal of Geophysics，2013，56（5）：1476－1486．（李珊珊，李志偉，胡俊，等．SBAS－InSAR技術(shù)監(jiān)測(cè)青藏高原季節(jié)性凍土形變［J］．地球物理學(xué)報(bào)，2013，56（5）：1476－1486．）

［9］ YIN Hongjie，ZHU Jianjun，LI Zhiwei，et al．Ground Subsidence Monitoring in Mining Area Using DInSAR SBAS Algorithm［J］．Acta Geodaetica et Cartographica Sinica，2011，40（1）：52－58．（尹宏杰，朱建軍，李志偉，等．基于SBAS的礦區(qū)形變監(jiān)測(cè)研究［J］．測(cè)繪學(xué)報(bào)，2011，40（1）：52－58．）

［10］ TIZZANI P，BERARDINO P，CASU F，et al．Surface Deformation of Long Valley Caldera and Mono Basin，California，Investigated with the SBAS－InSAR Approach［J］．Remote Sensing of Environment，2007，108（3）：277－289．

［11］ JIANG Mi，DING Xiaoli，LI Zhiwei．Hybrid Approach for Unbiased Coherence Estimation for Multi－temporal In－SAR［J］．IEEE Transactions on Geoscience and Remote Sensing，2014，52（5）：2459－2473．

［12］ ZEBKER H A，CH EN K．Accurate Estimation of Correlation in InSAR Observations［J］．IEEE Geoscience and Remote Sensing Letters，2005，2（2）：124－127．

［13］ MORA O，LANARI R，MALLORQUI J J，et al．A New Algorithm for Monitoring Localized Deformation Phenomena Based on Small Baseline Differential SAR Interferograms［C］∥Proceedings of the 2002 IEEE International Geoscience and Remote Sensing Symposium，2002：IGARSS＇02．Toronto：IEEE，2002：1237－1239．

［14］ HOOPER A．A Multi－temporal InSAR Method Incorporating both Persistent Scatterer and Small Baseline Approaches［J］．Geophysical Research Letters，2008，35（16）：L16302．

［15］ BUADES A，COLL B，MOREL J M．A Non－Local Algorithm for Image Denoising［C］∥Proceedings of IEEE Computer Society Conference on Computer Vision and Pattern Recognition：CVPR 2005．San Diego：IEEE，2005，2：60－65．

［16］ FERRETTI A，F(xiàn)UMAGALLI A，NOVALI F，et al．A New Algorithm for Processing Interferometric Data－stacks：SqueeSAR［J］．IEEE Transactions on Geoscience and Remote Sensing，2011，49（9）：3460－3470．

［17］ HOU Shumin，LI Yourong，LIU Guanglin．A New Method of Detecting Nonlinear for Time Series Based on KS Test［J］．Journal of Electronics ＆Information Technology，2007，29（4）：808－810．（侯澍旻，李友榮，劉光臨．一種基于KS檢驗(yàn)的時(shí)間序列非線性檢驗(yàn)方法［J］．電子與信息學(xué)報(bào)，2007，29（4）：808－810．）

［18］ BUADES A，COLL B，MOREL J M．Non－Local Image and Movie Denoising［J］．International Journal of Computer Vision，2008，76（2）：123－139．

［19］ HANSSEN R F．Radar Interferometry：Data Interpretation and Error Analysis［M］．Berlin：Springer，2001：197－244．

［20］ DA MOA H－AFARI P．Detecting Ground Settlement of Megacities Using InSAR Techniques［D］．Hong Kong：Hong Kong Polytechnic University，2009：91－173．

［21］ LU Lijun．Method for Land Deformation Detection Based on Time Series SAR Images and Its Application［D］．Wuhan：Wuhan University，2008：99－115．（盧麗君．基于時(shí)序SAR影像的地表形變檢測(cè)方法及其應(yīng)用［D］．武漢：武漢大學(xué)，2008：99－115．）

［22］ ZHOU Hui，F(xiàn)ENG Guangcai，LI Zhiwei，et al．The Fault Slip Distribution of the Myanmar Mw6．8 Earthquake Inferred from InSAR Measurements［J］．Chinese Journal of Geophysics，2013，56（9）：3011－3021．（周輝，馮光財(cái)，李志偉，等．利用InSAR資料反演緬甸Mw6．8地震斷層滑動(dòng)分布［J］．地球物理學(xué)報(bào)，2013，56（9）：3011－3021．）