建筑物的SAR图像仿真技术是根据建筑物的三维模型和雷达的成像参数，仿真电磁波在目标场景中的散射作用，生成逼真的SAR图像的过程^[1]。仿真SAR图像与真实SAR图像在几何特征和辐射特征方面保持一定的一致性^[2]。通过改变目标模型和成像参数，即可简单高效地预测出多角度、不同分辨率的SAR图像，为建筑物的分类、识别及三维重建等研究提供大量的SAR图像。因此对SAR图像仿真技术的研究具有重要的现实意义。

根据后向散射模型和模拟技术的不同，仿真方法可分为两种^[3]。一种是模拟SAR原始回波信号的仿真方法。Franceschett等^[4-6]通过物理光学(PO)^[7]和几何光学(GO)^[8]理论近似电磁波在目标表面的散射过程，得到回波信号；Xu和Jin^[9]采用积分方程模型(IEM)计算电磁散射矩阵，获得回波信号；Delliere等^[10]提出时域有限差分模型(FDTD)近似电磁波的后向散射过程以得到回波信号。这种仿真方法计算复杂，且仿真精度的提高强烈依赖目标建模精度，对建筑物中不可忽视的多次散射难以计算。另一种方法是基于几何特征的仿真方法，通过直接建立目标场景与二维SAR图像间的几何映射关系得到SAR仿真图像。Stefan^[11]设计的RaySAR仿真器用Phong Shading模型计算散射强度。Balz和Stilla^[12]提出的SARVIZ仿真器用光栅法计算几何映射关系。Bolter^[13]提出针对建筑物几何特征的仿真方法，并且对图像几何畸变和斑点噪声的特性进行深入研究。建筑物的外墙、窗户等部件易形成二、三面角结构，在SAR图像中形成明显的几何特征^[14-16]。基于特征的仿真方法不仅能够简洁高效地实现对SAR图像特征的仿真，还能建立目标空间结构与图像散射现象之间的对应关系，成为近年来的研究热点^[17-18]。

本文介绍一种基于图像几何特征的高分辨率SAR图像仿真方法。该方法首先由3D模型和雷达参数构建成像几何模型；其次，通过改进光学着色模型计算电磁波的散射特性；然后根据SAR成像原理改写Pov-Ray(persistence of vision ray)，在雷达与场景间展开射线追踪，建立三维场景与二维图像之间的映射关系，得到电磁波在场景中的散射点的位置、强度等信息；根据强度和位置进行二维成像，形成初步的SAR图像强度图；最后，对强度图进行系统卷积、添加噪声、插值等后处理过程，得到逼真的仿真SAR图像。

1 仿真模型 1.1 成像几何模型

成像几何模型在同一坐标系中描述SAR图像仿真中虚拟雷达的位置以及三维目标的几何位置。SAR系统采用主动侧视斜距成像模型，传感器在同一位置发射和接收电磁波。在同一场景范围内，雷达的航向和飞行高度不变。在远场的条件下，认为雷达在局部场景范围内以平行波入射，射线与水平面夹角均为θ。建立局部直角坐标系，如图 1所示。

以雷达飞行方向为Y轴(方位向)，水平面内垂直于飞行方向的为X轴(距离向)，高度向为Z轴。电磁波由虚拟雷达S发出，经路径r₁与目标场景产生第一个交点P₁，若电磁波直接反射回接收传感器，此时产生一次散射回波；若电磁波在P₁发生镜面反射，经路径r₁₂到达建筑物表面P₂，再由路径r₂到达接收器，则此时产生二次散射。

在传感器发射电磁波进行采样时，同一方位向的采样线数目n应满足n≥Xsinθ/ΔR，其中X表示距离向场景范围，ΔR表示距离向分辨率，取等号时，每个分辨单元有一个采样点。图像仿真中，可以根据需求增加每个分辨单元的采样点数，即增加每个单元的散射点数目，提高图像精度。在射线追踪的过程中，假设每条采样线照射在一个散射点内，相邻的分辨单元没有叠加。后续将得到的散射强度图与点扩散函数进行卷积处理，以接近雷达电磁波带宽有限的事实。

1.2 散射模型

由散射模型可以得到仿真SAR图像的强度。本文采用改进的光学着色模型Phong Shading作为电磁散射模型。在Phong Shading模型中，各点的光学着色通过顶点插值计算，逼真度较高，光强由环境光、镜面反射和漫反射光构成。在SAR成像中，雷达接收到环境散射的电磁波忽略不计，即对目标表面的漫反射只计算一次，不追踪漫反射的多次反射。在SAR图像仿真中，电磁波在目标表面的散射特性由镜面反射和朗伯体散射两部分近似，镜面反射和朗伯体散射强度的和值代表图像灰度大小。两种散射的计算公式为：

$ {I_{\rm{d}}} = {k_{\rm{d}}}{I_{{\rm{in}}}}\left( {\mathit{\boldsymbol{N}} \cdot \mathit{\boldsymbol{L}}} \right), $

(1)

$ {I_{\rm{s}}} = {p_i}{k_{\rm{s}}}{\left( {\mathit{\boldsymbol{N}} \cdot \mathit{\boldsymbol{H}}} \right)^n}. $

(2)

式中：I_d表示漫反射能量强度；k_d表示漫反射系数，取值在(0, 1)之间；I_in表示一个分辨单元内信号的入射强度；L表示入射向量；N是入射平面的法向量；I_s表示目标表面发生的镜面反射强度；H是入射向量与视点方向的半角向量，发生镜面反射时，H与N重合；k_s表示镜面反射系数，取值在(0, 1)之间；n是反射光的空间分布因子，n值越大，目标表面的散射空间分布越集中；p_i是控制因子，随着镜面反射次数i的增加，p_i减小，当i达到设定的最大镜面反射追踪次数(本文i最大值为5)时，p_i为0，认为反射回波能量太弱不能被接收到。当回波向量与镜面反射向量之间的夹角Δα小于1°时，认为发生的散射都是镜面散射。镜面散射中，反射向量S的计算公式为

$ \mathit{\boldsymbol{S}}{\rm{ = }}\mathit{\boldsymbol{R}}{\rm{ - 2}}\mathit{\boldsymbol{N}} \cdot \left( {\mathit{\boldsymbol{N}} \cdot \mathit{\boldsymbol{R}}} \right). $

(3)

式中：S表示反射向量，R表示由天线指向目标表面交点P方向的初始射线。

1.3 基于射线追踪的采样模型

射线追踪模拟雷达发射电磁波，跟踪电磁波与目标之间的相互作用过程。由射线追踪建立目标与二维SAR图像之间的映射关系。在仿真中，射线追踪部分借助Pov-Ray光学渲染软件实现。不同于光学成像原理，在SAR成像中，雷达发射与接收电磁波位于同一位置，且成像方式为斜距成像，因此要根据成像原理修改成像程序。具体流程如下：

1) 在同一坐标系下导入目标的三维模型及表面特征、雷达的位置及近距入射角。

2) 根据成像场景及仿真参数确定射线的数量。场景距离向幅宽为X，方位向为Y，分辨率为Δx、Δy，则射线数目为n_x=N·X/Δx, N=1, 2, 3，…，N的取值表示采样时每个分辨单元内散射点数目，N取整数使得仿真图像更加平滑逼真。方位向同理。

3) 在雷达与场景间进行射线追踪。以雷达的位置S为原点，雷达电磁波的入射方向为射线方向生成射线族，对每条射线的追踪结果用一个数组P_n=[x_nz_nr_nI_nb_nf_nX_nY_nZ_n]记录。y_n、z_n和r_n表示采样点的回波相位中心方位向、高度向及斜距距离坐标，I_n是雷达接收到的信号强度，由散射模型计算，b_n记录反射次数，f_n标记是否发生镜面散射，X_n, Y_n, Z_n为射线与场景目标第一次相交的交点坐标。以射线在目标表面发生二次散射为例，如图 1所示，射线在场景中的交点为P₁(x₁, y₁, z₁)和P₂(x₂, y₂, z₂)，根据斜距成像的原理，射线的回波相位中心距离向坐标r_n为射线路径r₁、r₁₂及r₂路径之和的一半，更多次散射计算同理。方位向x_n及高度向z_n的计算也采取各交点坐标平均的方式，即

$ \left\{ \begin{array}{l} {y_n} = \frac{{{y_1} + {y_2}}}{2}, \\ {z_n} = \frac{{{z_1} + {z_2}}}{2}, \\ {r_n} = \frac{{\left| {{\mathit{\boldsymbol{r}}_1}} \right| + \left| {{\mathit{\boldsymbol{r}}_{12}}} \right| + \left| {{\mathit{\boldsymbol{r}}_2}} \right|}}{2}. \end{array} \right. $

(4)

每发生一次镜面反射，f_n置1，b_n值加1。当射线不再与场景有交点或反射次数达到最大设定值、回波能量达到最小阈值，对本射线的追踪停止。

1.4 后处理模型

在采样的过程中，认为电磁波的能量照射到一个分辨单元内，信号被聚焦为一个像素点，相邻的分辨单元内没有相互影响，在这种假设下得到的是带宽无限的理想反射强度图。此时，需要引入点扩散函数与强度图进行卷积，点扩散函数选用Hamming窗函数：

$ \begin{array}{l} W\left( {n, \alpha } \right){\rm{ = }}\alpha {\rm{ + }}\left( {1{\rm{ - }}\alpha } \right)\cos \left( {\frac{{2\pi n}}{{N - 1}}} \right), \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;0 \le n \le N - 1. \end{array} $

(5)

式中：N表示信号长度；α的值决定窗的陡峭程度，工程实践中一般取α=0.54。将散射强度图进行傅里叶变换，在频域对二维的散射强度图和Hamming窗函数进行卷积，最后再将所得结果进行傅里叶逆变换，得到空域的散射强度图。

另外，采样过程选取的是像素点的散射模型，在将采样点的坐标按斜距成像映射到成像平面时，会出现多个采样点映射到单个分辨单元，也会出现有些分辨单元没有映射的情况。这便导致一次散射图像中易出现明暗相间的条纹特征。增加采样率可以减缓条纹，但采样率增加使得运算量大幅增加。故本文采用在一次反射矩阵中插值的方法来解决。插值选用三次样条插值，保证数据有良好的稳定性和收敛性，减少仿真误差的引入。

仿真过程采用的是非相干方式的SAR仿真方法，得到的散射强度图没有噪声，然而斑点噪声是SAR图像的一个显著特征，需要对图像添加噪声以接近仿真图像。相干斑在统计意义上可以用乘性噪声模型进行描述，本文采用瑞利相干斑模型为图像添加乘性噪声。瑞利噪声的二维分布密度函数为

$ {p_R}\left( {x, y} \right) = \left\{ \begin{array}{l} \frac{x}{{\mu _1^2}}{{\rm{e}}^{\frac{{ - {x^2}}}{{2\mu _1^2}}}} \cdot \frac{y}{{\mu _2^2}}{{\rm{e}}^{\frac{{ - {y^2}}}{{2\mu _2^2}}}}, \;\;\;\;\;x \ge 0, \\ 0, \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;x < 0. \end{array} \right. $

(6)

式中：μ₁、μ₂分别为图像在距离向和方位向的强度均值。

对于强度图中的每一个分辨单元，首先产生一个服从瑞利分布的均值为1的随机数，得到噪声的强度值。然后将分辨单元的强度值与噪声的值相乘，即可得到本分辨单元带有斑点噪声的散射强度值，最终得到含有斑点噪声的仿真SAR图像，计算公式为

$ \Delta \hat y\left( {x, y} \right) = \gamma \left( {x, y} \right) * n\left( {x, y} \right). $

(7)

式中：γ(x, y)表示散射强度图中(x, y)处的散射强度值；n(x, y)是对应的噪声强度；$\hat \gamma $(x, y)是添加噪声后的像素点强度。

2 仿真方法

基于图像几何特征的高分辨率SAR图像仿真流程如图 2所示。

1) 首先，建立仿真目标的三维模型，借助图形学的建模软件如Auto CAD、3DS Max等制作完成。然后根据雷达工作参数建立局部O-XYZ直角坐标系，明确雷达的飞行高度、飞行方向、姿态角、分辨率、目标的三维坐标，以及方位向距离向与X轴、Y轴的对应关系。至此，成像几何模型建立。

2) 其次，在雷达与目标场景间展开射线追踪，仿真电磁波的反射过程。射线追踪结束，输出一个多行的矩阵，矩阵的每一行代表一个采样点的信息，包含采样点的位置、强度、散射类型、散射次数等。根据散射次数的不同，将采样点信息矩阵分层存储为一次散射矩阵、二次散射矩阵等。

3) 然后，根据雷达成像原理，生成散射强度图。散射强度图的高度N_H和宽度N_W的计算公式为

$ \left\{ \begin{array}{l} {N_{\rm{H}}} = \frac{X}{{\Delta x}} + 1, \\ {N_{\rm{W}}} = \frac{Y}{{\Delta y}} + 1, \end{array} \right. $

(8)

$ \left\{ \begin{array}{l} {N_{{\rm{W}}n}} = \frac{{{N_{\rm{W}}}}}{{{r_{\max }} - {r_{\min }}}} \cdot \left( {{r_n} - {r_{\min }}} \right), \\ {N_{{\rm{H}}n}} = {y_n}. \end{array} \right. $

(9)

根据斜距r_n和方位y_n求得散射点在散射强度图中的位置，由于雷达工作在正侧视条件下，故方位向坐标N_Hn=y_n。将散射强度I_n映射到空白图像中对应的分辨单元，当有多个回波能量映射到同一个分辨单元时，对这些能量进行非相干叠加。最后对整幅图像进行灰度归一化处理，得到各层的散射强度图。也可以根据发生镜面反射与否作为标志，生成镜面散射图与非镜面散射图。

4) 最后，对散射强度图进行后处理。根据后处理模型，对一次散射矩阵进行三次样条插值，再将各散射层级的图像叠加，得到全层级的散射图像；接着在频域对点扩散函数与散射图进行卷积，傅里叶逆变换得到时域图像，在图像中添加乘性噪声。最终，得到高分辨率的SAR仿真图像。

3 仿真结果及分析

为验证本文所提方法的有效性，选取4组具有代表性的目标进行实验。前3组实验中，雷达的工作参数如表 1所示。

第1组实验以三面角反射器为目标，仿真模型及结果如图 3所示。角反射器对入射电磁波有汇聚作用，在建筑物目标中，窗台、窗框及与地面垂直的墙壁等结构极易构成三面角反射器。本文选取边长20 cm的角反射器作为测试目标。图 3(b)、3(c)分别是单个角反射器的仿真SAR图像与实测SAR图像。从仿真图像同实际SAR图像的对比可得，三面角反射器在SAR图像中以点目标的形式呈现，表现为一个边缘扩散的高亮度分辨单元，从像素级别验证了本文所提方法的有效性。

第2组实验选取简易的平顶房屋作为三维模型，并将仿真SAR图像与理论分析图像对比。图 4(a)表示一个长宽高为30 m×20 m×10 m的平顶屋，雷达飞行方向与屋顶的长边方向一致，仿真结果如图 4(b)所示，一次散射由屋顶和墙壁区域形成，由于斜距成像的原因，墙壁和屋顶产生叠掩。墙壁与地面垂直形成二面角，仿真图像上电磁波的相位中心位于墙角位置，回波叠加后出现二次散射亮线。电磁波照射不到的地方形成阴影区。对比图 4(c)可得，叠掩、阴影、顶底倒置及二次散射的亮线等几何特征均与理论分析结果一致，证明了本方法能较好的模拟简单建筑的几何特征。

第3组实验选取带窗户的房屋和有栏杆的筒仓模型进行实验。图 5(a)表示10 m×10 m×5 m的带窗房屋，包含窗户、窗框等元素，方位角为45°。图 5(b)表示筒仓的三维模型，外墙曲面以三角平面近似。图 6表示带窗户房屋的全散射和各次散射图。一次散射形成图像轮廓，二次散射发生在墙角和窗框部分，三次及以上的散射只发生在窗户部分。在多次散射中，斜距等于电磁波传播总路程的一半，故窗框四次、五次散射的位置其比一次散射的位置靠右。墙面每侧的3个亮区由3个窗户散射形成。

对筒仓模型的仿真结果如图 7所示，图 7(a)其中图 7(b)为散射层级分布图。无论是平面建筑还是曲面建筑，仿真结果都体现了叠掩、阴影、多次散射等几何特征。对于窗户、栏杆等建筑细节的仿真，说明本方法在高分辨率仿真中具有可行性。

第4组实验选取多伦多地区某组合建筑物为目标，其光学航拍图像及三维模型如图 8所示。在实验中，雷达的中心频率为9.65 GHz，带宽为0.15 GHz，俯仰角为46°，采样间隔为1.1 m，组合建筑的三维模型由Google Earth数据获取，雷达入射角与建筑物侧边的夹角为30°。

采用本文算法进行仿真，得到仿真图像如图 9(a)所示，由文献[1]算法得到仿真结果如图 9(b)所示，将仿真结果分别与实测图像图 9(c)进行对比。

仿真图像中建筑物的几何特征均能与实测SAR图像保持一致。由于文献[1]中对分辨单元内散射点的处理不充分，没有充分考虑散射模型、噪声及带宽等因素，仿真图像中出现了明暗相间的条纹特征，对比实测SAR图像可得，本文的仿真结果更加逼真。但真实SAR图像中有部分亮点在本文算法的仿真SAR图像中没有对应，是因为对建筑物的建模精度不够，阁楼、窗户、金属空调架以及屋顶的通风口等建筑构件在建模中被忽略了。因此，提高建模精度是得到逼真SAR图像的关键因素。

综合以上实验结果可得，针对规则的建筑物目标，本文的仿真方法可以得到几何特征一致的逼真SAR图像。

4 结论

本文提出一种基于建筑物几何特征的高分辨率SAR图像仿真方法，该方法首先根据三维场景模型和雷达的工作参数构建成像几何模型，然后基于改进的Phong Shading散射模型，在场景间展开射线追踪，得到采样点的散射信息矩阵，接着将得到的采样点数据进行后处理，最终得到高分辨率的仿真SAR图像。分析和仿真结果表明，该方法能够准确、高效的仿真建筑物在SAR图像中的几何特征，后处理过程严格考虑了SAR实际成像过程中的带宽、噪声等影响因素，另外，分层的数据处理方法为后续研究散射特征提供了基础。下一步将研究仿真图像在图像建筑物提取及目标重建中的作用。