CN105425258A - 一种惯导系统辅助的高动态微弱信号gps捕获方法 - Google Patents

Abstract

一种惯导系统辅助的高动态微弱信号GPS捕获方法，它有五大步骤：一、对比相干、非相干以及差分相干累积捕获性能，选取累积方式进行组合，加强信号能量累积；二、对输入中频信号进行相干累积，并利用惯导系统和卫星星历对接收信号的多普勒参数及码相位进行估算，以压缩频率及相位搜索空间；三、对多普勒参数及码相位估算值进行卡尔曼滤波平滑，并利用平滑结果计算本地信号与接收信号的动态频率/相位偏移量；四、采用循环平移算法补偿相干累积输出矩阵的动态偏移量；五、对补偿后的输出矩阵进行非相干累积，将输出结果与门限比较，完成捕获判决。本发明改善了GPS信号捕获的动态容忍性能，实现了高动态环境中GPS信号的高灵敏度捕获。

Description

一种惯导系统辅助的高动态微弱信号GPS捕获方法

一、技术领域

本发明涉及一种惯导系统辅助的高动态微弱信号全球定位系统(GPS,GlobalPositioningSystem)捕获方法，属于导航技术领域。

二、背景技术

GPS接收机对接收信号的处理包括捕获、跟踪、导航电文提取及导航解算等过程，信号捕获是首要也是关键的一步。随着GPS应用环境的复杂化，对高动态、高灵敏度GPS信号捕获技术的需求不断提升，其典型应用为地球同步轨道卫星、高轨卫星的自主定轨。高灵敏度捕获技术一般通过延长信号累积时间以获得高载噪比增益，但是在高动态环境中，载体速度及加速度变化范围较大，加剧了多普勒频移及多普勒频移率的不确定性，导致本地信号与接收信号的频率/相位匹配误差随着累积时间的延长而不断增大，降低了信号的检测性能。因此，高动态、高灵敏度信号捕获技术必须在提高载噪比增益的同时进行多普勒补偿。

为了在低载噪比环境下保持GPS信号捕获的灵敏度，需要通过相干累积获得增益以增加后处理载噪比，减弱噪声损耗。相干累积算法(COH,CoherentIntegration)利用了累积周期内信号的相关性，通过直接累积本地信号与接收信号的相关输出，减小信号功率损耗，提高捕获增益，但是由于处理过程中保留了相位信息，累积时长将受到比特翻转的限制。因此，为了延长累积时间以提高信号处理增益，通常采用能量累积算法进一步累积信号，例如非相干累积(NCH,Non-CoherentIntegration)和差分相干累积(DFC,DifferentialCoherentIntegration)等算法。非相干算法将信号相关结果平方后再累积，通过包络检波减小了由导航电文和残留载波多普勒频移带来的相位偏移，但由于噪声也被平方，导致处理增益存在平方损耗；差分相干算法是将相邻时刻的相关结果共轭相乘并累积，由于噪声在累积时间内相关性较弱，噪声能量不随累积时间累积，故差分相干抗干扰能力较强；但当累积时间较长时，将导致累积结果正负抵消衰减为零，因此，差分相干仅适用于累积时间较短的情况。

而针对高动态环境下的信号捕获问题，通常采取的解决方案包括基于快速傅里叶变换(FFT,FastFourierTransform)的频移捕获算法、基于序列估计的捕获算法、部分匹配滤波算法以及惯导系统辅助的捕获算法等。基于FFT的捕获法可搜索大范围多普勒频移，但由于FFT运算量庞大，捕获效率低下；基于序列估计的捕获法为了准确估计码序列需要足够长的观测时间，不利于信号的快速捕获；部分匹配滤波法将信号分段相关后进行FFT运算，虽然降低了频域运算量，但将导致相关器输出载噪比的衰落；惯导辅助的捕获法利用惯导输出的位置、速度信息与卫星星历，估算接收信号的多普勒频移，预先求取捕获的频率搜索范围，缩短了信号的捕获时间，但是在高动态、低载噪比环境下，载体与卫星间相对加速度引起的多普勒频移率较大，本地信号与接收信号偏移量不断扩大，使得惯导辅助的捕获法难以对信号进行长时间相干累积，导致捕获灵敏度大幅衰减。

为了满足高动态微弱信号的捕获要求，本发明提出了一种惯导系统辅助的高动态微弱信号GPS捕获方法。

三、发明内容

针对现有技术中存在的问题，本发明提出了一种惯导系统辅助的高动态微弱信号GPS捕获方法。它采用了相干/非相干组合累积方式以提高微弱信号在高动态环境中的载噪比处理增益；同时通过惯导系统和卫星星历信息估算信号的多普勒参数和码相位，以缩小搜索空间，提高捕获速度；利用滤波算法平滑多普勒参数和码相位估算值，并通过循环频移算法校正本地信号与接收信号之间的频率/相位动态偏移量，以增强信号捕获的动态性能和快速反应能力，从而实现高动态环境中GPS信号的高灵敏度捕获。

本发明一种惯导系统辅助的高动态微弱信号GPS捕获方法，具体包括以下步骤：

步骤一：对比相干、非相干以及差分相干累积捕获性能，选取累积方式进行组合，加强信号能量累积；

a.相干累积捕获性能

L1波段接收信号经过射频前端下变频及模数转换后得到数字中频信号，以采样频率f_s＝1/T_s进行采样，得到nT_s时刻采样信号以及本地相关信号模型分别为

$s_n^r={\mathrm{AC}}_n\left(\right(1+\mathrm{\η}\left)\right({\mathrm{nT}}_s-\mathrm{\τ}\left)\right)D_n\left({\mathrm{nT}}_s\right)\exp \left(2\mathrm{\π}\right(f_{IF}+f_d){\mathrm{nT}}_s+{\mathrm{\φ}}_0)+W_n---\left(1\right)$

$s_n^l=C_n\left(\right(1+\mathrm{\η}\left)\right({\mathrm{nT}}_s-\widehat{\mathrm{\τ}}\left)\right)\exp \left(2\mathrm{\π}\right(f_{IF}+{\hat{f}}_d){\mathrm{nT}}_s+{\widehat{\mathrm{\φ}}}_0)---\left(2\right)$

式中，下标n表示当前采样点；上标r表示采样信号，上标l表示本地信号；C(·)为粗捕获(C/A,Coarse/Acquisition)码序列；D(·)为导航数据；f_IF为载波中频频率；τ和分别为采样信号的真实码相位传播延迟以及本地信号的码相位估计值；f_d和分别为采样信号的真实多普勒频移和本地信号多普勒频移估计值；φ₀和分别为采样信号初始载波相位和本地信号初始载波相位估计值；多普勒频移f_d同时使C/A码速率产生偏移，偏移因子W_n为高斯白噪声。

假设由载体与GPS卫星相对加速度引起的多普勒频移率为则真实多普勒频移与多普勒频移估算值之间的关系为

$f_d-{\hat{f}}_d={\mathrm{\Δf}}_0+{\stackrel{\·}{f}}_d\·{\mathrm{nT}}_s/2---\left(3\right)$

式中，Δf₀为初始多普勒频移估计偏移量。

将采样信号与本地信号进行相关运算，则有

$\begin{array}{c}{s}_n^r\·s_n^l={\mathrm{AR}}_n(\mathrm{\τ}-\widehat{\mathrm{\τ}})\exp j\left(2\mathrm{\π}\right(f_d-{\hat{f}}_d){\mathrm{nT}}_s+{\mathrm{\φ}}_n-{\widehat{\mathrm{\φ}}}_n)\\ ={\mathrm{AR}}_n\left(\mathrm{\Δ}\mathrm{\τ}\right)\exp j\left(2\mathrm{\π}\right({\mathrm{\Δf}}_0+{\stackrel{\·}{f}}_d\·{\mathrm{nT}}_s/2){\mathrm{nT}}_s+{\mathrm{\Δ\φ}}_n)\\ ={\mathrm{AR}}_n\left(\mathrm{\Δ}\mathrm{\τ}\right)\exp j\left({\mathrm{\Δ\φ}}_n\right)\exp j\left(2{\mathrm{\π\Δf}}_0{\mathrm{nT}}_s\right)\exp j\[\mathrm{\π}{\stackrel{\·}{f}}_d{\left({\mathrm{nT}}_s\right)}^2\]\end{array}---\left(4\right)$

进一步，假设相干累积时间为T_coh，对式(4)中的相关结果号进行相干累积，则其输出检测量为

$\begin{array}{c}{Z}_{COH}=\underset{n=0}{\overset{N-1}{\Σ}}{\mathrm{AR}}_n\left(\mathrm{\Δ}\mathrm{\τ}\right)\exp j\left({\mathrm{\Δ\φ}}_n\right)\exp j\left(2{\mathrm{\π\Δf}}_0{\mathrm{nT}}_s\right)\exp j\[\mathrm{\π}{\stackrel{\·}{f}}_d{\left({\mathrm{nT}}_s\right)}^2\]\\ ={\mathrm{AR}}_n\left(\mathrm{\Δ}\mathrm{\τ}\right)N\sin c\left({\mathrm{\π\Δf}}_0{T}_{coh}\right)\exp \left({\mathrm{\π\Δf}}_0\right(N-1)T_s+{\mathrm{\Δ\φ}}_n)\underset{n=2}{\overset{N-1}{\Σ}}\exp j\[\mathrm{\π}{\stackrel{\·}{f}}_d{\left({\mathrm{nT}}_s\right)}^2\]\end{array}---\left(5\right)$

式中，N＝T_coh/T_s，表示相干累积的采样点数。

由式(5)可知，当不存在多普勒频移率，即时，相干累积输出量可简化为

$\begin{array}{c}{Z}_{COH}=\underset{n=0}{\overset{N-1}{\Σ}}{\mathrm{AR}}_n\left(\mathrm{\Δ}\mathrm{\τ}\right)\exp j\left({\mathrm{\Δ\φ}}_n\right)\exp j\left(2{\mathrm{\π\Δf}}_0{\mathrm{nT}}_s\right)\\ ={\mathrm{AR}}_n\left(\mathrm{\Δ}\mathrm{\τ}\right)N\sin c\left({\mathrm{\π\Δf}}_0{T}_{coh}\right)\exp \left({\mathrm{\π\Δf}}_0\right(N-1)T_s+{\mathrm{\Δ\φ}}_n)\end{array}---\left(6\right)$

此时，信号检测量的幅值仅由载波频率误差Δf₀和码相位误差Δτ两个变量决定，当Δf₀和Δτ均为零时，表示本地信号与接收信号的频率和相位均已对齐，此时输出检测量Z_COH即为捕获的相关峰。

而当多普勒频移率时，根据式(5)可知，的累积结果可表示为sinc函数，即Z_COH将沿着sinc函数主峰的下降梯度衰减。为了定量说明对检测量Z_COH幅值的损耗程度，将其定义为多普勒频移率损耗因子其中T为信号累积时间。对比式(5)和式(6)可知，相干累积算法的多普勒频移率损耗因子为

${\mathrm{\α}}_d({\stackrel{\·}{f}}_d,T_{coh})=\underset{n=0}{\overset{N-1}{\Σ}}\exp j\[\mathrm{\π}{\stackrel{\·}{f}}_d{\left({\mathrm{nT}}_s\right)}^2\]---\left(7\right)$

由于表现为sinc函数，因此多普勒频移率越大，幅度越小；同时随着累积时间T_coh的延长，幅度将进一步减小，导致信号累积能量的加速衰减。因此，为了在高动态环境中实现高灵敏度捕获，必须限制相干累积时间并进行多普勒频移率的补偿。

b.非相干、差分相干捕获性能

根据式(5)，可得相干累积后的同相、正交支路信号I_k和Q_k分别为

I_k＝AR_n(Δτ)sinc(πΔf₀T_coh)·cos(πΔf₀(N-1)T_s+Δφ_n)(8)

Q_k＝AR_n(Δτ)sinc(πΔf₀T_coh)·sin(πΔf₀(N-1)T_s+Δφ_n)(9)

而非相干和差分相干累积后得到的检测量可分别表示为

$Z_{NCH}=\underset{k=0}{\overset{T_{NCH}-1}{\Σ}}(I_k^2+Q_k^2)---\left(10\right)$

$Z_{DFC}=\underset{k=2}{\overset{T_{DFC}+1}{\Σ}}\left(\sqrt{{(I_k{I}_{k-1}+Q_k{Q}_{k-1})}^2+{(Q_k{I}_{k-1}-I_k{Q}_{k-1})}^2}\right)---\left(2\right)$

式中，Z_NCH，Z_DFC分别为非相干、差分相干输出矩阵；T_NCH和T_DFC分别为非相干和差分相干的累积周期。

将式(8)和式(9)分别代入式(10)和式(11)，得到非相干和差分相干的多普勒频移率衰减因子分别为

${\mathrm{\α}}_d^{NCH}({\stackrel{\·}{f}}_d,T_{NCH})=\underset{m=0}{\overset{T_{NCH}-1}{\Σ}}{(\underset{n=0}{\overset{N-1}{\Σ}}\exp j\[\mathrm{\π}{\stackrel{\·}{f}}_d{(m+{\mathrm{nT}}_s)}^2\])}^2---\left(12\right)$

${\mathrm{\α}}_d^{DFC}({\stackrel{\·}{f}}_d,T_{DFC})=\underset{m=1}{\overset{T_{DFC}-1}{\Σ}}(\underset{m=0}{\overset{N-1}{\Σ}}\exp j\[\mathrm{\π}{\stackrel{\·}{f}}_d{(m-1+{\mathrm{nT}}_s)}^2\]\·\underset{n=0}{\overset{N-1}{\Σ}}\exp j\[-\mathrm{\π}{\stackrel{\·}{f}}_d{(m+{\mathrm{nT}}_s)}^2\])---\left(13\right)$

在低动态环境中，非相干和差分相干的累积增益包括相干增益及非相干或差分相干损耗；但是在高动态环境中，由于多普勒频移率的存在，本地信号与接收信号的频率/相位偏移量随着累积时间的延长不断增大，导致相应累积增益的衰落。因此，高动态累积增益的计算模型需引入多普勒频移率损耗因子则非相干和差分相干累积增益模型分别为

$G_{NCH}\left(T_{NCH}\right)=G_{COH}\left(T_{NCH}\right)-L_{NCH}\left(T_{NCH}\right)-10\lg \[{\mathrm{\α}}_d^{NCH}({\stackrel{\·}{f}}_d,T_{NCH})\]---\left(14\right)$

$G_{DFC}\left(T_{DFC}\right)=G_{COH}\left(T_{DFC}\right)-L_{DFC}\left(T_{DFC}\right)-l0\lg \[{\mathrm{\α}}_d^{DFC}({\stackrel{\·}{f}}_d,T_{DFC})\]---\left(15\right)$

式中，G_COH(·)为相干累积增益；L_NCH(·)和L_DFC(·)分别为非相干和差分相干损耗。通过数值仿真可得，高动态条件下，随着信号累积时间的延长，非相干和差分相干的累积增益均不断下降，但是相比于非相干累积方式，差分相干增益的下降速度更快。因此，为了在高动态、低载噪比环境中获得较高的累积增益，应选择非相干方式对相干累积信号进行能量累积。

综上，组合累积方式选取相干和非相干累积方式进行组合。

步骤二：对输入中频信号进行相干累积(其结果见式(5)所示)，并利用惯导系统和卫星星历对接收信号的多普勒参数及码相位进行估算，以压缩频率及相位搜索空间；

a.多普勒参数估算

在信号捕获算法中引入惯导输出的载体位置、速度和加速度信息，结合卫星星历，则可得到由载体与卫星视线方向相对运动所产生的多普勒频移以及多普勒频移率，其计算式分别为

${\hat{f}}_d=-f_{L1}\frac{\left|\right|{\left({\hat{v}}_r\right)}_{LOS}\left|\right|}{c}=-f_{L1}\frac{({\hat{v}}_s-{\hat{v}}_u)\·{\hat{r}}_{LOS}}{c}---\left(16\right)$

$\widehat{\stackrel{\·}{f}}=-f_{L1}\frac{\left|\right|{\left({\hat{a}}_r\right)}_{LOS}\left|\right|}{c}=-f_{L1}\frac{({\hat{a}}_s-{\hat{a}}_u)\·{\hat{r}}_{LOS}}{c}---\left(17\right)$

式中，为GPS星历提供的卫星速度和加速度矢量；分别为惯导系统输出的载体速度和加速度矢量；为卫星与载体间的单位视线矢量，可由卫星与载体间的相对位置计算得到。

将由式(16)计算得到的多普勒频移估算值作为频率搜索中心，频率搜索范围根据多普勒频移的计算精度决定，其中多普勒频移误差由星历误差和惯导系统误差两部分组成。设由星历计算所得的多普勒频移计算误差为惯导输出的速度信息误差为则可以根据星历中给出的轨道误差项估算，而可通过下式估算

${\mathrm{\Δ}}_v^{SINS}=\underset{T_{COH}}{\∫}g\left(\underset{T_{COH}}{\∫}{C}_b^m\mathrm{\ϵ}dt\right)dt+\underset{T_{COH}}{\∫}{C}_b^m\▿dt---\left(18\right)$

式中，g为当地重力加速度；为本体系到导航系的姿态转换矩阵；ε和分别为陀螺仪和加速度计零偏。

因此，根据式(16)可得多普勒频移误差的方差为

${\mathrm{\σ}}_{Dopp}^2={\left(\frac{f_{L1}}{c}\right)}^2{\hat{r}}_{LOS}{\left({\mathrm{\Δ}}_v^{SINS}\right)}^2{\hat{r}}_{LOS}^T+{\left({\mathrm{\Δ}}_{Dopp}^{Eph}\right)}^2---\left(19\right)$

将式(19)计算得到的多普勒频移误差最大值作为搜索边界以压缩搜索范围，从而达到提高捕获效率的目的。

b.码相位估算

由于GPS卫星信号传输链路较长，相比于发射码，接收码将产生τ时长的传输延迟；另外，C/A码速率受到多普勒频移的影响而不断变化，导致接收码相位延时的不确定性。因此，随着信号累积时间的延长，相比于发射码相位，接收码相位将出现大幅超前或滞后。故在高动态、低载噪比环境下，需要对接收码相位进行预先估算，以生成较为准确的本地复制码，提高本地码与接收码的初始对准精度，达到缩小码相位搜索区间、提高捕获效率的目的。

为了估算接收码相位，建立t_k采样时刻GPS接收信号的C/A码序列模型为

$C(t_k,\mathrm{\τ},f_d,{\stackrel{\·}{f}}_d)=C\left(\right(1+\mathrm{\η}\left)\right(t_k-\mathrm{\τ}\left)\right)---\left(20\right)$

$\mathrm{\η}=\frac{f_d+{\stackrel{\·}{f}}_d\·t_k/2}{f_{L1}}---\left(21\right)$

式中，η表示由多普勒效应引起的C/A码速率收缩因子。

若载体到卫星的相对位置矢量为r_s,u＝r_s-r_u，则GPS信号的传播时间为

${\hat{t}}_{prop}=\left|\right|{\hat{r}}_{s,u}\left|\right|/c=\left|\right|{\hat{r}}_s-{\hat{r}}_u\left|\right|/c---\left(22\right)$

式中，为GPS星历提供的卫星位置矢量；为惯导系统解算得到的载体位置矢量。

因此，接收码相位的估算模型为

$\widehat{\mathrm{\τ}}=rem\left(\right(t_{GPS}+t_{corr}-{\hat{t}}_{prop})\·f_{CA},\frac{c}{f_{CA}(1+\mathrm{\η})})---\left(23\right)$

式中，t_GPS为当前C/A码的接收时刻(GPS时)；t_corr为卫星时钟修正量，可表示为

t_corr＝a_f0+a_f1(t_GPS-t_oc)+a_f2(t_GPS-t_oc)²(24)

式中，a_fi,i＝0,1,2为i阶卫星时钟修正系数，可从导航电文的第一子帧中获得；t_oc为第一子帧中第一数据块的参考时间。

码相位估算值误差包括载体-卫星相对位置误差及时间误差，可表示为

${\mathrm{\σ}}_{\mathrm{\τ}}^2=4{\mathrm{\σ}}_{pos}^2{\cos }^2\left({\mathrm{\φ}}_{el}\right)+{\mathrm{\σ}}_t^2---\left(25\right)$

式中，φ_el为卫星仰角；Δ_eph为由卫星星历计算得到的GPS卫星定轨误差，可通过导航电文中的轨道误差项估算得到，和分别为惯导系统位置输出在纬度和经度方向上的误差，可通过下式估算

$\begin{array}{c}{\mathrm{\σ}}_{pos}^2={\mathrm{\Δ}}_{eph}^2+{\mathrm{\Δ}}_{{\mathrm{SINS}}_L}^2+{\mathrm{\Δ}}_{{\mathrm{SINS}}_{\mathrm{\λ}}}^2\\ {\mathrm{\Δ}}_{{\mathrm{SINS}}_L}=\frac{1}{R}\underset{T_{COH}}{\∫}{\mathrm{\Δ}}_{v_y}^{SINS}\\ {\mathrm{\Δ}}_{{\mathrm{SINS}}_{\mathrm{\λ}}}=\frac{\sec L}{R}\underset{T_{COH}}{\∫}{\mathrm{\Δ}}_{v_x}^{SINS}\end{array}---\left(26\right)$

式中，和分别为导航坐标系x,y轴方向的惯导系统速度误差，可由式(18)计算获得。

根据式(23)得到码相位估算值，作为相位搜索空间中心；并根据式(25)计算相位误差的最大值，作为搜索边界，从而压缩搜索空间，提高高动态环境中的捕获效率。

步骤三：对多普勒参数及码相位估算值进行卡尔曼滤波平滑，并利用平滑结果计算本地信号与接收信号的动态频率/相位偏移量；

为了精确补偿不同相干累积周期间本地信号与接收信号间的动态偏移量，需要保证累积矩阵循环平移的精度，因此预先平滑处理信号捕获参数，选取卡尔曼滤波器的状态量和量测量分别为 $X={\left(\begin{array}{ccc}\mathrm{\τ}& f_d& {\stackrel{\·}{f}}_d\end{array}\right)}^T---\left(27\right)$

$Z={\left(\begin{array}{ccc}{(r_s-r_u)}_i& {(v_s-v_u)}_i& {(a_s-a_u)}_i\end{array}\right)}_{9\×1}^T,i=x,y,z---\left(28\right)$

则根据式(16)、式(17)和式(23)可得状态方程和量测方程为

$\left[\begin{array}{c}\stackrel{\·}{\mathrm{\τ}}\\ {\stackrel{\·}{f}}_d\\ {\stackrel{\·\·}{f}}_d\end{array}\right]=\left[\begin{array}{ccc}0& \frac{f_{C/A}}{f_{L1}}& 0\\ 0& 0& 1\\ 0& 0& 0\end{array}\right]\left[\begin{array}{c}\mathrm{\τ}\\ f_d\\ {\stackrel{\·}{f}}_d\end{array}\right]+\left[\begin{array}{ccc}2{\mathrm{\πf}}_{C/A}& 0& 0\\ 0& 2{\mathrm{\πf}}_{L1}& 0\\ 0& 0& 2\mathrm{\π}/{\mathrm{\λ}}_{L1}\end{array}\right]\left[\begin{array}{c}{w}_b\\ w_d\\ w_a\end{array}\right]---\left(29\right)$

Z＝HX+V(30)

$H={\left[\begin{array}{ccc}{\left(r_{LOS}\right)}_{3\×1}^T& 0_{3\×1}& 0_{3\×1}\\ 0_{3\×1}& {\left(r_{LOS}\right)}_{3\×1}^T& 0_{3\×1}\\ 0_{3\×1}& 0_{3\×3}& {\left(r_{LOS}\right)}_{3\×1}^T\end{array}\right]}_{9\×3}---\left(31\right)$

$V=diag\left(\begin{array}{ccc}{\mathrm{\σ}}_{\mathrm{\τ}}^2& {\mathrm{\σ}}_{Dopp}^2& {\mathrm{\σ}}_{Dopp\_Rate}^2\end{array}\right)---\left(32\right)$

式中，w_b为时钟偏差；w_d为时钟漂移；w_a为视线加速度的驱动噪声。

利用卡尔曼滤波器实时输出的信号参量平滑值，计算当前相干累积周期的累积矩阵输出与下一周期输出在载波频率搜索空间内的动态偏移量以及在码相位搜索空间内的动态偏移量其中

$\begin{array}{c}\mathrm{\δ}\hat{f}=(f_d-{\tilde{f}}_d)\·\mathrm{\Δ}t+{\stackrel{\·}{f}}_d\·{\mathrm{\Δt}}^2/2\\ \mathrm{\δ}\widehat{\mathrm{\τ}}=(\mathrm{\τ}-\widetilde{\mathrm{\τ}})\·\mathrm{\Δ}t+\frac{f_{CA}}{f_{L1}}\·{\stackrel{\·}{f}}_d\·{\mathrm{\Δt}}^2/2\end{array}---\left(33\right)$

式中，分别为当前搜索的多普勒频移值和码相位值；Δt为相干累积周期。

步骤四：采用循环平移算法补偿相干累积输出矩阵的动态偏移量；

对于第i个相干累积周期，由多普勒频移率导致的累积输出矩阵动态频率和相位偏移量分别为和在进行非相干累积之前，为了消除前i个周期的累积动态偏移量，第i个累积输出矩阵需沿频率搜索搜索空间循环平移个搜索单元，沿相位搜索搜索空间循环平移个搜索单元，若捕获搜索的频率分辨率和相位分辨率分别为f_Bin和τ_Bin，则和可表示为

$X_f^i=\[\frac{\underset{n=1}{\overset{i}{\Σ}}\mathrm{\δ}{\hat{f}}_i}{f_{Bin}}\],X_{\mathrm{\τ}}^i=\[\frac{\underset{n=1}{\overset{i}{\Σ}}\mathrm{\δ}{\widehat{\mathrm{\τ}}}_i}{{\mathrm{\τ}}_{Bin}}\]---\left(34\right)$

式中，[·]为取整运算。

步骤五：对补偿后的输出矩阵进行非相干累积(其算法见式(12)所示)，将输出结果与门限比较，完成捕获判决；

循环平移捕获周期内所有的相干输出矩阵后，对输出矩阵进行非相干累积，以进一步增强微弱信号的累积能量，提高载噪比处理增益，从而满足信号捕获最大峰值检测的灵敏度要求。

非相干累积输出的峰值为

$\begin{array}{c}{\mathrm{CP}}_{nch}=A^2{N}^2\underset{k=1}{\overset{n}{\&Sigma;}}{|R(f_{Bin}{X}_f^i-\underset{k=i}{\overset{n-1}{\&Sigma;}}\mathrm{\&delta;}{\widehat{\mathrm{\&tau;}}}_k)\exp \&lsqb;{\mathrm{\&pi;NT}}_s({\mathrm{\&tau;}}_{Bin}{X}_{\mathrm{\&tau;}}^i-\underset{k=i}{\overset{n-1}{\&Sigma;}}\mathrm{\&delta;}{\hat{f}}_k)\&rsqb;|}^2\\ =A^2{N}^2\underset{k=1}{\overset{n}{\&Sigma;}}{\left|\begin{array}{c}R({\mathrm{\&tau;}}_{Bin}\&CenterDot;{\&lsqb;\underset{k=i}{\overset{n-1}{\&Sigma;}}\mathrm{\&delta;}{\widehat{\mathrm{\&tau;}}}_k/{\mathrm{\&tau;}}_{Bin}\&rsqb;}_{round}-\underset{k=i}{\overset{n-1}{\&Sigma;}}\mathrm{\&delta;}{\widehat{\mathrm{\&tau;}}}_k)\&CenterDot;\\ \exp \&lsqb;{\mathrm{\&pi;NT}}_s(f_{Bin}\&CenterDot;{\&lsqb;\underset{k=i}{\overset{n-1}{\&Sigma;}}\mathrm{\&delta;}{\hat{f}}_k/f_{Bin}\&rsqb;}_{round}-\underset{k=i}{\overset{n-1}{\&Sigma;}}\mathrm{\&delta;}{\hat{f}}_k)\&rsqb;\end{array}\right|}^2\\ {\mathrm{\&omega;}}_{Bin}\&CenterDot;{\&lsqb;\underset{k=i}{\overset{n-1}{\&Sigma;}}\mathrm{\&delta;}{\widehat{\mathrm{\&omega;}}}_k/{\mathrm{\&omega;}}_{Bin}\&rsqb;}_{round}-\underset{k=i}{\overset{n-1}{\&Sigma;}}\mathrm{\&delta;}{\widehat{\mathrm{\&omega;}}}_k<<{\mathrm{\&omega;}}_{Bin}\mathrm{\&omega;}=\mathrm{\&tau;},f\end{array}---\left(35\right)$

若CP_nch＞2Threshold，则认为该颗卫星已成功捕获，其中，Threshold一般取非相干累积输出的第二峰值。

本发明的优点在于：

(1)本发明提出一种惯导系统辅助的高动态微弱信号GPS捕获方法，利用相干/非相干组合累积方式，一方面能够增强信号累积的处理增益，另一方面能够抑制载体动态应力的影响，以提高高动态条件下微弱信号的捕获性能。

(2)本发明提出一种惯导系统辅助的高动态微弱信号GPS捕获方法，利用惯导信息和卫星星历预先求取信号捕获参量，能够压缩频率/相位搜索空间，加快捕获速度。

(3)本发明提出一种惯导系统辅助的高动态微弱信号GPS捕获方法，利用循环平移算法修正不同周期相干累积输出的频率/相位动态偏移量，从而进一步提高信号累积的载噪比处理增益，改善了GPS信号捕获的动态容忍性能，实现了高动态环境中GPS信号的高灵敏度捕获。

四、附图说明

图1为本发明提出的一种惯导系统辅助的高动态微弱信号GPS捕获方法流程图；

图2为本发明提出的一种惯导系统辅助的高动态微弱信号GPS捕获方法结构示意图；

图3为本发明提出的循环频移算法原理图；

图2中：

N—相干累积周期n—非相干累积周期|·|²—取模的平方

图3中：

f—载波频率搜索空间τ—码相位搜索空间

—第k个相干累积输出矩阵载波频率动态偏差量

—第k个相干累积输出矩阵码相位动态偏差量

n—非相干累积周期

五、具体实施方式

下面将结合附图对本发明作进一步的详细说明。

本发明提出了一种惯导系统辅助的高动态微弱信号GPS捕获方法。采用了相干/非相干组合累积方式以提高微弱信号在高动态环境中的载噪比处理增益；同时通过惯导系统和卫星星历信息估算信号的多普勒参数和码相位，以缩小搜索空间，提高捕获速度；利用滤波算法平滑多普勒参数和码相位估算值并通过循环频移算法校正本地信号与接收信号之间的频率/相位动态偏移量以增强信号捕获的动态性能和快速反应能力，从而实现高动态环境中GPS信号的高灵敏度捕获。

见图1，本发明一种惯导系统辅助的高动态微弱信号GPS捕获方法，具体包括以下步骤

步骤一：对比相干、非相干以及差分相干累积捕获性能，选取累积方式进行组合，加强信号能量累积；

a.相干累积捕获性能

L1波段接收信号经过射频前端下变频及模数转换后得到数字中频信号，以采样频率f_s＝1/T_s(9.548MHz)进行采样，得到nT_s时刻采样信号以及本地相关信号模型分别为

$s_n^r={\mathrm{AC}}_n\left(\right(1+\mathrm{\&eta;}\left)\right({\mathrm{nT}}_s-\mathrm{\&tau;}\left)\right)D_n\left({\mathrm{nT}}_s\right)\exp \left(2\mathrm{\&pi;}\right(f_{IF}+f_d){\mathrm{nT}}_s+{\mathrm{\&phi;}}_0)+W_n---\left(36\right)$

$s_n^l=C_n\left(\right(1+\mathrm{\&eta;}\left)\right({\mathrm{nT}}_s-\widehat{\mathrm{\&tau;}}\left)\right)\exp \left(2\mathrm{\&pi;}\right(f_{IF}+{\hat{f}}_d){\mathrm{nT}}_s+{\widehat{\mathrm{\&phi;}}}_0)---\left(37\right)$

式中，下标n表示当前采样点；上标r表示采样信号，上标l表示本地信号；C(·)为C/A码序列；D(·)为导航数据；f_IF为载波中频频率；τ和分别为采样信号的真实码相位传播延迟以及本地信号的码相位估计值；f_d和分别为采样信号的真实多普勒频移和本地信号多普勒频移估计值；φ₀和分别为采样信号初始载波相位和本地信号初始载波相位估计值；多普勒频移f_d同时使C/A码速率产生偏移，偏移因子W_n为高斯白噪声。

假设由载体与GPS卫星相对加速度引起的多普勒频移率为则真实多普勒频移与多普勒频移估算值之间的关系为

$f_d-{\hat{f}}_d={\mathrm{\&Delta;f}}_0+{\stackrel{\&CenterDot;}{f}}_d\&CenterDot;{\mathrm{nT}}_s/2---\left(38\right)$

式中，Δf₀为初始多普勒频移估计偏移量。

将采样信号与本地信号进行相关运算，则有

$\begin{array}{c}{s}_n^r\&CenterDot;s_n^l={\mathrm{AR}}_n(\mathrm{\&tau;}-\widehat{\mathrm{\&tau;}})\exp j\left(2\mathrm{\&pi;}\right(f_d-{\hat{f}}_d){\mathrm{nT}}_s+{\mathrm{\&phi;}}_n-{\widehat{\mathrm{\&phi;}}}_n)\\ ={\mathrm{AR}}_n\left(\mathrm{\&Delta;}\mathrm{\&tau;}\right)\exp j\left(2\mathrm{\&pi;}\right({\mathrm{\&Delta;f}}_0+{\stackrel{\&CenterDot;}{f}}_d\&CenterDot;{\mathrm{nT}}_s/2){\mathrm{nT}}_s+{\mathrm{\&Delta;\&phi;}}_n)\\ ={\mathrm{AR}}_n\left(\mathrm{\&Delta;}\mathrm{\&tau;}\right)\exp j\left({\mathrm{\&Delta;\&phi;}}_n\right)\exp j\left(2{\mathrm{\&pi;\&Delta;f}}_0{\mathrm{nT}}_s\right)\exp j\&lsqb;\mathrm{\&pi;}{\stackrel{\&CenterDot;}{f}}_d{\left({\mathrm{nT}}_s\right)}^2\&rsqb;\end{array}---\left(39\right)$

进一步，假设相干累积时间为T_coh，对式(4)中的相关结果号进行相干累积，则其输出检测量为

$\begin{array}{c}{Z}_{COH}=\underset{n=0}{\overset{N-1}{\&Sigma;}}{\mathrm{AR}}_n\left(\mathrm{\&Delta;}\mathrm{\&tau;}\right)\exp j\left({\mathrm{\&Delta;\&phi;}}_n\right)\exp j\left(2{\mathrm{\&pi;\&Delta;f}}_0{\mathrm{nT}}_s\right)\exp j\&lsqb;\mathrm{\&pi;}{\stackrel{\&CenterDot;}{f}}_d{\left({\mathrm{nT}}_s\right)}^2\&rsqb;\\ ={\mathrm{AR}}_n\left(\mathrm{\&Delta;}\mathrm{\&tau;}\right)N\sin c\left({\mathrm{\&pi;\&Delta;f}}_0{T}_{coh}\right)\exp \left({\mathrm{\&pi;\&Delta;f}}_0\right(N-1)T_s+{\mathrm{\&Delta;\&phi;}}_n)\underset{n=2}{\overset{N-1}{\&Sigma;}}\exp j\&lsqb;\mathrm{\&pi;}{\stackrel{\&CenterDot;}{f}}_d{\left({\mathrm{nT}}_s\right)}^2\&rsqb;\end{array}---\left(40\right)$

式中，N＝T_coh/T_s，表示相干累积的采样点数。

由式(5)可知，当不存在多普勒频移率，即时，相干累积输出量可简化为

$\begin{array}{c}{Z}_{COH}=\underset{n=0}{\overset{N-1}{\&Sigma;}}{\mathrm{AR}}_n\left(\mathrm{\&Delta;}\mathrm{\&tau;}\right)\exp j\left({\mathrm{\&Delta;\&phi;}}_n\right)\exp j\left(2{\mathrm{\&pi;\&Delta;f}}_0{\mathrm{nT}}_s\right)\\ ={\mathrm{AR}}_n\left(\mathrm{\&Delta;}\mathrm{\&tau;}\right)N\sin c\left({\mathrm{\&pi;\&Delta;f}}_0{T}_{coh}\right)\exp \left({\mathrm{\&pi;\&Delta;f}}_0\right(N-1)T_s+{\mathrm{\&Delta;\&phi;}}_n)\end{array}---\left(41\right)$

此时，信号检测量的幅值仅由载波频率误差Δf₀和码相位误差Δτ两个变量决定，当Δf₀和Δτ均为零时，表示本地信号与接收信号的频率和相位均已对齐，此时输出检测量Z_COH即为捕获的相关峰。

而当多普勒频移率时，根据式(5)可知，的累积结果可表示为sinc函数，即Z_COH将沿着sinc函数主峰的下降梯度衰减。为了定量说明对检测量Z_COH幅值的损耗程度，将其定义为多普勒频移率损耗因子其中T为信号累积时间。对比式(5)和式(6)可知，相干累积算法的多普勒频移率损耗因子为

${\mathrm{\&alpha;}}_d({\stackrel{\&CenterDot;}{f}}_d,T_{coh})=\underset{n=0}{\overset{N-1}{\&Sigma;}}\exp j\&lsqb;\mathrm{\&pi;}{\stackrel{\&CenterDot;}{f}}_d{\left({\mathrm{nT}}_s\right)}^2\&rsqb;---\left(42\right)$

由于表现为sinc函数，因此多普勒频移率越大，幅度越小；同时随着累积时间T_coh的延长，幅度将进一步减小，导致信号累积能量的加速衰减。因此，为了在高动态环境中实现高灵敏度捕获，必须限制相干累积时间并进行多普勒频移率的补偿。

b.非相干、差分相干捕获性能

根据式(5)，可得相干累积后的同相、正交支路信号I_k和Q_k分别为

I_k＝AR_n(Δτ)sinc(πΔf₀T_coh)·cos(πΔf₀(N-1)T_s+Δφ_n)(43)

Q_k＝AR_n(Δτ)sinc(πΔf₀T_coh)·sin(πΔf₀(N-1)T_s+Δφ_n)(44)

而非相干和差分相干累积后得到的检测量可分别表示为

$Z_{NCH}=\underset{k=0}{\overset{T_{NCH}-1}{\&Sigma;}}(I_k^2+Q_k^2)---\left(45\right)$

$Z_{DFC}=\underset{k=2}{\overset{T_{DFC}+1}{\&Sigma;}}\left(\sqrt{{\left(I_k{I}_{k-1}{Q}_k{Q}_{k-1}\right)}^2+{(Q_k{I}_{k-1}-I_k{Q}_{k-1})}^2}\right)---\left(46\right)$

式中，Z_NCH，Z_DFC分别为非相干、差分相干输出矩阵；T_NCH和T_DFC分别为非相干和差分相干的累积周期。

将式(8)和式(9)分别代入式(10)和式(11)，得到非相干和差分相干的多普勒频移率衰减因子分别为

${\mathrm{\&alpha;}}_d^{NCH}({\stackrel{\&CenterDot;}{f}}_d,T_{NCH})=\underset{m=0}{\overset{T_{NCH}-1}{\&Sigma;}}{(\underset{n=0}{\overset{N-1}{\&Sigma;}}\exp \&lsqb;\mathrm{\&pi;}{\stackrel{\&CenterDot;}{f}}_d{(m+{\mathrm{nT}}_s)}^2\&rsqb;)}^2---\left(47\right)$

${\mathrm{\&alpha;}}_d^{DFC}({\stackrel{\&CenterDot;}{f}}_d,T_{DFC})=\underset{m=1}{\overset{T_{DFC}-1}{\&Sigma;}}(\underset{n=0}{\overset{N-1}{\&Sigma;}}\exp j\&lsqb;\mathrm{\&pi;}{\stackrel{\&CenterDot;}{f}}_d{(m-1+{\mathrm{nT}}_s)}^2\&rsqb;\&CenterDot;\underset{n=0}{\overset{N-1}{\&Sigma;}}\exp j\&lsqb;-\mathrm{\&pi;}{\stackrel{\&CenterDot;}{f}}_d{(m+{\mathrm{nT}}_s)}^2\&rsqb;)---\left(48\right)$

在低动态环境中，非相干和差分相干的累积增益包括相干增益及非相干或差分相干损耗；但是在高动态环境中，由于多普勒频移率的存在，本地信号与接收信号的频率/相位偏移量随着累积时间的延长不断增大，导致相应累积增益的衰落。因此，高动态累积增益的计算模型需引入多普勒频移率损耗因子则非相干和差分相干累积增益模型分别为

$G_{NCH}\left(T_{NCH}\right)=G_{COH}\left(T_{NCH}\right)-L_{NCH}\left(T_{NCH}\right)-10\lg \&lsqb;{\mathrm{\&alpha;}}_d^{NCH}({\stackrel{\&CenterDot;}{f}}_d,T_{NCH})\&rsqb;---\left(49\right)$

$G_{DFC}\left(T_{DFC}\right)=G_{COH}\left(T_{DFC}\right)-L_{DFC}\left(T_{DFC}\right)-l0\lg \&lsqb;{\mathrm{\&alpha;}}_d^{DFC}({\stackrel{\&CenterDot;}{f}}_d,T_{DFC})\&rsqb;---\left(50\right)$

式中，G_COH(·)为相干累积增益；L_NCH(·)和L_DFC(·)分别为非相干和差分相干损耗。采用累积周期分别为10ms、40ms、70ms以及100ms，得到非相干累积增益为16.0dB、12.5dB、10.0dB、8.1dB，差分相干累积增益为16.0dB、11.2dB、8.3dB、6.2dB。因此，高动态条件下，随着信号累积时间的延长，非相干累积方式的增益高于差分相干，故选择非相干方式对相干累积信号进行能量累积。

综上，组合累积方式选取相干和非相干累积方式进行组合。

步骤二：对输入中频信号进行10ms相干累积(其结果见式(5)所示)，并利用惯导系统和卫星星历对接收信号的多普勒参数及码相位进行估算，以压缩频率及相位搜索空间(具体实现结构见图2)；

a.多普勒参数估算

在信号捕获算法中引入惯导输出的载体位置、速度和加速度信息，结合卫星星历，则可得到由载体与卫星视线方向相对运动所产生的多普勒频移以及多普勒频移率，其计算式分别为

${\hat{f}}_d=-f_{L1}\frac{\left|\right|{\left({\hat{v}}_r\right)}_{LOS}\left|\right|}{c}=-f_{L1}\frac{({\hat{v}}_s-{\hat{v}}_u)\&CenterDot;{\hat{r}}_{LOS}}{c}---\left(51\right)$

${\widehat{\stackrel{\&CenterDot;}{f}}}_e=-f_{L1}\frac{\left|\right|{\left({\hat{a}}_r\right)}_{LOS}\left|\right|}{c}=-f_{L1}\frac{({\hat{a}}_s-{\hat{a}}_u)\&CenterDot;{\hat{r}}_{LOS}}{c}---\left(52\right)$

式中，为GPS星历提供的卫星速度和加速度矢量；分别为惯导系统输出的载体速度和加速度矢量；为卫星与载体间的单位视线矢量，可由卫星与载体间的相对位置计算得到。

将由式(16)计算得到的多普勒频移估算值作为频率搜索中心，频率搜索范围根据多普勒频移的计算精度决定，其中多普勒频移误差由星历误差和惯导系统误差两部分组成。设由星历计算所得的多普勒频移计算误差为(约为±150Hz)，惯导输出的速度信息误差为(若惯导系统在使用初期误差漂移较小，可认为估计误差在±5Hz以内)，则可以根据星历中给出的轨道误差项估算，而可通过下式估算

${\mathrm{\&Delta;}}_v^{SINS}=\underset{T_{COH}}{\&Integral;}g\left(\underset{T_{COH}}{\&Integral;}{C}_b^n\mathrm{\&epsiv;}dt\right)dt+\underset{T_{COH}}{\&Integral;}{C}_b^n\&dtri;dt---\left(53\right)$

式中，g为当地重力加速度；为本体系到导航系的姿态转换矩阵；ε和分别为陀螺仪和加速度计零偏。

因此，根据式(16)可得多普勒频移误差的方差为

${\mathrm{\&sigma;}}_{Dopp}^2={\left(\frac{f_{L1}}{c}\right)}^2{\hat{r}}_{LOS}{\left({\mathrm{\&Delta;}}_v^{SINS}\right)}^2{\hat{r}}_{LOS}^T+{\left({\mathrm{\&Delta;}}_{Dopp}^{Eph}\right)}^2---\left(54\right)$

将式(19)计算得到的多普勒频移误差最大值作为搜索边界以压缩搜索范围，从而达到提高捕获效率的目的。

b.码相位估算

由于GPS卫星信号传输链路较长，相比于发射码，接收码将产生τ时长的传输延迟；另外，C/A码速率受到多普勒频移的影响而不断变化，导致接收码相位延时的不确定性。因此，随着信号累积时间的延长，相比于发射码相位，接收码相位将出现大幅超前或滞后。故在高动态、低载噪比环境下，需要对接收码相位进行预先估算，以生成较为准确的本地复制码，提高本地码与接收码的初始对准精度，达到缩小码相位搜索区间、提高捕获效率的目的。

为了估算接收码相位，建立t_k采样时刻GPS接收信号的C/A码序列模型为

$C(t_k,\mathrm{\&tau;},f_d,{\stackrel{\&CenterDot;}{f}}_d)=C\left(\right(1+\mathrm{\&eta;}\left)\right(t_k-\mathrm{\&tau;}\left)\right)---\left(55\right)$

$\mathrm{\&eta;}=\frac{f_d+{\stackrel{\&CenterDot;}{f}}_d\&CenterDot;t_k/2}{f_{L1}}---\left(56\right)$

式中，η表示由多普勒效应引起的C/A码速率收缩因子。

若载体到卫星的相对位置矢量为r_s,u＝r_s-r_u，则GPS信号的传播时间为

${\hat{t}}_{prop}=\left|\right|{\hat{r}}_{s,u}\left|\right|/c=\left|\right|{\hat{r}}_s-{\hat{r}}_u\left|\right|/c---\left(57\right)$

式中，为GPS星历提供的卫星位置矢量；为惯导系统解算得到的载体位置矢量。

因此，接收码相位的估算模型为

$\widehat{\mathrm{\&tau;}}=rem\left(\right(t_{GPS}+t_{corr}-{\hat{t}}_{prop})\&CenterDot;f_{CA},\frac{c}{f_{CA}(1+\mathrm{\&eta;})})---\left(58\right)$

式中，t_GPS为当前C/A码的接收时刻(GPS时)；t_corr为卫星时钟修正量，可表示为

t_corr＝a_f0+a_f1(t_GPS-t_oc)+a_f2(t_GPS-t_oc)²(59)

式中，a_fi,i＝0,1,2为i阶卫星时钟修正系数，可从导航电文的第一子帧中获得；t_oc为第一子帧中第一数据块的参考时间。

码相位估算值误差包括载体-卫星相对位置误差及时间误差，可表示为

${\mathrm{\&sigma;}}_{\mathrm{\&tau;}}^2=4{\mathrm{\&sigma;}}_{pos}^2{\cos }^2\left({\mathrm{\&phi;}}_{el}\right)+{\mathrm{\&sigma;}}_t^2---\left(60\right)$

式中，φ_el为卫星仰角；Δ_eph为由卫星星历计算得到的GPS卫星定轨误差，可通过导航电文中的轨道误差项估算得到，和分别为惯导系统位置输出在纬度和经度方向上的误差，可通过下式估算

$\begin{array}{c}{\mathrm{\&sigma;}}_{pos}^2={\mathrm{\&Delta;}}_{eph}^2+{\mathrm{\&Delta;}}_{{\mathrm{SINS}}_L}^2+{\mathrm{\&Delta;}}_{{\mathrm{SINS}}_{\mathrm{\&lambda;}}}^2\\ {\mathrm{\&Delta;}}_{{\mathrm{SINS}}_L}=\frac{1}{R}\underset{T_{COH}}{\&Integral;}{\mathrm{\&Delta;}}_{v_y}^{SINS}\\ {\mathrm{\&Delta;}}_{{\mathrm{SINS}}_{\mathrm{\&lambda;}}}=\frac{\sec L}{R}\underset{T_{COH}}{\&Integral;}{\mathrm{\&Delta;}}_{v_x}^{SINS}\end{array}---\left(61\right)$

式中，和分别为导航坐标系x,y轴方向的惯导系统速度误差，可由式(18)计算获得。

根据式(23)得到码相位估算值，作为相位搜索空间中心；并根据式(25)计算相位误差的最大值，作为搜索边界，从而压缩搜索空间，提高高动态环境中的捕获效率。

步骤三：对多普勒参数及码相位估算值进行卡尔曼滤波平滑，并利用平滑结果计算本地信号与接收信号的动态频率/相位偏移量(具体实现结构见图2)；

为了精确补偿不同相干累积周期间本地信号与接收信号间的动态偏移量，需要保证累积矩阵循环平移的精度，因此预先平滑处理信号捕获参数，选取卡尔曼滤波器的状态量和量测量分别为

$X={\left(\begin{array}{ccc}\mathrm{\&tau;}& f_d& {\stackrel{\&CenterDot;}{f}}_d\end{array}\right)}^T---\left(62\right)$

$Z={\left(\begin{array}{ccc}{(r_s-r_u)}_i& {(v_s-v_u)}_i& {(a_s-a_u)}_i\end{array}\right)}_{9\&times;1}^T,i=x,y,z---\left(63\right)$

则根据式(16)、式(17)和式(23)可得状态方程和量测方程为

$\left[\begin{array}{c}\stackrel{\&CenterDot;}{\mathrm{\&tau;}}\\ {\stackrel{\&CenterDot;}{f}}_d\\ {\stackrel{\&CenterDot;\&CenterDot;}{f}}_d\end{array}\right]=\left[\begin{array}{ccc}0& \frac{f_{C/A}}{f_{L1}}& 0\\ 0& 0& 1\\ 0& 0& 0\end{array}\right]\left[\begin{array}{c}\mathrm{\&tau;}\\ f_d\\ {\stackrel{\&CenterDot;}{f}}_d\end{array}\right]+\left[\begin{array}{ccc}2{\mathrm{\&pi;f}}_{C/A}& 0& 0\\ 0& 2{\mathrm{\&pi;f}}_{L1}& 0\\ 0& 0& 2\mathrm{\&pi;}/{\mathrm{\&lambda;}}_{L1}\end{array}\right]\left[\begin{array}{c}{w}_b\\ w_d\\ w_a\end{array}\right]---\left(64\right)$

Z＝HX+V(65)

$H={\left[\begin{array}{ccc}{\left(r_{LOS}\right)}_{3\&times;1}^T& 0_{3\&times;1}& 0_{3\&times;1}\\ 0_{3\&times;1}& {\left(r_{LOS}\right)}_{3\&times;1}^T& 0_{3\&times;1}\\ 0_{3\&times;1}& 0_{3\&times;3}& {\left(r_{LOS}\right)}_{3\&times;1}^T\end{array}\right]}_{9\&times;3}---\left(66\right)$

$V=diag\left(\begin{array}{ccc}{\mathrm{\&sigma;}}_{\mathrm{\&tau;}}^2& {\mathrm{\&sigma;}}_{Dopp}^2& {\mathrm{\&sigma;}}_{Dopp\_Rate}^2\end{array}\right)---\left(67\right)$

式中，w_b为时钟偏差；w_d为时钟漂移；w_a为视线加速度的驱动噪声。

利用卡尔曼滤波器实时输出的信号参量平滑值，计算当前相干累积周期的累积矩阵输出与下一周期输出在载波频率搜索空间内的动态偏移量以及在码相位搜索空间内的动态偏移量其中

$\begin{array}{c}\mathrm{\&delta;}\hat{f}=(f_d-{\tilde{f}}_d)\&CenterDot;\mathrm{\&Delta;}t+{\stackrel{\&CenterDot;}{f}}_d\&CenterDot;{\mathrm{\&Delta;t}}^2/2\\ \mathrm{\&delta;}\widehat{\mathrm{\&tau;}}=(\mathrm{\&tau;}-\widetilde{\mathrm{\&tau;}})\&CenterDot;\mathrm{\&Delta;}t+\frac{f_{CA}}{f_{L1}}\&CenterDot;{\stackrel{\&CenterDot;}{f}}_d\&CenterDot;{\mathrm{\&Delta;t}}^2/2\end{array}---\left(68\right)$

式中，分别为当前搜索的多普勒频移值和码相位值；Δt为相干累积周期。

步骤四：见图3，采用循环平移算法补偿相干累积输出矩阵的动态偏移量；

对于第i个相干累积周期，由多普勒频移率导致的累积输出矩阵动态频率和相位偏移量分别为和在进行非相干累积之前，为了消除前i个周期的累积动态偏移量，第i个累积输出矩阵需沿频率搜索搜索空间循环平移个搜索单元，沿相位搜索搜索空间循环平移个搜索单元，若捕获搜索的频率分辨率和相位分辨率分别为f_Bin和τ_Bin，则和可表示为

$X_f^i=\&lsqb;\frac{\underset{n=1}{\overset{i}{\&Sigma;}}\mathrm{\&delta;}{\hat{f}}_i}{f_{Bin}}\&rsqb;,X_{\mathrm{\&tau;}}^i=\&lsqb;\frac{\underset{n=1}{\overset{i}{\&Sigma;}}\mathrm{\&delta;}{\widehat{\mathrm{\&tau;}}}_i}{{\mathrm{\&tau;}}_{Bin}}\&rsqb;---\left(69\right)$

式中，[·]为取整运算。

步骤五：对补偿后的输出矩阵进行10次非相干累积(其算法见式(12)所示)，将输出结果与门限比较，完成捕获判决(具体实现结构见图2)；

循环平移捕获周期内所有的相干输出矩阵后，对输出矩阵进行10次非相干累积，以进一步增强微弱信号的累积能量，提高载噪比处理增益，从而满足信号捕获最大峰值检测的灵敏度要求。

非相干累积输出的峰值为

$\begin{array}{c}{\mathrm{CP}}_{nch}=A^2{N}^2\underset{k=1}{\overset{10}{\&Sigma;}}{|R(f_{Bin}{X}_{\mathrm{\&tau;}}^k-\underset{k}{\overset{9}{\&Sigma;}}\mathrm{\&delta;}{\widehat{\mathrm{\&tau;}}}_k)\exp \&lsqb;{\mathrm{\&pi;NT}}_s({\mathrm{\&tau;}}_{Bin}{X}_f^k-\underset{k}{\overset{9}{\&Sigma;}}\mathrm{\&delta;}{\hat{f}}_k)\&rsqb;|}^2\\ =A^2{N}^2\underset{k=1}{\overset{10}{\&Sigma;}}{\left|\begin{array}{c}R({\mathrm{\&tau;}}_{Bin}\&CenterDot;{\&lsqb;\underset{k}{\overset{9}{\&Sigma;}}\mathrm{\&delta;}{\widehat{\mathrm{\&tau;}}}_k/{\mathrm{\&tau;}}_{Bin}\&rsqb;}_{round}-\underset{k}{\overset{9}{\&Sigma;}}\mathrm{\&delta;}{\widehat{\mathrm{\&tau;}}}_k)\&CenterDot;\\ \exp \&lsqb;{\mathrm{\&pi;NT}}_s(f_{Bin}\&CenterDot;{\&lsqb;\underset{k}{\overset{9}{\&Sigma;}}\mathrm{\&delta;}{\hat{f}}_k/f_{Bin}\&rsqb;}_{round}-\underset{k}{\overset{9}{\&Sigma;}}\mathrm{\&delta;}{\hat{f}}_k)\&rsqb;\end{array}\right|}^2\\ {\mathrm{\&omega;}}_{Bin}\&CenterDot;{\&lsqb;\underset{k}{\overset{9}{\&Sigma;}}\mathrm{\&delta;}{\widehat{\mathrm{\&omega;}}}_k/{\mathrm{\&omega;}}_{Bin}\&rsqb;}_{round}-\underset{k}{\overset{9}{\&Sigma;}}\mathrm{\&delta;}{\widehat{\mathrm{\&omega;}}}_k<<{\mathrm{\&omega;}}_{Bin}\mathrm{\&omega;}=\mathrm{\&tau;},f\end{array}---\left(70\right)$

若CP_nch＞2Threshold，则认为该颗卫星已成功捕获，其中，Threshold取非相干累积输出的第二峰值。

对于90％的检测概率，捕获处理信号的数据长度选取为100ms，相干累积10ms，非相干累积10次，则在高动态条件下，执行Monte-Carlo捕获仿真实验得到本发明方法的捕获灵敏度(有效捕获的最低载噪比)为29dB，相对于商用接收机的捕获灵敏度标准(35dB)提高了6dB。

Claims (1)

1.一种惯导系统辅助的高动态微弱信号GPS捕获方法，其特征在于：它具体包括以下步骤：

步骤一：对比相干、非相干以及差分相干累积捕获性能，选取累积方式进行组合，加强信号能量累积；

a.相干累积捕获性能

L1波段接收信号经过射频前端下变频及模数转换后得到数字中频信号，以采样频率f_s＝1/T_s进行采样，得到nT_s时刻采样信号以及本地相关信号模型分别为

$s_n^r={\mathrm{AC}}_n\left(\right(1+\mathrm{\&eta;}\left)\right({\mathrm{nT}}_s-\mathrm{\&tau;}\left)\right)D_n\left({\mathrm{nT}}_s\right)\exp \left(2\mathrm{\&pi;}\right(f_{IF}+f_d){\mathrm{nT}}_s+{\mathrm{\&phi;}}_0)+W_n---\left(1\right)$

$s_n^l=C_n\left(\right(1+\mathrm{\&eta;}\left)\right({\mathrm{nT}}_s-\widehat{\mathrm{\&tau;}}\left)\right)\exp \left(2\mathrm{\&pi;}\right(f_{IF}+\hat{f}){\mathrm{nT}}_s+{\widehat{\mathrm{\&phi;}}}_0)---\left(2\right)$

式中，下标n表示当前采样点；上标r表示采样信号，上标l表示本地信号；C(·)为粗捕获码序列；D(·)为导航数据；f_IF为载波中频频率；τ和分别为采样信号的真实码相位传播延迟以及本地信号的码相位估计值；f_d和分别为采样信号的真实多普勒频移和本地信号多普勒频移估计值；φ₀和分别为采样信号初始载波相位和本地信号初始载波相位估计值；多普勒频移f_d同时使C/A码速率产生偏移，偏移因子W_n为高斯白噪声；

假设由载体与GPS卫星相对加速度引起的多普勒频移率为则真实多普勒频移与多普勒频移估算值之间的关系为

$f_d-{\hat{f}}_d={\mathrm{\&Delta;f}}_0+{\stackrel{\&CenterDot;}{f}}_d\&CenterDot;{\mathrm{nT}}_s/2---\left(3\right)$

式中，Δf₀为初始多普勒频移估计偏移量；

将采样信号与本地信号进行相关运算，则有

$\begin{array}{c}{s}_n^r\&CenterDot;s_n^l={\mathrm{AR}}_n\left(\mathrm{\&tau;}-\widehat{\mathrm{\&tau;}}\right)\exp j\left(2\mathrm{\&pi;}\left(f_d-{\hat{f}}_d\right){\mathrm{nT}}_s+{\mathrm{\&phi;}}_n-{\widehat{\mathrm{\&phi;}}}_n\right)\\ ={\mathrm{AR}}_n\left(\mathrm{\&Delta;}\mathrm{\&tau;}\right)\exp j\left(2\mathrm{\&pi;}\left({\mathrm{\&Delta;f}}_0+{\stackrel{\&CenterDot;}{f}}_d\&CenterDot;{\mathrm{nT}}_s/2\right){\mathrm{nT}}_s+{\mathrm{\&Delta;\&phi;}}_n\right)\\ ={\mathrm{AR}}_n\left(\mathrm{\&Delta;}\mathrm{\&tau;}\right)\exp j\left({\mathrm{\&Delta;\&phi;}}_n\right)\exp j\left(2{\mathrm{\&pi;\&Delta;f}}_0{\mathrm{nT}}_s\right)\exp j\&lsqb;\mathrm{\&pi;}{\stackrel{\&CenterDot;}{f}}_d{\left({\mathrm{nT}}_s\right)}^2\&rsqb;\end{array}---\left(4\right)$

进一步，假设相干累积时间为T_coh，对式(4)中的相关结果号进行相干累积，则其输出检测量为

$\begin{array}{c}{Z}_{COH}=\underset{n=0}{\overset{N-1}{\&Sigma;}}{\mathrm{AR}}_n\left(\mathrm{\&Delta;}\mathrm{\&tau;}\right)\exp j\left({\mathrm{\&Delta;\&phi;}}_n\right)\exp j\left(2{\mathrm{\&pi;\&Delta;f}}_0{\mathrm{nT}}_s\right)\exp j\&lsqb;\mathrm{\&pi;}{\stackrel{\&CenterDot;}{f}}_d{\left({\mathrm{nT}}_s\right)}^2\&rsqb;\\ ={\mathrm{AR}}_n\left(\mathrm{\&Delta;}\mathrm{\&tau;}\right)N\sin c\left({\mathrm{\&pi;\&Delta;f}}_0{T}_{coh}\right)\exp \left({\mathrm{\&pi;\&Delta;f}}_0\left(N-1\right)T_s+{\mathrm{\&Delta;\&phi;}}_n\right)\underset{n=0}{\overset{N-1}{\&Sigma;}}\exp j\&lsqb;\mathrm{\&pi;}{\stackrel{\&CenterDot;}{f}}_d{\left({\mathrm{nT}}_s\right)}^2\&rsqb;\end{array}---\left(5\right)$

式中，N＝T_coh/T_s，表示相干累积的采样点数；

由式(5)知，当不存在多普勒频移率，即时，相干累积输出量简化为

$\begin{array}{c}{Z}_{COH}=\underset{n=0}{\overset{N-1}{\&Sigma;}}{\mathrm{AR}}_n\left(\mathrm{\&Delta;}\mathrm{\&tau;}\right)\exp j\left({\mathrm{\&Delta;\&phi;}}_n\right)\exp j\left(2{\mathrm{\&pi;\&Delta;f}}_0{\mathrm{nT}}_s\right)\\ ={\mathrm{AR}}_n\left(\mathrm{\&Delta;}\mathrm{\&tau;}\right)N\sin c\left({\mathrm{\&pi;\&Delta;f}}_0{T}_{coh}\right)\exp \left({\mathrm{\&pi;\&Delta;f}}_0\left(N-1\right)T_s+{\mathrm{\&Delta;\&phi;}}_n\right)\end{array}---\left(6\right)$

此时，信号检测量的幅值仅由载波频率误差Δf₀和码相位误差Δτ两个变量决定，当Δf₀和Δτ均为零时，表示本地信号与接收信号的频率和相位均已对齐，此时输出检测量Z_COH即为捕获的相关峰；

而当多普勒频移率时，根据式(5)知，的累积结果表示为sinc函数，即Z_COH将沿着sinc函数主峰的下降梯度衰减；为了定量说明对检测量Z_COH幅值的损耗程度，将其定义为多普勒频移率损耗因子其中T为信号累积时间；对比式(5)和式(6)知，相干累积算法的多普勒频移率损耗因子为

${\mathrm{\&alpha;}}_d({\stackrel{\&CenterDot;}{f}}_d,T_{coh})=\underset{n=0}{\overset{N-1}{\&Sigma;}}\exp j\&lsqb;\mathrm{\&pi;}{\stackrel{\&CenterDot;}{f}}_d{\left({\mathrm{nT}}_s\right)}^2\&rsqb;---\left(7\right)$

由于表现为sinc函数，因此多普勒频移率越大，幅度越小；同时随着累积时间T_coh的延长，幅度将进一步减小，导致信号累积能量的加速衰减，因此，为了在高动态环境中实现高灵敏度捕获，必须限制相干累积时间并进行多普勒频移率的补偿；

b.非相干、差分相干捕获性能

根据式(5)，得相干累积后的同相、正交支路信号I_k和Q_k分别为

I_k＝AR_n(Δτ)sinc(πΔf₀T_coh)·cos(πΔf₀(N-1)T_s+Δφ_n)(8)

Q_k＝AR_n(Δτ)sinc(πΔf₀T_coh)·sin(πΔf₀(N-1)T_s+Δφ_n)(9)

而非相干和差分相干累积后得到的检测量分别表示为

$Z_{NCH}=\underset{k=0}{\overset{T_{NCH}-1}{\&Sigma;}}\left(I_k^2+Q_k^2\right)---\left(10\right)$

$Z_{DFC}=\underset{k=2}{\overset{T_{DFC}+1}{\&Sigma;}}\left(\sqrt{{(I_k{I}_{k-1}+Q_k{Q}_{k-1})}^2+{(Q_k{I}_{k-1}-I_k{Q}_{k-1})}^2}\right)---\left(11\right)$

式中，Z_NCH，Z_DFC分别为非相干、差分相干输出矩阵；T_NCH和T_DFC分别为非相干和差分相干的累积周期；

将式(8)和式(9)分别代入式(10)和式(11)，得到非相干和差分相干的多普勒频移率衰减因子分别为

${\mathrm{\&alpha;}}_d^{NCH}({\stackrel{\&CenterDot;}{f}}_d,T_{NCH})=\underset{m=0}{\overset{T_{NCH}-1}{\&Sigma;}}{(\underset{n=0}{\overset{N-1}{\&Sigma;}}\exp j\&lsqb;\mathrm{\&pi;}{\stackrel{\&CenterDot;}{f}}_d{(m+{\mathrm{nT}}_s)}^2\&rsqb;)}^2---\left(12\right)$

${\mathrm{\&alpha;}}_d^{DFC}({\stackrel{\&CenterDot;}{f}}_d,T_{DFC})=\underset{m=1}{\overset{T_{DFC}-1}{\&Sigma;}}(\underset{n=0}{\overset{N-1}{\&Sigma;}}\exp j\&lsqb;\mathrm{\&pi;}{\stackrel{\&CenterDot;}{f}}_d{(m-1+{\mathrm{nT}}_s)}^2\&rsqb;\&CenterDot;\underset{n=0}{\overset{N-1}{\&Sigma;}}\exp j\&lsqb;-\mathrm{\&pi;}{\stackrel{\&CenterDot;}{f}}_d{(m+{\mathrm{nT}}_s)}^2\&rsqb;)---\left(13\right)$

在低动态环境中，非相干和差分相干的累积增益包括相干增益及非相干或差分相干损耗；但是在高动态环境中，由于多普勒频移率的存在，本地信号与接收信号的频率/相位偏移量随着累积时间的延长不断增大，导致相应累积增益的衰落；因此，高动态累积增益的计算模型需引入多普勒频移率损耗因子则非相干和差分相干累积增益模型分别为

$G_{NCH}\left(T_{NCH}\right)=G_{COH}\left(T_{NCH}\right)-L_{NCH}\left(T_{NCH}\right)-10\lg \&lsqb;{\mathrm{\&alpha;}}_d^{NCH}({\stackrel{\&CenterDot;}{f}}_d,T_{NCH})\&rsqb;---\left(14\right)$

$G_{DFC}\left(T_{DFC}\right)=G_{COH}\left(T_{DFC}\right)-L_{DFC}\left(T_{DFC}\right)-10\lg \&lsqb;{\mathrm{\&alpha;}}_d^{DFC}({\stackrel{\&CenterDot;}{f}}_d,T_{DFC})\&rsqb;---\left(15\right)$

式中，G_COH(·)为相干累积增益；L_NCH(·)和L_DFC(·)分别为非相干和差分相干损耗，通过数值仿真可得，高动态条件下，随着信号累积时间的延长，非相干和差分相干的累积增益均不断下降，但是相比于非相干累积方式，差分相干增益的下降速度更快，因此，为了在高动态、低载噪比环境中获得较高的累积增益，应选择非相干方式对相干累积信号进行能量累积；综上，组合累积方式选取相干和非相干累积方式进行组合；

步骤二：对输入中频信号进行相干累积，结果见式(5)所示，并利用惯导系统和卫星星历对接收信号的多普勒参数及码相位进行估算，以压缩频率及相位搜索空间；

a.多普勒参数估算

在信号捕获算法中引入惯导输出的载体位置、速度和加速度信息，结合卫星星历，则得到由载体与卫星视线方向相对运动所产生的多普勒频移以及多普勒频移率，其计算式分别为

${\hat{f}}_d=-f_{L1}\frac{\left|\right|{\left({\hat{v}}_r\right)}_{LOS}\left|\right|}{c}=-f_{L1}\frac{\left({\hat{v}}_s-{\hat{v}}_u\right)\&CenterDot;{\hat{r}}_{LOS}}{c}---\left(16\right)$

${\widehat{\stackrel{\&CenterDot;}{f}}}_d=-f_{L1}\frac{\left|\right|{\left({\hat{a}}_r\right)}_{LOS}\left|\right|}{c}=-f_{L1}\frac{({\hat{a}}_s-{\hat{a}}_u)\&CenterDot;{\hat{r}}_{LOS}}{c}---\left(17\right)$

式中，为GPS星历提供的卫星速度和加速度矢量；分别为惯导系统输出的载体速度和加速度矢量；为卫星与载体间的单位视线矢量，由卫星与载体间的相对位置计算得到；

将由式(16)计算得到的多普勒频移估算值作为频率搜索中心，频率搜索范围根据多普勒频移的计算精度决定，其中多普勒频移误差由星历误差和惯导系统误差两部分组成；设由星历计算所得的多普勒频移计算误差为惯导输出的速度信息误差为则根据星历中给出的轨道误差项估算，而通过下式估算

${\mathrm{\&Delta;}}_v^{SINS}=\underset{T_{COH}}{\&Integral;}g\left(\underset{T_{COH}}{\&Integral;}{C}_b^n\mathrm{\&epsiv;}dt\right)dt+\underset{T_{COH}}{\&Integral;}{C}_b^n\&dtri;dt---\left(18\right)$

式中，g为当地重力加速度；为本体系到导航系的姿态转换矩阵；ε和▽分别为陀螺仪和加速度计零偏；

因此，根据式(16)得多普勒频移误差的方差为

${\mathrm{\&sigma;}}_{Dopp}^2={\left(\frac{f_{L1}}{c}\right)}^2{\hat{r}}_{LOS}{\left({\mathrm{\&Delta;}}_v^{SINS}\right)}^2{\hat{r}}_{LOS}^T+{\left({\mathrm{\&Delta;}}_{Dopp}^{Eph}\right)}^2---\left(19\right)$

将式(19)计算得到的多普勒频移误差最大值作为搜索边界以压缩搜索范围，从而达到提高捕获效率的目的；

b.码相位估算

由于GPS卫星信号传输链路较长，相比于发射码，接收码将产生τ时长的传输延迟；另外，C/A码速率受到多普勒频移的影响而不断变化，导致接收码相位延时的不确定性，因此，随着信号累积时间的延长，相比于发射码相位，接收码相位将出现大幅超前或滞后，故在高动态、低载噪比环境下，需要对接收码相位进行预先估算，以生成准确的本地复制码，提高本地码与接收码的初始对准精度，达到缩小码相位搜索区间、提高捕获效率的目的；

为了估算接收码相位，建立t_k采样时刻GPS接收信号的C/A码序列模型为

$C(t_k,\mathrm{\&tau;},f_d,{\stackrel{\&CenterDot;}{f}}_d)=C\left(\right(1+\mathrm{\&eta;}\left)\right(t_k-\mathrm{\&tau;}\left)\right)---\left(20\right)$

$\mathrm{\&eta;}=\frac{f_d+{\stackrel{\&CenterDot;}{f}}_d\&CenterDot;t_k/2}{f_{L1}}---\left(21\right)$

式中，η表示由多普勒效应引起的C/A码速率收缩因子；

若载体到卫星的相对位置矢量为r_s,u＝r_s-r_u，则GPS信号的传播时间为

${\hat{t}}_{prop}=\left|\right|{\hat{r}}_{s,u}\left|\right|/c=\left|\right|{\hat{r}}_s-{\hat{r}}_u\left|\right|/c---\left(22\right)$

式中，为GPS星历提供的卫星位置矢量；为惯导系统解算得到的载体位置矢量；

因此，接收码相位的估算模型为

$\widehat{\mathrm{\&tau;}}=rem\left(\right(t_{GPS}+t_{corr}-{\hat{t}}_{prop})\&CenterDot;f_{CA},\frac{c}{f_{CA}(1+\mathrm{\&eta;})})---\left(23\right)$

式中，t_GPS为当前C/A码的接收时刻即GPS时；t_corr为卫星时钟修正量，表示为

t_corr＝a_f0+a_f1(t_GPS-t_oc)+a_f2(t_GPS-t_oc)²(24)

式中，a_fi,i＝0,1,2为i阶卫星时钟修正系数，从导航电文的第一子帧中获得；t_oc为第一子帧中第一数据块的参考时间；

码相位估算值误差包括载体-卫星相对位置误差及时间误差，表示为

${\mathrm{\&sigma;}}_{\mathrm{\&tau;}}^2=4{\mathrm{\&sigma;}}_{pos}^2{\cos }^2\left({\mathrm{\&phi;}}_{el}\right)+{\mathrm{\&sigma;}}_t^2---\left(25\right)$

式中，φ_el为卫星仰角；Δ_eph为由卫星星历计算得到的GPS卫星定轨误差，通过导航电文中的轨道误差项估算得到，和分别为惯导系统位置输出在纬度和经度方向上的误差，通过下式估算

${\mathrm{\&sigma;}}_{pos}^2={\mathrm{\&Delta;}}_{eph}^2+{\mathrm{\&Delta;}}_{{\mathrm{SINS}}_L}^2+{\mathrm{\&Delta;}}_{{\mathrm{SINS}}_{\mathrm{\&lambda;}}}^2$

${\mathrm{\&Delta;}}_{{\mathrm{SINS}}_L}=\frac{1}{R}\underset{T_{COH}}{\&Integral;}{\mathrm{\&Delta;}}_{v_y}^{SINS}$

(26)

${\mathrm{\&Delta;}}_{{\mathrm{SINS}}_{\mathrm{\&lambda;}}}=\frac{\sec L}{R}\underset{T_{COH}}{\&Integral;}{\mathrm{\&Delta;}}_{v_x}^{SINS}$

式中，和分别为导航坐标系x,y轴方向的惯导系统速度误差，由式(18)计算获得；

根据式(23)得到码相位估算值，作为相位搜索空间中心；并根据式(25)计算相位误差的最大值，作为搜索边界，从而压缩搜索空间，提高高动态环境中的捕获效率；

步骤三：对多普勒参数及码相位估算值进行卡尔曼滤波平滑，并利用平滑结果计算本地信号与接收信号的动态频率/相位偏移量；

为了精确补偿不同相干累积周期间本地信号与接收信号间的动态偏移量，需要保证累积矩阵循环平移的精度，因此预先平滑处理信号捕获参数，选取卡尔曼滤波器的状态量和量测量分别为

$X={\left(\begin{array}{ccc}\mathrm{\&tau;}& f_d& {\stackrel{\&CenterDot;}{f}}_d\end{array}\right)}^T---\left(27\right)$

$Z={\left(\begin{array}{ccc}{\left(r_s-r_u\right)}_i& {\left(v_s-v_u\right)}_i& {\left(a_s-a_u\right)}_i\end{array}\right)}_{9\&times;1}^T,i=x,y,z---\left(28\right)$

则根据式(16)、式(17)和式(23)得状态方程和量测方程为

$\left[\begin{array}{c}\stackrel{\&CenterDot;}{\mathrm{\&tau;}}\\ {\stackrel{\&CenterDot;}{f}}_d\\ {\stackrel{\&CenterDot;\&CenterDot;}{f}}_d\end{array}\right]=\left[\begin{array}{ccc}0& \frac{f_{C/A}}{f_{L1}}& 0\\ 0& 0& 1\\ 0& 0& 0\end{array}\right]\left[\begin{array}{c}\mathrm{\&tau;}\\ f_d\\ {\stackrel{\&CenterDot;}{f}}_d\end{array}\right]+\left[\begin{array}{ccc}2{\mathrm{\&pi;f}}_{C/A}& 0& 0\\ 0& 2{\mathrm{\&pi;f}}_{L1}& 0\\ 0& 0& 2\mathrm{\&pi;}/{\mathrm{\&lambda;}}_{L1}\end{array}\right]\left[\begin{array}{c}{w}_b\\ w_d\\ w_a\end{array}\right]---\left(29\right)$

Z＝HX+V(30)

$H={\left[\begin{array}{ccc}{\left(r_{LOS}\right)}_{3\&times;1}^T& 0_{3\&times;1}& 0_{3\&times;1}\\ 0_{3\&times;1}& {\left(r_{LOS}\right)}_{3\&times;1}^T& 0_{3\&times;1}\\ 0_{3\&times;1}& 0_{3\&times;3}& {\left(r_{LOS}\right)}_{3\&times;1}^T\end{array}\right]}_{9\&times;3}---\left(31\right)$

$V=diag\left(\begin{array}{ccc}{\mathrm{\&sigma;}}_{\mathrm{\&tau;}}^2& {\mathrm{\&sigma;}}_{Dopp}^2& {\mathrm{\&sigma;}}_{Dopp\_Rate}^2\end{array}\right)---\left(32\right)$

式中，w_b为时钟偏差；w_d为时钟漂移；w_a为视线加速度的驱动噪声；

利用卡尔曼滤波器实时输出的信号参量平滑值，计算当前相干累积周期的累积矩阵输出与下一周期输出在载波频率搜索空间内的动态偏移量以及在码相位搜索空间内的动态偏移量其中

$\mathrm{\&delta;}\hat{f}=(f_d-{\tilde{f}}_d)\&CenterDot;\mathrm{\&Delta;}t+{\stackrel{\&CenterDot;}{f}}_d\&CenterDot;{\mathrm{\&Delta;t}}^2/2$

$\mathrm{\&delta;}\widehat{\mathrm{\&tau;}}=(\mathrm{\&tau;}-\widetilde{\mathrm{\&tau;}})\&CenterDot;\mathrm{\&Delta;}t+\frac{f_{CA}}{f_{L1}}\&CenterDot;{\stackrel{\&CenterDot;}{f}}_d\&CenterDot;{\mathrm{\&Delta;t}}^2/2---\left(33\right)$

式中，分别为当前搜索的多普勒频移值和码相位值；Δt为相干累积周期；

步骤四：采用循环平移算法补偿相干累积输出矩阵的动态偏移量；

对于第i个相干累积周期，由多普勒频移率导致的累积输出矩阵动态频率和相位偏移量分别为和在进行非相干累积之前，为了消除前i个周期的累积动态偏移量，第i个累积输出矩阵需沿频率搜索搜索空间循环平移个搜索单元，沿相位搜索搜索空间循环平移个搜索单元，若捕获搜索的频率分辨率和相位分辨率分别为f_Bin和τ_Bin，则和表示为

$X_f^i=\&lsqb;\frac{\underset{n=1}{\overset{i}{\&Sigma;}}\mathrm{\&delta;}{\hat{f}}_i}{f_{Bin}}\&rsqb;,X_{\mathrm{\&tau;}}^i=\&lsqb;\frac{\underset{n=1}{\overset{i}{\&Sigma;}}\mathrm{\&delta;}{\widehat{\mathrm{\&tau;}}}_i}{{\mathrm{\&tau;}}_{Bin}}\&rsqb;---\left(34\right)$

式中，[·]为取整运算；

步骤五：对补偿后的输出矩阵进行非相干累积，其算法见式(12)所示，将输出结果与门限比较，完成捕获判决；

循环平移捕获周期内所有的相干输出矩阵后，对输出矩阵进行非相干累积，以进一步增强微弱信号的累积能量，提高载噪比处理增益，从而满足信号捕获最大峰值检测的灵敏度要求；

非相干累积输出的峰值为

$\begin{array}{c}{\mathrm{CP}}_{nch}=A^2{N}^2\underset{k=1}{\overset{n}{\&Sigma;}}|R\left(f_{Bin}{X}_f^i-\underset{k=i}{\overset{n-1}{\&Sigma;}}\mathrm{\&delta;}{\widehat{\mathrm{\&tau;}}}_k\right)\exp \&lsqb;{\mathrm{\&pi;NT}}_s\left({\mathrm{\&tau;}}_{Bin}{X}_{\mathrm{\&tau;}}^i-\underset{k=i}{\overset{n-1}{\&Sigma;}}\mathrm{\&delta;}{\hat{f}}_k\right)\&rsqb;{|}^2\\ =A^2{N}^2\underset{k=1}{\overset{n}{\&Sigma;}}|\begin{array}{c}R\left({\mathrm{\&tau;}}_{Bin}\&CenterDot;{\&lsqb;\underset{k=i}{\overset{n-1}{\&Sigma;}}\mathrm{\&delta;}{\widehat{\mathrm{\&tau;}}}_k/{\mathrm{\&tau;}}_{Bin}\&rsqb;}_{round}-\underset{k=i}{\overset{n-1}{\&Sigma;}}\mathrm{\&delta;}{\widehat{\mathrm{\&tau;}}}_k\right)\&CenterDot;\\ \exp \&lsqb;{\mathrm{\&pi;NT}}_s\left(f_{Bin}\&CenterDot;{\&lsqb;\underset{k=i}{\overset{n-1}{\&Sigma;}}\mathrm{\&delta;}{\hat{f}}_k/f_{Bin}\&rsqb;}_{round}-\underset{k=i}{\overset{n-1}{\&Sigma;}}\mathrm{\&delta;}{\hat{f}}_k\right)\&rsqb;\end{array}{|}^2\end{array}---\left(35\right)$

${\mathrm{\&omega;}}_{Bin}\&CenterDot;{\&lsqb;\underset{k=i}{\overset{n-1}{\&Sigma;}}\mathrm{\&delta;}{\widehat{\mathrm{\&omega;}}}_k/{\mathrm{\&omega;}}_{Bin}\&rsqb;}_{round}-\underset{k=i}{\overset{n-1}{\&Sigma;}}\mathrm{\&delta;}{\widehat{\mathrm{\&omega;}}}_k<<{\mathrm{\&omega;}}_{Bin}\mathrm{\&omega;}=\mathrm{\&tau;},f$

若CP_nch＞2Threshold，则认为该颗卫星已成功捕获，其中，Threshold取非相干累积输出的第二峰值。