基于高分三號SAR數(shù)據(jù)的城市建筑高分辨率高維成像

畢 輝 金 雙 王 瀟 李 勇 韓 冰 洪 文

①(南京航空航天大學電子信息工程學院 南京 211106)

②(南京工業(yè)大學計算機科學與技術學院 南京 211816)

③(中國科學院空天信息創(chuàng)新研究院 北京 100094)

1 引言

與傳統(tǒng)光學觀測手段不同，合成孔徑雷達(Synthetic Aperture Radar,SAR)具有全天時、全天候的工作能力，因而在國土資源勘測、自然災害監(jiān)測等領域得到廣泛應用[1]。然而傳統(tǒng)SAR成像只能獲取目標的方位-距離二維圖像，無法準確反映目標的三維散射特征，一定程度上影響了SAR圖像的進一步應用。層析合成孔徑雷達(Synthetic Aperture Radar Tomography,TomoSAR)是SAR成像技術進一步擴展。它將合成孔徑原理延伸至高程向，可基于多景二維復圖像獲得目標的方位-距離-高程信息，進而實現(xiàn)三維成像[2]。1999年，Reigber等人[3,4]首次展示了機載TomoSAR成像結果，實現(xiàn)了基于譜估計方法的森林區(qū)域三維成像。2006年，F(xiàn)ornaro和Serafino[5]基于長時間基線星載SAR數(shù)據(jù)，開展了TomoSAR成像星載實驗，證實了譜估計技術應用于高程向重構的可行性與有效性。2011年，Reale等人[6]證實了高分辨率數(shù)據(jù)和先進干涉處理技術相結合可以更好地重建建筑物三維結構。2013年，Shahzad和Zhu[7]提出了一種全新的建筑物立面重建方法，獲取了拉斯維加斯百樂宮酒店的三維點云，證實了TomoSAR點云在構建動態(tài)城市模型方面的優(yōu)越性。2017年，Wang等人[8]提出了迭代重加權的交替方向乘子算法，用于實現(xiàn)快速TomoSAR成像。2018年，Wang和Zhu[9]提出了一種基于核主成分分析的TomoSAR成像方法，使用極小代價即可分離同一方位-距離分辨單元中沿高程向分布的多個散射體。2019年，秦斐等人[10]針對TomoSAR成像中高程向分辨率較低、建筑物疊掩、提取建筑物目標特征效率較低等問題，提出了一種基于機器學習的建筑物目標識別和提取算法，提高了觀測目標的特征提取效率，并通過機載陣列三維SAR實驗數(shù)據(jù)驗證了方法的有效性。差分層析合成孔徑雷達(Differential Synthetic Aperture Radar Tomography,D-TomoSAR)是TomoSAR的進一步擴展。它基于多基線觀測數(shù)據(jù)可實現(xiàn)對目標的方位-距離-高程-時間四維成像，不僅解決了TomoSAR成像中的高度錯位與模糊問題，還可高精度獲取目標的形變信息。D-TomoSAR概念由Lombardini于2003年首次提出[11]。2007年，F(xiàn)ornaro等人[12]證明了D-TomoSAR可作為傳統(tǒng)永久散射體監(jiān)測技術的有效替代方案，實現(xiàn)對大場景形變的有效監(jiān)測。2008年，F(xiàn)ornaro等人[13]提出了一種可分離干擾散射體相關時間序列的技術，并應用于羅馬地區(qū)的D-TomoSAR成像中，證明了該技術可用于城市復雜場景的形變監(jiān)測。2009年，Zhu等人[14]基于多景TerraSAR-X數(shù)據(jù)，實現(xiàn)了對拉斯維加斯的D-TomoSAR成像，獲取了城市四維雷達地圖。2010年，F(xiàn)ornaro等人[15]使用ERS數(shù)據(jù)對羅馬市Grotta Perfetta地區(qū)進行了D-TomoSAR成像，結果顯示差分層析成像方法可以有效區(qū)分同一分辨單元中的不同散射體，克服了差分干涉的局限性，進一步提高了對城市基礎設施形變的監(jiān)測能力。2011年，Zhu和Bamler[16]提出了一種“時間扭曲”方法，并基于TerraSAR-X數(shù)據(jù)實現(xiàn)了城市區(qū)域的D-TomoSAR成像，有效獲取了城市建筑的線性運動和季節(jié)性運動速率。2015年，Siddique等人[17]將D-TomoSAR方法和永久散射體干涉技術相結合，基于50幅TerraSAR-X圖像獲得了巴塞羅那中高層建筑物立面散射體的時空反演結果。2020年，Wang和Liu[18]提出了一種廣義D-TomoSAR成像系統(tǒng)模型和一種基于擬極大似然的成像算法，同時反演出了建筑物線性運動、季節(jié)性運動等多個形變運動速率。目前，TomoSAR和D-TomoSAR成像技術已在城市建筑物和基礎設施的三維重建和長期形變監(jiān)測方面展現(xiàn)了極大的應用潛力。

壓縮感知(Compressive Sensing,CS)是一種重要的稀疏信號處理技術，它可以使用較少樣本實現(xiàn)對稀疏信號的高質(zhì)量恢復[19–21]。CS-TomoSAR成像的前提是觀測場景的高程向分布稀疏，而城市區(qū)域主要為人造建筑，其高程向分布都滿足稀疏性條件。因此，CS在城市區(qū)域三維、四維成像中具有廣闊的應用前景。起初，由于現(xiàn)代米級分辨率星載SAR系統(tǒng)軌道限制，TomoSAR成像高程向分辨率遠低于方位向和距離向，因而迫切需要超分辨算法來解決這一問題，2010年，Zhu和Bamler[22]介紹了一種基于CS理論的TomoSAR成像方法，相比于傳統(tǒng)譜估計方法，該算法實現(xiàn)了對目標高程向分布的超分辨重構。2012年，Zhu和Bamler[23]將所提出的SL1MMER算法應用到TerraSAR-X星載數(shù)據(jù)處理中，獲得了拉斯維加斯百樂宮酒店的高分辨率TomoSAR成像結果，并證明了該算法具有超分辨能力。2015年，Weiss等人[24]提出了一種適用于TomoSAR的自適應CS算法，準確識別了同一分辨單元中的兩個散射體位置。2017年，Li等人[25]研究了基于SPICE的TomoSAR成像方法，并利用8幅TerraSAR-X條帶影像，實現(xiàn)了對內(nèi)蒙古根河市某建筑的高精度三維重建。2010年，Zhu和Bamler[26]將CS技術應用在D-TomoSAR成像中，證明了在高程向多散射體分離上CS技術相比于傳統(tǒng)譜估計方法的優(yōu)越性，同時指出CS可自動識別散射體數(shù)量，非常適用于星載SAR系統(tǒng)的三維、四維成像。2010年，Zhu和Bamler[27]基于SL1MMER算法重構獲得了拉斯維加斯會議中心建筑群的D-TomoSAR成像結果，展示了CS技術在四維成像方面的優(yōu)勢和能力。2014年，Leng等人[28]將最小絕對收縮和選擇算子CS算法應用于建筑區(qū)域，并展示了巴塞羅那的高程重建和形變監(jiān)測結果。

高分三號衛(wèi)星是我國首顆分辨率達到1 m的C頻段多極化SAR衛(wèi)星，于2016年8月10日在中國太原衛(wèi)星發(fā)射中心由長征四號丙運載火箭發(fā)射升空[29]。高分三號衛(wèi)星是“國家高分辨率對地觀測系統(tǒng)重大專項”中唯一的民用微波遙感成像衛(wèi)星，具有高分辨率、大成像幅寬、多成像模式、長壽命運行等特點，可實現(xiàn)全天時、全天候的全球海洋與陸地監(jiān)測[30]。目前高分三號已成功應用于高精度測繪、自然災害監(jiān)測等多個領域[31,32]。然而，由于設計之初未考慮后續(xù)高維成像應用，現(xiàn)有高分三號獲取的SAR圖像存在有一定的空間、時間去相干問題，對其進一步應用于干涉SAR、差分干涉SAR、TomoSAR、D-TomoSAR等存在一定挑戰(zhàn)。2019年，余博等人[33]基于高分三號數(shù)據(jù)，對河南省登封市周圍地區(qū)進行了干涉測量實驗，通過與哨兵一號衛(wèi)星數(shù)據(jù)結果進行對比分析，驗證了高分三號的干涉能力以及可提取地表形變信息的能力。2021年7月，黃震等人[34]基于高分三號SAR數(shù)據(jù)進行了干涉測量實驗，成功提取了觀測區(qū)域的數(shù)字高程模型。

本文基于7景高分三號SAR復圖像數(shù)據(jù)，利用CS技術，開展了TomoSAR和D-TomoSAR成像實驗研究，獲取了北京市雁棲湖周圍建筑的高分辨率三維、四維SAR圖像，實現(xiàn)了建筑物的高質(zhì)量三維重建以及高精度形變監(jiān)測，為后續(xù)基于高分三號SAR數(shù)據(jù)的干涉系列應用及多維高分辨率成像提供了技術支撐。

本文后續(xù)結構如下：第2節(jié)主要介紹TomoSAR,D-TomoSAR成像模型，并給出了上述兩個模型的CS求解方案；第3節(jié)介紹了本文所使用的高分三號SAR復圖像數(shù)據(jù)集；第4節(jié)基于仿真數(shù)據(jù)，開展了TomoSAR,D-TomoSAR成像實驗，證明了CS技術在高分辨三維成像及高精度形變監(jiān)測方面的有效性；第5節(jié)基于7景高分三號SAR復圖像數(shù)據(jù)，對北京雁棲湖周圍建筑進行了TomoSAR,D-TomoSAR成像研究，獲取了代表性建筑和大觀測區(qū)域的三維、四維雷達圖像；第6節(jié)對文章進行了總結和展望。

2 成像模型

2.1 TomoSAR成像模型

TomoSAR利用對同一場景觀測獲取的多幅配準的二維SAR復圖像(多基線觀測數(shù)據(jù))在高程向上進行孔徑合成，以獲得高程向上的分辨能力，進而重構目標的三維散射信息[3,4,35,36]。TomoSAR成像幾何如圖1所示。設共有N條基線用于數(shù)據(jù)獲取，令bn(n=1,2,...,N)表示高程向孔徑分布，對于一個選定的方位-距離分辨單元，第n幅SAR圖像對應的聚焦測量值可以表示為

圖1 TomoSAR成像幾何Fig.1 TomoSAR imaging geometry

式中，ξn=-2bn/(λr)表 示高程向頻率，其中λ為波長，r為斜距；γ(s)表 示沿高程向s的復反射函數(shù)；Δs為 高程向跨度。沿高程向s對高程向復反射函數(shù)γ(s)進行離散化，則式(1)中的成像模型可近似表示為

式中，L為高程向離散化點數(shù)；g=[g1,g2,...,gN]T表示測量值向量；R=exp(-j2πξnsl)為根據(jù)Tomo-SAR成像幾何所構建的觀測矩陣；γ=[γ(s1),γ(s2),...,γ(sl)]T表示高程向離散復反射函數(shù)，其中sl(l=1,2,...,L)為離散高程向分布。從式(2)可以看出，TomoSAR成像模型可視為對γ(s)不規(guī)則采樣的離散傅里葉變換。因此，一個SAR測量值可以看作目標復反射函數(shù)沿高程向的一個譜參數(shù)。對于非參數(shù)化譜分析問題，高程向理論分辨率ρs依賴于高程孔徑大小 Δb，在高程向采樣密集的情況下，ρs可由式(3)進行計算。

2.2 D-TomoSAR成像模型

相比于TomoSAR,D-TomoSAR在三維的基礎上多了一個時間方向的維度。它沿高程向和形變速度向合成兩個孔徑，進而獲取被觀測目標的高程和形變速率的聯(lián)合分辨率，實現(xiàn)四維成像[11]。對于N個復圖像而言，當高程孔徑位置為bn、時間基線為tn時，第n幅圖像的聚焦測量值可以表示為

近年來，隨著我國社會經(jīng)濟的迅速發(fā)展，我國建筑工程項目數(shù)量日益增加，而工程測量是保證建筑工程項目正常進行的前提。傳統(tǒng)的工程測量必須要花費很多時間、人力和物力，并且不能確保測量精準度。地面三維激光掃描技術具有高效率、高精準度等優(yōu)點，將其在工程測量中應用，能夠明顯提升工程測量精準度。

式中，ηn=-2πtn/λ表 示形變速度頻率；V(s)表示形變速率。式(4)中的模型也可寫為

式中，Δv表 示被觀測目標的形變速率跨度；δ(·)是與形變項相關的譜分布。令aγ(s,v)=γ(s)δ(v-V(s))，則式(5)可以寫成

式(6)中的模型可以視為a(s,v)在高程-形變平面的二維傅里葉變換。因此，其在立面軸上的投影為反射率剖面γ(s)[37]。將式(6)中的s和v離散化后，D-TomoSAR成像模型可以表示為

其中，g=[g1,g2,...,gN]T表 示測量向量；R=exp(-j2π·(ξnsl+ηnvq))為D-TomoSAR成像觀測矩陣；sl(n=1,2,...,L)為 離散高程向分布；vq(q=1,2,...,Q)為 離散形變向分布；γ由離散化的a(s,v)組成。若時間孔徑大小 Δt，則形變分辨率ρv可由式(8)進行計算

2.3 CS重建

城市中被觀測目標主要為人造建筑，其高程向分布通常都是稀疏的，即每個方位-距離分辨單元中的散射體個數(shù)有限。因此，當測量矩陣R滿足有限等距性質(zhì)條件時，面向式(2)和式(7)中的模型，本文通過解決如下的最優(yōu)化問題分別實現(xiàn)基于CS的TomoSAR和D-TomoSAR成像。

式中，β為正則化參數(shù)，與噪聲水平和樣本數(shù)目有關。CS算法可以在短時間內(nèi)從獲取的樣本數(shù)據(jù)中實現(xiàn)高質(zhì)量信號恢復[38,39]?；谠撍惴ǖ膬?yōu)勢，本文采用CS算法進行TomoSAR和D-TomoSAR成像。

3 高分三號數(shù)據(jù)集簡介

本文所使用的高分三號數(shù)據(jù)集總共包含7景復圖像，具體參數(shù)如表1所示。該數(shù)據(jù)集的7景圖像的空間基線孔徑大小約為1417 m，時間基線跨度是從2018年6月到2019年9月，共464 d，表2給出其時空基線具體參數(shù)，其時空基線分布情況如圖2所示。本文以2019年3月1日獲取的SAR圖像為主影像，其余6景為輔圖像。輔圖像時空基線位置是相對于主影像計算得到的。本文第4節(jié)將基于表1的參數(shù)進行點目標仿真實驗，第5節(jié)將基于該高分三號數(shù)據(jù)集進行實驗。

圖2 高分三號數(shù)據(jù)集時空基線分布圖Fig.2 Spatial-temporal baseline distribution of GF-3 dataset

表1 高分三號數(shù)據(jù)集參數(shù)Tab.1 Parameters of GF-3 dataset

表2 高分三號數(shù)據(jù)集時空基線參數(shù)Tab.2 Spatial-temporal baseline parameters of GF-3 dataset

4 仿真實驗

基于表1中的實驗參數(shù)，本文設定高程向分布有兩個散射體，模擬產(chǎn)生了7景仿真數(shù)據(jù)，添加信噪比為20 dB的噪聲，并分別使用經(jīng)典譜估計方法和CS算法對高程向進行TomoSAR和D-TomoSAR成像[40]。本節(jié)將展示3種經(jīng)典譜估計算法即Beamforming (BF)[41],Adaptive beamforming (Capon)[42]和Multiple signal classification (MUSIC)[43,44]的仿真結果，用于與CS算法重構結果進行比較，以說明CS算法在TomoSAR,D-TomoSAR成像中的優(yōu)勢。圖3為高程向兩個散射點的TomoSAR成像結果，橫坐標為高程向位置分布，縱坐標為散射點的幅度值。本文設置兩個散射體的距離分別為11 m和50 m。由圖3可以看出，當距離為11 m時，由于小于高程向理論分辨率，3種譜估計算法的成像結果均存在分辨率較低和模糊嚴重等問題，無法準確分離兩個散射點，造成重建失??；而CS算法仍可有效識別兩個散射體，實現(xiàn)高程向的超分辨成像。當距離為50 m時，可以看出，BF算法重構結果具有嚴重的、不規(guī)則的旁瓣；Capon和MUSIC算法相較于BF算法，重構結果具有更低的旁瓣，提升了TomoSAR高程向重構質(zhì)量；而相較于3種譜估計算法，CS則可以更加有效地抑制旁瓣和噪聲，進一步提升了高程向散射體的可分辨能力。圖4為高程向3個散射點的TomoSAR成像結果，3個散射點的間隔是20 m，與兩個散射體的成像結果相似，可以看出CS算法對于多散射體分離也具有很好效果。

圖3 高程向兩個散射點的TomoSAR成像結果(左圖：兩個散射點之間的距離為11 m；右圖：兩個散射點之間的距離為50 m)Fig.3 TomoSAR reconstructed reflectivity profiles of two scattering points along the elevation direction (left image:the distance between two scattering points is 11 m;right image:the distance between two scattering points is 50 m)

圖4 高程向3個散射點的TomoSAR成像結果(3個散射點之間的間隔為20 m)Fig.4 TomoSAR reconstructed reflectivity profiles of three scattering points along the elevation direction (the distance between three scattering points is 20 m)

圖5 D-TomoSAR仿真結果(兩個散射體高程位置為–10 m,10 m；散射體形變速率分別為4毫米/年、–7毫米/年)Fig.5 D-TomoSAR simulation results (elevation position of two scatters are–10 m and 10 m;deformation velocity of two scatters are 4 mm/year and–7 mm/year,respectively)

5 高分三號實驗結果

5.1 代表性建筑物反演

本節(jié)將展示基于此高分三號SAR數(shù)據(jù)集的Tomo-SAR和D-TomoSAR重建結果。首先選取觀測場景內(nèi)兩處代表性建筑進行了TomoSAR和D-TomoSAR成像，分別為某生態(tài)農(nóng)業(yè)公司和雁棲湖會展中心。圖6(a)所示為北京某生態(tài)農(nóng)業(yè)公司的光學圖像，紅色虛線框出的區(qū)域是所關注的目標。圖6(b)為該生態(tài)農(nóng)業(yè)公司的二維SAR圖像。圖7(a)為TomoSAR重建結果，獲取了區(qū)域內(nèi)5棟建筑的高度信息。由圖7(a)可以看出，5個獨立建筑的高度均為35 m，這與實際該5棟完全相同的建筑高度相符，驗證了CS-TomoSAR成像技術的有效性，說明其可以用于高分三號SAR數(shù)據(jù)的三維高精度成像中。圖7(b)為基于CS算法的D-TomoSAR重建結果，可以看出，5棟建筑均有不同程度的形變，這可能與建筑材料隨著季節(jié)變化的熱脹冷縮緊密相關，在日常監(jiān)測中應該有所關注，防止出現(xiàn)相關危險。為更加直觀地展示所關注場景，圖8給出了生態(tài)農(nóng)業(yè)公司的三維點云，更加準確地反映出了建筑物的三維散射結構。

圖6 生態(tài)農(nóng)業(yè)公司Fig.6 Ecological agricultural company

圖7 生態(tài)農(nóng)業(yè)公司高程圖及形變速率圖Fig.7 Elevation and deformation velocity maps of ecological agricultural company

圖8 生態(tài)農(nóng)業(yè)公司三維點云圖Fig.8 3-D point cloud of ecological agricultural company

圖9(a)為北京雁棲湖國際會展中心的光學影像，該建筑真實高度約30 m。圖9(b)為雁棲湖國際會展中心的SAR圖像。圖10(a)給出了CS-TomoSAR成像結果，并在圖11給出了其三維點云圖。由圖10(a)和圖11可以看出，重建結果基本反映了建筑物的三維結構，重構的建筑物高度與實際相吻合，這也說明了本文基于高分三號SAR數(shù)據(jù)，已可實現(xiàn)對復雜建筑物較高精度的三維重建。圖10(b)展示了基于CS技術的雁棲湖國際會展中心D-TomoSAR重建結果。由圖10(b)可以看出，該建筑的形變大約在–10毫米/年～10毫米/年之間，且建筑的左半部分和右半部的線性形變速率正好相反，這反映了該建筑下方地面可能處于一邊抬高一邊塌陷的變化之中，應當著重關注。

圖9 北京雁棲湖國際會展中心Fig.9 Beijing Yanqi lake international convention and exhibition center

圖10 北京雁棲湖國際會展中心高程圖及形變速率圖Fig.10 Elevation and deformation velocity maps of Beijing Yanqi lake international convention and exhibition center

圖11 北京雁棲湖國際會展中心三維點云圖Fig.11 3-D point cloud of Beijing Yanqi lake international convention and exhibition center

5.2 大場景三維、四維成像

本節(jié)展示了基于高分三號數(shù)據(jù)集的大場景TomoSAR和D-TomoSAR成像結果。圖12分別給出了所選取的頂秀美泉小鎮(zhèn)區(qū)域的光學影像和SAR圖像。圖13(a)給出了基于CS算法的TomoSAR成像結果，重構圖像顯示該區(qū)域內(nèi)建筑高度均在15 m至20 m之間，準確反映了小鎮(zhèn)內(nèi)均為4～6層居民樓建筑的實際情況。圖13(b)給出了基于CS算法的D-TomoSAR成像結果。其顯示小鎮(zhèn)左下角區(qū)域和右上角區(qū)域的建筑形變約為–10毫米/年，可能由于這兩塊區(qū)域周圍的某些施工導致地面沉降，小鎮(zhèn)中心的形變速率約為5毫米/年，總體來說比較穩(wěn)定。該結果表明，基于高分三號SAR數(shù)據(jù)，已經(jīng)可以實現(xiàn)大場景的高質(zhì)量三維重建以及高精度形變監(jiān)測，進一步驗證了高分三號SAR衛(wèi)星在城市感知與監(jiān)測中的應用潛力。

圖13 頂秀美泉小鎮(zhèn)區(qū)域的高程圖及形變速率圖Fig.13 Elevation and deformation velocity maps of Dingxiumeiquan town

6 結論

本文開展了基于高分三號SAR復圖像數(shù)據(jù)的TomoSAR和D-TomoSAR成像研究，獲取了北京雁棲湖地區(qū)兩處代表性建筑的三維、四維成像結果，并給出了大面積觀測場景的三維、四維雷達圖像。研究結果顯示了我國高分三號SAR衛(wèi)星應用于TomoSAR和D-TomoSAR成像方面的潛力，為后續(xù)拓展高分三號的干涉系列應用領域提供了技術支撐。

后續(xù)本團隊將繼續(xù)開展相關研究工作，收集高分三號多個區(qū)域、多景數(shù)據(jù)信息，同時探索全新的高分辨TomoSAR,D-TomoSAR成像算法，實現(xiàn)面向復雜城市場景的大范圍三維、四維成像，進一步挖掘高分三號在干涉系列應用中的巨大潛力。