基于MEDLL的Loran-C天地波周期联合识别及信号跟踪方法及MATLAB仿真代码和实测信号处理

本文链接：https://blog.csdn.net/BoilingHotPot/article/details/142961275

引言

针对目前Loran-C接收机天地波识别方法的不足，提出了一种优化包络相关与MEDLL相结合的天地波识别方法，并与正峰值比检测算法结合，从而实现天地波周期的联合识别。并从优化包络峰值检测原理、MEDLL估计方法、正峰值比检测门限等方面描述了该方法。仿真分析了联合识别方法的抗噪性能和抗天波干扰性能，验证了天地波识别的有效性，最后将该方法应用到实际Loran-C信号的处理，进一步验证了该方法的天地波和周期识别性能。试验结果表明：该方法具有良好的抗天波干扰性能，当天波延迟大于20us时，该方法能够准确地进行地波和周期识别，满足Loran-C接收机大范围的使用要求。

天地波周期联合识别算法

（1）MEDLL天地波周期联合识别

实际接收到的罗兰C信号在这里插入图片描述为地波和天波的合成信号，由天地波参数构建合成信号表达式为：

式中，、为地波幅度和延迟，，（）为m次天波幅度和延迟。实测相关函数为：

其中，为罗兰C标准脉冲包络微分的自相关函数。因此，使估计均方误差最小的参数取值为：
在这里插入图片描述
由上式可知，参数间相互耦合，无法单独求出，需要通过迭代更新参数估计结果，实现流程如下:
1)令，，计算、，得到；
2）令，则，估计、，更新估计；
3），计算、，更新估计；
4）重复步骤2），3），直到估计均方误差收敛。
MEDLL（Multipath Estimation Delay Lock Loop）峰值检测是指基于多相关器原理，产生等间隔的标准罗兰C超前滞后包络微分（早迟码）分别与接收脉冲信号包络微分进行相关累积，从相关函数中提取地波和天波的幅度和延迟估计值，MEDLL峰值检测原理如下图所示：

在这里插入图片描述

图1 MEDLL峰值检测原理

（2）正峰值比

罗兰C脉冲各周期正向峰值在罗兰C包络峰值为1（归一化）时，都是确定的，如表1所示：

表1 Loran-C脉冲各周期归一化正峰值

零点(us)	10	20	30	40
归一化正峰值	0.0101	0.1860	0.4430	0.6795

各周期峰值与后一周期峰值比如表2所示：

表2 Loran-C脉冲周期正峰值比

零点(us)	10	20	30
归一化正峰值比	18.3790	2.3819	1.5338

因此，当比率接近 1.5538 时，相应的周期被认为是第三个周期。但在实际应用中，该值往往不等于 1.5338。一方面，在数字接收机中，信号是被采样的，因此不能保证采样点正好在峰值处。另一方面，由于噪声的影响，比值也会偏离理论值。在本文中，接收器的采样率为 1MHz，信噪比约为 18dB。经过长时间观测，确定第三个周期的峰值比在 1.53 和 1.80 之间。

（3）正向过零点估计

完成周期识别后，需要准确跟踪载波的第三周期正向过零点，由于实际中接收机的数据是经过采样后的信号，但在采样过程中，往往很难保证采样点恰好位于过零点处，因此在实际的信号处理中，一般采用极性检测的方法，即找到过零点附近的采样点，再通过插值的方法对过零点时间进行准确的估计，实际跟踪采样点的位置有以下三种情况：

在这里插入图片描述

图2 极性检测情况示意

在图2中，三个采样点的极性分别为“-++”，此时认为准确跟踪第三周正向过零点；在图2（b）中，三个采样点的极性分别为“+++”，没有负极性点，跟踪点左移一个采样点；在图2（c）中，三个采样点的极性分别为“–+”，虽然同时存在正负极性点，但为了统一测量标准，将跟踪点右移一个采样点，从而与（a）中情况保持一致，经过跟踪采样点的平移，得到两个相邻的正负极性点，然后在这两个采样点之间进行线性插值。假设前两个采样点的采样值为在这里插入图片描述和，采样点序号为和，则得到标准过零点的时刻为：

式中计算得到时刻的单位为：采样点，若需得到准确的时间，则用计算出的采样点时间乘采样间隔，即为绝对时刻。

关键代码及结果图

完整代码见文末

%% MEDLL天地波周期识别方法抗天波干扰仿真测试
clear
close all

%% 仿真生成标准Loran-C信号
y=zeros(1,1000);
j=zeros(1,1);
x1=j; x2=j;
fs=1e6;   %采样率1MHz
T=1e6/fs;
for i=0:T:300  %单位是us
    x1=((i)/65).^2;
    x2=exp(2*(0-(i)/65));
    x3=sin(0.2*pi*i);
    j=j+1;
    y(j)=x1*x2*x3;
end

%% 抗噪性能测试
for SNR=0:28
    swave=-10^0.6*[y(end-50+1:end) y(1:end-50)];%将原信号右移模拟天波信号
    r=y+swave;
    if SNR<14
        r1=awgn(r,SNR+30);
    else
        r1=awgn(r,50);
    end
    bze1(SNR+1)=TrackMainStation(fs,r1);
end
figure;plot(0:28,bze1);
title('抗噪性能测试结果');xlabel('SNR/dB');ylabel('周期识别输出/us');

%% 抗天波干扰性能测试天波
SCG=12;
m=1;
for delay=30:80;%天波延迟30~80微妙
     swave=-10^0.6*[y(end-delay+1:end) y(1:end-delay)];%将原信号右移模拟天波信号
     r2=y+swave;
     bze2(m)=TrackMainStation(fs,r2);
     m=m+1;
end
figure;plot(30:80,bze2);
title('抗天波干扰性能测试结果');xlabel('天波延迟/us');ylabel('周期识别输出/us');

%% 实际信号测试结果
load b3.mat;
r=b3;
load b.mat;
%fs=1e6;

me=abs(hilbert(r));%罗兰脉冲包络
yd=(diff(b))/max(diff(b));%y标准罗兰脉冲，yd标准罗兰脉冲微分
md=(diff(me))/max(diff(me));%实际接收罗兰脉冲微分

delta1=find(me==max(me))-find(b==max(b));  %接收信号包络峰值与标准包络的距离

%DLL峰值检测
delta2=fs*10e-6;%相关器间隔10us
Ec=delta1-delta2; %早码左移的距离
Lc=delta1+delta2;%迟码右移的距离
if Ec<0
   E=[b(abs(Ec)+1:end) zeros(1,abs(Ec))];
else if Ec>0
   E=[zeros(1,Ec) b(1:end-Ec)];     
    else E=b;
    end
end   %早码
if Lc<0
   L=[b(abs(Lc)+1:end) zeros(1,abs(Lc))];
else if Lc>0
   L=[zeros(1,Lc) b(1:end-Lc)];     
    else L=b;
    end
end   %迟码
et=xcorr(diff(me),diff(E))-xcorr(diff(me),diff(L));
%plot(et);
etd=diff(et);
right=find(et==max(et));
zz=length(et)/2;

R=etd(900:1099)/max(etd(900:1099));%实测
%R=R';
load dd.mat;
dd1=conversion(dd,1e6,fs);
R=[R zeros(1,length(dd1)-length(R))];
Rx=dd1;
Rx=Rx/max(Rx);%基准
%R1=zeros(1,length(R));
R2=zeros(1,length(R));
for i=1:40
    R1=R-R2;
    [a1,tao1,R1]=estimationR(R1,Rx);
    R2=R-R1;
    [a2,tao2,R2]=estimationR(R2,Rx);
    res=R-R1-R2;
    mres(i)=mean(res);
    stdres(i)=std(res);
end
figure;plot(R1);hold on;plot(R2);
title('实际信号测试结果');xlabel('延迟/us');ylabel('相关值');