时延估计是现代信号处理的重要内容，广泛应用于无源定位、雷达声纳、生物医学等领域.经过几十年的发展，形成了以广义互相关(GCC, generalized cross-correlation)方法为代表的时延估计方法^[1]，但这类方法对窄带弱信号估计精度较低.随后各种新的估计方法陆续被提出，如基于高阶累积量的估计方法^[2]；基于自适应滤波的估计方法^[3].在某些特定的背景下，这些方法能使时延估计精度提升，但对窄带弱信号估计精度低的问题仍未妥善解决.

近年来，压缩感知^[4] (CS, compressed sensing)理论兴起.待估参数在信号的某种表示域内具有稀疏性，是应用该理论的根本保证.随后，信道估计^[5]、频率估计^[6]等问题陆续被建模为信息缺失条件下信号的重构问题，应用该理论取得了较理想的结果.受到这些研究工作的启发，稀疏重构的思想也可以为改进窄带弱信号条件下的时延估计问题提供一种新的思路.

1 信号处理模型

假设空间分布的两观测站，采集到的信号y₁(n)和y₂(n)可以表示为

$ \left. \begin{array}{l} {y_1}\left( n \right) = {r_1}s\left( n \right) + {n_1}\left( n \right)\\ {y_2}\left( n \right) = {r_2}s\left( {n - \tau } \right) + {n_2}\left( n \right) \end{array} \right\} $

(1)

其中：s(n)表示目标源信号；r₁和r₂分别表示对应路径的衰减；τ表示信号到达两观测站的时延；n₁(n)和n₂(n)为加性高斯白噪声，n₁(n)和n₂(n)不相关且分别与源信号s(n)不相关，n=0, 1, …, N-1，N为信号采样点数.

对式(1)中两路信号分别做N点快速傅里叶变换(FFT, fast fourier transformation)，得其频谱：

$ \left. \begin{array}{l} {Y_1}\left( k \right) = {r_1}S\left( k \right) + {N_1}\left( k \right)\\ {Y_2}\left( k \right) = {r_2}S\left( k \right){{\rm{e}}^{ - {\rm{j}}\frac{{2{\rm{ \mathsf{ π} }}}}{N}k\tau }} + {N_2}\left( k \right) \end{array} \right\} $

(2)

其中，k=0, 1, …, N-1，S(k)=FFT[s(n)，Y_i(k)=FFT[y_i(n)](i=1, 2)，N_i(k)=FFT[n_i(n)](i=1, 2).

用X(k)=表示两信号的互相关谱(简称互谱)：

$ X\left( k \right) = {Y_1}\left( k \right)Y_2^ * \left( k \right) = {r_1}{r_2}{\left| {S\left( k \right)} \right|^2}{{\rm{e}}^{{\rm{j}}\frac{{2{\rm{ \mathsf{ π} }}}}{N}k\tau }} + N\left( k \right) $

(3)

其中：上标*表示取共轭操作；N(k)=为总噪声，根据加性高斯白噪声的特性，N(k)=也为加性高斯白噪声.

通过式(3)可以看出，待估时延参数τ存在于互谱的相位信息中.记ω=2πτ/N为待估中间参数，根据时频域的对偶关系，如从时域采样数据中估计出正弦信号频率一样^[6]，可以将X(k)=看作是信号在频域的采样数据，利用其对时延中间参数ω进行估计.进一步综合不同频点的数据，将式(3)表示成矢量形式：

$ \mathit{\boldsymbol{X}} = \mathit{\boldsymbol{Ad}} + \mathit{\boldsymbol{N}} $

(4)

其中：X=[X(1), X(2), …, X(N)]^T表示频域采样信号矢量数据，A=diag{r₁r₂[|S(1)|², |S(2)|², …, |S(N)|²]}表示信号幅度因子，d=[1, e^jω, …, e^j(N-1)ω]^T表示时延导向矢量，N =[N(1), N(2), …, N(N)]^T表示噪声矢量.

注意到待估参数在整个“参数域”具有天然的稀疏性，所以可以将式(4)稀疏表示.在允许一定量化误差的前提下，利用穷举方法得到一个完备集，使得待估参数构成该集合的一个较小子集.规模为N×L(1≪L)的超完备字典D可表示为

$ D = \left[ {{d_1}, \cdots ,{d_i}, \cdots ,{d_L}} \right] $

(5)

其中，d_i=[1, e^jω_i, …, e^j(N-1)ω_i]^T表示字典D中的第i个基.

根据式(5)，可将式(4)稀疏表示为

$ \mathit{\boldsymbol{X}} = \mathit{\boldsymbol{Dz}} + \mathit{\boldsymbol{N}} $

(6)

其中z为L×1维稀疏系数，在z中只有1个非零值，该非零值的位置代表信号待估参数，其大小对应待估信号幅度.此时，时延估计问题转化成在字典中挑选一个基函数并优化对应幅度，实现最佳拟合.由于1≪L，可以将式(6)表示为以下优化模型：

$ \mathop {\min {{\left\| \mathit{\boldsymbol{z}} \right\|}_0}}\limits_\mathit{\boldsymbol{z}} ,{\rm{s}}{\rm{.}}\;{\rm{t}}{\rm{.}}\;\;{\left\| {\mathit{\boldsymbol{X}} - \mathit{\boldsymbol{Dz}}} \right\|_2} < \varepsilon $

(7)

直接求解式(7)是一个NP-hard问题^[6]，可将目标函数中的l₀范数用l₁范数来代替，从而转化为一个如式(8)所示凸优化问题，该问题可以被方便地求解.

$ \mathop {\min {{\left\| \mathit{\boldsymbol{z}} \right\|}_1}}\limits_\mathit{\boldsymbol{z}} ,{\rm{s}}{\rm{.}}\;{\rm{t}}{\rm{.}}\;\;{\left\| {\mathit{\boldsymbol{X}} - \mathit{\boldsymbol{Dz}}} \right\|_2} < \varepsilon $

(8)

2 基于稀疏重构的时延估计算法

在信号模型(6)~(8)中，信号矢量X的长度为采集信号点数N，当N较大时，矩阵规模也随之变大，影响计算速度.考虑到待估参数维度较小，可设信号阶数为M(M的选择应大于待估参量个数)，将N点数据交叠划分为I=N-M+1个“快拍”数据按列排开，每个“快拍”为M×1维.此时信号模型变为

$ \mathit{\boldsymbol{\bar X}} = \mathit{\boldsymbol{\bar Dz}} + \mathit{\boldsymbol{N}} $

(9)

其中：X∈C^M×I表示多“快拍”数据矩阵，D∈C^M×L表示修正的稀疏字典，分别可表示为

$ \mathit{\boldsymbol{\bar X}} = \left[ {\begin{array}{*{20}{c}} {X\left( 1 \right)}&{X\left( 2 \right)}& \cdots &{X\left( {N - M + 1} \right)}\\ {X\left( 2 \right)}&{X\left( 3 \right)}& \cdots &{X\left( {N - M + 2} \right)}\\ \vdots&\vdots&\vdots&\vdots \\ {X\left( M \right)}&{X\left( {M + 1} \right)}& \cdots &{X\left( N \right)} \end{array}} \right] $

(10)

$ \mathit{\boldsymbol{\bar D}} = \left[ {\begin{array}{*{20}{c}} 1&1& \cdots &1\\ {{{\rm{e}}^{{\rm{j}}{\omega _1}}}}&{{{\rm{e}}^{{\rm{j}}{\omega _2}}}}& \cdots &{{{\rm{e}}^{{\rm{j}}{\omega _L}}}}\\ \vdots&\vdots&\vdots&\vdots \\ {{{\rm{e}}^{{\rm{j}}\left( {M - 1} \right){\omega _1}}}}&{{{\rm{e}}^{{\rm{j}}\left( {M - 1} \right){\omega _2}}}}& \cdots &{{{\rm{e}}^{{\rm{j}}\left( {M - 1} \right){\omega _L}}}} \end{array}} \right] $

(11)

2.1 数据预处理

多快拍稀疏分解随着快拍数增加运算量显著增大，对信号数据矩阵预先做奇异值分解(SVD, singular value decomposition)SVD分解，在信号子空间中进行后续操作，可以降低算法复杂度.对式(9)中构造的信号矩阵X做SVD分解：

$ \mathit{\boldsymbol{\bar X}} = \mathit{\boldsymbol{U \boldsymbol{\varSigma} }}{\mathit{\boldsymbol{V}}^{\rm{H}}} = \left[ {{\mathit{\boldsymbol{U}}_S},{\mathit{\boldsymbol{U}}_N}} \right]\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}{\mathit{\boldsymbol{V}}^{\rm{H}}} $

(12)

其中：U和V为酉矩阵，Σ表示以X的奇异值为对角线元素的矩阵.

Σ的前K个较大奇异值对应的特征向量所张成的子空间对应M×K维的信号子空间.令P_K=[I_K, O_K×(L-K)]，其中I_K为K×K维的单位矩阵，O_K×(L-K)为K×(L-K)维的零矩阵.需要说明的是，本模型没有考虑多径时延情况，待估时延参数有且仅有一个取值，所以可以设K=1.令X_S=XVP_K, z_S= zVP_K, N_S= NVP _K，则在信号子空间下可以得到：

$ {{\mathit{\boldsymbol{\bar X}}}_S} = \mathit{\boldsymbol{\bar D}}{\mathit{\boldsymbol{z}}_S} + {\mathit{\boldsymbol{N}}_S} $

(13)

对比式(6)和式(13)可以看出，SVD操作使得原问题的规模由M×L降为M×K.右乘矩阵为矩阵的初等列变换，因此z_S与z的行稀疏特性没有改变，因此不影响待估参数估计结果，且经过SVD得到的信号子空间中，多快拍数据中信号能量得到累积，有利于时延参数的估计.

2.2 算法步骤

正交匹配追踪(OMP, orthogonal matching pursuit)算法通过每次在稀疏字典中搜索与残余信号相关性最大的原子匹配出K个原子，求出稀疏系数确定待估参数.在第m次迭代过程中计算当前残差向量r^m-1与字典中各原子的相关性，选出最佳匹配原子，加入到已选出的原子集合Ω^m-1中，并计算原信号X_S在集合Ω^m中的正交投影，得到投影系数p^m，并更新下一次迭代的残差r^m.其过程可以通过式(14)所示模型表示：

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}^m} = {\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}^{m - 1}} \cup \mathit{\boldsymbol{\bar d}}_i^m,i \in \left[ {1,L} \right]}\\ {\left. \begin{array}{l} \mathop {\min }\limits_{{\mathit{\boldsymbol{p}}^m}} {\left\| {{{\mathit{\boldsymbol{\bar X}}}_S} - {\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}^m}{\mathit{\boldsymbol{p}}^m}} \right\|_2}\\ {\mathit{\boldsymbol{r}}^m} = {{\mathit{\boldsymbol{\bar X}}}_S} - {\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}^m}{\mathit{\boldsymbol{p}}^m} \end{array} \right\}} \end{array} $

(14)

综上，所提算法具体步骤如下：

步骤1  利用接收数据按照式(2)~(3)构造互谱；

步骤2  利用互谱按照式(10)和(11)构造数据矩阵X和稀疏字典D；

步骤3  对X做SVD分解，提取信号子空间；

步骤4  在信号子空间中通过式(14)迭代得到中间参数ω的估值；

步骤5  根据ω=2πτ/N的关系，获得时延τ的估计值.

2.3 计算复杂度分析

所提算法中用到了OMP算法，其计算复杂度为O(M²L)^[6]，由于引入SVD分解，将搜索空间压缩，运算量降至O(K²L).而SVD本身的复杂度为O(KM²).通常情况下，稀疏度K是一个很小的数，与待估计辐射源个数相关，如在本文中K=1，此时所提算法计算复杂度为O(L)，即在K确定时候，运算量与冗余字典规模L成正比.

直接在时域计算互相关^[1]的计算复杂度为O(N²)，其中N为信号采样点数，由于可将时域计算等效到频域计算，采用FFT快速算法，可将复杂度降至O(NlbN)；自适应时延估计的最小均方算法(LMS, least mean square)LMS^[3]计算复杂度为O(PN)，其中P为滤波器的阶数，P的选取需满足大于整数采样周期时延的限定条件，为保证算法的通用性，取P=N，则复杂度为O(N²)；四阶累积量算法CUM4^[2]计算复杂度为O(N⁴).

综上分析，所提算法计算复杂度略大于GCC方法，远小于自适应方法和高阶累积量方法.

3 仿真验证

首先定义时延估计的精度指标——均方根误差^[6](RMSE，root mean square error)，可表示为

$ R = \sqrt {\frac{1}{M}\sum\limits_{m = 1}^M {{{\left( {\hat \tau - \tau } \right)}^2}} } $

(15)

其中：M为蒙特卡洛仿真次数，${\hat{\tau }}$和τ分别为时延的估计值和真实值.

假设信源采用二进制相移键控(BPSK, binary phase shift keying)编码，中频70 kHz，采样率为200 kHz，真实时延设为31倍采样间隔，即τ=1.55×10^-4 s.分别设计实验，验证信号带宽、信噪比对时延估计精度的影响，并与克拉美罗下界(CRLB, Cramer-Rao lower bound)进行比较.

仿真实验1  仿真弱信号条件下，RMSE随带宽B的变化关系.信噪比(SNR, signal-to-noise ratio，)取2 dB，积累时间2 ms，带宽在2~20 kHz范围内，以2 kHz为间隔取值，每个带宽取值进行1 000次蒙特卡洛仿真实验，统计RMSE，结果如图 1所示.

“GCC”代表广义互相关方法；“LMS”代表自适应时延估计方法；“CUM4”代表四阶累积量方法；“CS”代表所提的基于稀疏重构时延估计方法；“CRLB”代表时延估计克拉美罗下界.从图 1中可以看出，在弱信号条件下，CUM4精度最差，LMS方法的估计精度与GCC方法相当，所提方法估计精度明显高于对比方法，且精度提升1倍以上.

仿真实验2  仿真窄带条件下，RMSE随信噪比的变化关系.基本条件设置同仿真实验1，带宽设为4 kHz，积累时间为2 ms，信噪比从0~20 dB，以2 dB为间隔取值，每个值进行1 000次蒙特卡洛实验，统计RMSE，结果如图 2所示.

在窄带条件下，CUM4精度最差，LMS方法的估计精度与GCC方法相当，所提方法估计精度明显高于对比方法，且精度提升1倍以上. 图 1和图 2所呈现的结果与理论分析相符.

此外，在计算复杂度方面，按照仿真条件，采样点数N=400，在LMS算法中设滤波器阶数P=N=400，对于所提算法，取L=10N，表示稀疏字典按照采样间隔的1/10划分.这样的条件下，计算复杂度按照从小到大依次是：GCC，CS，LMS和CUM4，所需的乘法运算次数分别为：3.46×10³，4×10³，1.6×10⁵和2.56×10¹⁰.该结果和2.3中的分析一致.

综合以上实验结果，相比于传统方法，该算法将窄带弱信号时延估计精度提高了约1倍.

4 结束语

随着信息技术的飞速发展，电磁环境日趋复杂，实际工程中往往会面临窄带信号、微弱信号等情况，以GCC为代表的传统方法在这些条件下时延估计精度低，从而限制了其应用.针对该问题，充分挖掘信号自身的结构信息，利用信号时延参数在“时延域”具有稀疏性这一隐含条件，提出了一种基于稀疏重构的窄带弱信号时延估计算法.理论分析和仿真实验验证了算法的正确性和有效性，与传统方法相比，该算法将窄带弱信号时延估计精度提高约1倍，具有一定的工程应用价值.