王愛春，向茂生，汪丙南

1. 中國科學(xué)院電子學(xué)研究所微波成像技術(shù)國家級重點實驗室,北京100190； 2. 中國科學(xué)院大學(xué),北京100049； 3. 中國資源衛(wèi)星應(yīng)用中心,北京100094

城區(qū)地表形變差分TomoSAR監(jiān)測方法

王愛春1,2,3，向茂生1，汪丙南1

1. 中國科學(xué)院電子學(xué)研究所微波成像技術(shù)國家級重點實驗室,北京100190； 2. 中國科學(xué)院大學(xué),北京100049； 3. 中國資源衛(wèi)星應(yīng)用中心,北京100094

壓縮感知技術(shù)(CS)的差分TomoSAR技術(shù)解決了中高分辨率SAR數(shù)據(jù)在城區(qū)出現(xiàn)的疊掩問題，實現(xiàn)了城區(qū)地表形變信息的重構(gòu)，但是該方法僅利用了目標(biāo)的稀疏特性并沒有考慮目標(biāo)的結(jié)構(gòu)特性，對具有這兩種特性的目標(biāo)進行重構(gòu)時其性能較差。針對這一問題，本文采用聯(lián)合Khatri-Rao子空間和塊壓縮感知(KRS-BCS)，提出了一種差分SAR層析成像方法。該方法依據(jù)目標(biāo)的結(jié)構(gòu)特性和重構(gòu)觀測矩陣具有的Khatri-Rao積性質(zhì)，將稀疏結(jié)構(gòu)目標(biāo)的差分TomoSAR問題轉(zhuǎn)化為Khatri-Rao子空間下的BCS問題，然后對目標(biāo)進行塊稀疏的l1/l2范數(shù)最優(yōu)化求解，最后通過理論分析和仿真試驗對分辨能力和重構(gòu)估計性能進行了定性和定量評價，仿真結(jié)果表明本文所采用的KRS-BCS方法不僅保持了高分辨率的優(yōu)點，而且有效地降低了虛假目標(biāo)出現(xiàn)的概率，大幅度提高了散射點準(zhǔn)確重構(gòu)概率，切實可行地解決了CS方法的不足。應(yīng)用實例研究中，利用34景Envisat衛(wèi)星ASAR時間序列影像對日本千葉縣茂原市城區(qū)進行地表形變監(jiān)測，并以一等水準(zhǔn)點和實時測量的GPS站點觀測數(shù)據(jù)作為參考形變結(jié)果進行驗證，試驗結(jié)果表明采用KRS-BCS方法反演的結(jié)果與參考形變結(jié)果保持了良好的一致性且形變速率整體偏差也較小，實現(xiàn)了較高精度的城區(qū)地表形變估計。

城區(qū)地表形變；差分TomoSAR技術(shù)；Khatri-Rao子空間；塊壓縮感知

隨著城市化和工業(yè)化的快速發(fā)展，大量抽取地下水所造成地面沉降的狀況呈現(xiàn)上升趨勢，已成為世界性的地質(zhì)災(zāi)害，給城市基礎(chǔ)設(shè)施造成了巨大的經(jīng)濟損失。星載合成孔徑雷達干涉測量技術(shù)(interferometric synthetic aperture radar，InSAR)作為新興的空間測量技術(shù)，全天時、全天候、覆蓋廣等優(yōu)點使其在地表形變的監(jiān)測中得到廣泛的關(guān)注；特別是近年來時序多基線InSAR技術(shù)的發(fā)展，大大提高了地表形變監(jiān)測精度和應(yīng)用范圍。時序多基線InSAR技術(shù)開始于20世紀(jì)90年代后期，經(jīng)過十多年的發(fā)展，國內(nèi)外學(xué)者提出了干涉圖堆棧技術(shù)(Stacking)[1]、最小二乘法(least square，LS)[2]、永久散射體法(permanent scatterer，PS)[3-6]、小基線子集法(small baseline subset，SBAS)[7-9]、相干目標(biāo)方法(coherence target，CT)[10-12]、多時相InSAR技術(shù)(multi temporal InSAR，MTI)[13]、多尺度InSAR時間序列分析法(multiscale InSAR time series，MInTS)[14]等多種技術(shù)方法。這些時序多基線InSAR技術(shù)的散射機制均以像素中占主導(dǎo)優(yōu)勢的散射體為主，當(dāng)應(yīng)用于城區(qū)建筑密集區(qū)域時，會產(chǎn)生嚴(yán)重的疊掩效應(yīng)使得單個雷達分辨率單元中包含了來自不同地面目標(biāo)的多個散射目標(biāo)信號，如果不能有效地將這些混疊的多個散射目標(biāo)信號分離開來，就無法獲取準(zhǔn)確的地表形變信息。

針對上述方法對城區(qū)建筑密集區(qū)中高分辨SAR數(shù)據(jù)由疊掩效應(yīng)帶來的無法正確解算地表形變量這一問題，文獻[15]提出了差分SAR層析成像技術(shù)：在SAR層析成像技術(shù)(tomography SAR，TomoSAR)基礎(chǔ)上，通過對時序多基線InSAR數(shù)據(jù)集在層析向-斜距形變速率向形成合成孔徑，獲得目標(biāo)方位向-距離向-層析向-形變速率向的四維空間信息。在該方法提出之初，文獻[15]采用Capon譜估計法實現(xiàn)了地表形變速率估計，但是該方法需要使用多視數(shù)據(jù)估計觀測協(xié)方差矩陣，降低了方位向-距離向的分辨率且旁瓣較高。為了解決上述問題，文獻[16]采用了SVD奇異值分解法，文獻[17] 采用了Backus-Gilber逆問題求解法，但它們在高度向-形變速率向的分辨率較低且還需要處理奇異值矩陣的病態(tài)問題。對于此問題，文獻[18]采用了RELAX算法，該方法雖然提高了層析向-形變速率向的分辨率，但需要進行目標(biāo)信號源數(shù)目的估計，當(dāng)估計的信號源數(shù)目與真實的信號源數(shù)目不一致時會產(chǎn)生嚴(yán)重的誤差。針對上述問題文獻[19—22]采用了壓縮感知技術(shù)(compressive sensing，CS)[23]，該方法不僅保持了方位向-距離向的分辨率而且還實現(xiàn)了層析向-斜距形變速率向的高分辨率聚焦，然而隨著層析向采樣間隔的減小，疊掩效應(yīng)中每個層析向上的散射目標(biāo)將會出現(xiàn)多個空間上呈聚集塊的采樣散射點，使得疊掩效應(yīng)中散射目標(biāo)的多個散射點具有塊結(jié)構(gòu)特性，但是該方法對塊結(jié)構(gòu)信號進行反演時其估計性能較差[24-26]。

為了解決CS方法對塊結(jié)構(gòu)散射目標(biāo)進行反演時的不足，本文利用疊掩效應(yīng)中散射目標(biāo)的塊結(jié)構(gòu)特性與差分TomoSAR中參數(shù)矩陣所具有的Khatri-Rao積性質(zhì)，提出采用聯(lián)合Khatri-Rao子空間與塊壓縮感知[27](Khatri-Rao subspace and block compressive sensing，KRS-BCS)的方法對城區(qū)地表進行反演，通過模擬仿真數(shù)據(jù)和Envisat衛(wèi)星ASAR數(shù)據(jù)與一等水準(zhǔn)點及GPS實時測量數(shù)據(jù)進行對比分析，結(jié)果表明本文所采用的KRS-BCS方法不僅保持了CS方法的優(yōu)點而且還更好地實現(xiàn)了城區(qū)地表形變信息的估計。

1 差分TomoSAR原理

差分TomoSAR成像示意圖如圖1所示。假設(shè)從不同空間和時間位置獲得了同一目標(biāo)區(qū)的M+1景航過SAR單視復(fù)圖像數(shù)據(jù)集，在對時序SAR單視復(fù)圖像數(shù)據(jù)集進行主圖影像選取、圖像序列配準(zhǔn)、去斜、相位補償去除大氣延遲以及基線估計等預(yù)處理步驟后，除主圖像外第m(m=1,2,…,M)航過影像中每個方位向-斜距向分辨率單元的復(fù)數(shù)值gm可看作是相同方位同一斜距下P個散射目標(biāo)信號在層析向-斜距形變速率向平面上的二維聯(lián)合譜，其表達式為

(1)

圖1 差分TomoSAR成像示意圖Fig.1 Differential TomoSAR imaging diagram

(2)

將式(2)代入方程式(1)，有

(3)

(4)

式中，vq(q=1,2,…,Q?M)為斜距形變速率向采樣點，并記斜距形變速率向采樣點間隔ρv=vq-vq-1；δ(·)為狄克拉函數(shù)。

G=Φ·γ+e

(5)

式中

差分TomoSAR技術(shù)主要解決的問題是根據(jù)方程式(5)中的觀測向量G，通過各種方法反演層析向-斜距形變速率向二維平面上目標(biāo)信號的γ值，進而不同方法采用不同準(zhǔn)則依據(jù)其值的位置確定散射點的層析向位置和斜距形變速率向的大小，從而實現(xiàn)差分SAR層析成像的四維成像。

2 KRS-BCS差分TomoSAR方法

2.1 散射目標(biāo)的塊結(jié)構(gòu)特性

從差分TomoSAR成像示意圖(圖1)可知，同一散射目標(biāo)中散射點數(shù)目或塊長度b與層析向采樣間隔ρs存在以下關(guān)系

(6)

由上面分析和式(6)可知，散射目標(biāo)由散射點組成，散射點的多少隨層析向采樣間隔的大小變化，層析向采樣間隔減小，散射目標(biāo)中獲得了更多的散射點，可以更全面地刻畫散射目標(biāo)的整體情況；層析向采樣間隔變大，散射目標(biāo)中獲得了較少的散射點，只能突出散射目標(biāo)的主要特性。塊結(jié)構(gòu)是散射目標(biāo)在壓縮感知中的另一種表述，疊掩中有幾個散射目標(biāo)就對應(yīng)幾個塊結(jié)構(gòu)，若層析向間隔ρs增大到以Δs采樣，則b=1，散射目標(biāo)中只有一個散射點，此時散射目標(biāo)不再具有塊結(jié)構(gòu)特性；另外，由于每一個塊結(jié)構(gòu)中多個散射點來自于同一個散射目標(biāo)，因此可認(rèn)為每一個塊結(jié)構(gòu)或散射目標(biāo)中多個散射點形變速率一樣。

2.2 差分TomoSAR的塊壓縮感知觀測模型

依據(jù)2.1節(jié)散射目標(biāo)的塊結(jié)構(gòu)特性，將散射目標(biāo)中散射點的反射能量矩陣γs按塊長度b分為N塊

γs=

(7)

式中，P=N×b；γT[n]表示第n(n=1,2,…N)個子塊，子塊中的值是同為零或非零的。

(8)

則目標(biāo)信號的γ值為

γ=γs⊙δv=[γT[1]?δT[1]

γT[2]?δT[2]…γT[N]?δT[N]]T

(9)

式中，“⊙”和“?”分別表示為Khatri-Rao積和Kronecker積;設(shè)

Khatri-Rao積定義為

按照與目標(biāo)信號γ相似的塊分塊，可將觀測矩陣Φ分塊定義為

Φ=Φs⊙Φv=

(10)式中

依據(jù)上述塊結(jié)構(gòu)特性，差分SAR層析成像的觀測模型可進一步表示成N個Khatri-Rao子空間下的塊壓縮感知觀測模型，其表達式為

(11)

根據(jù)稀疏重構(gòu)理論，欲在塊壓縮感知框架下對方程式(11)進行精確重構(gòu)，需要滿足塊結(jié)構(gòu)目標(biāo)信號γ的可稀疏表示性和觀測矩陣Φ的限制等距特性，文獻[26—28]已證明其滿足。因此，方程式(11)可通過混合l1/l2范數(shù)最小化求解將其準(zhǔn)確地重構(gòu)出來

(12)

2.3 性能分析

分辨能力和重構(gòu)估計性能是衡量成像方法的兩個重要方面，本文采用以Envisat衛(wèi)星34景ASAR數(shù)據(jù)集為實際參數(shù)(拍攝時間為2006年9月至2010年8月，主圖為2008年8月4日的垂直基線與時間基線分布如圖2所示)的仿真試驗(試驗對象為不受各種誤差影響的理想散射目標(biāo))對其進行評價。

圖2 垂直基線-時間基線分布圖Fig.2 Perpendicular baseline and time baseline distribution

2.3.1 分辨能力

圖3和圖4是塊稀疏數(shù)k=1、2、3和b=4時兩方法在層析向-斜距形變速率向三維平面上和二維平面上的重構(gòu)圖，其中綠色、青色和藍(lán)色的實矩形區(qū)是k=1、2、3時的真實位置，假定同一塊稀疏數(shù)中散射點的散射能量相同且γ1=1.0、γ2=0.4、γ3=0.7，紅色圓形是KRS-BCS方法重構(gòu)的位置，黑色實點是CS方法重構(gòu)的位置。從圖中可以看出，在塊稀疏數(shù)k=1、2、3時，KRS-BCS方法和CS方法都能清晰分辨出不同散射目標(biāo)中4個4.5m間隔的散射點，這充分展現(xiàn)了兩種方法具有相同的高分辨性能，但兩種方法在重構(gòu)估計性能方面差別較大，具體分析見2.3.2節(jié)。

圖3 b=4、不同k時的三維重構(gòu)位置Fig.3 3D reconstruction position of using b=4 and different k

圖4 b=4、不同k時的二維重構(gòu)位置Fig.4 2D reconstruction position of using b=4 and different k

2.3.2 重構(gòu)估計性能

虛假散射目標(biāo)出現(xiàn)概率和準(zhǔn)確重構(gòu)率是重構(gòu)估計性能的主要指標(biāo)，其由觀測矩陣Φ的RIP決定。RIP要求Φ對觀測向量G作用具有保距性，Φ的RIP性能越好保距性越好，保距性越好其重構(gòu)性能越好；由于實際應(yīng)用中RIP很難判斷，常用Φ中原子間的相似性μ代替，μ越小Φ的RIP性能越好；對于CS方法和KRS-BCS方法中的μ和μb其表達式分別為

(13)

(14)

依2.2節(jié)分析可知，隨著層析向采樣間隔的減少，目標(biāo)信號的總體采樣數(shù)目增加，Φ的維數(shù)也相應(yīng)地增加，但Φ的維數(shù)增加相當(dāng)于在各列向量間進行了插值，因此Φ的各列向量間差異隨著層析向采樣間隔的減少也逐漸減小。結(jié)合式(13)可知CS方法中Φ的相似性μ是列與列的相似性，因此μ的值將變大，其RIP性能將變差，故其重構(gòu)估計性能也將變差；而對于KRS-BCS方法中Φ的相似性μb，由式(14)可知μb是塊與塊的相似性，因此塊的大小變化不會對μb產(chǎn)生很大的影響，故其重構(gòu)估計性能穩(wěn)定。

為了更加合理地驗證上述理論分析，本文進行了仿真試驗。通過不同參數(shù)的設(shè)置，對KRS-BCS方法和CS方法的重構(gòu)估計性能從定性和定量兩個方面加以分析。

(1) 定性方面：采用2.3.1節(jié)的仿真，從圖3和圖4可以看出，在塊稀疏數(shù)k=1、2時，KRS-BCS方法和CS方法都沒有虛假散射目標(biāo)，但CS方法在每個散射目標(biāo)的周邊多出2～4個錯誤散射點，而KRS-BCS方法完整地恢復(fù)原始信號；在塊稀疏數(shù)k=3時，KRS-BCS方法還沒有虛假散射目標(biāo)出現(xiàn)，僅在周邊多出1個錯誤散射點，但CS方法出現(xiàn)大量虛假散射目標(biāo)且周邊都有錯誤散射點。

(2) 定量方面：保持2.3.1節(jié)的仿真參數(shù)不變，層析向采樣間隔ρs以18.5 m、8.5 m、4.5 m、2.5 m、1.2 m進行采樣，對應(yīng)子塊長度b=1、2、4、8、16，每個塊的斜距形變速率相同并以0.1 mm/a的分辨率在區(qū)間[-1.5 1.5]內(nèi)隨機產(chǎn)生。分析塊稀疏數(shù)k=1、2、3的情況下，不同子塊長度下的虛假散射目標(biāo)的出現(xiàn)概率和散射目標(biāo)中散射點的準(zhǔn)確重構(gòu)概率，結(jié)果如圖5和圖6所示。從圖5和圖6中可以看出，隨著塊稀疏數(shù)的增加，KRS-BCS方法和CS方法的虛假散射目標(biāo)的出現(xiàn)概率都相應(yīng)增加且散射點位置的準(zhǔn)確重構(gòu)率減小，但CS方法的增加和減小的幅度明顯大于KRS-BCS方法的增加和減小的幅度。另外，CS方法在塊稀疏數(shù)k=3和子塊長度b大于2時虛假散射目標(biāo)出現(xiàn)的概率就已大于0.8且準(zhǔn)確重構(gòu)率已衰減到0，在塊稀疏數(shù)k=1和子塊長度b大于4時虛假散射目標(biāo)出現(xiàn)概率大于0.2且準(zhǔn)確重構(gòu)率也已衰減到0，而KRS-BCS方法在子塊長度b大于4后才開始出現(xiàn)虛假散射目標(biāo)且最大虛假散射目標(biāo)出現(xiàn)概率不大于0.2，準(zhǔn)確重構(gòu)率也開始衰減但其最小值仍大于0.5，這一結(jié)果也表明KRS-BCS方法的性能明顯優(yōu)于CS方法。

圖5 不同k、b時虛假散射目標(biāo)出現(xiàn)概率Fig.5 Probability of false scattering target using different k and b

圖6 不同k、b時散射點準(zhǔn)確重構(gòu)率Fig.6 Reconstruction accurate of scattering point using different k and b

依據(jù)上述理論分析與仿真試驗的結(jié)果，可以得出本文所采用的KRS-BCS方法不僅保持了高分辨率的優(yōu)點，而且有效地降低了虛假目標(biāo)出現(xiàn)概并和大幅度地提高了散射點準(zhǔn)確重構(gòu)概率，切實可行地解決了CS方法的不足。

3 研究區(qū)域地表形變提取與驗證

3.1 研究區(qū)域概括

為了驗證KRS-BCS方法在城區(qū)地表形變監(jiān)測應(yīng)用中的可行性，本文選取日本千葉縣作為試驗研究區(qū)，該區(qū)域位于關(guān)東平原的東南部，具有豐富的天然氣資源和含鹽量很高的地下水，由于天然氣和地下水的過度開采，千葉縣成為日本地表形變最嚴(yán)重的區(qū)域。試驗選取了該區(qū)域2006年9月至2010年8月獲取的34景Envisat衛(wèi)星ASAR影像(2008年8月4日的垂直基線與時間基線分布如圖2所示)，拍攝范圍為35.34°N—35.55°N、140.22°E—140.47°E,如圖7所示。為進行針對性的分析對比，將試驗區(qū)進一步選定為高層建筑物較多的茂原市區(qū)，如圖8所示。圖9為試驗區(qū)的光學(xué)圖像，其中雨點狀標(biāo)記為一等水準(zhǔn)點、圖釘狀標(biāo)記為茂原區(qū)曾根尼之臺綜合公園GPS實時測量站點。

圖7 Envisat ASAR拍攝范圍Fig.7 Scene acquired by Envisat ASAR

圖8 試驗區(qū)的SAR幅度圖像Fig.8 SAR magnitude image of the test area

3.2 研究區(qū)域概括形變結(jié)果與精度驗證

試驗區(qū)在KRS-BCS差分TomoSAR成像處理中斜距形變速率采樣間隔的選取，依據(jù)2.3節(jié)的分析選擇為ρv=0.1 mm/a和塊稀疏的子塊長度選擇為b=4，此時虛假散射目標(biāo)出現(xiàn)概率為0和準(zhǔn)確重構(gòu)概率為1.0，可以充分保證重構(gòu)的準(zhǔn)確度，重構(gòu)結(jié)果如圖10和圖11所示。

圖9 試驗區(qū)的光學(xué)圖像Fig.9 Corresponding Optical image of the test area

圖10 KRS-BCS方法重構(gòu)的測試區(qū)高度與SAR幅度圖像疊加Fig.10 Reconstructing height from KRS-BCS overlapped SAR magnitude image of the test area

圖11 KRS-BCS方法重構(gòu)的測試區(qū)形變速率與SAR幅度圖像疊加Fig.11 Reconstructing deformation velocity from KRS-BCS overlapped SAR magnitude image of the test area

圖10是依據(jù)層析向與高度向轉(zhuǎn)換關(guān)系獲得的重構(gòu)高度與SAR幅度圖像的疊加圖，從中可以獲知測試區(qū)的高度從城區(qū)中心到沿海岸邊呈遞減趨勢，城區(qū)高層建筑物較多其高度大多分布在20～40 m區(qū)間內(nèi)，城區(qū)外以低矮建筑為主其高度分布在0～20 m區(qū)間內(nèi)，整個測試區(qū)的高度分布及趨勢與Google地圖中的三維建筑物高度呈整體一致性。圖11是斜距形變速率向的形變速率圖與SAR幅度圖像的疊加圖，從中可知測試區(qū)的小部分區(qū)域斜距形變速率大于0、呈現(xiàn)抬高趨勢、上升形變速率達20 mm/a，大部分區(qū)域斜距形變速率小于0、呈現(xiàn)下沉趨勢、下沉形變速率可達-15 mm/a，其整體趨勢與千葉縣環(huán)境生活部發(fā)布的形變趨勢基本吻合。

為進一步驗證KRS-BCS方法重構(gòu)的形變速率，本文選取了9個一等水準(zhǔn)點和1個GPS實時測量站點作為檢驗點，對比結(jié)果如圖12和圖13所示。圖12是2007年至2011年期間每年1月測量的9個水準(zhǔn)點形變量、線性擬合后的形變量與KRS-BCS方法重構(gòu)的線性形變量對比圖，為了便于分析將水準(zhǔn)點實測形變量轉(zhuǎn)化到斜距向或斜距形變速率向形變量，從圖12可知，KRS-BCS方法重構(gòu)的形變與水準(zhǔn)點測量的形變都呈下沉趨勢且保持了整體一致性，除3930水準(zhǔn)點外KRS-BCS方法重構(gòu)的形變速率均偏低于水準(zhǔn)點測量的形變速率且偏差也不大，在水準(zhǔn)點3927、3928、MB-12的偏差小于1.5 mm/a，在水準(zhǔn)點14、3230的偏差小于2 mm/a，與水準(zhǔn)點9、MB-15的偏差小于4 mm/a，與水準(zhǔn)點3929、MB-20的偏差在6 mm/a左右，所有點的整體偏差小于3.2 mm/a。

圖12 KRS-BCS方法重構(gòu)形變量與站點實測形變量的對比Fig.12 Comparison of reconstructing deformation from KRS-BCS and measuring deformation from first level point data

圖13是1996年5月至2011年1月間獲得的GPS實測形變量與KRS-BCS方法重構(gòu)的線性形變量對比圖，在分析前先將GPS實時測量數(shù)據(jù)進行濾波、ASAR數(shù)據(jù)成像時間重采樣、線性形變擬合以及高程向形變量到斜距向或斜距形變速率向形變量轉(zhuǎn)化等處理，然后比較GPS濾波后的成像時線性擬合形變量與KRS-BCS方法重構(gòu)的線性形變量，從圖13可知，以1996年5月數(shù)據(jù)采集起始點為參考點的GPS實時測量站點在ASAR數(shù)據(jù)成像期間呈現(xiàn)上升和下沉的交替沉降模式但其整體趨勢呈現(xiàn)下沉，KRS-BCS方法重構(gòu)的形變趨勢與GPS濾波后的成像時線性擬合形變趨勢保持了良好的一致性，都呈現(xiàn)下沉趨勢，且KRS-BCS方法重構(gòu)的線性形變速率-6.7 mm/a與GPS濾波后的成像時線性擬合形變速率-4.3 mm/a相差也較小,為2.4 mm/a。

圖13 KRS-BCS方法重構(gòu)形變量與GPS實測形變量的對比Fig.13 Comparison of reconstructing deformation from KRS-BCS and measuring deformation from GPS tracking data

4 結(jié) 論

本文提出了一種基于KRS-BCS差分TomoSAR技術(shù)的城區(qū)地表形變監(jiān)測方法，利用了城區(qū)地表疊掩效應(yīng)中散射目標(biāo)隨層析向采樣間隔變化而產(chǎn)生的塊結(jié)構(gòu)特性，構(gòu)建了Khatri-Rao子空間下的塊壓縮感知觀測模型，采用混合l1/l2范數(shù)最小化求解實現(xiàn)了反演參數(shù)的重構(gòu)，通過理論分析和仿真試驗說明了CS方法存在的不足，并進一步對比分析了KRS-BCS方法對CS方法不足之處的切實改進：分辨力方面KRS-BCS方法保持了CS方法的高分辨率優(yōu)點，可以實現(xiàn)0.1 mm/a級的地表形變監(jiān)測、重構(gòu)估計性能方面KRS-BCS方法，可以保證在不同塊稀疏和不同子塊長度下，最大虛假散射目標(biāo)出現(xiàn)概率不大于0.2和最小準(zhǔn)確重構(gòu)率仍大于0.5，而CS方法在塊稀疏數(shù)k=1和子塊長度b大于4時虛假散射目標(biāo)出現(xiàn)概率已大于0.2，且準(zhǔn)確重構(gòu)率也已衰減到0；在實例研究應(yīng)用中，利用Envisat衛(wèi)星ASAR數(shù)據(jù)集對日本千葉茂原區(qū)進行地表形變監(jiān)測，并利用一等水準(zhǔn)點和實時測量的GPS站點數(shù)據(jù)進行精度驗證，試驗結(jié)果表明本文采用KRS-BCS方法反演的結(jié)果不僅與一等水準(zhǔn)點及GPS結(jié)果保持了良好的一致性，而且形變速率整體偏差也較小，分別為3.2 mm/a和2.4 mm/a，進一步驗證了該方法的可靠性和精確性。在后續(xù)的研究中嘗試更加復(fù)雜的非線性形變模式和周期形變模式，使本文方法能滿足復(fù)雜條件下的地表形變監(jiān)測需求。

[1] SANDWELL D T,SICHOIX L. Topographic Phase Recovery from Stacked ERS Interferometry and A Low-resolution Digital Elevation Model[J].Journal of Geophysical Research:Atmospheres,2000,205(B12):28211-28222.

[2] USAI S. A New Approach for Long Term Monitoring of Deformations by Differential SAR Interferometry[D]. Delft: Delft University of Technology, 2001.

[3] FERRETTI A, PRATI C, ROCCA F. Permanent Scatterers in SAR Interferometry[J]. IEEE Transactions on Geoscience and Remote Sensing, 2001, 39(1): 8-20.

[4] HOOPER A, ZEBKER H, SEGALL P, et al. A New Method for Measuring Deformation on Volcanoes and Other Natural Terrains Using InSAR Persistent Scatterers[J]. Geophysical Research Letters, 2004, 31(23): L23611.

[5] KAMPES B M.Radar Interferometry: Persistent Scatterer Technique[M]. Netherlands: Springer, 2006.

[6] 李德仁, 廖明生, 王艷. 永久散射體雷達干涉測量技術(shù)[J]. 武漢大學(xué)學(xué)報(信息科學(xué)版), 2004, 29(8): 664-668. LI Deren, LIAO Mingsheng, WANG Yan. Progress of Permanent Scatterer Interferometry[J]. Geomatics and Information Science of Wuhan University, 2004, 29(8): 664-668.

[7] BERARDINO R, FORNARO 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.

[8] 鄧琳, 劉國祥, 張瑞, 等. 多平臺MC-SBAS長時序建模與形變提取方法[J]. 測繪學(xué)報, 2016, 45(2): 213-223. DOI: 10.11947/j.AGCS.2016.20140614. DENG Lin, LIU Guoxiang, ZHANG Rui, et al. A Multi-platform MC-SBAS Method for Extracting Long Gterm Ground Deformation[J]. Acta Geodaetica et Cartographica Sinica, 2016, 45(2): 213-223. DOI: 10.11947/j.AGCS.2016.20140614.

[9] 許文斌, 李志偉, 丁曉利, 等. 利用InSAR短基線技術(shù)估計洛杉磯地區(qū)的地表時序形變和含水層參數(shù)[J]. 地球物理學(xué)報, 2012, 55(2): 452-461. XU Wenbin, LI Zhiwei, DING Xiaoli, et al. Application of Small Baseline Subsets D-InSAR Technology to Estimate the Time Series Land Deformation and Aquifer Storage Coefficients of Los Angeles Area[J]. Chinese Journal of Geophysics, 2012, 55(2): 452-461.

[10] PERISSIN D, WANG Teng. Repeat-pass SAR Interferometry with Partially Coherent Targets[J]. IEEE Transactions on Geoscience and Remote Sensing, 2012, 50(1): 271-280.

[11] 張永紅, 龔文瑜, 張繼賢, 等. 基于SAR干涉點目標(biāo)分析技術(shù)的城市地表形變監(jiān)測[J]. 測繪學(xué)報, 2009, 38(6): 482-487. DOI: 10.3321/j.issn:1001-1595.2009.06.003. ZHANG Yonghong, GONG Wenyu, ZHANG Jixian, et al. Monitoring Urban Subsidence Based on SAR Interferometric Point Target Analysis[J].Acta Geodaetica et Cartographica Sinica, 2009, 38(6): 482-487. DOI: 10.3321/j.issn:1001-1595.2009.06.003.

[12] 王明洲, 李陶, 江利明, 等. 地表形變監(jiān)測的改進相干目標(biāo)法[J]. 測繪學(xué)報, 2016, 45(1): 36-43. DOI: 10.11947/j.AGCS.2016.20140617. WANG Mingzhou, LI Tao, JIANG Liming, et al. An Improved Coherent Targets Technology for Monitoring Surface Deformation[J]. Acta Geodaetica et Cartographica Sinica, 2016, 45(1): 36-43. DOI: 10.11947/j.AGCS.2016.20140617.

[13] HOOPER A. A Multi-temporal InSAR Method Incorporating Both Persistent Scatterer and Small Baseline Approaches[J]. Geophysical Research Letters, 2008, 35(16): L16302.

[14] HETLAND E A, MUSé P, SIMONS M, et al. Multiscale InSAR Time Series(MinTS) Analysis of Surface Deformation[J]. Journal of Geophysical Research: Solid Earth, 2012, 117(B2): B02404.

[15] LOMBARDINI F.Differential Tomography:A New Framework for SAR Interferometry[J]. IEEE Transactions on Geoscience and Remote Sensing, 2005, 43(1): 37-44.

[16] SERAFINO F, SOLDOVIERI F,LOMBARDINI F, et al. Singular Value Decomposition Applied to 4D SAR Imaging[C]∥Proceedings of International Geoscience and Remote Sensing Symposium. Seoul, Korea: IEEE, 2005: 2701-2704.

[17] 任笑真, 楊汝良. 一種基于逆問題的差分干涉SAR層析成像方法[J]. 電子與信息學(xué)報, 2010, 32(3): 582-586. REN Xiaozhen, YANG Ruliang. An Inverse Problem Based Approach for Differential SAR Tomography Imaging[J]. Journal of Electronics & Information Technology, 2010, 32(3): 582-586.

[18] 孫希龍, 余安喜, 董臻, 等. 一種差分SAR層析高分辨成像方法[J]. 電子與信息學(xué)報, 2012, 34(2): 273-278. SUN Xilong, YU Anxi, DONG Zhen, et al. A High Resolution Imaging Method for Differential SAR Tomography[J]. Journal of Electronics & Information Technology, 2012, 34(2): 273-278.

[19] ZHU Xiaoxiang,BAMLER R. Superresolving SAR Tomography for Multidimensional Imaging of Urban Areas: Compressive Sensing-based TomoSAR Inversion[J]. IEEE Signal Processing Magazine, 2014, 31(4): 51-58.

[20] WANG Yuanyuan, ZHU Xiaoxiang, BAMLER R. An Efficient Tomographic Inversion Approach for Urban Mapping Using Meter Resolution SAR Image Stacks[J]. IEEE Geoscience and Remote Sensing Letters, 2014, 11(7): 1250-1254.

[21] 廖明生, 魏戀歡, 汪紫蕓, 等. 壓縮感知在城區(qū)高分辨率SAR層析成像中的應(yīng)用[J]. 雷達學(xué)報, 2015, 4(2): 123-129. LIAO Mingsheng, WEI Lianhuan, WANG Ziyun, et al. Compressive Sensing in High-resolution 3D SAR Tomography of Urban Scenarios[J]. Journal of Radars, 2015, 4(2): 123-129.

[22] SIDDIQUE M A, HAJNSEK I, AERSOSPACE G, et al. Investigating the Combined Use of Differential SAR Tomography and PSI for Spatio-temporal Inversion[C]∥Proceedings of 2015 Joint Urban Remote Sensing Event(JURSE). Lausanne, Switzerland: IEEE, 2015: 1-4.

[23] DONOHO D L. Compressed Sensing[J]. IEEE Transactions on Information Theory, 2006, 52(4): 1289-1306.

[24] ZHU Xiaoxiang,BAMLER R. Super-resolution Power and Robustness of Compressive Sensing for Spectral Estimation with Application to Spaceborne Tomographic SAR[J]. IEEE Transactions on Geoscience and Remote Sensing, 2012, 50(1): 247-258.

[25] 張冰塵, 王萬影, 畢輝, 等. 基于壓縮多信號分類算法的森林區(qū)域極化SAR層析成像[J]. 電子與信息學(xué)報, 2015, 37(3): 625-630. ZHANG Bingchen, WANG Wanying, BI Hui, et al. Polarimetric SAR Tomography for Forested Areas Based on Compressive Multiple Signal Classification[J]. Journal of Electronics & Information Technology, 2015, 37(3): 625-630.

[26] 王愛春, 向茂生. 基于塊壓縮感知的SAR層析成像方法[J]. 雷達學(xué)報, 2016, 5(1): 57-64. WANG Aichun, XIANG Maosheng. SAR Tomography Based on Block Compressive Sensing[J]. Journal of Radars, 2016, 5(1): 57-64.

[27] ELDAR Y C,KUPPINGER P,BOLCSKEI H.Block-sparse Signals: Uncertainty Relations and Efficient Recovery[J]. IEEE Transactions on Signal Processing, 2010, 58(6): 3042-3054.

[28] FU Yuli, LI Haifeng, ZHANG Qiheng, et al. Block-sparse Recovery via Redundant Block OMP[J]. Signal Processing, 2014, 97: 162-171.

(責(zé)任編輯：張艷玲)

Method of Monitoring Urban Area Deformation Based on Differential TomoSAR

WANG Aichun1,2,3,XIANG Maosheng1,WANG Bingnan1

1. National Key Laboratory of Microwave Imaging Technology, Institute of Electronics, Chinese Academy of Sciences, Beijing 100190, China； 2. University of Chinese Academy of Sciences, Beijing 100049, China； 3. China Center for Resources Satellite Data and Application, Beijing 100094, China

While the use of differential TomoSAR based on compressive sensing (CS) makes it possible to solve the layover problem and reconstruct the deformation information of an observed urban area scene acquired by moderate-high resolution SAR satellite, the performance of the reconstruction decreases for a sparse and structural observed scene due to ignoring the structural characteristics of the observed scene. To deal with this issue, the method for differential SAR tomography based on Khatri-Rao subspace and block compressive sensing (KRS-BCS) is proposed. The proposed method changes the reconstruction of the sparse and structural observed scene into a BCS problem under Khatri-Rao subspace, using the structure information of the observed scene and Khatri-Rao product property of the reconstructed observation matrix for differential TomoSAR, such that the KRS-BCS problem is efficiently solved with a block sparse l1/l2norm optimization signal model, and the performance of resolution capability and reconstruction estimation is compared and analyzed qualitatively and quantitatively by the theoretical analysis and the simulation experiments, all of the results show the propose KRS-BCS method practicably overcomes the problems of CS method, as well as, quite maintains the high resolution characteristics, effectively reduces the probability of false scattering target and greatly improves the reconstruction accurate of scattering point. Finally, the application is taking the urban area of the Mobara(in Chiba, Japan) as the test area and using 34 ENVISAT-ASAR images, the accuracy is verifying with the reference deformations derived from first level point data and GPS tracking data, the results show the trend is consistent and the overall deviation is small between reconstruction deformations of the propose KRS-BCS method and the reference deformations, and the accuracy is high in the estimation of the urban area deformation.

urban area deformation; differential TomoSAR; Khatri-Rao subspace; block compressive sensingFoundation support： National Development and Reform Commission Satellite and Application Development Projects(No.[2012]2083)

WANG Aichun(1981—)，male，PhD candidate，engineer，majors in multi-baseline InSAR data process-ing and its applications

XIANG Maosheng

王愛春，向茂生，汪丙南.城區(qū)地表形變差分TomoSAR監(jiān)測方法[J].測繪學(xué)報，2016,45(12)：1413-1422.

10.11947/j.AGCS.2016.20160113. WANG Aichun, XIANG Maosheng, WANG Bingnan.Method of Monitoring Urban Area Deformation Based on Differential TomoSAR[J]. Acta Geodaetica et Cartographica Sinica,2016,45(12):1413-1422. DOI:10.11947/j.AGCS.2016.20160113.

P237

A

1001-1595(2016)12-1413-10

國家發(fā)改委衛(wèi)星及應(yīng)用產(chǎn)業(yè)發(fā)展專項([2012]2083)

2016-03-25

王愛春(1981—)，男，博士生，工程師，研究方向為多基線干涉SAR處理方法及應(yīng)用。

E-mail：wangaichun@cresda.com

向茂生

E-mail：xms@mail.ie.ac.cn

修回日期：2016-10-19