GPS信号的C/A码是码率为1.023 Mbps的Gold码，码长N=1 023。不同卫星的C/A码相互正交，同一卫星的C/A码及其循环移位序列也相互正交。因此，C/A码及其循环移位序列可以构成一个正交矩阵。将该正交矩阵作为变换基，就可以对GPS信号进行稀疏表示。

假设g₀=[g(0), g(1), ···, g(N−1)]^T∈R^N×1为任意满足上述条件的Gold码，将g₀循环移位m个码片后，得到序列g_m=[g(m), g(m+1), ···, g(N−1), g(0), g(1), ···, g(m−1)]^T∈R^N×1，m=0, 1, 2, ···, N−1, 则g_m可以表示为：

式中，G=[g₀, g₁, ···, g_N−1]∈R^N×N是由Gold码及其循环移位序列构成的正交变换基矩阵；γ_m=[γ (0), γ (1), ···, γ (N−1)]^T∈R^N×1是变换系数向量。显然，在γ_m中，除了γ (m)=1外，其余系数都为0，即γ_m是稀疏向量。因此，GPS中任意相位的Gold码序列都可以用式(1)进行稀疏表示。

据此可对捕获过程中的GPS信号进行稀疏表示。捕获的实质是进行多普勒频率和C/A码相位二维搜索，根据相关峰位置获得频率和相位估计值。在对某颗GPS卫星的采样信号进行载波(多普勒)剥离和半码片数据打包等处理后，一个C/A码周期的信号可以表示为：

式中，由于进行半码片打包，一个C/A码周期的序列长度为2N=2 046；ω_d是GPS信号的多普勒频率；${\hat \omega _d}$是多普勒频率估计值；Δω_d=${\hat \omega _d}$−ω_d是多普勒估计误差；a、d和ϕ分别表示GPS信号的幅度、导航电文和载波初相，考虑到导航电文码率较低、信号幅度变化缓慢，这里假设它们在捕获处理的时间段内恒定，表示为β=adexp(jϕ)；v(n)是噪声项；c(n; ω_d)表示多普勒为ω_d时的C/A码序列的半码片打包数据。

当多普勒估计准确时，有${\hat \omega _d}={\omega _d}$和Δω_d=0，式(2)表示为：

式中的近似处理考虑如下：在非高动态应用中，多普勒频率的范围通常为[−5, 5] KHz，折算到C/A码上约为[−3.24, 3.24] Hz。相对于1.023 M的C/A码码率而言，这个多普勒频率在1 ms内的影响可以忽略不计^[2]。因此，式中有c(n; ω_d)≈c(n; 0)。

定义c₀=[c(0; 0), c(1; 0), ···, c(2N−1; 0)]^T∈C^2N×1。根据上面的分析，c₀及其循环移位序列之间仍可以看作是近似正交的。令c_m为c₀循环移位m个数据后的序列，即c_m=[c(m; 0), c(m+1; 0), ···, c(2N−1; 0), c(0; 0), c(1; 0), ···, c(m−1; 0)]^T∈R^2N×1，m=0, 1, 2, ···, 2N−1。因此，在正确剥离载波和多普勒后，从任意码相位开始的一个C/A码周期的信号都可以构成一个向量，即r(ω_d)=[r(0; ω_d), r(1; ω_d), ···, r(2N−1; ω_d)]^T∈C^2N×1，它可以写成如下形式：

式中，C=[c₀, c₁, ···, c_2N−1]∈C^2N×2N为变换基矩阵；p(ω_d)=β[p(0; ω_d), p(1; ω_d), ···, p(2N−1; ω_d)]^T∈C^2N×1为多普勒估计准确时的码相位向量；v=[v(0), v(1), ···, v(2N−1)]^T∈C^2N×1为噪声向量。

显然，如果多普勒频率被成功剥离，则p(ω_d)是一个稀疏向量，仅有少数较大的非零值，其峰值出现的地方为C/A码相位估计值。与之相对，如果多普勒估计值存在误差，即${\hat \omega _d} \ne {\omega _d}$，则${{r}}({\hat \omega _d})$中仍残留有多普勒分量，${{p}}({\hat \omega _d})$的稀疏性就会受到影响。

基于式(1)中GPS信号的稀疏表示，可利用压缩感知理论^[16-17]完成对稀疏向量p(ω_d)的估计，从而取代传统捕获算法中的循环相关运算。

使用和变换基矩阵C不相关的观测矩阵Ω∈C^M×2N，M<<2N，对r(ω_d)进行压缩采样，得到观测序列：

式中，${{\varTheta}}={{\varOmega C}}$∈C^M×2N；u=Ωv。实际应用中，观测矩阵Ω的选择很多，常用的有随机高斯矩阵、贝努利矩阵、随机傅立叶矩阵等。

由于多普勒估计准确时，p(ω_d)为稀疏向量，可以通过求解如下优化问题：

获得p(ω_d)的估计值。

基于上述模型，本文提出了基于压缩感知的GPS信号捕获算法，其核心步骤为：

1)初始化捕获参数，获得随机观测矩阵Ω、基矩阵C和混合矩阵Θ=ΩC；

2) Repeat；

3)设置多普勒估计值${\hat \omega _d}$；

4)剥离载波和多普勒频率，半码片打包，获得${{r}}({\hat \omega _d})$；

5)压缩采样，获得${{y}}({\hat \omega _d})={{\varOmega r}}({\hat \omega _d})$；

6)$ \begin{array}{l} {{p}}({{\hat \omega }_d})=\mathop {\arg \min }\limits_{{{p}}({{\hat \omega }_d})} \;||{{p}}({{\hat \omega }_d})|{|_0} \\ \quad \quad \quad \quad {\rm{s}}{\rm{.t}}{\rm{.}}\;\quad \;{{y}}({{\hat \omega }_d})={{\varTheta p}}({{\hat \omega }_d}) \\ \end{array}$

7) Until遍历所有多普勒估计值${\hat \omega _d}$；

8)分析所有的${\bf{p}}({\hat \omega _d})$，判断是否捕获成功；

该算法与传统捕获算法最大的差别是不再需要循环相关运算，而是通过求解(P1)获得当前多普勒频率估计值${\hat \omega _d}$对应的相位估计结果${{p}}({\hat \omega _d})$，如步骤6)所示。与此同时，捕获过程仅需要存储维度为M的向量${{y}}({\hat \omega _d})$，而非维度为2N的向量${{r}}({\hat \omega _d})$，这里有M<<2N。当获得所有多普勒估计值对应的相位估计结果后，通过峰值检测和计算峰均比来判断是否成功捕获到某颗卫星的信号，并给出相应的多普勒频率和C/A码相位估计值。

由于l₀范数项是非凸的，(P1)的求解十分复杂，常见处理方法是用l₁范数来近似l₀范数，即：

由此得到的(P2)为凸优化问题，常见求解算法有基追踪法、(正交)匹配追踪法(orthogonal matching pursuit, OMP)、子空间追踪法^[16-17]等。除此之外，还可以通过软件包(如CVX^[24]等)直接求解。但是，上述一系列追踪算法本质上是贪婪搜索算法，计算复杂度偏高^[17]，同时难以并行执行(基于软件包的方法同样难以并行执行)。考虑到整个优化问题维数较高，分布式的并行求解算法显然更具吸引力。为此，本文将(P2)纳入ADMM的框架，提出一种高效的并行GPS信号捕获算法。

在下面的讨论中，考虑(P2)更一般的形式：

由于存在噪声v[n]，(P2a)可能没有可行解。为了完成捕获，对(P2a)做如下变形，得到：

式中，λ>0用于调节(P2a)求解过程中目标函数和限制条件的权重。

为将(P2b)分解成相对独立的子问题，引入冗余变量${{q}}({\hat \omega _d})={{p}}({\hat \omega _d})$，则(P2b)可以等效为：

使用增广拉格朗日(augmented Lagrangian)方法^[25]，(P2c)的最优解可以通过求解如下问题获得：

式中，ρ>0为惩罚因子；η=[η₀, η₁, ···, η_2N−1]^T是${{q}}({\hat \omega _d})={{p}}({\hat \omega _d})$的拉格朗日乘子(Lagrangian Multiplier)；${L_\rho }\left( {{{p}}({{\hat \omega }_d}),\;{{q}}({{\hat \omega }_d}),\;{{\eta}} } \right)$为增广拉格朗日函数，定义为：

于是，可将变量$\left\{ {{{p}}({{\hat \omega }_d}),\;{{q}}({{\hat \omega }_d})} \right\}$分成${{p}}({\hat \omega _d})$和${{q}}({\hat \omega _d})$，按照如下的ADMM框架^[25]迭代求解：

式中，t为迭代序号。尽管从表面上看，增加冗余变量${{q}}({\hat \omega _d})$使问题更加复杂，但更新${{p}}({\hat \omega _d})$和${{q}}({\hat \omega _d})$的问题可以分解成更多相对独立的子问题，从而进行分布式并行求解。

更新${{p}}({\hat \omega _d})$时，利用一阶最优条件，即：

求得${{{p}}^{(t + 1)}}({\hat \omega _d})$的最优解为：

注意${{{\varTheta }}^{\rm H}}{{y}}({\hat \omega _d})$和逆矩阵${({{{\varTheta }}^{\rm H}}{{\varTheta }} + \rho {{I}})^{ - 1}}$都可以事先计算，不妨将后者简记为：

即φ_m∈C^2N×1表示Φ的第m行，m=0, 1, ···, 2N−1。则${{{p}}^{(t + 1)}}({\hat \omega _d})$中的各元素可以并行计算，即：

根据式中各参数的维度，可知更新${p^{(t + 1)}}(m;\;{\hat \omega _d})$的运算量为O(2N)；更新整个${{{p}}^{(t + 1)}}({\hat \omega _d})$的运算量为O((2N)²)。

更新${{q}}({\hat \omega _d})$的问题可以表示为：

显然，(P4)可以分解为2N个独立的子问题，分别更新${q^{(t + 1)}}(m;\;{\hat \omega _d})$, m=0, 1, ···, 2N−1，即：

同样，利用一阶最优条件，有：

由于${\rm{|}}q(m;\;{\hat \omega _d})|$是非平滑的，这里需要使用它的次梯度(sub-gradient)，表示为^[26]：

将式代入式有：

式中，$\xi=\rho {p^{(t + 1)}}(m;\;{\hat \omega _d}) - \eta _m^{(t)}$。显然，更新${{{q}}^{(t + 1)}}(m;\;{\hat \omega _d})$的运算量为O(1)；更新整个${{{q}}^{(t + 1)}}({\hat \omega _d})$的运算量为O(2N)。

更新拉格朗日乘子η同样可以并行执行，即：

式中，更新$\eta _m^{(t + 1)}$的运算量为O(1)，更新整个${{{\eta}} ^{(t + 1)}}$的运算量为O(2N)。

因此，上述步骤6)可以利用基于ADMM的并行迭代算法来完成，最终获得C/A码相位估计值${{p}}({\hat \omega _d})$，算法的主要步骤如下所示。

1) 初始化ρ, λ, η⁽⁰⁾和${{{q}}^{(0)}}({\hat \omega _d})$, 计算${{\varPhi }}={({{{\varTheta }}^{\rm H}}{{\varTheta }} + \rho {{I}})^{ - 1}}$;

2) 重复

3) 更新${{{p}}^{(t + 1)}}({\hat \omega _d})$

for m=0, 1, 2, …, 2N−1, 并行执行

${p^{(t + 1)}}(m;\;{\hat \omega _d})={{\bf{\varphi }}_m}\left[ {{{\bf{\varTheta }}^{\rm H}}{{y}}({{\hat \omega }_d}) + \rho {{{q}}^{(t)}}({{\hat \omega }_d}) + {{\bf{\eta }}^{(t)}}} \right]$

end

4) 更新${{{q}}^{(t + 1)}}({\hat \omega _d})$

for m=0, 1, 2, ···, 2N−1, 并行执行

${q^{(t + 1)}}(m;\;{\hat \omega _d})=\left\{ \begin{array}{l} 0\quad \quad \quad \left| \xi \right| \leqslant \lambda \\ \dfrac{{\left| \xi \right| - \lambda }}{\rho } \dfrac{\xi }{{\left| \xi \right|}}\quad \text{其他} \\ \end{array} \right.$

end

5) 更新${{\bf{\eta }}^{(t + 1)}}$:

for m=1, 2, ···, 2N−1, 并行执行

$\eta _m^{(t + 1)}=\eta _m^{(t)} + \rho \left[ {{q^{(t + 1)}}(m;\;{{\hat \omega }_d}) - {p^{(t + 1)}}(m;\;{{\hat \omega }_d})} \right]$

end

6) Until迭代过程收敛；

ADMM算法可以确保求得(P2b)的全局最优解。由文献[25]可知，如果(P2b)具有可行解，且ADMM各个子问题都可以唯一地求解，则ADMM迭代过程的每一个聚集点都是原问题的最优解。由于(P2b)是一个无限制条件的凸优化问题，因而必然具有可行解。同时，由于ADMM迭代过程的每一个子问题均为强凸(strongly convex)的，因此，使用一阶最优条件可以唯一地获得其最优解。因此，本文提出的ADMM算法能够获得原问题的全局最优解。

同时，本文提出的ADMM算法具有可分解的结构，特别适合分布式并行实现，能够提高算法效率。如图1所示，整个捕获算法可以由2N个捕获通道并行执行，同时需要一个数据通道负责数据收集和分发。具体而言，每一轮ADMM迭代分为4个阶段。在阶段1，通道m更新${p^{(t + 1)}}(m;\;{\hat \omega _d})$；在阶段2，通道m更新${q^{(t + 1)}}(m;\;{\hat \omega _d})$；在阶段3，通道m更新$\eta _m^{(t + 1)}$；在阶段4，通道m将${q^{(t + 1)}}(m;\;{\hat \omega _d})$和$\eta _m^{(t + 1)}$传送给数据通道，数据通道收集整理这些数据，获得${{{q}}^{(t + 1)}}({\hat \omega _d})$和${{{\eta }}^{(t + 1)}}$，并将它们广播给2N个捕获通道。然后转入阶段1，开始下一轮ADMM迭代。

图  1  ADMM捕获算法的分布式并行执行

注意到在算法并行迭代的每一步都有简单闭合解，因而运算量很低。根据上述分析，ADMM算法每轮迭代的运算量为O((2N)²)，并行执行过程中单个通道的运算量为O(2N)。OMP算法每轮迭代的平均运算量为O(N³)，主要用于求解最小二乘问题，同时，OMP算法难以并行执行。因此，本文提出的ADMM算法比OMP算法更高效。

简便起见，使用仅包含1颗GPS卫星信号的数据进行仿真测试。该卫星为2号星，半码片的C/A码相位为201，多普勒频率为1 KHz，信号强度为−125 dBm。

图2为ADMM算法的典型收敛曲线，此时压缩比为M/(2N)=0.7，即一次捕获利用了2 046×0.7≈1 432个数据包，等效的采样频率为1.432 MHz。图中的结果表明，在各种惩罚因子取值下，ADMM算法的迭代过程都能快速收敛。一般经过50轮迭代后，ADMM捕获算法收敛。

图  2  ADMM捕获算法的典型收敛曲线

图3所示为ADMM算法的典型捕获结果。在多普勒频率和C/A相位的二维平面上，出现了明显的相关峰，表明捕获到了2号星的信号。同时，相关峰出现的位置对应的多普勒频率为1 000 Hz，C/A码相位为201，与预设值相同，实现正确捕获。

图  3  ADMM捕获算法的典型捕获结果

图4为各种捕获算法在不同信号强度下的捕获概率对比。参与比较的算法有：1) 传统时域循环相关算法；2) 基于OMP的捕获算法，它利用OMP算法求解压缩感知捕获问题(P1)，详见文献[19-20]；3) 基于SFT的频谱捕获算法，详见文献[23]；4) 本文的基于ADMM的捕获算法。图中结果表明，基于压缩感知的捕获算法能够在低于Nyquist采样频率时完成GPS信号的捕获，代价是捕获概率略有下降。与之相对，此时传统的时域循环相关算法几乎不能进行捕获。一般而言，采样频率越高(使用的数据量越多，压缩率越高)，捕获概率越高。在3种压缩感知捕获算法中，基于ADMM的捕获算法和基于OMP的捕获算法性能接近，优于基于SFT的捕获算法。与基于OMP的捕获算法相比，本文提出的基于ADMM的捕获算法运算量更低，且算法具有可分解结构，十分利于并行实现，使得算法的效率进一步提升，因而实用性更强。

图  4  各种算法的捕获概率对比

本文提出了一种基于ADMM的并行GPS信号捕获算法。基于GPS信号在相关域的稀疏特性，该算法利用C/A序列构造正交基，对GPS信号进行稀疏表示，据此完成了基于压缩感知的GPS信号捕获问题建模。为了高效地求解该问题，本文将该其纳入ADMM算法框架，提出一种高效的并行捕获算法。在该算法中，原压缩感知问题被分解成多个相对独立的子问题并行迭代求解，并且迭代的每一步都有简单闭合解，因而具有很低的运算量。仿真结果验证了该算法的正确性和有效性。