CN106597438B - 基于物理光学的太赫兹雷达回波成像方法 - Google Patents

基于物理光学的太赫兹雷达回波成像方法 Download PDF

Info

Publication number
CN106597438B
CN106597438B CN201611079232.6A CN201611079232A CN106597438B CN 106597438 B CN106597438 B CN 106597438B CN 201611079232 A CN201611079232 A CN 201611079232A CN 106597438 B CN106597438 B CN 106597438B
Authority
CN
China
Prior art keywords
surface element
triangular surface
slow
time
terahertz radar
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201611079232.6A
Other languages
English (en)
Other versions
CN106597438A (zh
Inventor
丁金闪
徐众
梁毅
王敏
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xidian University
Original Assignee
Xidian University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Xidian University filed Critical Xidian University
Priority to CN201611079232.6A priority Critical patent/CN106597438B/zh
Publication of CN106597438A publication Critical patent/CN106597438A/zh
Application granted granted Critical
Publication of CN106597438B publication Critical patent/CN106597438B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
    • G01S13/9004—SAR image acquisition techniques
    • G01S13/9017—SAR image acquisition techniques with time domain processing of the SAR signals in azimuth
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Physics & Mathematics (AREA)
  • Signal Processing (AREA)
  • Electromagnetism (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • General Physics & Mathematics (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明公开了一种基于物理光学的太赫兹雷达回波成像方法，其思路为：确定太赫兹雷达，所述太赫兹雷达包含N个天线，所述太赫兹雷达的检测范围内存在复杂目标，所述复杂目标为电大尺寸目标，然后对复杂目标进行三角面元划分，得到M个三角面元，对N个天线分别进行等间隔采样，得到N个天线的NP个采样点，进而得到NP个经过随机抖动的采样点以及NP个采样点对应的三角面元序列；预先设定W个慢时间，分别得到M个三角面元经过W个慢时间旋转后的可见三角面元区域，以及W个慢时间对应的被照亮可见三角面元区域的时延集合P；依次计算W个慢时间对应的太赫兹雷达回波数据和W个慢时间对应的太赫兹雷达回波数据的最终时域输出S(t)，并计算W个慢时间复杂目标对应的太赫兹雷达回波成像。

Description

基于物理光学的太赫兹雷达回波成像方法
技术领域
本发明涉及太赫兹雷达信号处理技术领域，特别涉及一种基于物理光学的太赫兹太赫兹雷达回波成像方法，是一种太赫兹频段电磁波与复杂目标相互作用的太赫兹雷达回波成像方法，适用于研究和模拟太赫兹频段下合成孔径太赫兹雷达(Synthetic ApertureRadar,SAR)成像的特征及其形成的物理过程。
背景技术
合成孔径雷达(SAR)具有全天时、全天候、远距离和高分辨率成像等特点，合成孔径雷达的一个重要功能是成像，并且这一重要功能已广泛应用于军事和民用领域。
随着合成孔径雷达技术越来越成熟，以及合成孔径雷达成像的分辨率不断提高，人们能够得到大量不同场景下不同分辨率的合成孔径雷达图像，且聚焦效果良好，而对这些良好聚焦的合成孔径雷达图像的解释工作显然成为了一项重要的工作；由于合成孔径雷达成像是实际三维物体在距离和方位维的二维投影，所以合成孔径雷达图像是不“直观”的，尤其对于人造的复杂物体如城市里此起彼伏的建筑等，电磁波在其中的相互作用往往会在合成孔径雷达成像中产生一些令人费解的现象，这将导致用肉眼识别和分析合成孔径雷达成像变得充满挑战并且非常耗时；如图1a所示的光学图像和图1b所示的图1a对应的星载合成孔径雷达图像。
发明内容
针对上述现有技术存在的问题，本发明的目的在于提出一种基于物理光学的太赫兹太赫兹雷达回波成像方法，该种基于物理光学的太赫兹雷达回波成像方法便于操作，是一种能够进行参数化输入的面向复杂目标的太赫兹雷达回波模拟方法，且该基于物理光学的太赫兹雷达回波成像方法能足够精确地模拟电磁波和任意复杂目标相互作用的物理过程并生成对应太赫兹雷达回波，然后进行合成孔径成像，使之能够很好地解释和理解对应的实际SAR成像。
为达到上述技术目的，本发明采用如下技术方案予以实现。
一种基于物理光学的太赫兹雷达回波成像方法，包括以下步骤：
步骤1，确定太赫兹雷达，所述太赫兹雷达包含N个天线，每个天线包含若干个阵元；所述太赫兹雷达的检测范围内存在复杂目标，所述复杂目标为电大尺寸目标，然后对复杂目标进行三角面元划分，得到M个三角面元，以及每个三角面元的三个顶点坐标和每个三角面元的法相向量；N、M分别为大于0的自然数；
步骤2，对N个天线分别进行等间隔采样，得到N个天线的NP个采样点，然后对第i个采样点进行随机抖动，得到第i个经过随机抖动的采样点以及第i个经过随机抖动的采样点对应的三角面元序列，所述第i个经过随机抖动的采样点对应的三角面元序列包含ni个三角面元；
然后令i分别取1至NP，进而得到第1个经过随机抖动的采样点至第NP个经过随机抖动的采样点，以及第1个经过随机抖动的采样点对应的三角面元序列至第NP个经过随机抖动的采样点对应的三角面元序列，并分别记为NP个经过随机抖动的采样点以及NP个采样点对应的三角面元序列，i＝1,2,…,NP，
n1+…+ni+…+nNP＝M，P表示每个天线进行等间隔采样后包含的采样点个数，且P个采样点为在每个天线中选取的阵元个数，P为大于0的自然数；
步骤3，预先设定W个慢时间，然后根据预先设定的W个慢时间对M个三角面元分别进行遮挡照亮判断，分别得到M个三角面元经过W个慢时间旋转后的可见三角面元区域，以及W个慢时间对应的被照亮可见三角面元区域的时延集合P；
步骤4，根据W个慢时间对应的被照亮可见三角面元区域的时延集合P，计算得到W个慢时间对应的太赫兹雷达回波数据；
步骤5，对W个慢时间对应的太赫兹雷达回波数据分别沿快时间t进行匹配滤波处理，得到W个慢时间对应的太赫兹雷达回波数据的最终时域输出S(t)，t表示快时间；
步骤6，根据W个慢时间对应的太赫兹雷达回波数据的最终时域输出S(t)，计算W个慢时间复杂目标对应的太赫兹雷达回波成像。
本发明的有益效果：本发明基于物理光学的太赫兹雷达回波成像方法能足够精确地模拟电磁波和任意复杂目标相互作用的物理过程并生成对应太赫兹雷达回波，然后进行合成孔径成像，使之能够很好地解释和理解对应的实际SAR成像；同时，借助本发明方法能够仿真电磁波在场景中多次反射在雷达图像中的效果，有助于人们更好地解释雷达图像。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1a为光学图像示意图；
图1b为图1a对应的星载合成孔径雷达图像示意图；
图2为本发明的一种基于物理光学的太赫兹雷达回波成像方法流程图；
图3a为3ds max软件下二面角的三角面元模型示意图；
图3b为使用本发明方法对3ds max软件下二面角的三角面元模型的成像结果示意图；
图4a为3ds max软件下二面角加遮挡面的三角面元示意图；
图4b为3ds max软件下二面角加遮挡面的三角面元原理示意图；
图4c为使用本发明方法对3ds max软件下二面角加遮挡面的模型的成像结果示意图；
图5a为3ds max软件下城市的三角面元模型示意图；
图5b为使用本发明方法对3ds max软件下城市三角面元模型的成像结果示意图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
参照图2，为本发明的一种基于物理光学的太赫兹雷达回波成像方法流程图；所述基于物理光学的太赫兹雷达回波成像方法，包括以下步骤：
步骤1，确定太赫兹雷达，所述太赫兹雷达包含N个天线，每个天线包含若干个阵元；所述太赫兹雷达的检测范围内存在复杂目标，所述复杂目标为电大尺寸目标，且所述复杂目标位于三维空间坐标系xyz中，复杂目标在y轴上的最大长度为ymax，复杂目标上的所有点在三维坐标系中唯一确定，并将其坐标分别记为列向量，复杂目标的中心在三维坐标系的原点，且复杂目标随着第n个慢时间tn旋转，得到第n个慢时间旋转后的角度θn，n＝1,2,…,W，W表示预先设定的慢时间总个数，
θn＝(n-1)×theta，theta表示太赫兹雷达所在转台角速度，
并且太赫兹雷达发射的电磁波在该电大尺寸目标内部发生至少一次反射后回到太赫兹雷达，所述电大尺寸目标的物理尺寸与太赫兹雷达发射的电磁波的比值大于10；然后利用电磁仿真的有限元技术按照预先设定的精度要求对复杂目标进行三角面元划分，得到M个三角面元，以及每个三角面元的三个顶点坐标和每个三角面元的法相向量，并保存为STL格式文件；N、M分别为大于0的自然数。
本实施例中N＝1，且使用主要利用电磁仿真的有限元技术的有限元分析软件Hypermesh对复杂目标进行三角面元划分；设定使用者要求的成像最小分辨力为fmin，则所述预先设定的精度要求为小于5倍成像最小分辨力为fmin。
步骤2，对N个天线分别进行等间隔采样，得到N个天线的NP个采样点，然后对第i个采样点进行随机抖动，得到第i个经过随机抖动的采样点以及第i个经过随机抖动的采样点对应的三角面元序列，所述第i个经过随机抖动的采样点对应的三角面元序列包含ni个三角面元；然后令i分别取1至NP，进而得到第1个经过随机抖动的采样点至第NP个经过随机抖动的采样点，以及第1个经过随机抖动的采样点对应的三角面元序列至第NP个经过随机抖动的采样点对应的三角面元序列，并分别记为NP个经过随机抖动的采样点以及NP个采样点对应的三角面元序列，i＝1,2,…,NP，n1+…+ni+…+nNP＝M，P表示每个天线进行等间隔采样后包含的采样点个数，且P个采样点为在每个天线中选取的阵元个数，P为大于0的自然数。
所述第i个经过随机抖动的采样点为：iT+ζ，i＝1,2,…,NP，T表示采样间隔，所述采样间隔为太赫兹雷达发射的电磁波载频波长的一半；ζ表示第i个采样点相邻两次采样间隔之间的均匀分布，并且第i个采样点相邻两次采样间隔是不相关的，该方式用以防止欠采样造成的频谱混叠。
步骤3，预先设定W个慢时间，然后根据预先设定的W个慢时间对M个三角面元分别进行遮挡照亮判断，分别得到M个三角面元经过W个慢时间旋转后的可见三角面元区域，以及W个慢时间对应的被照亮可见三角面元区域的时延集合P。
步骤3的子步骤为：
3.1令j表示第j个经过随机抖动的采样点，j∈{1,2,…,NP}，NP表示经过随机抖动的采样点总个数；同时j也表示第j个采样点对应的三角面元序列，NP表示NP个采样点对应的三角面元序列总个数；tn表示第n个慢时间，n的初始值为1，n∈{1,2,…,W}，W表示预先设定的慢时间总个数。
在所述NP个采样点分别对应的NP个可见三角面元序列中，第i个采样点对应若干个三角面元，该若干个三角面元包含不止一个的可见三角面元或不包含可见三角面元，所述不止一个的可见三角面元为第i个采样点对应的可见三角面元序列，i＝1,2,…,NP；同样，如果太赫兹雷达发射的电磁波照射到某个三角面元并发生至少一次反射后回到太赫兹雷达时的可见三角面元，为采样点对应的可见三角面元。
如果第j个采样点对应的三角面元序列为空，则第j个采样点没有对应的三角面元；如果第j个采样点对应的三角面元序列不为空，则第j个采样点存在对应的三角面元；设定第j个采样点对应的三角面元序列包含mj个三角面元，令r∈{1,2,…,mj}，r表示第r个三角面元，r的初始值为1；m1+…+mj+…+mNP＝M。
3.2将第j个经过随机抖动的采样点与太赫兹雷达天线相位中心的连线，作为第j个待判断光线，该第j个待判断光线的方向向量为Dj，第j个待判断光线的相位中心为Pj；并将第j个待判断光线的方向向量Dj，作为第j个太赫兹雷达发射的电磁波入射方向单位矢量sj；j∈{1,2,…,NP}，j的初始值为1。
3.3判断太赫兹雷达发射的电磁波是否照射到第r个三角面元经过第n个慢时间旋转后的三角面元不可见面，第r个三角面元经过第n个慢时间旋转后的三角面元的法相向量所在平面，为第r个三角面元经过第n个慢时间旋转后的三角面元可见面；反之为第r个三角面元经过第n个慢时间旋转后的三角面元不可见面。
将第r个三角面元经过第n个慢时间旋转后的三角面元法相向量nr与第j个太赫兹雷达发射的电磁波入射方向单位矢量sj相点乘，如果该点乘结果小于0，则进行3.4；否则进行3.8。
3.4判断第r个三角面元经过第n个慢时间旋转后的三角面元是否位于太赫兹雷达发射的电磁波传播方向的相反方向；
即将第j个待判断光线的相位中心Pj与第r个三角面元经过第n个慢时间旋转后的三角面元中点连接，得到第j个待判断光线、第r个三角面元经过第n个慢时间旋转后的三角面元的矢量mnjr，所述第j个待判断光线、第r个三角面元经过第n个慢时间旋转后的三角面元的矢量mnjr的方向为第j个待判断光线的相位中心Pj指向第r个三角面元经过第n个慢时间旋转后的三角面元中点，然后将第j个待判断光线、第r个三角面元经过第n个慢时间旋转后的三角面元的矢量mnjr与第j个太赫兹雷达发射的电磁波入射方向单位矢量sj相点乘，若此处得到的点乘结果大于0，则进行3.5；否则转至3.8。
3.5判断第j个太赫兹雷达发射的电磁波入射方向单位矢量sj是否偏离第r个三角面元经过第n个慢时间旋转后的三角面元：
计算第r个三角面元经过第n个慢时间旋转后的三角面元中心到第j个太赫兹雷达发射的电磁波入射方向单位矢量sj的距离，即将第j个待判断光线、第r个三角面元经过第n个慢时间旋转后的三角面元的矢量mnjr与第j个太赫兹雷达发射的电磁波入射方向单位矢量sj相叉乘，得到第j个待判断光线、第r个三角面元经过第n个慢时间旋转后的三角面元的向量模unjr。
如果第j个待判断光线、第r个三角面元经过第n个慢时间旋转后的三角面元的向量模unjr小于K倍第r个三角面元经过第n个慢时间旋转后的三角面元最长边长，则第j个待判断光线、第r个三角面元经过第n个慢时间旋转后的三角面元的向量模unjr没有偏离第r个可见光面元经过第n个慢时间旋转后的三角面元，并进行3.6；否则进行3.8；K表示设定的倍数值，为大于0的自然数。
3.6判断第j个太赫兹雷达发射的电磁波入射方向单位矢量sj是否最终照亮第r个三角面元经过第n个慢时间旋转后的三角面元：
即将第j个待判断光线的相位中心为Pj沿着第j个太赫兹雷达发射的电磁波入射方向单位矢量sj的方向在第r个三角面元经过第n个慢时间旋转后的三角面元所在平面进行投影，得到第j个待判断光线与第r个三角面元经过第n个慢时间旋转后的三角面元所在平面的交点Tnjr，然后根据几何算法的方法判断第j个待判断光线与第r个三角面元经过第n个慢时间旋转后的三角面元所在平面的交点Tnjr是否在第r个三角面元经过第n个慢时间旋转后的三角面元内部。
如果第j个待判断光线与第r个三角面元经过第n个慢时间旋转后的三角面元所在平面的交点Tnjr在第r个三角面元经过第n个慢时间旋转后的三角面元内部，则说明第r个三角面元经过第n个慢时间旋转后的三角面元被照亮，并得到第r个三角面元经过第n个慢时间旋转后的三角面元参数，所述第r个三角面元经过第n个慢时间旋转后的三角面元参数为第r个三角面元经过第n个慢时间旋转后的三角面元序号、第r个三角面元经过第n个慢时间旋转后的三角面元各自顶点坐标、第r个三角面元经过第n个慢时间旋转后的三角面元中点坐标、第r个三角面元经过第n个慢时间旋转后的三角面元法相向量，然后执行3.7。
如果第j个待判断光线与第r个三角面元经过第n个慢时间旋转后的三角面元所在平面的交点Tnjr没有在第r个三角面元经过第n个慢时间旋转后的三角面元内部，则说明第r个三角面元经过第n个慢时间旋转后的三角面元被遮挡，执行3.8。
具体地，设定第r个三角面元经过第n个慢时间旋转后的三角面元顶点分别为pr1、pr2和pr3，并根据第j个待判断光线与第r个三角面元经过第n个慢时间旋转后的三角面元所在平面的交点Tnjr，分别计算第一向量F1、第二向量F2和第三向量F3,F1＝pr2-pr1,F2＝Tjr-pr1,F3＝pr3-pr1，如果F2×F1与F3×F1的方向相反，即(F2×F1)·(F3×F1)＜0，则三个顶点坐标pnr1、pnr2、pnr3分别和Tnjr相交；如果(F2×F1)·(F3×F1)>0，则三个顶点坐标pnr1、pnr2、pnr3分别和Tnjr没有相交；第一向量F1和第三向量F3在第二向量F2两侧，(F2×F1)·(F3×F1)＝0，则Tnjr在线段pr1pr2上。
设定第j个待判断光线与第r个三角面元经过第n个慢时间旋转后的三角面元相交于点s，连接点s与第r个三角面元经过第n个慢时间旋转后的三角面元中点，得到线段ls，并设定第r个三角面元经过第n个慢时间旋转后的三角面元三条边分别为L1、L2、L3，如果点s在第r个三角面元经过第n个慢时间旋转后的三角面元内部，则线段ls与L1、L2、L3都不相交；如果点s在第r个三角面元经过第n个慢时间旋转后的三角面元外部，因此根据线段ls与L1、L2、L3的相交性确定点s是否在第r个三角面元经过第n个慢时间旋转后的三角面元内部，进而获知第r个三角面元经过第n个慢时间旋转后的三角面元是否被照亮。
3.7光线反射照亮判断：如果第j个待判断光线照射到第r个三角面元经过第n个慢时间旋转后的三角面元并发生反射，则分别计算第j个待判断光线的优化相位中心P'j和第j个太赫兹雷达发射的电磁波入射方向优化单位矢量s'j；所述第j个待判断光线的优化相位中心P'j为第r个三角面元经过第n个慢时间旋转后的三角面元中点，所述第j个太赫兹雷达发射的电磁波入射方向优化单位矢量s'j为：
s'j＝[(1+sj)·nr]×nr+sj，并将第j个待判断光线的优化相位中心P'j和第j个太赫兹雷达发射的电磁波入射方向优化单位矢量s'j分别替换第j个待判断光线的相位中心Pj和第j个太赫兹雷达发射的电磁波入射方向单位矢量sj，返回子步骤3.3；·表示点乘。
3.8令r加1，返回3.3，直到得到第mj'个可见三角面元经过第n个慢时间旋转后的三角面元的参数，并将此时得到第1个可见三角面元经过第n个慢时间旋转后的三角面元的参数至第mj'个可见三角面元经过第n个慢时间旋转后的三角面元的参数，作为第n个慢时间、第j个采样点对应的可见三角面元序列包含的mj'个被照亮的可见三角面元集合；mj'≤mj。
3.9令j加1，返回3.3，直到得到第n个慢时间、第NP个采样点对应的可见三角面元序列包含的mNP'个被照亮的可见三角面元集合,并将此时得到的第n个慢时间、第1个采样点对应的可见三角面元序列包含的m1'个被照亮的可见三角面元集合至第n个慢时间、第NP个采样点对应的可见三角面元序列包含的mNP'个被照亮的可见三角面元集合，作为第n个慢时间、NP个经过随机抖动的采样点对应的被照亮可见三角面元区域，m'1+…+m'j+…+m'NP≤M，并分别将第n个慢时间对应的被照亮可见三角面元区域记为Pn，将第n个慢时间对应的被照亮可见三角面元区域Pn中第i'个可见三角面元的路程记为则得到第n个慢时间对应的被照亮可见三角面元区域中第i'个可见三角面元的时延tni'，C表示光速，i'＝1,2,…,mNP'，进而得到第n个慢时间对应的被照亮可见三角面元区域Pn中的mNP'个时延τn。
3.10令n加1，返回3.3，直到得到第W个慢时间、NP个经过随机抖动的采样点对应的被照亮可见三角面元区域，以及第1个慢时间对应的被照亮可见三角面元区域P1中的m1'个时延至第W个慢时间对应的被照亮可见三角面元区域PW中的mNP'个时延，然后将此时得到的第1个慢时间、NP个经过随机抖动的采样点对应的被照亮可见三角面元区域至第W个慢时间、NP个经过随机抖动的采样点对应的被照亮可见三角面元区域，作为M个三角面元经过W个慢时间旋转后的可见三角面元区域；将第1个慢时间对应的被照亮可见三角面元区域P1中的mNP'个时延τ1至第W个慢时间对应的被照亮可见三角面元区域PW中的mNP'个时延τW，作为W个慢时间对应的被照亮可见三角面元区域的时延集合P，P＝{τ1,τ2,…,τW}。
步骤4，根据W个慢时间对应的被照亮可见三角面元区域的时延集合P，计算得到W个慢时间对应的太赫兹雷达回波数据。
4.1初始化：令τdn'表示第n'个慢时间对应的被照亮可见三角面元区域中第d个被照亮可见三角面元的时延，n'∈{1,2,…,W}，W表示被照亮可见三角面元区域的时延集合对应的慢时间个数；d∈{1,2,…,τn'}，τn'表示第n'个慢时间对应的被照亮可见三角面元区域的时延，Pn'∈{P1,P2,…,PW}；并且确定太赫兹雷达发射线性调频脉冲信号；i'＝1,2,…,mNP'，mNP'表示第NP个采样点对应的可见三角面元序列包含的被照亮可见三角面元的个数；n'的初始值为1。
4.2根据第n'个慢时间对应的被照亮可见三角面元区域的时延，计算得到第n'个慢时间对应的太赫兹雷达回波数据sn'(t)，其表达式为：
其中，τn'i'表示第n'个慢时间、第i'个采样点对应的可见三角面元序列包含的被照亮可见三角面元的时延，rect表示矩形窗函数，t表示时间变量，T表示线性调频脉冲信号的时间宽度，exp表示指数函数，K表示线性调频脉冲信号的调频率，fc表示线性调频脉冲信号载频。
4.3令n'加1,重复子步骤4.2，直到得到第W个慢时间对应的太赫兹雷达回波数据sW(t)，并将此时得到的第1个慢时间对应的太赫兹雷达回波数据s1(t)至第W个慢时间对应的太赫兹雷达回波数据sW(t)，记为W个慢时间对应太赫兹雷达回波数据。
步骤5，对W个慢时间对应的太赫兹雷达回波数据分别沿快时间t进行匹配滤波处理，得到W个慢时间对应的太赫兹雷达回波数据的最终时域输出S(t)，t表示快时间。
具体地，根据W个慢时间对应的太赫兹雷达回波数据，得到M个三角面元分别经过W个慢时间旋转后的可见三角面元区域在距离维的分布；所谓距离维即电磁波从太赫兹雷达到三角面元再反射回太赫兹雷达所走过的路程的一半，例如，第n'个慢时间、第i'个可见三角面元的距离维信息为 表示第n个慢时间对应的被照亮可见三角面元区域Pn中第i'个可见三角面元的路程，C表示光速；n'∈{1,2,…,W}，i'＝1,2,…,mNP'，W表示被照亮可见三角面元区域的时延集合对应的慢时间个数，mNP'表示第NP个采样点对应的可见三角面元序列包含的被照亮的可见三角面元个数，τn'i'表示第n'个慢时间、第i'个可见三角面元的时延，其得到过程为：
首先，计算第n'个慢时间、第i'个采样点对应的可见三角面元序列包含的被照亮的可见三角面元的太赫兹雷达回波数据sn'i'(t)，其表达式为：
该第n'个慢时间、第i'个采样点对应的可见三角面元序列包含的被照亮的可见三角面元的太赫兹雷达回波数据的频谱为Sn'i'(f)，其表达式为：
其次，确定W个慢时间对应的太赫兹雷达回波数据的频域匹配滤波器为H(f)，
其中，f表示W个慢时间对应的太赫兹雷达回波数据的频域匹配滤波器H(f)的频率，K表示设定的倍数值，为大于0的自然数；fc表示线性调频脉冲信号载频；T表示采样间隔，所述采样间隔为太赫兹雷达发射的电磁波载频波长的一半；exp表示指数函数，rect表示矩形窗函数。
将W个慢时间对应的太赫兹雷达回波数据的频域匹配滤波器H(f)和第n'个慢时间、第i'个采样点对应的可见三角面元序列包含的被照亮的可见三角面元的太赫兹雷达回波数据的频谱Sn'i'(f)相乘，得到第n'个慢时间、第i'个采样点对应的可见三角面元序列包含的被照亮的可见三角面元的太赫兹雷达回波数据的最终频谱输出Sn'i'out(f)，其表达式为：
然后，对第n'个慢时间、第i'个采样点对应的可见三角面元序列包含的被照亮的可见三角面元的太赫兹雷达回波数据的最终频谱输出Sn'i'out(f)进行快速傅里叶变换，得到第n'个慢时间、第i'个采样点对应的可见三角面元序列包含的被照亮的可见三角面元的太赫兹雷达回波数据的最终时域输出sn'i'out(t)，其表达式为：
sn'i'out(t)＝|K|Tsinc{KT(t-τn'i')}exp[-j2πfc(t-τn'i')]
第n'个慢时间、第i'个采样点对应的可见三角面元序列包含的被照亮的可见三角面元的太赫兹雷达回波数据的最终时域输出sn'i'out(t)是在快时间t＝τn'i'处出现的sinc函数；令i'分别取1至mNP'，对第n'个慢时间对应的太赫兹雷达回波数据sn'(t)进行匹配滤波，所述第n'个慢时间对应的太赫兹雷达回波数据sn'(t)包含mNP'个采样点对应的可见三角面元序列，进而得到第n'个慢时间对应的太赫兹雷达回波数据的最终时域输出sn'out(t)，在第n'个慢时间对应的太赫兹雷达回波数据的最终时域输出sn'out(t)中第n'个慢时间、第i'个采样点对应的可见三角面元序列包含的被照亮的可见三角面元，分别会形成快时间t＝τn'i'处出现的sinc函数，并且第n'个慢时间对应的太赫兹雷达回波数据的最终时域输出sn'out(t)在不同位置处出现的sinc函数峰体现了τn'i'的信息，从而实现了τn'i'的信息获取。
最后，令n'分别取1至W，分别得到第1个慢时间对应的太赫兹雷达回波数据的最终时域输出s1out(t)至第W个慢时间对应的太赫兹雷达回波数据的最终时域输出sWout(t)，并记为W个慢时间对应的太赫兹雷达回波数据的最终时域输出S(t)，t表示快时间。
步骤6，根据W个慢时间对应的太赫兹雷达回波数据的最终时域输出S(t)，计算W个慢时间复杂目标对应的太赫兹雷达回波成像。
具体地，假设复杂目标的旋转中心在坐标系原点，且复杂目标绕着z轴旋转，太赫兹雷达位于x轴上，且从x轴负半轴朝向复杂目标，并与复杂目标旋转中心即坐标系原点相距R，复杂目标的坐标为(R,0,0)。
6.1初始化：令tn”表示第tn”个慢时间，n”∈{1,2,…,W}，W表示被照亮可见三角面元区域的时延集合对应的慢时间个数，n”的初始值为1；第tn”个慢时间时复杂目标旋转的角度为θ(tn”)，θ(tn”)＝theta×tn”，fmin表示使用者要求的成像最小分辨力，ymax表示复杂目标在y轴上的最大长度。
选取复杂目标上任意一个三角面元的中点，记为点目标e，复杂目标为所有可见三角面元的中点形成的点目标的集合。
6.2计算第tn”个慢时间时点目标e到太赫兹雷达的径向距离为Re(tn”)，其表达式为：Re(tn”)＝R+xecos[θ(tn”)]-yesin[θ(tn”)]；其中，xe表示点目标e在x轴上的位置，ye表示点目标e在y轴上的位置，θ(tn”)表示第tn”个慢时间时复杂目标旋转的角度，θ(tn”)＝θ0+Ωtn”，R表示太赫兹雷达位于x轴负半轴时的位置到复杂目标旋转中心的距离，Ω表示第tn”个慢时间时复杂目标旋转的角速度，tn”表示第tn”个慢时间。
然后计算第tn”个慢时间时点目标e的基带太赫兹雷达回波信号SR(τn”e',tn”)，其表达式为：
其中，τn”e'表示点目标e在第tn”个慢时间对应的太赫兹雷达回波数据经过匹配滤波后得到的sinc函数峰值位置，ρ(xe,ye)表示点目标e的反射密度函数，ρ表示太赫兹雷达图像，xe表示点目标e在x轴上的位置，ye表示点目标e在y轴上的位置，表示太赫兹雷达发射电磁波的载频，Re(tn”)表示第tn”个慢时间时点目标e到太赫兹雷达的距离。
6.3对于复杂目标上点目标e的反射密度函数ρ(xe,ye)，其中点目标e在y轴上的位置ye通过确定，并确定点目标e在y轴上的位置ye包括个位置，且第e'个位置ye'通过确定，τn”e'表示第e'个位置ye'在第tn”个慢时间处对应的sinc函数峰值，yimage为设定的太赫兹雷达图像在y轴大小，fmin表示使用者要求的成像最小分辨力；然后计算点目标e在第e'个位置ye'处的反射密度函数ρ(xe,ye')，分别令进而分别得到点目标e在第1个位置y1处的反射密度函数ρ(xe,y1)至点目标e在第个位置处的反射密度函数记为点目标e最终的太赫兹雷达回波数据成像。
6.4令快时间t＝τn”e'，τn”e'表示第e'个位置ye'在第tn”个慢时间处对应的sinc函数峰值，得到W个慢时间对应的太赫兹雷达回波数据在第e'个位置、第tn”个慢时间处的最终时域输出S(τn”e',tn”)，其表达式为：
其中， 表示太赫兹雷达发射电磁波的载频，Re(tn”)表示第tn”个慢时间时点目标e到太赫兹雷达的距离，ρ(xe,ye')表示点目标e在第e'个位置ye'处的反射密度函数，exp表示指数函数。
令
结合第tn”个慢时间时点目标e到太赫兹雷达的径向距离为Re(tn”)，Re(tn”)＝R+xecos[θ(tn”)]-yesin[θ(tn”)]，得到：
对第tn”个慢时间，由于且该式中的值非常小，因此在不影响精度的情况下令ksinθ(tn”)＝0，k cosθ(tn”)＝k×tn”。
其中， 表示太赫兹雷达发射电磁波的载频；进而计算得到第tn”个慢时间时复杂目标的优化基带太赫兹雷达回波信号其表达式为：
其中，是由于太赫兹雷达位于x轴负半轴时的位置到复杂目标旋转中心的距离R引起的相位项，所以使用的复共轭将去除，该过程在信号处理中为径向运动补偿；表示点目标e在第e'个位置ye'处的优化反射密度函数，其表达式为：
其中，xe表示点目标e在x轴上的位置，IFT表示逆傅里叶变换。
6.5令tn”分别取t1至tW，依次重复子步骤6.2至子步骤6.4，直到得到第t1个慢时间时复杂目标的优化基带太赫兹雷达回波信号至第tW个慢时间时复杂目标的优化基带太赫兹雷达回波信号并对此时得到的第t1个慢时间时复杂目标的优化基带太赫兹雷达回波信号至第tW个慢时间时复杂目标的优化基带太赫兹雷达回波信号进行成像处理，得到W个慢时间复杂目标对应的太赫兹雷达回波成像。
所述W个慢时间复杂目标对应的太赫兹雷达回波成像为复杂目标在距离维和方位维的分布，且所述W个慢时间复杂目标对应的太赫兹雷达回波成像为二维矩阵，组成该二维矩阵的每一列元素分别是复杂目标在距离维的分布，组成该二维矩阵的每一行元素分别是复杂目标在方位维的分布，而复杂目标在方位维的分布垂直于距离维方向的分布。
通过下述仿真实验对本发明效果作进一步验证说明。
(一)仿真参数：
太赫兹雷达的载频fc为1200Thz，太赫兹雷达发射的电磁波带宽为10Thz，太赫兹雷达所在转台角速度theta为0.000001度/慢时间，预先设定的慢时间总个数W为0-300，采样间隔T为2×10-9s，太赫兹雷达位于x轴负半轴时的位置到复杂目标旋转中心的距离R为848m，使用者要求的成像最小分辨力fmin为1.5×10-5m，复杂目标在y轴上的最大长度ymax为1.5×10-1m。
(二)仿真过程及结果分析：
图3a为3ds max软件下二面角的三角面元模型示意图，图4a为3ds max软件下二面角加遮挡面的三角面元示意图，图4b为3ds max软件下二面角加遮挡面的三角面元原理示意图；图4b的d为入射射线的方向，虚线L是垂直入射方向的平面，即等相位面；图5a为3dsmax软件下城市的三角面元模型示意图；使用本发明方法对图3a、图4a和图5a分别进行成像，结果如图3b、图4c和图5b所示，图3b为使用本发明方法对3ds max软件下二面角的三角面元模型的成像结果示意图；图4c为使用本发明方法对3ds max软件下二面角加遮挡面的模型的成像结果示意图，图5b为使用本发明方法对3ds max软件下城市三角面元模型的成像结果示意图。
从图3b中可以看到许多排清晰的亮点，这些亮点的排布正好能反映图3a中二面角上三角面元的排布，所以，亮点代表的是一个三角面元的回波；另外，在图3b中图不难发现，图3b中左边用红线圈住的1区域能正好和面元的排布一一对应，而图3b中右侧的2区域高于正常亮度的直线区域却在模型中找不到直接与之对应的部分，因此，图3b中右侧的2区域即是由于二次反射造成的。
从图4c中可以看到，图4b中1区域的反射点分布和二面角经过一次反射后的回波分布相同，这是因为图4a中蕴含着二面角，所以形成图4a成像的原理和图3b所述的一致。而图4c中的2区域：入射射线经过图4b原理图所示的4次反射之后回到等相位面位置的距离为图4b中射线所经过的路程，从图4b的几何结构可以看出该距离约为原理图中路程2*l的二倍(距离2*l对应的为距离最远的三角面元的直接照射回波，其在图4c对应于左起第八排的亮点)；而这点在仿真图像图4c中得到了印证；因此认为图4c最右侧出现的亮点即是由于4次反射造成。
从图5b中可以看到图中有多处长短不一的明亮线条结构，这些明亮线条结构即是仿真的二次反射现象，在图5b中还有另一种明显的阴影现象，该阴影的形状和建筑的形状是一致的，这说明图像中阴影的形成是由于建筑物的遮挡造成的，其原理和光学图片中阴影的形成原理一致，这一点在实际的SAR图像中也是可以体现的。
综上所述，仿真实验验证了本发明的正确性，有效性和可靠性。
以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。

Claims (9)

1.一种基于物理光学的太赫兹雷达回波成像方法，其特征在于，包括以下步骤：
步骤1，确定太赫兹雷达，所述太赫兹雷达包含N个天线，每个天线包含若干个阵元；所述太赫兹雷达的检测范围内存在复杂目标，所述复杂目标为电大尺寸目标，然后对复杂目标进行三角面元划分，得到M个三角面元，以及每个三角面元的三个顶点坐标和每个三角面元的法相向量；N、M分别为大于0的自然数；
步骤2，对N个天线分别进行等间隔采样，得到N个天线的NP个采样点，然后对第i个采样点进行随机抖动，得到第i个经过随机抖动的采样点以及第i个经过随机抖动的采样点对应的三角面元序列，所述第i个经过随机抖动的采样点对应的三角面元序列包含ni个三角面元；
然后令i分别取1至NP，进而得到第1个经过随机抖动的采样点至第NP个经过随机抖动的采样点，以及第1个经过随机抖动的采样点对应的三角面元序列至第NP个经过随机抖动的采样点对应的三角面元序列，并分别记为NP个经过随机抖动的采样点以及NP个采样点对应的三角面元序列，i＝1，2，…，NP，
n1+…+ni+…+nNP＝M，P表示每个天线进行等间隔采样后包含的采样点个数，且户个采样点为在每个天线中选取的阵元个数，户为大于0的自然数；
步骤3，预先设定W个慢时间，然后根据预先设定的W个慢时间对M个三角面元分别进行遮挡照亮判断，分别得到M个三角面元经过W个慢时间旋转后的可见三角面元区域，以及W个慢时间对应的被照亮可见三角面元区域的时延集合P；
步骤4，根据W个慢时间对应的被照亮可见三角面元区域的时延集合P，计算得到W个慢时间对应的太赫兹雷达回波数据；
步骤5，对W个慢时间对应的太赫兹雷达回波数据分别沿快时间t进行匹配滤波处理，得到W个慢时间对应的太赫兹雷达回波数据的最终时域输出S(t)，t表示快时间；
步骤6，根据W个慢时间对应的太赫兹雷达回波数据的最终时域输出S(t)，计算W个慢时间复杂目标对应的太赫兹雷达回波成像。
2.如权利要求1所述的一种基于物理光学的太赫兹雷达回波成像方法，其特征在于，在步骤1中，所述复杂目标，还包括：复杂目标位于三维空间坐标系xyz中，且随着第n个慢时间tn旋转，得到第n个慢时间旋转后的角度θn，n＝1，2，…，W，W表示预先设定的慢时间总个数，
θn＝(n-1)×theta，ymax表示复杂目标在y轴上的最大长度，fmin表示设定的成像最小分辨力。
3.如权利要求1所述的一种基于物理光学的太赫兹雷达回波成像方法，其特征在于，在步骤2中，所述第i个经过随机抖动的采样点为：iT+ζ，i＝1，2，…，NP，T表示采样间隔，所述采样间隔为太赫兹雷达发射的电磁波载频波长的一半；ζ表示第i个采样点相邻两次采样间隔之间的均匀分布。
4.如权利要求1所述的一种基于物理光学的太赫兹雷达回波成像方法，其特征在于，步骤3的子步骤为：
3.1令j表示第j个经过随机抖动的采样点，j∈{1，2，…，NP}，NP表示经过随机抖动的采样点总个数；同时j也表示第j个采样点对应的三角面元序列，NP表示NP个采样点对应的三角面元序列总个数；tn表示第n个慢时间，n的初始值为1，n∈{1，2，…，W}，W表示预先设定的慢时间总个数；
在所述NP个采样点分别对应的NP个可见三角面元序列中，第i个采样点对应若干个三角面元，该若干个三角面元包含不止一个的可见三角面元或不包含可见三角面元，所述不止一个的可见三角面元为第i个采样点对应的可见三角面元序列，i＝1，2，…，NP；同样，如果太赫兹雷达发射的电磁波照射到某个三角面元并发生至少一次反射后回到太赫兹雷达时的可见三角面元，为采样点对应的可见三角面元；
如果第j个采样点对应的三角面元序列为空，则第j个采样点没有对应的三角面元；如果第j个采样点对应的三角面元序列不为空，则第j个采样点存在对应的三角面元；设定第j个采样点对应的三角面元序列包含mj个三角面元，令r∈{1，2，…，mj}，r表示第r个三角面元，r的初始值为1；m1+…+mj+…+mNP＝M；
3.2将第j个经过随机抖动的采样点与太赫兹雷达天线相位中心的连线，作为第j个待判断光线，该第j个待判断光线的方向向量为Dj，第j个待判断光线的相位中心为Pj；并将第j个待判断光线的方向向量Dj，作为第j个太赫兹雷达发射的电磁波入射方向单位矢量sj；j∈{1，2，…，NP)，j的初始值为1；
3.3判断太赫兹雷达发射的电磁波是否照射到第r个三角面元经过第n个慢时间旋转后的三角面元不可见面，第r个三角面元经过第n个慢时间旋转后的三角面元的法相向量所在平面，为第r个三角面元经过第n个慢时间旋转后的三角面元可见面；反之为第r个三角面元经过第n个慢时间旋转后的三角面元不可见面；
将第r个三角面元经过第n个慢时间旋转后的三角面元法相向量nr与第j个太赫兹雷达发射的电磁波入射方向单位矢量sj相点乘，如果该点乘结果小于0，则进行3.4；否则进行3.8；
3.4判断第r个三角面元经过第n个慢时间旋转后的三角面元是否位于太赫兹雷达发射的电磁波传播方向的相反方向；
即将第j个待判断光线的相位中心Pj与第r个三角面元经过第n个慢时间旋转后的三角面元中点连接，得到第j个待判断光线、第r个三角面元经过第n个慢时间旋转后的三角面元的矢量mnjr，所述第j个待判断光线、第r个三角面元经过第n个慢时间旋转后的三角面元的矢量mnjr的方向为第j个待判断光线的相位中心Pj指向第r个三角面元经过第n个慢时间旋转后的三角面元中点，然后将第j个待判断光线、第r个三角面元经过第n个慢时间旋转后的三角面元的矢量mnjr与第j个太赫兹雷达发射的电磁波入射方向单位矢量sj相点乘，若此处得到的点乘结果大于0，则进行3.5；否则转至3.8；
3.5判断第j个太赫兹雷达发射的电磁波入射方向单位矢量sj是否偏离第r个三角面元经过第n个慢时间旋转后的三角面元：
计算第r个三角面元经过第n个慢时间旋转后的三角面元中心到第j个太赫兹雷达发射的电磁波入射方向单位矢量sj的距离，即将第j个待判断光线、第r个三角面元经过第n个慢时间旋转后的三角面元的矢量mnjr与第j个太赫兹雷达发射的电磁波入射方向单位矢量sj相叉乘，得到第j个待判断光线、第r个三角面元经过第n个慢时间旋转后的三角面元的向量模unjr；
如果第j个待判断光线、第r个三角面元经过第n个慢时间旋转后的三角面元的向量模unjr小于K倍第r个三角面元经过第n个慢时间旋转后的三角面元最长边长，则第j个待判断光线、第r个三角面元经过第n个慢时间旋转后的三角面元的向量模unjr没有偏离第r个可见光面元经过第n个慢时间旋转后的三角面元，并进行3.6；否则进行3.8；K表示设定的倍数值，为大于0的自然数；
3.6判断第j个太赫兹雷达发射的电磁波入射方向单位矢量sj是否最终照亮第r个三角面元经过第n个慢时间旋转后的三角面元：
即将第j个待判断光线的相位中心为Pj沿着第j个太赫兹雷达发射的电磁波入射方向单位矢量sj的方向在第r个三角面元经过第n个慢时间旋转后的三角面元所在平面进行投影，得到第j个待判断光线与第r个三角面元经过第n个慢时间旋转后的三角面元所在平面的交点Tnjr，然后根据几何算法的方法判断第j个待判断光线与第r个三角面元经过第n个慢时间旋转后的三角面元所在平面的交点Tnjr是否在第r个三角面元经过第n个慢时间旋转后的三角面元内部；
如果第j个待判断光线与第r个三角面元经过第n个慢时间旋转后的三角面元所在平面的交点Tnjr在第r个三角面元经过第n个慢时间旋转后的三角面元内部，则说明第r个三角面元经过第n个慢时间旋转后的三角面元被照亮，并得到第r个三角面元经过第n个慢时间旋转后的三角面元参数，所述第r个三角面元经过第n个慢时间旋转后的三角面元参数为第r个三角面元经过第n个慢时间旋转后的三角面元序号、第r个三角面元经过第n个慢时间旋转后的三角面元各自顶点坐标、第r个三角面元经过第n个慢时间旋转后的三角面元中点坐标、第r个三角面元经过第n个慢时间旋转后的三角面元法相向量，然后执行3.7；
3.7光线反射照亮判断：如果第j个待判断光线照射到第r个三角面元经过第n个慢时间旋转后的三角面元并发生反射，则分别计算第j个待判断光线的优化相位中心P′j和第j个太赫兹雷达发射的电磁波入射方向优化单位矢量s′j；所述第j个待判断光线的优化相位中心P′j为第r个三角面元经过第n个慢时间旋转后的三角面元中点，所述第j个太赫兹雷达发射的电磁波入射方向优化单位矢量s′j为：
s′j＝[(1+sj)·nr]×nr+sj，并将第j个待判断光线的优化相位中心P′j和第j个太赫兹雷达发射的电磁波入射方向优化单位矢量s′j分别替换第j个待判断光线的相位中心Pj和第j个太赫兹雷达发射的电磁波入射方向单位矢量sj，返回子步骤3.3；·表示点乘；
3.8令r加1，返回3.3，直到得到第mj′个可见三角面元经过第n个慢时间旋转后的三角面元的参数，并将此时得到第1个可见三角面元经过第n个慢时间旋转后的三角面元的参数至第mj′个可见三角面元经过第n个慢时间旋转后的三角面元的参数，作为第n个慢时间、第j个采样点对应的可见三角面元序列包含的mj′个被照亮的可见三角面元集合；mj′≤mj；
3.9令j加1，返回3.3，直到得到第n个慢时间、第NP个采样点对应的可见三角面元序列包含的mNP′个被照亮的可见三角面元集合，并将此时得到的第n个慢时间、第1个采样点对应的可见三角面元序列包含的m1′个被照亮的可见三角面元集合至第n个慢时间、第NP个采样点对应的可见三角面元序列包含的mNP′个被照亮的可见三角面元集合，作为第n个慢时间、NP个经过随机抖动的采样点对应的被照亮可见三角面元区域，m′1+…+m′j+…+m′NP≤M，并分别将第n个慢时间对应的被照亮可见三角面元区域记为Pn，将第n个慢时间对应的被照亮可见三角面元区域Pn中第i′个可见三角面元的路程记为则得到第n个慢时间对应的被照亮可见三角面元区域中第i′个可见三角面元的时延tni′，
C表示光速，i′＝1，2，…，mNP′，进而得到第n个慢时间对应的被照亮可见三角面元区域Pn中的mNP′个时延τn；
3.10令n加1，返回3.3，直到得到第W个慢时间、NP个经过随机抖动的采样点对应的被照亮可见三角面元区域，以及第1个慢时间对应的被照亮可见三角面元区域P1中的m1′个时延至第W个慢时间对应的被照亮可见三角面元区域PW中的mNP′个时延，然后将此时得到的第1个慢时间、NP个经过随机抖动的采样点对应的被照亮可见三角面元区域至第W个慢时间、NP个经过随机抖动的采样点对应的被照亮可见三角面元区域，作为M个三角面元经过W个慢时间旋转后的可见三角面元区域；将第1个慢时间对应的被照亮可见三角面元区域P1中的mNP′个时延τ1至第W个慢时间对应的被照亮可见三角面元区域PW中的mNP′个时延τW，作为W个慢时间对应的被照亮可见三角面元区域的时延集合P，P＝{τ1，τ2，…，τW}。
5.如权利要求4所述的一种基于物理光学的太赫兹雷达回波成像方法，其特征在于，所述根据几何算法的方法判断第j个待判断光线与第r个三角面元经过第n个慢时间旋转后的三角面元所在平面的交点Tnjr是否在第r个三角面元经过第n个慢时间旋转后的三角面元内部，过程为：
设定第r个三角面元经过第n个慢时间旋转后的三角面元顶点分别为pr1、pr2和pr3，并根据第j个待判断光线与第r个三角面元经过第n个慢时间旋转后的三角面元所在平面的交点Tnjr，分别计算第一向量F1、第二向量F2和第三向量F3，F1＝pr2-pr1，F2＝Tnjr-pr1，F3＝pr3-pr1，如果F2×F1与F3×F1的方向相反，即(F2×F1)·(F3×F1)＜0，则三个顶点坐标pnr1、pnr2、pnr3分别和Tnjr相交；
如果(F2×F1)·(F3×F1)＞0，则三个顶点坐标pnr1、pnr2、pnr3分别和Tnjr没有相交；第一向量F1和第三向量F3在第二向量F2两侧，(F2×F1)·(F3×F1)＝0，则Tnjr在线段pr1 pr2上；设定第j个待判断光线与第r个三角面元经过第n个慢时间旋转后的三角面元相交于点s，连接点s与第r个三角面元经过第n个慢时间旋转后的三角面元中点，得到线段ls，并设定第r个三角面元经过第n个慢时间旋转后的三角面元三条边分别为L1、L2、L3，如果点s在第r个三角面元经过第n个慢时间旋转后的三角面元内部，则线段ls与L1、L2、L3都不相交；如果点s在第r个三角面元经过第n个慢时间旋转后的三角面元外部，因此根据线段ls与L1、L2、L3的相交性确定点s是否在第r个三角面元经过第n个慢时间旋转后的三角面元内部，进而获知第r个三角面元经过第n个慢时间旋转后的三角面元是否被照亮。
6.如权利要求4所述的一种基于物理光学的太赫兹雷达回波成像方法，其特征在于，在步骤4中，所述得到W个慢时间对应的太赫兹雷达回波数据，其过程为：
4.1初始化：令τdn′表示第n′个慢时间对应的被照亮可见三角面元区域中第d个被照亮可见三角面元的时延，n′∈{1，2，…，W}，W表示被照亮可见三角面元区域的时延集合对应的慢时间个数；d∈{1，2，…，τn′}，τn′表示第n′个慢时间对应的被照亮可见三角面元区域的时延，τn′∈{τ1，τ2，…，τW)；并且确定太赫兹雷达发射线性调频脉冲信号；i′＝1，2，…，mNP′，mNP′表示第NP个采样点对应的可见三角面元序列包含的被照亮可见三角面元的个数；n′的初始值为1；
4.2根据第n′个慢时间对应的被照亮可见三角面元区域的时延，计算得到第n′个慢时间对应的太赫兹雷达回波数据sn′(t)，其表达式为：
其中，τn′i′表示第n′个慢时间、第i′个采样点对应的可见三角面元序列包含的被照亮可见三角面元的时延，rect表示矩形窗函数，t表示时间变量，T表示线性调频脉冲信号的时间宽度，exp表示指数函数，K表示线性调频脉冲信号的调频率，fc表示线性调频脉冲信号载频；
4.3令n′加1，重复子步骤4.2，直到得到第W个慢时间对应的太赫兹雷达回波数据sW(t)，并将此时得到的第1个慢时间对应的太赫兹雷达回波数据s1(t)至第W个慢时间对应的太赫兹雷达回波数据sW(t)，记为W个慢时间对应太赫兹雷达回波数据。
7.如权利要求1所述的一种基于物理光学的太赫兹雷达回波成像方法，其特征在于，在步骤5中，所述得到W个慢时间对应的太赫兹雷达回波数据的最终时域输出S(t)，其过程为：
首先，计算第n′个慢时间、第i′个采样点对应的可见三角面元序列包含的被照亮的可见三角面元的太赫兹雷达回波数据sn′i′(t)，其表达式为：
该第n′个慢时间、第i′个采样点对应的可见三角面元序列包含的被照亮的可见三角面元的太赫兹雷达回波数据的频谱为Sn′i′(f)，其表达式为：
其次，确定W个慢时间对应的太赫兹雷达回波数据的频域匹配滤波器为H(f)，
其中，τn′i′表示第n′个慢时间、第i′个采样点对应的可见三角面元序列包含的被照亮可见三角面元的时延，f表示W个慢时间对应的太赫兹雷达回波数据的频域匹配滤波器H(f)的频率，K表示设定的倍数值，为大于0的自然数；fc表示线性调频脉冲信号载频；T表示采样间隔，所述采样间隔为太赫兹雷达发射的电磁波载频波长的一半；exp表示指数函数，rect表示矩形窗函数；
将W个慢时间对应的太赫兹雷达回波数据的频域匹配滤波器H(f)和第n′个慢时间、第i′个采样点对应的可见三角面元序列包含的被照亮的可见三角面元的太赫兹雷达回波数据的频谱Sn′i′(f)相乘，得到第n′个慢时间、第i′个采样点对应的可见三角面元序列包含的被照亮的可见三角面元的太赫兹雷达回波数据的最终频谱输出Sn′i′out(f)，其表达式为：
然后，对第n′个慢时间、第i′个采样点对应的可见三角面元序列包含的被照亮的可见三角面元的太赫兹雷达回波数据的最终频谱输出Sn′i′out(f)进行快速傅里叶变换，得到第n′个慢时间、第i′个采样点对应的可见三角面元序列包含的被照亮的可见三角面元的太赫兹雷达回波数据的最终时域输出sn′i′out(t)，其表达式为：
sn′i′out(t)＝|K|Tsin c{KT(t-τn′i′)}exp[-j2πfc(t-τn′i′)]
第n′个慢时间、第i′个采样点对应的可见三角面元序列包含的被照亮的可见三角面元的太赫兹雷达回波数据的最终时域输出sn′i′out(t)是在快时间t＝τn′i′处出现的sinc函数；令i′分别取1至mNP′，对第n′个慢时间对应的太赫兹雷达回波数据sn′(t)进行匹配滤波，所述第n′个慢时间对应的太赫兹雷达回波数据sn′(t)包含mNP′个采样点对应的可见三角面元序列，进而得到第n′个慢时间对应的太赫兹雷达回波数据的最终时域输出sn′out(t)；
最后，令n′分别取1至W，分别得到第1个慢时间对应的太赫兹雷达回波数据的最终时域输出slout(t)至第W个慢时间对应的太赫兹雷达回波数据的最终时域输出swout(t)，并记为W个慢时间对应的太赫兹雷达回波数据的最终时域输出S(t)，t表示快时间。
8.如权利要求1所述的一种基于物理光学的太赫兹雷达回波成像方法，其特征在于，在步骤6中，所述W个慢时间复杂目标对应的太赫兹雷达回波成像，其得到过程为：
6.1初始化：令tn″表示第tn″个慢时间，n″∈{1，2，…，W}，W表示被照亮可见三角面元区域的时延集合对应的慢时间个数，n″的初始值为1；第tn″个慢时间时复杂目标旋转的角度为θ(tn″)，θ(tn″)＝theta×tn″，fmin表示设定的成像最小分辨力，ymax表示复杂目标在y轴上的最大长度；
选取复杂目标上任意一个三角面元的中点，记为点目标e，复杂目标为所有可见三角面元的中点形成的点目标的集合；
6.2计算第tn″个慢时间时点目标e到太赫兹雷达的径向距离为Re(tn″)，其表达式为：
Re(tn″)＝R+xecos[θ(tn″)]-yesin[θ(tn″)]
其中，xe表示点目标e在x轴上的位置，ye表示点目标e在y轴上的位置，θ(tn″)表示第tn″个慢时间时复杂目标旋转的角度，θ(tn″)＝θ0+Ωtn″，R表示太赫兹雷达位于x轴负半轴时的位置到复杂目标旋转中心的距离，Ω表示第tn″个慢时间时复杂目标旋转的角速度，tn″表示第tn″个慢时间；
然后计算第tn″个慢时间时点目标e的基带太赫兹雷达回波信号SR(τn″e′，tn″)，其表达式为：
其中，τn″e′表示点目标e在第tn″个慢时间对应的太赫兹雷达回波数据经过匹配滤波后得到的sinc函数峰值位置，p(xe，ye)表示点目标e的反射密度函数，ρ表示太赫兹雷达图像，xe表示点目标e在x轴上的位置，ye表示点目标e在y轴上的位置，表示太赫兹雷达发射电磁波的载频，Re(tn″)表示第tn″个慢时间时点目标e到太赫兹雷达的距离；
6.3对于复杂目标上点目标e的反射密度函数ρ(xe，ye)，其中点目标e在y轴上的位置ye通过确定，并确定点目标e在y轴上的位置ye包括个位置，且第e′个位置ye′通过确定，τn″e′表示第e′个位置ye′在第tn″个慢时间处对应的sinc函数峰值，yimage为设定的太赫兹雷达图像在y轴大小，fmin表示使用者要求的成像最小分辨力；然后计算点目标e在第e′个位置ye′处的反射密度函数ρ(xe，ye′)，分别令进而分别得到点目标e在第1个位置y1处的反射密度函数ρ(xe，y1)至点目标e在第个位置处的反射密度函数记为点目标e最终的太赫兹雷达回波数据成像；
6.4令快时间t＝τn″e′，τn″e′表示第e′个位置ye′在第tn″个慢时间处对应的sinc函数峰值，得到W个慢时间对应的太赫兹雷达回波数据在第e′个位置、第tn″个慢时间处的最终时域输出S(τn″e′，tn″)，其表达式为：
其中， 表示太赫兹雷达发射电磁波的载频，Re(tn″)表示第tn″个慢时间时点目标e到太赫兹雷达的距离，ρ(xe，ye′)表示点目标e在第e′个位置ye′处的反射密度函数，exp表示指数函数；
结合第tn″个慢时间时点目标e到太赫兹雷达的径向距离为Re(tn″)，
Re(tn″)＝R+xecos[θ(tn″)]-yesin[θ(tn″)]，得到：
对第tn″个慢时间，由于且该式中的值非常小，因此令k sinθ(tn″)＝0，k cosθ(tn″)＝k×tn″；
其中， 表示太赫兹雷达发射电磁波的载频；进而计算得到第tn″个慢时间时复杂目标的优化基带太赫兹雷达回波信号其表达式为：
其中，表示点目标e在第e′个位置ye′处的优化反射密度函数，其表达式为：
其中，xe表示点目标e在x轴上的位置，IFT表示逆傅里叶变换；
6.5令tn″分别取t1至tW，依次重复子步骤6.2至子步骤6.4，直到得到第t1个慢时间时复杂目标的优化基带太赫兹雷达回波信号至第tW个慢时间时复杂目标的优化基带太赫兹雷达回波信号并对此时得到的第t1个慢时间时复杂目标的优化基带太赫兹雷达回波信号至第tW个慢时间时复杂目标的优化基带太赫兹雷达回波信号进行成像处理，得到W个慢时间复杂目标对应的太赫兹雷达回波成像。
9.如权利要求8所述的一种基于物理光学的太赫兹雷达回波成像方法，其特征在于，所述W个慢时间复杂目标对应的太赫兹雷达回波成像，还包括：所述W个慢时间复杂目标对应的太赫兹雷达回波成像为复杂目标在距离维和方位维的分布，且所述W个慢时间复杂目标对应的太赫兹雷达回波成像为二维矩阵，组成该二维矩阵的每一列元素分别是复杂目标在距离维的分布，组成该二维矩阵的每一行元素分别是复杂目标在方位维的分布，而复杂目标在方位维的分布垂直于距离维方向的分布。
CN201611079232.6A 2016-11-30 2016-11-30 基于物理光学的太赫兹雷达回波成像方法 Active CN106597438B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201611079232.6A CN106597438B (zh) 2016-11-30 2016-11-30 基于物理光学的太赫兹雷达回波成像方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201611079232.6A CN106597438B (zh) 2016-11-30 2016-11-30 基于物理光学的太赫兹雷达回波成像方法

Publications (2)

Publication Number Publication Date
CN106597438A CN106597438A (zh) 2017-04-26
CN106597438B true CN106597438B (zh) 2019-05-21

Family

ID=58595734

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201611079232.6A Active CN106597438B (zh) 2016-11-30 2016-11-30 基于物理光学的太赫兹雷达回波成像方法

Country Status (1)

Country Link
CN (1) CN106597438B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111090102B (zh) * 2018-10-08 2021-10-15 中国科学院沈阳自动化研究所 一种超分辨反射式太赫兹三维目标重建成像方法
CN109655841B (zh) * 2018-12-13 2020-08-14 北京遥测技术研究所 太赫兹多模实时成像系统

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104698456A (zh) * 2015-03-26 2015-06-10 西安电子科技大学 基于预调制aic的太赫兹雷达成像方法
CN105842698A (zh) * 2016-05-18 2016-08-10 西安电子科技大学 一种高分辨sar图像的快速仿真方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2012000075A1 (en) * 2010-06-28 2012-01-05 Institut National D'optique Synthetic aperture imaging interferometer

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104698456A (zh) * 2015-03-26 2015-06-10 西安电子科技大学 基于预调制aic的太赫兹雷达成像方法
CN105842698A (zh) * 2016-05-18 2016-08-10 西安电子科技大学 一种高分辨sar图像的快速仿真方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
基于快速物理光学法的太赫兹目标RCS计算;江月松 等;《光学学报》;20141231;第34卷(第12期);第1-7页

Also Published As

Publication number Publication date
CN106597438A (zh) 2017-04-26

Similar Documents

Publication Publication Date Title
Stroe et al. Discovery of spectral curvature in the shock downstream region: CIZA J2242. 8+ 5301
CN102812379B (zh) 超分辨率成像雷达
CN104007440A (zh) 一种加速分解后向投影聚束合成孔径雷达成像方法
KR20110097277A (ko) 합성 개구면 레이더 영상 형성 시뮬레이션 시스템 및 방법
CN112612024A (zh) 微波阵列快速成像方法
CN106597438B (zh) 基于物理光学的太赫兹雷达回波成像方法
Deming et al. Polar format algorithm for SAR imaging with Matlab
CN105044720A (zh) 一种基于直角坐标系的后向投影成像方法
Jones et al. Marine vehicles simulated SAR imagery datasets generation
Zherdev et al. Object recognition using real and modelled SAR images
Cochin et al. MOCEM-An'all in one'tool to simulate SAR image
Miller et al. Impact of ground mover motion and windowing on stationary and moving shadows in synthetic aperture radar imagery
Kelly et al. A fast decimation-in-image back-projection algorithm for SAR
Demirci et al. Wide‐field circular SAR imaging: An empirical assessment of layover effects
Anglberger et al. Fast ISAR image generation through localization of persistent scattering centers
CN109061645A (zh) 综合孔径辐射计虚拟相关稀疏成像方法及系统
CN112612026B (zh) 基于双雷达距离像融合的目标角分辨方法
Dumont et al. DIONISOS. An end-to-end SAR simulator
Perälä et al. ESPRESS—On efficient bistatic characterization of radar targets
CN112764030A (zh) 基于电磁与统计混合模型的海面sar图像杂波生成方法
CN114624666B (zh) 基于时域回波散射模型的雷达图像快速获取方法和装置
Mackay et al. Millimetre wave ray tracing simulator with phase and beam effects using the Wigner distribution function
Zhang et al. The video SAR simulation algorithm based on the DEM and SBR
Haapaniemi Generating synthetic SAR images to train target classification algorithms
Yu et al. A Rapid Single-view Radar Imaging Method with Window Functions

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant