1 引言

太赫兹频段作为电磁波谱中尚未被充分开发的区域，具有与微波和红外频段不同的特性，在生物、医疗、通信、成像等领域显现出广阔应用前景^[1–4]，因此太赫兹正受到越来越多研究人员的关注。太赫兹雷达作为一项重要的太赫兹应用技术，在民用与军事领域都具有巨大应用潜力。当前太赫兹雷达技术正处在如火如荼的研究热潮之中，美国喷气推进实验室(Jet Propulsion Laboratory, JPL)报道了其研制出的世界上第1部高分辨率太赫兹成像系统^[5]，并一直致力于对已有系统的改进与升级^[6–8]，德国弗朗科夫(Fraunhofer Institute, FGAN)实验室研制了一部工作频率为220 GHz的太赫兹成像雷达COBRA-220，能获取更多目标细节信息^[9,10]。国内方面，中科院研制了基于220 GHz雷达的成像系统^[11]，电子科技大学研发了330 GHz频段逆合成孔径(Inverse Synthetic Aperture Radar, ISAR)成像系统^[12,13]，工程物理研究院分别研发了140 GHz频段ISAR成像雷达^[14]和340 GHz站开式扫描成像系统^[15,16]。以上结果均展示了太赫兹雷达在成像应用中的巨大优势与潜力。

更高的频率使太赫兹雷达具有更高的分辨能力，同时更短的波长使目标散射表现出新的特点规律。美国亚毫米波实验室(Submillimeter-wave Techniques Laboratory, STL)^[17]研究了不同粗糙度与表面细微结构表面金属圆柱对太赫兹雷达图像的影响，国防科大^[18–20]利用数值仿真手段研究了周期粗糙表面和随机粗糙表面目标在太赫兹雷达图像中的表现规律。这些研究均表明，目标太赫兹雷达图像与微波雷达图像的特征具有显著差异。造成这一差异的一项重要原因是，更短的波长使目标表面粗糙起伏的电尺寸大大增加，目标表面由光滑过渡为粗糙，从而造成散射特性的变化。正是这种粗糙特性使得太赫兹雷达图像中目标的细节和轮廓信息更加丰富，为基于成像的后续应用提供了更多可能。但另一方面，这也给对目标进行散射建模和回波仿真造成了更大的困难。

当前，关于太赫兹频段目标粗糙表面对散射特性的影响已引起了学者们的广泛关注。天津大学以实验方式对粗糙目标在太赫兹频段的散射特性展开了系列研究^[21–24]，文献[25]从通信的角度研究了室内墙壁粗糙面对太赫兹波散射的影响。这些研究成果充分说明了微波频段下往往忽略不计的目标粗糙表面能对太赫兹波散射产生显著影响。近年来在太赫兹频段目标散射特性的计算方面亦取得了不少进展，文献[26]利用高频近似方法对太赫兹频段目标散射进行建模计算，正如其作者指出，目标表面的粗糙是影响目标散射特性的重要因素，而文献[26]中尚未涉及这一问题。北京航空航天大学亦提出了多种目标THz频段雷达散射截面(Radar Cross Section, RCS)的快速计算方法^[27–29]，并进一步对表面粗糙和表面涂覆等因素进行了研究^[30]。文献[31]亦对太赫兹频段粗糙目标的散射计算问题进行了研究。然而，在现有研究中，考虑粗糙表面影响时，研究的对象均为目标RCS，即表征目标散射回波能力强弱的实数值。从雷达探测的角度看，为了成像以及后续应用的需要，除了需要散射回波的强度信息，更关键的是要获取回波的相位信息。但现有针对太赫兹频段粗糙目标的研究尚无提供相位信息的能力，因而亦无法支撑后续成像的需求。

要解决这一问题，为太赫兹雷达成像体制设计、算法开发提供带有相位信息的回波数据。一方面需要有对超电大尺寸目标的处理能力，另一方面要将目标表面粗糙特性考虑在内。综合考虑这两方面因素，本文提出了一种基于面片分级的半确定性建模方法，并利用粗糙面全波法(Full-Wave Algorithm, FWA)进行散射场计算，通过小尺寸目标计算结果与高频数值方法的对比，验证了本文方法的有效性，并给出了全尺寸粗糙锥体的成像结果。初步解决了THz频段全尺寸凸体粗糙目标散射建模及回波生成问题。

2 回波散射建模方法 2.1 面片分级建模

太赫兹频段的波长范围是30 μm～3 mm，对于波长为1 mm的太赫兹波，一个1 m×1 m的目标，电尺寸已达到1000×1000 ${\lambda ^2}$ 。从电磁计算的角度，这一规模不论对于精确数值方法还是高频近似方法都是难以承受的。面片的数量直接决定了矩量法中未知数的个数，也直接决定了基于射线追踪的高频近似方法射线追踪的复杂度，可见巨大的面片数量是造成计算难以实现的直接原因。因此传统的面片剖分成为制约太赫兹频段全尺寸目标散射计算的重要因素。针对太赫兹频段的散射建模问题，受到激光RCS计算方法的启迪，我们提出了太赫兹频段下基于面片分级的半确定性建模方法。

面片分级建模的思路是：首先利用宏观小微观大面片对目标进行建模并构成第1级面片，然后根据目标表面粗糙参数的先验信息，在每个宏观小微观大面片上进行进一步划分，形成更精细的第2级面片。上面提到的“宏观小微观大”需要满足以下两个要求：在宏观上能够较精确地描述目标的形状特征；在微观上能够体现出表面粗糙的统计特征。图1展示了圆锥、球、圆柱3种全尺寸简单体的一级面片划分结果。

二级面片划分的目的是使全尺寸目标体现表面粗糙特征，为了降低算法的内存开销，可以采用如下方法建立二级面片：首先根据建模对象表面粗糙度的先验信息生成粗糙模板样片，该模板面片尺寸须大于所有一级面片尺寸的最大值。其次，依据一级面片的形状，从粗糙样片上抠出与其形状相符的区域。二级面片的建模过程如图2所示。

利用两级面片的建模方法，可以根据目标表面粗糙的先验信息，对不同粗糙参数的区域生成相应的二级面片模板，只需存储少量的模板面片信息，从而大大降低了计算所需的存储量。简单地说，第1级面片用于勾勒目标形状，而不考虑目标表面粗糙，第2级面片则用于在一级面片之上精细地描绘粗糙面起伏形状，从而使目标具有表面粗糙的特征。通过这两级面片模型，既能充分描绘目标形状的几何信息，又能反映目标表面的粗糙特征。

在建模过程中，一级面片的建模(图1)是确定性的，可利用经典电磁计算方法中所用的网格剖分方法，只是此时的面片尺寸并非小于 $\lambda /8$ ，而是远大于波长，本文要求一级面片边长大于 $50\lambda $ 。二级面片建模是半确定性的，可以将粗糙面片看作一个2维随机过程。其与真实目标表面具有相同的统计参数，而生成二级模板面片的过程即是对该随机过程进行了一次实例采样，生成后的二级模板面片在后续建模中不再变化，以此来保证所生成雷达回波的相干性。

2.2 粗糙面散射场计算

通过前面介绍的面片分级建模方法，对整个全尺寸目标的散射计算问题可分解为对每个一级面片的散射计算问题。计算得到一级面片的散射场后经过坐标变换与相位补偿可以得到该面片在目标全局坐标系中的散射场。将所有一级面片的散射场相干叠加便得到了目标总散射场。因此，本节针对单个一级面片研究其散射场计算方法。

粗糙面散射问题其实并不是一个新问题，经典的粗糙面理论^[32]就是针对粗糙面散射问题的一整套理论与方法。其中典型的粗糙面散射理论包括基尔霍夫近似法、微扰法等^[33]。尽管如此，从雷达回波仿真的角度来看，经典粗糙面理论却难以胜任这项任务。其中核心问题在于粗糙面理论的研究目标是提供粗糙面的散射系数，即从能量角度刻画粗糙面的散射能力强弱，而如何提供雷达回波仿真中所需的相位信息并不在其研究范围。

粗糙面全波法最初由Ezekiel Bahar在上世纪80年代提出^[34]，该方法从Maxwell方程组和边界条件出发，较严格地推导了粗糙面散射场(包含相位信息)的表达式，且该方法对粗糙面高度起伏均方根和相关长度具有更大的适用范围。本文我们选取该方法计算粗糙面的散射场。为讨论方便，定义一级面片坐标系如图3所示。关于后文公式中的变量书写，本文做如下规定：标量用斜体表示，如k；矢量由斜体加矢量标识表示，如 ${ \mathop E \limits^{\raisebox{-10pt}{\rightharpoonup}}} $ ， $ \mathop r\limits^{\rightharpoonup} $ ，单位矢量加尖角标识，如 $\hat e$ ；矩阵用斜体表示，如F。

根据粗糙面全波法^[34]，粗糙面散射场可表示为：

(1)

其中，T ⁱ, T ^s分别为入射场、散射场的极化旋转矩阵，F为散射系数矩阵， ${\hat e^{\rm{i}}}$ , ${\hat e^{\rm{s}}}$ 分别为入射、散射方向的单位矢量， ${k_0} = \omega \sqrt {{\mu _0}/{\varepsilon _0}} $ 为自由空间波数， $\omega $ 为电磁波角频率， $ U\left( {\mathop r\limits^{\rightharpoonup}}' \right)$ 为遮蔽函数，S为当前一级面片中包含的全体二级面片集合。入射与散射场分别用极化分量表示如下：

场的极化矢量和极化旋转矩阵定义依赖于一、二级面片坐标系的定义，一、二级面片坐标轴的关系如下：

其中， $ \hat x$ , $ \hat y$ , $ \hat z$ 分别为图3中各坐标轴的单位矢量， $ \hat n$ 为二级面片的单位法矢量，其指向随着二级面片位置矢量 $ {{\mathop r\limits^{\rightharpoonup}}} '$ 的不同而不同，下标n代表与法矢 $ \hat n$ 对应的二级面片，利用式(3)的二级面片坐标系定义，入射与散射场的水平极化单位矢量在一、二级面片坐标系中的定义为：

根据式(4)中对一、二级面片坐标系中极化矢量的定义，可以得出极化旋转矩阵T ⁱ, T ^s的表达式为：

其中， $C_n\,\!\!^{\rm{i}} = {\hat h^{\rm{i}}} \cdot \hat h_n^{\rm{i}}$ , $S_n^{\rm{i}} = {\hat e^{\rm{i}}} \cdot {\hat h^{\rm{i}}} \times \hat h_n^{\rm{i}}$ , $C_n\,\!\!^{\rm{s}} = {\hat h^{\rm{s}}} \cdot \hat h_n^{\rm{s}}$ , $S_n^{\rm{s}} = {\hat e^{\rm{s}}} \cdot {\hat h^{\rm{s}}} \times \hat h_n^{\rm{s}}$ ，可以看出T ⁱ, T ^s对不同的二级面片是不同的。式(1)中散射系数矩阵F的定义如下：

其中，n_r, ${\eta _{\rm{r}}}$ , ${\varepsilon _{\rm{r}}}$ , ${\mu _{\rm{r}}}$ 分别为反射率、相对波阻抗、相对介电系数和相对磁导率，它们可分别由介质0与介质1中电参数表示如下：

$C_n\,\!\!^{\rm{i}}$ , $S_n^{\rm{i}}$ , $C_n\,\!\!^{\rm{s}}$ , $S_n^{\rm{s}}$ 分别代表入射与散射角的余弦与正弦，上标0或1用以区分介质0和介质1，两种介质中入射与散射角的关系可以通过折射定律确定。式(6)中 $\cos \left( {\phi _n^{\rm{s}} - \phi _n^{\rm{i}}} \right)$ , $\sin \left( {\phi _n^{\rm{s}} - \phi _n^{\rm{i}}} \right)$ 分别代表入射与散射平面间夹角的余弦与正弦。

2.3 全尺寸凸体粗糙目标回波仿真

利用上一节介绍的方法，可以得到考虑表面粗糙的单个一级面片的散射场。为得到目标的总散射场，须将一级面片坐标系中的计算结果转换至目标坐标系并将所有一级面片的散射场相干叠加。为此，分别假定目标坐标系中与x, y, z轴平行的单位矢量为 ${\hat x_{\rm{t}}}$ , ${\hat y_{\rm{t}}}$ , ${\hat z_{\rm{t}}}$ ，其中下标t代表目标。与式(4)相似的，可以得到目标坐标系中水平极化方向的单位矢量为：

为了便于在目标坐标系中区分不同的一级面片坐标系，令 $m = 1,2, ·\!·\!· ,M$ 为一级面片的编号，M为一级面片的总数量。于是，将上一节中与第m个一级面片相关的入射场、散射场、入射水平极化矢量、散射水平极化矢量加上下标m，重写为 $ {\mathop E\limits^{\rightharpoonup}} _m^{\rm{i}}$ , $ {\mathop E\limits^{\rightharpoonup}} _m^{\rm{s}}$ , $\hat h_m^{\rm{i}}$ , $\hat h_m^{\rm{s}}$ 。与式(5)相似的，定义目标坐标系到一级面片坐标系的极化旋转矩阵为：

其中， $C_m\,\!\!\!\!^{\rm{i}} = \hat h_{\rm{t}}^{\rm{i}} \cdot \hat h_m^{\rm{i}}$ , $S_m^{\rm{i}} = {\hat e^{\rm{i}}} \cdot \hat h_{\rm{t}}^{\rm{i}} \times \hat h_m^{\rm{i}}$ , $C_m^{\rm{s}} = \hat h_{\rm{t}}^{\rm{s}} \cdot \hat h_m^{\rm{s}}$ , $S_m^{\rm{s}} = {\hat e^{\rm{s}}} \cdot \hat h_{\rm{t}}^{\rm{s}} \times \hat h_m^{\rm{s}}$ 。再次，为了书写方便，将式(1)表示为算子形式：

其中， $A\left( \cdot \right)$ 代表在一级面片坐标系中根据入射场计算散射场的算子。基于以上表述，目标的总散射场可以表示为：

(11)

其中， $ {{\mathop r\limits^{\rightharpoonup}} _m}$ 为第m个一级面片坐标系原点在目标坐标系中的位置矢量。由式(11)，我们就获得了考虑目标表面粗糙的且带有相位信息的散射回波。在雷达成像仿真中，该计算结果可以直接作为雷达的谱域回波使用。

需要说明的是，以上建模过程中为了简化问题，忽略了面片间的互遮挡问题。对于形状较复杂的目标，互遮挡判断往往是不可或缺的一步。为此，本文在后续内容中假定目标是不存在部件互遮挡的简单凸体。下面给出全尺寸目标的散射建模与回波仿真计算流程，如图4所示。

3 回波生成与成像仿真结果分析

本节利用所提方法对不同目标进行散射计算、回波生成与成像仿真，由于内存与计算复杂度的限制，现有电磁计算软件难以在普通个人电脑上处理太赫兹频段全尺寸目标模型，因此本节先以小尺寸模型为研究对象，并通过与电磁计算软件中高频数值方法的计算结果进行对比以验证本文方法的有效性。

为了便于检验所提方法的正确性，本文涉及的成像仿真利用方位俯仰两个方向的转角进行单站方位俯仰成像。在这种观测几何条件下，获得的雷达图像与光学图像具有相同的视角，从而便于对成像结果进行分析。设置仿真频率为300 GHz，对应波长为1 mm。以下仿真中，均假定目标为处于真空背景中的理想导电体(Perfect Electric Conductor, PEC)。与式(7)对应， ${\varepsilon _0}$ , ${\mu _0}$ 分别为真空介电常数与磁导率。假定目标电导率为 $\sigma = {10^{20}}\;{\rm{S}}/{\rm{m}}$ , ${\mu _1} = {\mu _0}$ , ${\varepsilon _1} = {\varepsilon _0} - j\sigma /\omega $ ，其中 $\omega $ 为电磁波角频率。

3.1 小尺寸粗糙矩形板

我们首先选取一个小尺寸矩形板为目标，其尺寸为3 cm×4 cm，对应电尺寸为 $30 \times 40{\lambda ^2}$ 。对2维高斯随机过程进行一次实例采样获得粗糙面片数据，粗糙面高度起伏均方根为 $\lambda /4$ ，相关长度为 $\lambda$ 。我们首先从RCS计算的准确性角度对所提方法进行分析，固定方位角为 $\varphi {\rm{ = }}{30^ \circ }$ ，令俯仰角范围为 $\theta {\rm{ = }}{0^ \circ } \sim {85^ \circ }$ 。基于本文FWA方法的计算结果与利用CST软件中弹跳射线算法(Shooting and Bouncing Ray, SBR)的计算结果分别绘于图5中。

图5展示了4种极化关系下的粗糙矩形板RCS计算结果。由图5(a)与图5(c)可见，HH极化与VV极化在入射角角度小于50°的区间，两者计算的RCS曲线基本重合，接近掠入射时，FWA与SBR计算结果出现一定偏差。由图5(b)和图5(d)可见，FWA计算结果与SBR计算结果的趋势基本相符，在小入射角时偏差稍大。对比同极化和交叉极化散射分量可以发现，同极化分量的能量明显大于交叉极化分量，因此对交叉极化分量计算的偏差对总体RCS计算结果的影响较小。我们认为造成FWA与SBR计算结果在若干角度存在一定偏差的原因有以下几点：①当入射角接近掠入射时，二级面片间的互遮挡现象将十分普遍，而当前实现的FWA算法中并未考虑面片间的互遮挡问题，从而造成了一定的计算误差；②对于交叉极化散射分量，面片间的相互耦合是去极化现象形成的重要原因，而由式(1)可以看出，FWA算法中并未包含射线追踪操作，因此对交叉极化分量的计算存在一定误差。尽管如此，从雷达回波仿真的角度本文所提方法的计算结果仍是令人满意的。

下面进一步对该小尺寸粗糙板进行方位俯仰成像仿真。回波生成及成像仿真参数列于表1中。获得目标的方位俯仰单站散射回波后利用2维积分法实现成像，成像结果如图6所示，图像的两个坐标轴分别代表方位向(Azimuth)和俯仰向(Elevation)。

由图6可以看出，依据FWA和SBR计算数据回波的成像结果基本一致，对比图6(a)与图6(c)、图6(b)与图6(d)可发现，两种方法的同极化成像结果在强点位置、粗糙纹理上都十分相似。交叉极化的成像结果中的纹理细节虽不一样，但仍具有相似的粗糙特征。VV极化和VH极化的成像结果分别与HH极化和HV极化的成像结果是相似的，在此不再展示。

表2展示了利用CST高频求解器和本文方法分别生成图6所用回波数据的内存和时间消耗。计算硬件平台是一台配有Intel i7-6500处理器和8 GB内存的个人电脑。电磁计算软件版本为CST2014，所用电磁求解器为A求解器。本文方法的软件实现平台为Matlab2016A。其中CST软件内存占用仅考虑求解器进程的内存消耗，本文方法的内存消耗以程序运行前与运行峰值时Matlab软件的内存占用差为准。可以看出，两种方法的内存消耗差别不大，这是因为对于该小尺寸目标，两种方法均存储了所有面片信息。而在耗时方面本文方法则具有明显优势，仅为SBR方法的0.8%。图5、图6和表2共同说明了利用本文方法进行粗糙面散射场计算与雷达回波仿真的有效性和高效性。

3.2 全尺寸粗糙锥体

下面对全尺寸粗糙凸体目标进行回波生成与成像仿真。所用锥体模型同图1中所示，锥体高1 m，底面直径为0.5 m。回波生成及成像仿真的具体参数见表3。

利用本文回波生成方法，获得目标的方位俯仰单站散射回波并利用2维积分法成像，成像结果如图7、图8所示。图7、图8中假设锥体表面为高斯粗糙面，相关长度均为 $\lambda $ ，图7中高度起伏均方根为 $\lambda /8$ ，图8中高度起伏均方根为 $\lambda /4$ 。由图可以看出，由于粗糙面散射的影响，太赫兹雷达图像与传统微波图像中目标往往表现为少数孤立的散射点不同，目标体现为相连的面快且具有更丰富的纹理信息。且与同极化图像相比，交叉极化图像中目标的轮廓形状更加凸显。对比图7和图8，当粗糙度增加时，目标的轮廓亦表现得更为清晰，上述现象与文献[19]中的报道一致。

表4中列出了对上述全尺寸锥体进行回波生成的资源消耗，其数据获取方式与表2中数据记录的方法相同。由于该锥体在300 GHz下电尺寸已达 $1.3 \times {10^6}{\lambda ^2}$ ，在前面所提的硬件平台上，受到8 GB内存的限制，已无法在CST软件中建立如此庞大的模型，因而无法进行计算。与其相比，本文采用了面片分级的建模方法，计算中的峰值内存消耗不超过3 GB，这使得在普通计算机上进行太赫兹频段全尺寸目标回波仿真成为可能。利用FWA方法，生成图7、图8成像所需回波数据耗时约为88 h，可以满足后续成像研究的需求。以上结果说明了本文方法用于太赫兹频段全尺寸粗糙目标散射建模与回波生成的有效性。

4 结论

本文提出了一种基于面片分级与粗糙面全波法的THz全尺寸凸体粗糙目标散射建模与雷达回波仿真方法，并获得了全尺寸简单体的雷达仿真图像。在面片分级建模中，一级面片用于勾勒目标形状，二级面片则用于在一级面片之上描绘目标表面粗糙起伏。随后利用粗糙面全波法对目标进行散射场计算，生成带有相位信息的雷达回波。在小尺寸目标上，通过与电磁计算软件中高频数值方法的对比，验证了本文所提方法的计算精度、相位生成能力以及在计算空间与时间消耗方面的优势。后给出了全尺寸锥体的算例，验证了本文方法在有限存储空间和合理时间内提供雷达仿真回波的能力。为进一步研究太赫兹雷达成像体制和成像算法打下了良好基础。另外，本文方法目前将目标形状限定为凸体，后续将通过加入互遮挡判断等操作进一步扩展该方法的适用范围。