本发明属于雷达信号处理领域，特别涉及一种基于位置与波形参数联合估计的多发射源被动定位方法。
背景技术：
被动定位技术由于雷达自身不需要发射源发射电磁波，相比于主动雷达，其具有低截获率、低功耗等优点。目前，正处在研制当中的无源雷达有两种：一种是依靠目标辐射器探测和定位目标，而另一种则利用其它辐射器发射的电磁波来进一步定位目标。这里考虑前一种情况，其中目标本身就是辐射源，它们可能是机载雷达、干扰机和通信设备等，因此也可称之为发射机定位。
现有技术中，有一种联合处理各基站接收机观测数据的直接定位方法(dpd)，该方法没有传统定位方法的参数提取过程，尽可能地保留了目标信息。仿真显示其较传统两步定位方法定位精度大大提高。这种方法涉及了两种情况，一种是目标发射机发射的信号已知，如训练信号或同步信号，我们称这种情况的dpd定位算法为dpd-known算法。另一种较为常见的情况是目标发射的信号完全未知，主要针对非合作的发射机。由于不知道目标信号的形式，dpd定位方法通过使目标函数特征值最大的方法实现了发射机定位，仿真显示其定位性能仍然优于传统的定位方法，称这种方法为dpd-unknown算法。但是，该方法忽略了发射信号的信号特征，其定位精度受限，不能适应低信噪比下的目标的高精度定位。
此外，还有一种将dpd定位技术与短时傅里叶变换相结合，在估计信号参数的同时定位发射源。但该方法考虑的是单个发射源的定位问题，无法解决实际场景更为常见且更具有挑战性的多目标的定位问题。
技术实现要素：
为了解决现有定位技术定位精度差且无法定位多个发射机目标的问题，本发明提供了一种基于位置与波形参数联合估计的多发射源被动定位方法，能实现低信噪比情况下的多发射源的高精度定位。
本发明采用的方案为：
一种基于位置与波形参数联合估计的多发射源被动定位方法，包括以下步骤：
步骤1，初始化参数，其中，所述参数包括接收机个数、各接收机位置、采样间隔、各接收机通道上的噪声协方差矩阵、短时傅里叶变换的窗长度和步长以及霍夫变换参数；
步骤2，读取各接收机的量测，对量测进行采样，得到数个离散的量测向量；
步骤3，根据所述短时傅里叶变换的窗长度和步长，对所述量测向量进行短时傅里叶变换，得到各接收机量测的功率谱；
步骤4，对所述功率谱进行霍夫变换，得到发射源个数估计值、接收信号中每个线性调频信号成分在每个频率上的最大功率，估计出对应发射源发射信号的最大功率谱和对应的瞬时频率；
步骤5，应用基本分类算法，以所述发射源个数估计值为类别数，将所述瞬时频率分为所述类别数类；
步骤6，划分目标位置网格搜索区间、信号发射时间的网格搜索区间及信号长度的网格搜索区间；
步骤7，采用网格搜索的方法对发射源进行定位。
进一步地，所述步骤4中的所述霍夫变换应用hough和houghlines函数。
进一步地，所述步骤5中的所述基本分类算法为k-means分类算法。
进一步地，所述步骤7包括以下流程：
步骤71，初始化已定位发射源个数；
步骤72，采用网格搜索的方法定位一个发射源；
步骤73，更新所述已定位发射源个数；
步骤74，判断所述已定位发射源个数与所述发射源个数估计值的大小；
步骤75，当所述已定位发射源个数不小于所述发射源个数估计值时，结束定位。
进一步地，所述步骤7还包括以下流程：
步骤76，当所述已定位发射源个数小于所述发射源个数估计值时，流程回到所述步骤72。
进一步地，所述步骤72包括以下流程：
步骤721，设置代价初值；
步骤722，选择一个网格点作为发射信号的位置坐标、发射时间及发射信号长度，计算网格点的信号起始时间和终止时间，生成对应的时间区间；
步骤723，根据生成的所述时间区间，从所述瞬时频率中选择对应所述时间区间的频率值作为线性调频信号的瞬时频率估计值，从所述最大功率谱中选择对应所述时间区间的最大功率谱值；
步骤724，从数个所述量测向量中估计出定位的发射源的相同数量个线性调频信号的瞬时频率的估计值的均值；
步骤725，采用最小二乘法，通过参数拟合从所述均值中估计出定位的发射源发射的线性调频信号的初始频率和调频斜率；
步骤726，根据网格点和估计的线性调频信号参数，计算定位的发射源的代价函数值，其中，线性调频信号参数包括线性调频信号的初始频率和调频斜率；
步骤727，当所述代价函数值大于所述代价初值时，以所述代价函数值替代所述代价初值，将对应的网格点的位置坐标作为发射源的位置估计；
步骤728，判断数据平面上所有的网格点是否已经被遍历；
步骤729，当所有的网格未被全部遍历，流程回到步骤722；
步骤730，当所有的网格点被遍历，最终的发射源的位置估计为该发射源的定位结果。
进一步地，所述步骤726包括以下流程：
根据发射信号的起始时间、信号长度及估计的线性调频信号参数构建匹配函数，其中，线性调频信号参数包括线性调频信号的初始频率和调频斜率；
将匹配函数在时间上离散化，得到离散后的匹配函数；
计算对应网格点的代价函数值。
本发明的有益效果：针对未知数目目标提出了可以有效定位所有目标的算法，利用霍夫变换处理短时傅里叶变换得到的功率谱，结合基本分类算法将多目标的定位问题分解为多个单目标的定位问题。同时针对单个目标的定位问题，通过充分挖掘发射信号的波形信息，利用线性调频信号的时频特征，结合短时傅里叶变换联合估计信号参数和位置，使定位性能在信号参数的帮助下有极大提高。
附图说明
图1为本发明的流程图。
图2为图1中步骤7的流程图。
图3图2中步骤72的流程图。
图4为图3中步骤726的流程图。
具体实施方式
本发明主要采用计算机仿真的方法进行验证，所有步骤、结论都在matlab-r2016a上验证正确。具体实施步骤如下：
下面结合附图对本发明的实施例做进一步的说明。
请参阅图1，图1为本发明的流程图。
本发明的一种基于位置与波形参数联合估计的多发射源被动定位方法，具体实施步骤如下：
步骤1，初始化参数，其中，所述参数包括接收机个数、各接收机位置、采样间隔、各接收机通道上的噪声协方差矩阵、短时傅里叶变换的窗长度和步长以及霍夫变换参数。
本实施例中，初始化的参数包括：目标发射机的位置、雷达基站接收机的个数l＝4、雷达基站接收机位置、采样周期ts＝10-7s，短时傅里叶变换的窗长度40ts，步长为1ts，噪声协方差矩阵
步骤2，读取各接收机的量测，对量测进行采样，得到数个离散的量测向量。
本实施例中，得到的l个用于仿真的离散的量测向量rl，rl＝[rl[0],rl[1],...,rl[ns-1]]，l＝1,2,...,l。其中，ns表示量测样本个数。
这里用于仿真的三个调频信号的初始频率分别为20mhz、16mhz和18mhz，调频斜率分别为0.1mhz/μs、0.15mhz/μs和-0.2mhz/μs。
步骤3，根据所述短时傅里叶变换的窗长度和步长，对所述量测向量进行短时傅里叶变换，得到功率谱。
本实施例中，根据设置的短时傅里叶变换的窗长度和步长，对量测向量rl作短时傅里叶变换，得到其功率谱
步骤4，对所述功率谱进行霍夫变换，得到发射源个数估计值、接收信号中每个线性调频信号成分在每个频率上的最大功率，估计出对应发射源发射信号的最大功率谱和对应的瞬时频率。
本实施例中，对短时傅里叶变换得到的功率谱进行霍夫变换，得到目标个数的估计值，即发射源个数估计值q，得到接收信号中每个线性调频信号成分在每个频率上的最大功率，即发射信号的最大功率谱及对应的频率，估计出对应目标发射信号的瞬时频率
本实施例中，霍夫变换的应用主要需要调用hough及houghlines两个函数，分别用于处理功率谱以及提取最大功率谱对应的瞬时频率。
步骤5，应用基本分类算法，以所述发射源个数估计值为类别数，将所述瞬时频率分为所述类别数类。
本实施例中，应用基本分类算法，将从接收信号中分离得到的各发射信号的瞬时频率
本实施例中，应用的基本分类算法为k-means分类算法。
步骤6，划分目标位置网格搜索区间、信号发射时间的网格搜索区间及信号长度的网格搜索区间。
本实施例中，划分目标位置网格搜索区间(xgrid,ygrid)，信号发射时间的网格搜索区间tgrid及信号长度的网格搜索区间tgrid。
步骤7，采用网格搜索的方法对发射源进行定位。
请参阅图2，步骤7通过以下流程实现：
步骤71，初始化已定位发射源个数。
本实施例中，初始化已定位的发射源个数q＝0。
步骤72，采用网格搜索的方法定位一个发射源。
本实施例中，采用网格搜索的方法定位q个发射源中的一个。
步骤73，更新所述已定位发射源个数。
本实施例中，更新已定位的发射源个数q＝q+1，即每一次对一个发射源的定位完成，更新一次已定位发射源的个数。
步骤74，判断所述已定位发射源个数与所述发射源个数估计值的大小。
本实施例中，判断q和q的大小关系。
步骤75，当所述已定位发射源个数不小于所述发射源个数估计值时，结束定位。
本实施例中，当q≥q时，即已定位发射源个数不小于发射源个数估计值时，对发射源的定位完成。
步骤76，当所述已定位发射源个数小于所述发射源个数估计值时，流程回到所述步骤72。
本实施例中，当q＜q时，即已定位发射源个数不足发射源个数估计值，回到步骤72继续对其它发射源进行定位。
请参阅图3，步骤72通过以下流程实现：
步骤721，设置代价初值。
本实施例中，设置代价初值value＝0。
步骤722，选择一个网格点作为发射信号的位置坐标、发射时间及发射信号长度，计算网格点的信号起始时间和终止时间，生成对应的时间区间。
本实施例中，选取一个网格点η1grid＝(xgrid,ygrid,tgrid,tgrid)作为发射信号的位置坐标、发射时间及发射信号长度，根据网格点的参数，计算信号的起始时间tls,q(η1grid)和终止时间tle，q(η1grid)，生成对应的时间区间[tls,q(η1grid),tle,q(η1grid)]。
步骤723，根据生成的所述时间区间，从所述瞬时频率中选择对应所述时间区间的频率值作为线性调频信号的瞬时频率估计值，从所述最大功率谱中选择对应所述时间区间的最大功率谱值。
本实施例中，根据生成的时间区间，从瞬时频率
步骤724，从数个所述量测向量中估计出定位的发射源的相同数量个线性调频信号的瞬时频率的估计值的均值。
本实施例中，对从l个量测中估计出定位的发射源的l个线性调频信号的瞬时频率的估计值
步骤725，采用最小二乘法，通过参数拟合从所述均值中估计出定位的发射源发射的线性调频信号的初始频率和调频斜率。
本实施例中，采用最小二乘法，通过参数拟合从均值
步骤726，根据网格点和估计的线性调频信号参数，计算定位的发射源的代价函数值，其中，线性调频信号参数包括线性调频信号的初始频率和调频斜率。
本实施例中，根据网格点η1grid＝(xgrid,ygrid,tgrid,tgrid)和估计的线性调频信号参数
步骤727，当所述代价函数值大于所述代价初值时，以所述代价函数值替代所述代价初值，将对应的网格点的位置坐标作为发射源的位置估计。
本实施例中，比较一个网格点得出的所定位的发射源的代价函数值
步骤728，判断数据平面上所有的网格点是否已经被遍历。
本实施例中，直到全部网格点被遍历，从而估计的发射源位置为
步骤729，当所有的网格未被全部遍历，流程回到步骤722。
本实施例中，若未被全部遍历，回到步骤722继续对网格点进行遍历。
步骤730，当所有的网格点被遍历，最终的发射源的位置估计为该发射源的定位结果。
本实施例中，当所有的网格点被遍历，此时该发射源的定位结果为能使代价初值最大的那个代价函数值所对应的网格点的位置坐标，保存定位结果
请参阅图4，步骤726通过以下流程实现：
步骤61，根据发射信号的起始时间、信号长度及估计的线性调频信号参数构建匹配函数，其中，线性调频信号参数包括线性调频信号的初始频率和调频斜率。
本实施例中，根据信号的起始时间tls,q(η1grid)、信号长度tgrid及估计的线性调频信号参数
式中，l＝1,2,...,l，
步骤62，将匹配函数在时间上离散化，得到离散后的匹配函数。
本实施例中，得到离散后的匹配函数
步骤63，计算对应网格点的代价函数值。
本实施例中，计算对应网格点η1grid的代价函数值：
式中，h表示对矩阵共轭装置，rl代表第l个接收机通道上目标附近噪声及干扰的协方差矩阵。
本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的权利要求范围之内。