国防科大太赫兹测量声致水面微幅波

相位补偿

在仿真的时候，水面干扰超过一定数值后，脉压之后，最强散射点出现在不同的距离门。
声致水面微幅波检测的原理是相位检测，重要的是相位提取这一过程。
最强散射点表示目标信息，然而在两个距离门过渡的阶段，提取相位，会出现相位跳变的情况。



由于相位提取出现错误，后续的小波分解会在跳变处分解出很多频率分量，也就是最后一张图片出现频率条纹的原因。
对其相位补偿，得到的结果：

卡尔曼滤波用来去除加性噪声，如雷达相位噪声。
相位解缠的方法：相邻采样点相位比较

文章内容

一：没有海面干扰，交替发送100Hz/150Hz

问题：发送频率大于350Hz以后，振幅非常微弱。

这和毫米波雷达chirp脉冲宽度有关系。脉宽越窄，可测量的频率范围越大。

实验

图6：单频500Hz，无干扰，

图7：

作为对比，22.5GHz能检测200Hz，但不能检测出500Hz声源。
结论：相同算法处理下，太赫兹雷达可以探测的频率范围更宽。

图8：

小波降噪后，太赫兹雷达最大可探测频率1100Hz，而微波最大可探测频率为400Hz
验证了：太赫兹在测量水面振动时精度更高。

平静水面依然存在微米量级的自然扰动，（没有说明微米量级扰动的数学模型）

图10：

0.2m深，1100Hz振幅97.6nm
0.5m深，500Hz振幅85.9nm
1m深，350Hz振幅61.3nm

图12：cm量级小尺度波浪进行干扰

图14：一维距离向内信号幅值差异大，信号连续性差。

即使滤除低频信号，原始水面振动包含大量噪声，在功率谱无法观察到待测振动的频率信号
使用本文算法，检测出的信号连续性差且强度较弱，信号中混有较多杂波

图15，二级海况，波浪起伏60cm，没有位于同一距离门。

0.5米深声源持续发射300Hz的振动信号，用本文算法可以提取有效信号。

半物理法

半物理法是指针对仿真研究内容，将被仿真对象系统的一部分以实物方式引入仿真回路；被仿真对象系统的其余部分以数学模型描述，并转化为仿真计算模型。
常用来验证控制系统方案的正确性和可行性，进行故障模式的仿真以及对个研制阶段的控制系统进行闭路动态验收试验。

半物理法：在原始振动位移加入与原始待测信号频率相同，幅值不同的正弦振动信号，
确定实际测量到300Hz的声源的幅度为130nm，而理论振幅为544nm

仿真

x = -5:0.1:5;
P = 180;    %点声压
f = 0:10:1000;      %声波频率
omega = 2*pi*f; %角频率:/s
rho = 1e3;
c = 1500;   %声波在水介质的传播速度
lamda = c./f;
k = 2*pi./lamda; %水面波动的波数:/m
mu = 1;%水介质粘度：Pa s
g = 9.8;%重力加速度:N/kg
sigma = 1;%表面张力系数:N/m
alpha = 4*mu*k.^2.*sqrt(k*g+sigma*k.^3/rho)./(rho*g+3*sigma*k.^2/rho);
ymax = 2*P./omega/rho/c;
figure(1);
subplot(1,2,1);plot(f,alpha);title('振幅衰减系数随频率变化');
subplot(1,2,2);plot(f,ymax);title('振幅随频率变化');

close all
x = -50:0.1:50;
N = length(x);
X = zeros(N,N);
for i=1:N
    for j=1:N
        X(i,j) = sqrt(x(i)^2+x(j)^2);
    end
end
P = 140;    %点声压
f = 130;      %声波频率
omega = 2*pi*f; %角频率:/s
rho = 1e3;
c = 1450;   %声波在水介质的传播速度
lamda = c/f;
k = 2*pi/lamda; %水面波动的波数:/m
mu = 1;%水介质粘度：Pa s
g = 9.8;%重力加速度:N/kg
sigma = 72.75e-3;%表面张力系数:N/m
alpha = 4*mu*k^2*sqrt(k*g+sigma*k^3/rho)/(rho*g+3*sigma*k^2/rho);
ymax = 2*P/omega/rho/c;
t = 0.5/c;
y = ymax * exp(-alpha*X) .* cos(k*X-omega*t);
figure;imagesc(x,x',real(y));
figure; plot(x,real(y(round(N/2),:)))
[X,Y] = meshgrid(x,x');
figure; mesh(X,Y,abs(y))

close all;
%{
不同频率，分时激励，水面振动持续2s
%}
f0 = [50:50:500,600:100:2000];
n1 = length(f0);
n2 = 5000;
%相当于采样率，大于频率的2倍，才能得到正负值
f1 = zeros(1,n1*n2*2);
T1 = n1*2*2;
t1 = linspace(0,T1,n1*n2*2);
for i = 0:n1-1
    f1(1,i*2*n2+1:i*2*n2+n2) = f0(i+1);
end
%--------------------------计算某时刻的振动y--------------------
rho = 1e3;
P = 160;
c = 1500;   %声波在水介质的传播速度
mu = 1;%水介质粘度：Pa s
g = 9.8;%重力加速度:N/kg
sigma = 1;%表面张力系数:N/m
y = zeros(1,n1*n2*2);
x0 = 0.2;
for i=1:n1
    f = f0(i);
    lamda = c/f;
    k = 2*pi/lamda; %水面波动的波数:/m
    omega = 2*pi*f; %角频率:/s
    alpha = 4*mu*k^2*sqrt(k*g+sigma*k^3/rho)/(rho*g+3*sigma*k^2/rho);
    ymax = 2*P/omega/rho/c;
    %每个单频产生的振动持续2s
    if (i-1)*n2*2+n2*2>n1*n2*2
        right_ind = n1*n2*2;
    else
        right_ind = (i-1)*n2*2+n2*2;
    end
    tmp = ymax * exp(-alpha*x0) .* cos(k*x0-omega*t1((i-1)*n2*2+1:right_ind));
    y((i-1)*n2*2+1:right_ind) = y((i-1)*n2*2+1:right_ind) + tmp;
end
figure;
subplot(1,2,1),plot(t1,f1/1e3,'.'),ylim([1,2500]/1e3);title('（1）单频');
subplot(1,2,2),plot(t1,real(y)*1e7);title('单频声源，雷达接收信号');ylabel('理论水面振动');