一種三維森林場景極化SAR數據的快速模擬方法

孫晗偉 胡 程 曾 濤

(北京理工大學信息與電子學院 北京 100081)

1 引言

森林是地球生態物質的主體，區域和全球森林生物量的估測對于檢測生態系統碳循環和氣候變化具有重要意義[1]。近年來，合成孔徑雷達(SAR)遙感技術在森林生物信息探測方面的作用得到了廣泛的關注，已成為大區域森林參數估測的有力工具[2?4]。但是參數的估測精度經常受雷達系統不穩定性和自然環境復雜性的影響[5]。為了對各種作用因素進行定量化分析，促進反演算法的發展，需要大量不同雷達系統參數、不同自然環境條件下的遙感數據作支撐。而這些數據往往難以通過實際遙感系統獲取，獲取的數據維度有限，不足以進行全面地分析和研究。因此，針對遙感應用系統的定量化模擬成為獲取大量不同物理條件和客觀環境下遙感數據的重要手段，為反演算法的研究提供了豐富的數據源。在此驅動下，森林微波后向散射模型逐步地發展[6,7]，SAR森林模擬數據得到驗證和應用[8?10]。

目前模擬 SAR森林遙感數據采用兩種圖像直接模擬的方式，一是通過將森林冠層后向散射系數與點擴展函數直接作用得到SAR圖像[8]，二是根據后向散射系數直接得到散射系數圖[9,10]。這兩種方式都避免了復雜的回波模擬和成像處理過程，應用效果明顯，在一定程度上能夠反映SAR遙感系統與森林場景的作用特點。但是，前者在空間中點擴展特點單一，后者不具備點擴展特點，與真實SAR影像數據存在差異。特別是當遙感系統存在非理想因素時(如天線指向出現誤差)，SAR點擴展函數將不再是理想的形式，而且具有空變特性，導致觀測場景內無法用統一的點擴展函數描述[11]。如果仍采用圖像直接模擬方式，將使森林遙感影像數據過于理想，與真實遙感系統的響應結果不符，影響對森林參數反演精度的分析和評估。解決這一問題的有效途徑是在數據模擬過程中引入 SAR原始回波信號模擬的過程，遙感系統各種非理想因素首先作用于回波信號，回波模擬過程可以充分考慮這些因素的影響[11]，使模擬數據能夠反映遙感應用系統的真實特性。

但是，引入回波模擬過程的突出問題是大規模計算量的問題，產生該問題的主要原因是：目前森林雷達后向散射模型已從描述2維概率分布的整體結構模型[6]發展為針對單木建模的 3維結構模型[7?10]，考慮的森林生物學結構越來越豐富；在模擬過程中將3維森林冠層看成由大量的散射介電粒子組成，森林面積越大、密度越大，散射粒子的數量越多，使回波模擬的計算量越大。例如，1 m分辨率1.0×106m2中等密度的森林冠層將離散化為5億個散射粒子，SAR回波數據模擬的時間將超過500 h。

本文針對3維森林場景極化SAR回波數據模擬運算量較大的問題，提出了一種3維森林場景極化SAR數據的快速模擬方法。該方法基于等效散射處理，通過構造一系列虛擬散射粒子，在斜距誤差較小的情況下將森林冠層內多個散射粒子對 SAR回波的綜合貢獻用單個虛擬散射粒子等效，大幅度降低散射粒子的數量，將模擬效率提高1個數量級以上。本文從SAR回波信號模型出發，根據多普勒項誤差限制、線性調頻項誤差限制以及其它限制推導了等效散射處理各控制參數被約束關系，分析了各控制參數的誤差敏感度，得到了模擬精度和效率的優化組合準則。數據域的復相關性驗證表明，經本方法模擬得到的3維森林場景極化SAR回波數據和影像數據在幅度和相位上與傳統回波模擬方法的模擬結果具有較強的一致性，保持了較高的數據精度，同時模擬效率提高了近10倍。

本文首先簡要介紹了3維森林場景極化SAR數據模擬的主要模型和方法；然后給出了等效散射處理原理，提出了基于等效散射的快速模擬方法，詳細給出了3維森林場景切分步驟，推導了控制參數優化組合準則并進行了誤差敏感度分析；最后，通過與傳統回波模擬方法的模擬結果進行數據域復相關驗證了快速模擬方法的精度和效率，同時給出了快速模擬方法用于大范圍森林場景模擬與森林高度反演應用的實例。

2 3維森林場景SAR數據模擬模型與方法

本文討論的3維森林場景極化SAR數據模擬基于SAR回波模擬過程，包括4個步驟，如圖1所示。首先，綜合考慮地形、林木分布、林木大小以及單木生長結構，建立3維森林冠層的幾何模型；其次，根據3維森林冠層內各個散射介電粒子的大小、指向角、雷達的頻率以及電磁波入射和出射方向并采用合適的散射模型計算每個散射介電粒子的多極化后向散射系數；然后，建立SAR遙感系統模型，充分考慮遙感系統的各項特性，包括運動特性、天線特性、接收機特性、極化特性等，模擬生成SAR回波數據；最后，回波數據經過成像處理得到SAR森林遙感模擬影像數據。

圖1 3維森林場景極化SAR數據模擬過程圖

2.1 單木3維結構模型

本文采用整體幾何結構模型構造單木 3維結構[12]。該方法提供從樹根到樹葉的任意路徑，通過樹木分叉的數目、分叉的角度、樹枝長度、半徑變化和各種擾動來構造模型，從整體結構出發用若干參數對樹木造型進行控制。樹木關鍵參數(胸徑、冠層厚度、枝葉大小等)的計算均通過生物學公式或經驗值獲取，在一定程度上保持了所構造樹木的生物特性，同時采用遞歸迭代的計算方法保證了計算機處理的高效性。單木建模效果如圖2所示。

2.2 雷達后向散射模型

3維森林場景可以看成由大量的散射介電粒子組成，這些散射粒子包括圓柱體(樹干和樹枝)、圓盤體(葉片)等[7]。森林的雷達后向散射模型需同時考慮地表的表面散射和樹木冠層的體散射，包括：雷達-地表直接散射，雷達-樹冠直接散射，雷達-樹干直接散射，雷達-樹冠-地表-雷達二次散射以及雷達-樹干-地表-雷達二次散射。采用有限長介電圓柱體(finite dielectric cylinder)散射模型近似計算樹干和樹枝的雷達后向散射[7,10]，采用 GRG(Generalized Rayleigh-Gans)近似的圓盤體散射模型計算葉片的雷達后向散射[7,10]。

圖2 闊葉樹3維建模效果

2.3 SAR回波信號模型

根據后向散射模型，將3維森林冠層看成由大量散射介電粒子組成的復雜介質，其中第k個散射粒子的回波信號模型如式(1)所示。

式中p和q分別表示入射極化方式和接收極化方式，u和t分別代表雷達方位慢時間和距離快時間，Tp為脈沖寬度，l為電磁波波長，c為光速，kr為信號線性調頻率，為后向散射系數，Rk(u)為u時刻的雷達與散射介電粒子的瞬時斜距。

方位向u時刻接收的森林場景回波信號可以表示為該時刻波束照射范圍內K個散射介電粒子的SAR回波信號相干疊加[13]，如式(2)所示。可以看到，回波信號的模擬效率直接決定于散射粒子的數量。

3 等效散射處理原理

由于 SAR回波信號模擬效率與散射粒子的數量有直接關系，因此設法減少散射粒子數量是提高模擬效率的有效途徑。但并不是忽略對某些散射粒子的計算，因為這樣會造成散射信息的缺失導致模擬精度大幅下降。本文提出一種針對3維森林場景的等效散射處理方法，通過構造虛擬的散射粒子，將多個真實散射粒子對 SAR回波信號的綜合貢獻等效為單個的虛擬散射粒子對 SAR回波信號的貢獻。

如圖3所示，設有M個散射介電粒子在空間中位置分布較為接近，其斜距歷程用Rm(u)表示，根據式(1)回波信號模型，M個散射粒子的總回波信號可以表示為

圖3 等效散射處理示意圖

圖 3中E點為虛擬散射粒子，其斜距歷程用近似M個散射粒子的斜距歷程，并用二者SAR孔徑中心的斜距差補償對回波多普勒項的近似。于是M個散射粒子的總回波信號可以近似表示為

其中RE和Rm分別為虛擬散射粒子和真實散射粒子在SAR孔徑中心時刻的斜距，定義為虛擬散射粒子的等效散射系數，表示為

式(4)表示M個散射粒子的總回波信號用位于E點的虛擬散射粒子的回波信號進行等效。由于虛擬散射粒子的個數小于真實散射粒子的個數，因此SAR回波模擬的效率幾乎提高了M倍，并且可以發現，散射粒子空間分布越集中、密度越大，效率提高越顯著。

4 基于等效散射的快速模擬方法

由于3維森林冠層內散射介電粒子成簇聚集分布，十分符合等效散射處理的特點。因此，本文根據等效散射處理方法，通過在3維森林冠層內構造一系列虛擬散射粒子，將3維森林冠層總回波信號用虛擬散射粒子的回波信號等效，達到提高SAR森林遙感回波模擬效率的目的。該方法在雷達后向散射計算后、SAR回波數據生成前進行，與電磁波散射層面相對獨立，不影響森林電磁波散射的物理描述，這意味著森林冠層完整的散射信息將得到保留，它們是被“等效”而不是簡單的“忽略”。

4.1 3維森林場景切分步驟

本小節將給出等效散射處理步驟，重點在于虛擬等效散射粒子的位置確定方法。將整個3維森林場景按照如下步驟進行切分，切分效果如圖4所示。

圖4 等效散射處理效果示意圖

(1)將3維森林場景沿方位向(Y軸)劃分子場景，子場景的寬度為Dy。

(2)對于每一個子場景沿地距向和高度向劃分“場景塊”，“場景塊”的寬度為Dx，高度為Dz。

(3)在孔徑中心時刻，將雷達與“場景塊”中心的視線方向定義為r方向，r與雷達到“場景塊”的視角q隨“場景塊”的空間位置不同而發生變化。在“場景塊”內沿r方向的垂線方向劃分“場景帶”，寬度為Dr，每個“場景帶”定義為一個虛擬散射單元。

(5)針對所有的子場景按照(2)-(4)步驟進行處理，得到全部虛擬散射粒子的中心位置和等效散射系數。

(6)剔除所有等效散射系數為“零”的虛擬散射粒子。若某個虛擬散射粒子等效散射系數為“零”，則意味著該“場景帶”內不存在真實散射粒子，該虛擬散射粒子無效，對回波信號沒有貢獻，應予以剔除。

(7)用所有有效(等效散射系數“非零”)的虛擬散射粒子代替全部的真實散射粒子，代入回波模擬和成像處理過程，從而得到森林3維場景模擬回波數據和影像數據。

4.2 控制參數的約束條件

SAR回波信號模型包含線性調頻項和多普勒項，等效散射處理后，等效回波的線性調頻項和多普勒項均將引入斜距誤差。為了保持等效散射處理后的回波信號精度，使等效回波信號能夠發揮原始回波信號相同的作用，需要對線性調頻項和多普勒項的誤差進行約束。

從對3維森林場景的等效處理過程可以看出，子場景的寬度Dy、“場景塊”的寬度Dx和高度Dz、“場景帶”寬度Dr是控制虛擬散射粒子位置和等效范圍的重要參數(稱為“控制參數”)，也是控制斜距誤差的決定性參數。因此，本小節將從SAR成像幾何與3維森林場景的等效散射處理過程出發，分別從多普勒項約束、線性調頻項約束以及其它約束3個方面推導控制參數的約束條件

圖5中長方體為場景劃分的某個“場景塊”，其中心Oi的坐標為,qi為“場景塊”中心的雷達視角，ri為孔徑中心時刻雷達與“場景塊”中心的視線方向，si為ir的垂線方向；雷達在方位向(Y正向)隨慢時間的位置變化用,L為合成孔徑長度)表示，雷達斜視角變化為j(u),qa為雷達方位向波束寬度)；A為該“場景塊”內某“場景帶”的中心點(即虛擬等效散射粒子中心)，在rs坐標系中的坐標為(r0,0),B為該“場景帶”內某一真實散射粒子，在rs坐標系中的坐標為，沿方位向與Oi偏差Δy。則有，其中Ds為“場景塊”內一點到ri軸的最遠距離。

圖5 算法控制參數推導示意圖

4.2.1 多普勒項約束條件 根據第3節等效散射處理的原理，將位于B點的真實散射粒子用位于A點的虛擬散射粒子等效后，等效回波信號多普勒項殘余斜距誤差可以寫為

則最大殘余斜距誤差可以寫為

定義kr,ky,ks分別為Dr,Dy,Ds的敏感因子，則有

Dx和Dz的大小與Ds存在直接關系，可以得到

于是，控制參數與等效回波信號多普勒項的殘余斜距誤差應滿足式(11)條件，其中Dopplerr為多普勒項的殘余斜距誤差要求，根據波長和精度需求確定。

由式(9)和式(11)分析可知，由于雷達方位向波束寬度一般取值較小，因此ky相對較大，kr和ks相對較小，圖6給出了3個敏感因子隨波束寬度的變化曲線。可以看出，對斜距誤差貢獻最大的項為Dy，次大的項為Dr,Ds的貢獻最小。式(10)導出的Dx和Dz與局部入射角有關，該角度隨“場景塊”的不同位置而發生變化，由于Ds對誤差很不敏感，因此Ds對Dx和Dz的限制可以統一采用場景中心視角，由此也可以得到均勻的場景切分，使3維場景切分更為簡便。即可將式(11)重新寫為式(12)，式中q0為雷達觀測3維森林場景的中心視角。

圖6 敏感因子分析曲線

4.2.2 線性調頻項約束條件 根據第3節等效散射處理的原理，將位于B點的真實散射粒子用位于A點的虛擬散射粒子等效后，等效回波信號線性調頻項殘余斜距誤差可以寫為

則最大殘余斜距誤差可以寫為

于是，從精度上考慮，控制參數與等效回波信號線性調頻項的斜距誤差應滿足式(15)條件，其中為線性調頻項的殘余斜距誤差要求，根據距離向采樣間隔和精度需求確定。

由式(15)可以看出，由于qa一般較小，Dy項一般可以忽略，因此等效回波信號線性調頻項斜距誤差主要與Dr有直接約束關系，這也與圖5中Dr表示的物理含義(斜距向的切分長度)一致。

4.2.3 其它約束條件 除了受回波信號多普勒項和線性調頻項斜距誤差的限制，切分3維森林場景的各控制參數還受方位向采樣和場景大小的限制。

SAR方位向采樣間隔為V/ PRF ，其中V為飛行速度，PRF為脈沖重復頻率。對森林場景等效散射處理后，若真實散射粒子與虛擬散射粒子在回波信號的方位采樣點不同，不僅會影響回波信號的精度，而且影響干涉配準，降低了本方法在干涉模擬中的應用精度。因此，應限制子場景的寬度Dy在半個方位采樣間隔之內，使虛擬散射粒子與被其等效的真實散射粒子位于相同的方位采樣內，即

另外，對于“場景塊”的寬度Dx和高度Dz，其取值應限制在整個3維森林場景的范圍內，即

式中Wx表示整個森林場景的地距向寬度，Wz表示整個森林場景最大高度。

4.3 控制參數優化組合準則

綜合多普勒項、線性調頻項以及其它條件，得出3維森林場景等效散射處理的控制參數與斜距誤差、方位向采樣和場景大小的約束關系：

其中多普勒項斜距誤差對Dy的約束相對較嚴格，線性調頻項斜距誤差對Dr的約束相對較嚴格。

5 快速模擬方法驗證

本節從數據域復相關的角度對本文提出的快速模擬方法進行驗證，并從模擬效率和模擬精度兩方面對應用效果進行評估。復相關通過式(19)來表示，式中s1和s2分別表示經等效散射的快速模擬方法和傳統模擬方法得到的數據(回波數據或影像數據)。數據相關性驗證在表1的參數下進行，選擇隨機分布的闊葉林作為模擬樣地。3維森林場景的切分控制參數如表2所示。

根據表2參數對3維森林場景進行等效散射處理，模擬得到等效回波信號(圖7(a), 7(b), 7(c))，并通過CS成像算法[14]處理得到模擬影像(圖8(a), 8(b),8(c))。為了對比應用效果，本文同時給出未經過等效散射處理(以下稱為傳統方法)的回波信號模擬結果(圖7(d), 7(e), 7(f))和影像模擬結果(圖8(d), 8(e),8(f))。通過將等效模擬數據和傳統模擬數據進行對比發現，相同極化方式的數據從定性角度看幾乎沒有差別。

表1 快速模擬方法驗證模擬參數表

表2 等效散射處理控制參數表

為了進一步驗證本文方法的應用效果，以下從定量角度對模擬效率和模擬精度進行分析。在模擬效率方面，分別統計本文提出方法與傳統方法的散射粒子數量和模擬時間，統計結果如表3所示；在模擬精度方面，將本文提出方法與原始方法的模擬結果分別在回波數據域和影像數據域進行復相關處理，相關系數越接近 1，則說明二者一致性越好，復相關結果如表4所示。

表3 模擬效率對比結果

表4 數據復相關結果

模擬效率的對比結果表明，經等效散射處理后，散射粒子數量降低了一個數量級，使SAR回波信號模擬耗時降低了一個數量級，模擬效率提高了近10倍；復相關結果表明，無論從回波數據域還是影像數據域，本文提出方法的模擬結果與傳統模擬結果復相關程度十分接近 1，二者從幅度和相位兩方面都具有一致性，這說明經等效散射處理后的模擬結果保持了較高的精度。因此，從效率和精度兩方面看，本文提出方法應用效果明顯，實現了3維森林場景的極化SAR數據快速模擬。鑒于作為比較對象的原始方法計算復雜度較高，因此本文的應用實例是在低密度林分給出的，當林分密度進一步增大，森林冠層散射粒子密度也隨之增高，則本文提出的快速模擬方法更加有效。

6 應用實例

本節利用本文提出的快速模擬方法模擬1.0×106m2闊葉林極化干涉影像對，采用典型的森林高度提取算法對模擬影像進行處理，反演得到森林高度。模擬參數如表5所示。

圖7 SAR模擬回波信號，(a)(b)(c)經等效散射處理結果，(d)(e)(f)未經等效散射處理結果，(a)(d)HH極化，(b)(e)HV極化，(c)(f)VV極化

圖8 SAR模擬影像，(a)(b)(c)經等效散射處理結果，(d)(e)(f)未經等效散射處理結果，(a)(d)HH極化，(b)(e)HV極化，(c)(f)VV極化

模擬數據與樹高反演結果如圖9所示，反演樹高統計直方圖9(c)表明，反演平均樹高為9.9479 m，十分接近平均樹高的模擬設置值10 m。這說明本文提出的快速模擬方法可用于大范圍森林場景的SAR遙感數據模擬，能夠將森林生物學結構信息忠實地反映在模擬數據中，模擬過程是有效的，模擬數據可用于森林生物學參數反演算法的研究。

7 小結

SAR森林遙感數據模擬引入SAR回波模擬過程后計算量較大成為阻礙其應用的瓶頸。本文提出了一種基于等效散射處理的3維森林場景極化SAR數據的快速模擬方法，通過將3維森林場景劃分為“子場景”、“場景塊”和“場景帶”構造出一系列虛擬散射粒子，在斜距誤差較小的情況下將多個真實散射粒子對 SAR回波信號的綜合貢獻用單個虛擬散射粒子等效。本文推導了各等效控制參數之間的被約束關系，分析了各參數的誤差敏感度，得到了模擬精度和效率的優化組合準則。驗證結果表明，該方法在保持較高模擬精度的同時可以將模擬效率提高1個數量級以上，可以有效地解決3維森林場景極化SAR回波模擬計算較大的問題。模擬結果可用于SAR林業遙感系統性能分析、成像驗證以及森林生物信息參數反演算法的相關研究。另外，本方法在玉米、水稻等農作物的極化SAR數據快速模擬方面也具有一定的應用前景。

表5 樹高反演應用模擬參數表

圖9 模擬數據及樹高反演結果

[1] 陳爾學. 合成孔徑雷達森林生物量估測研究進展[J]. 世界林業研究, 1999, (6): 18-23.Chen Er-xue. Development of forest biomass estimation using SAR data[J].World Forestry Research, 1999, (6): 18-23.

[2] Garestier F and Le Toan T. Forest modeling for height inversion using single-baseline InSAR/Pol-InSAR data[J].IEEE Transactions on Geoscience and Remote Sensing, 2010,48(3): 1528-1539.

[3] Frey O and Meier E. 3-D time-domain SAR imaging of a forest using airborne multibaseline data at L- and P-bands [J].IEEE Transactions on Geoscience and Remote Sensing, 2011,49(10): 3660-3664.

[4] Garestier F and Le Toan T. Estimation of the backscatter vertical profile of a pine forest using single baseline P-band(Pol-)InSAR data[J].IEEE Transactions on Geoscience and Remote Sensing, 2010, 48(9): 3340-3348.

[5] 孫國清. 雷達后向散射模型及其在雷達圖像地形影響糾正中的應用[J]. 遙感學報, 2002, (6): 406-411.Sun Guo-qing. Radar backscattering model and its application in terrain-effect reduction of radar imagery[J].Journal of Remote sensing, 2002, (6): 406-411.

[6] Ulaby F T, Sarabandi K, McDonald K,et al.. Michigan microwave canopy scattering model[J].International Journal of Remote Sensing, 1990, 11(7): 1223-1253.

[7] Sun G and Ranson K J. A three-dimensional radar backscatter model of forest canopies[J].IEEE Transactionson Geoscience and Remote Sensing, 1995, 33(2): 372-382.

[8] Williams M L. The theory for a forward SAR model:implementation, applications and challenges[C]. European Conference on Synthetic Aperture Radar, Dresden Germany,May 16-18, 2006.

[9] Xu Feng and Jin Ya-Qiu. Imaging simulation of polarimetric SAR for a comprehensive terrain scene using the mapping and projection algorithm[J].IEEE Transactions on Geoscience and Remote Sensing, 2006, 44(11): 3219-3234.

[10] Liu D W, Sun G Q, Guo Z F,et al.. Three-dimensional coherent radar backscatter model and simulations of scattering phase center of forest canopies[J].IEEE Transactions on Geoscience and Remote Sensing, 2010, 48(1):349-357.

[11] Franceschetti G, Iodice A, Perna S,et al.. SAR sensor trajectory deviations: Fourier domain formulation and extended scene simulation of raw signal[J].IEEE Transactions on Geoscience and Remote Sensing, 2006, 44(9):2323-2334.

[12] Weber J and Penn J. Creation and rendering of realistic trees[C]. Proceedings of the 22nd Annual Conference on Computer Graphics and Interactive Techniques, New York USA, 1995: 119-128.

[13] 張順生. 合成孔徑雷達自然場景回波仿真技術研究[D]. [博士論文], 北京理工大學, 2007.

[14] Raney R K, Runge H, Bamler R,et al.. Precision SAR processing using Chirp scaling[J].IEEE Transactions on Geoscience and Remote Sensing, 1994, 32(4): 786-799.