一种具有空间相关性的水声阵列海洋环境噪声仿真方法与流程

本发明涉及水声阵列信号仿真技术领域，特别是涉及一种具有空间相关性的水声阵列海洋环境噪声仿真方法。

背景技术：

声纳等水声系统利用水听器阵列接收水中声信号，实现水中目标检测、定位、识别等功能。随着现代计算机性能的不断提高，越来越多的水声技术研究人员利用计算机仿真方式生成水声阵列接收信号，对水声信号处理算法以及系统进行测试评估，以降低海上试验成本，提高测试评估效率。为了达到这一目的，水声阵列接收信号仿真系统需要尽可能地还原海上试验情景，生成的仿真信号需要具有与海上实测信号良好的吻合度。水声阵列接收信号是随机信号，因此仿真生成的水声阵列接收信号的统计特性需要尽可能地与海上实测信号的统计特性相吻合。其中，水声阵列各阵元接收的海洋环境噪声的时间和空间统计特性就是水声阵列信号仿真的重要部分。

为提高信号处理的增益，现有的声纳系统通常利用水听器阵列接收水下声信号，并采用波束形成方法进行目标方位测量和目标跟踪。在波束形成方法中，一般假设各阵元接收的背景噪声服从高斯分布，且是不相关的，这样波束形成可获得最大的空间增益[1]。因此，大多数现有水声阵列信号仿真系统忽略了各阵元背景噪声之间的相关性[1][2][3]。但是，作为水声阵列主要背景噪声源之一的海洋环境噪声具有一定的空间相关性[4]，这使得阵列信号处理算法的仿真测试和实际使用之间存在性能上的差异。

海洋环境噪声是水声阵列背景噪声的重要组成部分，是一种宽带噪声。一方面，由于海洋环境噪声自身的多变性以及水声阵列工作环境的复杂性，建立起一套易于信号仿真生成使用的宽带噪声理论模型和应用模型难度很大。另一方面，受到计算能力的限制，对于海洋环境噪声仿真，在满足符合时频特征要求的情况下，很难兼顾阵元间的海洋环境噪声空间相关性。因此，建立既满足符合时频特征要求，又符合空间相关特性要求的宽带海洋环境噪声仿真生成方法，对于提升水声阵列信号仿真的逼真性和有效性具有重要意义。

本发明构建了水声阵列宽带海洋环境噪声阵元间空间相关序列，根据空间相关序列的特性对离散频率进行了频带分解，并构造了子频带内恒定空间相关系数矩阵。利用分频带方法生成符合各子频带相关特性的带限白噪声，再合成空间相关宽带高斯白噪声；然后将相关宽带高斯白噪声作为ar(autoregressive，自回归)滤波器激励信号，采用ar建模方式生成宽带有色海洋环境噪声生成，利用ar滤波器的线性特性，将高斯白噪声激励间的相关性线性映射为不同阵元海洋环境噪声间的相关性，实现具有空间相关性的水声阵列海洋环境噪声时域序列仿真生成。

技术实现要素：

针对现有水声阵列信号仿真技术中存在的问题，本发明提供一种具有空间相关性的水声阵列海洋环境噪声仿真方法，该方法先求解分频带空间相关系数矩阵，生成符合各子频带相关特性的带限白噪声，然后将各子频带带限白噪声合成得到相关宽带高斯白噪声，并激励ar滤波器，最终仿真生成具有空间相关特性的水声阵列宽带海洋环境噪声，为达此目的，本发明提供一种具有空间相关性的水声阵列海洋环境噪声仿真方法，包括如下步骤：

(1)设定水声阵列海洋环境噪声仿真参数；

(2)根据水声阵列参数求解m个阵元两两之间随频率变化的海洋环境噪声空间相关函数序列quv(fp)，1≤u<v≤m，m为阵元数，p＝0,1,…,nf-1，nf为离散频率总点数；并依照水声阵列中间距最小的两个阵元对应的空间相关函数序列(fp)，对离散频率序列fp进行频带分解，得到h个子频带bj，j＝1,2…,h；

(3)根据空间相关函数序列quv(fp)，计算各子频带内相关系数的参考值

(4)生成m个长度为n，且互不相关的高斯白噪声序列xl(n)，l＝1,2,…,m，n＝0,1,2,…,n-1；

(5)依照步骤2中的频带分解结果，对步骤4中m个互不相关的宽带高斯白噪声序列xl(n)进行带通滤波，生成h组带限高斯白噪声序列，每组中的m个带限高斯白噪声序列互不相关；

(6)依照步骤3中得到的各子频带内空间相关系数参考值，生成h组相关带限高斯白噪声序列，第j组中的m个带限高斯白噪声序列分别对应于m个阵元在bj频带内的带限相关高斯白噪声序列，且两两之间具有特定的相关系数；

(7)将与每个阵元对应的h个相关带限高斯白噪声序列相加，生成与m个阵元对应的相关宽带高斯白噪声序列yl(n)；

(8)根据设定的海洋环境噪声功率谱参数，求解海洋环境噪声ar滤波器系数；

(9)用步骤7中得到的与各阵元对应的相关宽带高斯白噪声序列yl(n)作为输入，激励步骤8中得到的ar滤波器，生成具有空间相关性的水声阵列海洋环境噪声。

作为本发明进一步改进，步骤1设定如下参数：

声速c，采样频率fs，离散频率点数nf，频率间隔δf＝fs/nf，水声阵列阵元数量m，依照附图2所示的坐标系设定各阵元坐标(xl，yl，zl)，阵元坐标用矩阵p表示，向量pl表示第l个阵元的坐标，取第一个阵元位置为坐标原点，即矩阵p可用式表示；设定空间相关函数划分精度δq，0＜δq≤1；

设定海洋环境噪声功率谱参数σ、fm、fm0和λ。

作为本发明进一步改进，步骤2具体包括如下步骤：

(2.1)依据式设定的水声阵列阵元坐标，代入～式计算水声阵列阵元两两之间的间距序列suv，阵元连线俯仰角序列γuv以及阵元连线的方位角序列ζuv，再代入式计算对应的海洋环境噪声空间相关函数序列quv(fp)，其中zuv＝zv-zu，xuv＝xv-xu；

式中表示方向概率密度函数，i为虚数单位，θ表示俯仰角，表示方位角；

(2.2)依据得到的海洋环境噪声空间相关函数序列quv(fp)，选取suv的最小值对应的空间相关函数序列(fp)作为参考序列进行频带划分；

首先，依照空间相关函数序列(fp)的单调性，将离散频率序列fp划分为k个区间；前k-1个区间内，各离散频率对应的相关函数值单调变化，且最大值和最小值的差不小于δq；第k个区间内，各离散频率对应的相关函数最大值和最小值的差小于δq；第1个离散频率区间表示为第2个离散频率区间表示为第k个离散频率区间表示为第k个离散频率区间表示为且有sk＝nf-1，k＝1,2,…,k；

然后，在每个频率区间内按照式所示的离散频率序列下标间隔δnk对各离散频率区间再进行等间隔划分，每个频率区间划分为nk个子频带，nk取值如式所示，总共划分为h个子频带：

式和式中round表示四舍五入取整，qkmax为第k个区间内相关函数序列的最大值，qkmin为第k个子区间内相关函数序列的最小值，h和k的满足式的关系；

作为本发明进一步改进，步骤3具体包括如下步骤：

(3.1)根据步骤2中得到的频带划分结果，对每个quv(fp)进行频带划分，得到各子频带内的空间相关函数序列(fq)，其中

(3.2)根据式求解空间各子频带内的相关系数参考值

作为本发明进一步改进，步骤4具体包括如下步骤：

(4.1)设定高斯白噪声数学期望μ，方差σg²；

(4.2)采用box-muller算法生成m个点数为n的宽带高斯白噪声序列xl(n)，l＝1,2,…,m，n＝0,1,2,…,n-1，并组成矩阵xn×m＝[x1,x2,…,xm]，该矩阵的每一列为n点的宽带高斯白噪声序列构成的列向量，矩阵各列之间不相关。

作为本发明进一步改进，步骤5具体包括如下步骤：

(5.1)根据步骤2的频带划分结果，采用窗函数法设计fir滤波器，生成h个fir带通滤波器系数序列hj，hj的通带为子频带bj；

(5.2)将步骤4中生成的矩阵xn×m的各列分别与带通滤波器系数序列hj进行卷积运算，如式所示，得到h个带限高斯白噪声矩阵每个带限高斯白噪声矩阵的各列之间是不相关的；

式中xl^j表示矩阵的第l列，xl表示矩阵xn×m的第l列，表示卷积运算。

作为本发明进一步改进，步骤6具体包括如下步骤：

(6.1)根据步骤3中得到的各子频带内相关系数参考值构成如式所示的子频带内空间相关系数矩阵

(6.2)对空间相关系数矩阵进行cholesky分解，将空间相关系数矩阵分解为上三角矩阵的转置与其自身的乘积，如式所示；

(6.3)将不相关带限高斯白噪声矩阵与上三角矩阵相乘，如式所示，得到相关带限高斯白噪声矩阵的第i1列和i2列之间的相关系数等于空间相关系数矩阵的第i1行、i2列的元素

作为本发明进一步改进，步骤7包括如下步骤：

对步骤6得到的h个相关带限高斯白噪声矩阵依照式进行相加，得到相关宽带高斯白噪声矩阵yn×m；

作为本发明进一步改进，步骤8包括如下步骤：

(8.1)采用三参数模型法模拟海洋环境噪声的功率谱序列e(fp)，如所示；

(8.2)建立qn阶ar滤波器的yule-walker方程，如所示；

其中a(ln)和b0为滤波器系数，δ(kn)为冲击函数；rx(kn)为e(fp)对应的自相关函数序列，ln＝1,2,…,qn，kn＝0,1,…,n-1；

(8.3)采用levison-durbin算法求解式，得到qn阶ar滤波器系数a(ln)和b0。

作为本发明进一步改进，步骤9包括如下步骤：

利用步骤7得到的相关宽带高斯白噪声yn×m激励步骤8求得的滤波器系数a(ln)和b0，将yn×m中的每一列yl(n)分别代入式，得到输出序列zl(n)；

zl(n)组成的矩阵zn×m＝[z1(n),z2(n),…,zm(n)]即为具有空间相关性的宽带海洋环境噪声，zn×m的每一列对应为水声阵列的一个阵元接收的海洋环境噪声序列，各列间具有指定的相关特性。

本发明公开的一种具有空间相关性的宽带海洋环境噪声仿真生成的方法，首先计算水声阵列宽带海洋环境噪声阵元间空间相关序列，根据空间相关序列的特性对离散频率进行了频带分解，生成子频带内恒定空间相关系数矩阵。利用分频带方法生成符合各子频带相关特性的带限白噪声，再合成相关宽带高斯白噪声。宽带有色海洋环境噪声生成采用ar建模方式，将相关宽带高斯白噪声作为ar滤波器激励信号，利用ar滤波器的线性特性，将高斯白噪声激励间的相关性线性映射为海洋环境噪声间的相关性，实现具有空间相关性的水声阵列海洋环境噪声时域序列生成。

与现有技术相比，本发明公开的方法具有以下优点：本发明公开的方法将高斯白噪声激励间的相关性线性映射为海洋环境噪声间的相关性，使得海洋环境噪声空间相关模型清晰，易于仿真系统编程实现；采用频带分解与合成的方法，降低计算量，应用简单直接；本发明公开方法生成的海洋环境噪声通道间具有指定的相关性，更贴近真实的水声阵列海洋环境噪声。

附图说明

图1为本发明方法的实施流程图；

图2为本发明方法建立的空间坐标系；

图3为实施例中的海洋环境噪声空间相关函数序列；

图4为实施例中的空间相关系数参考值；

图5为实施例中的海洋环境噪声功率谱；

图6至图10分别为实施例中仿真得到的水声阵列阵元1至阵元5空间相关宽带海洋环境噪声功率谱；

图11为实施例中仿真得到的空间相关系数与理论值对比图；

图12为实施例中仿真得到的空间相关系数多次平均值与理论值对比图。

具体实施方式

下面结合附图与具体实施方式对本发明作进一步详细描述：

本发明提供一种具有空间相关性的水声阵列海洋环境噪声仿真方法，该方法先求解分频带空间相关系数矩阵，生成符合各子频带相关特性的带限白噪声，然后将各子频带带限白噪声合成得到相关宽带高斯白噪声，并激励ar滤波器，最终仿真生成具有空间相关特性的水声阵列宽带海洋环境噪声。

实施例1：

下面结合附图和具体实施方式，进一步阐明本发明。

一种具有空间相关性的水声阵列宽带海洋环境噪声仿真生成方法，如图1所示，包括如下步骤：

步骤1，设定仿真参数

设定声速c＝1500m/s，采样频率fs＝15khz，离散频率点数nf＝150000，频率间隔δf＝fs/nf＝0.1hz，水声阵列阵元数m＝5，δq取值为0.1。本实施例设定水声阵列为等间距线列阵，阵元坐标矩阵p表示如下：

设定海洋环境噪声功率谱参数fm＝200，f0＝500，λ＝0，σ＝0.01。

步骤2，根据水声阵列阵元坐标，求解m个阵元两两之间随频率变化的海洋环境噪声空间相关函数序列quv(fp)，并选取序列suv最小值对应的空间相关函数序列(fp)作为参考序列，进行频带划分，包括如下步骤：

(2.1)依据设定的水声阵列阵元坐标p，代入(2)式～(4)3计算水声阵列阵元间距序列suv，阵元连线俯仰角序列γuv以及阵元连线的方位角序列ζuv，结果如表1所示。

表1阵元间距、连线俯仰角、连线方位角求解结果

表1结果代入(5)式求得阵元间空间相关函数序列quv(fp)。本实施例中，概率密度函数如(19)式所示。

其中π/4≤θ1≤π/2，u1>0，α1以及a1满足(20)式，i0(·)是零阶第一类贝塞尔函数。

本例中，取θ1＝π/3,u1＝100，α1＝10。联合(2)～(4)、(19)～(20)式求解得到水声阵列空间相关函数序列quv(fp)，组如图3所示。

(2.2)依据得到的海洋环境噪声空间相关函数序列quv(fp)，选取suv的最小值对应的空间相关函数序列(fp)作为参考序列进行频带划分。首先，依照空间相关函数序列(fp)的单调性，将离散频率序列fp划分为k个区间；前k-1个区间内，各离散频率对应的相关函数值单调变化，且最大值和最小值的差不小于δq；第k个区间内，各离散频率对应的相关函数最大值和最小值的差小于δq。第1个离散频率区间表示为第2个离散频率区间表示为第k个离散频率区间表示为第k个离散频率区间表示为且有sk＝nf-1，k＝1,2,…,k。本实施例中求得的k＝3，h＝15，频带划分结果如表2所示。

表2频带划分结果

步骤3，求解各阵元间各子频带内相关系数参考值，包括如下步骤：

(3.1)根据表2所示的频带划分结果，对每个quv(fp)进行频带划分，得到子得到各子频带内的空间相关函数序列(fq)。

(3.2)根据(10)式求解空间相关函数序列各子频带内的相关系数参考值本实施例求得的空间相关系数参考值如表3所示，作图如图4所示。

表3空间相关系数参考值

步骤4包括如下步骤：

(4.1)设定高斯白噪声数学期望μ＝0，方差σg²＝1.

(4.2)采用box-muller算法，模拟生成5个点数为15000的宽带高斯白噪声序列xl(n)，l＝1,2,…,5，n＝0,1,2,…,14999，组成矩阵x15000×5＝[x1,x2,…,xm]，该矩阵的每一列为15000点的宽带高斯白噪声序列构成的列向量，矩阵各列之间不相关。

步骤5包括如下步骤：

(5.1)根据表2的频带划分结果，采用函数法设计fir滤波器，窗函数选用汉明窗，每个滤波器为512阶，生成15个fir带通滤波器hj。

(5.2)将步骤(4)中生成的矩阵x15000×5的各列分别与带通滤波器系数序列hj进行卷积运算，得到15个带限高斯白噪声矩阵每个带限高斯白噪声矩阵的各列之间是不相关的。

步骤6，生成多通道相关带限高斯白噪声激励信号，包括如下步骤：

(6.1)根据表3各子频带内相关系数参考值生成子频带内空间相关系数矩阵例如本实施例中第3个子频带内空间相关系数矩阵如(21)所示。

(6.2)对空间相关系数矩阵进行cholesky分解，将空间相关系数矩阵分解为上三角矩阵的转置与其自身的乘积，例如本实施例中第3个子频带内上三角矩阵如(22)所示。

(6.3)将不相关带限高斯白噪声矩阵与上三角矩阵相乘，得到相关带限高斯白噪声矩阵

步骤7包括如下步骤：

对步骤6得到的15个相关带限高斯白噪声矩阵依照(15)式进行相加，得到相关宽带高斯白噪声矩阵y15000×5。

步骤8包括如下步骤：

(8.1)采用三参数模型法模拟海洋环境噪声的功率谱序列e(fp)，功率谱e(fp)如图5所示。

(8.2)建立qn阶ar滤波器的yule-walker方程，rx(kn)为e(fp)对应的自相关函数序列，可利用功率谱e(fp)通过维纳-辛钦定理求得，ln＝1,2,…,qn，kn＝0,1,…,15000-1；本实施例中，qn＝8。

(8.3)采用levison-durbin算法求解yule-walker方程，得到qn阶ar滤波器系数如表4所示。

表4.ar模型系数

步骤9包括如下步骤：

利用步骤7得到的相关宽带高斯白噪声激励y15000×5以及表4的系数a(ln)，ln＝1,2,…,8和b0，将y15000×5中的每一列yl(n)分别代入(18)式，得到输出序列zl(n)。

本例仿真生成的具有空间相关特性的海洋环境噪声序列长度为1秒，对5个阵元的功率谱进行多次平均，结果如图6至图10所示。

为验证水声阵列宽带海洋环境噪声的相关特性，对海洋环境噪声依据表2频带划分进行分频带相关系数求解。仿真共进行10次，仿真结果如图11所示，曲线为理论的相关系数，散点为每次仿真的空间相关系数求解结果，可以看出每次仿真结果在理论值附近波动。对10次仿真的每个子频带内的相关系数进行平均，对比理论的相关系数与仿真的平均海洋环境噪声相关系数如图12所示，可以看出二者之间的相关性一致。

仿真结果表明，采用以上方法能够有效地实现具有相关性的水声阵列宽带海洋环境噪声仿真生成，当空间相关系数仿真精度取为0.1时，仿真生成的海洋环境噪声的阵元间相关特性以及频谱特性与理论值一致。

以上所述，仅是本发明的较佳实施例而已，并非是对本发明作任何其他形式的限制，而依据本发明的技术实质所作的任何修改或等同变化，仍属于本发明所要求保护的范围。