简介一套完整的MATLAB实现覆盖北斗B1C信号从捕获、跟踪到定位的闭环处理流程面向卫星导航、通信与信号处理方向的研究生和工程师。代码按信号仿真、捕获FFT峰值检测与相关器、跟踪PLL/DLL环路、帧同步与解码、最小二乘定位等阶段组织适合从算法原理到工程实现进行系统性学习。包内共51个文件包含42个m脚本和9个mat数据文件脚本负责各模块算法与主流程数据文件提供仿真中频/基带信号和中间结果资源包约25.56MB目录按流程阶段划分便于快速定位关键函数。目前已有530人学习下载。通过逐段运行和修改代码可直观理解B1C的BOC调制、多普勒频移补偿、伪距测量及PVT解算还能利用附带的电离层/对流层修正函数完善定位精度并附有中间结果可对照验证是一份既能入门又能深入的实践资料。1. 从北斗B1C信号的捕获到定位Matlab全链路该怎么搭拿到北斗B1C信号处理任务时最容易被低估的一点是它和GPS L1 C/A的差别远不止换了个频点B1C在1.023Mcps码率下使用10ms码周期和QMBOC调制伪码长度为10230捕获、跟踪、定位三个阶段各自都有需要单独处理的工程细节。在Matlab里把这条链路完整跑通不只是把算法按顺序摆出来更重要的是把副载波多峰、环路滤波器带宽、码周期到伪距的换算这些容易被忽略的参数变成可以量化、可以回查的对象。这篇内容适合做接收机算法预研、信号质量评估和教学仿真的人你手头有一批B1C中频数据或者你想在Matlab里从信号样本一直解算到接收机位置。2. 备料B1C信号样本的生成、读取与Matlab数据组织2.1 B1C信号结构和捕获跟踪前必须确定的三个参数B1C的公开参数里和信号处理链路直接相关的是这几个载波频率1575.42MHz伪码速率1.023Mcps码长10230个码片码周期10ms。调制方式上数据分量是BOC(1,1)导频分量是QMBOC(6,1,4/33)也就是说信号里同时存在1.023MHz和6.138MHz两路子载波分量。工程上处理B1C时捕获和跟踪通常都选导频分量因为导频支路没有电文调制相干积分不担心数据位跳变PLL鉴相器可以直接用atan2。数据分量留给后续B-CNAV1电文解析。开始写Matlab代码前有三个参数必须定下来否则后面捕获和跟踪的每一个单位换算都会出问题。第一是采样率。只处理BOC(1,1)主瓣时2.046MHz带宽就够但B1C里还有BOC(6,1)分量完整保留它需要更高的采样率。我一般用20.46MHz也就是码率的20倍BOC(1,1)主瓣完整保留BOC(6,1)的主瓣中心6.138MHz在带内但高频段被截掉一些对捕获和伪距测量影响很小。如果要做精细的信号质量分析采样率提到30MHz以上更稳妥。第二是数据格式。Matlab里处理B1C中频数据最常见的来源是软件接收机前端的bin文件两路int16交织存储先I后Q。读取时按交织顺序拆成复数注意文件里到底是不是这个顺序读反了信号会变成镜像频谱在频偏搜索时表现为多普勒频率符号异常。第三是码表。B1C主码长10230需要一个完整的本地码序列。可以从公开ICD文档的描述生成也可以从现成接收机工程里导出二进制码表Matlab里存成double数组即可值域取±1。注意B1C一个码周期是10ms不是GPS C/A的1ms码表长度对应的时间跨度直接决定后续伪距换算因子。2.2 读取数字中频IQ样本的Matlab入口函数下面是读取交织int16格式IQ数据的函数兼容复数单精度格式。function iq readB1cIq(filePath, nSamples, dtype) % 读取B1C前端的IQ数据文件 % filePath : bin文件路径 % nSamples : 需要读取的样本点数复数 % dtype : i16 表示交织int16cf 表示单精度复数I/Q fid fopen(filePath, rb); if strcmp(dtype, i16) raw fread(fid, 2*nSamples, *int16); iq complex(raw(1:2:end), raw(2:2:end)); elseif strcmp(dtype, cf) raw fread(fid, 2*nSamples, *single); iq complex(raw(1:2:end), raw(2:2:end)); else error(不支持的dtype: %s, dtype); end fclose(fid); end这个函数做了两件事按交织顺序拆出I路和Q路然后组成Matlab复数数组。后面捕获和跟踪都在这个复数序列上进行。读取逻辑之外有个容易忽略的点文件读取多少数据由捕获策略决定。捕获至少需要一个完整码周期也就是10ms在20.46MHz采样率下对应204600个样本。实际建议一次读够100ms以上这样捕获做完后跟踪模块可以直接用同一批数据不必重复读文件。如果要模拟弱信号场景文件读取后乘一个衰减因子再叠加高斯白噪声比在硬件上改射频增益更容易控制信噪比。2.3 采样率、码周期和FFT长度的换算关系捕获代码里处处是换算最常见的错误是把码周期当成1ms去算FFT长度和伪距。B1C的码周期是10ms那么单个码片时间宽度 1 / 1.023e6 秒 单个码周期时间宽度 10230 / 1.023e6 0.01 秒 20.46MHz采样率下一个码片对应20个样本 一个码周期对应 20.46e6 * 0.01 204600 个样本这段换算同时决定了捕获FFT的大小如果对10ms数据直接做FFT相关FFT长度至少取204600。这样一次FFT会比较大但Matlab还能接受。如果嫌慢可以先把信号降到2.046MHz再捕获代价是BOC(6,1)分量被滤掉捕获灵敏度和峰形判断会受影响这个取舍放到下一章说。码多普勒也是一个在跟踪前就要算好的参数。多普勒频率对码速率的影响按载波频率的比例折算载波偏移3000Hz时码速率偏移大约是1.023e6乘以3000除以1575.42e6约等于1.95Hz。这个偏移单独看很小但10ms的码周期里会积累约千分之二个码片跟踪跑几秒钟后相关峰会明显塌陷所以跟踪环路里必须把载波多普勒按比例补偿到码NCO上。参数取值说明伪码速率1.023 McpsB1C主码速率码长10230 chips一个主码周期码周期10 ms非GPS L1的1ms采样率20.46 MHz每码片20个样本10ms样本数204600一个码周期对应样本3. 捕获B1C副载波多峰下的码相位搜索与验证3.1 为什么B1C捕获不能直接套用BPSK接收机的峰值判决BPSK信号的自相关函数在主峰两侧单调下降捕获时找到相关幅值最大点就可以把码相位和频率交给跟踪环路。BOC(1,1)不是这样它的自相关函数在主峰两侧约±0.5码片处有幅度接近主峰一半的副峰非相干累加后副峰幅度约为主峰的1/4。这个特性导致的直接后果是如果捕获代码只做了一次“找全局最大”就返回结果跟踪环路很可能落在错锁的副峰上。码环会锁到一个偏了0.5码片的延迟上伪距直接产生约146米的偏差而且环路带宽足够窄时还很难自己挣脱出来。所以在B1C捕获的判决逻辑里除了找峰值还要检查峰值形态。我常用的办法是在检测到的峰值位置附近搜索是否存在离主峰约0.5码片、幅度在主峰60%以上的另一个峰。如果存在且多普勒频率bin相同说明这是BOC副载波的多峰现象而不是真正的信号到达时刻。更稳妥的做法是直接用ASPeCT类算法在相关运算里对BOC(1,1)的副峰做抑制但先把多峰特征暴露出来用于调试对教学和算法验证更有价值。3.2 基于FFT并行码相位搜索的捕获代码下面给出一个可直接运行的B1C导频捕获函数使用FFT并行搜索码相位多普勒频率网格外层循环。function [codePhase, fDop, peakMetric] b1cAcquisition(iq, code, fs, fIf, freqGrid, nNoncoh) % iq : 复中频或基带IQ长度至少 nNoncoh*10ms % code : B1C导频主码长度10230±1 % fs : 采样率 % fIf : 中频频率基带时填0 % freqGrid : 多普勒搜索网格 % nNoncoh : 非相干累加次数每个频点累加多少个10ms块 segLen round(fs * 0.01); % 单个码周期样本数 Nfft 2^nextpow2(segLen); % 补零到2的幂 n (0:segLen-1).; bestMetric 0; codePhase 0; fDop freqGrid(1); for k 1:numel(freqGrid) acc zeros(Nfft, 1); for m 1:nNoncoh segStart (m-1)*segLen 1; seg iq(segStart : segStartsegLen-1); lo exp(1j*2*pi*(fIf freqGrid(k))*n/fs); baseband seg .* lo; % 载波剥离 S fft(baseband, Nfft); C conj(fft(code, Nfft)); corr ifft(S .* C); acc acc abs(corr).^2; % 非相干累加 end [peakVal, idx] max(acc); if peakVal bestMetric bestMetric peakVal; fDop freqGrid(k); codePhase (idx - 1) / Nfft * 10230; % FFT索引回落到码片 end end peakMetric bestMetric / (mean(acc) eps); % 峰均比用于门限判断 end代码逻辑分成三段多普勒剥离、FFT相关、峰值提取。外层循环的每个频点都做一次完整的码相位搜索内层循环把每个10ms子块的相干相关结果做幅值平方累加提高弱信号下的检测概率。参数上面有几个点需要说明。freqGrid的步长要根据相干积分时长取10ms积分对应的频率分辨率约100Hz所以多普勒网格步长取80到100Hz比较合适步长再大边缘频率的积分损耗会明显增加。nNoncoh取1对强信号就够取4左右能多几dB处理增益但要注意非相干累加本身也有平方损耗并不是线性叠加增益。peakMetric用峰均比做判决本地噪声统计稳定时峰均比超过8到10通常可以判为捕获成功实际门限应在仿真数据上先做统计标定不要直接抄一个固定值。codePhase从FFT索引等比例回落到码片序号这个换算只给出码相位的粗估分辨率受限于采样间隔20.46MHz下约0.05码片更细的小数部分留给跟踪环路去逼近。3.3 捕获参数、副峰校验和弱信号时的非相干累加捕获相关参数的推荐取值按我的常用配置整理如下。参数推荐值说明多普勒搜索范围±5 kHz覆盖MEO卫星动态和本地晶振漂移多普勒搜索步长100 Hz与10ms相干积分匹配相干积分长度10 ms一个完整B1C主码周期非相干累加次数1~4弱信号时加大注意平方损耗检测门限峰均比8~10需按噪声统计实际标定码相位输出精度1个采样间隔约0.05码片捕获输出在交到跟踪之前建议做一次副峰校验。校验逻辑很简单在捕获得到的码相位附近以0.5码片为间隔搜索是否存在另一个幅度超过主峰60%的相关峰。如果有把本次检测标记为可疑可以先用BPSK-like方式重新捕获也就是把本地副载波先乘进输入信号再做单峰搜索损失约0.5dB的BOC(6,1)分量但峰形干净。弱信号场景下有一种比单纯增加nNoncoh更省时间的做法先用2.046MHz采样率的降采样数据粗捕获码相位定位到1码片以内然后在原始采样率上用窄多普勒网格精捕获。粗捕获阶段BOC(6,1)分量被滤掉但BOC(1,1)主瓣还在粗捕获的多普勒搜索网格可以放宽到250Hz整体计算量能降一个数量级。这个两级捕获策略在Matlab里跑真实长时间数据时非常实用因为直接把每一个频点的10ms数据做204600点FFT十几个频点几十次FFT下来单次捕获几十秒就过去了。4. 跟踪导频分量上的DLL/PLL环路与观测量生成4.1 跟踪环路结构与B1C导频分量的用途捕获给出的是粗略的码相位和多普勒频率跟踪环路负责把这两个值细化并持续锁定同时输出伪距和伪距率观测量。B1C跟踪的标准结构是码环加载波环码环用非相干超前减滞后功率鉴相器载波环用锁相环在导频分量上做就可以省略Costas环的180度相位模糊处理直接对同相支路做atan2鉴相。码环的典型相关器间隔取0.5码片超前和滞后支路分别位于码NCO相位两侧。载波环跟踪载波相位和频率码环跟踪码相位两个环路的误差信号都来自同一个即时支路的相关结果但鉴相方式不同环路带宽也不同。码环带宽窄载波环带宽宽这是接收机里最常见的不对称配置。B1C的10ms码周期对跟踪有一个实际好处积分时间可以覆盖整个主码周期而没有导航电文跳变相干积分的信噪比基础比1ms码周期的信号好。但这也意味着环路更新率不是1kHz而是100Hz对高动态场景下的载波环参数要求更高。如果载体动态大载波环带宽就要适当加宽或者改用二阶锁频环辅助三阶锁相环不能在固定带宽参数上一路走到黑。4.2 环路滤波器带宽、阻尼和系数的取值方法环路滤波器带宽和阻尼是跟踪调试里最值得花时间标定的参数。下面给出一个计算二阶环路滤波器系数的函数可以同时用于PLL和DLL的环路系数初始化。function [alpha, beta] calcLoopCoeff(Bn, zeta, T, loopType) % Bn : 环路噪声带宽单位Hz % zeta : 阻尼比通常取0.707 % T : 环路更新周期单位s % loopType : pll 或 dll if strcmp(loopType, pll) || strcmp(loopType, dll) wn Bn * 8 * zeta / (4*zeta^2 1); alpha 2 * zeta * wn * T; beta (wn * T)^2; else error(未知环路类型); end end这个函数的依据是典型二阶环路的噪声带宽与自然频率的关系式阻尼0.707时环路在噪声抑制和动态响应之间比较均衡。alpha对应比例项beta对应积分项。跟踪主循环里码环和载波环各自维护一组alpha、beta每次鉴相器输出的误差乘上这两个系数后累加到NCO频率上。我常用的带宽起点是载波环15到20Hz码环1Hz。这个组合对静态和低动态场景稳定。载波环带宽太窄时捕获残差里的多普勒估计误差会超过环路牵入范围太宽时热噪声引起的相位抖动变大。码环带宽则直接决定伪距噪声1Hz的码环带宽下C/N0为45dBHz时伪距噪声通常在几十厘米量级带宽收到0.5Hz能进一步压噪声但对动态应力的容忍度下降。4.3 跟踪主循环Matlab实现相关、鉴相与环路更新跟踪主循环按码周期逐个处理每个周期生成本地载波和本地码做三个支路的相关然后更新NCO。下面是一个可运行的教学版跟踪循环重点放在环路骨架和观测量输出上。function [codePhaseMesh, dopplerMesh, cn0Est] trackB1c(iq, code, fs, fDop0, codePhase0, params) % 跟踪单个B1C信号 % fDop0 : 捕获得到的多普勒频率 % codePhase0 : 捕获得到的码相位码片 % params : 结构体包含环路带宽等参数 segLen round(fs * 0.01); % 10ms样本数 epochLen 10230; % 每个码周期码片数 chipPerSec 1.023e6; carrierFreq fDop0; codeFreq chipPerSec * (1 fDop0 / 1575.42e6); carrierPhase 0; codeNco codePhase0; % 码NCO相位单位码片 T 0.01; [alphaPll, betaPll] calcLoopCoeff(params.pllBn, 0.707, T, pll); [alphaDll, betaDll] calcLoopCoeff(params.dllBn, 0.707, T, dll); carrierFreqAcc 0; codeFracout zeros(1, numel(iq)/segLen); dopplerOut zeros(1, numel(iq)/segLen); for ep 1 : floor(numel(iq) / segLen) segStart (ep-1)*segLen 1; seg iq(segStart : segStartsegLen-1); n (0:segLen-1).; lo exp(1j*(2*pi*carrierFreq*n/fs carrierPhase)).; s seg .* lo; % E/P/L三路码相位索引 halfChip 0.25; % 相关器间距0.5码片时超前滞后各偏0.25 pIdx mod(round(codeNco / epochLen * segLen), segLen) 1; eIdx mod(round((codeNco - halfChip) / epochLen * segLen), segLen) 1; lIdx mod(round((codeNco halfChip) / epochLen * segLen), segLen) 1; codeE code([eIdx:end 1:eIdx-1]); codeP code([pIdx:end 1:pIdx-1]); codeL code([lIdx:end 1:lIdx-1]); E mean(s .* codeE); P mean(s .* codeP); L mean(s .* codeL); % 码环鉴相非相干超前减滞后功率 dllDisc 0.5 * (abs(E)^2 - abs(L)^2) / (abs(E)^2 abs(L)^2 eps); % 载波环鉴相导频分量直接用atan2 pllDisc atan2(imag(P), real(P)); % 环路滤波与NCO更新 carrierFreqAcc carrierFreqAcc betaPll * pllDisc; carrierFreq carrierFreq alphaPll * pllDisc carrierFreqAcc; carrierPhase mod(carrierPhase 2*pi*carrierFreq*T, 2*pi); codeFreq chipPerSec * (1 carrierFreq / 1575.42e6); codeNco codeNco codeFreq * T; codeFracout(ep) mod(codeNco, epochLen); dopplerOut(ep) carrierFreq; end end代码的核心顺序是载波剥离、三支路相关、鉴相、环路滤波、NCO推进。载波剥离把输入信号从残余中频搬到基带三支路相关分别输出E、P、L三个相关值。码环鉴相器输出的是码相位误差单位码片载波环鉴相器输出的是载波相位误差单位弧度。这个教学版本做了一些简化码环的E、P、L相关值直接用码片索引换算成样本索引没有做插值环路滤波器用的是一阶比例加积分形式。工程实现里码环通常还要在相关器输出上做精细插值载波环在动态场景下会升级到三阶。直接用这个版本上真实数据时如果发现跟踪稳定但有小的固定偏差先检查E/L支路索引的舍入方式再看码环鉴相器是不是被BOC副峰带偏了。4.4 从码NCO和载波NCO生成伪距、伪距率跟踪环路稳定后观测量从两个NCO里读出来。码NCO的整周期计数和当前相位表示信号到达时刻在一个码周期内的位置载波NCO的相位累加值用于载波相位观测量。伪距生成的第一个坑是码周期基准。B1C一个码周期是10ms对应的时间是0.01秒而不是0.001秒。如果照搬GPS C/A的1ms基准伪距会直接差出约2700公里定位结果完全发散。伪距的粗测值可以这样算码相位延迟秒 codeFrac / 1.023e6 发射时刻 接收时刻 - 整数码周期 * 0.01 - 码相位延迟 伪距 光速 * (接收时刻 - 发射时刻)这里有一个整数码周期模糊度。码NCO只能给出当前周期内的相位无法直接告诉你信号是第几个周期发出的。真实接收机通过B-CNAV1电文里的时间信息解掉这个模糊度Matlab做定位算法验证时通常可以用已知接收时刻的粗时先验或者直接注入真值来绕开。伪距率则直接来自载波NCO的频率跟踪环路里的carrierFreq就是载波多普勒转换为伪距率只需乘上光速除以载波频率也就是多普勒频移与径向速度的经典关系。做定位解算时伪距率可以作为接收机速度观测方程的输出不过在只验位置解的流程里这项可以留到后续扩展。5. 定位最小二乘伪距解算与全链路一致性验证5.1 B1C单频伪距观测方程里哪些误差项不能省伪距观测方程的标准形式是伪距等于几何距离、接收机钟差、卫星钟差、电离层延迟、对流层延迟和噪声误差之和。B1C单频数据处理时电离层延迟没有双频组合可以用必须依靠模型或外部改正。在Matlab仿真的理想场景里卫星钟差和星历误差可以用注入真值的方式消掉对流层和电离层则可以按源数据里是否存在延迟来决定要不要加模型。试验里常见的做法是先跑一组不含电离层和对流层的干净数据确认定位算法本身没有偏差再叠加上模型改正去看误差预算。如果直接用真实采集数据而不做任何误差处理最小二乘残差会大得离谱这时优先检查伪距本身是否连续而不是急着调定位算法。卫星位置的计算也要注意B1C对应的卫星类型。BDS的MEO卫星用开普勒参数算位置GEO和IGSO卫星的位置计算里包含5度轨道倾角旋转这个旋转项不能丢。用RINEX导航文件时B1C对应的星历参数组在较新的RINEX版本里已经有定义读取时注意版本兼容性。5.2 带加权的最小二乘定位Matlab函数下面是定位解算的核心函数输入是至少4颗卫星的ECEF坐标和对应伪距输出是接收机位置、钟差和GDOP。function [pos, clkBias, gdop, residual] solveB1cPos(satPos, pr, pos0, b0, sigmaPr) % satPos : n×3 ECEF卫星坐标单位米 % pr : n×1 伪距单位米 % pos0 : 接收机初始位置 % b0 : 初始钟差单位米 % sigmaPr: 伪距噪声标准差单位米加权最小二乘用 pos pos0; b b0; n length(pr); W diag(1 ./ (sigmaPr.^2)); for it 1:8 r vecnorm(satPos - pos, 2, 2); H [(pos - satPos) ./ r, ones(n,1)]; res pr - (r b); dx (H * W * H) \ (H * W * res); pos pos dx(1:3); b b dx(4); if norm(dx(1:3)) 1e-3 break; end end residual pr - (vecnorm(satPos - pos, 2, 2) b); gdop sqrt(trace(inv(H * H))); end几何矩阵H的行向量是接收机到卫星的单位方向向量与1的拼接对应位置和钟差两个未知数列。残差更新时方程的线性化点在当前估计位置上所以每次迭代都要重算H和r不能只在循环外算一次。加权矩阵W来自伪距噪声方差的先验估计直接用单位阵就是普通最小二乘用信噪比或仰角换算的方差做权定位结果在高仰角卫星上会明显更稳。GDOP可以当做一个快速判据如果解算出来的GDOP大于10说明卫星几何结构很差定位结果即使残差很小也不可信。调试定位问题时先看GDOP再看残差顺序不要反。5.3 用已知时延注入检验捕获、跟踪、定位三段链路一致性全链路跑通后最值得做的一项验证是注入测试在Matlab里构造一个假的B1C信号注入已知的码时延和载波多普勒然后看捕获输出、跟踪输出和定位结果能不能还原出真值。构造方法是在本地生成B1C导频码后按预设的码相位延迟做循环移位乘上已知多普勒的载波再叠加带宽匹配的高斯白噪声。接收机位置和卫星位置都设为已知值伪距由几何关系直接算出。跑完捕获、跟踪、定位后把每一级的输出和注入值做对比。检查项注入值示例处理链输出允许误差捕获码相位5000.23码片5000.2左右半个采样间隔捕获多普勒-1300Hz-1300±50Hz搜索步长一半跟踪平均码相位5000.23码片5000.23±0.03码环噪声定位坐标已知ECEF差10m噪声与几何共同决定如果定位结果偏差超出预期按链路顺序排查先看捕获码相位是否离注入值在半个采样间隔以内再看跟踪输出相对捕获值有没有系统性偏移最后看伪距观测值是否和几何距离吻合。三者逐步收紧问题一定出在漏掉的那一级。更进一步的验证里可以对跟踪输出的载波相位残差做频谱分析对锁定后的即时支路相位做几十秒的FFT如果出现明显单频分量说明载波环存在稳态偏差或者有多径这时候回查捕获峰值形态和多峰校验记录。位置域的输出如果已经正确再考虑用卡尔曼滤波做位置平滑和视觉目标跟踪里的特征点跟踪用的是一套状态估计思路滤波本身掩盖不了伪距残差的问题所以先把最小二乘残差压干净再加滤波器。本文还有配套的精品资源点击获取