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

摘要

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

主权项

一种惯导系统辅助的高动态微弱信号GPS捕获方法，其特征在于：它具体包括以下步骤：步骤一：对比相干、非相干以及差分相干累积捕获性能，选取累积方式进行组合，加强信号能量累积；a.相干累积捕获性能L1波段接收信号经过射频前端下变频及模数转换后得到数字中频信号，以采样频率f_s＝1/T_s进行采样，得到nT_s时刻采样信号以及本地相关信号模型分别为$s_n^r={\mathrm{AC}}_n\left(\right(1+\eta \left)\right({\mathrm{nT}}_s-\tau \left)\right)D_n\left({\mathrm{nT}}_s\right)\exp \left(2\pi \right(f_{IF}+f_d){\mathrm{nT}}_s+{\phi }_0)+W_n---\left(1\right)$$s_n^l=C_n\left(\right(1+\eta \left)\right({\mathrm{nT}}_s-\hat{\tau }\left)\right)\exp \left(2\pi \right(f_{IF}+\hat{f}){\mathrm{nT}}_s+{\hat{\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{·}{f}}_d·{\mathrm{nT}}_s/2---\left(3\right)$式中，Δf₀为初始多普勒频移估计偏移量；将采样信号与本地信号进行相关运算，则有$\begin{array}{c}{s}_n^r·s_n^l={\mathrm{AR}}_n\left(\tau -\hat{\tau }\right)\exp j\left(2\pi \left(f_d-{\hat{f}}_d\right){\mathrm{nT}}_s+{\phi }_n-{\hat{\phi }}_n\right)\\ ={\mathrm{AR}}_n\left(\Delta \tau \right)\exp j\left(2\pi \left({\mathrm{\Delta f}}_0+{\stackrel{·}{f}}_d·{\mathrm{nT}}_s/2\right){\mathrm{nT}}_s+{\mathrm{\Delta \phi }}_n\right)\\ ={\mathrm{AR}}_n\left(\Delta \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\left[\pi {\stackrel{·}{f}}_d{\left({\mathrm{nT}}_s\right)}^2\right]\end{array}---\left(4\right)$进一步，假设相干累积时间为T_coh，对式(4)中的相关结果号进行相干累积，则其输出检测量为$\begin{array}{c}{Z}_{COH}=\underset{n=0}{\overset{N-1}{\Sigma }}{\mathrm{AR}}_n\left(\Delta \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\left[\pi {\stackrel{·}{f}}_d{\left({\mathrm{nT}}_s\right)}^2\right]\\ ={\mathrm{AR}}_n\left(\Delta \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\left[\pi {\stackrel{·}{f}}_d{\left({\mathrm{nT}}_s\right)}^2\right]\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(\Delta \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(\Delta \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)知，相干累积算法的多普勒频移率损耗因子为${\alpha }_d({\stackrel{·}{f}}_d,T_{coh})=\underset{n=0}{\overset{N-1}{\Sigma }}\exp j\left[\pi {\stackrel{·}{f}}_d{\left({\mathrm{nT}}_s\right)}^2\right]---\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)，得到非相干和差分相干的多普勒频移率衰减因子分别为${\alpha }_d^{NCH}({\stackrel{·}{f}}_d,T_{NCH})=\underset{m=0}{\overset{T_{NCH}-1}{\Sigma }}{\left(\underset{n=0}{\overset{N-1}{\Sigma }}\exp j\right[\pi {\stackrel{·}{f}}_d{(m+{\mathrm{nT}}_s)}^2\left]\right)}^2---\left(12\right)$${\alpha }_d^{DFC}({\stackrel{·}{f}}_d,T_{DFC})=\underset{m=1}{\overset{T_{DFC}-1}{\Sigma }}\left(\underset{n=0}{\overset{N-1}{\Sigma }}\exp j\right[\pi {\stackrel{·}{f}}_d{(m-1+{\mathrm{nT}}_s)}^2]·\underset{n=0}{\overset{N-1}{\Sigma }}\exp j[-\pi {\stackrel{·}{f}}_d{(m+{\mathrm{nT}}_s)}^2\left]\right)---\left(13\right)$在低动态环境中，非相干和差分相干的累积增益包括相干增益及非相干或差分相干损耗；但是在高动态环境中，由于多普勒频移率的存在，本地信号与接收信号的频率/相位偏移量随着累积时间的延长不断增大，导致相应累积增益的衰落；因此，高动态累积增益的计算模型需引入多普勒频移率损耗因子则非相干和差分相干累积增益模型分别为$G_{NCH}\left(T_{NCH}\right)=G_{COH}\left(T_{NCH}\right)-L_{NCH}\left(T_{NCH}\right)-10\lg \left[{\alpha }_d^{NCH}({\stackrel{·}{f}}_d,T_{NCH})\right]---\left(14\right)$$G_{DFC}\left(T_{DFC}\right)=G_{COH}\left(T_{DFC}\right)-L_{DFC}\left(T_{DFC}\right)-10\lg \left[{\alpha }_d^{DFC}({\stackrel{·}{f}}_d,T_{DFC})\right]---\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)·{\hat{r}}_{LOS}}{c}---\left(16\right)$${\widehat{\stackrel{·}{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)·{\hat{r}}_{LOS}}{c}---\left(17\right)$式中，为GPS星历提供的卫星速度和加速度矢量；分别为惯导系统输出的载体速度和加速度矢量；为卫星与载体间的单位视线矢量，由卫星与载体间的相对位置计算得到；将由式(16)计算得到的多普勒频移估算值作为频率搜索中心，频率搜索范围根据多普勒频移的计算精度决定，其中多普勒频移误差由星历误差和惯导系统误差两部分组成；设由星历计算所得的多普勒频移计算误差为惯导输出的速度信息误差为则根据星历中给出的轨道误差项估算，而通过下式估算${\Delta }_v^{SINS}=\underset{T_{COH}}{\int }g\left(\underset{T_{COH}}{\int }{C}_b^nϵdt\right)dt+\underset{T_{COH}}{\int }{C}_b^n▿dt---\left(18\right)$式中，g为当地重力加速度；为本体系到导航系的姿态转换矩阵；ε和▽分别为陀螺仪和加速度计零偏；因此，根据式(16)得多普勒频移误差的方差为${\sigma }_{Dopp}^2={\left(\frac{f_{L1}}{c}\right)}^2{\hat{r}}_{LOS}{\left({\Delta }_v^{SINS}\right)}^2{\hat{r}}_{LOS}^T+{\left({\Delta }_{Dopp}^{Eph}\right)}^2---\left(19\right)$将式(19)计算得到的多普勒频移误差最大值作为搜索边界以压缩搜索范围，从而达到提高捕获效率的目的；b.码相位估算由于GPS卫星信号传输链路较长，相比于发射码，接收码将产生τ时长的传输延迟；另外，C/A码速率受到多普勒频移的影响而不断变化，导致接收码相位延时的不确定性，因此，随着信号累积时间的延长，相比于发射码相位，接收码相位将出现大幅超前或滞后，故在高动态、低载噪比环境下，需要对接收码相位进行预先估算，以生成准确的本地复制码，提高本地码与接收码的初始对准精度，达到缩小码相位搜索区间、提高捕获效率的目的；为了估算接收码相位，建立t_k采样时刻GPS接收信号的C/A码序列模型为$C(t_k,\tau ,f_d,{\stackrel{·}{f}}_d)=C\left(\right(1+\eta \left)\right(t_k-\tau \left)\right)---\left(20\right)$$\eta =\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星历提供的卫星位置矢量；为惯导系统解算得到的载体位置矢量；因此，接收码相位的估算模型为$\hat{\tau }=rem\left(\right(t_{GPS}+t_{corr}-{\hat{t}}_{prop})·f_{CA},\frac{c}{f_{CA}(1+\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为第一子帧中第一数据块的参考时间；码相位估算值误差包括载体‑卫星相对位置误差及时间误差，表示为${\sigma }_{\tau }^2=4{\sigma }_{pos}^2{\cos }^2\left({\phi }_{el}\right)+{\sigma }_t^2---\left(25\right)$式中，φ_el为卫星仰角；Δ_eph为由卫星星历计算得到的GPS卫星定轨误差，通过导航电文中的轨道误差项估算得到，和分别为惯导系统位置输出在纬度和经度方向上的误差，通过下式估算${\sigma }_{pos}^2={\Delta }_{eph}^2+{\Delta }_{{\mathrm{SINS}}_L}^2+{\Delta }_{{\mathrm{SINS}}_{\lambda }}^2$${\Delta }_{{\mathrm{SINS}}_L}=\frac{1}{R}\underset{T_{COH}}{\int }{\Delta }_{v_y}^{SINS}$                  (26)${\Delta }_{{\mathrm{SINS}}_{\lambda }}=\frac{\sec L}{R}\underset{T_{COH}}{\int }{\Delta }_{v_x}^{SINS}$式中，和分别为导航坐标系x,y轴方向的惯导系统速度误差，由式(18)计算获得；根据式(23)得到码相位估算值，作为相位搜索空间中心；并根据式(25)计算相位误差的最大值，作为搜索边界，从而压缩搜索空间，提高高动态环境中的捕获效率；步骤三：对多普勒参数及码相位估算值进行卡尔曼滤波平滑，并利用平滑结果计算本地信号与接收信号的动态频率/相位偏移量；为了精确补偿不同相干累积周期间本地信号与接收信号间的动态偏移量，需要保证累积矩阵循环平移的精度，因此预先平滑处理信号捕获参数，选取卡尔曼滤波器的状态量和量测量分别为$X={\left(\begin{array}{ccc}\tau & f_d& {\stackrel{·}{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{·}{\tau }\\ {\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}\tau \\ f_d\\ {\stackrel{·}{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\pi /{\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}{\sigma }_{\tau }^2& {\sigma }_{Dopp}^2& {\sigma }_{Dopp\_Rate}^2\end{array}\right)---\left(32\right)$式中，w_b为时钟偏差；w_d为时钟漂移；w_a为视线加速度的驱动噪声；利用卡尔曼滤波器实时输出的信号参量平滑值，计算当前相干累积周期的累积矩阵输出与下一周期输出在载波频率搜索空间内的动态偏移量以及在码相位搜索空间内的动态偏移量其中$\delta \hat{f}=(f_d-{\tilde{f}}_d)·\Delta t+{\stackrel{·}{f}}_d·{\mathrm{\Delta t}}^2/2$$\delta \hat{\tau }=(\tau -\tilde{\tau })·\Delta t+\frac{f_{CA}}{f_{L1}}·{\stackrel{·}{f}}_d·{\mathrm{\Delta t}}^2/2---\left(33\right)$式中，分别为当前搜索的多普勒频移值和码相位值；Δt为相干累积周期；步骤四：采用循环平移算法补偿相干累积输出矩阵的动态偏移量；对于第i个相干累积周期，由多普勒频移率导致的累积输出矩阵动态频率和相位偏移量分别为和在进行非相干累积之前，为了消除前i个周期的累积动态偏移量，第i个累积输出矩阵需沿频率搜索搜索空间循环平移个搜索单元，沿相位搜索搜索空间循环平移个搜索单元，若捕获搜索的频率分辨率和相位分辨率分别为f_Bin和τ_Bin，则和表示为$X_f^i=\left[\frac{\underset{n=1}{\overset{i}{\Sigma }}\delta {\hat{f}}_i}{f_{Bin}}\right],X_{\tau }^i=\left[\frac{\underset{n=1}{\overset{i}{\Sigma }}\delta {\hat{\tau }}_i}{{\tau }_{Bin}}\right]---\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 }}\delta {\hat{\tau }}_k\right)\exp \left[{\mathrm{\pi NT}}_s\left({\tau }_{Bin}{X}_{\tau }^i-\underset{k=i}{\overset{n-1}{\Sigma }}\delta {\hat{f}}_k\right)\right]{|}^2\\ =A^2{N}^2\underset{k=1}{\overset{n}{\Sigma }}|\begin{array}{c}R\left({\tau }_{Bin}·{\left[\underset{k=i}{\overset{n-1}{\Sigma }}\delta {\hat{\tau }}_k/{\tau }_{Bin}\right]}_{round}-\underset{k=i}{\overset{n-1}{\Sigma }}\delta {\hat{\tau }}_k\right)·\\ \exp \left[{\mathrm{\pi NT}}_s\left(f_{Bin}·{\left[\underset{k=i}{\overset{n-1}{\Sigma }}\delta {\hat{f}}_k/f_{Bin}\right]}_{round}-\underset{k=i}{\overset{n-1}{\Sigma }}\delta {\hat{f}}_k\right)\right]\end{array}{|}^2\end{array}---\left(35\right)$${\omega }_{Bin}·{\left[\underset{k=i}{\overset{n-1}{\Sigma }}\delta {\hat{\omega }}_k/{\omega }_{Bin}\right]}_{round}-\underset{k=i}{\overset{n-1}{\Sigma }}\delta {\hat{\omega }}_k<<{\omega }_{Bin}\omega =\tau ,f$若CP_nch＞2Threshold，则认为该颗卫星已成功捕获，其中，Threshold取非相干累积输出的第二峰值。