基于模型迭代匹配的單幅高分辨SAR圖像遮擋建筑物高度估計方法

蔣李兵，王壯，雷琳，郁文賢

(1. 國防科學技術大學 ATR重點實驗室，湖南 長沙 410073；2. 上海交通大學 電子信息與電氣工程學院，上海 200240)

1 引言

SAR是一種以主動模式工作的微波遙感設備，具有全天時、全天候工作的特點，在遙感對地觀測體系中扮演著重要的角色。近年來，災害監測和評估成為高分辨SAR遙感應用的新熱點[1]。例如，在2008的四川汶川地震和2009年的意大利拉奎拉地震中，TerraSAR-X和Cosmo-SkyMed衛星在震后及時獲取了受災地區的高分辨 SAR圖像，為災害評估提供數據支持。建筑物作為一類重要的人造目標，其高度是災后評估中的一項重要指標。

雷達立體測量[1]（radargrammetry）和干涉測量(interferometry)[2]是2種經典的SAR建筑物高度估計方法。然而，這2種方法對SAR圖像數據均有一定的前提要求：立體測量法需要事先對已配好的2幅或多幅SAR圖像，干涉測量法則常需借助于額外的輔助數據（如地面高程數據、光學圖像數據等）來解決高度估計中的關鍵問題——相位解纏，這令它們的應用范圍受到一定的限制。在一些時效性要求很強的場合（如地震災害應急救援），受衛星重返周期、航空飛行氣象條件等因素的影響，在某個時間段內往往只能可靠地獲得單幅 SAR圖像。如何從單幅 SAR圖像中估計建筑物高度信息已成為一個現實而亟待解決的問題。特別地，在城區等建筑物密集分布區域，相鄰建筑物在沿雷達視線方向（RLOS, radar line of sight）上易形成遮擋而導致其圖像特征發生變化，給高度估計問題帶來新的挑戰。

圍繞著單幅SAR圖像的建筑物高度估計問題，相關領域學者開展了一系列的研究工作[3~12]。目前已有的成果主要集中在孤立情況下的建筑物高度估計，按照處理思路的不同，可以歸納為直接測量法[3~5]、電磁反演法[6,7]和模型匹配法[8~12]3大類。直接測量法人工或半自動地提取疊掩和陰影特征，通過測量它們的尺寸來推算建筑物高度，其估計誤差受特征提取精度的制約。電磁反演法則依據電磁散射理論，構建 SAR圖像中的可觀測量與建筑物電磁參數、幾何參數之間的確定性方程關系來求解高度。此類方法不僅要求構建準確的后向散射模型，而且還需要目標的粗糙度、介電常數等多個電磁參數作為先驗，前提條件苛刻，實用性不高。模型匹配法從一系列模型假設中搜索與觀測圖像最為匹配的解，是目前高度估計研究的一大熱點。其中，用于模型匹配可以發生在圖像域[8,9]，也可以選擇在強散射點、邊緣等低層特征域[10,11]進行。在圖像域上進行的模型匹配需要通過電磁仿真來模擬 SAR圖像，算法的計算復雜度偏高。而在低層特征域上進行模型匹配時，除了要根據模型假設計算建筑物的理論低層特征外，還需要額外地從 SAR觀測圖像中提取相應的實際低層特征，易受相干斑和特征提取算子的影響而導致匹配錯誤。而針對遮擋情況下的建筑物高度估計研究目前尚處于探索階段，由于遮擋情況的建筑物 SAR圖像特征與孤立建筑物有顯著不同，直接沿用現有的孤立建筑物高度估計方法并不合理。

本文從遮擋情況下的建筑物圖像特征出發，提出了一種針對遮擋情況的平頂建筑物高度估計方法。該方法首先通過本文構建的遮擋情況建筑物特征信號計算模型，將建筑物幾何模型假設變換為特征信號。接著設計了度量特征信號與SAR觀測圖像一致程度的似然函數，在模型假設空間中搜索最優解。與經典的單幅SAR圖像建筑物高度估計方法相比，本文方法既不需要從觀測SAR圖像中提取低層特征，也無須進行復雜的SAR圖像模擬計算。

2 受遮擋建筑物圖像特征分析

建筑物在 SAR圖像中呈現為疊掩、二面角、陰影等多類區域的集合，它們是 SAR側視距離成像的結果[13]，其形狀、尺寸和位置與建筑物幾何參數直接相關。圖1給出了2種孤立平頂建筑物的特征區域分布情況。來自墻面、地面及屋頂的等距單次反射疊加于同一分辨單元內，形成較亮的疊掩反射區（圖1中a區）；同時，正對RLOS方向的建筑物墻面和地面構成二面角反射器，在圖像上形成一條亮細線（圖1中b區）；屋頂的單次反射構成了次亮的屋頂反射區（圖1中c區）；而在背向RLOS方向的一側，則因無反射回波而形成陰影（圖1中d區）。本文將上述區域的集合定義為建筑物特征信號，它不包含建筑物的散射強度信息，但卻刻畫了建筑物在斜距平面上全體特征區域的形狀、尺寸和位置，因而可以理解為是表征建筑物結構的二維掩模。孤立建筑物特征信號與高度的關系如式(1)所示。

圖1 2類孤立平頂建筑物的特征信號

然而，當建筑物分布密集時，相鄰建筑物會在沿RLOS方向上產生遮擋，其特征信號也隨之發生改變。圖2以2個相鄰的建筑物B0和B1為例，給出了不同高度配置時建筑物特征信號剖面圖。B0的高度固定為h0，當B1的高度為h1時，兩建筑物間不存在遮擋，它們的斜距剖面均完整地包含著疊掩、二面角、屋頂反射和陰影。隨著B1的高度增加至 h2，兩建筑物之間一部分對 B0的疊掩有貢獻的地面被 B1所遮擋，導致兩者的特征信號均有所變化：B1的陰影被縮短，B0的疊掩則被劃分為強度不一的2部分，且其二面角反射強度也有所減弱。這種效應隨著 B1高度的增加而加劇。當 B1的高度增加至 h3時，B1對 B0的遮擋不僅使B1自身的陰影被進一步縮短，還導致B0的二面角徹底消失。

圖2中的例子說明，受遮擋建筑物的特征信號不僅與自身的幾何參數有關，而且還與引起遮擋建筑物和被其遮擋的建筑物有關。因此，要估計受遮擋建筑物高度，必須將受遮擋影響的所有相鄰建筑物目標同時進行考慮。

圖2 相鄰平頂建筑物在不同高度配置時的特征信號

3 遮擋情況的建筑物高度估計問題建模

以往在研究建筑物高度估計問題時，一般均假定建筑物除高度之外的其他幾何參數（如長、寬、方位角、空間位置）是已知的[6,7,9,10]。這些幾何參數可以通過災前的光學圖像或地理信息系統直接獲知，也可以經由特征提取步驟從當前的單幅SAR觀測圖像中估計，文獻[10]和文獻[14]對此有詳細的介紹。本文依然沿用這一假定，并以常見的一類建筑物——平頂建筑物為例，研究遮擋情況的建筑物高度估計問題。

將平頂建筑物用參數模型表示，記為B(l, w,h, α, C )，其中，l,w,h,α分別代表建筑物的長、寬、高和方位角， C = ( d x, d y)是建筑物底面中心坐標。SAR傳感器的高度和入射角分別是H和β。為方便表述，將建筑物和SAR的已知參數用矢量形式表示，分別記為 θ =(l, w, α, C)和 φ =(H , β)。依據模型匹配方法[9]的一般思路，基于單幅SAR圖像的建筑物高度估計問題可以表示為如下形式

式(2)中的 S ( B(h, θ), φ)（在下文中簡記為 S (·)）是在成像條件φ下，從建筑物幾何參數模型B(h, θ)到某一特征域的計算模型。I是已知的SAR觀測圖像， L{·}是度量 S(·)與觀測圖像I一致程度的似然函數。

考慮到受遮擋建筑物高度估計問題的特點，需要同時估計其周圍受到遮擋影響的所有建筑物高度，式(2)相應地改寫為

其中， h = (h1, h2,… ,hn)是受遮擋影響的相鄰建筑物目標集{B(hi, θi) ,i = 1,… ,n }的高度矢量。

求解式(3)需要確定用于模型匹配的特征域，并設計相應的計算模型 S (·)和似然函數L{·}。如引言中所述，圖像域和低層特征域是用于模型匹配的 2種常用特征域，但當面對遮擋情況時，它們均存在著局限性：圖像域上的模型匹配過程涉及電磁仿真計算，其復雜度將隨著建筑物數目的增加而急劇上升；而在低層特征域進行模型匹配時，相干斑現象會直接影響低層特征提取的精度。考慮到遮擋情況建筑物高度估計問題的實際特點，本文并沒有選擇在上述 2種特征域上進行模型匹配，而是將匹配過程作用了建筑物特征信號域，并定義了與之相適應的似然函數。

4 遮擋情況的建筑物高度估計方法

本文提出了一種從單幅 SAR圖像中估計受遮擋平頂建筑物高度的方法。算法從 SAR成像參數和建筑物幾何參數模型出發，首先根據建筑物高度假設計算對應的特征信號。接著通過本文所設計的似然函數，與 SAR觀測圖像進行匹配。該“高度假設——特征信號計算——匹配驗證”過程中的最優高度假設求解通過遺傳算法實現。圖3歸納了算法的主要流程。

圖3 遮擋情況的建筑物高度估計方法流程

4.1 遮擋情況的建筑物特征信號計算模型

雷達散射機理表明，建筑物的疊掩、陰影等特征區域的形狀和尺寸主要取決于建筑物及其所處環境與雷達之間的照射關系，而它們在斜距平面上的位置則由起主導作用的后向反射決定。對于本文所研究的平頂建筑物（如圖 2所示），其疊掩是由受雷達回波直接照射的建筑物表面及地面的單次反射相疊加而成，陰影取決于被遮擋的建筑物表面及地面，而二面角則固定地位于墻—地交線處。結合圖1和圖2中的例子和前文的特征信號定義可以得出如下結論：平頂建筑物在斜距平面上的特征信號可以由單次反射完全確定。于是，根據斜距成像原理，建筑物表面上直接被雷達照射的點(,,)x y z經單次反射后在 SAR斜距平面上的位置(r,a)可以由正交投影來近似，投影方程為

于是，遮擋情況的建筑物特征信號計算模型可以分為 2步：1）計算直接被雷達照射的建筑物表面及周圍地面；2）將 1）中的點投影至斜距平面。

上式中的(x, y, z)滿足方程組(5)的條件為

建筑物B周圍未被遮擋的地面unG 也可根據式(5)~式(7)類似地計算得到，此處不再贅述。

對于給定的建筑物高度假設h，按式(5)~式(7)逐一計算其中每一個建筑物及周圍地面未被遮擋的部分并相加，代入式(4)中，就得到了相鄰建筑物的特征信號(,)r a∪。它雖然不包含反射強度，但卻描述了疊掩、陰影和單反射區域的形狀、尺寸及位置信息。

4.2 似然函數

4.1節給出了遮擋情況下建筑物特征信號的計算模型。盡管該特征信號并不包括散射強度信息，但卻完整地刻畫了建筑物在斜距平面上的疊掩、陰影等不同特征區域的形狀、尺寸及位置。將該特征信號疊加于 SAR觀測圖像上，就可以比較兩者之間的相似程度，從而間接地實現高度參數假設與SAR觀測圖像之間的匹配。

記 PI(h)為特征信號 ∪ ( r, a)疊加于SAR觀測圖像I之上所得到的圖像區域，它可以用全體特征區域的集合{Ri} 表示，也可以等價地用特征區域邊緣段的集合{Ej}表示。

其中，N、M分別是特征區域和邊緣段的個數。

建筑物的特征區域是由具有不同局部入射角和反射類型的后向反射疊加形成的，一般來說，屬于同一特征區域的像素強度趨于一致，而不同特征區域間的像素強度則差異明顯。從這個性質出發，分別構建了區域均勻度函數region()C ·和邊緣段對比度函數edge()C ·。

其中，I( k)是圖像I在k處的像素值，μi是區域Ri的對數均值，grad(·)是梯度算子。

從而， PI(h)的似然函數定義為區域均勻度函數和邊緣段對比度函數的加權和。

其中，w是預定義的加權因子。

4.3 基于遺傳算法的高度參數求解

式(3)是一個目標函數極大值求解問題，通常可由確定性方法或隨機優化方法進行求解。由于似然函數 ()Lh的空間曲面性質未知，為了避免陷入局部極值，本文采用隨機優化算法中的遺傳算法求解式(3)。遺傳算法是一種自適應全局優化概率搜索算法，其求解主要包括染色體編碼、種群初始化、適應度函數、遺傳算子和終止條件5個方面的內容。針對本文的似然函數極值求解問題，遺傳算法的具體參數設置如下。

1) 染色體編碼

由于待求解的高度參數h是由n個實值變量構成，因此采用實數編碼方式，即染色體由長度為n的實數向量 ( q1, q2,… ,qn)表示。

2) 種群初始化

初始種群采用指定和隨機相結合的方式生成。設初始種群規模為 p，其中，一定比例的染色體由沿建筑物不同剖面上的疊掩和陰影尺寸從式(1)計算得到，剩余的染色體則在參數取值范圍內隨機生成。盡管大部分直接指定的染色體受測量誤差的影響而與真值相差甚遠，但從概率上來說，也會存在位于真值附近的染色個體。采用這種方式對種群進行初始化，可以在保證種群多樣性的同時提高種群收斂到真值附近的概率。

3) 適應度函數

遺傳算法中的適應度標識了個體生存的概率，適應度較高的個體遺傳到下一代的可能性較大。這與本文所構建的用于度量高度假設h的似然函數L(h)是一致的，因此直接以式(11)定義的似然函數作為適應度函數。

4) 遺傳算子

遺傳算子主要包括選擇、交叉和變異3部分。本文采用比例選擇算子的方法來進行選擇運算，采用算術交叉的方法來進行交叉運算，采用均勻變異的方法進行變異運算。這3種算子均是遺傳算法中的常用算子，此處不再贅述。

5) 終止條件

遺傳算法的終止條件設置為迭代次數達到預設的最大迭代次數g或連續c代平均適應度不發生變化。

5 實驗結果與分析

為驗證本文方法的有效性，首先對模擬圖像進行了實驗，分析了高度估計的誤差，并與直接測量法的估計結果進行了比較。隨后，在模擬圖像實驗的基礎上，將本文方法用于實測圖像。其中，模擬圖像是基于Lambertian反射模型，根據不同的建筑物幾何參數模型、噪聲強度及成像條件仿真生成的。實測圖像來源于德國 DLR公開發布的 Noerdlingen地區TerraSAR-X衛星遙感圖像。

圖4以包含2個相鄰平頂建筑物的場景為例，給出了高度估計的中間過程和最終結果。圖4(a)是仿真場景的三維示意，兩相鄰建筑物的參數模型分別為 B1(h1,θ1)= (60,(80,30,0°,(25,100)))和 B2(h2,θ2)=(40,(80,30,0°,(95,80)))。圖 4(b)是根據圖 4(a)的場景按Lambertian反射模型計算得到的模擬SAR圖像，成像參數θs(H,β)=(5 000,45°)，圖像尺寸90×350，距離和方位分辨率均為 0.3m，等效視數（ENL,equivalent number of looks）為3。遺傳算法種群數設為 20，依直接測量法得到的 2組高度估計值(65.0,43.1)和(24.0,43.1)作為初始種群中的高度假設初值。圖 4(c)中的灰色細輪廓代表高度(24.0,43.1)所對應的建筑物特征信號。隨著遺傳進化的進行，每一代最優高度假設個體所對應的特征信號不斷逼近 SAR觀測圖像中實際建筑物各個特征區域的真實邊界。圖4(d)中的灰色細輪廓是對應著高度估計值=(59.7,39.1)的建筑物特征信號，其與建筑物的實際特征區域輪廓相吻合。從數值上看，估計結果(59.7,39.1)與高度真值(60,40)非常接近，相較于直接測量法的2個估計結果(65.0,43.1)和(24.0,43.1)，精度有了較大提高，表明本文方法用于估計遮擋情況的建筑物高度是有效的。圖4(e)是最優個體適應度隨遺傳代數的演化曲線。由于具有可靠的高度初值，使得遺傳進化過程能較快地以大概率收斂于真值附近。

進一步選取不同方位角、成像入射角和圖像噪聲方差的SAR仿真圖像共40幅來測試本文方法的性能。仿真場景中包含2個相互平行的建筑物，建筑物方位角范圍選取為 0°、10°、20°、30°、40°、50°、60°、70°、80°和 90°共 10 組，每組包含 40°和50°這2種成像入射角，噪聲方差則取為0.1和0.3（對應的ENL分別為3和1）。圖4(f)給出了本文算法得到的高度估計絕對值誤差統計結果，其中，h1和h2分別為近雷達側和遠雷達側建筑物的高度。經統計得到的高度估計絕對誤差均值約為0.5m，最大估計誤差在 1m左右。圖4(f)的統計結果曲線說明，本文方法能有效估計受遮擋程度不一的建筑物高度。同時，在不同方位角下的 h1和 h2誤差曲線起伏不明顯，表明本文的高度估計方法對建筑物方位角不敏感。不同噪聲方差下的誤差曲線也較為接近，說明估計誤差并未隨噪聲水平的增大而顯著變大，算法具有一定的頑健性。

針對上述同樣的數據集，還選用了文獻[10]中經典的單幅 SAR圖像孤立建筑物高度估計方法與本文方法進行比較。其中，文獻[10]算法涉及的圖像先驗統計分布選用Gamma分布，且假定特征信號的最優匹配位置已知。仿真實驗得到的絕對估計誤差統計曲線如圖4(g)所示。高度估計h1的誤差曲線與本文算法得到的結果相近，但由于受建筑物部分區域被遮擋的影響，絕對誤差均值上升為1.2m。而h2的誤差曲線則遠高于本文算法所得到的結果，絕對誤差均值達到 8.6m。同時，從 h2誤差曲線的走勢可以發現，文獻[10]方法的估計誤差隨著受遮擋程度的加劇而趨于增大。上述比較實驗進一步印證了本文算法針對遮擋建筑物目標高度估計的有效性。

圖5給出了基于實測SAR圖像的建筑物高度估計結果。圖5(a)是德國Noerdlingen地區的TerraSAR-X衛星在高分辨聚束（HighRes SpotLight 300 MHz）模式下所獲得的SAR圖像切片，分辨率為1m，局部入射角為32°；圖5(b)是對應的光學圖像（@Google Earth），其中包含了3棟相鄰的平頂建筑物（分別記為 A、B、C）。圖 5(c)是由本文算法估計得到的高度參數所對應的特征信號（灰色細輪廓）與原圖疊加的結果。可以看出，該特征信號與 3棟相鄰建筑物的整體特征區域輪廓基本吻合，較準確地反映了建筑物目標A和C之間由于遮擋所引起的特征信號變化。圖5(d)是依據高度估計結果反演的三維場景示意。上述基于實測 SAR圖像的實驗結果進一步說明了本文算法的有效性。

圖4 兩相鄰建筑物遮擋情況的高度估計結果

圖5 實測SAR圖像建筑物高度估計結果

6 結束語

對城區中建筑物進行高度估計在雷達遙感領域有著重要的應用。遮擋效應導致建筑物圖像發生變化，使得現有的方法難以對遮擋建筑物特征進行準確的描述，降低了高度估計的精度。本文從遮擋效應產生的幾何機理出發，通過分析受遮擋建筑物的圖像特征，構建了其特征信號的計算模型，進而提出了一種基于模型迭代匹配的單幅 SAR圖像受遮擋建筑物高度估計方法。基于模擬和實測 SAR圖像的實驗結果表明，本文方法用于估計遮擋情況下的建筑物高度是可行的。

需要指出的是，本文提出的方法目前僅限于平頂建筑物。對于外形不規則的一類建筑物，其與雷達波的相互作用類型較為復雜，僅用單次和二次反射無法完整描述。下一步的研究將致力于解決該問題。

[1] SOERGEL U. Radar Remote Sensing of Urban Areas[M]. Berlin:Springer, 2010. 1-8.

[2] 馬侖,廖桂生,李真芳. 利用分布式小衛星InSAR系統獲取寬域、高分辨、高精度三維地形[J]. 電子學報, 2009, 37(9):1900-1906.MA L, LIAO G S, LI Z F. Using constellation InSAR system to reconstruct three-dimensional topography with wide swath, high resolution and high precision[J]. Acta Electronica Sinica, 2009,37(9): 1900-1906.

[3] TUPIN F. Extraction of 3D information using overlay detection on SAR images[A]. 2nd GRSS/ISPRS Joint Workshop on Data Fusion and Remote Sensing over Urban Areas[C]. Berlin, Germany, 2003.72-76.

[4] BENNETT A J, BLACKNELL D. The extraction of building dimensions from high resolution SAR imagery[A]. IEEE Radar Conference[C].Huntsville, Alabama, USA, 2003. 182-187.

[5] THIELE A, CADARIO E, SCHULZ K, et al. Analysis of gable-roofed building signature in multiaspect InSAR data[J]. IEEE Geoscience and Remote Sensing Letters, 2010, 7(1):83-87.

[6] FRANCESCHETTI G, GUIDA R, IODICE A, et al. Building feature extraction via a deterministic approach: application to real high resolution SAR images[A]. Proceedings of IEEE International Geoscience and Remote Sensing Symposium (IGARSS2007)[C]. Barcelona, Spain,2007. 2681-2684

[7] GUIDA R, IODICE A, RICCIO D. Height retrieval of isolated buildings from single high-resolution SAR images[J]. IEEE Trans on Geoscience and Remote Sensing, 2010, 48(7):2967-2979.

[8] QUARTULLI M, DATCU M. Stochastic geometrical modeling for built-up area understanding from a single SAR intensity image with meter resolution[J]. IEEE Trans on Geoscience and Remote Sensing,2004, 42(9):1996-2003.

[9] BRUNNER D, LEMOINE G, BRUZZONE L, et al. Building height retrieval from VHR SAR imagery based on an iterative simulation and matching technique[J]. IEEE Trans on Geoscience and Remote Sensing, 2010, 48(3):1487-1504.

[10] SPORTOUCHE H, TUPIN F, DENISE L. Building detection and height retrieval in urban areas in the framework of high resolution optical and SAR data fusion[A]. 2010 IEEE International Geoscience& Remote Sensing Symposium (IGARSS2010)[C]. Honolulu, Hawaii,USA, 2010.3660-3663.

[11] 趙凌君. 高分辨率 SAR圖像建筑物提取方法研究[D]. 長沙: 國防科技大學, 2009. 107-128.ZHAO L J. Building Extraction from High Resolution SAR Imagery[D]. Changsha: National University of Defense Technique, 2009.107-128.

[12] BARTHELET E, MERCIER G, DENISE L. Building change detection in a couple of optical and SAR high resolution images[A]. 2011 IEEE International Geoscience & Remote Sensing Symposium (IGARSS2011)[C].Vancouver, Canada, 2011. 2393-2396.

[13] XIA Z G, HENDERSON F M. Understanding the relationship between radar response patterns and the bio-and geophysical parameters of urban areas[J]. IEEE Transactions on Geoscience and Remote Sensing,1997, 35(1):93-101.

[14] FERRO A, BRUNNER D, BRUZZONE L. Detection and reconstruction of building footprints from single VHR SAR images[A]. European Conference on Synthetic Aperture Radar (EUSAR)[C]. Aachen,Germany, 2010. 1094-1097.