CN103926579A - 一种开放式阵列信号的仿真方法 - Google Patents

Abstract

本发明公开了一种开放式阵列信号的仿真方法，该方法包括以下步骤：第一步：获取水声目标的实时的目标运动轨迹；第二步：获取水声目标辐射噪声信号；第三步：测算等间隔排列的线阵阵形畸变后各个阵元的坐标；第四步：生成实时的海洋信道脉冲响应序列；第五步：无阵元随机抖动情况下，生成阵元信号；第六步：出现随机抖动时，第i₇号阵元实际生成的阵元信号；第七步：获取阵元的实时海洋环境噪声信号；第八步：获得被动声纳阵列信号。该仿真方法使得阵列信号仿真系统具有开放性。

Description

一种开放式阵列信号的仿真方法

技术领域

本发明属于信号处理领域，具体来说，涉及一种开放式阵列信号的仿真方法。

背景技术

传统的水声阵列信号仿真系统结构单一。一方面，传统的水声目标运动轨迹仿真大都是对典型运动方式的仿真，如：直线航路、周(圆弧)航路和蛇形追踪航路。其中，直线航路是鱼雷攻击目标最常见的方式；圆周(圆弧)航路是目标调准航迹的主要方式之一，其主要特点是目标在水平面内作等半径的圆弧运动；蛇形追踪航路是潜艇跟踪和侦查敌方目标的一种主要方式，它包括两个部分：蛇形运动和追踪。但是这几种典型的运动方式都局限于平稳运动，运动方式单一，不能仿真非平稳目标运动轨迹；一方面，我国的海域辽阔海洋环境复杂多变，同一海域的海洋信道也会随气候洋流等的变换而不同，而传统的水声阵列信号仿真系统只能做到对特定海洋环境下海洋信道的仿真，显然这是不够的。一方面，一般的水声目标辐射噪声信号仿真根据谱特性将目标辐射噪声划分为线谱、时变调制谱和平稳连续谱三大谱分量分别进行建模仿真。若想使仿真信号逼近某一具体舰船的辐射信号则需对仿真参数进行精确设置，但有时并不具备这些参数信息；另一方面，随着现代水声计算的发展，对能够兼容实测数据并保留实测数据的特征的被动声纳阵列信号仿真系统的渴望越来越强烈。

发明内容

技术问题：本发明所要解决的技术问题是：一种开放式阵列信号的仿真方法，该仿真方法使得阵列信号仿真系统具有开放性。

技术方案：为解决上述技术问题，本发明采用的技术方案是：

一种开放式阵列信号的仿真方法，该仿真方法包括以下步骤：

第一步：获取水声目标的实时的目标运动轨迹：首先建立三维直角坐标系：以等间隔排列的线阵末端所在位置为圆心O，以船体拖着等间隔排列的线阵行驶的方向为X轴方向，等间隔排列的线阵扰动的方向为Y轴方向，垂直于海平面、且指向地心的方向为Z轴方向；设置开放式阵列信号仿真系统的抽样时间是T_s，采样率为F_s，水声目标的运动速度为V，实时的目标运动轨迹为n1表示整型变量；

获取水声目标的实时的目标运动轨迹的方法有：第一种方法：基于直线航路、圆周航路或者蛇形追踪航路的运动方式模拟水声目标的实时目标运动轨迹；或者第二种方法：基于B样条拟合的非平稳目标运动轨迹实时仿真方法，获得实时的目标运动轨迹；

第一种方法的过程为：当水声目标的运动方式是直线航路时，设置直线航路的初始位置为(x[0],y[0],z[0])^T，上标T表示转置；运动速度为V，运动的水平方向角为α，垂直仰角为θ，获得如式(1)所示的直线航路的实时的目标运动轨迹

$\stackrel{\→}{P}\left[n1\right]=\left(\begin{array}{c}x\left[n1\right]\\ y\left[n1\right]\\ z\left[n1\right]\end{array}\right)=\left(\begin{array}{c}x\left[0\right]+n1\·{\mathrm{VT}}_s\cos \left(\mathrm{\θ}\right)\cos \left(\mathrm{\α}\right)\\ y\left[0\right]+n1\·{\mathrm{VT}}_s\cos \left(\mathrm{\θ}\right)\sin \left(\mathrm{\α}\right)\\ z\left[0\right]+n1\·{\mathrm{VT}}_s\sin \left(\mathrm{\θ}\right)\end{array}\right)$    式(1)

其中，x[n1]表示实时的目标运动轨迹在X轴方向的坐标，y[n1]表示实时的目标运动轨迹在Y轴方向的坐标，z[n1]表示实时的目标运动轨迹在Z轴方向的坐标；

当水声目标的运动方式是圆周航路时，设置圆周航路的初始位置为(x[0],y[0],z[0])^T，速度为V，圆心上标T表示转置；圆周运动半径水声目标的盘旋角速度ω满足式(2)：

$\mathrm{\ω}=\frac{V}{r}=\frac{V}{\sqrt{{(X_0-x{\left[0\right]}^2)}^2+{(Y_0-y[0\left]\right)}^2}}$    式(2)

圆周航路的实时的目标运动轨迹如式(3)所示：

$\stackrel{\→}{p}\left[n1\right]=\left(\begin{array}{c}x\left[n1\right]\\ y\left[n1\right]\\ z\left[n1\right]\end{array}\right)=\left(\begin{array}{c}{X}_0+\left(x\right[n1-1]-X_0)\cos \left({\mathrm{\ωT}}_s\right)+\left(y\right[n1-1]-Y_0)\sin \left({\mathrm{\ωT}}_s\right)\\ Y_0+\left(y\right[n1-1]-Y_0)\cos \left({\mathrm{\ωT}}_s\right)-\left(x\right[n1-1]-X_0)\sin \left({\mathrm{\ωT}}_s\right)\\ Z_0\end{array}\right)$    式(3)

其中，x[n1-1]表示n1-1时刻的实时的目标运动轨迹在X轴方向的坐标，y[n1-1]表示n1-1时刻的实时的目标运动轨迹在Y轴方向的坐标；

当水声目标的运动方式是蛇形追踪航路时，设置蛇形追踪航路的运动速度为V，轨迹的水平跨度为L₀，拱高H，蛇形追踪航路的实时的目标运动轨迹如式(4)所示：

$\stackrel{\→}{P}\left[n1\right]=\stackrel{\→}{P}[n1-1]+\mathrm{\Δ}\stackrel{\→}{P}[n1-1]$    式(4)

其中，表示蛇形追踪航路n1-1时刻的运动坐标，表示蛇形追踪航路当前时刻和前一时刻，即在n1时刻和n1-1时刻的运动坐标的差值，满足式(5)：

$\mathrm{\Δ}\stackrel{\→}{P}[n1-1]={\mathrm{VT}}_s\left(\begin{array}{c}\frac{1}{\sqrt{1+{\left(\frac{2\mathrm{\πH}}{L_0}\cos \left(\frac{2\mathrm{\πx}[n1-1]}{L_0}\right)\right)}^2}}\\ \frac{\frac{2\mathrm{\πH}}{L_0}\cos \left(\frac{2\mathrm{\πx}[n1-1]}{L_0}\right)}{\sqrt{1+{\left(\frac{2\mathrm{\πH}}{L_0}\cos \left(\frac{2\mathrm{\πx}[n1-1]}{L_0}\right)\right)}^2}}\end{array}\right)$    式(5)

其中，n1表示整型变量，π表示圆周率；

第二种方法的过程为：获取采样率为f_s的实测水声目标运动轨迹序列m(n)：该运动轨迹序列m(n)包含N个运动轨迹点，待拟合的目标运动轨迹的采样率为F_s；利用基于B样条拟合的非平稳目标运动轨迹实时仿真方法，获得实时的目标运动轨迹；

第二步：获取水声目标辐射噪声信号：设水声目标辐射噪声信号由线谱分量、调制谱分量和连续谱分量组成，其中，线谱的数量为M3，第1—M3中的第q₄条线谱的幅度为第q₄条线谱的频率为和第q₄条线谱的随机相位为设调制谱分量为由M4个具有随机幅度、同样形状、相同重复周期的脉冲性随机过程组成，其中单个脉冲的形状取为高斯型，设置第1—M4中的第q₅个脉冲随机过程的调制深度为第q₅个脉冲随机过程的脉冲宽度为螺旋桨旋转的叶频为f_b，并且连续谱分量利用三参数模型模拟，然后进行AR建模，计算出AR滤波器的分子系数b₀和M5个AR滤波器分母系数中的第q₆个分母系数a_p(q₆)，依据式(6)获得实时的水声目标辐射噪声信号s(n2)为：

   式(6)

其中，x(n2)表示利用高斯白噪声进行AR建模，获得的n2时刻的水声目标辐射噪声信号的连续谱分量，x(n2-q₆)表示n2-q₆时刻的水声目标辐射噪声信号的连续谱分量；q₄、q₅和q₆为标号，q₄取1—M3，q₅取1—M4，q₆取1—M5，T_b＝1/f_b，v(n2)表示高斯白噪声，n2为整型变量，表示第q₅个脉冲随机过程，满足式(7)：

${\mathrm{\Λ}}_{q_5}(T_s\·n2-(q_5-1)\·T_b)=\frac{E_{q_5}}{\sqrt{2\mathrm{\π}}}{e}^{(-\frac{{(T_s\·n2-(q_5-1)\·T_b)}^2}{{2\mathrm{\σ}}_{q_5}})}$    式(7)

设水声目标辐射噪声信号仿真采用用于水声目标连续信号仿真的有限长数据循环拼接方法，获取有限长的实测水声目标数据序列s(n3)，序列长度为N0；利用用于水声目标连续信号仿真的有限长数据循环拼接方法获得实时的水声目标辐射噪声信号，n3为整型变量，n3取1—N0；

第三步：测算等间隔排列的线阵阵形畸变后各个阵元的当前坐标：首先获取等间隔排列的线阵的阵元间距d，阵元个数n10和声纳布放的深度z，从而测算等间隔排列的线阵无阵形畸变时各个阵元的坐标，记无阵形畸变时第1—n10号阵元中的第i₇号阵元的坐标为i₇为阵元标号，i₇取值为1—n10；其中，阵元的标号从坐标系的原点开始，按照X轴方向，顺序标号，从1标至n10号；然后，测算阵形畸变后各个阵元的当前坐标：在等间隔排列的线阵出现阵形畸变时，当阵形畸变是正弦随机扰动时，设置正弦随机扰动的扰动参数，正弦随机扰动的扰动参数包括正弦随机扰动的随机频率f₇、正弦随机扰动的随机幅度amp₇和标号为1的阵元的初始相位当阵形畸变是高斯随机扰动时，设置高斯随机扰动的扰动参数，高斯随机扰动的扰动参数包括高斯随机扰动的均值μ₇和高斯随机扰动的方差从而获得阵形畸变后第1—n10号阵元中的第i₇号阵元的当前坐标

第四步：生成实时的海洋信道脉冲响应序列：获取实测海洋环境参数，包括声速曲线、信号源的水平位置和垂直深度、波束的出射开角、波束传播的最远距离和最大深度、海底系数、海面系数和海洋环境的中心频率，利用基于射线模型的海洋信道仿真方法生成实时的海洋信道脉冲响应序列；

第五步：无阵元随机抖动情况下，生成阵元信号：将实时的水声目标辐射噪声信号和从水声目标的当前位置到达各个阵元当前坐标的实时的海洋信道脉冲响应序列相卷积，生成无阵元随机抖动情况下水声目标辐射噪声信号从水声目标位置传播到各个阵元的阵元信号k7为整型变量；实时的水声目标辐射噪声信号从第二步获得，水声目标的当前位置是从第一步获得的当前时刻的实时的目标运动轨迹，各个阵元当前坐标从第三步获得；

第六步：出现随机抖动时，第i₇号阵元实际生成的阵元信号：第1—n10号阵元中的第i₇号阵元出现随机抖动时，第i₇号阵元实际生成的阵元信号为：

${s^{\′}}_{i_7}\left(k7\right)=s_{i_7}\left(k7\right)(1-a_{i_7})+M_{i_7}$    式(8)

其中，表示第i₇号阵元随机抖动的均值偏移，表示第i₇号阵元随机抖动的幅度衰减；

当第i₇号阵元无随机抖动时，第i₇号阵元实际生成的阵元信号 ${s^{\′}}_{i_7}\left(k7\right)=s_{i_7}\left(k7\right);$

第七步：获取阵元的实时海洋环境噪声信号：假定海洋环境噪声的功率谱为p_e(f)，将p_e(f)通过傅里叶反变换测算出自相关序列，再利用Levinson-Durbin递归算法求解出AR滤波器的分子系数b_e0和M8个AR滤波器分母系数中的第q₈个分母系数a_ep(q₈)，依据式(9)获得第1—n10号阵元中的第i₇号阵元的实时海洋环境噪声信号

${\mathrm{se}}_{i_7}\left(n8\right)=-\underset{q_8=1}{\overset{M8}{\mathrm{\Σ}}}{a}_{\mathrm{ep}}\left(q_8\right)\·{\mathrm{se}}_{i_7}(n8-q_8)+b_{e0}v\left(n8\right)$    式(9)

其中，为n8-q₈时刻的实时海洋环境噪声信号，v(n8)表示第1—n10号阵元中的第i₇号阵元的实时海洋环境噪声信号对应的高斯白噪声，第1—n10号阵元的各个阵元的实时海洋环境噪声信号对应的高斯白噪声相互独立；

第八步：获得被动声纳阵列信号：被动声纳阵列信号的第1—n10号阵元中的第i7号阵元端接收到的信号是第i7号阵元实际生成的阵元信号和第i₇号阵元的实时海洋环境噪声信号的叠加。

进一步，所述的第一步中，基于B样条拟合的非平稳目标运动轨迹实时仿真方法，获得实时的目标运动轨迹的过程如下：

步骤101)，获得实测水声目标运动轨迹子序列：化简f_s/F_s为最简分子式，该最简分子式的分子为M，当M大于10时，则M0＝M，当M小于或等于10时，令Q＝2gM，g为从1开始的正整数，且Q为第一次大于10的数时，则M0＝Q；然后将实测水声目标运动轨迹序列m(n)从起始位置按长度M0截取，从而获得L段截取的实测水声目标运动轨迹子序列，最后一段子序列长度为M1，M1小于或等于M0，从第1段到第L-1段中的第i段子序列包含的运动轨迹点为m(M0(i-1)+1),m(M0(i-1)+2),…,m(M0(i-1)+M0)，第L段子序列包含的运动轨迹点为m(M0(L-1)+1),m(M0(L-1)+2),…,m(M0(L-1)+M1)；i表示子序列的段号，i为1至L；0<f_s<F_s；

步骤102)：获取各子序列反求三次B样条曲线的控制顶点：对每段目标运动轨迹子序列反求三次B样条曲线的控制顶点，当第L段子序列包含M0个运动轨迹点时，每段子序列的控制顶点为p_i(1),p_i(2),…,p_i(M0+1),p_i(M0+2)；当第L段子序列包含小于M0个运动轨迹点时，前L-1段子序列的控制顶点为p_i(1),p_i(2),…,p_i(M0+1),p_i(M0+2)，第L段子序列的控制顶点为p_L(1),p_L(2),…,p_L(M1+1),p_L(M1+2)；

步骤103)：获取各子序列在采样率F_s下的目标运动轨迹：利用获得的各段子序列的控制顶点分别作三次均匀B样条拟合，获得从第1段到第L-1段中的第i段子序列在采样率F_s下的目标运动轨迹m′_s(t₁)：

$m_i^{\&prime;}\left(t_1\right)=\frac{1}{6}\left[\begin{array}{cccc}1& t_1& t_1^2& t_1^3\end{array}\right]\left[\begin{array}{cccc}1& 4& 1& 0\\ -3& 0& 3& 0\\ 3& -6& 3& 0\\ -1& 3& -3& 1\end{array}\right]\left[\begin{array}{c}{p}_i\left(j_1\right)\\ p_i(j_1+1)\\ p_i(j_1+2)\\ p_i(j_1+3)\end{array}\right]$    式(10)

式(10)中，j₁＝1、2、…、M0-3,t₁＝k×f_s/F_s,t₁≤M0，k表示整型变量，j₁为大于t₁的最小整数；

当第L段子序列包含M0个运动轨迹点时，第L段子序列在采样率F_s下的目标运动轨迹m′_L,M0(t₂)：

$m_{L,M0}^{\&prime;}\left(t_2\right)=\frac{1}{6}\left[\begin{array}{cccc}1& t_2& t_2^2& t_2^3\end{array}\right]\left[\begin{array}{cccc}1& 4& 1& 0\\ -3& 0& 3& 0\\ 3& -6& 3& 0\\ -1& 3& -3& 1\end{array}\right]\left[\begin{array}{c}{p}_L\left(j_2\right)\\ p_L(j_2+1)\\ p_L(j_2+2)\\ p_L(j_2+3)\end{array}\right]$    式(11)

式(11)中，j₂＝1、2、…、M0-3,t₂＝k×f_s/F_s,t₂≤M0，k表示整型变量，j₂为大于t₂的最小整数；

当第L段子序列包含小于M0个运动轨迹点时，第L段子序列在采样率F_s下的目标运动轨迹m′_L,M1(t₃)：

$m_{L,M1}^{\&prime;}\left(t_3\right)=\frac{1}{6}\left[\begin{array}{cccc}1& t_3& t_3^2& t_3^3\end{array}\right]\left[\begin{array}{cccc}1& 4& 1& 0\\ -3& 0& 3& 0\\ 3& -6& 3& 0\\ -1& 3& -3& 1\end{array}\right]\left[\begin{array}{c}{p}_L\left(j_3\right)\\ p_L(j_3+1)\\ p_L(j_3+2)\\ p_L(j_3+3)\end{array}\right]$    式(12)

式(12)中，j₃＝1、2、…、M1-3,t₃＝k×f_s/F_s,t₃≤M1，k表示整型变量，j₃为大于t₃的最小整数；

步骤104)：获得实时的目标运动轨迹：将第1段到第L段的各段子序列在采样率F_s下的目标运动轨迹按段号由小到大合并，获得采样率F_s下的非平稳目标运动轨迹，从而获得实时的目标运动轨迹。

进一步，所述的步骤102)包括以下步骤：假设目标运动轨迹子序列的长度为M2，第1段到第L-1段中的第i段子序列对应的M2＝M0，当最后一段子序列长度为M1等于M0时，M2＝M0，当最后一段子序列长度M1小于M0时，M2＝M1，利用式(13)求解子序列的控制顶点p(1),p(2),…,p(M2+1),p(M2+2)：

(p(q),p(q+1),p(q+2))/6＝m(q)   式(13)

其中，q表示整型变量，q＝1,2,…,M2；m(q)表示子序列包含的目标运动轨迹点；式(13)方程组有M2个方程，M2+2个未知数，根据式(14)设定两个边界条件：

$\left\{\begin{array}{c}p\left(1\right)=p\left(2\right)\\ p(M2+1)=p(M2+2)\end{array}\right.$    式(14)

将式(13)和式(14)相结合，得到如式(15)所示的M2+2个线性方程：

$\left[\begin{array}{cccccc}-1& 1& & & & \\ 1& 4& 1& & & 0\\ & 1& 4& 1& & \\ & \&CenterDot;\&CenterDot;\&CenterDot;& & & & \\ & 0& & 1& 4& 1\\ & & & & 1& -1\end{array}\right]\left[\begin{array}{c}p\left(1\right)\\ p\left(2\right)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ p\left(M2\right)\\ p(M2+1)\\ p(M2+2)\end{array}\right]=6\left[\begin{array}{c}0\\ m\left(1\right)\\ m\left(2\right)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ m\left(M2\right)\\ 0\end{array}\right]$    式(15)

每段子序列通过上述方法得到A，从而得到该段子序列中A的伴随矩阵A^*，1≤q₁≤M2+2，1≤q₂≤M2+2，表示矩阵A的第q₁行第q₂列元的代数余子式，A^*满足式(16)和式(17)：

$\left\{\begin{array}{c}{A}_{1M2+2}={(-1)}^{M2+3}\\ A_{2M2+2}={(-1)}^{M2+5}\\ A_{q_{3}M2+2}=-4A_{q_3-1M2+2}-A_{q_3-2M2+2}\end{array}\right.$    式(16)

|A|＝-A_M2+2M2+2+A_M2+1M2+2   式(17)

其中，A_1M2+2表示矩阵A的第1行第M2+2列元的代数余子式，A_2M2+2表示矩阵A的第2行第M2+2列元的代数余子式，表示矩阵A的第q₃行第M2+2列元的代数余子式，表示矩阵A的第q₃-1行第M2+2列元的代数余子式，表示矩阵A的第q₃-2行第M2+2列元的代数余子式，A_M2+2M-2表示矩阵A的第M2+2行第M2+2列元的代数余子式，A_M2+1M2+2表示矩阵A的第M2+1行第M2+2列元的代数余子式；3≤q₃≤M2+2；

利用式(16)和式(17)得到A^*的第M2+2行元素和|A|，从而依据式(18)得到p(M2+2)：

$p(M2+2)=\left(6\underset{q=1}{\overset{M2}{\mathrm{\&Sigma;}}}{A}_{q+1M2+2}m\left(q\right)\right)/\left|A\right|$    式(18)

其中，A_q+1M2+2表示矩阵A的第q+1行第M2+2列元的代数余子式；

最后递推式(13)，得到M2+2个控制顶点p(1),p(2),…,p(M2+1),p(M2+2)，从而获得各段子序列的控制顶点。

进一步，所述的第二步中，三参数模型是用户通过设置三个参数f₀，f_m和λ来精确定义平稳连续谱p_c(f)的形状，p_c(f)满足：

$p_c\left(f\right)\frac{{\mathrm{\&sigma;}}_1^2}{2\mathrm{\&pi;}}[\frac{f_m+\mathrm{\&lambda;}(f+f_0)}{f_m^2+{(f+f_0)}^2}+\frac{f_m-\mathrm{\&lambda;}(f-f_0)}{f_m^2+{(f-f_0)}^2}]$    式(19)

其中，表示用来控制平稳连续谱信号的能量，f表示功率谱的线频率；

然后，将p_c(f)通过傅里叶反变换计算出自相关序列，再利用Levinson-Durbin递归算法求解AR滤波器的分子系数b₀和M5个AR滤波器分母系数中的第q₆个分母系数a_p(q₆)。

进一步，所述的第二步中，用于水声目标连续信号仿真的有限长数据循环拼接方法包括以下步骤：

首先：在实测数据段s(n3)首尾两端各找到一段能量最低的区域：从实测数据s(n3)的首尾部开始，在N_s点的搜索范围内以移窗的方式分别从i₁点开始搜索首部能量最低的区域，从i₂点开始搜索尾部能量最低的区域，窗长WL，记首尾部每段窗内的数据段的能量记为ph(i₁)和ph(i₂)：

$\mathrm{ph}\left(i_1\right)=\underset{n4=0}{\overset{\mathrm{WL}-1}{\mathrm{\&Sigma;}}}{s}^2(n4+i_1)$    式(20)

$\mathrm{ph}\left(i_2\right)=\underset{n5=0}{\overset{\mathrm{WL}-1}{\mathrm{\&Sigma;}}}{s}^2(n5+i_2)$    式(21)

其中，i₁表示搜索首部能量最低区域过程中，首部每次搜索的起点；i₂表示搜索尾部能量最低区域过程中，尾部每次搜索的起点；i₁∈[1,N_s-1]，i₂∈[N0-N_s,N0-WL]，WL表示每次搜索的长度；记窗头位置移到N₁时找到数据段首部能量最低的区域，记首部能量最低区域为W₁；记窗头位置移到N₂时找到数据段尾部能量最低的区域，记尾部能量最低区域为W₂；N0表示实测水声目标数据序列s(n3)的序列长度，s²(n₄+i₁)表示n4+i₁点对应的实测水声目标辐射噪声值的平方，s²(n₅+i₂)表示n5+i₂点对应的实测水声目标辐射噪声值的平方；

其次：在首尾两段能量最低区域中找出频谱特征相似的两段区域：以移窗的方式在区域W₁和W₂中分别从i₃点开始搜索首部频谱特性相似的的数据段W3，从i₄点开始搜索尾部频谱特性相似的数据段W₄，记窗长为WL₁，且WL₁<WL，最低能量区域W₁和W₂内每段窗长的数据段的频谱分别记为和

${\mathrm{Sh}}_{i_3}\left(e^{\mathrm{j\&omega;}}\right)=\underset{n5=0}{\overset{{\mathrm{WL}}_1-1}{\mathrm{\&Sigma;}}}s(n5+i_3)\&CenterDot;e^{-\mathrm{j\&omega;n}5}$    式(22)

${\mathrm{St}}_{i_4}\left(e^{\mathrm{j\&omega;}}\right)=\underset{n6=0}{\overset{{\mathrm{WL}}_1-1}{\mathrm{\&Sigma;}}}s(n6+i_4)\&CenterDot;e^{-\mathrm{j\&omega;n}6}$    式(23)

其中，i₃表示搜索首部最低能量区域中频谱特征相似的过程中，首部每次搜索的起点，i₄表示搜索尾部最低能量区域中频谱特征相似的过程中，尾部每次搜索的起点，i₃∈[N₁,N₁+WL-WL₁-1]，i₄∈[N₂,N₂+WL-WL₁-1]，j表示虚数单位，ω表示频谱的角频率，ω取值范围为-π—π；s(n5+i₃)表示n5+i₃点对应的实测水声目标辐射噪声值，s(n6+i₄)表示n6+i₄点对应的实测水声目标辐射噪声值；设W₁和W₂区域内窗头位置分别移到N₃和N₄，和最为接近时，则找到最低能量区域内频谱最为相似的两个数据段W₃和W₄；N₃表示和最为接近时，首部搜索的起点；N₄表示和最为接近时，尾部搜索的起点；

再次：在频谱特征相似的两段区域里找出导数同号的两组过零点：分别在区域W₃、W₄中找到符合条件的两组过零点s(N₅),s(N₅+1)和s(N₆),s(N₆+1)，满足：

s(N₅)·s(N₅+1)＜0,s(N₆)·s(N₆+1)＜0   式(24)

(s(N₅+1)-s(N₅))·(s(N₆+1)-s(N₆))＞0   式(25)

其中，N₅、N₅+1、N₆和N₆+1分别表示整型变量，s(N₅)表示N₅对应的实测水声目标辐射噪声值，s(N₅+1)表示N₅+1对应的实测水声目标辐射噪声值，s(N₆)表示N₆对应的实测水声目标辐射噪声值，s(N₆+1)表示N₆+1对应的实测水声目标辐射噪声值；

然后：确定循环数据段：实测信号中s(N₅)至s(N₆+1)之间的信号段，记为s_r(n7)，即满足局部相似条件的循环数据段，n7为整型变量；

最后：将循环数据段s_r(n7)不断重复拼接，获得实时的连续水声目标信号数据。

进一步，所述的第三步中，依据式(26)得到无阵形畸变时第1—n10号阵元中的第i₇号阵元的坐标

$\left\{\begin{array}{c}{x}_{i_7}=(i_7-1)d\\ y_{i_7}=0\\ z_{i_7}=z\end{array}\right.$    式(26)

其中，d表示等间隔排列的线阵的阵元间距，i₇表示阵元标号，z表示声纳布放的深度；

依据式(27)得到正弦随机扰动时第1—n10号阵元中的第i₇号阵元的当前坐标 $(x_{i_7}^{\&prime;},y_{i_7}^{\&prime;},z_{i_7}^{\&prime;}):$

   式(27)

依据式(28)得到高斯随机扰动时第1—n10号阵元中的第i₇号阵元的当前坐标 $(x_{i_7}^{\&prime;},y_{i_7}^{\&prime;},z_{i_7}^{\&prime;}):$

$\left\{\begin{array}{c}{x}_{i_7}^{\&prime;}(i_7-1)d\\ f\left(y_{i_7}^{\&prime;}\right)=\frac{1}{\sqrt{2\mathrm{\&pi;}{\mathrm{\&sigma;}}_7^2}}{e}^{(-\frac{{(y_{i_7}^{\&prime;}-{\mathrm{\&mu;}}_7)}^2}{2{\mathrm{\&sigma;}}_7^2})}\\ z_{i_7}^{\&prime;}=z\end{array}\right.$    式(28)

其中，表示第i₇号阵元的坐标的概率密度函数，满足高斯分布，e表示数学常数，为2.7182818。

进一步，所述的第四步中，利用基于射线模型的海洋信道仿真方法生成实时的海洋信道脉冲响应序列的步骤包括：

首先，基于实测海洋环境参数，在主机端利用射线模型的Bellhop软件包生成实测海洋环境下的海洋信道脉冲响应序列库，海洋信道脉冲响应序列库是水声目标辐射噪声信号从不同的水声目标位置传播到不同的阵元的本征声线，包括本征声线的数目，每条本征声线的幅度衰减、相移、传播时延、反射次数、折射次数和本征声线到达阵元时的入射角；

其次：获得实时的海洋信道脉冲响应序列：主机从DSP板读取水声目标实时的目标运动轨迹，每次读取水声目标实时的目标运动轨迹后，主机根据水声目标的当前坐标和各个阵元的当前坐标获取水声目标辐射噪声信号到达各个阵元的本征声线，测算当前的海洋信道脉冲响应序列，然后更新水声目标辐射噪声信号到达各个阵元的海洋信道脉冲响应序列为当前的海洋信道脉冲响应序列，再通过总线将当前的海洋信道脉冲响应序列加载进DSP板；通过主机每隔L1个采样间隔从DSP板读取一次水声目标运动坐标，在DSP板上获得实时的海洋信道脉冲响应序列；L1的取值范围为128—512之间的整数。

进一步，所述的第四步中，根据水声目标的当前坐标和各个阵元的坐标，主机获取水声目标辐射噪声信号到达各个阵元的本征声线，当从不同水声目标位置传播到不同阵元的本征声线小于m1条时，选出所有本征声线作为主机获取的水声目标辐射噪声信号到达各个阵元的本征声线，当本征声线大于或等于m1条时，将所有的本征声线按照幅度衰减降序排列，筛选出前m1条本征声线作为主机获取的水声目标辐射噪声信号到达各个阵元的本征声线，m1为10—20之间的整数。

进一步，所述的第四步中，测算当前的海洋信道脉冲响应序列的过程是：首先选取1号阵元为参考阵元，假设水声目标在运动过程中以球面波的形式不断的向周围发射辐射噪声信号，则水声目标辐射噪声信号到达各个阵元的传播时延ti₇为：

$t_{i_7}=({\mathrm{positon}}_{{\mathrm{array}}_{i7}}-{\mathrm{positon}}_{t\arg \mathrm{et}})/c$    式(29)

其中，为第i₇号阵元的坐标，positon_target为水声目标的当前坐标，c为海洋中声波传播速度；

各个阵元的信号到达时间相对于参考阵元的信号到达时间的时延差为：

$t_{i_7}-t_{i_0}=({\mathrm{positon}}_{{\mathrm{array}}_{i7}}-{\mathrm{positon}}_{t\arg \mathrm{et}})/c-({\mathrm{positon}}_{{\mathrm{array}}_{i0}}-{\mathrm{positon}}_{t\arg \mathrm{et}})/c$    式(30)

其中，下标i₀表示参考阵元的序号，表示参考阵元的坐标；

水声目标辐射噪声信号从当前坐标到达各个阵元的海洋信道脉冲响应序列 $h_{i_7}\left(n8\right)$ 为：

   式(31)

其中，i₇表示第i₇号阵元，k1表示第k1条本征声线，为第i₇号阵元对应的本征声线数目，表示第i₇号阵元的第k1条本征声线的幅度衰减，表示第i₇号阵元的第k1条本征声线的相移，表示第i₇号阵元的第k1条本征声线的传播时延；e表示数学常数，为2.7182818，j表示虚数单位，n8表示整型变量， $\mathrm{\&sigma;}(n8-{\mathrm{\&tau;}}_{i_{7}k1}-(t_{i_7}-t_{i_0}))$ 表示单位采样信号， $\mathrm{\&sigma;}(n8-{\mathrm{\&tau;}}_{i_{7}k1}-(t_{i_7}-t_{i_0}))$ 满足式(32)：

$\mathrm{\&sigma;}(n8-{\mathrm{\&tau;}}_{i_{7}k1}-(t_{i_7}-t_{i_0}))=\left\{\begin{array}{cc}1& n8={\mathrm{\&tau;}}_{i_{7}k1}-(t_{i_7}-t_{i_0})\\ 0& n8\&NotEqual;{\mathrm{\&tau;}}_{i_{7}k1}-(t_{i_7}-t_{i_0})\end{array}\right.$    式(32)。

进一步，所述的第六步中，第i₇号阵元随机抖动的幅度衰减服从均匀分布：

$f\left(a_{i_7}\right)=\frac{1}{B_7-A_7}$    式(33)

其中，表示幅度衰减的概率密度函数，A₇和B₇是均匀分布的两个参数，且 $A_7\&le;a_{i_7}\&le;B_7,$ 均匀分布的均值为 $\frac{A_7+B_7}{2},$ 均匀分布的方差为 $\frac{{(B_7-A_7)}^2}{12}.$

有益效果：与现有技术相比，本发明具有以下有益效果：

1、本发明的方法除了实现一般声纳阵列信号仿真，还基于实测数据的阵列信号仿真并保留实测数据的特征，具有开放性。本发明的仿真方法利用基于B样条拟合的非平稳目标运动轨迹实时仿真方法、用于水声目标连续信号仿真的有限长数据循环拼接方法、基于射线模型的海洋信道仿真方法以及正弦随机扰动和高斯阵形扰动表现声纳信号仿真的阵形畸变，基于均匀分布的阵元信号随机抖动表现声纳信号仿真的随机抖动。整体上做到阵列信号仿真系统的开放性：第一，目标运动的方式可选，包括利用直线航路、圆周航路和蛇形追踪航路的目标运动模型模拟目标运动和根据实测的目标运动轨迹拟合目标运动；第二，辐射噪声的信号源可选，包括应用辐射噪声的线谱分量、调制谱分量和连续谱分量实时产生噪声和利用实测海洋数据做低失真的实测噪声循环回放；第三，根据设置的实测海洋环境参数的不同，生成实测海洋环境下的海洋信道脉冲响应序列，这本身就是开放的；第四，声纳阵列嵌入阵形畸变和阵元随机抖动模块，这也与实际情况相符。

2、本发明的方法采用B样条曲线拟合，可以用于对非平稳目标运动轨迹的拟合，使得因目标运动状态改变而可能产生“拐点”的运动轨迹更加平滑。基于B样条拟合的非平稳目标运动轨迹实时仿真通过建模的方法拟合出符合实际要求的水声目标运动轨迹，特别是当实际记录的目标运动数据的采样率不高时，实现采样率F_s下的水声目标变速和变曲线运动轨迹的模拟，使得因目标运动状态改变而可能产生“拐点”的运动轨迹更加平滑，仿真系统将有更好的扩展性。目标运动状态改变就是指非平稳目标运动轨迹。产生“拐点”的运动轨迹越平滑，仿真的精度越高。

3、本发明的方法通过主机和DSP的交互完成对水声传播信道的实时仿真，且同时确保仿真精度。本发明的仿真方法选取n条本征声线条数来模拟海洋信道，以L个采样间隔为间隔更新海洋信道，实时生成的阵列信号基本保证了信号的能量和信号的连续性。即，在保证仿真精度的同时提高了海洋信道在DSP上实时仿真的速度，且保证实时仿真系统的采样率不低于6KHz，提升了海洋信道仿真的采样率，加快了仿真速度。

4、与现有技术相比，本发明的仿真方法更加精确。现有的声纳信号仿真不考虑阵形畸变和阵元随机抖动，忽视了在实际操作过程中，声纳基阵出现阵形畸变，阵元信号出现随机抖动的情况。本专利的仿真方法利用正弦随机扰动和高斯阵形扰动表现声纳信号仿真的阵形畸变，基于均匀分布的阵元信号随机抖动表现声纳信号仿真的随机抖动，使得声纳信号仿真方法更加精确。本发明的仿真方法实现对声纳信号仿真中阵形畸变和阵元抖动的仿真，丢弃声纳信号仿真中阵形不畸变和阵元不抖动的假设，使得声纳信号仿真更加真实和精确。

5、本发明的方法利用在TigerSharc TS201高速阵列信号处理板实时仿真阵列信号，具有工程实用性。

附图说明

图1所示为本发明的流程框图。

图2为实施例1中采样率f_s下的目标运动轨迹和基于三次均匀B样条算法拟合的采样率F_s下的目标运动轨迹。

图3为实施例1中采样率f_s和F_s下目标运动对应的阵元信号实部序列。

图4为实施例1中目标运动调整过程中2号阵元的阵元信号实部序列。

图5为实施例1中采样率fs和Fs下目标运动对应的2号阵元的阵元信号实部序列的功率谱。

图6为实施例1中等间隔排列的线阵加入正弦随机扰动的基阵坐标。

图7为实施例1中等间隔排列的线阵加入正弦随机扰动但无阵元随机抖动的2号阵元的阵元信号实部序列的功率谱。

图8为实施例1中等间隔排列的线阵加入正弦随机扰动和阵元随机抖动的1号阵元的阵元信号。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其它实施例，都属于本发明保护的范围。

如图1所示，本发明的一种开放式阵列信号的仿真方法，包括以下步骤：

第一步：获取水声目标的实时的目标运动轨迹：首先建立三维直角坐标系：以等间隔排列的线阵末端所在位置为圆心O，以船体拖着等间隔排列的线阵行驶的方向为X轴方向，等间隔排列的线阵扰动的方向为Y轴方向，垂直于海平面、且指向地心的方向为Z轴方向；设置开放式阵列信号仿真系统的抽样时间是T_s，采样率为F_s，水声目标的运动速度为V，实时的目标运动轨迹为n1表示整型变量；

获取水声目标的实时的目标运动轨迹的方法有：第一种方法：基于直线航路、圆周航路或者蛇形追踪航路的运动方式模拟水声目标的实时目标运动轨迹；或者第二种方法：基于B样条拟合的非平稳目标运动轨迹实时仿真方法，获得实时的目标运动轨迹；

第一种方法的过程为：当水声目标的运动方式是直线航路时，设置直线航路的初始位置为(x[0],y[0],z[0])^T，上标T表示转置；运动速度为V，运动的水平方向角为α，垂直仰角为θ，获得如式(1)所示的直线航路的实时的目标运动轨迹 $\stackrel{\&RightArrow;}{P}\left[n1\right]:$

$\stackrel{\&RightArrow;}{P}\left[n1\right]=\left(\begin{array}{c}x\left[n1\right]\\ y\left[n1\right]\\ z\left[n1\right]\end{array}\right)=\left(\begin{array}{c}x\left[0\right]+n1\&CenterDot;{\mathrm{VT}}_s\cos \left(\mathrm{\&theta;}\right)\cos \left(\mathrm{\&alpha;}\right)\\ y\left[0\right]+n1\&CenterDot;{\mathrm{VT}}_s\cos \left(\mathrm{\&theta;}\right)\sin \left(\mathrm{\&alpha;}\right)\\ z\left[0\right]+n1\&CenterDot;{\mathrm{VT}}_s\sin \left(\mathrm{\&theta;}\right)\end{array}\right)$    式(1)

其中，x[n1]表示实时的目标运动轨迹在X轴方向的坐标，y[n1]表示实时的目标运动轨迹在Y轴方向的坐标，z[n1]表示实时的目标运动轨迹在Z轴方向的坐标；

当水声目标的运动方式是圆周航路时，设置圆周航路的初始位置为(x[0],y[0],z[0])^T，速度为V，圆心上标T表示转置；圆周运动半径水声目标的盘旋角速度ω满足式(2)：

$\mathrm{\&omega;}=\frac{V}{r}=\frac{V}{\sqrt{{(X_0-x{\left[0\right]}^2)}^2+{(Y_0-y[0\left]\right)}^2}}$    式(2)

圆周航路的实时的目标运动轨迹如式(3)所示：

$\stackrel{\&RightArrow;}{p}\left[n1\right]=\left(\begin{array}{c}x\left[n1\right]\\ y\left[n1\right]\\ z\left[n1\right]\end{array}\right)=\left(\begin{array}{c}{X}_0+\left(x\right[n1-1]-X_0)\cos \left({\mathrm{\&omega;T}}_s\right)+\left(y\right[n1-1]-Y_0)\sin \left({\mathrm{\&omega;T}}_s\right)\\ Y_0+\left(y\right[n1-1]-Y_0)\cos \left({\mathrm{\&omega;T}}_s\right)-\left(x\right[n1-1]-X_0)\sin \left({\mathrm{\&omega;T}}_s\right)\\ Z_0\end{array}\right)$    式(3)

其中，x[n1-1]表示n1-1时刻的实时的目标运动轨迹在X轴方向的坐标，y[n1-1]表示n1-1时刻的实时的目标运动轨迹在Y轴方向的坐标；

当水声目标的运动方式是蛇形追踪航路时，设置蛇形追踪航路的运动速度为V，轨迹的水平跨度为L₀，拱高H，蛇形追踪航路的实时的目标运动轨迹如式(4)所示：

$\stackrel{\&RightArrow;}{P}\left[n1\right]=\stackrel{\&RightArrow;}{P}[n1-1]+\mathrm{\&Delta;}\stackrel{\&RightArrow;}{P}[n1-1]$    式(4)

其中，表示蛇形追踪航路n1-1时刻的运动坐标，表示蛇形追踪航路当前时刻和前一时刻，即在n1时刻和n1-1时刻的运动坐标的差值，满足式(5)：

$\mathrm{\&Delta;}\stackrel{\&RightArrow;}{P}[n1-1]={\mathrm{VT}}_s\left(\begin{array}{c}\frac{1}{\sqrt{1+{\left(\frac{2\mathrm{\&pi;H}}{L_0}\cos \left(\frac{2\mathrm{\&pi;x}[n1-1]}{L_0}\right)\right)}^2}}\\ \frac{\frac{2\mathrm{\&pi;H}}{L_0}\cos \left(\frac{2\mathrm{\&pi;x}[n1-1]}{L_0}\right)}{\sqrt{1+{\left(\frac{2\mathrm{\&pi;H}}{L_0}\cos \left(\frac{2\mathrm{\&pi;x}[n1-1]}{L_0}\right)\right)}^2}}\end{array}\right)$    式(5)

其中，n1表示整型变量，π表示圆周率。

第二种方法的过程为：获取采样率为f_s的实测水声目标运动轨迹序列m(n)：该运动轨迹序列m(n)包含N个运动轨迹点，待拟合的目标运动轨迹的采样率为F_s；利用基于B样条拟合的非平稳目标运动轨迹实时仿真方法，获得实时的目标运动轨迹。

在第一步中，基于B样条拟合的非平稳目标运动轨迹实时仿真方法，获得实时的目标运动轨迹的过程如下：

步骤101)，获得实测水声目标运动轨迹子序列：化简f_s/F_s为最简分子式，该最简分子式的分子为M，当M大于10时，则M0＝M，当M小于或等于10时，令Q＝2gM，g为从1开始的正整数，且Q为第一次大于10的数时，则M0＝Q；然后将实测水声目标运动轨迹序列m(n)从起始位置按长度M0截取，从而获得L段截取的实测水声目标运动轨迹子序列，最后一段子序列长度为M1，M1小于或等于M0，从第1段到第L-1段中的第i段子序列包含的运动轨迹点为m(M0(i-1)+1),m(M0(i-1)+2),…,m(M0(i-1)+M0)，第L段子序列包含的运动轨迹点为m(M0(L-1)+1),m(M0(L-1)+2),…,m(M0(L-1)+M1)；i表示子序列的段号，i为1至L；0<f_s<F_s；

步骤102)：获取各子序列反求三次B样条曲线的控制顶点：对每段目标运动轨迹子序列反求三次B样条曲线的控制顶点，当第L段子序列包含M0个运动轨迹点时，每段子序列的控制顶点为p_i(1),p_i(2),…,p_i(M0+1),p_i(M0+2)；当第L段子序列包含小于M0个运动轨迹点时，前L-1段子序列的控制顶点为p_i(1),p_i(2),…,p_i(M0+1),p_i(M0+2)，第L段子序列的控制顶点为p_L(1),p_L(2),…,p_L(M1+1),p_L(M1+2)；

步骤103)：获取各子序列在采样率F_s下的目标运动轨迹：利用获得的各段子序列的控制顶点分别作三次均匀B样条拟合，获得从第1段到第L-1段中的第i段子序列在采样率F_s下的目标运动轨迹m′_i(t₁)：

$m_i^{\&prime;}\left(t_1\right)=\frac{1}{6}\left[\begin{array}{cccc}1& t_1& t_1^2& t_1^3\end{array}\right]\left[\begin{array}{cccc}1& 4& 1& 0\\ -3& 0& 3& 0\\ 3& -6& 3& 0\\ -1& 3& -3& 1\end{array}\right]\left[\begin{array}{c}{p}_i\left(j_1\right)\\ p_i(j_1+1)\\ p_i(j_1+2)\\ p_i(j_1+3)\end{array}\right]$    式(10)

式(10)中，j₁＝1、2、…、M0-3,t₁＝k×f_s/F_s,t₁≤M0，k表示整型变量，j₁为大于t₁的最小整数；

当第L段子序列包含M0个运动轨迹点时，第L段子序列在采样率F_s下的目标运动轨迹m′_L,M0(t₂)：

$m_{L,M0}^{\&prime;}\left(t_2\right)=\frac{1}{6}\left[\begin{array}{cccc}1& t_2& t_2^2& t_2^3\end{array}\right]\left[\begin{array}{cccc}1& 4& 1& 0\\ -3& 0& 3& 0\\ 3& -6& 3& 0\\ -1& 3& -3& 1\end{array}\right]\left[\begin{array}{c}{p}_L\left(j_2\right)\\ p_L(j_2+1)\\ p_L(j_2+2)\\ p_L(j_2+3)\end{array}\right]$    式(11)

式(11)中，j₂＝1、2、…、M0-3,t₂＝k×f_s/F_s,t₂≤M0，k表示整型变量，j₂为大于t₂的最小整数；

当第L段子序列包含小于M0个运动轨迹点时，第L段子序列在采样率F_s下的目标运动轨迹m′_L,M1(t₃)：

$m_{L,M1}^{\&prime;}\left(t_3\right)=\frac{1}{6}\left[\begin{array}{cccc}1& t_3& t_3^2& t_3^3\end{array}\right]\left[\begin{array}{cccc}1& 4& 1& 0\\ -3& 0& 3& 0\\ 3& -6& 3& 0\\ -1& 3& -3& 1\end{array}\right]\left[\begin{array}{c}{p}_L\left(j_3\right)\\ p_L(j_3+1)\\ p_L(j_3+2)\\ p_L(j_3+3)\end{array}\right]$    式(12)

式(12)中，j₃＝1、2、…、M1-3,t₃＝k×f_s/F_s,t₃≤M1，k表示整型变量，j₃为大于t₃的最小整数；

步骤104)：获得实时的目标运动轨迹：将第1段到第L段的各段子序列在采样率F_s下的目标运动轨迹按段号由小到大合并，获得采样率F_s下的非平稳目标运动轨迹，从而获得实时的目标运动轨迹。

上述的步骤102)包括以下步骤：假设目标运动轨迹子序列的长度为M2，第1段到第L-1段中的第i段子序列对应的M2＝M0，当最后一段子序列长度为M1等于M0时，M2＝M0，当最后一段子序列长度M1小于M0时，M2＝M1，利用式(13)求解子序列的控制顶点p(1),p(2),…,p(M2+1),p(M2+2)：

(p(q),p(q+1),p(q+2))/6＝m(q)   式(13)

其中，q表示整型变量，q＝1,2,…,M2；m(q)表示子序列包含的目标运动轨迹点；式(13)方程组有M2个方程，M2+2个未知数，根据式(14)设定两个边界条件：

$\left\{\begin{array}{c}p\left(1\right)=p\left(2\right)\\ p(M2+1)=p(M2+2)\end{array}\right.$    式(14)

将式(13)和式(14)相结合，得到如式(15)所示的M2+2个线性方程：

$\left[\begin{array}{cccccc}-1& 1& & & & \\ 1& 4& 1& & & 0\\ & 1& 4& 1& & \\ & \&CenterDot;\&CenterDot;\&CenterDot;& & & & \\ & 0& & 1& 4& 1\\ & & & & 1& -1\end{array}\right]\left[\begin{array}{c}p\left(1\right)\\ p\left(2\right)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ p\left(M2\right)\\ p(M2+1)\\ p(M2+2)\end{array}\right]=6\left[\begin{array}{c}0\\ m\left(1\right)\\ m\left(2\right)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ m\left(M2\right)\\ 0\end{array}\right]$    式(15)

每段子序列通过上述方法得到A，从而得到该段子序列中A的伴随矩阵A^*，1≤q₁≤M2+2，1≤q₂≤M2+2，表示矩阵A的第q₁行第q₂列元的代数余子式，A^*满足式(16)和式(17)：

$\left\{\begin{array}{c}{A}_{1M2+2}={(-1)}^{M2+3}\\ A_{2M2+2}={(-1)}^{M2+5}\\ A_{q_{3}M2+2}=-4A_{q_3-1M2+2}-A_{q_3-2M2+2}\end{array}\right.$    式(16)

|A|＝-A_M2+2M2+2+A_M2+1M2+2   式(17)

其中，A_1M2+2表示矩阵A的第1行第M2+2列元的代数余子式，A_2M2+2表示矩阵A的第2行第M2+2列元的代数余子式，表示矩阵A的第q₃行第M2+2列元的代数余子式，表示矩阵A的第q₃-1行第M2+2列元的代数余子式，表示矩阵A的第q₃-2行第M2+2列元的代数余子式，A_M2+2M-2表示矩阵A的第M2+2行第M2+2列元的代数余子式，A_M2+1M2+2表示矩阵A的第M2+1行第M2+2列元的代数余子式；3≤q₃≤M2+2；

利用式(16)和式(17)得到A^*的第M2+2行元素和|A|，从而依据式(18)得到p(M2+2)：

$p(M2+2)=\left(6\underset{q=1}{\overset{M2}{\mathrm{\&Sigma;}}}{A}_{q+1M2+2}m\left(q\right)\right)/\left|A\right|$    式(18)

其中，A_q+1M2+2表示矩阵A的第q+1行第M2+2列元的代数余子式；

最后递推式(13)，得到M2+2个控制顶点p(1),p(2),…,p(M2+1),p(M2+2)，从而获得各段子序列的控制顶点。

第二步：获取水声目标辐射噪声信号：设水声目标辐射噪声信号由线谱分量、调制谱分量和连续谱分量组成，其中，线谱的数量为M3，第1—M3中的第q₄条线谱的幅度为第q₄条线谱的频率为和第q₄条线谱的随机相位为设调制谱分量为由M4个具有随机幅度、同样形状、相同重复周期的脉冲性随机过程组成，其中单个脉冲的形状取为高斯型，设置第1—M4中的第q₅个脉冲随机过程的调制深度为第q₅个脉冲随机过程的脉冲宽度为螺旋桨旋转的叶频为f_b，并且连续谱分量利用三参数模型模拟，然后进行AR建模，计算出AR滤波器的分子系数b₀和M5个AR滤波器分母系数中的第q₆个分母系数a_p(q₆)，依据式(6)获得实时的水声目标辐射噪声信号s(n2)为：

   式(6)

其中，x(n2)表示利用高斯白噪声进行AR建模，获得的n2时刻的水声目标辐射噪声信号的连续谱分量，x(n2-q₆)表示n2-q₆时刻的水声目标辐射噪声信号的连续谱分量；q₄、q₅和q₆为标号，q₄取1—M3，q₅取1—M4，q₆取1—M5，T_b＝1/f_b，v(n2)表示高斯白噪声，n2为整型变量，表示第q₅个脉冲随机过程，满足式(7)：

${\mathrm{\&Lambda;}}_{q_5}(T_s\&CenterDot;n2-(q_5-1)\&CenterDot;T_b)=\frac{E_{q_5}}{\sqrt{2\mathrm{\&pi;}}}{e}^{(-\frac{{(T_s\&CenterDot;n2-(q_5-1)\&CenterDot;T_b)}^2}{{2\mathrm{\&sigma;}}_{q_5}})}$    式(7)

设水声目标辐射噪声信号仿真采用用于水声目标连续信号仿真的有限长数据循环拼接方法，获取有限长的实测水声目标数据序列s(n3)，序列长度为N0；利用用于水声目标连续信号仿真的有限长数据循环拼接方法获得实时的水声目标辐射噪声信号，n3为整型变量，n3取1—N0。

在第二步中，三参数模型是用户通过设置三个参数f₀，f_m和λ来精确定义平稳连续谱p_c(f)的形状，p_c(f)满足：

$p_c\left(f\right)\frac{{\mathrm{\&sigma;}}_1^2}{2\mathrm{\&pi;}}[\frac{f_m+\mathrm{\&lambda;}(f+f_0)}{f_m^2+{(f+f_0)}^2}+\frac{f_m-\mathrm{\&lambda;}(f-f_0)}{f_m^2+{(f-f_0)}^2}]$    式(19)

其中，表示用来控制平稳连续谱信号的能量，f表示功率谱的线频率；

然后，将p_c(f)通过傅里叶反变换计算出自相关序列，再利用Levinson-Durbin递归算法求解AR滤波器的分子系数b₀和M5个AR滤波器分母系数中的第q₆个分母系数a_p(q₆)。

在第二步中，用于水声目标连续信号仿真的有限长数据循环拼接方法包括以下步骤：

首先：在实测数据段s(n3)首尾两端各找到一段能量最低的区域：从实测数据s(n3)的首尾部开始，在N_s点的搜索范围内以移窗的方式分别从i₁点开始搜索首部能量最低的区域，从i₂点开始搜索尾部能量最低的区域，窗长WL，记首尾部每段窗内的数据段的能量记为ph(i₁)和ph(i₂)：

   式(20)

$\mathrm{ph}\left(i_2\right)=\underset{n5=0}{\overset{\mathrm{WL}-1}{\mathrm{\&Sigma;}}}{s}^2(n5+i_2)$    式(21)

其中，i₁表示搜索首部能量最低区域过程中，首部每次搜索的起点；i₂表示搜索尾部能量最低区域过程中，尾部每次搜索的起点；i₁∈[1,N_s-1]，i₂∈[N0-N_s,N0-WL]，WL表示每次搜索的长度；记窗头位置移到N₁时找到数据段首部能量最低的区域，记首部能量最低区域为W₁；记窗头位置移到N₂时找到数据段尾部能量最低的区域，记尾部能量最低区域为W₂；N0表示实测水声目标数据序列s(n3)的序列长度，s²(n₄+i₁)表示n4+i₁点对应的实测水声目标辐射噪声值的平方，s²(n₅+i₂)表示n5+i₂点对应的实测水声目标辐射噪声值的平方；

其次：在首尾两段能量最低区域中找出频谱特征相似的两段区域：以移窗的方式在区域W₁和W₂中分别从i₃点开始搜索首部频谱特性相似的的数据段W₃，从i₄点开始搜索尾部频谱特性相似的数据段W₄，记窗长为WL₁，且WL₁<WL，最低能量区域W₁和W₂内每段窗长的数据段的频谱分别记为和

${\mathrm{Sh}}_{i_3}\left(e^{\mathrm{j\&omega;}}\right)=\underset{n5=0}{\overset{{\mathrm{WL}}_1-1}{\mathrm{\&Sigma;}}}s(n5+i_3)\&CenterDot;e^{-\mathrm{j\&omega;n}5}$    式(22)

${\mathrm{St}}_{i_4}\left(e^{\mathrm{j\&omega;}}\right)=\underset{n6=0}{\overset{{\mathrm{WL}}_1-1}{\mathrm{\&Sigma;}}}s(n6+i_4)\&CenterDot;e^{-\mathrm{j\&omega;n}6}$    式(23)

其中，i₃表示搜索首部最低能量区域中频谱特征相似的过程中，首部每次搜索的起点，i₄表示搜索尾部最低能量区域中频谱特征相似的过程中，尾部每次搜索的起点，i₃∈[N₁,N₁+WL-WL₁-1]，i₄∈[N₂,N₂+WL-WL₁-1]，j表示虚数单位，ω表示频谱的角频率，ω取值范围为-π—π；s(n5+i₃)表示n5+i₃点对应的实测水声目标辐射噪声值，s(n6+i₄)表示n6+i₄点对应的实测水声目标辐射噪声值；设W₁和W₂区域内窗头位置分别移到N₃和N₄，和最为接近时，则找到最低能量区域内频谱最为相似的两个数据段W3和W4；N₃表示和最为接近时，首部搜索的起点；N₄表示和最为接近时，尾部搜索的起点；

再次：在频谱特征相似的两段区域里找出导数同号的两组过零点：分别在区域W₃、W₄中找到符合条件的两组过零点s(N₅),s(N₅+1)和s(N₆),s(N₆+1)，满足：

s(N₅)·s(N₅+1)＜0,s(N₆)·s(N₆+1)＜0   式(24)

(s(N₅+1)-s(N₅))·(s(N₆+1)-s(N₆))＞0   式(25)

其中，N₅、N₅+1、N₆和N₆+1分别表示整型变量，s(N₅)表示N₅对应的实测水声目标辐射噪声值，s(N₅+1)表示N₅+1对应的实测水声目标辐射噪声值，s(N₆)表示N₆对应的实测水声目标辐射噪声值，s(N₆+1)表示N₆+1对应的实测水声目标辐射噪声值；

然后：确定循环数据段：实测信号中s(N₅)至s(N₆+1)之间的信号段，记为s_r(n7)，即满足局部相似条件的循环数据段，n7为整型变量；

最后：将循环数据段s_r(n7)不断重复拼接，获得实时的连续水声目标信号数据。

第三步：测算等间隔排列的线阵阵形畸变后各个阵元的当前坐标：首先获取等间隔排列的线阵的阵元间距d，阵元个数n10和声纳布放的深度z，从而测算等间隔排列的线阵无阵形畸变时各个阵元的坐标，记无阵形畸变时第1—n10号阵元中的第i₇号阵元的坐标为i₇为阵元标号，i₇取值为1—n10；其中，阵元的标号从坐标系的原点开始，按照X轴方向，顺序标号，从1标至n10号；然后，测算阵形畸变后各个阵元的当前坐标：在等间隔排列的线阵出现阵形畸变时，当阵形畸变是正弦随机扰动时，设置正弦随机扰动的扰动参数，正弦随机扰动的扰动参数包括正弦随机扰动的随机频率f₇、正弦随机扰动的随机幅度amp7和标号为1的阵元的初始相位当阵形畸变是高斯随机扰动时，设置高斯随机扰动的扰动参数，高斯随机扰动的扰动参数包括高斯随机扰动的均值μ₇和高斯随机扰动的方差从而获得阵形畸变后第1—n10号阵元中的第i₇号阵元的当前坐标 $(x_{i_7}^{\&prime;},y_{i_7}^{\&prime;},z_{i_7}^{\&prime;}).$

在第三步中，依据式(26)得到无阵形畸变时第1—n10号阵元中的第i₇号阵元的坐标 $(x_{i_7},y_{i_7},z_{i_7}):$

$\left\{\begin{array}{c}{x}_{i_7}=(i_7-1)d\\ y_{i_7}=0\\ z_{i_7}=z\end{array}\right.$    式(26)

其中，d表示等间隔排列的线阵的阵元间距，i₇表示阵元标号，z表示声纳布放的深度；

依据式(27)得到正弦随机扰动时第1—n10号阵元中的第i₇号阵元的当前坐标 $(x_{i_7}^{\&prime;},y_{i_7}^{\&prime;},z_{i_7}^{\&prime;}):$

   式(27)

依据式(28)得到高斯随机扰动时第1—n10号阵元中的第i₇号阵元的当前坐标 $(x_{i_7}^{\&prime;},y_{i_7}^{\&prime;},z_{i_7}^{\&prime;}):$

$\left\{\begin{array}{c}{x}_{i_7}^{\&prime;}(i_7-1)d\\ f\left(y_{i_7}^{\&prime;}\right)=\frac{1}{\sqrt{2\mathrm{\&pi;}{\mathrm{\&sigma;}}_7^2}}{e}^{(-\frac{{(y_{i_7}^{\&prime;}-{\mathrm{\&mu;}}_7)}^2}{2{\mathrm{\&sigma;}}_7^2})}\\ z_{i_7}^{\&prime;}=z\end{array}\right.$    式(28)

其中，表示第i₇号阵元的坐标的概率密度函数，满足高斯分布，e表示数学常数，为2.7182818。

第四步：生成实时的海洋信道脉冲响应序列：获取实测海洋环境参数，包括声速曲线、信号源的水平位置和垂直深度、波束的出射开角、波束传播的最远距离和最大深度、海底系数、海面系数和海洋环境的中心频率，利用基于射线模型的海洋信道仿真方法生成实时的海洋信道脉冲响应序列。

在第四步中，利用基于射线模型的海洋信道仿真方法生成实时的海洋信道脉冲响应序列的步骤包括：

首先，基于实测海洋环境参数，在主机端利用射线模型的Bellhop软件包生成实测海洋环境下的海洋信道脉冲响应序列库，海洋信道脉冲响应序列库是水声目标辐射噪声信号从不同的水声目标位置传播到不同的阵元的本征声线，包括本征声线的数目，每条本征声线的幅度衰减、相移、传播时延、反射次数、折射次数和本征声线到达阵元时的入射角；

其次：获得实时的海洋信道脉冲响应序列：主机从DSP板读取水声目标实时的目标运动轨迹，每次读取水声目标实时的目标运动轨迹后，主机根据水声目标的当前坐标和各个阵元的当前坐标获取水声目标辐射噪声信号到达各个阵元的本征声线，测算当前的海洋信道脉冲响应序列，然后更新水声目标辐射噪声信号到达各个阵元的海洋信道脉冲响应序列为当前的海洋信道脉冲响应序列，再通过总线将当前的海洋信道脉冲响应序列加载进DSP板；通过主机每隔L1个采样间隔从DSP板读取一次水声目标运动坐标，在DSP板上获得实时的海洋信道脉冲响应序列；L1的取值范围为128—512之间的整数。

在第四步中，根据水声目标的当前坐标和各个阵元的坐标，主机获取水声目标辐射噪声信号到达各个阵元的本征声线，当从不同水声目标位置传播到不同阵元的本征声线小于m1条时，选出所有本征声线作为主机获取的水声目标辐射噪声信号到达各个阵元的本征声线，当本征声线大于或等于m1条时，将所有的本征声线按照幅度衰减降序排列，筛选出前m1条本征声线作为主机获取的水声目标辐射噪声信号到达各个阵元的本征声线，m1为10—20之间的整数。

在第四步中，测算当前的海洋信道脉冲响应序列的过程是：首先选取1号阵元为参考阵元，假设水声目标在运动过程中以球面波的形式不断的向周围发射辐射噪声信号，则水声目标辐射噪声信号到达各个阵元的传播时延为：

$t_{i_7}=({\mathrm{positon}}_{{\mathrm{array}}_{i7}}-{\mathrm{positon}}_{t\arg \mathrm{et}})/c$    式(29)

其中，为第i₇号阵元的坐标，positon_target为水声目标的当前坐标，c为海洋中声波传播速度；

各个阵元的信号到达时间相对于参考阵元的信号到达时间的时延差为：

$t_{i_7}-t_{i_0}=({\mathrm{positon}}_{{\mathrm{array}}_{i7}}-{\mathrm{positon}}_{t\arg \mathrm{et}})/c-({\mathrm{positon}}_{{\mathrm{array}}_{i0}}-{\mathrm{positon}}_{t\arg \mathrm{et}})/c$    式(30)

其中，下标i₀表示参考阵元的序号，表示参考阵元的坐标；

水声目标辐射噪声信号从当前坐标到达各个阵元的海洋信道脉冲响应序列 $h_{i_7}\left(n8\right)$ 为：

   式(31)

其中，i₇表示第i₇号阵元，k1表示第k1条本征声线，为第i₇号阵元对应的本征声线数目，表示第i₇号阵元的第k1条本征声线的幅度衰减，表示第i₇号阵元的第k1条本征声线的相移，表示第i₇号阵元的第k1条本征声线的传播时延；e表示数学常数，为2.7182818，j表示虚数单位，n8表示整型变量， $\mathrm{\&sigma;}(n8-{\mathrm{\&tau;}}_{i_{7}k1}-(t_{i_7}-t_{i_0}))$ 表示单位采样信号， $\mathrm{\&sigma;}(n8-{\mathrm{\&tau;}}_{i_{7}k1}-(t_{i_7}-t_{i_0}))$ 满足式(32)：

$\mathrm{\&sigma;}(n8-{\mathrm{\&tau;}}_{i_{7}k1}-(t_{i_7}-t_{i_0}))=\left\{\begin{array}{cc}1& n8={\mathrm{\&tau;}}_{i_{7}k1}-(t_{i_7}-t_{i_0})\\ 0& n8\&NotEqual;{\mathrm{\&tau;}}_{i_{7}k1}-(t_{i_7}-t_{i_0})\end{array}\right.$    式(32)。

第五步：无阵元随机抖动情况下，生成阵元信号：将实时的水声目标辐射噪声信号和从水声目标的当前位置到达各个阵元当前坐标的实时的海洋信道脉冲响应序列相卷积，生成无阵元随机抖动情况下水声目标辐射噪声信号从水声目标位置传播到各个阵元的阵元信号k7为整型变量；实时的水声目标辐射噪声信号从第二步获得，水声目标的当前位置是从第一步获得的当前时刻的实时的目标运动轨迹，各个阵元当前坐标从第三步获得。

第六步：出现随机抖动时，第i₇号阵元实际生成的阵元信号：第1—n10号阵元中的第i₇号阵元出现随机抖动时，第i₇号阵元实际生成的阵元信号为：

${s^{\&prime;}}_{i_7}\left(k7\right)=s_{i_7}\left(k7\right)(1-a_{i_7})+M_{i_7}$    式(8)

其中，表示第i₇号阵元随机抖动的均值偏移，表示第i₇号阵元随机抖

动的幅度衰减；

当第i7号阵元无随机抖动时，第i7号阵元实际生成的阵元信号 ${s^{\&prime;}}_{i_7}\left(k7\right)=s_{i_7}\left(k7\right).$

在第六步中，第i₇号阵元随机抖动的幅度衰减服从均匀分布：

$f\left(a_{i_7}\right)=\frac{1}{B_7-A_7}$    式(33)

其中，表示幅度衰减的概率密度函数，A₇和B₇是均匀分布的两个参数，且 $A_7\&le;a_{i_7}\&le;B_7,$ 均匀分布的均值为 $\frac{A_7+B_7}{2},$ 均匀分布的方差为 $\frac{{(B_7-A_7)}^2}{12}.$

第七步：获取阵元的实时海洋环境噪声信号：假定海洋环境噪声的功率谱为p_e(f)，将p_e(f)通过傅里叶反变换测算出自相关序列，再利用Levinson-Durbin递归算法求解出AR滤波器的分子系数b_e0和M8个AR滤波器分母系数中的第q₈个分母系数a_ep(q₈)，依据式(9)获得第1—n10号阵元中的第i₇号阵元的实时海洋环境噪声信号

${\mathrm{se}}_{i_7}\left(n8\right)=-\underset{q_8=1}{\overset{M8}{\mathrm{\&Sigma;}}}{a}_{\mathrm{ep}}\left(q_8\right)\&CenterDot;{\mathrm{se}}_{i_7}(n8-q_8)+b_{e0}v\left(n8\right)$    式(9)

其中，为n8-q₈时刻的实时海洋环境噪声信号，v(n8)表示第1—n10号阵元中的第i₇号阵元的实时海洋环境噪声信号对应的高斯白噪声，第1—n10号阵元的各个阵元的实时海洋环境噪声信号对应的高斯白噪声相互独立。

第八步：获得被动声纳阵列信号：被动声纳阵列信号的第1—n10号阵元中的第i7号阵元端接收到的信号是第i7号阵元实际生成的阵元信号和第i₇号阵元的实时海洋环境噪声信号的叠加。

本发明利用开放式目标轨迹仿真、开放式目标辐射噪声仿真、开放式环境噪声仿真、开放式水声传播信道仿真和开放式阵元信号仿真构建了开放式阵列信号仿真方法。

下面例举一实施例。

现有实测浅海海洋环境参数：海底声速：1650m/s，海底密度1.6g/cm³，海底吸收系数0.59dB/Wavelength，海深100m。在该实测海洋环境下，某等间隔排列的线阵，阵元间距0.25m，阵元个数36个，布放在距离海平面30m深的海水中。由于受到拖船、海流等的影响，阵形发生畸变；由于声纳基阵在生产工艺和在安装过程中出现了误差，第1号阵元出现随机抖动，第1号阵元的阵元信号出现0.2v的均值偏移，幅度衰减在0—0.2左右。水声目标辐射噪声信号假定为单频信号，采样率为f_s＝6KHz，单频信号的频率为300Hz。在上述条件下，以采样率为f_s＝10Hz模拟这一水声目标运动，目标首先沿45°方向做直线加速运动，继而转变运动方向沿270°方向做直线加速运动，轨迹序列长度为6，如图2所示。另外待拟合的目标运动轨迹的采样率为F_s＝6KHz。

首先，简化f_s/F_s为最简分子式，最简式的分子为1，然后按第一步介绍的方法选取M＝16，然后将实测水声目标运动轨迹序列m(i)按长度M＝16截取，第一段取第1～16采样点的运动轨迹，第二段取16～31采样点的运动轨迹，依次分段截取，最后一段序列的长度为6。

接着，对每段截取的运动轨迹补充边界条件方程，给出反求控制顶点的线性方程组。归纳总结A的伴随矩阵的递推公式和A行列式的计算公式。根据递推公式求出各段轨迹序列对应的控制顶点。

接着，利用各段轨迹序列的控制顶点作三次均匀B样条拟合，模拟出采样率F_s＝6KHz的目标运动轨迹，如图2所示。

利用采样率f_s和采样率F_s下的目标运动轨迹分别仿真阵元信号，整个运动过程的阵元信号如图3所示，其调整运动方向过程的阵元信号如图4所示。从图2可以明显的看出，目标调整运动方向过程中采样率f_s下的目标运动轨迹有“拐点”，而采样率F_s下的目标运动轨迹比较平滑。从图4可以看出运动调整过程的采样率F_s下的阵元信号比采样率f_s下的阵元信号更光滑，没有出现采样率f_s下出现的信号跳变。如图5中看出，采样率f_s下的功率谱图毛刺较多，而采样率F_s下的阵元信号实部的功率谱图相比采样率f_s下的功率谱图更加光滑。产生“拐点”的运动轨迹越平滑，说明仿真的精度越高。

接着，在主机加载实测海洋环境参数，包括信号源的水平位置500m和垂直深度60m，波束出射开角为-20°—20°，波束传播的最远距离为50km，最大深度为100m，海洋环境的中心频率3KHz，以及海底声速、海底密度、海底吸收系数和海洋声速曲线等。同时，在DSP端基于采样率F_s下的目标运动轨迹序列实时仿真水声目标的运动轨迹，在主机端利用Bellhop软件包生成实测海洋环境下的海洋信道脉冲响应序列库。

同时，测算等间隔排列的线阵无阵形畸变时第1～36号阵元的各个阵元的坐标，测算出无阵形畸变时各个阵元的坐标如图6所示。接着，选取阵形畸变为正弦随机扰动，设置正弦随机扰动的扰动参数，正弦随机扰动的扰动参数包括正弦随机扰动的随机频率f₇为0.05Hz、正弦随机扰动的随机幅度A为1v和标号为1的阵元的初始相位为0°。利用阵形畸变的参数测算出阵形畸变后各个阵元的当前坐标，如图5所示。

接着，设置声纳基阵的标号为1的阵元出现随机抖动，随机抖动的均值偏移为0.2v，幅度衰减a₁服从均匀分布：0≤a₁≤0.2；标号为2—36的阵元无随机抖动。

然后，主机每隔256个采样点从DSP板读取一次目标运动坐标，主机根据水声目标的当前坐标和各个阵元的当前坐标获取水声目标辐射噪声信号到达各个阵元的本征声线，当从不同水声目标位置传播到不同阵元的本征声线小于10条时，选出所有本征声线，当本征声线大于或等于10条时，将所有的本征声线按照幅度衰减降序排列，筛选出前10条本征声线，测算当前的海洋信道脉冲响应序列，然后更新水声目标辐射噪声信号到达各个阵元的海洋信道脉冲响应序列为当前的海洋信道脉冲响应序列。再通过总线将求得的海洋信道脉冲响应序列加载进DSP板。将水声目标实时发出的目标辐射噪声信号和海洋信道脉冲响应序列相卷积，获得实时仿真的被动声纳各个阵元的阵元信号，图7是无阵元随机抖动的2号阵元的阵元信号的功率谱图。对实时仿真的2号阵元的阵元信号分析研究，发现选取10条本征声线来模拟海洋信道，以256个采样间隔更新海洋信道，生成的阵列信号基本保证了信号的能量和信号的连续性，阵元信号实部的功率谱如图7所示。从图7看出阵元信号实部的功率谱只在300Hz频率附近有峰值，在其他频率功率谱几乎为0。也就是说，阵元信号实部的功率谱在300Hz频率附近以外的频率没有毛刺。这恰恰反映了阵元信号的连续性。另外，由于标号为1的阵元出现随机抖动，标号为1的阵元实际生成的阵元信号依据式(8)获得，标号为1的阵元在随机抖动前后，1号阵元的阵元信号如图8所示。从图8看出，第1号阵元出现随机抖动时，第1号阵元的阵元信号出现0.2v的均值偏移和0.1的幅度衰减。该数值符合本实施例开始设定的第1号阵元出现随机抖动时，第1号阵元的阵元信号出现0.2v的均值偏移，幅度衰减在0—0.2左右。第1号阵元的仿真效果达到阵元抖动的实际效果。

最后，设置海洋环境参数，计算出各个阵元的海洋环境噪声，同时叠加各个阵元实际生成的阵元信号，获得的就是被动声纳阵列信号。在DSP板利用DSP的0号定时器测试单位时间内生成的被动声纳阵列信号的序列长度为6000，满足采样率f_s＝6KHz的仿真速度要求。

以上对本发明实施例所提供的一种开放式阵列信号仿真系统，进行了详细介绍，本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。

Claims (10)

1.一种开放式阵列信号的仿真方法，其特征在于，该仿真方法包括以下步骤：

第一步：获取水声目标的实时的目标运动轨迹：首先建立三维直角坐标系：以等间隔排列的线阵末端所在位置为圆心O，以船体拖着等间隔排列的线阵行驶的方向为X轴方向，等间隔排列的线阵扰动的方向为Y轴方向，垂直于海平面、且指向地心的方向为Z轴方向；设置开放式阵列信号仿真系统的抽样时间是T_s，采样率为F_s，水声目标的运动速度为V，实时的目标运动轨迹为n1表示整型变量；

获取水声目标的实时的目标运动轨迹的方法有：第一种方法：基于直线航路、圆周航路或者蛇形追踪航路的运动方式模拟水声目标的实时目标运动轨迹；或者第二种方法：基于B样条拟合的非平稳目标运动轨迹实时仿真方法，获得实时的目标运动轨迹；

第一种方法的过程为：当水声目标的运动方式是直线航路时，设置直线航路的初始位置为(x[0],y[0],z[0])^T，上标T表示转置；运动速度为V，运动的水平方向角为α，垂直仰角为θ，获得如式(1)所示的直线航路的实时的目标运动轨迹

$\stackrel{\&RightArrow;}{P}\left[n1\right]=\left(\begin{array}{c}x\left[n1\right]\\ y\left[n1\right]\\ z\left[n1\right]\end{array}\right)=\left(\begin{array}{c}x\left[0\right]+n1\&CenterDot;{\mathrm{VT}}_s\cos \left(\mathrm{\&theta;}\right)\cos \left(\mathrm{\&alpha;}\right)\\ y\left[0\right]+n1\&CenterDot;{\mathrm{VT}}_s\cos \left(\mathrm{\&theta;}\right)\sin \left(\mathrm{\&alpha;}\right)\\ z\left[0\right]+n1\&CenterDot;{\mathrm{VT}}_s\sin \left(\mathrm{\&theta;}\right)\end{array}\right)$    式(1)

其中，x[n1]表示实时的目标运动轨迹在X轴方向的坐标，y[n1]表示实时的目标运动轨迹在Y轴方向的坐标，z[n1]表示实时的目标运动轨迹在Z轴方向的坐标；

当水声目标的运动方式是圆周航路时，设置圆周航路的初始位置为(x[0],y[0],z[0])^T，速度为V，圆心上标T表示转置；圆周运动半径水声目标的盘旋角速度ω满足式(2)：

$\mathrm{\&omega;}=\frac{V}{r}=\frac{V}{\sqrt{{(X_0-x{\left[0\right]}^2)}^2+{(Y_0-y[0\left]\right)}^2}}$    式(2)

圆周航路的实时的目标运动轨迹如式(3)所示：

$\stackrel{\&RightArrow;}{p}\left[n1\right]=\left(\begin{array}{c}x\left[n1\right]\\ y\left[n1\right]\\ z\left[n1\right]\end{array}\right)=\left(\begin{array}{c}{X}_0+\left(x\right[n1-1]-X_0)\cos \left({\mathrm{\&omega;T}}_s\right)+\left(y\right[n1-1]-Y_0)\sin \left({\mathrm{\&omega;T}}_s\right)\\ Y_0+\left(y\right[n1-1]-Y_0)\cos \left({\mathrm{\&omega;T}}_s\right)-\left(x\right[n1-1]-X_0)\sin \left({\mathrm{\&omega;T}}_s\right)\\ Z_0\end{array}\right)$    式(3)

其中，x[n1-1]表示n1-1时刻的实时的目标运动轨迹在X轴方向的坐标，y[n1-1]表示n1-1时刻的实时的目标运动轨迹在Y轴方向的坐标；

当水声目标的运动方式是蛇形追踪航路时，设置蛇形追踪航路的运动速度为V，轨迹的水平跨度为L₀，拱高H，蛇形追踪航路的实时的目标运动轨迹如式(4)所示：

$\stackrel{\&RightArrow;}{P}\left[n1\right]=\stackrel{\&RightArrow;}{P}[n1-1]+\mathrm{\&Delta;}\stackrel{\&RightArrow;}{P}[n1-1]$    式(4)

其中，表示蛇形追踪航路n1-1时刻的运动坐标，表示蛇形追踪航路当前时刻和前一时刻，即在n1时刻和n1-1时刻的运动坐标的差值，满足式(5)：

$\mathrm{\&Delta;}\stackrel{\&RightArrow;}{P}[n1-1]={\mathrm{VT}}_s\left(\begin{array}{c}\frac{1}{\sqrt{1+{\left(\frac{2\mathrm{\&pi;H}}{L_0}\cos \left(\frac{2\mathrm{\&pi;x}[n1-1]}{L_0}\right)\right)}^2}}\\ \frac{\frac{2\mathrm{\&pi;H}}{L_0}\cos \left(\frac{2\mathrm{\&pi;x}[n1-1]}{L_0}\right)}{\sqrt{1+{\left(\frac{2\mathrm{\&pi;H}}{L_0}\cos \left(\frac{2\mathrm{\&pi;x}[n1-1]}{L_0}\right)\right)}^2}}\end{array}\right)$    式(5)

其中，n1表示整型变量，π表示圆周率；

第二种方法的过程为：获取采样率为f_s的实测水声目标运动轨迹序列m(n)：该运动轨迹序列m(n)包含N个运动轨迹点，待拟合的目标运动轨迹的采样率为F_s；利用基于B样条拟合的非平稳目标运动轨迹实时仿真方法，获得实时的目标运动轨迹；

第二步：获取水声目标辐射噪声信号：设水声目标辐射噪声信号由线谱分量、调制谱分量和连续谱分量组成，其中，线谱的数量为M3，第1—M3中的第q₄条线谱的幅度为第q₄条线谱的频率为和第q₄条线谱的随机相位为设调制谱分量为由M4个具有随机幅度、同样形状、相同重复周期的脉冲性随机过程组成，其中单个脉冲的形状取为高斯型，设置第1—M4中的第q₅个脉冲随机过程的调制深度为第q₅个脉冲随机过程的脉冲宽度为螺旋桨旋转的叶频为f_b，并且连续谱分量利用三参数模型模拟，然后进行AR建模，计算出AR滤波器的分子系数b₀和M5个AR滤波器分母系数中的第q₆个分母系数a_p(q₆)，依据式(6)获得实时的水声目标辐射噪声信号s(n2)为：

   式(6)

其中，x(n2)表示利用高斯白噪声进行AR建模，获得的n2时刻的水声目标辐射噪声信号的连续谱分量，x(n2-q₆)表示n2-q₆时刻的水声目标辐射噪声信号的连续谱分量；q₄、q₅和q₆为标号，q₄取1—M3，q₅取1—M4，q₆取1—M5，T_b＝1/f_b，v(n2)表示高斯白噪声，n2为整型变量，表示第q₅个脉冲随机过程，满足式(7)：

${\mathrm{\&Lambda;}}_{q_5}(T_s\&CenterDot;n2-(q_5-1)\&CenterDot;T_b)=\frac{E_{q_5}}{\sqrt{2\mathrm{\&pi;}}}{e}^{(-\frac{{(T_s\&CenterDot;n2-(q_5-1)\&CenterDot;T_b)}^2}{{2\mathrm{\&sigma;}}_{q_5}})}$    式(7)

设水声目标辐射噪声信号仿真采用用于水声目标连续信号仿真的有限长数据循环拼接方法，获取有限长的实测水声目标数据序列s(n3)，序列长度为N0；利用用于水声目标连续信号仿真的有限长数据循环拼接方法获得实时的水声目标辐射噪声信号，n3为整型变量，n3取1—N0；

第三步：测算等间隔排列的线阵阵形畸变后各个阵元的当前坐标：首先获取等间隔排列的线阵的阵元间距d，阵元个数n10和声纳布放的深度z，从而测算等间隔排列的线阵无阵形畸变时各个阵元的坐标，记无阵形畸变时第1—n10号阵元中的第i₇号阵元的坐标为i₇为阵元标号，i₇取值为1—n10；其中，阵元的标号从坐标系的原点开始，按照X轴方向，顺序标号，从1标至n10号；然后，测算阵形畸变后各个阵元的当前坐标：在等间隔排列的线阵出现阵形畸变时，当阵形畸变是正弦随机扰动时，设置正弦随机扰动的扰动参数，正弦随机扰动的扰动参数包括正弦随机扰动的随机频率f₇、正弦随机扰动的随机幅度amp₇和标号为1的阵元的初始相位当阵形畸变是高斯随机扰动时，设置高斯随机扰动的扰动参数，高斯随机扰动的扰动参数包括高斯随机扰动的均值μ₇和高斯随机扰动的方差从而获得阵形畸变后第1—n10号阵元中的第i₇号阵元的当前坐标

第四步：生成实时的海洋信道脉冲响应序列：获取实测海洋环境参数，包括声速曲线、信号源的水平位置和垂直深度、波束的出射开角、波束传播的最远距离和最大深度、海底系数、海面系数和海洋环境的中心频率，利用基于射线模型的海洋信道仿真方法生成实时的海洋信道脉冲响应序列；

第五步：无阵元随机抖动情况下，生成阵元信号：将实时的水声目标辐射噪声信号和从水声目标的当前位置到达各个阵元当前坐标的实时的海洋信道脉冲响应序列相卷积，生成无阵元随机抖动情况下水声目标辐射噪声信号从水声目标位置传播到各个阵元的阵元信号k7为整型变量；实时的水声目标辐射噪声信号从第二步获得，水声目标的当前位置是从第一步获得的当前时刻的实时的目标运动轨迹，各个阵元当前坐标从第三步获得；

第六步：出现随机抖动时，第i₇号阵元实际生成的阵元信号：第1—n10号阵元中的第i₇号阵元出现随机抖动时，第i₇号阵元实际生成的阵元信号为：

${s^{\&prime;}}_{i_7}\left(k7\right)=s_{i_7}\left(k7\right)(1-a_{i_7})+M_{i_7}$    式(8)

其中，表示第i₇号阵元随机抖动的均值偏移，表示第i₇号阵元随机抖动的幅度衰减；

当第i₇号阵元无随机抖动时，第i₇号阵元实际生成的阵元信号 ${d^{\&prime;}}_{i_7}\left(k7\right)=s_{i_7}\left(k7\right);$

第七步：获取阵元的实时海洋环境噪声信号：假定海洋环境噪声的功率谱为p_e(f)，将p_e(f)通过傅里叶反变换测算出自相关序列，再利用Levinson-Durbin递归算法求解出AR滤波器的分子系数b_e0和M8个AR滤波器分母系数中的第q₈个分母系数a_ep(q₈)，依据式(9)获得第1—n10号阵元中的第i₇号阵元的实时海洋环境噪声信号

${\mathrm{se}}_{i_7}\left(n8\right)=-\underset{q_8=1}{\overset{M8}{\mathrm{\&Sigma;}}}{a}_{\mathrm{ep}}\left(q_8\right)\&CenterDot;{\mathrm{se}}_{i_7}(n8-q_8)+b_{e0}v\left(n8\right)$    式(9)

其中，为n8-q₈时刻的实时海洋环境噪声信号，v(n8)表示第1—n10号阵元中的第i₇号阵元的实时海洋环境噪声信号对应的高斯白噪声，第1—n10号阵元的各个阵元的实时海洋环境噪声信号对应的高斯白噪声相互独立；

第八步：获得被动声纳阵列信号：被动声纳阵列信号的第1—n10号阵元中的第i₇号阵元端接收到的信号是第i7号阵元实际生成的阵元信号和第i7号阵元的实时海洋环境噪声信号的叠加。

2.按照权利要求1所述的开放式阵列信号的仿真方法，其特征在于，所述的第一步中，基于B样条拟合的非平稳目标运动轨迹实时仿真方法，获得实时的目标运动轨迹的过程如下：

步骤101)，获得实测水声目标运动轨迹子序列：化简f_s/F_s为最简分子式，该最简分子式的分子为M，当M大于10时，则M0＝M，当M小于或等于10时，令Q＝2gM，g为从1开始的正整数，且Q为第一次大于10的数时，则M0＝Q；然后将实测水声目标运动轨迹序列m(n)从起始位置按长度M0截取，从而获得L段截取的实测水声目标运动轨迹子序列，最后一段子序列长度为M1，M1小于或等于M0，从第1段到第L-1段中的第i段子序列包含的运动轨迹点为m(M0(i-1)+1),m(M0(i-1)+2),…，m(M0(i-1)+M0)，第L段子序列包含的运动轨迹点为m(M0(L-1)+1),m(M0(L-1)+2),…,m(M0(L-1)+M1)；i表示子序列的段号，i为1至L；0<f_s<F_s；

步骤102)：获取各子序列反求三次B样条曲线的控制顶点：对每段目标运动轨迹子序列反求三次B样条曲线的控制顶点，当第L段子序列包含M0个运动轨迹点时，每段子序列的控制顶点为p_i(1),p_i(2),…,p_i(M0+1),p_i(M0+2)；当第L段子序列包含小于M0个运动轨迹点时，前L-1段子序列的控制顶点为p_i(1),p_i(2),…,p_i(M0+1),p_i(M0+2)，第L段子序列的控制顶点为p_L(1),p_L(2),…,p_L(M1+1),p_L(M1+2)；

步骤103)：获取各子序列在采样率F_s下的目标运动轨迹：利用获得的各段子序列的控制顶点分别作三次均匀B样条拟合，获得从第1段到第L-1段中的第i段子序列在采样率F_s下的目标运动轨迹m′_i(t₁)：

$m_i^{\&prime;}\left(t_1\right)=\frac{1}{6}\left[\begin{array}{cccc}1& t_1& t_1^2& t_1^3\end{array}\right]\left[\begin{array}{cccc}1& 4& 1& 0\\ -3& 0& 3& 0\\ 3& -6& 3& 0\\ -1& 3& -3& 1\end{array}\right]\left[\begin{array}{c}{p}_i\left(j_1\right)\\ p_i(j_1+1)\\ p_i(j_1+2)\\ p_i(j_1+3)\end{array}\right]$    式(10)

式(10)中，j₁＝1、2、…、M0-3,t₁＝k×f_s/F_s,t₁≤M0，k表示整型变量，j₁为大于t₁的最小整数；

当第L段子序列包含M0个运动轨迹点时，第L段子序列在采样率F_s下的目标运动轨迹m′_L,M0(t₂)：

$m_{L,M0}^{\&prime;}\left(t_2\right)=\frac{1}{6}\left[\begin{array}{cccc}1& t_2& t_2^2& t_2^3\end{array}\right]\left[\begin{array}{cccc}1& 4& 1& 0\\ -3& 0& 3& 0\\ 3& -6& 3& 0\\ -1& 3& -3& 1\end{array}\right]\left[\begin{array}{c}{p}_L\left(j_2\right)\\ p_L(j_2+1)\\ p_L(j_2+2)\\ p_L(j_2+3)\end{array}\right]$    式(11)

式(11)中，j₂＝1、2、…、M0-3,t₂＝k×f_s/F_s,t₂≤M0，k表示整型变量，j₂为大于t₂的最小整数；

当第L段子序列包含小于M0个运动轨迹点时，第L段子序列在采样率F_s下的目标运动轨迹m′_L,M1(t₃)：

$m_{L,M1}^{\&prime;}\left(t_3\right)=\frac{1}{6}\left[\begin{array}{cccc}1& t_3& t_3^2& t_3^3\end{array}\right]\left[\begin{array}{cccc}1& 4& 1& 0\\ -3& 0& 3& 0\\ 3& -6& 3& 0\\ -1& 3& -3& 1\end{array}\right]\left[\begin{array}{c}{p}_L\left(j_3\right)\\ p_L(j_3+1)\\ p_L(j_3+2)\\ p_L(j_3+3)\end{array}\right]$ 式(12)

式(12)中，j₃＝1、2、…、M1-3,t₃＝k×f_s/F_s,t₃≤M1，k表示整型变量，j₃为大于t₃的最小整数；

步骤104)：获得实时的目标运动轨迹：将第1段到第L段的各段子序列在采样率F_s下的目标运动轨迹按段号由小到大合并，获得采样率F_s下的非平稳目标运动轨迹，从而获得实时的目标运动轨迹。

3.按照权利要求2所述的开放式阵列信号的仿真方法，其特征在于，所述的步骤102)包括以下步骤：假设目标运动轨迹子序列的长度为M2，第1段到第L-1段中的第i段子序列对应的M2＝M0，当最后一段子序列长度为M1等于M0时，M2＝M0，当最后一段子序列长度M1小于M0时，M2＝M1，利用式(13)求解子序列的控制顶点p(1),p(2),…,p(M2+1),p(M2+2)：

(p(q),p(q+1),p(q+2))/6＝m(q)   式(13)

其中，q表示整型变量，q＝1,2,…,M2；m(q)表示子序列包含的目标运动轨迹点；式(13)方程组有M2个方程，M2+2个未知数，根据式(14)设定两个边界条件：

$\left\{\begin{array}{c}p\left(1\right)=p\left(2\right)\\ p(M2+1)=p(M2+2)\end{array}\right.$    式(14)

将式(13)和式(14)相结合，得到如式(15)所示的M2+2个线性方程：

$\left[\begin{array}{cccccc}-1& 1& & & & \\ 1& 4& 1& & & 0\\ & 1& 4& 1& & \\ & \&CenterDot;\&CenterDot;\&CenterDot;& & & & \\ & 0& & 1& 4& 1\\ & & & & 1& -1\end{array}\right]\left[\begin{array}{c}p\left(1\right)\\ p\left(2\right)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ p\left(M2\right)\\ p(M2+1)\\ p(M2+2)\end{array}\right]=6\left[\begin{array}{c}0\\ m\left(1\right)\\ m\left(2\right)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ m\left(M2\right)\\ 0\end{array}\right]$   式(15)

每段子序列通过上述方法得到A，从而得到该段子序列中A的伴随矩阵A^*，1≤q₁≤M2+2，1≤q₂≤M2+2，表示矩阵A的第q₁行第q₂列元的代数余子式，A^*满足式(16)和式(17)：

$\left\{\begin{array}{c}{A}_{1M2+2}={(-1)}^{M2+3}\\ A_{2M2+2}={(-1)}^{M2+5}\\ A_{q_{3}M2+2}=-4A_{q_3-1M2+2}-A_{q_3-2M2+2}\end{array}\right.$    式(16)

|A|＝-A_M2+2M2+2+A_M2+1M2+2   式(17)

其中，A_1M2+2表示矩阵A的第1行第M2+2列元的代数余子式，A_2M2+2表示矩阵A的第2行第M2+2列元的代数余子式，表示矩阵A的第q₃行第M2+2列元的代数余子式，表示矩阵A的第q₃-1行第M2+2列元的代数余子式，表示矩阵A的第q₃-2行第M2+2列元的代数余子式，A_M2+2M-2表示矩阵A的第M2+2行第M2+2列元的代数余子式，A_M2+1M2+2表示矩阵A的第M2+1行第M2+2列元的代数余子式；3≤q₃≤M2+2；

利用式(16)和式(17)得到A^*的第M2+2行元素和|A|，从而依据式(18)得到p(M2+2)：

$p(M2+2)=\left(6\underset{q=1}{\overset{M2}{\mathrm{\&Sigma;}}}{A}_{q+1M2+2}m\left(q\right)\right)/\left|A\right|$    式(18)

其中，A_q+1M2+2表示矩阵A的第q+1行第M2+2列元的代数余子式；

最后递推式(13)，得到M2+2个控制顶点p(1),p(2),…,p(M2+1),p(M2+2)，从而获得各段子序列的控制顶点。

4.按照权利要求1所述的开放式阵列信号的仿真方法，其特征在于，所述的第二步中，三参数模型是用户通过设置三个参数f₀，f_m和λ来精确定义平稳连续谱p_c(f)的形状，p_c(f)满足：

$p_c\left(f\right)\frac{{\mathrm{\&sigma;}}_1^2}{2\mathrm{\&pi;}}[\frac{f_m+\mathrm{\&lambda;}(f+f_0)}{f_m^2+{(f+f_0)}^2}+\frac{f_m-\mathrm{\&lambda;}(f-f_0)}{f_m^2+{(f-f_0)}^2}]$    式(19)

其中，表示用来控制平稳连续谱信号的能量，f表示功率谱的线频率；

然后，将p_c(f)通过傅里叶反变换计算出自相关序列，再利用Levinson-Durbin递归算法求解AR滤波器的分子系数b₀和M5个AR滤波器分母系数中的第q₆个分母系数a_p(q₆)。

5.按照权利要求1所述的开放式阵列信号的仿真方法，其特征在于，所述的第二步中，用于水声目标连续信号仿真的有限长数据循环拼接方法包括以下步骤：

首先：在实测数据段s(n3)首尾两端各找到一段能量最低的区域：从实测数据s(n3)的首尾部开始，在N_s点的搜索范围内以移窗的方式分别从i₁点开始搜索首部能量最低的区域，从i₂点开始搜索尾部能量最低的区域，窗长WL，记首尾部每段窗内的数据段的能量记为ph(i₁)和ph(i₂)：

$\mathrm{ph}\left(i_1\right)=\underset{n4=0}{\overset{\mathrm{WL}-1}{\mathrm{\&Sigma;}}}{s}^2(n4+i_1)$    式(20)

$\mathrm{ph}\left(i_2\right)=\underset{n5=0}{\overset{\mathrm{WL}-1}{\mathrm{\&Sigma;}}}{s}^2(n5+i_2)$    式(21)

其中，i₁表示搜索首部能量最低区域过程中，首部每次搜索的起点；i₂表示搜索尾部能量最低区域过程中，尾部每次搜索的起点；i₁∈[1,N_s-1]，i₂∈[N0-N_s,N0-WL]，WL表示每次搜索的长度；记窗头位置移到N₁时找到数据段首部能量最低的区域，记首部能量最低区域为W₁；记窗头位置移到N₂时找到数据段尾部能量最低的区域，记尾部能量最低区域为W₂；N0表示实测水声目标数据序列s(n3)的序列长度，s²(n₄+i₁)表示n4+i₁点对应的实测水声目标辐射噪声值的平方，s²(n₅+i₂)表示n5+i₂点对应的实测水声目标辐射噪声值的平方；

其次：在首尾两段能量最低区域中找出频谱特征相似的两段区域：以移窗的方式在区域W₁和W₂中分别从i₃点开始搜索首部频谱特性相似的的数据段W₃，从i₄点开始搜索尾部频谱特性相似的数据段W₄，记窗长为WL₁，且WL₁<WL，最低能量区域W₁和W₂内每段窗长的数据段的频谱分别记为和

${\mathrm{Sh}}_{i_3}\left(e^{\mathrm{j\&omega;}}\right)=\underset{n5=0}{\overset{{\mathrm{WL}}_1-1}{\mathrm{\&Sigma;}}}s(n5+i_3)\&CenterDot;e^{-\mathrm{j\&omega;n}5}$    式(22)

${\mathrm{St}}_{i_4}\left(e^{\mathrm{j\&omega;}}\right)=\underset{n6=0}{\overset{{\mathrm{WL}}_1-1}{\mathrm{\&Sigma;}}}s(n6+i_4)\&CenterDot;e^{-\mathrm{j\&omega;n}6}$    式(23)

其中，i₃表示搜索首部最低能量区域中频谱特征相似的过程中，首部每次搜索的起点，i₄表示搜索尾部最低能量区域中频谱特征相似的过程中，尾部每次搜索的起点，i₃∈[N₁,N₁+WL-WL₁-1]，i₄∈[N₂,N₂+WL-WL₁-1]，j表示虚数单位，ω表示频谱的角频率，ω取值范围为-π—π；s(n5+i₃)表示n5+i₃点对应的实测水声目标辐射噪声值，s(n6+i₄)表示n6+i₄点对应的实测水声目标辐射噪声值；设W₁和W₂区域内窗头位置分别移到N₃和N₄，和最为接近时，则找到最低能量区域内频谱最为相似的两个数据段W₃和W₄；N₃表示和最为接近时，首部搜索的起点；N₄表示和最为接近时，尾部搜索的起点；

再次：在频谱特征相似的两段区域里找出导数同号的两组过零点：分别在区域W₃、W₄中找到符合条件的两组过零点s(N₅),s(N₅+1)和s(N₆),s(N₆+1)，满足：

s(N₅)·s(N₅+1)＜0,s(N₆)·s(N₆+1)＜0   式(24)

(s(N₅+1)-s(N₅))·(s(N₆+1)-s(N₆))＞0   式(25)

其中，N₅、N₅+1、N₆和N₆+1分别表示整型变量，s(N₅)表示N₅对应的实测水声目标辐射噪声值，s(N₅+1)表示N₅+1对应的实测水声目标辐射噪声值，s(N₆)表示N₆对应的实测水声目标辐射噪声值，s(N₆+1)表示N₆+1对应的实测水声目标辐射噪声值；

然后：确定循环数据段：实测信号中s(N₅)至s(N₆+1)之间的信号段，记为s_r(n7)，即满足局部相似条件的循环数据段，n7为整型变量；

最后：将循环数据段s_r(n7)不断重复拼接，获得实时的连续水声目标信号数据。

6.按照权利要求1所述的开放式阵列信号的仿真方法，其特征在于，所述的第三步中，依据式(26)得到无阵形畸变时第1—n10号阵元中的第i₇号阵元的坐标 $(x_{i_7},y_{i_7},z_{i_7}):$

$\left\{\begin{array}{c}{x}_{i_7}=(i_7-1)d\\ y_{i_7}=0\\ z_{i_7}=z\end{array}\right.$    式(26)

其中，d表示等间隔排列的线阵的阵元间距，i₇表示阵元标号，z表示声纳布放的深度；

依据式(27)得到正弦随机扰动时第1—n10号阵元中的第i₇号阵元的当前坐标 $(x_{i_7},y_{i_7},z_{i_7}):$

   式(27)

依据式(28)得到高斯随机扰动时第1—n10号阵元中的第i₇号阵元的当前坐标

$\left\{\begin{array}{c}{x}_{i_7}^{\&prime;}(i_7-1)d\\ f\left(y_{i_7}^{\&prime;}\right)=\frac{1}{\sqrt{2\mathrm{\&pi;}{\mathrm{\&sigma;}}_7^2}}{e}^{(-\frac{{(y_{i_7}^{\&prime;}-{\mathrm{\&mu;}}_7)}^2}{2{\mathrm{\&sigma;}}_7^2})}\\ z_{i_7}^{\&prime;}=z\end{array}\right.$    式(28)

其中，表示第i₇号阵元的坐标的概率密度函数，满足高斯分布，e表示数学常数，为2.7182818。

7.按照权利要求1所述的开放式阵列信号的仿真方法，其特征在于，所述的第四步中，利用基于射线模型的海洋信道仿真方法生成实时的海洋信道脉冲响应序列的步骤包括：

首先，基于实测海洋环境参数，在主机端利用射线模型的Bellhop软件包生成实测海洋环境下的海洋信道脉冲响应序列库，海洋信道脉冲响应序列库是水声目标辐射噪声信号从不同的水声目标位置传播到不同的阵元的本征声线，包括本征声线的数目，每条本征声线的幅度衰减、相移、传播时延、反射次数、折射次数和本征声线到达阵元时的入射角；

其次：获得实时的海洋信道脉冲响应序列：主机从DSP板读取水声目标实时的目标运动轨迹，每次读取水声目标实时的目标运动轨迹后，主机根据水声目标的当前坐标和各个阵元的当前坐标获取水声目标辐射噪声信号到达各个阵元的本征声线，测算当前的海洋信道脉冲响应序列，然后更新水声目标辐射噪声信号到达各个阵元的海洋信道脉冲响应序列为当前的海洋信道脉冲响应序列，再通过总线将当前的海洋信道脉冲响应序列加载进DSP板；通过主机每隔L1个采样间隔从DSP板读取一次水声目标运动坐标，在DSP板上获得实时的海洋信道脉冲响应序列；L1的取值范围为128—512之间的整数。

8.按照权利要求7所述的开放式阵列信号的仿真方法，其特征在于，所述的第四步中，根据水声目标的当前坐标和各个阵元的坐标，主机获取水声目标辐射噪声信号到达各个阵元的本征声线，当从不同水声目标位置传播到不同阵元的本征声线小于m1条时，选出所有本征声线作为主机获取的水声目标辐射噪声信号到达各个阵元的本征声线，当本征声线大于或等于m1条时，将所有的本征声线按照幅度衰减降序排列，筛选出前m1条本征声线作为主机获取的水声目标辐射噪声信号到达各个阵元的本征声线，m1为10—20之间的整数。

9.按照权利要求8所述的一种开放式阵列信号仿真系统，其特征在于，所述的第四步中，测算当前的海洋信道脉冲响应序列的过程是：首先选取1号阵元为参考阵元，假设水声目标在运动过程中以球面波的形式不断的向周围发射辐射噪声信号，则水声目标辐射噪声信号到达各个阵元的传播时延为：

$t_{i_7}=({\mathrm{positon}}_{{\mathrm{array}}_{i7}}-{\mathrm{positon}}_{t\arg \mathrm{et}})/c$    式(29)

其中，为第i₇号阵元的坐标，positon_target为水声目标的当前坐标，c为海洋中声波传播速度；

各个阵元的信号到达时间相对于参考阵元的信号到达时间的时延差为：

$t_{i_7}-t_{i_0}=({\mathrm{positon}}_{{\mathrm{array}}_{i7}}-{\mathrm{positon}}_{t\arg \mathrm{et}})/c-({\mathrm{positon}}_{{\mathrm{array}}_{i0}}-{\mathrm{positon}}_{t\arg \mathrm{et}})/c$    式(30)

其中，下标i₀表示参考阵元的序号，表示参考阵元的坐标；

水声目标辐射噪声信号从当前坐标到达各个阵元的海洋信道脉冲响应序列 $h_{i_7}\left(n8\right)$ 为：

   式(31)

其中，i₇表示第i₇号阵元，k1表示第k1条本征声线，为第i₇号阵元对应的本征声线数目，表示第i₇号阵元的第k1条本征声线的幅度衰减，示第i₇号阵元的第k1条本征声线的相移，表示第i₇号阵元的第k1条本征声线的传播时延；e表示数学常数，为2.7182818，j表示虚数单位，n8表示整型变量， $\mathrm{\&sigma;}(n8-{\mathrm{\&tau;}}_{i_{7}k1}-(t_{i_7}-t_{i_0}))$ 表示单位采样信号， $\mathrm{\&sigma;}(n8-{\mathrm{\&tau;}}_{i_{7}k1}-(t_{i_7}-t_{i_0}))$ 满足式(32)：

$\mathrm{\&sigma;}(n8-{\mathrm{\&tau;}}_{i_{7}k1}-(t_{i_7}-t_{i_0}))=\left\{\begin{array}{cc}1& n8={\mathrm{\&tau;}}_{i_{7}k1}-(t_{i_7}-t_{i_0})\\ 0& n8\&NotEqual;{\mathrm{\&tau;}}_{i_{7}k1}-(t_{i_7}-t_{i_0})\end{array}\right.$    式(32)。

10.按照权利要求1所述的开放式阵列信号的仿真方法，其特征在于，所述的第六步中，第i₇号阵元随机抖动的幅度衰减服从均匀分布：

$f\left(a_{i_7}\right)=\frac{1}{B_7-A_7}$    式(33)

其中，表示幅度衰减的概率密度函数，A₇和B₇是均匀分布的两个参数， $A_7\&le;a_{i_7}\&le;B_7,$ 均匀分布的均值为 $\frac{A_7+B_7}{2},$ 均匀分布的方差为 $\frac{{(B_7-A_7)}^2}{12}.$