CN104280737B - 加权宽带时反算子分解声成像方法 - Google Patents

Abstract

一种加权宽带时反算子分解声成像方法：将收发合置阵接收到的信号进行快速傅里叶变换，经过窄带滤波器划分为各个子带，将各个子带对应的时反算子奇异值分解，将最大特征值对应的特征向量与基于声传播模型的传递向量相关后获得模糊度函数，提取各子带模糊度函数的最大值及所对应的二维声场空间，同时以最大特征值设计加权系数，加权模糊度函数，并将所有子带的结果相干累加后呈现为三维图像。本发明的优点在于将目标的散射特性通过加权的方式应用于宽带时反算子分解成像方法，充分利用了宽带信息，真实地显示扩展性目标的空间特征，并通过实验验证了该发明方案的有效性和可靠性。

Description

加权宽带时反算子分解声成像方法

技术领域

本发明属于海洋资源开发和利用的声探测领域，具体的说是一种采用阵列对点目标和扩展性目标进行声成像的方法。

背景技术

由于声波在水中传播的性能优于光波和电磁波，所以，应用于海洋资源的探测基本上是以声学手段为主。目前实际应用于声纳系统的空间处理技术仍然是以上世纪五、六十年代兴起并发展起来的以白噪声背景和自由场传播条件下的平面波波束形成。然而，在浅海环境下，传统的基于平面波模型的处理方法不再适用，由于声信号在海洋中的传播呈现多路径/多简正模状态，很容易造成检测过程中的伪目标或虚影现象，影响声纳系统的检测性能。时反处理从上世纪末开始应用于水声领域，是目前颇有应用前景的主动目标探测手段，时反算子分解是通过分析目标散射矩阵达到聚焦发射、定位以及成像目标的一种方法。

然而，目前时反算子分解方法是每个频点对应一个目标散射矩阵，每个频带单独处理的方式，存在如下缺点：(1)宽带信号的处理通常采用仅仅考虑中心频率或者整个频段的处理结果非相干累加的方式，没有有效地利用宽带信号的信息；(2)特别是波导环境中的扩展性目标，由于受目标本身的散射特性、波导传播效应以及换能器与频率的关系等因素的影响，最终表现出的目标散射特性是频率的函数，而目前处理过程中并没有考虑和利用这种散射特性。

发明内容

本发明要克服现有技术的上述缺点，提供一种加权宽带时反算子分解声成像方法，在现有的时反算子分解方法的基础上，有效地利用目标散射特性和宽带信息，从而达到消除伪目标以及有效显示扩展性目标空间特征的目的。

本发明提出的加权宽带时反算子分解声成像方法，提取每个子带模糊度函数的最大值以及相应的二维声场空间，根据最大特征值与频率的函数关系设计加权系数，加权子带的模糊度函数，并将加权后的模糊度函数相干累加，将最终结果显示为三维声场图像，确定目标的二维声场空间。

下面对本发明方法做进一步说明。

对环境的声学参数如水深、声速剖面、密度、沉积层参数等有足够了解的情况下，按如下步骤完成目标的三维声学图像及目标空间位置的识别：

(1)建立坐标系。设P个换能器组成的收发合置线阵，垂直放置在水中，以线阵为z轴，水平方向为r轴建立坐标系，水面为坐标原点。第一个换能器距水面的深度为z₁,第P个换能器距水面的深度为z_P；

(2)将换能器阵列接收到数据快拍排列成列向量，y_i(t)＝[y_1i(t)…y_pi(t)…y_Pi(t)]^T，p＝1,…,P，i＝1,…,T，T表示数据快拍数；

(3)根据预定采样频率对每个数据快拍进行采样，

y_i(n)＝[y_1i(n)…y_pi(n)…y_Pi(n)]^T；

(4)对采样后的信号进行短时傅里叶变换

$y_i\left(\mathrm{\ω}\right)=\left[\begin{array}{c}{\mathrm{\Σ}}_{m=1}^M{y}_{1i}\left(m\right)w(n-m)e^{-\mathrm{j\ωm}}\\ \·\\ \·\\ \·\\ {\mathrm{\Σ}}_{m=1}^M{y}_{\mathrm{pi}}\left(m\right)w(n-m)e^{-\mathrm{j\ωm}}\end{array}\right]---\left(\text{1}\right)$

其中w(n-m)为窗口函数序列；

(5)将所有的T个数据快拍排列成矩阵

$Y\left(\mathrm{\ω}\right)=\left[\begin{array}{ccc}{y}_{11}\left(\mathrm{\ω}\right)& \·\·\·& y_{1T}\left(\mathrm{\ω\ω}\right)\\ \·\·\·& \·\·\·& \·\·\·\\ y_{P1}\left(\mathrm{\ω}\right)& \·\·\·& y_{\mathrm{PT}}\left(\mathrm{\ω}\right)\end{array}\right]---\left(2\right)$

则时反算子为K(ω)＝Y(ω)Y^H(ω)；

(6)将时反算子划分为多个子带，进行奇异值分解：

K(ω)＝U(ω)Λ(ω)V*(ω) (3)

(7)输出各子带非相干累加后的特征值与水平距离r的函数曲线，根据最大特征值与水平距离的函数曲线判断包含目标的时间窗，找出所对应的时反算子K_o(ω)，特征值Λ_o(ω)以及特征向量U_o(ω)和V_o*(ω)；

(8)根据特征值Λ_o(ω)设计加权系数，比较所有子带的最大特征值λ₁(ω)，找出其最大峰值λ_max，其余各子带的最大特征值

与峰值λ_max的比值即为加权系数：

η(ω)＝λ₁(ω)/λ_max (4)

(9)将感兴趣的空间区域划分为网格，网格的交点即为假定目标所在的位置(r,z)，其中r表示假定目标距收发合置阵的水平距离，z表示假定目标的水深；

(10)根据成像方法所面临的水声环境，确定所适用的传播模型，分别得到假定目标与收发合置阵之间的传递向量

g(r,z,ω)＝[g₁(z₁,r,z,ω)…g_p(z_p,r,z,ω)…g_P(z_P,r,z,ω)]^T，

其中g₁(z₁,r,z,ω)表示假定目标与第一个换能器之间的传递函数，例如，简正模传播模型的传递函数为：

$g_1(z_1,r,z,\mathrm{\ω})=\mathrm{j\π}{e}^{\mathrm{j\ωt}}{\mathrm{\Σ}}_l{Z}_l\left(z_l\right)Z_l\left(z\right)H_0^{\left(2\right)}\left({\mathrm{\κ}}_{1}r\right)---\left(5\right)$

其中Z_l(·)表示第l号传播模对应的特征函数，H₀ ⁽²⁾(·)为汉克尔函数，κ_l表示第l号传播模对应的波数；

(11)各个子带的模糊度函数为：

I(r,z,ω)＝|g^H(r,z,ω)u₁(ω)|² (6)

其中u₁(ω)表示Λ_o(ω)的最大特征值λ₁(ω)对应的特征向量；

(12)提取各个子带模糊度函数的最大值I_max(r,z,ω)及其对应的声场空间位置(r_max(ω),z_max(ω))；

(13)将所有子带的声场空间(r_max(ω),z_max(ω))整理为一个集合Ω(r_max,z_max)，并重新定义模糊度函数：

$I(r,z,\mathrm{\ω})=\left\{\begin{array}{c}0,(r,z)\∉\mathrm{\Ω}(r_{\max },z_{\max })\\ I(r,z,\mathrm{\ω}),(r,z)\∈\mathrm{\Ω}(r_{\max },z_{\max })\end{array}\right.---\left(7\right)$

(14)对所有子带的模糊度函数进行加权累加：

I_end(r,z,ω)＝Ση(ω)I(r,z,ω) (8)

将最终获得的模糊度函数显示为距离r和深度z的三维图像，可确定目标所对应的二维空间区域。

本发明的优点在于将目标的散射特性通过加权的方式应用于宽带时反算子分解成像方法，真实地显示扩展性目标的空间特征，充分利用了宽带信息，并通过实验验证了该发明方案的有效性和可靠性。

附图说明

图1本发明的实验布置及坐标设置图。

图2本发明的特征值与距离的函数关系图。

图3本发明根据最大特征值设计的加权系数曲线。

图4为本发明目标在二维波导空间的定位图。

图5为本发明目标空间特征放大展示图。

具体实施方案

下面通过具体实施例子对本发明作进一步地描述。在一个深度为1.44m的波导实验室水池中，水池的三面贴消声瓦，水池的底部铺了一层0.22m厚的沙。水池的环境参数如下：水中声速为常数，通过测量水体温度计算得知为c₁＝1493m/s，水体密度ρ₁＝1000kg/m3，沉积层参数密度ρ₂＝1800kg/m3,声速为c₂＝1650m/s,衰减系数α₂＝0.67dB/λ，基底的密度为ρ₃＝1800kg/m3，声速c₃＝1580m/s，衰减系数α₃＝0.8dB/λ。

(1)本实施例子中，实验布置及坐标系的建立如图1所示，收发合置阵32个阵元，全场布阵，第一个阵元距水面0.04m，阵元间距0.04m，垂直放置在水池中。一个直径为0.21m，长度为0.51m的圆柱型目标放置在距阵约8m的水池底部。发射信号为脉冲宽度为1ms的中心频率为12kHz，带宽为4kHz的线性调频信号。

(2)将换能器阵列接收到数据快拍排列成列向量，

y_i(t)＝[y_1i(t)…y_pi(t)…y_Pi(t)]^T，p＝1,…,32，i＝1,…,8；

(3)根据采样频率为50kHz对每个数据快拍进行采样，

y_i(n)＝[y_1i(n)…y_pi(n)…y_Pi(n)]^T；

(4)选择的短时窗为发射脉宽的两倍，即点数为100的矩形窗，对采样后的信号进行短时傅里叶变换

$y_i\left(\mathrm{\ω}\right)=\left[\begin{array}{c}{\mathrm{\Σ}}_{m=1}^M{y}_{1i}\left(m\right)w(n-m)e^{-\mathrm{j\ωm}}\\ \·\\ \·\\ \·\\ {\mathrm{\Σ}}_{m=1}^M{y}_{\mathrm{pi}}\left(m\right)w(n-m)e^{-\mathrm{j\ωm}}\end{array}\right]---\left(\text{1}\right)$

(5)将所有的8个数据快拍排列成矩阵

$Y\left(\mathrm{\ω}\right)=\left[\begin{array}{ccc}{y}_{11}\left(\mathrm{\ω}\right)& \·\·\·& y_{1T}\left(\mathrm{\ω\ω}\right)\\ \·\·\·& \·\·\·& \·\·\·\\ y_{P1}\left(\mathrm{\ω}\right)& \·\·\·& y_{\mathrm{PT}}\left(\mathrm{\ω}\right)\end{array}\right]---\left(2\right)$

则时反算子为K(ω)＝Y(ω)Y^H(ω)；

(6)将4kHz的带宽的时反算子划分为401个子带进行奇异值分解：

K(ω)＝U(ω)Λ(ω)V*(ω) (3)

(7)输出各子带非相干累加后的特征值与水平距离r的函数曲线，根据最大特征值与水平距离的函数曲线判断包含目标的时间窗，找出所对应的时反算子K_o(ω)，特征值Λ_o(ω)以及特征向量U_o(ω)和V_o*(ω)；图2是前5个特征值与距离的函数关系图，以所有子带最大特征值的最大值规一化后以dB的形式显示，由上往下依次为最大、第二、第三、第四和第五特征值的函数关系曲线。由图2可知，最大特征值在目标所在的距离(8m)附近有一个局部峰值，从而可以确定目标所在的时间窗。

(8)根据特征值Λ_o(ω)设计加权系数，比较所有子带的最大特征值λ₁(ω)，找出其最大峰值λ_max，由图3所见，在中心频率12kHz附近，纵坐标显示为0dB，其余各子带的最大特征值与峰值λ_max的比值即为加权系数：

η(ω)＝λ₁(ω)/λ_max (4)

图3是根据最大特征值设计的加权系数，结果以dB形式显示。

(9)将感兴趣的空间区域划分为网格，水平起始距离为0.1m，搜索步长为0.05m，截止距离为15m，垂直起始深度为0m，搜索步长为0.0075m，截止深度为1.5m。网格的交点即为假定目标所在的位置(r,z)，其中r表示假定目标距收发合置阵的水平距离，z表示假定目标的水深；

(15)本例采用简正模传播模型，传递向量为

g(r,z,ω)＝[g₁(z₁,r,z,ω)…g_p(z_p,r,z,ω)…g_P(z_P,r,z,ω)]^T，

其中g₁(z₁,r,z,ω)表示假定目标与第一个换能器之间的传递函数，表达式为：

$g_1(z_1,r,z,\mathrm{\ω})=\mathrm{j\π}{e}^{\mathrm{j\ωt}}{\mathrm{\Σ}}_l{Z}_l\left(z_l\right)Z_l\left(z\right)H_0^{\left(2\right)}\left({\mathrm{\κ}}_{1}r\right)---\left(5\right)$

其中Z_l表示第l号传播模对应的特征函数，H₀ ⁽²⁾(·)为汉克尔函数，κ_l表示波数；

(10)各个子带的模糊度函数为：

I(r,z,ω)＝|g^H(r,z,ω)u₁(ω)|² (6)

其中u₁(ω)表示Λ_o(ω)的最大特征值λ₁(ω)对应的特征向量；

(11)提取各个子带模糊度函数的最大值I_max(r,z,ω)及其对应的声场空间位置(r_max(ω),z_max(ω))；

(12)将所有子带的声场空间(r_max(ω),z_max(ω))整理为一个集合Ω(r_max,z_max)，并重新定义模糊度函数：

$I(r,z,\mathrm{\ω})=\left\{\begin{array}{c}0,(r,z)\∉\mathrm{\Ω}(r_{\max },z_{\max })\\ I(r,z,\mathrm{\ω}),(r,z)\∈\mathrm{\Ω}(r_{\max },z_{\max })\end{array}\right.---\left(7\right)$

(13)对所有子带的模糊度函数进行加权累加：

I_end(r,z,ω)＝Ση(ω)I(r,z,ω) (8)

图4是目标在二维波导空间的定位图，即最终获得的模糊度函数与距离r和深度z的三维图像，展示了目标在整个水池所处的空间位置，图5目标空间特征放大展示图。

本说明书实施例所述的内容仅仅是对发明构思的实现形式的列举，本发明的保护范围不应当被视为仅限于实施例所陈述的具体形式，本发明的保护范围也及于本领域技术人员根据本发明构思所能够想到的等同技术手段。

Claims (1)

1.加权宽带时反算子分解声成像方法，包括如下步骤：

(1)建立坐标系；设P个换能器组成的收发合置线阵，垂直放置在水中，以线阵为z轴，水平方向为r轴建立坐标系，水面为坐标原点；第一个换能器距水面的深度为z₁,第P个换能器距水面的深度为z_P；

(2)将换能器阵列接收到数据快拍排列成列向量，

y_i(t)＝[y_1i(t)…y_pi(t)…y_Pi(t)]^T，p＝1,…,P，i＝1,…,T，T表示数据快拍数；

(3)根据预定采样频率对每个数据快拍进行采样，

y_i(n)＝[y_1i(n)…y_pi(n)…y_Pi(n)]^T；

(4)对采样后的信号进行短时傅里叶变换

$y_i\left(\mathrm{\ω}\right)=\left[\begin{array}{c}{\Σ}_{m=1}^M{y}_{1i}\left(m\right)w(n-m)e^{-j\mathrm{\ω}m}\\ .\\ .\\ .\\ {\Σ}_{m=1}^M{y}_{Pi}\left(m\right)w(n-m)e^{-j\mathrm{\ω}m}\end{array}\right]---\left(1\right)$

其中M表示窗的长度，w(n-m)为窗口函数序列；

(5)将所有的T个数据快拍排列成矩阵

$Y\left(\mathrm{\ω}\right)=\left[\begin{array}{ccc}{y}_{11}\left(\mathrm{\ω}\right)& ...& y_{1T}\left(\mathrm{\ω}\right)\\ ...& ...& ...\\ y_{P1}\left(\mathrm{\ω}\right)& ...& y_{PT}\left(\mathrm{\ω}\right)\end{array}\right]---\left(2\right)$

则时反算子为K(ω)＝Y(ω)Y^H(ω)；

(6)将时反算子划分为多个子带，进行奇异值分解：

K(ω)＝U(ω)Λ(ω)V*(ω) (3)

(7)输出各子带非相干累加后的特征值与水平距离r的函数曲线，根据最大特征值与水平距离的函数曲线判断包含目标的时间窗，找出所对应的时反算子K_o(ω)，特征值Λ_o(ω)以及特征向量U_o(ω)和V_o*(ω)；

(8)根据特征值Λ_o(ω)设计加权系数，比较所有子带的最大特征值λ₁(ω)，找出其最大峰值λ_max，其余各子带的最大特征值与峰值λ_max的比值即为加权系数：

η(ω)＝λ₁(ω)/λ_max (4)

(9)将感兴趣的空间区域划分为网格，网格的交点即为假定目标所在的位置(r,z)，其中r表示假定目标距收发合置阵的水平距离，z表示假定目标的水深；

(10)根据成像方法所面临的水声环境，确定所适用的传播模型，分别得到假定目标与收发合置阵之间的传递向量

g(r,z,ω)＝[g₁(z₁,r,z,ω)…g_p(z_p,r,z,ω)…g_P(z_P,r,z,ω)]^T，

其中g₁(z₁,r,z,ω)表示假定目标与第一个换能器之间的传递函数，简正模传播模型的传递函数为：

$g_1(z_1,r,z,\mathrm{\ω})={\mathrm{j\πe}}^{j\mathrm{\ω}t}{\Σ}_l{Z}_l\left(z_1\right)Z_l\left(z\right)H_0^{\left(2\right)}\left({\mathrm{\κ}}_{l}r\right)---\left(5\right)$

其中Z_l(·)表示第l号传播模对应的特征函数，H₀ ⁽²⁾(·)为汉克尔函数，κ_l表示第l号传播模对应的波数；

(11)各个子带的模糊度函数为：

I(r,z,ω)＝|g^H(r,z,ω)u₁(ω)|² (6)

其中u₁(ω)表示Λ_o(ω)的最大特征值λ₁(ω)对应的特征向量；

(12)提取各个子带模糊度函数的最大值I_max(r,z,ω)及其对应的声场空间位置(r_max(ω),z_max(ω))；

(13)将所有子带的声场空间(r_max(ω),z_max(ω))整理为一个集合Ω(r_max,z_max)，并重新定义模糊度函数：

$I(r,z,\mathrm{\ω})=\left\{\begin{array}{cc}0,& (r,z)\∉\mathrm{\Ω}(r_{\max },z_{\max })\\ I(r,z,\mathrm{\ω}),& (r,z)\∈\mathrm{\Ω}(r_{\max },z_{\max })\end{array}\right.---\left(7\right)$

(14)对所有子带的模糊度函数进行加权累加：

I_end(r,z,ω)＝Ση(ω)I(r,z,ω) (8)

将最终获得的模糊度函数显示为距离r和深度z的三维图像，可确定目标所对应的二维空间区域。