本发明涉及一种将水听器镶嵌在电缆上形成线列阵，由拖曳电缆拖在舰尾后水中探测目标的声纳。具体涉及一种声纳拖曳线列阵进行被动合成孔径时存在子阵位移误差情况下，同时对子阵位阵列移误差和波达方向估计方法。
背景技术：
拖曳声纳(towed sonar)是将换能器基阵拖曳在运载平台尾后水中探侧目标的声纳。拖线阵声纳基阵是柔性的。拖曳过程中由于舰艇机动和海流影响，以及自身抖动共振，其阵形难以保持稳定，阵型畸变将使得拖线阵声纳难以达到理论性能，这个问题尤其在采用自适应信号处理，空间谱估计等现代阵处理方法中更为严重。因为常规阵处理方法只对阵元能量进行累加，而现代阵处理方法还对阵元信号相关性，信号协方差矩阵特征值等进行计算，虽然现代阵处理方法可以大大提高精度和目标分辨能力，但对阵型稳定更为严格。此外，拖线阵没有垂直孔径，不能区分水面水下目标，在航运密集海域或编队作战时大量目标对其性能产生严重干扰，会产生到处是目标的“满天星”现象。低频主动拖曳线列阵声呐是目前探测安静型潜艇最有效的手段。由于拖曳线列阵声呐基阵远离本舰,且为柔性阵,探测的方位存在较大误差。现代拖曳线列阵声纳系统逐渐趋向于工作在低频段，探测距离不断提高，向着更高的作用距离和探测精度发展。在这样环境下，若要获得较好的空间分辨能力，需要使用孔径较大的阵列，这在实际中，通常意味着更大的系统复杂度和更高的设备成本。在不改变拖曳阵列参数的情况下，利用拖曳线列阵进行被动合成孔径(PSA)，增加了有效孔径的可能性。PSA技术是利用信号时间和空间上的相关性，以及阵列的运动信息、位置信息，构造一个含有多个虚拟子阵的阵列来获得更高的方位分辨力。目前，常见的运动直线阵PSA方法有： Yen、Carey被动合成孔径方法、基于快速傅立叶变换的被动合成孔径，以及扩展拖曳阵测量方法等。虽然具有多个虚拟子阵的PSA已经被广泛地应用到了声纳波束形成，以及波达方位(DOA)估计等领域，但在目标数目较多，以及目标舰船与本舰运动参数与真实条件不一致的情况下，性能会急剧下降。具体地，只存在单目标时，PSA方法可以在相对速度出现误差情况下工作，直接从数据中得到真实的相位校正因子，合成阵元位置不重叠时，仍可准确地进行DOA估计；但在多目标时，若相对速度无误差，PSA方法可以工作，若相对速度存在误差，则此时的DOA结果，甚至会差于单阵常规波束形成的DOA结果，阵列扩展不但没有带来好处，反而增加了定位误差。除此之外，还有其它因素影响DOA的估计精度,例如阵元间距(通常取工作波长的一半)、目标来波之间的角度差、信号的相干性等等。在PSA的基础上，通常利用高分辨DOA方法可以得到目标的精确角度估计。许多高分辨方法，例如多信号分类 MUSIC(Multiple Signal Classification)多重信号分类方法、旋转不变子空间方法、稀疏空间谱估计等可以在阵列参数精确已知的情况，提高阵列的空间分辨率。MUSIC算法是一种基于矩阵特征空间分解的方法，从几何角度讲，从几何角度，信号处理的观测空间可以分解为两个空间是正交的信号子空间和噪声子空间。MUSIC算法利用这两个互补空间之间的正交特性来估计空间信号的方位。噪声子空间的所有向量被用来构造谱，所有空间方位谱中的峰值位置对应信号的来波方位，其基本思想是对任意阵列输出数据的协方差矩阵进行特征分解，从而得到与信号分类相对应的信号子空间和与信号分量相正交的噪声。MUSIC算法处理任务是设法估计出入射到阵列的空间信号的个数D以及空间信号源的强度及其来波方向。在实际处理中，观测到的阵列输出数据复向量Y得到的数据是有限时间段内的有限次数的样本，也称快拍或快摄，但是原型MUSIC算法要求来波信号是不相干的。以MUSIC为代表的算法存在一个缺点，即对相干信号处理的不理想。在针对相干信号源的一系列处理方案中，比较经典的是空间平滑技术，如空间平滑(SS)和修正的空间平滑(MSS)算法。然而，空间平滑技术是以损失阵列有效孔径为代价的，而且只适用于等距均匀线阵(ULA)。当存在阵元误差时，假定的杂波模型与实际接收数据不匹配，在模型参数失配的情况下，原本属于信号子空间的目标信号会被错误划分到噪声子空间内，破坏MUSIC算法目标函数的正交性，导致定位失效和。同理，其余各类高分辨DOA算法也对阵列流形的误差很敏感，且在低信噪比、少快拍数、声源相干情况出现时，高分辨算法的性能将会大为退化。由于重构杂波协方差矩阵依赖于杂波模型的准确性，若在假定的杂波模型与实际的数据不匹配时，将导致算法性能下降甚至失效。因此在拖曳线列阵PSA应用中，阵元误差，信号相干、阵元互耦、通道失配等均会使理想模型不再成立。如何更有效地利用接收数据来估计实际的真实阵列流行向量，并进行准确方位估计，成为一个十分重要的课题。
目前，国内外在解决这类问题中已进行了一系列的探索。这些文献根据各自应用领域中存在的实际误差来源，从不同侧面对阵列误差进行校正。例如Weiss和Friedlander提出了一种基于最大似然估计阵列自校正的方法，该方法对目标方位和位移失配量进行分块交替优化迭代。由于此类方法是对所有阵元进行建模，在应用到长基线、大规模、长时间PSA中时，会导致维数急剧增多，目标函数收敛极为缓慢、计算效率急剧下降。因此，需要发展针对长基线、大规模、长时间情况拖曳线列阵PSA的阵列流行自校正和方位估计方法。
拖曳线列阵声呐也称“拖曳阵声呐”(可简称拖线阵)。它是将水听器镶嵌在电缆上形成线列阵，由拖曳电缆拖在舰艇尾后水中探测目标的声呐。主要用于听测潜艇辐射噪声，进行远程监视、测向和识别，有的也可用于测距。由线列基阵、拖曳电缆、收放装置和绞盘、电子机柜等组成。拖曳线列阵又由前导段、仪器段、基阵段、后导段和尾段构成，阵长数十米至数百米，工作深度可变。具有基阵尺寸大、工作频率低、利于线谱检测，能远距离隐蔽地发现目标等优点；但对拖带舰艇的旋回和倒车等机动有不利影响。声源、海洋信道与水听器阵是水声学研究中的三个基本要素。声源在水中辐射声信号，是形成声场的源头；海洋信道则决定了声波在海洋中的传播特点；水听器阵用以接收声信号，对水中的声场分布进行采样。这三者之间相互密切联系，构成不可分割的统一整体。已知其中两者，就可以推断第三者，这就是水声匹配场处理的基本依据。如果已知水听器阵接收信号和海洋信道信息，待求解的是包括声源位置在内的声源信息，这就是匹配场被动定位。如果已知的是水听器阵接收信号和包括声源位置在内的声源信息，待求解的是海洋信道信息，这就是匹配场反演(MFI： MatchedField Inversion)。它们都是水声匹配场处理研究的重要内容。
匹配场处理技术近年来在水下目标检测、被动定位、海洋环境参数反演等方面的应用受到广泛的关注。水声阵列信号的匹配场匹配场处理(MFP：Matched Field Processing)是利用海洋环境参数和声传播信道特性，通过水下声场模型计算得到接收基阵的声场幅度和相位，形成拷贝场向量，并与基阵接收数据进行匹配，从而实现水下目标的被动定位和海洋环境参数的精确估计。基于远近场声传播特性的拖线阵声纳平台，远场目标信号的平面波传播特性将匹配场定位技术和平面波目标方位估计技术结合，在远场平面波假设下，有N个角频率为ω的声源信号入射到具有P个子阵的线阵上。假设每个子阵内部有Mp个阵元，整个阵列的阵元数为子阵内的阵元位置精确已知，第p个子阵的模型可以表示为
xp(t)＝Aps(t)+ep(t), (1)
其中，第p个子阵阵列流形矩阵Ap＝[ap(θ1),ap(θ2),...,ap(θN)]，θn是第n个来波的方位；向量ap(θn)＝[1,exp(-jωdcos(θn)/c),...,exp(-jω(M-1)dcos(θn)/c)]T是第p个子阵θn方向对应的阵列流形向量，c为声速，上标T表示转置；向量s(t)＝[s1(t),s2(t),...,sN(t)]T表示t时刻对应的信号波形向量；向量ep(t)＝[e1(t),e2(t),...,eM(t)]T代表t时刻第p个子阵对应的噪声。通常情况下，声源数量N小于阵元个数M，来波在空域具有稀疏性。假设标量rp是第p个子阵的第一个阵元到第一个子阵的第一个阵元的距离，全阵的阵列流形向量可以表示为
aw(θ)＝V(θ)h(θ), (2)
式中为子阵内阵列流行矩阵，ap(θ)为第p个子阵的阵列流行向量，h为子阵间近似阵列流行向量，θ为来波的方位角。尽管ap(θ)精确已知，实际测量的子阵相对位置向量一般不等于预先设置的子阵相对位置向量r，所以真实的全阵阵列流形向量与预置阵列流形向量aw(θ)间会存在偏差。存在子阵间位移失配的阵列模型和一种目标测向方法。
技术实现要素：
本发明针对现有技术存在的不足之处，提供一种能够减小定位误差，并能得到准确的方位估计值，更高角度分辨率，可以同时提高稳健性的声纳拖曳线列阵子阵误差失配模型的目标定位的方法。
本发明的上述目的可以通过以下技术方案予以实现，一种拖曳线列阵声纳子阵误差失配估计方法，其特征在于包括如下步骤：在存在子阵间位移失配的阵列流形矩阵模型中，将真实的全阵列流形矩阵表示为每个子阵的位移误差量和每个子阵误差量对全阵列流形矩阵贡献的线性组合；将子阵间的位置误差引入测向模型，对含有子阵间位移失配的全阵列模型运用贝叶斯法则建立数据融合模型，对子阵位置误差向量β和目标真实方位角α-1进行同时求解，根据融合传感器节点采集到的数据，使用贝叶斯算法同时对子阵位移误差和波达方向进行估计，计算得到多快拍观测值的似然函数；用后验函数得到均方根误差随信噪比变化的方位估计值和位移误差估计值。
本发明相比于现有技术具有如下有益效果：
本发明将存在子阵位移误差情况下真实的全阵阵列矩阵近似为每个子阵的位移误差量和每个子阵误差量对全阵阵列流形矩阵的贡献的线性组合，然后使用贝叶斯框架进行求解，可以在低信噪比情况下对子阵位置误差向量和目标真实方位进行同时求解，当存在子阵位移误差情况下可同时实现阵列的自校正与目标准确测向，通过将子阵间的位置误差引入测向模型，通过贝叶斯算法实现对子阵位移误差和波达方向的同时估计，对阵列流形的误差不敏感，在减小了定位误差的同时提高了稳健性。
本发明使用贝叶斯算法计算得到如图2所示的位移误差估计值的均方根误差随信噪比的变化，可以看到位移误差估计量的均方根误差随信噪比的增加而减小。方位估计值的均方根误差随如图3所示的信噪比的变化，比使用多重信号分类MUSIC算法得到的均方根误差稳定1度附近。
本发明使用贝叶斯算法计算得到的均方根误差均随信噪比的增加而减小，且均小于1度，性能较MUSIC算法有较大提高。信噪比为0dB时使用全阵CBF算法、全阵MUSIC算法以及本发明方法的定位结果如图4所示，通过比较可以看出，本发明不仅能够得到准确的方位估计值，还具有比其他方法更高的角度分辨率。
在使用具有子阵位移误差的阵列对相干声源的方位角进行估计中取得了明显实施效果。与直接使用全阵CBF方法和全阵MUSIC方法进行方位估计相比，其优势主要在于：
(1)通过将子阵误差引入阵列模型，可以同时对子阵位移误差和目标波达方向进行估计；
(2)通过将含有子阵位移误差的阵列模型引入贝叶斯框架，可以得到更稳健、精度更高的定位性能。
本发明适用于传感器子阵处理、射频天线阵列扩展、声纳被动合成孔径等场景，主要用于听测潜艇辐射噪声，进行远程监视、测距、测向和识别阵列信号处理、声纳信号处理方法，
附图说明
图1是本发明具有子阵位移误差的阵列扩展示意图。
图2是本发明位移误差估计值的均方根误差随信噪比的变化曲线示意图。
图3是本发明方位估计值的均方根误差随信噪比的变化曲线示意图。
图4是本发明信噪比为0dB时的定位结果曲线示意图。
下面结合附图对本发明作进一步的描述。
具体实施方式
参阅图1-图4。根据本发明，在存在子阵间位移失配的阵列流形矩阵模型中，将真实的全阵列流形矩阵表示为每个子阵的位移误差量和每个子阵误差量对全阵列流形矩阵贡献的线性组合；将子阵间的位置误差引入测向模型，对含有子阵间位移失配的全阵列模型，运用贝叶斯法建立数据融合模型和融合传感器节点采集到的数据来计算目标的位置和速度，对子阵位置误差向量β和目标真实方位角α-1进行同时求解，计算得到多快拍观测值的似然函数；对子阵位移误差和波达方向进行同时估计，用后验函数得到均方根误差随信噪比变化的方位估计值和位移误差估计值。具体步骤为
步骤1：失配后的真实子阵间阵列流行向量在真实子阵相对位置向量的预设子阵相对位置向量r处进行一阶泰勒展开逼近，得到子阵间近似阵列流行向量
式中，向量h为子阵间近似阵列流行向量，e为欧拉常数，j为虚数单位常数，k为波数，上标T表示转置，P为子阵数，θ为表示来波的方位角，向量为真实的子阵相对位置向量，为第p个子阵的真实相对位置，向量r＝[r1,...,rP]T为预设的子阵相对位置向量，rp为第p个子阵的预设相对位置，β＝[β1,...,βP]T为子阵位置误差向量，为第 p个子阵位移误差，diag(β)表示以向量β的元素作为对角线元素的对角矩阵。
根据一阶泰勒展开逼近得到的子阵间近似阵列流行向量，将真实的全阵阵列流形向量 aw近似为
子阵间阵列流行矩阵
式中，是来波方位为θ时第p列为vp(θ)其余列为零向量的矩阵，向量vp为子阵间阵列流行矩阵V(θ)的第p列，ap(θ)为第p个子阵的阵列流行向量。
在全阵列模型中，将全阵的阵列流形向量沿N个θ1-θN的来波方位进行扫描，组合成全阵的近似阵列流行矩阵A1，再将每个子阵的位移误差量和每个子阵误差量对全阵列流形矩阵的贡献的线性组合为:
且无子阵间位置误差时预设的全阵阵列流形矩阵Aw＝[aw(θ1),...,aw(θN)]，第p个子阵误差投影矩阵Bp＝-jk[Vp(θ1)h(θ1,r),...,Vp(θN)h(θN,r)]，
式中，N为扫描方位网格数量，aw为预设的全阵阵列流形向量，βpBp是第p个子阵位置误差引起的全阵阵列流形矩阵误差的一阶逼近乘积。
步骤2.在基于子阵间位移失配阵列模型的贝叶斯定位中，使用贝叶斯算法对全阵列模型子阵间含子阵位置误差向量β和目标真实方位α-1进行同时求解，得到表示近似的全阵阵列流行矩阵和近似的全阵的阵列模型全阵接收信号向量x(t)＝Φ(β)s(t)+e(t)，根据全阵接收信号向量x(t)＝[x1(t),x2(t),...,xN(t)]Tt时刻阵列接收到的信号、t时刻对应目标信号波形的声源向量s(t)＝[s1(t),s2(t),...,sN(t)]T和第t时刻全阵噪声的噪声向量e(t)＝[e1(t),e2(t),...,eM(t)]T，在多快拍情况下，将全阵的阵列模型改写为多个快拍数下的阵列接收信号矩阵X＝Φ(β,θ)S+E，其中，T是快拍数，矩阵X＝[x(1),x(2),...,x(T)]， x(t)为阵列t时刻的阵列接收信号向量，矩阵S＝[s(1),s(2),...,s(T)]表示声源信号矩阵，s(t)为 t时刻的声源信号向量，E＝[e(1),e(2),...,e(T)]表示噪声矩阵，e(t)表示t时刻的噪声向量，θ表示扫描来波的方位角向量。
来波的扫描方位角向量为L个方位扫描网格向量，且L＞＞N，t时刻所有扫描方向的信号向量为矩阵是L个扫描网格条件下的T个时刻的目标信号矩阵，按照多快拍数下的阵列接收信号矩阵公式 X＝Φ(β,θ)S+E，将存在子阵间位移失配的阵列流形矩阵的L个方位扫描网格条件下的全阵列模型可以表示为
式中，矩阵为存在子阵间误差全阵扫描阵列流行矩阵，向量表示扫描来波的扫描方位角向量，为第l个扫描方位，为无子阵间位置误差全阵扫描阵列流形矩阵，向量为第l个扫描方位时的真实的全阵扫描阵列流形向量，矩阵为第p个子阵误差的扫描投影矩阵，矩阵为扫描方位为θl时第p列为其余列为零向量的矩阵，β为子阵位置误差向量，是L个扫描网格条件下的T个时刻的目标信号矩阵，为t时刻所有扫描方向的信号向量，E表示噪声矩阵。
在每个阵元的噪声独立且满足均值为0方差为目标真实初始方位角的复高斯分布时，求得多快拍观测值的似然函数
式中，I是L维单位矩阵，det()表示矩阵求行列式，exp()表示求指数函数，||()||2表示求向量二范数。噪声精度α0服从参数为a和b的伽马分布
p(α0|a,b)＝Gamma(α0|a,b)
式中，函数Γ(a)表示变量为a的伽马函数。
信源向量服从复高斯分布，目标信号矩阵的概率密度函数为
式中，
协方差矩阵
α2为第l个方位的精度。信号精度向量参数α＝[α1,α2,...,αL]服从参数为c和d的伽马分布
目标信号矩阵的后验概率分布为
式中，α0为初始噪声精度参数，μ(t)表示后验均值向量，Σ表示的后验协方差矩阵。
根据贝叶斯模型更新方法求得t时刻初始噪声精度参数α0的更新后验均值向量μ(t)和后验协方差矩阵Σ，其中
后验均值向量
后验协方差矩阵
式中，为子阵间误差时全阵扫描阵列流行矩阵，H为共轭转置，β为子阵位置误差向量，为t时刻未知的接收信号向量，矩阵Λ为先验协方差矩阵，矩阵Λ-1为矩阵Λ的逆。
根据贝叶斯模型更新方法求得t时刻信号精度参数αi和噪声精度参数α0的更新噪声精度参数更新噪声精度参数计算公式
更新方差σ的
式中，σ为方差，||()F||代表向量的Frobenius范数算子，均值矩阵H＝[μ(1),μ(2),...,μ(T)]，Σii是协方差矩阵Σ的第i个对角元素，标量γi＝1-αiΣii。
根据贝叶斯模型迭代更新方法求得子阵位置误差向量β的估计值
第一中间矩阵T＝G+Q，
第二中间矩阵
第三中间矩阵
第四中间矩阵
第一中间向量
第二中间向量
整个贝叶斯模型更新迭代过程可以总结如下：叶斯模型迭代更新方法对初始噪声精度参数α0、信号精度向量α以及子阵位置误差向量β赋予初值，使用更新后验均值向量μ(t)公式(8)和后验协方差矩阵Σ公式(9)更新均值向量μ和协方差矩阵Σ，然后使用按照公式更新噪声精度参数计算公式(10)、更新方差σ的(σ2)new计算公式(11)和子阵位置误差向量β的估计值计算公式(12)更新初始噪声精度参数α0，信号精度向量α以及子阵位置误差向量β，重复以上过程，直到收敛。迭代完成后的初始噪声精度参数α0、信号精度向量α和子阵位置误差向量β分别表示了噪声能量、特定方位的信号能量和每个虚拟子阵的位移误差。
下面以具体范例说明：
参阅图1。子阵位移误差的阵列扩展的预设位置包含虚拟子阵1、虚拟子阵2、虚拟子阵3，虚拟子阵2的实际位置偏离β＝0.11，虚拟子阵3的实际位置偏离β＝0.2，有目标数K＝2个远场窄带信号入射到阵元数M＝4元均匀线阵上，中心频率f＝250Hz，声源信号入射角度θ1和θ2分别为60°和65°，线阵阵元间距为0.68米。P＝3个虚拟子阵首阵元的预设位置分别距离初始位置0m、3.4m和6.8m，子阵位移误差β1,β2和β3分别为0米、0.11m和0.2m。每个位置阵列采集的快拍数T＝200，获得目标信号矩阵X。根据公式(5)计算第p个子阵位置误差引起的全阵阵列流形矩阵矩阵Bp。
初始噪声精度参数α0的超参数分别设为a＝b＝1×10-4，信号精度向量的超参数分别设为c＝1，d＝0.01；迭代过程初始噪声精度参数α0的初始值设置为信号精度向量的初始值设置为子阵位置误差向量β的初始值设置为β＝0。
对每个阵元接收的信号添加方差σ为σ2的独立高斯白噪声,定义信噪比仿真的信噪比范围为0-10dB，在每个信噪比下的仿真次数 R＝200。
按照后验均值向量μ(t)公式(8)和后验协方差矩阵Σ公式(9)更新后验均值向量μ和后验协方差矩阵Σ，然后按照公式更新噪声精度参数计算公式(10)、更新方差σ的(σ2)new计算公式(11)和子阵位置误差向量β的估计值计算公式(12)更新更新初始噪声精度参数α0、信号精度向量α以及子阵位置误差向量β，重复以上过程，直到收敛，得到第i次仿真的位移误差向量估计值和方位估计值其中是信号精度向量α中值最小的K个分量的倒数组合而成的向量。
在给定信噪比下波达方向均方根误差，根据波达方向均方根误差计算公式计算，
位移误差均方根误差根据位移误差均方根误差计算公式计算，式中R表示仿真次数，P表示子阵数。
使用均方根误差计算公式计算得到如图2所示的位移误差估计值的均方根误差随信噪比的变化，可以看到位移误差估计量的均方根误差随信噪比的增加而减小。方位估计值的均方根误差随如图3所示的信噪比的变化，比使用多重信号分类MUSIC算法得到的均方根误差稳定1度附近。使用本发明方法计算得到的均方根误差均随信噪比的增加而减小，且均小于1 度，性能较MUSIC算法有较大提高。信噪比为0dB时使用全阵CBF算法、全阵MUSIC算法以及本发明方法的定位结果如图4所示，通过比较可以看出，本发明方法不仅能够得到准确的方位估计值，还具有比其他方法更高的角度分辨率。