BDS-3 B1C/B2a双频SDR接收机MATLAB实现：从捕获跟踪到PVT

简介面向北斗三号信号处理学习与科研实践的 MATLAB 开源工程聚焦 B1C、B2a 频点的软件接收机实现适合具备一定 GNSS 基础、希望打通从中频数据到定位结果完整链路的高年级本科生、研究生与工程人员使用。压缩包共 87 个文件以 82 个 .m 脚本为核心辅以 md 说明、txt 许可与元数据、1 份 pdf 文档整体约 3.72MB代码按初始化、参数配置、捕获、跟踪、导航电文解码、伪距测量与 PVT 解算等模块组织并包含通用功能与 B1C/B2a 分目录实现。已有 1936 人学习说明其在同类开源接收机资料中具有一定参考价值。读者可据此复现数据与导频双通道跟踪流程结合 IF 数据集及元数据完成离线后处理也可参考文档理解各 SDR 接收机之间的共性设计与测试思路用于课程实验、算法验证和二次开发。1. B1C/B2a 双频 SDR 接收机一段中频数据能跑到定位吗去年帮同事调试一段 B1C 中频采集数据用几十行 MATLAB 代码就找到了相关峰他以为剩下的只是按老办法堆代码。结果一上手就卡住了B1C 是 BOC 调制主码周期 10 ms导频和数据分量还分开捕获那套老参数直接套用会失效B2a 码率 10.23 Mcps采样率一低相关峰就会被削平。这套 BDS-3 B1C/B2a SDR 接收机的开源 MATLAB 实现把从 IQ 采样到 PVT 解算的整条链路完整写了出来适合两类人一类是想弄明白双频接收机内部到底怎么运转的学生和刚入行的工程师另一类是手里有采集数据、想拿它当基线做二次开发的人。2. 工程结构与信号体制B1C/B2a 双频链路全景2.1 从 IQ 文件到定位结果的六个环节这套代码的骨架是一条单线程、文件驱动的处理链没有实时性要求也不依赖任何硬件驱动喂一段复基带采样就能从头跑到尾。整条链路拆成六段读采样文件并把中频搬移到零中频、按 PRN 做二维捕获、对捕获到的卫星开跟踪通道、比特同步与帧同步、电文解析出星历和钟差、观测量组帧后做最小二乘定位。每一段都是一个独立函数段与段之间只通过结构体传递这样改其中一段不会牵连其他部分。前两段最容易被轻视。采样文件通常是 int16 交织的 I/Q 或者 complex float读进来第一件事是确认字节序和交织方式第二件事是确认采样率与标称值是否一致。很多采集设备标称 8.184 MHz实际晶振偏几十 ppm落到 B1C 的 10 ms 主码上就是几十纳秒的漂移短时间看不出来跑几分钟跟踪就慢慢滑掉。代码里把采样率、中频、信号带宽、码率全部收在一个rx_config.m里就是为了让这种系统性偏差只有一个源头可以改。捕获和跟踪的边界也要划清楚。捕获只负责给出「这颗星大概在哪个码相位、多普勒大概多少」精度不需要很高码相位误差半个码片、频率误差几十赫兹都能被跟踪环拉回来。跟踪才是把精度做上去的地方它每个积分周期输出一次载波频率、码相位和 I/Q 累加值后面的比特同步、电文解析、伪距生成全部建立在这串观测量上。把这两段的职责混在一起写是新手最容易犯的结构性错误。最后一段是观测量到 PVT。伪距的物理含义是「接收时刻减去发射时刻乘光速」发射时刻要从电文的参考时间加上码相位折算接收时刻来自采样点索引两个时间必须落在同一个时间系统里。代码里用了一个t_rx数组统一管理每个采样点对应的接收时刻所有通道共用避免每个通道各自算一套时间导致伪距之间互相不一致。2.2 B1C 与 B2a 的信号参数对比双频接收机不是把两个单频接收机简单堆在一起两个频点的信号体制差别很大代码里必须分开处理。下面这张表是我按公开的接口文件整理的关键参数写代码时对着它改常数基本不会错项目B1CB2a中心频率1575.42 MHz1176.45 MHz调制方式数据分量 BOC(1,1)导频分量 BOC(6,1)合成为 QMBOC(6,1,4/33)BPSK(10)主码速率1.023 Mcps10.23 Mcps主码长度1023010230主码周期10 ms1 ms分量构成数据分量 导频分量数据分量 导频分量导航电文B-CNAV1B-CNAV2典型采样率8.184 MHz 及以上40.92 MHz 及以上B1C 的难点在调制。BOC(1,1) 的功率谱主瓣被搬到了 ±1.023 MHzBOC(6,1) 更远合成之后的 QMBOC 相关峰比 BPSK 窄得多主峰宽度只有 BPSK 的一半左右。窄主峰意味着测距精度更高代价是对码相位误差更敏感跟踪环的 DLL 间隔如果还按 BPSK 的 1 个码片去设鉴别器曲线会直接跑出线性区。代码里超前减滞后间隔默认取 0.5 个码片B1C 通道可以再压到 0.3 左右。B2a 的难点在采样率。码率 10.23 Mcps主瓣宽度 ±10.23 MHz采样率低于 40.92 MHz 时高次谐波被砍掉相关峰展宽、峰值下降捕获的归一化门限要相应调低否则明明有信号却判决不出来。另外 B2a 主码周期只有 1 ms相干积分时长最多堆到 1 ms 的整数倍想靠长时间积分换信噪比就得走导频分量加次级码那条路。两个频点的电文结构也不一样B-CNAV1 和 B-CNAV2 的编码方式、帧长、纠错方式都不同解析函数必须分开写。这里不要图省事写一个通用解析器两套电文的字段位置和编码规则差异太大强行合并只会让调试变成玄学。2.3 目录结构与各模块职责代码按功能分目录一个目录只干一件事。常见做法是下面这种组织方式我按它拆过几个类似的接收机改起来最顺手BDS3-SDR/ ├── main_acq.m % 捕获入口输出 acqResult.mat ├── main_track.m % 跟踪入口读捕获结果输出 trkResult.mat ├── main_pvt.m % 定位入口读跟踪结果输出 PVT ├── config/ │ └── rx_config.m % 采样率、中频、带宽、环路带宽、门限 ├── acq/ │ ├── acquire_b1c.m │ ├── acquire_b2a.m │ └── gen_local_code.m % 主码 次级码生成 ├── trk/ │ ├── track_channel.m % 单通道跟踪主循环 │ ├── dll_disc.m │ └── pll_disc.m ├── nav/ │ ├── bcnav1_decode.m │ ├── bcnav2_decode.m │ └── eph2pos.m % 星历转卫星位置 ├── pvt/ │ └── ls_pvt.m ├── data/ % 采样文件放这里 └── utils/ % FFT、坐标转换、绘图三个main_*.m各自独立可跑中间结果存成.mat好处是调跟踪的时候不用每次重跑捕获。捕获一段 100 ms 的数据在普通笔记本上要几十秒B2a 因为采样率高还要更久每次改跟踪参数都重跑捕获一天下来大部分时间都耗在等结果上。gen_local_code.m是共用模块两个频点的主码生成都从这里出。B1C 和 B2a 的主码都是 10230 长但生成多项式和初相不同函数入参里带一个频点标识内部走不同分支。次级码单独放在这个文件里因为导频分量的次级码周期长跟踪阶段要用它做电文剥离和主码不是一回事。config/rx_config.m建议只放常数不要放逻辑。我见过有人把环路带宽的自适应切换写在配置里结果调试时改一个参数要翻三个文件才能确认最终生效的是哪个值血泪经验。3. 捕获模块拆解并行码相位搜索的 MATLAB 实现3.1 捕获主函数逐段讲捕获做的是二维搜索码相位一维多普勒一维。最直接的做法是每个多普勒频点上都做一次时域相关计算量随多普勒网格数线性增长B2a 这种 40 MHz 采样率的数据根本跑不动。代码里用的是频域并行码相位搜索把本地码的 FFT 只算一次每个频点只做一次 FFT 和一次 IFFT把码相位那一维的搜索全部压进频域乘法里。function acq acquire_b1c(x, fs, prn, cfg) % x : 复基带采样列向量 % fs : 采样率 Hz % prn : 卫星 PRN 号 % cfg : 捕获配置多普勒范围、步进、相干积分时长、判决门限 codeRate 1.023e6; % B1C 主码速率 Nc 10230; % 主码长度 % 1) 生成本地主码并按采样率重采样到与信号同一时间栅格 code gen_local_code(prn, B1C); % 1 x 10230取值 ±1 splPerCode round(fs * Nc / codeRate); % 每个主码周期的采样点数 tLocal (0 : splPerCode-1) / fs; localCode interp1((0:Nc-1)/codeRate, code, mod(tLocal, Nc/codeRate), linear); localCode localCode(:); % 2) 相干积分长度B1C 导频分量可做 10 ms 整数倍数据分量压到 1 ms 量级 Ncoh round(fs * cfg.cohMs * 1e-3); nBlk floor(length(x) / Ncoh); % 3) 多普勒搜索网格 dopGrid cfg.fMin : cfg.fStep : cfg.fMax; % 典型 -6000:250:6000 % 4) 本地码 FFT 只算一次长度取线性卷积不混叠的下界 L Ncoh splPerCode - 1; codeF fft(localCode, L); best struct(peak, 0, codePhase, 0, doppler, 0, metric, []); for k 1:nBlk seg x((k-1)*Ncoh (1:Ncoh)); for d 1:numel(dopGrid) % 载波剥离把多普勒假设搬掉 n (0:Ncoh-1).; baseband seg .* exp(-1j*2*pi*dopGrid(d)*n/fs); % 频域并行码相位搜索 corr ifft(fft(baseband, L) .* conj(codeF)); corr corr(1:splPerCode); % 只取一个主码周期内的相关值 [v, idx] max(abs(corr)); if v best.peak best.peak v; best.codePhase idx - 1; % 单位采样点 best.doppler dopGrid(d); best.metric abs(corr); end end end acq best; acq.prn prn; acq.splPerCode splPerCode; acq.detRatio best.peak / mean(abs(best.metric)); acq.detected acq.detRatio cfg.thrRatio; end代码里有几个点值得展开说。第一本地码重采样用的是interp1线性插值对 ±1 的方波码来说插值会引入一点幅度损失但比起码相位量化的误差这点损失可以接受。如果想做得更干净可以按采样率直接生成码表省掉插值。第二FFT 长度L取Ncoh splPerCode - 1保证循环卷积不产生混叠取完corr的前splPerCode点才是有效的线性相关结果这段长度算错是新手最常见的翻车点相关峰会莫名其妙地分裂成两个。第三多普勒网格是外层循环、FFT 是内层这样本地码 FFT 可以复用但如果多普勒范围很大可以把 IFFT 那一维也并行化用pagefun或者直接循环展开。detRatio用的是峰值比均值不是峰值比次峰。峰值比均值对单颗强星效果好但在强信号压制弱信号的时候均值会被强星拉高弱星就漏了。更稳的做法是峰值比次峰把主峰附近 ±2 个码片、±1 个频点以外的区域取最大值当次峰代码里可以加一个掩膜。这一步多十几行但能明显降低虚警。3.2 多普勒网格、相干积分与门限怎么定多普勒范围取决于应用场景。地面静态接收B1C 的最大多普勒在 ±5 kHz 量级B2a 因为频率低一些量级相近但符号可能不同车载场景要放到 ±10 kHz机载还要更大。宁可在捕获阶段把范围给宽一点也不要在跟踪阶段因为初值偏太远拉不回来。代价是网格点变多计算时间线性增长。频率步进和相干积分时长是一对。二阶环的频率牵引范围大约是1/(2*Tcoh)的量级10 ms 相干积分对应 50 Hz 左右的半带宽步进取 250 Hz 时相邻网格之间会有重叠不会漏掉真正的峰。如果为了省时间把步进取到 1 kHz牵引范围覆盖不住跟踪环大概率失锁。我一般按fStep 1/(2*Tcoh)先算一遍再往上放一档取整250 Hz 和 500 Hz 是最常用的两个值。相干积分时长受电文符号率限制。B1C 导频分量没有电文翻转可以把相干积分堆到 10 ms 甚至 20 msB1C 数据分量和 B2a 数据分量有电文符号符号率分别是 1000 sps 和 500 sps 量级相干积分超过符号长度就会因为符号跳变而能量对消。所以捕获阶段如果拿不准用的是哪个分量保守取 1 ms 相干积分用非相干累加补信噪比累加次数 10 到 20 次。判决门限thrRatio没有普适值和积分时长、采样率、噪声带宽都有关。经验值是 2.5 到 3.5 之间先用一段确认有信号的数据把detRatio打出来看真峰和噪声底的实际差距再定门限。B2a 因为采样率高、噪声带宽大同样信噪比下detRatio会比 B1C 低一截门限要相应下调。3.3 捕获结果怎么验证捕获结果不能只看detected是 true 就放心。最直接的办法是把acq.metric画成二维图横轴码相位、纵轴多普勒看主峰是不是一个干净的尖峰。如果主峰旁边有明显旁瓣或者主峰被削平多半是采样率不够或者本地码重采样出了问题。第二个验证方法是对同一时刻的多颗卫星做码相位对比。同一段数据里捕获到的几颗星码相位折算成时间之后两两之差应该和它们到接收机的几何距离差对得上量级在毫秒以内。差得离谱说明码相位到时间的换算系数写错了常见错误是把splPerCode当成码周期用了。第三个方法是反算伪距率。捕获给出的多普勒是载波多普勒乘上光速除以载波频率就是伪距率和卫星运动引起的距离变化率应该在同一个量级。如果某颗星的伪距率符号反了说明载波剥离时用了j而不是-j这个错误在跟踪阶段会表现为载波频率一直往一个方向漂。4. 跟踪与电文解调从环路鉴别器到比特流4.1 跟踪主循环的完整代码跟踪是整条链路里最容易出问题的一段。捕获给的码相位和多普勒只是粗值跟踪要在每个积分周期里把这个粗值收紧同时输出载波相位和码相位的小数部分供后面的伪距计算使用。代码按标准的二阶环写DLL 和 PLL 各一路共用同一个积分周期。function trk track_channel(x, fs, acq, cfg) % 单通道跟踪输入已捕获的码相位与多普勒输出每积分周期的观测量 codeRate 1.023e6; Nc 10230; splPerCode acq.splPerCode; code gen_local_code(acq.prn, B1C); T splPerCode / fs; % 积分周期B1C 为 10 ms [w_dll, a_dll] loop_gain(cfg.Bn_dll, T, 0.707); [w_pll, a_pll] loop_gain(cfg.Bn_pll, T, 0.707); codePhase acq.codePhase; % 单位采样点可带小数 carrFreq acq.doppler; % Hz carrPhase 0; nIter floor((length(x) - 3*splPerCode) / splPerCode); trk struct(prn, acq.prn, I_P, zeros(nIter,1), Q_P, zeros(nIter,1), ... carrFreq, zeros(nIter,1), codePhase, zeros(nIter,1), ... CN0, zeros(nIter,1)); for k 1:nIter % ---- 1) 本地三路码超前/即时/滞后间隔 0.5 码片 ---- offE -0.5 * fs / codeRate; offL 0.5 * fs / codeRate; idxE mod(round(codePhase offE) (0:splPerCode-1), Nc) 1; idxP mod(round(codePhase) (0:splPerCode-1), Nc) 1; idxL mod(round(codePhase offL) (0:splPerCode-1), Nc) 1; cE code(idxE); cP code(idxP); cL code(idxL); % ---- 2) 载波剥离与相关累加 ---- n (0:splPerCode-1).; localCarr exp(1j*(2*pi*carrFreq*n/fs carrPhase)); seg x((k-1)*splPerCode (1:splPerCode)); yE sum(seg .* localCarr .* cE.); yP sum(seg .* localCarr .* cP.); yL sum(seg .* localCarr .* cL.); I_E real(yE); Q_E imag(yE); I_P real(yP); Q_P imag(yP); I_L real(yL); Q_L imag(yL); % ---- 3) 鉴别器 ---- pE I_E^2 Q_E^2; pL I_L^2 Q_L^2; dCode (pE - pL) / max(pE pL, eps); % 归一化超前减滞后功率 dCarr atan2(Q_P, I_P); % 四象限相位鉴别器 % ---- 4) 环路滤波与 NCO 更新 ---- codeFreq codeRate w_dll * dCode; carrFreq carrFreq w_pll * dCarr; carrPhase mod(carrPhase 2*pi*carrFreq*T, 2*pi); codePhase mod(codePhase codeFreq * T, Nc); % ---- 5) 记录观测量与载噪比 ---- trk.I_P(k) I_P; trk.Q_P(k) Q_P; trk.carrFreq(k) carrFreq; trk.codePhase(k) codePhase; trk.CN0(k) 10*log10(max(I_P^2 Q_P^2, eps) / ...); end endloop_gain按二阶环的标准形式写自然角频率wn 2*Bn/(zeta 1/(4*zeta))两个滤波系数分别是a 4*zeta*wn*T/(12*zeta*wn*T(wn*T)^2)和与之配套的另一个系数。这里不需要手推只要保证传入的带宽Bn和阻尼系数zeta是物理量纲系数由公式算出来改带宽时只改一个数。三路码的间隔用固定的 0.5 码片B1C 因为相关峰窄可以改到 0.3。间隔越小鉴别器增益越低但线性区越窄跟踪精度反而可能更好。实际调的时候先把间隔设 0.5 跑通再往下压每压一次看CN0有没有掉掉了就说明压过头了。4.2 鉴别器选型和环路带宽的经验值DLL 用的是归一化超前减滞后功率这个鉴别器在 ±0.5 码片内近似线性超过之后曲线会翻转所以跟踪一旦失锁鉴别器输出会乱跳环路自己拉不回来。判断失锁的简单办法是看CN0连续几个周期掉到 25 dB-Hz 以下就可以认为这个通道废了标记一下不再输出观测量别让它污染后面的定位解算。PLL 用四象限atan2线性区 ±90 度比二象限的atan稳。B1C 导频分量没有电文翻转可以直接用纯 PLL数据分量因为符号跳变PLL 鉴别器会在符号边界上产生 ±180 度的跳变这时候要么先用 FLL 把频率拉稳再切 PLL要么用判决反馈把符号剥掉。代码里最简单的做法是前 20 个积分周期只用 FLL频率误差收敛到几十赫兹以内再切 PLL切的时候相位从 0 开始不要沿用 FLL 阶段的相位。环路带宽的经验值静态场景 PLL 取 10 到 15 HzDLL 取 0.5 到 1 Hz车载场景 PLL 放到 25 到 30 HzDLL 放到 2 Hz高动态场景还要再宽但宽到一定程度噪声就压不住了。带宽改窄会让跟踪更平滑但牵引范围变小捕获给出的多普勒初值如果偏了 200 Hz 以上窄带宽拉不回来表现为频率一直朝一个方向缓慢漂。积分周期也不是越长越好。B1C 主码周期 10 ms数据分量的电文符号率是 1000 sps一个主码周期里有 10 个符号跨符号积分会互相抵消。导频分量没有这个问题可以老老实实积 10 ms。代码里如果两个分量都跟踪数据通道的积分周期要压到 1 ms 量级导频通道保持 10 ms最后在观测量层面把两个通道的结果合起来。4.3 比特同步、帧同步与电文解析跟踪出来的I_P是带符号的累加值符号里承载着电文比特但比特边界在哪里还不知道。比特同步就是在若干个可能的边界位置里找出正确的那个。直方图法最省事把一个比特周期内的I_P值按位置分箱累加能量最大的那个箱就是比特边界。% 比特同步直方图法 % 假设一个电文符号跨 Nsym 个主码周期每个周期输出一个 I_P Nsym cfg.symPerCode; % B1C 数据分量约为 1 ms 符号 / 10 ms 主码 acc zeros(1, Nsym); cnt zeros(1, Nsym); for k 1:numel(trk.I_P) bin mod(k-1, Nsym) 1; acc(bin) acc(bin) abs(trk.I_P(k)); cnt(bin) cnt(bin) 1; end acc acc ./ max(cnt, 1); [~, symBoundary] max(acc); % 比特边界所在位置找到边界之后把每个符号周期内的I_P按边界对齐相加得到符号判决值取符号位就是比特。这里有一个容易忽略的点极性可能整体反了。I_P的正负取决于本地载波和信号载波的相位关系如果本地载波取了共轭整串比特会全部翻转。解决办法是在帧同步阶段用已知的帧头或前导码做相关相关峰为正说明极性对为负就把整串比特取反。帧同步之后是电文解析。B-CNAV1 和 B-CNAV2 的帧结构、纠错编码都不一样解析函数分开写。B2a 的 B-CNAV2 用了 LDPC 编码解析前要先做译码译码失败就丢弃这一帧不要硬着头皮往下解否则星历参数里会混进垃圾值最后定位结果飞出几百公里都算正常。星历解析出来之后至少要做一次合理性检查轨道半长轴应该在 26000 km 量级倾角在 55 度左右偏心率是个小数任何一项对不上就说明帧对齐或者译码出了问题。电文解析的产出是每颗星的星历、钟差参数、健康标志。健康标志一定要用不健康的卫星即使观测量质量很好也要剔掉这是最省事的一道保险。5. 定位解算与常见问题排查5 个在实际调试中遇到的坑5.1 从伪距到 PVT 的最小二乘伪距的生成依赖两个时刻。接收时刻从采样点索引和采样率算出来发射时刻从电文里的参考时间加上码相位折算两者相减乘光速就是伪距。代码里用一个统一的t_rx数组管理接收时刻所有通道共用同一套时间基准这是保证通道间一致性的前提。function [pos, dtr, nIter] ls_pvt(satPos, pr, satClk, cfg) % satPos : N x 3ECEF 卫星位置单位米 % pr : N x 1伪距单位米已扣除卫星钟差 % satClk : N x 1卫星钟差等效距离单位米 % cfg : 初始位置与迭代上限 c 299792458; x [cfg.initPos(:); 0]; % [x y z c*dt] nIter 0; for it 1:cfg.maxIter nIter it; rho sqrt(sum((satPos - x(1:3).).^2, 2)); pred rho x(4) - satClk; H [-(satPos - x(1:3).)./rho, ones(size(rho))]; dx (H. * H) \ (H. * (pr - pred)); x x dx; if norm(dx(1:3)) 1e-3 break; end end pos x(1:3); dtr x(4) / c; end这段代码本身不复杂难的是前置条件。卫星位置必须和伪距在同一时刻卫星在信号发射时刻的位置要用迭代的方式逼近直接用接收时刻的位置会引入几十米的误差。另外至少要 4 颗卫星才能解出三维位置加钟差只有 3 颗的时候要么固定高度要么放弃。迭代收敛的判据用位置增量小于 1 mm 就够了但初始位置不能给得太离谱。初始位置偏差几百公里本文还有配套的精品资源点击获取