低采样点数下的太赫兹高分辨率成像方法与流程

本发明属于成像
技术领域：
，涉及一种太赫兹成像方法，具体涉及一种低采样点数下的太赫兹高分辨率成像方法。
背景技术：
：太赫兹波是一种频率范围在0.1～10THz的电磁波，介于微波与红外波段之间。与微波、红外线等类似，太赫兹波也可以应用在成像系统中。并且，太赫兹波具有与其他电磁波不相同的特性，例如能量低、穿透力强、具有高方向性等，这些特性使得太赫兹成像系统具有其他成像系统不具备的优势，因而具有良好的应用前景。但是，现有的太赫兹成像技术均受到分辨率的制约。为了提高太赫兹成像的分辨率，目前较为普遍的方法是提高雷达合成孔径的宽度，使采样得到的太赫兹信号源增多。上述方法是通过改进太赫兹成像系统的硬件条件来提高分辨率的。然而，合成孔径的宽度越高，太赫兹成像系统的制造成本也就越高，不利于太赫兹成像技术的推广应用。并且，现有技术中针对成像系统的成像方法进行的改进较为少见。基于逆傅里叶变换的成像法是目前常用的一种成像方法。这种方法较为简单，通过直接对方向位坐标进行傅里叶逆变换即可得到目标各个点的坐标，从而得到目标图像。这种方法的步骤较为简单，在样本数量大的情况下对目标点坐标的估计精度比较高，但是在采样数少的情况下估计精度则较低。并且，这种算法需要将目标区域进行网格化划分，然后只能估计出落在网格上的目标点坐标；对于网格外的坐标点，该方法采取最接近的网格点坐标作为目标点的估计结果，因此这种方法进行成像的估计误差较大，导致图像分辨率难以提升。技术实现要素：为解决上述问题，提出一种改进的太赫兹成像方法，使低采样条件下也能够对目标点的坐标作出准确的估计从而得到高分辨率的图像，本发明采用了如下技术方案：本发明提供一种太赫兹高分辨率成像方法，用于根据低采样点数条件下得到的太赫兹回波来获取目标点坐标，从而实现太赫兹高分辨率成像，其特征在于，包括如下步骤：步骤A，使用太赫兹成像系统向待成像物体发射太赫兹波，获取目标点的太赫兹回波，得到Q次采样的回波矩阵Sr；步骤B，根据回波矩阵，用Akaike信息量准则估计出信号源数量；步骤C，对步骤A得到的回波矩阵Sr进行x维及y维的傅里叶变换，得到以波数域分量表示的回波矩阵Sr(kx’(m)，ky’(n))，其中，kx’(m)对应观测空间域x维的波数域分量，ky’(n)对应观测空间域y维的波数域分量；步骤D，对步骤C中得到的回波矩阵Sr(kx’(m)，ky’(n))进行相位补偿，得到相位补偿后的回波矩阵步骤E，通过增广矩阵束算法估算得到目标点的一系列坐标值，并进行匹配得到目标点坐标。进一步地，本发明提供的太赫兹高分辨率成像方法，还可以具有如下技术特征：其中，步骤E包括如下子步骤：子步骤E1，将步骤D中得到的相位补偿后的回波矩阵所形成的矩阵记为H，对于H的每一行向量进行矩阵扩展，得到矩阵Hm；子步骤E2，对子步骤E1得到的矩阵Hm进行矩阵平移，得到矩阵Hm的增广矩阵He；子步骤E3，对子步骤E2得到的增广矩阵He进行特征值分解，提取与x维和y维相关的两维目标信息相关的特征向量；子步骤E4，根据子步骤E3得到的特征向量进行广义特征值计算，得到信号的x维频率信息以及y维频率信息；子步骤E5，根据子步骤E4得到的x维频率信息以及y维频率信息，计算得到目标点的x维坐标以及y维坐标；子步骤E6，对子步骤E5得到的x维坐标以及y维坐标进行配对，得到目标点坐标。发明作用与效果根据本发明的太赫兹高分辨率成像方法，由于采用了MEMP算法，该方法在低采样点数的条件下仍然具有较高的目标点估计精度，并且不需要对目标区域进行网格化划分，在任意位置的目标点都能够较好地估计得出其坐标，因此能够在低采样点数的条件下实现高分辨率成像。附图说明图1是本发明的太赫兹高分辨率成像方法的流程图；图2是本发明的太赫兹成像系统的观测空间示意图；图3是本发明的太赫兹高分辨率成像方法与基于逆傅里叶变换的太赫兹成像结果对比图。具体实施方式以下结合附图及实施例来说明本发明的具体实施方式。<实施例>图1是本发明的太赫兹高分辨率成像方法的流程图。如图1所示，本发明的太赫兹高分辨率成像方法主要包括步骤A～步骤E这五个步骤，具体过程如下所述。步骤A，使用太赫兹成像系统向待成像物体发射太赫兹波，获取目标点的太赫兹回波，得到Q次采样的回波矩阵Sr，主要包括如下子步骤：子步骤A1，使用太赫兹成像系统向待成像物体发射太赫兹波，并从待成像物体表面返回的回波，进行一次采样过程。图2是本发明的太赫兹成像系统的观测空间示意图。如图2所示，观测空间为一个分布有多个采样点的平面，发出太赫兹波后接收相应的回波。例如，图中的太赫兹成像系统从采样平面方位发出太赫兹波，到达待成像物体(即、成像目标)的一个目标点(x，y)以后，产生一个回波返回采样点(x’,y’,Z0)。由此，采样点就得到了关于该目标点的回波信息。在本实施例中，收集的回波为散射较强的目标点产生的回波，这些散射较强的点在下文中称为散射点。另外，由于本实施例中，收集的回波为散射点所产生的回波，因此散射点的数量也就等于信号源的数量。将上述回波信息按下式(1)的形式进行记录，sr(xm′,yn′,tq)=Σl=1Lσl(xl,yl)ej2πf(tq-2(xm′-xl)2+(yn′-yl)2+Z02c)---(1)]]>式(1)中，x'm为观测空间收发机沿x维的空间采样位置，m＝0,…,M-1，M为满足均匀采样时在x维的空间采样点数量；y'n为沿y维的空间采样位置，n＝0,…,N-1，N为满足均匀采样时在y维的空间采样点数量。σl(xl，yl)表示第l个散射点的散射系数，L为散射点总个数(即、信号源总个数)，f是信号频率，(xl，yl)为散射点坐标，Z0是目标中心到探测器观测平面的距离，c是光速。子步骤A2，将上述过程中的采样时刻记为tq，收集tq时刻全部采样位置的的回波，并将其整理成如式(2)所示的向量形式，得到tq时刻回波：sr(tq)=sr(x0′,y0′,tq)...sr(x0′,yN-1′,tq)...sr(xM-1′,y0′,tq)...sr(xM-1′,yN-1′,tq)---(2)]]>由式(2)可知，由于满足均匀采样时，x维的空间采样点数量为M，y维的空间采样点数量为N，因此tq时刻的回波中总共包含M*N(以下记为MN)个采样点的采样信息。子步骤A3，用太赫兹成像系统进行Q次采样，每次采样均按照子步骤A1～子步骤A2所描述的过程收集每次采样时刻下的回波，从而得到Q次采样的回波。将上述Q次采样的回波整理成如式(3)所示的矩阵形式：式(3)中，矩阵Sr为回波矩阵，即Q次采样下，MN个采样点所得到的回波所构成的矩阵。步骤B，根据上述步骤A得到的回波矩阵Sr，用Akaike信息量准则估计出信号源数目，具体包括如下子步骤：子步骤B1，根据式(4)获取回波的协方差矩阵R，R=Sr×SrH---(4)]]>式(4)中，为Sr的共轭转置矩阵。子步骤B2，基于Akaike信息量准则(以下简称AIC准则)，根据式(5)估算出信号源数量L0：AIC(L)＝2T(MN-L)lnΛ(L)+2L(2MN-L)(5)式(5)中，T为采样数，Λ(L)为似然函数，具体计算式如式(6)所示：Λ(L)=1/MN-LΣi=L+1MNλi/(Πi=L+1MNλi)1/MN-L---(6)]]>式(6)中，λi为R的特征值。根据式(5)及式(6)，得到0,1,...,MN-1之间使AIC准则取值最小的L值，该L值即为信号源数量L0。步骤C，对步骤A得到的回波矩阵Sr进行x维及y维的傅里叶变换，得到以波数域分量表示的回波矩阵Sr(kx’(m)，ky’(n))，其中，kx’(m)对应观测空间域x维的波数域分量，ky’(n)对应观测空间域y维的波数域分量，具体包括如下子步骤：子步骤C1，对发射信号和接收信号作混频处理，消去时间变量tq，过程如式(7)所示：sr(xm′,yn′)=sr(xm′,yn′,tq)×e-j2πftq=Σl=1Lσl(xl,yl)ej2πf(tq-2(xm′-xl)2+(yn′-yl)2+Z02c)×e-j2πftq=Σl=1Lσl(xl,yl)e-j4πf2(xm′-xl)2+(yn′-yl)2+Z02c---(7)]]>式(7)中，e为自然常数，j为复信号的表示。子步骤C2，将式(7)中的球面波经过平面波叠加展开，得到波数域的回波表达式，即：sr(xm′,yn′)=Σl=1Lσl(xl,yl)e-j2k(xm′-xl)2+(yn′-yl)2+Z02=Σl=1Lσl(xl,yl){Σkx′(m)Σky′(n)e-jkx′(m)(xl-xm′)-jky′(n)(yl-yn′)+jkz′(m,n)Z0}---(8)]]>其中，k＝2πf/c是波数，kx'(m)和ky'(n)分别对应观测空间域x和y维的波数域分量，为观测空间域z维的波数域分量。由于从波数域回波的表达式中可以看出目标散射系数σl(xl,yl)与回波的傅里叶变换关系，因此利用傅里叶变换关系，式(8)可以简化为式(9)的形式：sr(kx(m),ky′(n))=ΔFT2D{sr(xm′,yn′)}e-jkz·(m,n)Z0=FT2D{σl(xl,yl)}---(9)]]>因此，通过对回波sr(xm',yn')的x、y方向维回波进行如上所述的傅里叶变换，可以得到波数域回波与散射点坐标的直接关系。步骤D，对步骤C中得到的回波矩阵Sr(kx’(m)，ky’(n))进行相位补偿，得到相位补偿后的回波矩阵具体过程如下：由式(9)可以得到，第(m,n)个观测点的观测回波为如下式(10)所述的形式：sr(kx′(m),ky′(n))=Σl=1Lσl(xl,yl)e-jkx′(m)xi-jky′(n)yle-jkz′(m,n)Z0---(10)]]>式(10)中，得到的回波信息中包含了和目标信息无关的相位信息为了消去该相位信息，需要根据式(11)对回波进行相位补偿：s~r(kx′(m),ky′(n))=Δsr(kx′(m),ky′(n))ejkz′(m,n)Z0=Σl=1Lσl(xl,yl)e-jkx′(m)xl-jky′(n)yl---(11)]]>由此，得到的相位补偿后的回波矩阵中不再包含相位信息。步骤E，通过增广矩阵束算法估算得到目标点的一系列坐标值，并进行匹配得到目标点坐标，包括如下子步骤：子步骤E1，将步骤D中得到的相位补偿后的回波矩阵所形成的矩阵记为H，对于H的每一行向量进行矩阵扩展，得到矩阵Hm，过程如下：令对H的每一行向量进行矩阵扩展，得到的矩阵Hm如下式(12)所示：式(12)中，Hm是I×(N-I+1)的Hankel矩阵，Hm的每一列长度为I，并且满足I≥L。子步骤E2，对子步骤E1得到的矩阵Hm进行矩阵平移，得到矩阵Hm的增广矩阵He，过程如下式(13)所示：式(13)中，He是K×(M-K+1)的Hankel分块矩阵，He的每一列有K个Hankel块矩阵，并且满足K≥L。子步骤E3，对子步骤E2得到的增广矩阵He进行特征值分解，提取与x维和y维相关的两维目标信息相关的特征向量，过程如下：首先，对矩阵He进行如式(14)所示的特征值分解，He=UsΣsVsH+UnΣnVnH---(14)]]>式(14)中，Us是包含了和目标信息相关的特征向量，Un是和目标信息无关的特征向量。然后，提取构建和两维信号信息相关的特征向量：令Us1＝Us(1:KL-L,:),Us2＝Us(L+1:end,:)，并且令Usp＝P·Us，提取和x维、y维目标信息相关的特征向量，即Usp1＝Usp(1:KL-L,:)，Usp2＝Usp(L+1:end,:)。其中，P的表达式如下式(15)所示:P=pT(1)pT(1+I)...pT(1+(K-1)I)...pT(I)pT(I+I)...pT(I+(K-1)I)---(15)]]>式(15)中，p(I+(K-1)I)表示第I+(K-1)I行元素为1，其余元素为0的KI×1矩阵；pT为p的转置矩阵。子步骤E4，根据子步骤E3得到的特征向量进行广义特征值计算，得到信号的x维频率信息以及y维频率信息，过程如下：计算Us1-λUs2的广义特征值，得到信号的x维频率信息fx,l；计算Usp1-λUsp2的广义特征值，得到信号的y维频率信息fy,l。子步骤E5，根据子步骤E4得到的x维频率信息以及y维频率信息，计算得到目标点的x维坐标以及y维坐标，过程如下：根据下式(16)得到目标点的x维坐标，根据下式(17)得到目标点的y维坐标：fx,l＝Δfxxl(16)fy,l＝Δfyyl(17)式(16)中，为x方向波束域采样间隔；式(17)中，为y方向波束域采样间隔。子步骤E6，由于子步骤E5中，得到的坐标点并不一定就是目标的坐标点，因此还需将x维坐标信息和y维坐标信息进行配对，得到实际上的目标点的坐标，该配对过程包括如下子步骤：子步骤E6.1，令i＝1，计算找到Js(i,j)取最大值时的j并组成配对(i,j(i))。子步骤E6.2，令i＝i+1，计算Js(i,j)，j＝1,2,...,L，找到Js(i,j)取最大值时的j并组成配对(i,j(i))，但是j≠j(k),其中k＝1,2,...,i-1。子步骤E6.3，重复子步骤E6.2，直到i＝L-1。子步骤E6.4，根据(i,j(i))，完成(yi,zj(i))的配对。经过上述过程，即可得到目标点的坐标。图3是本发明的太赫兹高分辨率成像方法与基于逆傅里叶变换的太赫兹成像结果对比图。图3(a)和图(b)中，白色原点指示真实的散射点位置。图(a)示出在采样点数为50×50，x、y维目标点坐标不在网格上时，通过基于逆傅里叶变换法得出的太赫兹成像结果图；图(b)示出在采样点数为50×50，x、y维目标点坐标不在网格上时通过本发明的成像方法得出的太赫兹成像结果图。由图3(a)和图3(b)的对比可知，基于逆傅里叶变换的成像方法的估计误差为0.0513，经过这种算法仿真出来的一些目标点连在一起无法区分开。而本发明的成像方法的估计误差为0.0224，并且没有出现目标点连在一起的情况。可以看出，相比于基于逆傅里叶变换的成像方法，本发明的成像方法在较低的样本数下仍然具有较高的样本估计精度；并且，本发明的成像方法不需要对目标区域进行网格划分，因此对于没有分布在网格上的目标点也可以较好的估计出来，其结果就是在进行太赫兹成像时有更高的分辨率。实施例作用与效果根据本实施例的太赫兹高分辨率成像方法，由于采用了MEMP算法，该方法在低采样点数的条件下仍然具有较高的目标点估计精度，并且不需要对目标区域进行网格化划分，在任意位置的目标点都能够较好地估计得出其坐标，因此能够在低采样点数的条件下实现高分辨率成像。当前第1页1 2 3