1.算法概述
载波跟踪环是传统独立式GPS接收机最脆弱的环节，针对弱信号环境下其比伪码跟踪环路更容易失锁的问题，给出一种基于矢量频率锁定环（vector-frequency lock loop，VFLL）的载波跟踪方法。给出VFLL理论推导及实现过程，并以最小二乘估计方法证明VFLL在载波跟踪性能上优于频率锁定环（frequency lock loop，FLL）。静止场景时9颗卫星实验结果显示，本文给出的方法能够实现14 dB/Hz微弱GPS信号的载波跟踪。 矢量型GPS信号跟踪算法（矢量延迟锁定环VDLL）
而在VDLL中，仅仅在DLL中对码跟踪进行改进，使其通过中心滤波器，而载波跟踪和传统的算法相同。所以，下面将重点对延迟锁定环进行改进，也就是你的课题的VDLL延迟锁定环。而VDFLL则是对码和载波分别进行改进。其基本结构如下所示：
即使用卡尔曼替代PLL，使用EKF替代DLL。这个是VDFLL，而VDLL则使用扩展卡尔曼滤波替代原DLL即可。
提出的VDLL(vector delay lock loop)方法直接估计用户位置信息,由于用户物理动态有限,与传统的独立通道码环相比,跟踪的维度和带宽都更小,所以该方法具有更强的鲁棒性.阐述了VDLL与传统独立通道码跟踪环的本质区别,建立了VDLL的非线性系统模型,推导了系统观测量与传输延迟估计误差的具体线性化关系,确立了观测误差方差矩阵的计算公式;然后对非线性系统模型进行线性化,给出了多卫星联合跟踪下用户位置更新的EKF(extended Kalman filte-ring)滤波算法.
2.仿真效果预览
matlab2022a仿真结果如下：
3.MATLAB部分代码预览 `...............................................................
%参数初始化
time = 1000*(10^(-3)); % 数据发送时间
time_unit = 20*(10^(-3)); % 数据跳变时间单位
time_cyc = 1*(10^(-3)); % 一个完整扩频码周期
fs = 5*(10^6);
nn = time_cyc*fs;
kk = (time/time_cyc)*nn;% 数据总采样点
F_if = 1.25*(10^6); % 载波中频
F_Carrier = 1575.42*(10^6); % L1波段载波频率
CA_freq = 1.023*(10^6); % CA码速率
%%
%%
%生成C/A码
PN = func_CAcodegen(svnum);
CA = [];
k = 5;
for n = 1:length(PN)
CA((1+k*(n-1)):k*n) = PN(n)*ones(1,k);%CA码扩展
end
tc = 1/(k*CA_freq);
loop_time = time/time_cyc;
%%
%%
%模拟产生测试信号源
[Signal_Source,Phase_signal,buffer_bit_data]=func_CreateSource(iniphcode,inifd,iniph,snr);
%%
%%
%在模拟之前，首先需要进行捕获
[fd_ac,f_ac_code,Corr_value] = func_acquire(Signal_Source);
figure
mesh(Corr_value);title('捕获结果');
%信号中断起始时间
break_start = 400;
break_end = 800;
P0 = [0 0
0 1];
P = [P0 zeros(2,2*(loop_time-1))];
T = 0.1;
LL = loop_time;
Y0 = [0;1];
data_out(:,1) = Y0; %Y的第一列等于Y0
A = [1 T
0 1];
B = [1/2*(T)^2 T]';
H = [1 0];
Q = (0.25)^2;
R = (0.25)^2;
X = zeros(1,loop_time);
%%
%%
%跟踪参数设置
tracking_parameter();
%进行跟踪
for i = 1:1:loop_time
i
%开始循环，每次循环去除一段数据
%开始循环，每次循环去除一段数据
%模拟信号突然中断
if i > break_start & i < break_end
Signal = 0.0001*rand(1,nn);%中断的时候，该段数据为随机的噪声干扰
else
Signal = Signal_Source((i-1)nn+1:inn);
end
%产生本地载波
t = [0:nn-1]*ts;
track_dopplar = fd_ac + track_freq_pll;
Track_Freq_Buffer = [Track_Freq_Buffer track_dopplar];
track_dopplar2 = [track_dopplar2 track_freq_pll];
Local_I = cos(2pi(F_if+track_dopplar)*t + Last_Phase);
Local_Q = sin(2pi(F_if+track_dopplar)*t + Last_Phase);
Iph = 2pi(F_if+track_dopplar)*t + Last_Phase;
Local_Ph_Buffer = [Local_Ph_Buffer Iph];
Last_Phase = Last_Phase + 2pi(F_if+track_dopplar)*time_cyc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Carrier_I = Local_I;%产生本地的载波
Carrier_Q = Local_Q;%产生本地的载波
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%产生本地相位码
%利用DLL的思路
%当前
ph_code_p = offside;
fd_code_p = track_dopplar;
CA_Code_p = func_CA(ph_code_p,fd_code_p,i);
lc_p = CA_Code_p.*Signal;
%早
ph_code_e = offside+diffoffside;
fd_code_e = track_dopplar;
CA_Code_e = func_CA(ph_code_e,fd_code_e,i);
lc_e = CA_Code_e.*Signal;
%迟
ph_code_l = offside-diffoffside;
fd_code_l = track_dopplar;
CA_Code_l = func_CA(ph_code_l,fd_code_l,i);
lc_l = CA_Code_l.*Signal;
%下变频
Local_P_I = lc_p.*Carrier_I;
Local_P_Q = lc_p.*Carrier_Q;
Local_E_I = lc_e.*Carrier_I;
Local_E_Q = lc_e.*Carrier_Q;
Local_L_I = lc_l.*Carrier_I;
Local_L_Q = lc_l.*Carrier_Q;
%积分运算
IPSum = sum(Local_P_I);
QPSum = sum(Local_P_Q);
IESum = sum(Local_E_I);
QESum = sum(Local_E_Q);
ILSum = sum(Local_L_I);
QLSum = sum(Local_L_Q);
%码相位环路控制
%鉴想器
theta_code = ((IESum.^2+QESum.^2)-(ILSum.^2+QLSum.^2))/((IESum.^2+QESum.^2)+(ILSum.^2+QLSum.^2));
I2_Q2(i) = IESum.^2 + QESum .^2;
%kalman
data(:,i) = theta_code;
j = (i-1)*2+1;
K = P(:,j:j+1)*H'inv(HP(:,j:j+1)*H'+R);%滤波增益
data_out(:,i) = data_out(:,i)+K*(data(1,1)-H*data_out(:,i)); %估计
data_out(:,i+1) = A*data_out(:,i); %预测
P(:,j:j+1) = (eye(2,2)-K*H)*P(:,j:j+1); %误差
P(:,j+2:j+3) = AP(:,j:j+1)A'+BQB'; %kalman滤波
CodeErr = data_out(1,i)/20;
%码环NCO
offside = offside_old+k1*CodeErr; %码NCO的输出
theta_code_old = theta_code; %将当前结果保存，用于下一个循环的码跟踪
CodeErr_old = CodeErr; %将当前结果保存，用于下一个循环的码跟踪
offside_old = offside; %将当前结果保存，用于下一个循环的码跟踪
Bk_DLL = [Bk_DLL theta_code]; %记录跟踪过程中的码环鉴想器的输出
Track_Code_Buffer = [Track_Code_Buffer offside]; %记录跟踪过程中的码环NCO的数出
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%载波跟踪
%载波跟踪
%载波跟踪
theta_pll = atan(QPSum/IPSum);
PLLinput = theta_pll/(2pitime_cyc);
Bk_PLL = [Bk_PLL theta_pll];
%LoopFilter
PLLoutput = func_CarLoopFilter(carrierw,carrierpllb/2,PLLinput,PLLinput_old,PLLoutput_old);
track_freq_pll = -PLLoutput;
PLLinput_old = PLLinput;
PLLoutput_old = PLLoutput;
adj_flag = track_dopplar - track_dopplar_old;
track_dopplar_old = track_dopplar;
adj_buffer = [adj_buffer adj_flag];
outdata = sign(real(IPSum));
ALL_Buffer_Data = [ALL_Buffer_Data outdata];
if adj_flag < 1
add = add+1;
else
add = 0;
end
if add >= 2
dem_flag = 1;
end
if dem_flag == 1
count_time = i;
count_buffer = [count_buffer count_time];
Buffer_Data = [Buffer_Data outdata];
end
end
%%
%%
%位同步与数据解调
Buffer_Data_out = func_bitssync(Buffer_Data,count_buffer);
l_i_d = time/time_unit;
l_o_d = length(Buffer_Data_out);
l_zeros = l_i_d - l_o_d;
Buffer_Data_out = [zeros(1,l_zeros) Buffer_Data_out];
%跟踪误差
l_dll = length(Track_Code_Buffer);
l_fll = length(Track_Freq_Buffer);
diata_dll = (Track_Code_Buffer(40:l_dll)-iniphcode);
break_start = 400;
break_end = 800;
%多普勒频率跟踪
figure;
plot(Track_Freq_Buffer);
xlabel('时间(ms)');
ylabel('多普勒频率跟踪结果(Hz)')
title('多普勒频率跟踪结果');
grid on
hold on
plot(break_start,min(Track_Freq_Buffer):0.01:max(Track_Freq_Buffer),'r-*','LineWidth',3);hold on
plot(break_end,min(Track_Freq_Buffer):0.01:max(Track_Freq_Buffer),'r-*','LineWidth',3);hold off
figure;
plot(I2_Q2(1:end),'LineWidth',3);
xlabel('时间(ms)');
ylabel('I^2+Q^2(Hz)')
title('I^2+Q^2');
grid on
hold on
plot(break_start,min(I2_Q2(1:end)):100000:max(I2_Q2(1:end)),'r-*','LineWidth',3);hold on
plot(break_end,min(I2_Q2(1:end)):100000:max(I2_Q2(1:end)),'r-*','LineWidth',3);hold off
%码相位跟踪
figure;
subplot(211);
plot(Track_Code_Buffer);
xlabel('时间(ms)');
ylabel('码相位跟踪结果');
title('码相位跟踪结果');
grid on
axis([0,length(Track_Code_Buffer),0.8iniphcode,1.2iniphcode]);
subplot(212);
plot(diata_dll);
xlabel('时间(ms)');
ylabel('码相位跟踪误差');
title('码相位跟踪误差');
grid on
axis([0,length(diata_dll),-10,10]);
hold on
plot(break_start,-10:0.1:10,'r-*','LineWidth',3);hold on
plot(break_end,-10:0.1:10,'r-*','LineWidth',3);hold off
01_048_m`
