一种三维森林场景极化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.