本发明属于水声工程、海洋工程和声纳技术等领域，涉及一种深海低信噪比条件下声源深度自主实时跟踪方法，是一种在深海大深度接收环境中，低信噪比条件下，实现声源深度自主实时跟踪的方法，适用于深海环境中低声源级宽带运动声源深度的自主跟踪。
背景技术：
深海环境中，深海等温层的声速随着深度的增加而增加，至临界深度(典型值为4500m-5500m)，声速与海表面附近的最大声速相同。近海面声源至临界深度以下水听器之间的直达波(和海面反射波)被称为可靠声路径。广义上讲，当不存在临界深度但海深较深时(>2000m)，近海底水听器也可接收到中等距离近海面声源的直达波信号，也可称为可靠声路径，但这种条件下的低频海洋环境噪声级会增加。可靠声路径由于其传播损失低，近海面覆盖范围广且不存在盲区(5倍的海深以内)，将潜标布放在近海底附近，用于水下声源探测，具有重要的军事应用价值。由于这种潜标是无人值守的，声源检测、距离估计和深度估计法均需要自主完成。声源检测和距离估计可以通过同步垂直阵的方位估计结合射线理论自主完成，而声源深度的自主估计难度较大。近海面声源的直达波和海面反射波的时延差随声源深度增加而增加。例如在近距离处，水下50m和5m声源时延差相差近10倍。因此，通过时延差是估计声源深度的一个重要途径。对于信号形式未知的宽带信号，自相关函数会在时延差位置形成峰值(称为时延峰)，从而获得两种路径的时延差。但是，在低信噪比环境下，一次自相关函数的时延峰的幅度和噪声峰值的幅度相差不大甚至更低，因此需要利用时间累积增益。
通过滤波方法，以每一个时刻的自相关函数为一次观测，是利用时间累积增益实现时延差估计的一类优化算法。常用的有卡尔曼滤波类方法和粒子滤波类方法。基于粒子滤波的方法参见“particlefilteringforarrivaltimetrackinginspaceandsourcelocalization”，该文2012年11月发表在《thejournaloftheacousticsocietyofamerican》第132期，起始页码为3041。但这些滤波类方法需要建立声源运动状态方程，并且对噪声协方差矩阵、观测误差协方差矩阵等参数的选取均比较敏感，往往需要人工调整参数来保证算法的性能。此外，这一类方法的运算量大，硬件实现后导致系统功耗、存储需求过高，无法应用于长期无人值守的低功耗潜标系统。因此，目前仍缺少可应用于大深度接收环境的声源深度自主实时跟踪方法。
技术实现要素：
要解决的技术问题
为了避免现有技术的不足之处，本发明提出一种深海低信噪比条件下声源深度自主实时跟踪方法，克服在实际低信噪比环境下，现有基于时间累积效应的声源深度跟踪算法无法自主工作的缺点。
技术方案
一种深海低信噪比条件下声源深度自主实时跟踪方法，其特征在于：接收同步垂直阵列布放在近海底，辐射宽带信号的声源与垂直阵的水平距离小于35km；实时跟踪步骤如下：
步骤1：令c(z)为本地的声速剖面，z表示深度；将声速剖面0-300m，均匀分割为i层，每层下界面的深度为zb,i，i＝1，2，…，i；其中i为从上往下声速层的编号；定义zb,0为海面的深度，设置为0m；利用下式计算每一个声速层的线性回归系数：
利用下式计算每一个声速层的参考声速：
设定近海面声速平均值c0为近海面任意选取的一个深度范围内声速的平均值，c(zr)为接收阵中心位置的声速；
步骤2：垂直阵采集长度为1s-10s的信号，信号的采样率为fs，取值范围1khz-100khz；对每一个阵元采集的信号进行傅里叶变换，得到第m个阵元处频点fl处的频谱xm(fl)，l＝1,2,…l；f1和fl为宽带波束形成所选用频率范围的上下界，频点fl处的信号向量为x(fl)＝[x1(fl),x2(fl),…xm(fl)]t，上标t表示转置；利用下式计算频点fl处的信号协方差矩阵rx(f)：
rx(fl)＝x(fl)xh(fl)
其中x(f)为垂直阵接收到的数据向量，上标h表示共轭转置；
利用公式计算宽带方位谱：
记录p(θ)最大值对应的方向θmax；
所述
步骤3：利用下式构造频点fl处的空域滤波器
利用下式估计直达波和海面反射波叠加的波形y(t)
其中，
然后利用以下自相关处理方法估计接收信号的自相关函数：
其中，t为空域滤波后输出信号的长度，τ为时延坐标；
由于信号的采样率为fs，则τ的分辨率为1/fs；记离散的时延坐标值为{τp}，p＝1,2,….；
步骤4：设置初始条件：i＝1，参数γ0＝0，与时延坐标集合{τp}对应的声源深度集合设置为{zs(τp)}＝nan，其中nan表示无效值；
步骤5：利用下式计算参考出射角θref,i
依次计算下列公式，得到参数κi，ai，bi，ξi：
ai＝2qiκi+2cscθref,i
bi＝4piκic0
步骤6：利用下式计算每个时延坐标值τp对应的声源深度
当时延坐标τp满足zs,i(τp)≥zb,i-1，则对应的声源深度集合更新为
zs(τp)＝zs,i(τp).
若i为最后一个声速层，进入下一步骤；否则，计算下式
同时，令i＝i+1，返回至步骤5，在下一个声速层中计算声源深度；
步骤7：将接收信号的自相关函数利用下式表示为深度相关函数
e(zs)＝c(τ(zs)),
其中τ(zs)是zs(τ)的反函数，直接由声源深度集合{zs(τp)}给出；
步骤8：每隔5s～30s，重复步骤2至步骤7，得到深度相关函数序列{e1(zs),…,ek(zs),…,eq(zs)}，其中，q为平滑窗的长度，取值范围为1-50，e1(zs)为起始时刻的深度相关函数，ek(zs)为k时刻的深度相关函数；
利用下式计算过去q个深度相关函数的平均值函数：
则q时刻的声源深度估计值为平均值函数最小值点对应的深度，即
步骤9：每隔5s-30s，获得一个新的深度相关函数，将其加入函数序列，并将该序列中时刻最早的深度相关函数剔除，此时的序列更新为{ek-q+1(zs),…,ek(zs)}，k为当前时刻；
根据步骤8的公式计算该时刻序列的平均值函数，然后估计声源深度zs,k；
依次类推，得到实时的声源深度序列{zs,q,…,zs,k,…,zs,k}，k为总的观测时长。
所述接收同步垂直阵在近海底为深度2000-6000m。
所述声源深度范围小于1000m。
所述接收同步垂直阵阵元数为m，取值范围为8-128，阵元编号从浅到深依次增加。
有益效果
本发明提出的一种深海低信噪比条件下声源深度自主实时跟踪方法，该方法的系统配置为同步垂直阵。本发明首先将阵列布放在深海大深度位置，在某一时刻检测到声源信号，接收声源辐射的宽带信号，计算接收信号到达角；之后利用空域滤波获得该方向波束内的信号，并求得信号的自相关函数；将自相关函数的时延自变量表示为声源深度的函数，则原自相关函数转换为以声源深度为因变量的复合函数，称为深度相关函数；将当前时刻的深度相关函数作为平滑滤波器的输入，则输出为过去一段时间长度内深度相关函数的平均值函数；深度相关函数的背景噪声由于具有随机性而无法累加，而真实声源深度附近的峰值会在该平均值函数内累积，其对应的深度值即为该时刻估计的声源深度；下一时刻的深度相关函数进入平滑滤波器可估计下一时刻的声源深度，依此类推实现声源深度跟踪。本发明利用了下述性质解决了声源深度实时自主跟踪问题：在深海大深度接收环境下，海面反射-直达波时延可以显式的表示为声源深度、到达角和声速剖面为参量的方程。因此，将时延作为中间变量，将该方程带入自相关函数中，便得到了相关能量(自相关函数的因变量)随假设声源深度变化的函数。本发明的关键是时延方程，其所需的声速剖面参数可以实时测量，到达角参数可以由垂直阵的宽带波束形成实现自主估计，无需人工干预，从而实现自主的声源深度跟踪。另外，在任意时刻，仅存储了最近一段时间内的深度相关函数，存储量小，便于实时处理。
有益效果是：该方法首先在步骤1中对声速剖面进行了分段线性近似处理，在步骤2中得到信号在垂直阵上的到达角信息，然后利用该信息在步骤3中获得了信号自相关函数，其中包含声源深度信息；步骤4至步骤7利用低运算量的若干方程获得了自相关函数的每个时延坐标所对应的声源深度，从而将自相关函数变换为以假设声源深度为自变量的函数，称为深度相关函数；最后在步骤8中将平滑滤波器应用于深度相关函数的时间序列，削弱背景噪声，提取信号的相关峰，跟踪声源深度。本方法的有益效果主要基于以信号到达角、直达-海面波时延、声速剖面为参数的声源深度计算公式，该公式将自相关函数变换为随假设声源深度变化的深度相关函数，与传统基于射线模型的声源深度估计方法相比，计算量得到大幅度降低；此外，简单并且可以实时运行的平滑滤波器适用于该函数的时间序列，避免使用复杂的滤波跟踪算法，从而用很小的计算量和存储代价便获得了时间累积增益，提高了低信噪比下本发明的稳健性。
附图说明
图1是深海低信噪比条件下声源深度自主实时跟踪方法工作环境示意图。(a)声源-水听器在大深度接收环境中的几何关系示意图；(b)仿真环境的声速剖面。
图2是本发明方法声源运动路线的俯视图
图3是本发明方法阵元接收信号信噪比的平均值随声源距离变化的曲线。
图4是本发明方法多段线性近似的声速剖面(虚线)与仿真使用的声速剖面(实线)的对比图。
图5是本发明方法自相关函数随距离/信号到达角变化的灰度图。
图6是本发明方法深度相关函数随距离/信号到达角变化的灰度图。
图7是本发明方法应用于仿真实验的结果。(a)图6中的深度相关函数的时间序列通过平滑滤波器后的输出序列。(b)本发明方法声源深度跟踪结果(虚线)和实际声源深度(实线)的对比图。
图8是本发明方法应用于跟踪一艘海上航船深度的实验结果。(a)深度相关函数随航船信号到达角变化的灰度图。(b)子图(a)中的深度相关函数的时间序列通过平滑滤波器后的输出序列，及本发明方法声源深度跟踪结果(虚线)。
具体实施方式
现结合实施例、附图对本发明作进一步描述：
同步采集垂直阵列布放在近海底附近，接收近海面运动声源辐射的宽带信号，并通过波束形成、空域滤波、自相关处理和函数变换方法，实现声源深度的跟踪，其过程分为以下步骤：
声源进入大深度接收环境范围内，同步垂直阵接收到声源辐射的宽带信号，此时声源和垂直阵的水平距离为0km～35km，声源运动状态(包括距离、深度和速度等)未知。垂直阵阵元数为m，取值范围为8-128，阵元编号从浅到深依次增加。
步骤1：令c(z)为本地的声速剖面，z表示深度。将声速剖面0-300m，均匀分割为i层，每层下界面的深度为zb,i，i＝1，2，…，i，其中i为从上往下声速层的编号。定义zb,0为海面的深度，设置为0m。利用下式计算每一个声速层的线性回归系数
利用下式计算每一个声速层的参考声速
设定近海面声速平均值c0为近海面任意选取的一个深度范围内声速的平均值，c(zr)为接收阵中心位置的声速。
步骤2：垂直阵采集长度为1s-10s的信号。信号的采样率为fs，取值范围1khz-100khz。对每一个阵元采集的信号进行傅里叶变换，得到第m个阵元处频点fl处的频谱xm(fl)，l＝1,2,…l。f1和fl为宽带波束形成所选用频率范围的上下界。频点fl处的信号向量为x(fl)＝[x1(fl),x2(fl),…xm(fl)]t，上标t表示转置。利用下式计算频点fl处的信号协方差矩阵rx(f)
rx(fl)＝x(fl)xh(fl)(3)
其中x(f)为垂直阵接收到的数据向量，上标h表示共轭转置。利用公式
计算宽带方位谱，记录p(θ)最大值对应的方向θmax。
所述
步骤3：利用下式构造频点fl处的空域滤波器
利用下式估计直达波和海面反射波叠加的波形y(t)
其中，
其中，t为空域滤波后输出信号的长度，τ为时延坐标。由于信号的采样率为fs，则τ的分辨率为1/fs。记离散的时延坐标值为{τp}，p＝1,2,….
步骤4：设置初始条件：i＝1，参数γ0＝0，与时延坐标集合{τp}对应的声源深度集合设置为{zs(τp)}＝nan，其中nan表示无效值。
步骤5：利用下式计算参考出射角θref,i
依次计算下列公式，得到参数κi，ai，bi，ξi：
ai＝2qiκi+2cscθref,i(10)
bi＝4piκic0(11)
步骤6：利用下式计算每个时延坐标值τp对应的声源深度
当时延坐标τp满足zs,i(τp)≥zb,i-1，则对应的声源深度集合更新为
zs(τp)＝zs,i(τp).(14)
若i为最后一个声速层，进入步骤八。否则，计算下式
同时，令i＝i+1，返回至步骤六，在下一个声速层中计算声源深度。
步骤7：将公式(7)得到的自相关函数利用下式表示为深度相关函数
e(zs)＝c(τ(zs)),(16)
其中τ(zs)是zs(τ)的反函数，直接由声源深度集合{zs(τp)}给出。
步骤8：每隔5s～30s，重复步骤三至步骤八，得到深度相关函数序列{e1(zs),…,ek(zs),…,eq(zs)}，其中，q为平滑窗的长度，取值范围为1-50，e1(zs)为起始时刻的深度相关函数，ek(zs)为k时刻的深度相关函数。利用下式计算过去q个深度相关函数的平均值函数
则q时刻的声源深度估计值为平均值函数最小值点对应的深度，即
步骤9：每隔5s-30s，获得一个新的深度相关函数，将其加入函数序列，并将该序列中时刻最早的深度相关函数剔除，此时的序列更新为{ek-q+1(zs),…,ek(zs)}，k为当前时刻。参照公式(17)和(18)计算该时刻序列的平均值函数，然后估计声源深度zs,k。依次类推，得到实时的声源深度序列{zs,q,…,zs,k,…,zs,k}，k为总的观测时长。
具体实施例：
参照图1。水声信道建模：仿真环境的声速剖面为典型的munk声速剖面，仿真环境的海深为5500m，临界深度为4800m，海底参数为声速1600m/s，密度1.65g/cm3，吸收损失0.1db/λ，其中λ为波长。接收阵为16元垂直线列阵，阵元间距为5m，阵列中心的深度为3700m。
参照图2。声源运动建模：声源由远及近运动，在15km处(a点)开始声源深度跟踪。声源最近通过点距离潜标的水平距离为5km。之后，声源由近及远运动，至15km处(d点)跟踪结束。在a至b点声源深度100m；b至c点，声源深度匀速变为50m；之后保持不变。声源运动速度恒定为5m/s。声源辐射50hz至500hz的宽带随机信号，在该带宽内，信号的功率谱是平坦的。
参照图3。接收信号建模：在某一声源-接收点位置下，取声源辐射的宽带随机信号时长为4s，采样率为5khz，对其做傅里叶变换，得到fs(f)，然后采用射线模型软件bellhop计算所有频点处的信道响应，记为fc(f)，最后fs(f)和fc(f)对应频点相乘并做反傅里叶变换，得到水听器接收到的无干扰信号。将高斯白噪声通过通带为50hz-500hz的带通滤波器。然后将噪声进行相应放大或衰减并与无干扰信号叠加，以获得相应低信噪比的信号。仿真得到信号的信噪比均低于0db，并且随着距离增加，信噪比进一步降低。
具体跟踪过程如下：
声源进入大深度接收环境范围内，水听器接收到声源辐射的宽带信号，此时声源距离接收阵的水平距离为15km，声源运动状态(包括距离、深度和速度等)未知。
步骤1参照图4。将声速剖面0-300m，均匀分割为5层，则声速层上下边界深度zb,i＝i*60。由公式(1)和公式(2)计算每一个声速层的回归系数{pi,qi}和参考声速{cref,i}。声速平均值c0为1545m/s，接收阵中心位置的声速c(zr)为1530.7m/s。
步骤2将阵元接收的4s长度的信号，进行傅里叶变换。然后利用公式(3)和(4)计算信号到达的方向θmax，其中频率上下界分别f1＝50hz，fl＝500hz。
步骤3参考图5。利用公式(5)构造空域滤波器，然后利用公式(6)获得波束输出，最后利用公式(7)计算接收信号的自相关函数c(τ)，将τ的范围限制在0-0.08s。由于采样率为5khz，则τ的分辨率为0.2ms，记离散的时延坐标值为{τ1,…,τp,…，τ401}，其中τ1＝0s，τ401＝0.08s。
步骤4设置初始条件：i＝1，参数γ0＝0，声源深度集合{zs(τp)}＝nan。
步骤5利用公式(8)计算参考出射角θref,i，利用公式(9)至(12)计算参数κi，ai，bi，ξi。
步骤6利用公式(13)计算每个时延坐标值τp对应的声源深度zs,i(τp)，若满足条件zs,i(τp)≥zb,i-1，则利用公式(14)更新对应的声源深度集合。若i＝5，则进入步骤八。否则利用公式(15)计算参数γi，同时，令i＝i+1，返回至步骤六，在下一个声速层中计算声源深度。
步骤7参考图6利用公式(16)计算深度相关函数e1(zs)。
步骤8每隔10s，重复步骤三至步骤八，得到深度相关函数序列{e1(zs),…,ek(zs),…,eq(zs)}，其中q选取为12。利用公式(17)计算该时刻深度相关函数序列的平均值函数，然后利用公式(18)估计声源深度zs,12。
步骤9参考图7。每隔10s，获得一个新的深度相关函数，将其加入函数序列，并将序列中时刻最早的深度相关函数剔除，此时的序列为{ek-q+1(zs),…,ek(zs)}，k为当前时刻。参照公式(17)计算该时刻深度相关函数序列的平均值函数，然后参照公式(18)估计声源深度zs,k。由此，得到实时的声源深度序列{zs,q,…,zs,k,…,zs,k}，k为总的观测时长。
参照图8。本发明应用于一次海上试验，被跟踪的声源为沿直线匀速运动的一艘航船，海深为3950m。接收系统为布放在近海底附近的均匀垂直线列阵，阵元数16，阵元间距5m，阵列中心深度3715m。在航船接近过程中，距离潜标水平距离8km时开始跟踪；在远离过程中，水平距离13km时，停止跟踪。航船的最近通过距离为4km。跟踪过程共采集50组信号，信号带宽为700-1500hz。本方法将自相关函数变换为深度相关函数，如子图(a)所示，可以看出在11.6m深度附近出现幅度较高的条纹。当平滑窗q的长度设置为12时，本方法的声源深度跟踪结果如子图(b)所示，与条纹吻合较好。由航船的尺寸、声源级等参数分析，该航船的螺旋桨中心的深度较深，在10m-15m之间，与本方法的跟踪结果相符。