辐射源航迹矢量直接估计方法

本发明属于辐射源定向定位领域，特别涉及一种辐射源航迹矢量直接估计方法。

背景技术：

辐射源的航迹估计是指：多个观测站同步获取辐射源的观测数据，基于某种算法，确定辐射源的运动轨迹。传统方法中，将辐射源估计问题分解为三个步骤：第一步获得单次观测参数，例如tdoa(到达时间差)，toa(到达时间)，aoa(到达方向角)，rssi(接收信号强度指示)等；第二步基于单次观测参数，建立优化模型，估计辐射源位置；第三步基于滤波方法估计目标轨迹。该方法将原本的优化问题人为地拆解为3个优化子问题，从最优化的角度来说，优化过程中没有充分利用三个子问题之间的约束条件，容易受到优化门限效应的限制，损失了优化性能。直接航迹矢量优化方法中将决定航迹的参数矢量建模为高维空间的一个点，并基于接收波形和信号传播模型构建航迹矢量的最大似然估计模型，运用最优化算法，实现从观测波形直接到航迹矢量估计的一步优化，充分利用了多个观测站观测数据之间，以及同一个观测站的不同时刻的观测数据之间的自然冗余关系，可以获得比传统三步航迹估计方法更优的性能，尤其是在辐射源运动规律明确，低信噪比条件下，性能提升显著。但是直接航迹矢量估计方法需要更大的计算量，在计算能力受限的应用场景中难以发挥作用。

技术实现要素：

为此，本发明提供一种辐射源航迹矢量直接估计方法，能够在目标运动规律明显、信噪比低、采集快拍数少的条件下，获得较高的航迹估计。

按照本发明所提供的设计方案，一种辐射源航迹矢量直接估计方法，用于水声通信信号识别，包含如下内容：

基于多个观测站采集的接收波形和电磁波传播模型，构建相邻阶段的目标运动模型；

依据目标运动模型和相邻两个阶段观测目标状态转移概率构建包含接收波形数据和航迹矢量的似然函数；

通过对似然函数进行优化求解获取目标航迹矢量的最大似然估计。

作为本发明辐射源航迹矢量直接估计方法，进一步地，运动辐射源航迹估计的状态矢量表示为其中，为第k次观测的状态矢量，(xk,yk)为第k次观测时目标位置坐标，vk,φk分别为第k次观测时目标运动的速度大小和速度角度，为角加速度，k表示等间距观测次数；假设目标辐射源运动过程中状态转移是一阶马尔科夫过程，相邻阶段的目标运动模型表示为xk+1＝φ(xk)+γ(wk)，φ(·)表示状态转移关系，wk表示状态转移扰动项，γ(·)表示扰动转移关系。

作为本发明辐射源航迹矢量直接估计方法，进一步地，对于复杂运动目标，采用通过不同运动模型之间的切换来逐段估计的混合模型作为目标运动模型；对于非马尔科夫过程的运动目标，采用通过历史数据进行深度学习的lstm模型作为运动模型。

作为本发明辐射源航迹矢量直接估计方法，进一步地，依据目标运动模型获取相邻两个阶段的状态转移，建立航迹状态矢量与接收波形之间的关系。

作为本发明辐射源航迹矢量直接估计方法，进一步地，观测站i在第k次观测时获得的波形向量表示为：ri,k＝h(ai,k,fi,k,ti,k,fk,tk,sk)+nk，其中，ai,k表示第k次观测时从运动辐射源到观测站i的路径衰减复系数，fi,k表示第k次观测时从运动辐射源到观测站i的频率漂移，ti,k表示第k次观测时从运动辐射源到接收站的传播延迟，tk表示第k次观测运动辐射源与观测站的钟差，fk表示第k次观测中运动辐射源的频率漂移，sk表示发送波形向量，nk表示第k次观测的观测噪声向量。

作为本发明辐射源航迹矢量直接估计方法，进一步地，针对接收波形，对观测函数进行建模，表示为：ri,k＝d(fi,k)f^-1t(ti,k)t(tk)fd(fk)skai,k+nk，其中，d(f)表示频率漂移f赫兹矩阵，t(t)表示信号延迟t秒矩阵，f表示傅里叶变换矩阵。

作为本发明辐射源航迹矢量直接估计方法，进一步地，针对目标运动模型，依据航迹状态矢量和发送波形，并根据目标运动速度计算接收机的多普勒频偏；依据多普勒频偏将发送信号在时域移动相应频率，根据钟差和航迹状态矢量计算接收站延迟；将信号傅里叶变换到频域，并在频域内移动与延迟对应的频率，将信号通过傅里叶反变换转换到时域；将信号乘以路径衰减系数，获取航迹状态矢量与接收波形的关系。

作为本发明辐射源航迹矢量直接估计方法，进一步地，假设各观测站观测噪声为高斯噪声，依据第k次观测时观测站i接收波形的概率密度函数，在给定辐射源第k次观测状态xk的条件下，获取第k+1次观测时转移到xk+1的概率密度函数；联合两个概率密度函数，得到基于观测站i的航迹状态矢量估计的对数似然函数。

作为本发明辐射源航迹矢量直接估计方法，进一步地，运用双向矢量估计方法对航迹状态矢量估计的似然函数进行求解，获取用于获取目标航迹矢量的最大似然估计值。

作为本发明辐射源航迹矢量直接估计方法，进一步地，似然函数求解包含如下内容：首先，基于第k次观测数据，使用直接定位方法估计目标位置，并基于crlb作为参数估计计算直接定位方法的定位误差协方差矩阵；然后，基于从0到k-1次直接定位结果采用扩展kalman滤波方法获得航迹矢量的kalman滤波估计，计算kalman滤波估计的误差协方差矩阵，其中，k表示等间距观测次数；再者，基于kalman滤波估计的结果，逆向计算航迹矢量。

本发明的有益效果：

本发明通过构建包含波形数据和航迹矢量的似然函数，直接获得航迹矢量的最大似然估计，充分利用同一目标的不同观测站的观测数据之间、以及不同时刻的观测数据之间的自然冗余关系，能够在目标运动规律明显、信噪比低、采集快拍数少的条件下，获得较高的航迹估计，显著提高航迹估计性能，具有较好的应用前景。

附图说明：

图1为实施例中辐射源航迹矢量直接估计流程示意；

图2为实施例中辐射源的航迹估计原理示意；

图3为实施例中辐射源航迹估计算法流程示意；

图4为实施例中构建波形与航迹矢量关系流程示意；

图5为实施例中前向efk滤波与逆向矢量平滑估计流程示意；

图6为实施例中目标航迹与估计结果示意；

图7为实施例中目标航迹与估计结果局部示意；

图8为实施例中目标航迹估计误差示意。

具体实施方式：

为使本发明的目的、技术方案和优点更加清楚、明白，下面结合附图和技术方案对本发明作进一步详细的说明。

辐射源的航迹估计是指：多个观测站同步获取辐射源的观测数据，基于某种算法，参见图2所示，确定辐射源的运动轨迹。本发明实施例，提供一种辐射源航迹矢量直接估计方法，参见图1所示，包含如下内容：

s101、基于多个观测站采集的接收波形和电磁波传播模型，构建相邻阶段的目标运动模型；

s102、依据目标运动模型和相邻两个阶段观测目标状态转移概率构建包含接收波形数据和航迹矢量的似然函数；

s103、通过对似然函数进行优化求解获取目标航迹矢量的最大似然估计。

通过构建相邻两个阶段的目标运动模型，建立起状态向量与观测波形之间的关系；最后基于状态转移概率波形和观测概率模型构建航迹矢量的似然函数，并通过求解获取航迹矢量的最大似然估计值，能够在目标运动规律明显、信噪比低、采集快拍数少的条件下获得较高的航迹估计。

假设只有1个未知航迹的运动辐射源，辐射源发出已知波形的信号。运动辐射源在运动过程中，等间隔观测k次，每次观测持续时间为td，每次观测获得n个快拍,第k次观测中发送波形记为其中tn＝n/fs,fs是观测站的采样间隔。本发明实施例，进一步地，运动辐射源航迹估计的状态矢量表示为其中，为第k次观测的状态矢量，(xk,yk)为第k次观测时目标位置坐标，vk,φk分别为第k次观测时目标运动的速度大小和速度角度，为角加速度，k表示等间距观测次数；参见图3所示，假设辐射源运动时遵循一定的规律，即下一时刻的运动状态与之前的运动状态有关，目标辐射源运动过程中状态转移是一阶马尔科夫过程，相邻阶段的目标运动模型表示为xk+1＝φ(xk)+γ(wk)，φ(·)表示状态转移关系，wk表示状态转移扰动项，γ(·)表示扰动转移关系。扰动转移关系表示系统扰动与状态扰动的关系。假设系统扰动项wk是一个高斯随机过程，单次扰动的协方差矩阵为典型的状态转移模型有cv(匀速直线运动模型)，ca(匀加速运动模型)，ctrv(匀转弯和匀速运动模型)等。进一步地，对于复杂运动目标，例如受人工控制的运动目标，可以采用混合模型，逐段估计运动模型，并进行不同运动模型的切换。对于非马尔可夫过程的运动目标，可以采用深度学习lstm模型等，运用历史数据对运动模型进行深度学习。

作为本发明实施例中辐射源航迹矢量直接估计方法，进一步地，依据目标运动模型获取相邻两个阶段的状态转移，建立航迹状态矢量与接收波形之间的关系。

假设有m个已实现同步采集数据的静止观测站，第i个观测站的位置已知并记为pi，观测站i在第k次观测时，按照fs的采样速率进行采样，获得波形向量为

ri,k＝h(ai,k,fi,k,ti,k,fk,tk,sk)+nk

其中ai,k是第k次观测时，从运动信源到观测站i的路径衰减复系数，是一个标量，表征链路衰减和收发系统相位差；fi,k是第k次观测，从辐射源到观测站i的频率漂移，产生原因是目标与观测站的径向运动导致的多普勒频移，ti,k是从辐射源到接收站的传播延迟，tk是第k次观测辐射源与观测站的钟差；fk是第k次观测中辐射源的频率漂移，假设观测站的时钟和频率通过gnss驯服时钟授时，能够获得足够高的时钟同步经度和频率经度。sk是根据发送波形向量。nk是第k次观测的观测噪声向量，其协方差矩阵为

一个典型的观测函数可以建模为：

ri,k＝d(fi,k)f^-1t(ti,k)t(tk)fd(fk)skai,k+nk

其中d(f)是信号源频率漂移f赫兹的矩阵

其中t⁽ⁱ⁾＝i/fs。

t(t)是信号延迟t秒的矩阵

其中f⁽ⁱ⁾＝fs(i+1-n/2).

f是傅里叶变换矩阵

作为本发明实施例中辐射源航迹矢量直接估计方法，进一步地，针对目标运动模型，依据航迹状态矢量和发送波形，并根据目标运动速度计算接收机的多普勒频偏；依据多普勒频偏将发送信号在时域移动相应频率，根据钟差和航迹状态矢量计算接收站延迟；将信号傅里叶变换到频域，并在频域内移动与延迟对应的频率，将信号通过傅里叶反变换转换到时域；将信号乘以路径衰减系数，获取航迹状态矢量与接收波形的关系。

对发送时域波形矢量sk在时域移动频率fk赫兹，表示辐射源的实际频率与标称频率的漂移量fk赫兹，fk是未知量，需要在航迹估计过程中进行估计。通过傅里叶变化矩阵f获得移频后的频域矢量，并在频域移动时间tk秒，表示辐射源与接收系统的钟差为tk秒，tk是未知量，需要在航迹估计过程中进行估计。在频率域移动时间ti,k秒，表示第k次观测从辐射源到接收站i的延迟ti,k秒，ti,k与第k次观测时辐射源到观测站i的距离有关，即与航迹矢量x有关。通过傅里叶逆变换f^-1，将频域矢量变换到时域矢量，并在时域将信号频率移动频率fi,k赫兹，表示因辐射源和接收站的径向速度导致的多普勒频移为fi,k赫兹，fi,k与第k次观测时辐射源与观测站i的径向速度有关，即与航迹矢量x有关。将时域波形矢量乘以路径衰减系数标量ai,k，表示因路径传播衰减和收发系统载波相位差导致的复衰减系数ai,k，ai,k是一个未知量，需要在航迹矢量估计过程中估计。参见图4所示，依据上述步骤可以构建观测矢量和状态矢量x的关系，可以根据实际情况进行删减和添加。

作为本发明辐射源航迹矢量直接估计方法，进一步地，假设各观测站观测噪声为高斯噪声，依据第k次观测时观测站i接收波形的概率密度函数，在给定辐射源第k次观测状态xk的条件下，获取第k+1次观测时转移到xk+1的概率密度函数；联合两个概率密度函数，得到基于观测站i的航迹状态矢量估计的对数似然函数。

假设各站观测噪声为高斯噪声，则第k次观测，观测站i的接收波形矢量ri,k的概率密度函数为：

其中

在给定辐射源第k次观测状态xk的条件下，第k+1次观测时转移到xk+1的概率密度函数为：

联合上述两个似然函数，并取对数，去掉常数项，得到基于观测站i的航迹估计的对数似然函数是：

上述似然函数是x，ai,k,tk,fk的函数。

作为本发明实施例中辐射源航迹矢量直接估计方法，进一步地，运用双向矢量估计方法对航迹状态矢量估计的似然函数进行求解，获取用于获取目标航迹矢量的最大似然估计值。进一步地，求解过程可设计为包含如下内容：首先，基于第k次观测数据，使用直接定位方法估计目标位置，并基于crlb作为参数估计计算直接定位方法的定位误差协方差矩阵；然后，基于从0到k-1次直接定位结果采用扩展kalman滤波方法获得航迹矢量的kalman滤波估计，计算kalman滤波估计的误差协方差矩阵，其中，k表示等间距观测次数；再者，基于kalman滤波估计的结果，逆向计算航迹矢量。

依据最大化似然函数获得航迹状态矢量x的最优估计上述优化问题的决策变量包括x,ai,k,tk,fk，决策矢量的维度为：5k+2mk+k+k，直接求解该优化问题复杂度很大。本案实施例中可利用优化方法快速求解算法。参见图5所示，具体优化算法可设计为如下：

步骤1、基于第k次观测数据，使用直接定位方法估计目标位置，即求解以下优化问题：

记并基于crlb分析方法计算直接定位方法的定位误差协方差矩阵rk。

步骤2、基于从0到k-1次直接定位结果zk采用扩展kalman滤波(ekf)的方法，获得航迹矢量的kalman滤波估计计算kalman估计的误差协方差矩阵pk。

步骤3、基于kalman滤波的结果，逆向计算航迹矢量：

其中

上述优化算法可以快速求解高维矢量优化问题，但不限于上述算法，rts平滑方法、动态规划方法、分布式搜索算法(例如遗传算法、粒子群算法)等都可以实现高维矢量优化的目的。

为验证本案方案有效性，下面结合试验数据做进一步解释说明：

假设有4个观测站，位于边长为100km的四个角，目标从[0,50]km水平向东飞行，速度为720km/s，每隔2秒观测一次，一共观测1000秒。目标运动模型为ctrv。基于接收波形，分别使用dpd方法、ekf方法和本案提出的dtd方法对航迹进行估计，结果如图6～8所示。通过仿真数据可以看出，dpd方法的平均定位误差为1239米，ekf方法的平均误差为423米，本案dtd方法平均误差为241米。相对与直接定位和ekf，本案能够充分利用同一目标的不同观测站的观测数据之间、以及不同时刻的观测数据之间的自然冗余关系，航迹估计性能和定位精度获得了较显著的提升。

除非另外具体说明，否则在这些实施例中阐述的部件和步骤的相对步骤、数字表达式和数值并不限制本发明的范围。

在这里示出和描述的所有示例中，任何具体值应被解释为仅仅是示例性的，而不是作为限制，因此，示例性实施例的其他示例可以具有不同的值。

应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步定义和解释。

最后应说明的是：以上所述实施例，仅为本发明的具体实施方式，用以说明本发明的技术方案，而非对其限制，本发明的保护范围并不局限于此，尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，其依然可以对前述实施例所记载的技术方案进行修改或可轻易想到变化，或者对其中部分技术特征进行等同替换；而这些修改、变化或者替换，并不使相应技术方案的本质脱离本发明实施例技术方案的精神和范围，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应所述以权利要求的保护范围为准。