水下运动声源的接收信号仿真

提示：文章写完后，目录可以自动生成，如何生成可参考右边的帮助文档

---

前言

提示：这里可以添加本文要记录的大概内容：

项目地址：https://gitee.com/perry12341234/csdn-learn/tree/master/02SimRecSgnMvUWSdSrcs
时域信号的模拟技术目前使用较多的是频域卷积法，使用声场模型进行分频计算，与信号的傅里叶变换结果即频域在相应的频点上相乘，然后通过频域反傅里叶变换方法得到时域输出。
但是这样做在深海上往往计算量过大，计算速度慢，精度低，而且不稳定。更简单的方法是，使用BELLHOP设定带宽内的一个频点作为这个信号的中心频率，计算声场到达结构，忽略不同频点上到达结构的细微差距，使用到达结构幅度、相位和时延对信号做滤波，得到输出信号，这样做的速度快、计算精度高、稳定。

---

提示：以下是本篇文章正文内容，下面案例可供参考

一、理想情况下运动声源信号仿真技术

声源运动可以追溯到最简单的匀速直线运动，声波近似为平面波的情况。

假设发射信号为中心频率$f_c$的正弦信号，信号时长为$T$，其信号满足下面形式
$s(t)=\sin (2\pi f_{c}t),t\in (0,T)$
接收器静止位于$x=0$的位置，声源在做速度为$v$的匀速直线运动，初始位置位于$x=L$，在$t=0$时刻发射信号。
根据多普勒定律
接收信号满足
$r(t)=\sin (2\pi f_c\left(t-\frac{L}{c}-\frac{v}{c}t\right))=\sin (2\pi f_c\left(1-\frac{v}{c}\right)\left(t-\frac{L}{c}\right))$
从上面的公式可以看出，由于声源的位置和运动状态，时间因子$t$出现了时延和缩放。
在matlab中，可以使用重采样的方式实现。
归一化插值方法是比较文件的重采样方法，具体操作如下

1. 根据数字采样，得到发射信号的时间序列，并归一化
2. 按照多普勒的时间变化，得到接收信号时长，根据数字采样，得到接收信号时间序列，并将接收信号的时间序列归一化
3. 按照发射和接收时间序列归一化结果对信号进行插值
4. 将接收信号时间序列加上声线到达声源初始位置的时延$L/c$，并和插值后信号对应起来，得到接收信号和对应时间

clc;
% 运动声源的接收信号模拟
fs = 4000;
T = 1; % 信号时长
ts = 0: 1/fs: T; % 发射时间序列
fc = 100; % 信号中心频率
st = sin(2* pi * fc * ts); % 信号 s(t)
v = 10; % 声源运动速度(m/s)
L = 100; % 声源和接收的初始距离(m)
c = 1500; % 声速(m/s)
Tr = T * (1 + v/c); % 经过了多普勒后的信号时长
tau0 = L / c; % 信号到达的时间
tr = 0: 1/fs: Tr; % 重采样
trNor = tr / max(tr); % 归一化插值方法
tsNor = ts / max(ts);
sr = interp1(tsNor, st, trNor, 'linear', 'extrap'); % 接收信号
figure;
plot(ts, st, 'DisplayName', '发射信号');
hold on; 
plot(tr+tau0, sr, 'DisplayName', '接收信号');
legend();
ylim([-2 2]) 
xlabel('时间 (s)')
ylabel('幅度')

最终结果为：

二、水下运动声源信号仿真技术

1.原理

根据声场互易理论，可以以接收器为中心，使用射线理论计算水下声线，如上图所示，bellhop计算的声线可以根据下面6组参数确定，在每次步进中，射线可以近似为直线。因此对于运动声源，根据其运动状态，可以确定坐标随时间的关系。

$x_s=(r(t),z(t))$
进而得到声线与声源法相距离和径向长度随时间的关系
$m_s=(s(t),n(t))$
根据时延插值和声线贡献理论，可以得到时延和幅度随时间的关系
$A_s=(A(t),\tau (t))$
进而得到接收信号
$y(t)=A(t)\cdot s(t-\tau (t))$
因此，根据归一化插值理论可以求得该条声线下运动声源的接收信号。
最终的接收信号为所有声线的接收信号的和
$y(t)={\sum }_{n=1}^{n=Nbeams}{A}_n(t)\cdot s(t-{\tau }_n(t))$

2.代码

bellhop源代码中，.ray文件只输出了声线坐标，但是没有输出上表的其他参数，而且写入txt格式，比较浪费磁盘空间，因此我把bellhop整编了一下，使其输出二进制的声线参数。
代码如下（示例）：

    WRITE( BINRAYFile, REC = iRec1 ) alpha0, N2, ray2D( Nsteps1 )%NumTopBnc, ray2D( Nsteps1 )%NumBotBnc
    
    DO is = 1, N2
       iRec1 = iRec1 + 1
       WRITE( BINRAYFile, REC = iRec1 ) REAL(ray2D( is )%x(1), KIND = 4), REAL(ray2D( is )%x(2), KIND = 4), REAL(ray2D( is )%q(1), KIND = 4), &
           REAL(ray2D( is )%c, KIND = 4), real(ray2D( is )%tau, KIND = 4), REAL(imag(ray2D( is )%tau), KIND = 4), REAL(ray2D( is )%Amp, KIND = 4), &
           REAL(ray2D( is )%phase, KIND = 4)
    END DO

matlab读取ray参数后，仿真运动声源的信号代码如下，项目地址： 水下运动声源接收信号仿真

代码如下（示例）：

clc;

fs = 4000; % 采样率
T = 20; % 信号时长
ts = 0: 1/fs: T; % 发射时间序列
fc = 100; % 信号中心频率
st = sin(2* pi * fc * ts); % 信号 s(t)
v = 10; % 声源运动速度(m/s)

omega = 2 * pi * fc;
zs0 = 200; % 声源深度
rs0 = 10000; % 声源和接收的初始距离(m)

% rayfil = 'MunkB_ray_f.ray';
% [ PlotTitle, PlotType, freqVec, freq0, ray ] = readray(rayfil);
binrayfil = 'MunkB_ray_f.binray';
ray = readbinray(binrayfil);
Dalpha = ray(2).alpha0 - ray(1).alpha0; % 每次步进的角度

q0      = ray(1).c(1) / Dalpha;   % Reference for J = Q0 / Q
caustic = 1;   % caustic is multiplied by i after each caustic
TT = 40; % 经过海洋信道扩展后的信号长度一般不超过20s
reall = zeros(1,fs*TT);
for i = 1: length(ray)
    alpha = ray(i).alpha0;
    Rat1 = sqrt( abs( cos( alpha ) ) );
    % *** step along beams ***

    for is = 1 : max( ray(i).nsteps ) - 1
        if ( is > 1 )
            if ( ( ( ray( i ).q( is, 1 ) <= 0.0 ) && ( ray( i ).q( is - 1, 1 ) > 0.0 ) ) || ...
                    ( ( ray( i ).q(  is, 1 ) >= 0.0 ) && ( ray( i ).q( is - 1, 1 ) < 0.0 ) ) )
                caustic = 1i * caustic;   % apply the caustic phase change
            end
        end
        if ray( i ).x( is+1, 1 ) > rs0
            is0 = is;
            flag = 1;
            break;
        end
        flag = 0;
    end
    if ~flag % 如果找不到比起始声源距离大的声线，则跳过这条声线
        continue;
    end
    for is = is0: max( ray(i).nsteps ) - 1
        xray = ray( i ).x(is, :);
        rvec = v * ts.'; % 声源运动向量

        xrcvr  = [rs0 + rvec ones( length(rvec), 1 )*zs0];
        irr1 = find((xrcvr(:,1)<=ray( i ).x( is+1, 1 )) .* (xrcvr(:,1)>ray( i ).x( is, 1 ))); % 找到在该声线步进距离范围内的所有声源距离
        if isempty(irr1)
            break;
        end
        tray = ray( i ).x(is + 1, :) - ray( i ).x(is, :);
        rlen = norm( tray );

        % IF ( rlen < TINY( xV( 1, I ) ) ) CYCLE   % if duplicate point in ray, skip to next step along the ray
        tray        = tray / rlen;
        tray_scaled = tray / rlen;

        nray = [ -tray( 2 ) tray( 1 ) ];      % unit normal  to ray
        alpha1 = atan(tray(2)/ tray(1));



        % xrcvr  = [ xrcvr( ones( length( Pos.r.z ), 1 ) ) Pos.r.z' ];
        s      =      ( xrcvr(irr1,:) - xray ) * tray_scaled';  % proportional distance along ray
        n      = abs( ( xrcvr(irr1,:) - xray ) * nray' );       % normal distance to ray

        q      = ray( i ).q( is, 1 ) + s * ( ray( i ).q( is+1, 1 ) -  ray( i ).q( is, 1 ) );  % interpolated amplitude
        RadMax = abs( q / q0 );          	            % beam radius
        irr = find ( n < RadMax );
        if ~isempty(irr)
            delay = ray( i ).tau(is, :) + s * ( ray( i ).tau(is + 1, :) - ray( i ).tau(is, :) ); % interpolated delay 求时延
            amp = zeros(length(s), 1);
            % if (n < RadMax)
            A     = 1 ./ RadMax( irr );


            % shift phase for rays that have passed through a caustic
            caust = caustic * ones( 1, length( irr ) );
            ii    = ( ( q( irr ) <= 0.0 ) & ( ray( i ).q( is, 1 ) > 0.0 ) ) | ...
                ( ( q( irr ) >= 0.0 ) & ( ray( i ).q( is, 1 ) < 0.0 ) );
            caust( ii ) = 1i * caust( ii );   % apply the caustic phase change

            % const = Rat1 * ray( i ).Amp(is,1) * realsqrt( ray( i ).c(is) ./ abs( q( irr ) ) ) .* A .* caust.';
            const = Rat1 * ray( i ).Amp(is, 1) * realsqrt( ray( i ).c(is)./ abs( q( irr ) ) ) .* A .* caust.';

            amp(irr) = abs(const) .* ( RadMax( irr ) - n( irr ) )   .* exp( omega * delay(irr,2)) ;
            % end
            tr = ts(irr1).' + delay(:,1); % 虚部应该加到幅度上去
            tr1 = ceil(min(tr)*fs)/fs: 1/fs: floor(max(tr)*fs)/fs; % 接收的采样点，需要按fs进行采样
            % trNor = (tr - min(tr)) / max(tr - min(tr));
            % trNor = tr / max(tr); % 归一化插值方法
            tsNor = ts / max(ts);
            sr = interp1(tr, amp.*st(irr1).', tr1, 'linear', 'extrap'); % 接收信号
            reall(ceil(min(tr)*fs):  floor(max(tr)*fs)) = reall(ceil(min(tr)*fs):  floor(max(tr)*fs)) + sr;
            % yt =  sr;
        end
    end
end
tt = (0: TT*fs-1)/fs;
figure;
plot(tt, reall)
figure;
stft(reall, fs);

仿真的时域和时频谱结果如下，可以看出信号在某些部分还是有些不连贯（算法有待优化），能看出干涉效应。

三、 连续信号仿真

由上面的仿真结果可以看出，仿真的信号存在不连续，从理论中分析原因，不连续主要是由声线步进之间的方向向量不连续导致的，声线步进过程，在上面的程序中，使用直线近似的，那么在步进的边界点处不光滑，也就是声线轨迹不是处处可导的，声线轨迹导数是方向，方向向量不是连续的。
根据snell定律，可以证明，在声速梯度为常数，即声速剖面是直线的条件下，声线的轨迹是圆
根据步进思想
$\{\begin{array}{ll}\frac{dr}{ds}=\cos \theta =c\xi & \\ \frac{dz}{ds}=\sin \theta =c\zeta & \end{array}$
只需要在每一个步进声线内，把声线当作圆来处理，即可得到光滑的声线轨迹。
在每一个步进声线内，根据声线两端的方向向量和声线弧长，可以计算得到圆心坐标和半径
$R=\frac{1}{\kappa }=dS/d\phi$
${\mathbf{x}}_{\mathbf{c}}={\mathbf{x}}_0+R\ast \mathbf{n}$
根据目标位置、圆心、半径、以及声线起点坐标，可以计算得到目标位置在声线上的百分比映射
${\phi }_x=angle(\mathbf{x}-{\mathbf{x}}_{\mathbf{c}})$
${\phi }_{x0}=angle({\mathbf{x}}_0-{\mathbf{x}}_{\mathbf{c}})$
${\phi }_{x1}=angle({\mathbf{x}}_1-{\mathbf{x}}_{\mathbf{c}})$
$s_x=\frac{{\phi }_x-{\phi }_{x_0}}{{\phi }_{x_1}-{\phi }_{x0}}$

其他代码与上文保持一致。

clc;
% ! bellhop2D.exe MunkB_ray_f
fs = 4000; % 采样率
T = 20; % 信号时长
TT = 40; % 经过海洋信道扩展后的信号长度一般不超过20s
ts = 0: 1/fs: T; % 发射时间序列
fc = 100; % 信号中心频率
B = 50; % 带宽
% st = sin(2* pi *  (fc - B/2) * ts + pi *  B / T * ts.^2); % 信号 s(t)
% st = ones(1, length(ts));
st = sin(2* pi * fc * ts); % 信号 s(t)
v = 5; % 声源运动速度(m/s)


zs0 = 200; % 声源深度
rs0 = 50000; % 声源和接收的初始距离(m)
PT = rs0/1800;

omega = 2 * pi * fc;
% rayfil = 'MunkB_ray_f.ray';
% [ PlotTitle, PlotType, freqVec, freq0, ray ] = readray(rayfil);
binrayfil = 'MunkB_ray_f.binray';
ray = readbinray(binrayfil);
Dalpha = ray(2).alpha0 - ray(1).alpha0; % 每次步进的角度

q0      = ray(1).c(1) / Dalpha;   % Reference for J = Q0 / Q
caustic = 1;   % caustic is multiplied by i after each caustic

reall = zeros(1,fs*TT);
for i = 1: length(ray)
    alpha = ray(i).alpha0;
    Rat1 = sqrt( abs( cos( alpha ) ) );
    x = ray(i).x;
    t = ray(i).t;
    
    dx = x(2:end,:) - x(1:end-1,:);
    L = sqrt(sum(dx.^2, 2));
    dxnorm = dx./sqrt(sum(dx.^2,2));
    tnorm = t./sqrt(sum(t.^2,2));
    nnorm = [-dxnorm(:,2) dxnorm(:,1)];
    ntnorm = [-tnorm(:,2) tnorm(:,1)];
    
    phi = atan(tnorm(:, 2)./ tnorm(:, 1));
    dphi = phi(2:end) - phi(1:end-1);
    k = dphi ./ L;
    R = 1./k;
    xo = x(1:end-1,:) + R*ones(1,2) .* ntnorm(1:end-1,:);
    % *** step along beams ***

    for is = 1 : max( ray(i).nsteps ) - 1
        if ( is > 1 )
            if ( ( ( ray( i ).q( is, 1 ) <= 0.0 ) && ( ray( i ).q( is - 1, 1 ) > 0.0 ) ) || ...
                    ( ( ray( i ).q(  is, 1 ) >= 0.0 ) && ( ray( i ).q( is - 1, 1 ) < 0.0 ) ) )
                caustic = 1i * caustic;   % apply the caustic phase change
            end
        end
        if ray( i ).x( is+1, 1 ) > rs0
            is0 = is;
            flag = 1;
            break;
        end
        flag = 0;
    end
    if ~flag % 如果找不到比起始声源距离大的声线，则跳过这条声线
        continue;
    end
    for is = is0: max( ray(i).nsteps ) - 1
        xray = ray( i ).x(is, :);
        rvec = v * ts.'; % 声源运动向量

        xrcvr  = [rs0 + rvec ones( length(rvec), 1 )*zs0];
        
        tray = ray( i ).x(is + 1, :) - ray( i ).x(is, :);
        rlen = norm( tray );

        % IF ( rlen < TINY( xV( 1, I ) ) ) CYCLE   % if duplicate point in ray, skip to next step along the ray
        tray        = tray / rlen;
        tray_scaled = tray / rlen;

        % tray_scaled = ray(i).t(is,:);
        % tray = tray_scaled ./ norm(tray_scaled);

        nray = [ -tray( 2 ) tray( 1 ) ];      % unit normal  to ray
        alpha1 = atan(tray(2)/ tray(1));



        % xrcvr  = [ xrcvr( ones( length( Pos.r.z ), 1 ) ) Pos.r.z' ];
        % s      =      ( xrcvr(:,:) - xray ) * tray_scaled';  % proportional distance along ray
        xstep = x(is,:);
        xstep1 = x(is+1,:);
        xo1 = xo(is,:);
        rsvec = xrcvr - xo1;
        r0vec = xstep - xo1;
        r1vec = xstep1 - xo1;
        Rc = R(is);
        if abs(Rc) > 1e7
            tnorm_scale = dxnorm(is,:) / L(is);
            s = (xrcvr - xstep)* tnorm_scale.';
        else
            phis = angle(rsvec(:,1)+ 1j* rsvec(:,2));
            phi0 = angle(r0vec(1)+ 1j* r0vec(2));
            phi1 = angle(r1vec(1)+ 1j* r1vec(2));
            s = (phis - phi0) / (phi1 - phi0);
        end
        irr1 = find(s>0& s<=1); % 找到在该声线步进距离范围内的所有声源距离
        if isempty(irr1)
            break;
        end
        n      = abs( ( xrcvr(irr1,:) - xray ) * nray' );       % normal distance to ray
        % s = ( xrcvr(irr1,1) - ray( i ).x(is, 1) ) / (ray( i ).x(is + 1, 1) - ray( i ).x(is, 1));
        % n = abs(( xrcvr(irr1,2) - ray( i ).x(is, 2) ) + s * (ray( i
        % ).x(is + 1, 2) - ray( i ).x(is, 2)));1
        q      = ray( i ).q( is, 1 ) + s(irr1) * ( ray( i ).q( is+1, 1 ) -  ray( i ).q( is, 1 ) );  % interpolated amplitude
        RadMax = abs( q / q0 );          	            % beam radius
        irr = find ( n < RadMax );
        if ~isempty(irr)
            delay = ray( i ).tau(is, :) + s(irr1) * ( ray( i ).tau(is + 1, :) - ray( i ).tau(is, :) ); % interpolated delay 求时延
            amp = zeros(length(s(irr1)), 1);
            % if (n < RadMax)
            A     = 1 ./ RadMax( irr );


            % shift phase for rays that have passed through a caustic
            caust = caustic * ones( 1, length( irr ) );
            ii    = ( ( q( irr ) <= 0.0 ) & ( ray( i ).q( is, 1 ) > 0.0 ) ) | ...
                ( ( q( irr ) >= 0.0 ) & ( ray( i ).q( is, 1 ) < 0.0 ) );
            caust( ii ) = 1i * caust( ii );   % apply the caustic phase change

            % const = Rat1 * ray( i ).Amp(is,1) * realsqrt( ray( i ).c(is) ./ abs( q( irr ) ) ) .* A .* caust.';
            const = Rat1 * ray( i ).Amp(is, 1) * realsqrt( ray( i ).c(is)./ abs( q( irr ) ) ) .* A .* caust.';

            amp(irr) = abs(const) .* ( RadMax( irr ) - n( irr ) )   .* exp( omega * delay(irr,2)) ;
            % end
            tr = ts(irr1).' + delay(:,1) - PT; % 虚部应该加到幅度上去
            tr1 = ceil(min(tr)*fs)/fs: 1/fs: floor(max(tr)*fs)/fs; % 接收的采样点，需要按fs进行采样
            % trNor = (tr - min(tr)) / max(tr - min(tr));
            % trNor = tr / max(tr); % 归一化插值方法
            tsNor = ts / max(ts);
            sr = interp1(tr, amp.*st(irr1).', tr1, 'linear', 'extrap'); % 接收信号
            reall(ceil(min(tr)*fs):  floor(max(tr)*fs)) = reall(ceil(min(tr)*fs):  floor(max(tr)*fs)) + sr;
            % yt =  sr;
        end
    end
end
tt = (0: TT*fs-1)/fs;
figure;
plot(tt, reall)
figure;
stft(reall, fs);

接收信号的结果为
与上文对比，可以发现明显少了很多毛刺

四、 被动定位baseline

1. 四种典型海深

浅海（200m） 过渡海域（200-1000m） 南海（4000m） 西太深海（6000m）

2. 声纳或者基阵参数

舰艇上常用的声纳
球阵、线阵、面阵
布放或者吊放阵
浅海坐底圆环阵、垂直阵、圆柱阵、深海吊放垂直阵、线阵

3. 常见的声源或者目标

实验用的爆炸声源（高声源级）、换能器发射的各种信号（CW、LFM等，信噪比一般，但是可以使用脉冲压缩等方法）、水面的船（声源级比较小，但频带宽，连续的信号）、水下的潜艇（声源级非常小，但是一般在水下）

4. 常见的定位场景

在浅海0-30km作DOA和距离深度估计（浅海深度定位更难）
在深海直达区或者汇聚区DOA和距离深度估计（直达区和汇聚区的传播损失更小）
深海声影区的定位（作为附加题作为研究）

4.1 浅海

注意现在用的是bellhop模型，得注意一下频段
$f>c/H=1500/200=7.5Hz$
基本上满足

1. 地形
   平坦 200m
2. 声速剖面
   WOA读取一个3月份
3. 声纳
   坐底圆环阵，半径7.34m，37阵元，均分，采样率4000
4. 声源
   换能器 75-150Hz LFM信号，脉冲时间1s 140dB
5. 场景
   大地坐标系11.1°，深度88.8m，速度2.2m/s，方向33.1°
6. 任务
   方位角、距离、深度
7. 参数级粗计算
   声纳方程
   $TL=15lgr=60dB$
   $NL = 68 + lg(150-75) \approx 87dB $
   $DI = 10lgM \approx 16 $
   $DT=SL-TL-NL+DI\approx 140-60-87+16=9dB$

4.2 过渡海域

200-1000m从etopo1中选择一个比较曲线的地形，比较真实

4.3 南海

4000m 平底
声速使用WOA读取一个7月份的

4.4 西太深海

5000m 平底
Munk声速剖面（西太经典声速剖面）

1. 地形
   5000m 平坦

2. 声速剖面
   Munk声速剖面
   $c=c_0\ast (1.0+0.00737\ast (z_t-1.0+\exp (-z_t)))$
   $z_t=2.0\ast (z-1300.0)/1300.0$
   其中$c_0=1500.0$

3. 声纳
   垂直阵 全海深 M = 35
   88 138 188 238 288 338 388 438 488 538 588 638 688 738 788 838 888 938 988 1033 1133 1233 1333 1433 1533 1633 1733 1833 1933 2077 2577 3077 3577 4077 4577
   球阵 阵元数 10*10 深度 200m

4. 声源
   换能器 LFM信号 83-127Hz 2.8s 145dB
   舰艇 目标辐射噪声，150dB@100Hz，20s

5. 场景
   大地坐标系25.3°，深度13.2m，速度4.1m/s，8.2km接近

6. 任务
   方位角、距离、深度

7. 参数级粗计算
   声纳方程
   $TL=20lgr\approx 78dB$
   $NL = 68 + lg(150-75) \approx 87dB $
   $DI1 = 10lgM \approx 16 $
   $DI2 = 10lgM \approx 20 $
   $DT1=SL-TL-NL+DI1\approx 150-78-87+16=1dB$
   $DT2=SL-TL-NL+DI2\approx 150-78-87+20=5dB$
   
   使用LFM得到的声纳输出信号如图所示。

总结

提示：这里对文章进行总结：
以上就是今天要讲的内容，本文仅仅简单介绍了水下运动声源的接收信号仿真，使用射线法计算声场，使用射线参数，对幅度、时延等参数进行计算，得到了较为准确的水下运动声源的接收信号模拟。