基于优化包络相关的Loran-C（罗兰C）信号天地波识别算法及MATLAB仿真代码和实测信号处理

引言

Loran-C接收信号处理中，完成信号的捕获后需要进行信号的跟踪，Loran-C信号的跟踪按照信号处理流程分为天地波识别、周期识别和正向过零点跟踪三个步骤。Loran-C信号通过天波和地波两种方式到达接收点，导致接收到信号为地波和天波的合成信号，而实际到达时间测量是通过测量地波信号的到达时间实现的，这就需要从接收合成信号中提取地波信号，这个过程称为天地波识别，也是信号跟踪主要的环节。本文主要经典的基于优化包络相关的天地波识别方法原理，并给出算法仿真及结果。

优化包络相关天地波识别方法

基于相关系数的相关匹配算法是经典的匹配滤波技术。它是指通过计算参考信号与被测信号之间的相关系数来识别两个信号参数之间的相似度的过程。相似度可以反映信号的形状与时间序列之间的关系。在Loran-C信号的处理方法中，通常将信号的包络及其改进序列作为重要特征形状，而脉冲包络的微分是最为常用的优化包络之一。Loran-C信号的包络表达式为：

式中，为与天线电流峰值有关的归一化常数，单位为安培；为时间，单位为μs；为包络的时间起点与载波时间起点的差值，称为包周差（Envelope Cycle Difference，ECD），单位为μs。实际上，由于很小，远小于天地波相对延迟，在进行天地波识别中可以忽略，因此为了简化公式推导，令等于0，得到脉冲包络的微分为：

令等于0，得到为0或65μs，所以脉冲包络的上升时间为65μs，然而根据USCG发布的《Loran-C接收机最低性能标准》，实际中需要考虑的天波延迟为35μs至500μs。所以脉冲包络的上升时间大于最小的天波延迟，实际接收的地波信号峰值已经被天波干扰，不能作为特征形状。脉冲包络的二阶微分为：

令等于0，得到为19μs，所以脉冲包络一阶微分的上升时间为19μs，不会受到天波干扰的影响，并且按照包络微分峰值时间再向后推一周（Loran-C脉冲载波周期为10μs）即为第三周期，因此，只要找到峰值的位置，对应的时间加上11μs就是第三个正的过零点。因此在优化包络相关识别方法中，将标准脉冲包络的微分作为本地参考信号，其波形如图1所示：

图1 Loran-C脉冲包络微分

假设接收信号的脉冲包络微分为，而本地标准包络微分为，可以得到它们的相关函数并表示如下：

图2 Loran-C包络微分相关函数

如图2所示，当两个信号完全匹配时，相关函数具有一个特别突出的峰值。峰值的位置是时间匹配点，即峰值的TOA延迟。这种匹配滤波技术被认为是最好的、最有效的时延估计方法，它可以获得最大输出信噪比，与输入信号的波形和加性噪声的分布特性无关。因此，这些特性可以用来估计Loran-C天地波信号的时延。

其他性能与实测数据结果：

关键代码及仿真结果

完整代码见文末
运行代码需要安装matlab通信工具箱，运行时根据报错信息操作即可。

%Loran-C信号仿真
j=zeros(1,1);
x1=j; x2=j;
fs=1e6;
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-100+1:end) y(1:end-100)];%将原信号右移模拟天波信号
    r=y+swave;
    if SNR<23
        r1=awgn(r,SNR+18);
    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微妙
    if delay<36
        swave=-10^0.6*[y(end-delay+1:end) y(1:end-delay)];%将原信号右移模拟天波信号
    else
        swave=-10^0.6*[y(end-100+1:end) y(1:end-100)];%将原信号右移模拟天波信号
    end
     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;
load b.mat;
%fs=1e6;
r=b3;
me=abs(hilbert(r));%罗兰脉冲包络
yd=(diff(b))/max(diff(b));%y标准罗兰脉冲，yd标准罗兰脉冲微分
md=(diff(me))/max(diff(me));%实际接收罗兰脉冲微分
x=xcorr(md,yd);%微分相关
figure;plot(x);
title('实际信号测试结果');xlabel('延迟/us');ylabel('相关值');

%抗天波干扰性能测试天波
SCG=12;
m=1;
for delay=30:80;%天波延迟30~80微妙
    if delay<36
        swave=-10^0.6*[y(end-delay+1:end) y(1:end-delay)];%将原信号右移模拟天波信号
    else
        swave=-10^0.6*[y(end-100+1:end) y(1:end-100)];%将原信号右移模拟天波信号
    end
     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;
load b.mat;
%fs=1e6;
r=b3;
me=abs(hilbert(r));%罗兰脉冲包络
yd=(diff(b))/max(diff(b));%y标准罗兰脉冲，yd标准罗兰脉冲微分
md=(diff(me))/max(diff(me));%实际接收罗兰脉冲微分

x=xcorr(md,yd);%微分相关
figure;plot(x);
title('实际信号测试结果');xlabel('延迟/us');ylabel('相关值');

完整代码分为多个文件，数量较多，无法放在文章中。完整代码在公众号（沸腾的火锅资源号）中自取。