0 概述

在非合作直接序列扩频(Direct Sequence Spread Spectrum, DSSS)信号通信系统中, 接收者无法截获信号的先验知识, 因此, DSSS信号伪随机码的估计对信息的获取是一个关键的、重要的技术。对于短码DSSS信号的PN码估计, 已经有较多的研究。长码DSSS信号一个扩频周期需调制多位信息符号, 截获长码DSSS信号的PN码的盲估计变得更加困难和具有挑战性。

2000年, 文献[1]提出利用特征值分解(EVD)方法估计DSS信号的扩频码, 该方法已广泛地应用于工程上, 并有许多改进的EVD方法^[1-3], 其中的一些方法已经应用于长码DSSS信号的PN码估计, 其扩频码长为m序列^[4-6]。文献[7]利用爬山算法实现对同步DSSS信号的扩频码估计, 通过将扩频波形估计问题视为一个低秩缺失数据矩阵近似问题。文献[8]提出一个适用于长码DSSS信号的扩频波形估计方法。文献[9]提出2种基于自相关和主分量分析的方法来估计DSSS信号的序列长度。三阶统计数据对缺少样本不敏感, 可以公开更多的信号信息, 因此广泛应用于信号检测和参数估计。文献[10]将三阶相关函数(TCF)应用于m序列估计。文献[11]提出三阶统计过程可以用于非合作环境下DSSS信号检测甚至确定码特征。文献[12]提出对于周期长码信号, 可以通过将其等效为周期时变信道来解决。文献[12]提出一种新的方法解决信号方差矩阵。该方法首先对接收信号进行码片延迟相乘, 然后分析相关矩阵, 得到用于长伪随机码DSSS系统的伪随机码转换序列的估计序列。文献[13]提出对于周期长码直扩信号, 可以等效为同步多用户短码直扩信号来解决。文献[14]提出将周期长码直扩信号等效为多用户短码直扩CDMA系统, 该方法减弱了对信号类型的限制, 并且可以获得非同步信号的扩频序列估计。文献[15]引入了一种由主分量分析实现自适应特征提取的在线无监督学习神经网络的方法, 来实现直接序列码分多址(DS-CDMA)信号的伪码序列盲估计。文献[16]针对批处理方法在实现非等功率同步DS-CDMA信号伪码序列盲估计时, 存在复杂度高、收敛速度慢的问题, 引入了Sanger NN、LEAP NN和APEX NN 3种多主分量神经网络, 实验结果表明, 该方法具有较低的复杂度。但是, 上述方法同样存在不足, 即所处理的信号模型都是线性的, 当信号模型为非线性时, 上述方法均失效。

无迹卡尔曼滤波(UKF)方法在处理非线性问题时具有明显的优势, 因此, 本文针对长码直扩信号的扩频码及信息序列盲估计问题, 结合重叠分段和MCMC算法, 提出一种基于重叠分段MCMC-UKF的盲估计算法。

1 信号模型

持续时间为N个扩频周期的接收信号为:

$ \begin{array}{l} r\left( t \right) = s\left( t \right) \times g\left( t \right) + w\left( t \right) = \sum\limits_{i = 0}^{{N_0} - 1} {b\left( i \right)q\left( {t - i{T_s} - \tau } \right)} \times \\ \;\;\;\;\;\;\;\;\;\sum\limits_{m = 0}^{N - 1} {h\left( {t - mT - \tau } \right){{\rm{e}}^{{\rm{j}}2{\rm{ \mathsf{ π} }}ft - \varphi }}} + w\left( t \right) \end{array} $

(1)

其中, s(t)表示持续时间为N个扩频周期的发射信号模型:

$ s\left( t \right) = A\sum\limits_{i = 0}^{{N_0} - 1} {b\left( i \right)q\left( {t - i{T_s}} \right)} \sum\limits_{m = 0}^{N - 1} {h\left( {t - mT} \right){{\rm{e}}^{{\rm{j}}2{\rm{ \mathsf{ π} }}ft}}} $

(2)

其中, f为载波频率, N₀为N个扩频周期长度包含的信息符号的数目。信号经过信道传播, 其表达式为:

$ g\left( t \right) = \delta \left( {t - \tau } \right){{\rm{e}}^{{\rm{j}}2{\rm{ \mathsf{ π} }}ft}} $

(3)

其中, δ(·)为狄拉克函数, τ为时延。

初相φ=2πfτ通过载波同步^[16]估计出信号载频和初相, 对接收信号r(t)进行下变频, 得到等效基带信号为:

$ \begin{array}{l} r\left( t \right) = \sum\limits_{i = 0}^{{r_0} - 1} {b\left( i \right)} \times \sum\limits_{l = 1}^L {{{\tilde \beta }_l}q\left( {t - i{T_s} - {\tau _l}} \right)} \times \\ \;\;\;\;\;\;\;\;\;\sum\limits_{m = 0}^{N - 1} {h\left( {t - mT - {\tau _l}} \right)} + \tilde w\left( t \right) \end{array} $

(4)

其中, $\tilde{w}(t)=w(t) \mathrm{e}^{-\mathrm{j} 2 {\rm{ \mathsf{ π} }} f t+\varphi}$为均值为0、方差为σ²的加性复高斯白噪声。以码片速率对接收信号进行采样, 定义$r(k) \triangleq r\left(k T_{c}\right)$, 则可得接收信号的离散形式为:

$ \begin{array}{l} r\left( k \right) = \sum\limits_{i = 0}^{{N_0} - 1} {b\left( i \right)q\left( {k - iP - \frac{\tau }{{{T_c}}}} \right)} \times \\ \;\;\;\;\;\;\;\;\;\;\sum\limits_{m = 0}^{N - 1} {h\left( {k - mR - \frac{\tau }{{{T_c}}}} \right)} + \tilde w\left( k \right) \end{array} $

(5)

将式(5)以扩频周期进行分段, 得到的分段观测数据如图 1所示。

得到的第n个分段信号表示为:

$ \begin{array}{l} r\left( n \right) = \left[ {{b_0}\left( n \right){c_0}\left( n \right) + {b_1}\left( n \right){c_1}\left( n \right) + \cdots + } \right.\\ \;\;\;\;\;\;\;\;\;\;\left. {{b_M}\left( n \right){c_M}\left( n \right)} \right] + w\left( n \right),1 \le n \le N \end{array} $

(6)

其中, M=R/P为扩频比, $r(n) \triangleq[r(n R+1), r(n R+2)$, $\cdots, r(n R+R) ], b_{m}(n) \in \mathbb{C}^{1}, c_{m}(n) \in \mathbb{R}^{1 \times R}, 0 \leqslant m \leqslant M$, 且有:

$ \left\{ \begin{array}{l} {b_0}\left( n \right) = b\left( {Mn} \right)\\ {b_1}\left( n \right) = b\left( {Mn + 1} \right)\\ \;\;\;\;\;\;\; \vdots \\ {b_M}\left( n \right) = b\left( {Mn + M} \right) \end{array} \right. $

(7)

$ \begin{array}{l} {\mathit{\boldsymbol{c}}_0}\left( n \right) = {\left[ {c\left( {\tau /{T_c} + 1} \right)c\left( {\tau /{T_c} + 2} \right) \cdots c\left( P \right)\underbrace {\begin{array}{*{20}{c}} 0& \cdots &0 \end{array}}_{R - P + \tau /{T_c}}} \right]^{\rm{T}}}\\ {\mathit{\boldsymbol{c}}_M}\left( n \right) = {\left[ {\underbrace {\begin{array}{*{20}{c}} 0& \cdots &0 \end{array}}_{R - \tau /{T_c}}c\left( 1 \right)c\left( 2 \right) \cdots c\left( {\tau /{T_c}} \right)} \right]^{\rm{T}}}\\ {\mathit{\boldsymbol{c}}_m}\left( n \right) = {\left[ {\underbrace {\begin{array}{*{20}{c}} 0& \cdots &0 \end{array}}_{P - \tau /{T_c} + \left( {m - 1} \right)P}c\left( 1 \right)c\left( 2 \right) \cdots c\left( P \right)\underbrace {\begin{array}{*{20}{c}} 0& \cdots &0 \end{array}}_{R - mP - P + \tau /{T_c}}} \right]^{\rm{T}}} \end{array} $

(8)

将式(6)进一步写成下面的向量形式:

$ \mathit{\boldsymbol{r}}\left( n \right) = \mathit{\boldsymbol{b}}\left( n \right)\mathit{\boldsymbol{c}}\left( n \right) + \mathit{\boldsymbol{w}}\left( n \right),1 \le n \le N $

(9)

其中, b(n)=[b₀(n)   b₁(n)   …   b_M(n)], c(n)=[c₀(n)  c₁  (n) … c_M(n)]^T。结合图 1并分析式(6)~式(9)可以发现, 在n取不同值时, 各时间窗内第一个符号起始时刻相对于该时间窗起始时刻的偏移不相同, 从而导致c(n)为一个时变的向量, 因此无法通过建立统一的数学模型, 采用MCMC方法对扩频码和信息序列进行联合估计。

针对上述问题, 下文将一个完整的扩频周期分成小的片段, 使扩频码在每个片段的持续时间内对应一个符号, 从而使c(n)保持恒定, 这样就可以通过MCMC方法估计出每个片段的扩频码及信息序列。对于长码直扩信号, 为进一步提高算法的估计性能, 避免每一段估值运算时产生的扩频码极性不一致问题, 将一个扩频周期向量分成S个的相互重叠。图 2所示为重叠长度为H₀、长度为H的信号片段, 其中, 信号片段的长度H远小于扩频增益P。

则式(5)可重新表示成下式的形式:

$ R_{s}(k)=A(\boldsymbol{b}(s)) D(\boldsymbol{c}(s), \tau)+W(k) $

(10)

即:

$ \begin{array}{l} \left[ {\begin{array}{*{20}{c}} {r_1^s\left( k \right)}\\ {r_2^s\left( k \right)}\\ \vdots \\ {r_N^s\left( k \right)} \end{array}} \right]{\rm{ = }}\left[ {\begin{array}{*{20}{c}} {b\left( {i_1^s} \right)}\\ {b\left( {i_2^s} \right)}\\ \vdots \\ {b\left( {i_N^s} \right)} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {h_1^s\left( k \right)}\\ {h_2^s\left( k \right)}\\ \vdots \\ {h_L^s\left( k \right)} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {w_1^s\left( k \right)}\\ {w_2^s\left( k \right)}\\ \vdots \\ {w_N^s\left( k \right)} \end{array}} \right],\\ 1 \le k \le H,1 \le s \le S \end{array} $

(11)

其中:

$ \begin{gathered} r_n^s\left( k \right) \triangleq \left\{ {r\left[ {\left( {n - 1} \right)R + \left( {s - 1} \right)\left( {H - {H_0}} \right) + k} \right]} \right\}_{1 \leqslant k \leqslant H}^{1 \leqslant n \leqslant N} \in \hfill \\ \;\;\;\;\;\;\;\;\;\;\;{\mathbb{C}^{1 \times H}} \hfill \\ \end{gathered} $

$ \mathit{\boldsymbol{b}}\left( s \right) = {\left[ {b\left( {i_1^s} \right),b\left( {i_2^s} \right), \cdots ,b\left( {i_N^s} \right)} \right]^{\text{T}}} $

$ i_n^s = \left\lfloor {\frac{{\left( {n - 1} \right)R + \left( {s - 1} \right)\left( {H - {H_0}} \right)}}{P}} \right\rfloor ,1 \leqslant n \leqslant N $

$ \mathit{\boldsymbol{c}}\left( s \right) = {\left\{ {c\left[ {\left( {s - 1} \right)\left( {H - {H_0}} \right) + k} \right]} \right\}_{1 \leqslant k \leqslant H}} $

$ {h^s}\left( k \right) \triangleq {\left\{ {h\left[ {\left( {s - 1} \right)\left( {H - {H_0}} \right) + k - \frac{\tau }{{{T_c}}}} \right]} \right\}_{1 \leqslant k \leqslant H}} $

本文研究目标是在知道每一小段数据R_s的前提下, 采用MCMC-UKF方法对参数空间Θ $\triangleq${b(s), c(s), τ, σ²}进行联合估计, 重复进行S次独立的运算之后, 进而用所有的样本采样值对扩频码及信息序列进行拼接和恢复。通过分析参数空间Θ中的元素可以发现, 对于参数σ², 可采用UKF方法进行估计, 而b(s)、c(s)和τ的全条件后验分布未知且不易抽取, 因此, 采用混合MH算法进行抽取。

2 算法原理 2.1 无迹卡尔曼滤波

无迹卡尔曼滤波(UKF)是一种递归式贝叶斯估计方法, 但是无须进行非线性模型的求解(即不需要求解雅可比矩阵), 其基本思想是利用UT(Unscented Transformation)转换, 用一组确定的样本点近似求解测量条件下系统状态的后验概率密度p(x_k|z_k)的均值和方差, 实现系统状态递推均值和方差(一阶矩、二阶矩)的估计。

设具有如下状态空间模型的非线性系统状态方程和量测方程分别为:

$ {x_{k + 1}} = f\left( {{x_k}} \right) + {n_k} $

(12)

$ {z_{k + 1}} = h\left( {{x_{k + 1}}} \right) + {v_{k + 1}} $

(13)

其中, n_k是协方差为Q_k的系统噪声, v_k+1是协方差为R_k+1的测量噪声。基于UT变换的UKF算法基本过程如下:

1) 初始化$\hat{x}_{0}=E\left(x_{0}\right), \boldsymbol{P}_{0}=E\left[\left(x_{0}-\hat{x}_{0}\right) \times\right.$$\left(x_{0}-\hat{x}_{0}\right)^{\mathrm{T}} ]$;

2) 确定样本点Sigma集合s_i(k)及相应的权值w_i(k), i=0, 1, …, 2n;

3) 根据系统状态方程求取样本点传递值:

$ {x_{{s_{i\left( {k + 1,k} \right)}}}} = f\left( {{s_{i\left( k \right)}}} \right) $

(14)

4) 系统状态均值和方差的一步预测, 为:

$ {{\hat x}_{k + 1,k}} = \sum\limits_{i = 0}^{2n} {w_i^m{x_{{s_{i\left( {k + 1,k} \right)}}}}} $

(15)

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{P}}_{xx\left( {k + 1,k} \right)}} = {Q_{k + 1}} + \sum\limits_{i = 0}^{2n} {w_i^c\left( {{x_{{s_{i\left( {k + 1,k} \right)}}}} - {{\hat x}_{k + 1,k}}} \right)} \times }\\ {{{\left( {{x_{{s_{i\left( {k + 1,k} \right)}}}} - {{\hat x}_{k + 1,k}}} \right)}^{\rm{T}}}} \end{array} $

(16)

5) 根据系统量测方程进一步求取状态一步预测的传递值:

$ {z_{{s_{i\left( {k + 1,k} \right)}}}} = h\left( {{x_{{s_{i\left( {k + 1,k} \right)}}}}} \right) $

(17)

6) 预测量均值和协方差:

$ {\hat z_{k + 1,k}} = \sum\limits_{i = 0}^{2n} {w_i^mz\left( {{s_{i\left( {k + 1,k} \right)}}} \right)} $

(18)

$ \begin{array}{l} {\mathit{\boldsymbol{P}}_{zz}} = {R_{k + 1}} + \sum\limits_{i = 0}^{2n} {{w_i}\left( {z\left( {{s_{i\left( {k + 1,k} \right)}}} \right) - {{\hat z}_{k + 1,k}}} \right)} \times \\ \;\;\;\;\;\;\;\;{\left( {z\left( {{s_{i\left( {k + 1,k} \right)}}} \right) - {{\hat z}_{k + 1,k}}} \right)^{\rm{T}}} \end{array} $

(19)

$ \begin{array}{l} {\mathit{\boldsymbol{P}}_{xz}} = \sum\limits_{i = 0}^{2n} {w_i^c\left( {x\left( {{s_{i\left( {k + 1,k} \right)}}} \right) - {{\hat x}_{k + 1,k}}} \right)} \times \\ \;\;\;\;\;\;\;\;{\left( {z\left( {{s_{i\left( {k + 1,k} \right)}}} \right) - {{\hat z}_{k + 1,k}}} \right)^{\rm{T}}} \end{array} $

(20)

其中, P_zz是量测方差矩阵, P_xz是状态向量与量测向量的协方差矩阵。

7) 计算UKF增益, 更新状态向量与方差, 即:

$ {\mathit{\boldsymbol{K}}_{k + 1}} = {\mathit{\boldsymbol{P}}_{xz}}\mathit{\boldsymbol{P}}_{_{zz}}^{ - 1} $

(21)

$ {\hat x_{k + 1,k + 1}} = {\hat x_{k + 1,k}} + {\mathit{\boldsymbol{K}}_{k + 1}}\left( {{z_{k + 1}} - {{\hat z}_{k + 1,k}}} \right) $

(22)

$ {\mathit{\boldsymbol{P}}_{xx\left( {k + 1,k + 1} \right)}} = {\mathit{\boldsymbol{P}}_{xx\left( {k + 1,k} \right)}} - {\mathit{\boldsymbol{K}}_{k + 1}}{\mathit{\boldsymbol{P}}_{zz}}\mathit{\boldsymbol{K}}_{k + 1}^{\rm{T}} $

(23)

2.2 MCMC算法

MCMC算法的基本思想是构造一条Markov链, 使其平稳分布为待估参数的后验分布。利用这条链上的各个样本值就可以估计参数, 进行各种统计推断, 该算法包括2个部分:

1) Monte Carlo部分利用平均近似一个期望, 如需计算期望E[f(θ)]=∫f(θ)p(θ|X)dθ, 则一个好的近似为$E[f(\theta)] \approx \frac{1}{I} \sum\limits_{i=0}^{I-1} f\left(\theta^{(i)}\right)$, 当I→∞时一定收敛于真实期望, 即取等号, 这里θ⁽ⁱ⁾是依据p(θ|X)产生的一系列样本;

2) Markov部分产生具有不变分布p(θ|X)(一般用π(θ)表示不变分布)的Markov链{θ⁽⁰⁾, θ⁽¹⁾, …}, 一个Markov链由θ⁽⁰⁾初始分布和所谓的转移核l(θ⁽ⁱ⁾, dθ⁽ⁱ⁺¹⁾)完全定义。经典的MCMC方法有Gibbs采样和MH算法, 它们定义的转移核使得无论θ⁽⁰⁾是何种初始分布, 都一定收敛于代价分布p(θ|X)。

下面给出2种算法的具体实现步骤。

1) Gibbs采样

参数空间θ={θ₁, θ₂, …, θ_d}, 为了从后验分布p(θ|X)中抽取样本值, 利用前一时刻的抽样值, 对各参数的全条件后验分布进行迭代抽取, 得到后一时刻各参数的抽样值。假设在n－1次迭代时刻的空间值为θ^(n－1)={θ₁^(n－1), θ₂^(n－1), …, θ_d^(n－1)}, 于是n时刻的空间值通过以下抽取步骤得到:

步骤1  从p(θ₁|θ₂^(n－1), θ₃^(n－1), …, θ_d^(n－1), X)中抽取θ₁⁽ⁿ⁾。

步骤2  利用更新值θ₁⁽ⁿ⁾, 从p(θ₂|θ₁⁽ⁿ⁾, θ₂^(n－1)…, θ_d^(n－1), X)中抽取θ₂⁽ⁿ⁾。

步骤3  利用更新值θ₁⁽ⁿ⁾、θ₂⁽ⁿ⁾, 从p(θ₃θ₁⁽ⁿ⁾, θ₂⁽ⁿ⁾, θ₄^(n－1), …, θ_d^(n－1), X)中抽取θ₃⁽ⁿ⁾。

步骤4  依此类推, 最后利用θ₁⁽ⁿ⁾, θ₂⁽ⁿ⁾, …, θ_d－1⁽ⁿ⁾, 从p(θ_d|θ₁⁽ⁿ⁾, θ₂⁽ⁿ⁾, …, θ_d－1⁽ⁿ⁾, X)中抽取θ_d⁽ⁿ⁾。

2) MH算法

MH算法是Gibbs采样的扩展, 它主要是针对那些全条件后验分布概率不易抽取的参数, 其具体步骤如下:

步骤1  根据先验分布选择初始值θ⁽⁰⁾。

步骤2  在第n次迭代时刻, 从建议分布q(θ^*|θ^(n－1))中抽取备择采样值θ^*。

步骤3  计算接收概率:

$ r = \frac{{p\left( {{\theta ^ * }\left| X \right.} \right)q\left( {{\theta ^{\left( {n - 1} \right)}}\left| {{\theta ^ * }} \right.} \right)}}{{p\left( {{\theta ^{\left( {n - 1} \right)}}\left| X \right.} \right)q\left( {{\theta ^ * }\left| {{\theta ^{\left( {n - 1} \right)}}} \right.} \right)}} $

(24)

步骤4  以接收概率min(1, r)接收θ^*作为θ⁽ⁿ⁾的更新值, 否则θ⁽ⁿ⁾=θ^(n－1)。

步骤5  重复步骤2~步骤4进行多次迭代, 从而完成对p（θ|X)的抽取。

2.3 贝叶斯模型推导

由第1节的信号模型可知, 时延等参数的状态空间模型和观测方程分别为:

$ {\tau _n} = {\tau _{n - 1}} $

(25)

$ {R_s} = A\left( {\mathit{\boldsymbol{b}}\left( s \right)} \right)D\left( {\mathit{\boldsymbol{c}}\left( s \right),\tau } \right) + W\left( s \right) $

(26)

1) 参数的先验假设

先验分布的选取一般基于2个原则:简化计算与对后验分布的影响最小(即无信息先验分布)。在简化计算上一般采用共轭先验分布, 它能保证参数的后验分布与先验分布相同。

由贝叶斯公式可知, 参数空间的后验分布函数等价为:

$ \begin{array}{l} p\left( {\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}\left| {{\mathit{\boldsymbol{R}}_s}\left( k \right)} \right.} \right) = p\left( {\mathit{\boldsymbol{ \boldsymbol{\varTheta} }},{\mathit{\boldsymbol{R}}_s}\left( k \right)} \right)/p\left( R \right) \propto \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;p\left( {{R_s}\left( k \right)/\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}} \right)p\left( \mathit{\boldsymbol{ \boldsymbol{\varTheta} }} \right) \end{array} $

(27)

其中, p(R_s(k)|Θ)为似然函数, p(Θ)为联合先验概率密度函数。根据式(10)可得:

$ \begin{array}{l} p\left( {{\mathit{\boldsymbol{R}}_s}\left( k \right)\left| {\mathit{\boldsymbol {\varTheta} } }\right.} \right) = {\left( {2{\rm{ \mathsf{ π} }}{\sigma ^2}} \right)^{ - NH/2}} \times \\ \exp \left( { - \frac{1}{{2{\sigma ^2}}}\sum\limits_{i = 1}^N {{{\left\| {r_{i,:}^s - A{{\left( {\mathit{\boldsymbol{b}}\left( s \right)} \right)}_{i,:}}D\left( {\mathit{\boldsymbol{c}}\left( s \right),\tau } \right)} \right\|}^2}} } \right) \end{array} $

(28)

其中, (·)_{i, :}表示相应矩阵的第i行, R_{i, :}^s(k)为矩阵R_s(k)的第i行。

$ \begin{array}{l} p\left( \boldsymbol{\varTheta } \right) = p\left( {\mathit{\boldsymbol{b}}\left( s \right),\mathit{\boldsymbol{c}}\left( s \right),\tau ,{\sigma ^2}} \right) = \\ \;\;\;\;\;\;\;\;\;\;\;p\left( {\mathit{\boldsymbol{b}}\left( s \right)\left| {\mathit{\boldsymbol{c}}\left( s \right)} \right.,\tau ,{\sigma ^2}} \right)p\left( {\mathit{\boldsymbol{c}}\left( s \right)\left| {\tau ,{\sigma ^2}} \right.} \right)p\left( {\tau \left| {{\sigma ^2}} \right.} \right)p\left( {{\sigma ^2}} \right) = \\ \;\;\;\;\;\;\;\;\;\;\;p\left( {\mathit{\boldsymbol{b}}\left( s \right)} \right)p\left( {\mathit{\boldsymbol{c}}\left( s \right)} \right)p\left( \tau \right)p\left( {{\sigma ^2}} \right) \end{array} $

(29)

由于噪声方差σ²的先验分布可认为服从逆伽玛分布Ig(v₀, γ₀), 因此概率密度函数表达式为:

$ p\left( {{\sigma ^2}} \right) = \frac{{{{\left( {{\gamma _0}} \right)}^{{\nu _0}}}}}{{\mathit{\Gamma }\left( {{v_0}} \right)}}{\left( {{\sigma ^2}} \right)^{ - {\nu _0} - 1}}\exp \left( { - \frac{{{\gamma _0}}}{{{\sigma ^2}}}} \right) $

(30)

本文选取ν₀=γ₀=0, 逆伽玛分布则成为Jeffrey无信息先验分布p(σ²)∝1/σ²。

扩频码向量c(s)的每个元素相互独立, 可认为服从$\mathbb{S} \triangleq\{-1, +1\}$上的离散均匀分布$D U(\mathbb{S})$, 则其概率密度函数为:

$ p\left( {\mathit{\boldsymbol{c}}\left( s \right)} \right) = {\left( {1/2} \right)^R} $

(31)

信息序列向量b(s)中的元素相互独立, 可认为服从离散有限字符集$\mathbb{C}$上的均匀分布, 若采用BPSK调制, 则$\mathbb{C} \triangleq\{-1, +1\}$。本文中假设调制样式已知, 则其概率密度函数为:

$ p(\mathit{\boldsymbol{b}}(s)) = {\left( {\frac{1}{2}} \right)^N} $

(32)

时延τ服从[0, T_s]内的均匀分布U[0, T_s], 则其概率密度函数为:

$ p(\tau)=\frac{1}{T_{s}} $

(33)

根据上述先验分布, 式(29)可写成:

$ \begin{array}{l} p(\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}) = p(\mathit{\boldsymbol{b}}(s))p(\mathit{\boldsymbol{c}}(s))p(\tau )p\left( {{\sigma ^2}} \right) \propto \\ \;\;\;\left( {\frac{1}{{{T_s}}}} \right){\left( {\frac{1}{2}} \right)^N}{\left( {\frac{1}{2}} \right)^R} \times \frac{{{{\left( {{\gamma _0}} \right)}^{{v_0}}}}}{{\mathit{\Gamma }\left( {{v_0}} \right)}}{\left( {{\sigma ^2}} \right)^{ - {\nu _0} - 1}}\exp \left( { - \frac{{{\gamma _0}}}{{{\sigma ^2}}}} \right) \end{array} $

(34)

由此, 可得:

$ \begin{array}{l} p\left( {\boldsymbol{\varTheta }\left| {{R_s}} \right.} \right) \propto p\left( {{R_s}\left| \boldsymbol{\varTheta } \right.} \right)p\left( \boldsymbol{\varTheta } \right) = {\left( {2{\rm{ \mathsf{ π} }}{\sigma ^2}} \right)^{ - NH/2}} \times \\ \;\;\;\exp \left( { - \frac{1}{{2{\sigma ^2}}}\sum\limits_{i = 1}^N {{{\left\| {R_{i,:}^s - A{{\left( {\mathit{\boldsymbol{b}}\left( s \right)} \right)}_{i,:}}D\left( {\mathit{\boldsymbol{c}}\left( s \right),\tau } \right)} \right\|}^2}} } \right) \times \\ \;\;\;\left( {\frac{1}{{{T_s}}}} \right){\left( {\frac{1}{2}} \right)^N}{\left( {\frac{1}{2}} \right)^R} \times \frac{{{{\left( {{\gamma _0}} \right)}^{{\nu _0}}}}}{{\mathit{\Gamma }\left( {{v_0}} \right)}}{\left( {{\sigma ^2}} \right)^{ - {\nu _0} - 1}}\exp \left( { - \frac{{{\gamma _0}}}{{{\sigma ^2}}}} \right) \end{array} $

(35)

结合式(28)并对σ²积分, 可得:

$ \begin{gathered} p\left( {\mathit{\boldsymbol{b}}\left( s \right),\mathit{\boldsymbol{c}}\left( s \right),\tau \left| {{R_s}} \right.} \right) \propto \hfill \\ {\left( {\frac{{2{\gamma _0} + \sum\limits_{i = 1}^N {{{\left\| {R_{i,:}^s - A{{\left( {\mathit{\boldsymbol{b}}\left( s \right)} \right)}_{i,:}}D\left( {\mathit{\boldsymbol{c}}\left( s \right),\tau } \right)} \right\|}^2}} }}{2}} \right)^{\left( { - \frac{{NH + 2{\nu _0}}}{2}} \right)}} \times \hfill \\ DU\left( \mathbb{C} \right)DU\left( \mathbb{S} \right)U\left[ {0,{T_s}} \right] \hfill \\ \end{gathered} $

(36)

下文给出各参数的全条件后验分布。

1) σ²的全条件后验分布

根据式(26)可知, 在第n个观测周期, σ²服从逆伽玛分布, 为了方便公式推导, 令β=σ², 即:

$ \begin{array}{l} p\left( {\hat \beta _{n + 1}^s\left| {\mathit{\boldsymbol{b}}\left( s \right),\mathit{\boldsymbol{c}}\left( s \right),{\tau _{n + 1}},\Delta {f_{n + 1}},\sigma _{n + 1}^2,R_{n + 1}^s} \right.} \right)\sim \\ \;\;\;\;{N_c}\left( {\hat \beta _{n + 1}^s,{P_{\beta \beta \left( {n + 1} \right)}}} \right) \end{array} $

(37)

其中, 均值$\hat{\beta}_{n+1}^{s}$和方差P_ββ(n+1)通过无迹卡尔曼滤波方法, 在给出$\hat{\beta}_{n+1}^{s}$和P_ββ(n)时进行迭代更新, 其过程如下:

初始化:

$ \hat \beta _n^s = E\left( {\beta _n^s} \right) $

(38)

$ {\mathit{\boldsymbol{P}}_{\beta \beta \left( n \right)}} = E\left[ {\left( {\beta _n^s - \hat \beta _n^s} \right){{\left( {\beta _n^s - \hat \beta _n^s} \right)}^{\rm{T}}}} \right] $

(39)

根据UKF算法, 确定第n个周期的样本点Sigma集合及相应的权值:

$ \left\{ \begin{array}{l} s_{0\left( n \right)}^s = \hat \beta _n^s\\ s_{1\left( n \right)}^s = \hat \beta _n^s + {\left( {\sqrt {\left( {L + \lambda } \right){\mathit{\boldsymbol{P}}_{\beta \beta }}} } \right)_i},i = 1,2, \cdots ,L\\ s_{i\left( n \right)}^s = \hat \beta _n^s - {\left( {\sqrt {\left( {L + \lambda } \right){\mathit{\boldsymbol{P}}_{\beta \beta }}} } \right)_{i - L}},i = L + 1, \cdots ,2L \end{array} \right. $

(40)

$ \left\{ \begin{array}{l} w_0^m = \lambda /\left( {L + \lambda } \right),w_0^c = \lambda /\left( {L + \lambda } \right) + \left( {1 - {\alpha ^2} + \beta } \right)\\ w_i^m = w_i^c = \lambda /\left[ {2\left( {L + \lambda } \right)} \right],i = 1,2, \cdots ,2L \end{array} \right. $

(41)

$ \beta _{{s_{i\left( {n + 1,n} \right)}}}^s = Fs_{i\left( n \right)}^s $

(42)

$ \hat \beta _{n + 1,n}^s = \sum\limits_{i = 0}^{2L} {w_i^m\beta _{{s_{i\left( {n + 1,n} \right)}}}^s} $

(43)

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{P}}_{\beta \beta \left( {n + 1,n} \right)}} = {{\Sigma }_V} + \sum\limits_{i = 0}^{2L} {w_i^c\left( {\beta _{{s_{i\left( {n + 1,n} \right)}}}^s - \hat \beta _{n + 1,n}^s} \right)} \times }\\ {{{\left( {\beta _{{s_{i\left( {n + 1,n} \right)}}}^s - \hat \beta _{n + 1,n}^s} \right)}^{\rm{T}}}} \end{array} $

(44)

$ R_{{s_{i\left( {n + 1,n} \right)}}}^s = A\left( {\mathit{\boldsymbol{b}}\left( s \right)} \right)D\left( {\mathit{\boldsymbol{c}}\left( s \right),{\tau _n}} \right) $

(45)

$ \hat R_{n + 1,n}^s = \sum\limits_{i = 0}^{2L} {w_i^mR\left( {{s_{i\left( {n + 1,n} \right)}}} \right)} $

(46)

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{P}}_{RR}} = {{\Sigma }_W} + \sum\limits_{i = 0}^{2L} {{w_i}\left( {R\left( {{s_{i\left( {n + 1,n} \right)}}} \right) - \hat R_{n + 1,n}^s} \right)} \times }\\ {{{\left( {R\left( {{s_{i\left( {n + 1,n} \right)}}} \right) - \hat R_{n + 1,n}^s} \right)}^{\rm{T}}}} \end{array} $

(47)

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{P}}_{\beta R}} = \sum\limits_{i = 0}^{2L} {w_i^c\left( {\beta \left( {{s_{i\left( {n + 1,n} \right)}}} \right) - \hat \beta _{n + 1,n}^s} \right)} \times }\\ {{{\left( {\beta \left( {{s_{i\left( {n + 1,n} \right)}}} \right) - \hat \beta _{n + 1,n}^s} \right)}^{\rm{T}}}} \end{array} $

(48)

$ {K_{n + 1}} = {\mathit{\boldsymbol{P}}_{\beta R}}\mathit{\boldsymbol{P}}_{RR}^{ - 1} $

(49)

$ \hat \beta _{n + 1,n + 1}^s = \hat \beta _{n + 1,n}^s + {K_{n + 1}}\left( {R_{n + 1}^s - \hat R_{n + 1,n}^s} \right) $

(50)

$ {P_{\beta \beta \left( {n + 1,n + 1} \right)}} = {\mathit{\boldsymbol{P}}_{\beta \beta \left( {n + 1,n} \right)}} - {K_{n + 1}}{\mathit{\boldsymbol{P}}_{RR}}\mathit{\boldsymbol{K}}_{n + 1}^{\rm{T}} $

(51)

其中, L为多样本点数, $\hat{\beta}_{n+1, n}^{s}$为状态预测均值, P_{ββ(n+1, n)}为状态预测值的协方差矩阵, $\hat{R}_{n+1, n}^{s}$为预测的观测均值, P_RR是观测方差矩阵, P_βR是状态向量与观测向量的协方差矩阵, K_n+1为卡尔曼增益。由此, 可得$\tilde{\beta}_{n+1}^{s}$的全条件后验分布为:

$ \begin{array}{l} p\left( {\tilde \beta _{n + 1}^s\left| {\mathit{\boldsymbol{b}}\left( s \right),\mathit{\boldsymbol{c}}\left( s \right),\tau ,{R_s}} \right.} \right) \sim \\ Ig\left( {{\nu _0} + \frac{{NH}}{2},\frac{1}{2}\sum\limits_{i = 1}^N {{{\left\| {R_{i,:}^s - A{{\left( {\mathit{\boldsymbol{b}}\left( s \right)} \right)}_{i,:}}D\left( {\mathit{\boldsymbol{c}}\left( s \right),\tau } \right)} \right\|}^2}} + {\gamma _0}} \right) \end{array} $

即σ²的全条件后验分布为:

$ \begin{array}{l} p\left( {{\sigma ^2}\left| {\mathit{\boldsymbol{b}}\left( s \right),\mathit{\boldsymbol{c}}\left( s \right),\tau ,{R_s}} \right.} \right) \sim \\ Ig\left( {{\nu _0} + \frac{{NH}}{2},\frac{1}{2}\sum\limits_{i = 1}^N {{{\left\| {R_{i,:}^2 - A{{\left( {\mathit{\boldsymbol{b}}\left( s \right)} \right)}_{i,:}}D\left( {\mathit{\boldsymbol{c}}\left( s \right),\tau } \right)} \right\|}^2}} + {\gamma _0}} \right) \end{array} $

(52)

2) b(s)、c(s)和τ的全条件后验分布

由式(36)可知, b(s)、c(s)和τ的全条件后验分布表达式均为:

$ {\left( {\frac{{2{\gamma _0} + \sum\limits_{i = 1}^N {{{\left\| {R_{i,:}^s - A{{\left( {\mathit{\boldsymbol{b}}\left( s \right)} \right)}_{i,:}}D\left( {\mathit{\boldsymbol{c}}\left( s \right),\tau } \right)} \right\|}^2}} }}{2}} \right)^{\left( { - \frac{{NH + 2{\nu _0}}}{2}} \right)}} $

(53)

在计算后验分布时内部的已知条件取值不同。

下文给出采用混合MH算法抽样b(s)、c(s)和τ的过程, 式(53)已经给出扩频码向量c(s)、信息序列向量b(s)以及时延τ中元素的全条件后验概率分布, 并将其分别用作所需抽样参数的不变分布。针对扩频码向量c(s), 首先取[0, 1]之间的均匀分布值μ, 如果μ<λ(λ为仿真参数, 文中取λ=0.5), 则建议概率分布$q_{1}\left(\boldsymbol{c}^{*}(s) | \boldsymbol{c}(s)\right) \sim D U(\mathbb{S})$随机更新扩频码片, 其中$\mathbb{S} \triangleq\{-1, +1\}$, 否则以q₂(c^*(s)|c(s))~N(c(s), σ_s²)更新扩频码片; 同理对于信息向量b(s)分别取建议概率分布$q_{1}\left(\boldsymbol{b}^{*}(s) | \boldsymbol{b}(s)\right) \sim D U(\mathbb{C})$和q₂(b^*(s)|b(s))~N(b(s), σ_s²); 时延τ的建议分布取q₁(τ^*|τ)~U[0, T_s]和q₂(τ_s^*|τ_s)~N(τ_s, σ_s²), 且3种建议概率分布的接受概率分别为:

$ A\left( {{\mathit{\boldsymbol{c}}_n}\left( s \right),\mathit{\boldsymbol{c}}_n^ * \left( s \right)} \right) = \min \left\{ {1,{{\left( {\frac{{2{\gamma _0} + \sum\limits_{i = 1}^N {{{\left\| {R_{i,:}^s - A{{\left( {\mathit{\boldsymbol{b}}\left( s \right)} \right)}_{i,:}}D\left( {\mathit{\boldsymbol{c}}\left( s \right),\tau } \right)} \right\|}^2}} }}{{2{\gamma _0} + \sum\limits_{i = 1}^N {{{\left\| {R_{i,:}^s - A{{\left( {\mathit{\boldsymbol{b}}\left( s \right)} \right)}_{i,:}}D\left( {\mathit{\boldsymbol{c}}_n^ * \left( s \right),\tau } \right)} \right\|}^2}} }}} \right)}^{\left( {\frac{{NH + 2{\nu _0}}}{2}} \right)}}} \right\} $

(54)

$ A\left( {b\left( {i_n^s} \right),{b^ * }\left( {i_n^s} \right)} \right) = \min \left\{ {1,{{\left( {\frac{{2{\gamma _0} + {{\left\| {R_{i,:}^s - A{{\left( {b\left( {i_n^s} \right)} \right)}_{i,:}}D\left( {\mathit{\boldsymbol{c}}\left( s \right),\tau } \right)} \right\|}^2} + {\Sigma }}}{{2{\gamma _0} + {{\left\| {R_{i,:}^s - A{{\left( {{b^ * }\left( {i_n^s} \right)} \right)}_{i,:}}D\left( {\mathit{\boldsymbol{c}}\left( s \right),\tau } \right)} \right\|}^2} + {\Sigma }}}} \right)}^{\left( {\frac{{NH + 2{\nu _0}}}{2}} \right)}}} \right\} $

(55)

$ A\left( {\tau ,{\tau ^ * }} \right) = \min \left\{ {1,{{\left( {\frac{{2{\gamma _0} + \sum\limits_{i = 1}^N {{{\left\| {R_{i.:}^s - A{{\left( {\mathit{\boldsymbol{b}}\left( s \right)} \right)}_{i,:}}D\left( {\mathit{\boldsymbol{c}}\left( s \right),\tau } \right)} \right\|}^2}} }}{{2{\gamma _0} + \sum\limits_{i = 1}^N {{{\left\| {R_{i.:}^s - A{{\left( {\mathit{\boldsymbol{b}}\left( s \right)} \right)}_{i,:}}D\left( {\mathit{\boldsymbol{c}}\left( s \right),{\tau ^ * }} \right)} \right\|}^2}} }}} \right)}^{\left( {\frac{{NH + 2{\nu _0}}}{2}} \right)}}} \right\} $

(56)

其中, $\varSigma=\sum\limits_{l=1, l \neq m}^{N}\left\|R_{i, :}^{s}-A(\boldsymbol{b}(s))_{i, :} D(\boldsymbol{c}(s), \tau)\right\|^{2}$为常数, c^*(s)、b^*(s)和τ^*与更新前的c(s)、b(s)和τ完全相似, 只是用c_k^*(s)代替扩频码向量c(s)的第k个值c_k(s), b_n^*(s)代替信息序列向量b(s)的第n个值b(i_n^s), τ^*代替时延τ。

2.4 序列拼接算法

由于算法是在对S个分段数据做独立抽取后获得的结果, 因此需要将获得的b(s)和c(s)进行拼接和重构(1≤s≤S), 得到完整的扩频码及信息序列。针对不同的片段, 在算法的执行过程中, 首先根据估计到的时延τ对基带信号进行码片同步, 得到用户的扩频向量c(s)。在扩频码的拼接过程中, 利用片段间重叠部分的相似性来提高序列的估计性能。

对于信息序列b(s), 在s为不同值时, 由于在一个完整的扩频周期内包含多个信息序列, 致使每个向量对应的信息符号不是完整的信息序列。可以通过对抽取得到的b(s)按照下标进行排序, 得到完整的信息序列。在算法的实现过程中, 同样可以对重复抽取得到的样本值做均值处理来提高估计性能。

由上文的推导过程可以看出, 本节的算法对于短码和长码直扩信号都能适应, 其具体实现步骤如下:

步骤1  根据已知的扩频周期, 将接收到的基带信号按照图 2所示进行重叠分段, 得到S个分段数据R_s(1≤s≤S)。

步骤2  初始化数据。根据各参数的先验分布, 得到Θ$\triangleq${b(s), c(s), τ, σ²}的初始值Θ⁽⁰⁾, 设置σ²的初始方差, 并令迭代变量n=1。

步骤3  针对每一个片段数据, 在第n次迭代过程中, 首先利用混合MH算法对b(s), c(s)和τ进行抽样更新; 然后采用UKF方法对σ²进行估计。

步骤4  第n次迭代结束, 进入n+1次迭代, 重复步骤3, 直至完成所有迭代。

步骤5  重复步骤3~步骤4直到获得所有的序列片段。

最后根据给出的序列拼接方法, 估计出完整的扩频码和信息序列。

3 仿真实验结果与分析

本文算法是在长码直扩信号条件下推导出的, 但由分段模型和推导过程可知, 该算法对短码的直扩信号同样适用。下文通过2个仿真实验分别对短码和长码直扩信号的参数估计性能进行仿真验证。算法共进行500次迭代运算, 由于算法收敛需要一个过程, 为了增加仿真结果的可信度和准确性, 取后400次的迭代结果累计进行相应的运算。

仿真实验1   对短码直扩信号的参数估计性能验证

仿真条件:扩频码为R=63位随机序列, 调制样式为BPSK, 码片速率10 MHz, 符号速率为10 MHz/63=158.7 kHz, 扩频增益P=63, 仿真数据为N=300个扩频周期, 分段长度H=8, 重叠长度H₀=H－1, 则每个扩频周期分为S=56个片段, 时延τ从[0, T_s]上随机选取, 并以码片速率进行采样, 噪声为零均值加性复高斯白噪声, 输入信噪比变化范围为－15 dB到15 dB。

图 3直接给出SNR为－9 dB时, 由MCMC算法所得的扩频码。由图 3可以看出, 该算法能在低信噪比条件下较好地适应短码直扩信号。

图 4给出在不同分段长度H和不同周期数目N时, 算法得出的扩频码及信息序列估计性能随输入信噪比变化的曲线。参数设置为:N分别取100、200、300;分段长度H分别取8、16, 且重叠部分均为H₀=H－1。从图 4可以看出, 对于短码直扩信号的扩频码及信息序列盲估计, 重叠分段的MCMC-UKF算法利用300个扩频周期的数据和8个码片周期的分段长度, 就能在－12dB的信噪比条件下, 使得估计扩频码与真实扩频码的相似度超过0.95。该算法性能与使用的周期数目N成正比, 在具有足够数据窗的情况下, 即使具有更低的信噪比, 该算法依然适应, 但随之而来的是成倍增加的计算量。另外, 由图 4还可以看出, 随着分段长度的减少, 算法性能会随之提高, 但算法的计算量也同样会随着增加。因此, 在实际的应用中, 需要考虑性能与计算复杂度的折中。

仿真实验2   对长码直扩信号的参数估计性能验证

仿真条件:扩频码为R=63位随机序列, 调制样式为BPSK, 码片速率10 MHz, 符号速率为10 MHz/30=333.3 kHz, 扩频增益P=30, 仿真数据为N=300个扩频周期, 分段长度H=8, 重叠长度H₀=H－1, 则每个扩频周期分为S=56个片段, 其他参数与仿真实验1相同。

图 5给出了SNR为－9 dB时, 采用后400次迭代值进行的时延参数τ的估计, 并以直方图的形式表示。从图 5可以看出, 时延估计在真值的位置上具有最大的后验概率, 因此本文算法能够较好地在低信噪比条件下对多径时延参数进行估计, 仿真时τ的真值为{0.6T_s}。

图 6给出在不同分段长度H和不同周期数目N时, 长码信号扩频码及信息序列的估计性能随输入信噪比变化的曲线。参数设置为:N分别取100、200、300;分段长度H分别取8、16, 且重叠部分均为H₀=H－1。由图 6可以看出, 对于长码直扩信号, 信息符号在每段持续时间H内都有一定发生跳变的可能, 这样就会形成干扰, 所以当参数条件相同时, 长码信号的扩频码及信息序列估计性能明显低于短码信号。随着分段长度的增加, 信息符号发生跳变的概率也随之增大, 扩频码及信息序列的估计性能就会相应降低。利用300个扩频周期的数据和8个码片周期的分段长度, 只能在－10 dB的信噪比条件下使得估计扩频码与真实扩频码的相似度超过0.95。因此, 针对长码直扩信号中的信息调制影响, 如要得到与短码相当的估计性能, 就需要更高的信噪比或更多的时间窗。

4 结束语

本文针对长码直扩信号的扩频码及信息序列盲估计问题, 结合重叠分段和MCMC算法的思想, 提出基于重叠分段MCMC-UKF的盲估计算法。该算法能够实现对短码和长码信号的有效估计, 并且不受扩频序列类型的限制。仿真实验结果表明, 该算法具有适应较低信噪比的能力, 且能获得较好的估计性能。