北斗三号双频B1C/B2a MATLAB软件接收机捕获跟踪定位

简介这套面向北斗三号软件接收机学习与二次开发的 MATLAB 工程聚焦 B1C/B2a 双频信号处理适合具备 GNSS 基础、希望从理论走向实现的高校学生、科研人员与 SDR 开发者。压缩包共 87 个文件以 82 个 .m 脚本与函数为主体辅以 2 个 Markdown 说明、2 个 txt 和 1 份 PDF 文档整体 3.72MB代码与文档配套便于按模块阅读和复现。已有 1936 人学习说明其在北斗三号开源接收机方向具有一定参考价值。内容覆盖数据与导频通道的捕获跟踪、导航电文解码、伪距测量生成和位置计算支持 RF 二进制中频文件后处理并给出基于 NUT4NT 采样器实测 IF 信号的验证路径目录按通用功能、BDS-3 B1C/B2a 及文档分区适合课程设计、算法验证与接收机链路拆解。1. 双频 B1C/B2a 软件接收机跑通这套 MATLAB 工程能拿到什么一段冷冰冰的 IQ 采样文件跑完这套系统开源的 MATLAB 工程之后会变成每颗卫星的多普勒频移、码相位、伪距以及解算出来的终端坐标。它处理的是北斗三号的 B1C 和 B2a 两个民用信号B1C 在 1575.42 MHz码率 1.023 Mcps主码 10230 长、十个毫秒一个周期B2a 在 1176.45 MHz码率 10.23 Mcps主码同样 10230 长、但只有一毫秒。工程本身不含射频前端输入是已经下变频采好的二进制文件所以只要手上有一段录制的数据或者在电脑里造一段仿真数据整条捕获、跟踪、电文、定位链路就能跑起来。适合两类人一类是想弄明白 BOC 调制到底把捕获难在哪里的在校学生另一类是要拿双频观测量做电离层组合、或者在把算法搬到 C 和嵌入式之前先做验证的工程人员。2. 工程骨架参数面板、数据源与主流程怎么串起来2.1 文件清单与调用链MATLAB 工程最容易写成一锅粥任何时候从任何一个脚本点运行都能跑但跑出来的东西对不上。这套代码的骨架是刻意收敛过的只有一个参数入口只有一个主脚本其余全是纯函数靠入参出参通信不依赖工作区里的隐式变量。这种写法在system里调evalin调用别处的变量虽然方便但换个数据集就崩血泪经验就是不碰。路径作用改动频率main.m总入口读参数、分块、逐块送捕获和跟踪、收尾解算低settings_init.m全局参数面板返回一个 struct高gen_prn.m按信号体制生成 10230 长的本地主码低read_iq_block.m从交错 int16 文件里读一段复数样点低acq_fft.m并行码相位捕获多普勒串行扫描中track_channel.m单通道单块跟踪输出即时支路中nav_decode.m帧同步、电文解析、星历参数提取中pvt_solve.m伪距组装、最小二乘定位中调用链是单向的main.m先调settings_init.m拿到 cfg打开数据文件按cfg.msPerBlock一块一块读第一块交给acq_fft.m扫出所有可见卫星的码相位和多普勒初值之后每一块把上一块的跟踪状态带进去继续调track_channel.m。等某个通道的位同步锁定main.m才开始往nav_decode.m里灌积分结果。整个流程没有任何一个函数回头去改 cfg这样做的好处是参数一旦定死同样的数据文件在任何机器上跑出来的结果一致排查问题的时候变量只有一个。实际项目里最容易被忽略的是块与块之间的状态延续。track_channel.m必须把载波相位、码相位、环路滤波器里攒的那点量全部带出去再带回来否则每一块都相当于重新入锁伪距会出现周期性锯齿。这套代码把跟踪状态统一塞进一个trk结构体通道号做索引主循环里就是trk(ch) track_channel(block, acq(ch), cfg, band, trk(ch))这么一行看着平淡但少了它后面全是玄学。2.2 参数面板采样率、中频与相干积分长度的三角关系参数面板是整个工程的地基改错一个数后面的现象会离谱到你怀疑代码有 bug。下面这份是能直接用的起点每一项都建议先照抄再微调% settings_init.m —— 全工程唯一参数入口其余函数只读不改 function cfg settings_init() cfg.fs 40.92e6; % 采样率, 取 B2a 码率的 4 倍, 兼顾 B1C 的 BOC(6,1) 分量 cfg.fIF 0; % 0 表示零中频复数基带; 实数中频需在载波复现里改符号 cfg.dataType int16; % 每样点两路 int16, I 在前 Q 在后 cfg.msPerBlock 10; % 单块时长, 10 ms 正好一个 B1C 主码周期 cfg.b1c.fc 1575.42e6; % 载波频率 cfg.b1c.rc 1.023e6; % 码率 cfg.b1c.Nc 10230; % 主码长度, 10 ms 一周期 cfg.b2a.fc 1176.45e6; cfg.b2a.rc 10.23e6; cfg.b2a.Nc 10230; % 主码长度, 1 ms 一周期 % 捕获参数 cfg.acq.dopplerSpan 5000; % 搜索范围 ±5 kHz, 静态接收机够用 cfg.acq.dopplerStep 250; % 步长, 要小于 1/(2T) 才能不漏峰 cfg.acq.threshold 2.5; % 峰值/噪声均值, 低于此判为未捕获 % 跟踪参数 cfg.trk.spacing 0.5; % 早迟相关器间隔, 单位码片 cfg.trk.bnPll 18; % PLL 噪声带宽 Hz cfg.trk.bnDll 1.5; % DLL 噪声带宽 Hz cfg.trk.zeta 0.707; % 阻尼系数 end采样率这里有个必须说清的约束。B2a 码率 10.23 MHz主瓣宽度是码率的两倍也就是 20.46 MHz按奈奎斯特定理fs至少要到 20.46 MHz 才不混叠取 40.92 MHz 是四倍过采样好处是每个码片正好 4 个采样点早迟相关器的半码片间隔能落到整数样点上不用每块都做分数插值。B1C 那边更麻烦一点导频分量里除了 BOC(1,1) 还有 BOC(6,1)后者的主瓣要宽到 6.138 MHz 附近两边加起来带宽也得十几个兆所以同一个fs兼顾两个频点是合理的省得切换频点时重采样。多普勒搜索步长和相干积分长度是绑死的。频率搜索的栅格间隔要小于1/(2T)T是相干积分时长B2a 用 1 ms 的时候这个数就是 500 Hz。你按 500 Hz 走最坏情况卫星真实多普勒落在两个格点正中间相关峰损失大约 2 dB勉强能捕到但跟踪入锁会很吃力按 250 Hz 走损失压到 0.5 dB 以内代价是搜索时间翻倍。我一般先用 250 Hz 跑一遍看看有多少颗星长期跑再考虑用码相位并行加频率 FFT 的两维搜索压时间。参数改动后的直接后果fs调低到 20.46e6每码片只剩 2 个采样点早迟间隔 0.5 码片正好 1 样点DLL 抖动明显变大fIF填错捕获峰整体偏移表现为所有卫星多普勒都在同一方向偏了一大截msPerBlock超过次级码符号长度相干积分跨符号翻转相关值被抵消跟踪直接失锁dopplerStep大于 500 Hz部分卫星能捕到但峰很矮跟踪阶段频繁丢锁最后一行那个坑值得单独提一句。B2a 导频分量的次级码长度是 100 个符号符号率 1 kHz也就是 100 ms 才翻一次1 ms 的相干积分完全不受影响可一旦你为了提信噪比把积分拉到 10 ms 又没剥次级码就等于把十个可能变号的符号硬加在一起运气不好直接加到零。B1C 那边反过来主码周期就是 10 ms次级码 1800 长、18 秒才翻一次10 ms 相干积分是安全的。两个频点的安全积分长度不一样千万别用同一套参数硬套。2.3 复数基带数据的读取与分块数据源这一层最不起眼但翻车率不低。录制设备吐出来的通常是 I、Q 交错的 int16也就是文件里按I0 Q0 I1 Q1 ...排用fread的时候维度给错读出来的波形看着像噪声你会以为是前端坏了。function [iq, nRead] read_iq_block(fid, nSamp) % 从交错 int16 文件读一段复数样点 % nSamp: 期望读取的复数样点个数 raw fread(fid, [2, nSamp], int16double); % 按列填充: 第1行全是 I, 第2行全是 Q nRead size(raw, 2); if nRead 0 iq complex(zeros(0,1)); return; end iq complex(raw(1,:), raw(2,:)).; % 转成列向量, 与后续索引习惯一致 endfread给[2, nSamp]这个尺寸MATLAB 是按列优先填充的所以文件里前两个数I0 Q0落在第一列第三四个数落在第二列最后得到的raw第一行是完整的 I 路、第二行是完整的 Q 路。这个细节写反了不会报错只会让你对着一堆看不出相关峰的图发呆。转成列向量是为了跟后面(0:N-1).这种时间轴对齐行向量和列向量在 MATLAB 里做.*会隐式扩展成矩阵几百兆的内存瞬间吃满这个坑我踩过不止一次。分块读取的另一个作用是控制内存。40.92 MHz 采样率、int16 双通道一秒数据就是 160 MB 出头直接fread整个文件在普通笔记本上很容易把内存打爆。按 10 ms 一块读单块 40 万个复数样点占用不到 7 MB跑几百秒的数据也很从容。文件末尾那块通常读不满nRead会小于nSamp跟踪代码里所有跟长度相关的计算都得用nRead而不是nSamp否则最后一块会索引越界。3. 捕获FFT 并行码相位搜索怎么按住 BOC 的多峰3.1 载波剥离与循环相关捕获要解决的是二维搜索码相位一个维度多普勒一个维度。朴素做法是每个多普勒格点、每个码相位都做一次相关B2a 一个码周期 10230 个码片、40.92 MHz 采样下四万个样点一轮下来计算量能让人等到睡着。FFT 循环相关的意义就是把码相位这一维一次性算完。function acq acq_fft(iq, prn, cfg, band) % 单颗卫星: 多普勒串行扫描 码相位 FFT 并行搜索 L numel(iq); t (0:L-1). / cfg.fs; code gen_prn(prn, band, cfg); % 本地主码, ±1 序列, 长度 Nc code code(1:min(numel(code), L)); % 与数据等长, 不足则循环补齐 codeF conj(fft(code, L)); % 本地码共轭谱, 循环外只算一次 bestPeak 0; bestCp 0; bestFd 0; bestCorr []; for fd -cfg.acq.dopplerSpan : cfg.acq.dopplerStep : cfg.acq.dopplerSpan carr exp(-1j*2*pi*(cfg.fIF fd)*t); % 复本振, 剥离残余载波 x iq(:) .* carr; % 混频到基带 R ifft(fft(x, L) .* codeF); % 循环相关, 一次给出全部码相位 P abs(R); [pk, idx] max(P); if pk bestPeak bestPeak pk; bestCp idx-1; bestFd fd; bestCorr P; end end acq.peak bestPeak; acq.codePhase bestCp; acq.doppler bestFd; acq.corr bestCorr; end核心就三行混频、正变换、乘共轭谱再反变换。ifft(fft(x).*conj(fft(code)))这个式子在数学上等价于把本地码在时域上循环移位后逐一做相关只不过全部用频域乘法一次算完复杂度从O(L·Nc)降到O(L log L)。L取一个码周期的整数倍最省事B2a 取 1 ms 也就是 40920 个样点B1C 取 10 ms 也就是 409200 个样点。后者明显更慢实际调试时可以先按 1 ms 粗捕捕到了再用 10 ms 精化码相位。gen_prn.m内部按信号体制走不同分支B2a 用的是分层码结构由一段基础码和扩展序列拼出来B1C 用的是两个 m 序列截断组合的 Weil 码导频分量还要叠一层子码。这些码表第一次生成之后建议存成mat或者直接persistent缓存别每次调用都重算一颗星 10230 个码片、几十颗星跑下来光码生成就能占掉可观的时间。3.2 BOC(1,1) 主副峰判据与双峰确认普通 BPSK 信号的相关峰是个干净的三角形取最大值就完事。BOC 调制的自相关函数因为副载波的存在被劈成了多个峰中心主峰两侧各约半码片处还有次峰幅度大概只有主峰的一半。麻烦在于跟踪环的早迟间隔只有半码片如果捕获阶段锁到了次峰上跟踪环会非常稳定地锁在一个错误的码相位上伪距直接偏出去上百米而且环路一切正常你看不出任何异常——这是最阴的一种翻车。% 在最大峰附近找次峰, 用峰间距判断是否落在 BOC 副峰上 [pks, locs] findpeaks(acq.corr, SortStr, descend, NPeaks, 4); if numel(locs) 2 dSample abs(locs(1) - locs(2)); % 主次峰间隔, 单位样点 dChip dSample * band.rc / cfg.fs; % 换算成码片 ratio pks(2) / pks(1); % 次峰与主峰幅度比 % BOC(1,1) 的次峰大致落在半码片外, 幅度约为主峰的一半 if abs(dChip - 0.5) 0.2 ratio 0.35 acq.warn sub-peak-nearby; % 标记可疑, 交给跟踪段双峰复核 end end这段判据的作用不是当场做决定而是打一个标记。真实场景里主峰和次峰的采样点间隔会因为你fs和码率的比值而偏移dChip算出来可能是 0.45 也可能是 0.55容差给到 ±0.2 码片比较稳。幅度比同理受噪声和残余多普勒影响会浮动0.35 这个下限是宁可错杀不可放过。标记之后交给跟踪段处理常见做法是双峰确认让跟踪环先在捕获给出的码相位上入锁稳定几百毫秒后人为把码相位拨动一个副峰间距看相关值是不是掉下去又回来。如果拨动之后相关值反而升高说明之前锁的就是副峰这时候把码相位整体搬一个间隔就行。另一种更省事的思路是在捕获阶段直接对相关函数做副峰消除把 BOC(1,1) 的自相关跟它的移位版本相减主峰保留、副峰压制代价是主峰也会变窄一点DLL 的牵引范围变小。两种做法各有取舍工程上我更倾向双峰确认因为它不改动捕获的相关函数形状出问题时现象好解释。3.3 门限、多星并行与搜索步长捕获输出的判决量用「峰值比噪声均值」比用绝对阈值靠谱得多因为不同录制设备的增益、不同天线的噪声底完全不一样写死一个数换台机器就废。判决量建议门限说明峰值 / 均值 2.5均值取整张相关图的平均含噪声不含峰峰值 / 次大峰 1.5防止多颗星互相干扰时误判有效卫星数≥ 4少于四颗直接放弃本次解算多星并行在 MATLAB 里没必要上parfor一颗星一轮捕获也就几百毫秒三十多颗星顺序跑完十几秒调试阶段反而更好定位问题。真要提速优先做的是把多普勒扫描改成频域并行对混频后的数据做 FFT在频域上搬移一次反变换就能覆盖多个多普勒格点这个改造在第 6 章还会提到。搜索范围dopplerSpan也不是越大越好。地面静态接收机卫星最大多普勒在 B1 频段大概 ±5 kHz 量级范围开到 ±10 kHz 只是白白多花一倍时间还会把一些杂散峰纳进候选。真要做高动态场景多普勒变化率也得进参数表否则跟踪环一开始就跟不上。4. 跟踪DLL/PLL 鉴别器与环路滤波器的落地写法4.1 早即时晚相关器与积分清零跟踪段每一块做的事情说穿了就是「用上一轮的估计值生成本地复本跟输入混一遍看误差有多大再把估计值往正确方向推一点」。三个支路早、即时、晚间隔半码片对应码相位的三个采样点。function trk track_channel(iqBlock, cfg, band, trk) N numel(iqBlock); n (0:N-1).; % 1) 载波复现: 用上一轮估计的多普勒与相位 carr exp(-1j*(2*pi*(cfg.fIF trk.fd)*n/cfg.fs trk.phi)); base iqBlock(:) .* carr; % 2) 三路码相关, 间隔 spacing 码片 corr zeros(1,3); for k 1:3 off (k-2) * trk.spacing; % -d, 0, d code gen_prn_shift(band.prn, band, cfg, trk.codePhase off); corr(k) sum(base .* code); % 相干积分, 一个码周期 end E corr(1); P corr(2); L corr(3); % 3) 鉴别器 trk.errCode (abs(E)-abs(L)) / (abs(E)abs(L)eps); % 归一化早迟功率 trk.errCarr atan2(imag(P), real(P)); % Costas 环, 对数据翻转不敏感 % 4) 环路滤波 NCO 更新 trk.fd trk.fd trk.k2 * trk.errCarr; trk.phi trk.phi trk.k1 * trk.errCarr trk.fd*2*pi*N/cfg.fs; trk.codePhase trk.codePhase trk.kc * trk.errCode * cfg.fs/band.rc; end几个参数的单位必须较真。trk.fd是 Hztrk.phi是弧度trk.codePhase是码片。gen_prn_shift内部要做的是按fs把 10230 长的码表上采样成四万个样点的序列然后按整数偏移加分数偏移取值分数部分用线性插值这一步偷懒直接用round取整DLL 的稳态误差会明显变大。eps加在分母上是防除零静态场景下噪声再低也可能出现E和L同时为零的瞬间。band.prn存在trk里而不是每次从 cfg 里翻是因为跟踪函数要足够纯粹给同样的输入和状态输出必须一致这样出问题可以单独把某一块数据抠出来复现。4.2 鉴别器公式与环路系数环路系数不是拍脑袋定的二阶环的系数由噪声带宽、阻尼系数和积分时间三个量决定公式在 GNSS 教材里翻来覆去就那一个但代入的时候容易把T用错。T cfg.msPerBlock * 1e-3; % 积分时间, 秒 Bn cfg.trk.bnPll; % PLL 噪声带宽, Hz zeta cfg.trk.zeta; theta Bn*T / (zeta 1/(4*zeta)); % 归一化环路增益 trk.k1 4*zeta*theta / (1 2*zeta*theta theta^2); % 比例支路 trk.k2 4*theta^2 / (1 2*zeta*theta theta^2); % 积分支路 % 一阶 DLL, 系数按指数衰减形式给, 比线性近似稳 trk.kc 1 - exp(-4 * cfg.trk.bnDll * T);T这里必须是实际参与相干积分的时长不是主循环的块长。如果 B2a 用了 1 ms 积分、主循环 10 ms 一块那T就是 1e-3写成 10e-3 会让环路带宽等效变窄跟踪看起来更稳实际上对动态响应变迟钝高动态场景直接丢锁。这个错我早期犯过现象是静止场景效果好得离谱一上车就全崩。噪声带宽的选取是典型折中。PLL 带宽压到 10 Hz 以下载波相位噪声小了载波平滑伪距的质量高但接收机稍微一动环路就跟不上多普勒变化率容易失锁放到 25 Hz 以上又太抖伪距噪声变大。18 Hz 配 0.707 阻尼静态和低速场景都能凑合动态大的场景建议做成带宽自适应用鉴别器输出的绝对值做滑动平均超过阈值就临时把带宽提上去。DLL 带宽比 PLL 低一个数量级是常识1.5 Hz 对静态足够了动态再往上加。4.3 位同步、子码剥离与相干积分长度的取舍载波和码都锁住之后还得知道电文比特从哪开始。B2a 一毫秒一个主码周期二十个周期凑一个比特比特边界落在哪二十个里面得靠统计。% 位同步: 在即时支路同相分量上找符号跳变最集中的位置 spb 20; % 每比特对应的积分块数, B2a 取 20 hist zeros(1, spb); for k 1:numel(Ip)-1 if sign(Ip(k1)) ~ sign(Ip(k)) % 相邻积分块符号是否翻转 hist(mod(k, spb)1) hist(mod(k, spb)1) 1; end end [~, bitEdge] max(hist); % 跳变最密集的位置即比特边界这个方法靠的是电文本身有足够多的 0/1 跳变数据量不够的时候直方图会很平max取到哪都差不多这时候别急着判定锁定继续攒。判据可以加一条最大计数要显著高于平均计数比如超过均值两倍以上。B1C 那边主码周期就是 10 ms比特率和次级码的关系不一样spb这个数得单独配直接照搬 B2a 的值就是白跑。相干积分长度是这章最后一个要权衡的点。积分越长信噪比越高但受限于符号翻转B2a 数据分量 1 ms 就可能翻1 ms 是安全上限导频分量有 100 ms 周期的次级码剥掉之后能拉到 20 ms 甚至更长信噪比能提十几 dB。工程上的常见做法是双分量都用、数据分量走 1 ms 保证鲁棒导频分量在次级码已知的情况下走长积分提精度最后把两路的载波相位估计融合。这套代码里预留了分量开关但默认只跑一路改造起来不算费劲。5. 踩坑排查捕获到了跟踪不住、伪距跳变这几件事5.1 捕获与跟踪段的翻车记录现象一捕获图上有明显峰跟踪环却始终不入锁鉴别器输出在零附近乱跳。原因多半是多普勒格点没对准。按 500 Hz 步长扫出来的峰最坏情况跟真实频率差 250 Hz1 ms 相干积分下这个偏差带来的相关损失还算可接受但一进跟踪环PLL 的牵引范围只有几十赫兹根本拉不回来。解决办法是把dopplerStep降到 250 Hz 甚至 125 Hz 重捕一次或者干脆在捕获结果上做抛物线插值用峰值和左右两个格点的幅度估出亚格点的真实频率代价只是三行代码。现象二锁到副峰上伪距稳定地偏出去一段且环路指标一切正常。这是 BOC 调制独有的坑。判断依据是伪距残差在整段数据里保持一个常值偏移而且这个偏移换算成距离大致是几十米量级。确认方法是手动把码相位拨动半个码片看相关值是不是反而上升。解决就是第 3.2 节说的双峰确认把码相位整体搬一个副峰间隔再重新入锁。血泪经验是这种错误在小样本数据上极难发现必须跑足够长的数据看伪距随时间的变化曲线平不平。现象三跟踪稳了几秒之后突然失锁然后重新捕获又稳几秒循环往复。先把msPerBlock和实际积分长度对一遍。B2a 数据分量用 10 ms 相干积分又没剥次级码正好撞上次级码翻转的位置就会被抵消一次表现出来就是周期性失锁。另一个可能的原因是载波环路的相位累加用了single精度几十秒之后相位误差累积到不可忽略换成double立刻就好。这两种原因的现象很像排查顺序是先看积分长度再看数据类型前者改参数就能验证。现象四相关峰的峰值随数据块推进缓慢下降最后掉到门限以下。码多普勒没跟。卫星相对运动带来的多普勒不止作用在载波上码率也会按fd·rc/fc的比例缩放B1C 的码率只有 1.023 MHz比值小忽视它问题不大B2a 码率本文还有配套的精品资源点击获取