一种基于改进后的SAGE算法的信道参数估计方法

摘要

一种基于改进后的SAGE算法的信道参数估计方法，包括以下步骤：利用连续干扰消除原理对接收信号进行处理后，使用非相干最大似然估计法进行初始化估计；将初始化的参数结果合成模拟的接收信号；将原接收信号与得到的模拟接收信号之差定义为欧式距离，计算令欧式距离最小的部分干扰因子；代入部分干扰因子，利用部分干扰消除的方式对接收信号进行处理后，使用非相干最大似然估计法得到最终初始化结果；计算Z函数，代入最终初始化结果进行迭代更新，找到使Z函数取最大值的点作为参数估计结果。本发明的方法计算精度高，在提高整体收敛性能的条件下对实际信道信道特征作以较为准确的估计，计算效果好，有效解决了现有技术的不足。

说明书

技术领域

本发明属于无线通信技术领域，具体涉及一种基于改进后的SAGE算法的信道参数估计方法，可用于在多径信号情况下以较好的收敛性能估计信道参数，表征实际信道特性。

背景技术

随着无线通信技术的不断发展，无线信道作为无线通信系统的基础，对通信系统标准制定、无线新系统设计、无线设备以及网络布局的选择具有重要意义，因此需要提取尽可能逼近实际传播环境的准确的信道参数。由于无线信道的动态特性，无线传播环境中的参数变化很快。为了保证高质量通信信号的高效传输，基于一定算法的信道参数估计起着非常关键的作用，因此研究信道参数估计的方法具有一定的意义。

Schmidt教授提出了基于多重信号分类(MUSIC)算法的信道参数估计方法，但是MUSIC算法对来波信号的要求较高，不适合对相干信号的估计。(R.O.Schmidt,Multiple emitter location and signal parameter estimation[J],IEEE Trans.AP,1986,34(2):276-280)。Roy教授提出了基于旋转子空间不变法(ESPRIT)的参数估计方法，常用于对入射波到达角或者谐波频率的估计，但是该算法的误差较大(Roy R,Kailath T.ESPRIT-estimation of signal parameters via rotational Invariance techniques[J].IEEE Trans.AP,1989,37(7):984-995P)。Fleury教授提出了基于空间交替通用期望最大化(SAGE)算法的信道估计方法，但是此算法的精确度、收敛速度以及迭代次数等整体性能有待提高 (Fleury B H,Tschudin M,Heddergott R,et al.Channel parameter estimation in mobile radio environments using the SAGE algorithm[J].Selected Areas in Communications IEEE Journal on,1999,17(3):434-450)。

发明内容

针对现有技术的不足，本发明旨在提供一种基于改进后的SAGE算法的信道参数估计方法，在提高整体收敛性能的条件下对实际信道信道特征作以较为准确的估计。

为了实现上述目的，本发明采用的技术方案如下：

S1、利用连续干扰消除原理对接收信号进行处理后，使用非相干最大似然估计法进行初始化估计；

S2、将初始化的参数结果合成模拟的接收信号；

S3、将原接收信号与S2得到的模拟接收信号之差定义为欧式距离，计算令欧式距离最小的部分干扰因子；

S4、代入部分干扰因子，利用部分干扰消除的方式对接收信号进行处理后，使用非相干最大似然估计法得到最终初始化结果；

S5、计算Z函数，代入最终初始化结果进行迭代更新，找到使Z函数取最大值的点作为参数估计结果。

需要说明的是，步骤S1中，利用连续干扰消除原理对接收信号进行处理后，使用非相干最大似然估计法进行初始化估计按以下进行：

首先，我们将发射天线阵的发射信号定义为由此得到接收天线阵接收到的第l簇响应为：

$s(t;{\theta }_l)=\exp \left\{j2{\mathrm{\pi \upsilon }}_{l}t\right\}{C}_2\left({\Omega }_{2,l}\right)A_1{C}_1^H\left({\Omega }_{1,l}\right)u(t-{\tau }_l)$

其中包含参数集θ_l＝[Ω_1,l,Ω_2,l,τ_l,ν_l,A_l]表示第l条簇散射波的参数矢量集合。分别为离开角、到达角、传播时延、多普勒频率和极化幅度矩阵。矢量矩阵C_k(Ω)表示发/收天线的导向矢量，其计算表达式如下：

$C_k\left(\Omega \right)={\left[f_{k,m}\exp \right\{j2{\mathrm{\pi \lambda }}_0^{-1}(\Omega ·r_{m,k})\};m=1,...,M_k]}^T(k=1,2)$

Ω代表波发或波达角，r_m,k代表天线阵元的位置向量，Ω·r_m,k体现了波程差的计算。为了在实测中使用ULA天线，必须采用特殊的天线摆放方式，这样可得波程差以及收发端导向矢量计算式为：

C_k＝[1,fm*exp(j2π*Δ/λ),fm*exp(j2π*2Δ/λ),...,fm*exp(j2π*(N_k-1)Δ/λ)]其中，fm是天线厂方向图。通过上式的计算，可以将得到的导向矢量带入SAGE，计算流程中进行参数提取。仿真时使用全向天线，即fm＝1。

接收端的接收信号即为L条多径信号的加和，并且叠加高斯噪声：

$Y\left(t\right)={[Y_1\left(t\right),...,Y_{M_2}\left(t\right)]}^T=\underset{l=1}{\overset{L}{\Sigma }}s(t;{\theta }_l)+\sqrt{\frac{N_0}{2}}W\left(t\right)$

其中，为随机白噪声，满足E[W_m(t)·W_m(t)′]＝2δ(t-t′)。定义时间窗函数为：

$q_T\left(t\right)=\left(\begin{array}{c}1,t\in [0,T)\\ 0,otherwise\end{array}\right)$

I个循环周期为一个快拍，在一个快拍内，令q_k,m(t)(k＝1代表发射端，k＝2代表接收端)表示发/收端的时间窗函数，由第m个发射阵元发射信号的时间窗为：

$q_{1,m}\left(t\right)=\underset{i=1}{\overset{I}{\Sigma }}{q}_{T_t}(t-t_{i,m})$

其中，t_i,m＝(i-1)T_cy+(m-1)T_t,m＝1,...,M₁,i＝1,...,I,发射天线阵发射的信号则为u(t)＝q₁u(t),其中

同理可以计算接收端第n个阵元接收信号的时间窗函数为：

$q_{2,n}\left(t\right)=\underset{i=1}{\overset{I}{\Sigma }}\underset{m=1}{\overset{M_1}{\Sigma }}{q}_{T_{sc}}(t-t_{i,n,m})$

其中，t_i,n,m＝(i-1)T_cy+(m-1)T_t,m＝1,...,M₁,i＝1,...I,那么它实际上是一个值域为{0,1}的时间开关函数。则接收天线接收到的第l簇散射波形成的响应可以进一步计算为：

s(t；θ_l)＝α_lexp(j2πv_lt)q₂(t)^Tc₂(Ω_2,l)c₁(Ω_1,l)^Tq₁(t-τ_l)u(t-τ_l)

＝α_lexp(j2πv_lt)c₂(Ω_2,l)^TU(t；τ_l)c₁(Ω_1,l)

其中U(t；τ_l)＝q₂(t)q₁(t-τ_l)^Tu(t-τ_l)是维的信号探测矩阵，表征了开关切换情况，用以确定子信道，该探测矩阵内的元素与为：

$u_{n,m}(t;{\tau }_l)=q_{2,n}\left(t\right)q_{1,m}(t-{\tau }_l)u(t-{\tau }_l)=\underset{i=1}{\overset{I}{\Sigma }}{q}_{Tsc}(t-t_{i,n,m})u(t-{\tau }_l)$

其中n＝1,...,M₂,m＝1,...,M₁,i＝1,...,I。

利用连续干扰消除(Successive Interference Cancellation,SIC)方法对第l条径参数进行初始化时，首先将除第l条径的其它径的信号进行滤除。其步骤为：

对τ_l进行最大似然估计，τ_l的初始化值为：

$\widehat{{\tau }_l}\left(0\right)\underset{{\tau }_l}{\arg \max }\left\{\underset{i=1}{\overset{I}{\Sigma }}\underset{n=1}{\overset{Nr}{\Sigma }}\underset{m=1}{\overset{Nt}{\Sigma }}{\left|{\int }_0^{Tsc}{y}^{\left(l\right)}(t+t_{i,n,m})u^{*}(t-{\tau }_l)dt\right|}^2\right\}$

得到Ω_2,l的非相干最大似然估计过程如下：

${\hat{\phi }}_{2,l}\left(0\right)=\underset{{\phi }_{2,l}\left(0\right)}{\arg \max }\left\{\underset{i=1}{\overset{I}{\Sigma }}\underset{m=1}{\overset{N_t}{\Sigma }}{|\underset{n=1}{\overset{N_t}{\Sigma }}{\tilde{c}}_{2,n}{\left({\phi }_{2,l}\right)}^{*}·{\int }_0^{Tsc}{y}^{\left(l\right)}(t+t_{i,n,m})u(t-\hat{\tau }(0\left)\right)dt|}^2\right\}$

得到后，对Ω_1,l进行非相关最大似然估计过程如下：

得到后，对υ_l进行非相干最大似然估计步骤为：

$\widehat{{\upsilon }_l}\left(0\right)=\underset{{\upsilon }_l}{\arg \max }\left\{{\left|\underset{i=1}{\overset{I}{\Sigma }}\underset{m=1}{\overset{N_t}{\Sigma }}{\tilde{c}}_{1,m}{\left({\hat{\phi }}_{1,l}\right)}^{*}\left(\underset{n=1}{\overset{N_r}{\Sigma }}\begin{array}{c}\exp (-j2{\mathrm{\pi \upsilon }}_l{t}_{i,n,m})·{\tilde{c}}_{2,n}{\left({\hat{\phi }}_{2,l}\right(0\left)\right)}^{*}\\ ·{\int }_0^{Tsc}{y}^{\left(l\right)}(t+t_{i,n,m})u^{*}(t-{\hat{\tau }}_l(0\left)\right)dt\end{array}\right)\right|}^2\right\}$

得到后，复幅度α_l的初始化值可以由下式计算得到：

${\hat{\alpha }}_l\left(0\right)={\left[\right|c_2\left({\hat{\phi }}_{2,l}\right(0\left)\right)\left|\right|c_1\left({\hat{\phi }}_{1,l}\right(0\left)\right)\left|IPTsc\right]}^{-1}.$

$\{\underset{i=1}{\overset{I}{\Sigma }}\underset{m=1}{\overset{Nt}{\Sigma }}{\tilde{c}}_{1,m}{\left({\hat{\phi }}_{1,l}\right)}^{*}{\left(\underset{n=1}{\overset{Nr}{\Sigma }}\exp \right(-j2\pi \widehat{{\upsilon }_l}\left(0\right)t_{i,n,m})·{\tilde{c}}_{2,n}({\hat{\phi }}_{2,l}\left)\right)}^{*}·{\int }_0^{Tsc}{y}^{\left(l\right)}(t+t_{i,n,m})u^{*}(t-\widehat{{\tau }_l}(0\left)\right)dt\}$

这样依次将上一次得到的初始化值代入下一个参数的估计计算中，得到一组五维的参数即为信道参数初始化估计值。

需要说明的是，步骤S3中，将原接收信号与步骤S2得到的模拟接收信号之差定为欧式距离，计算令欧式距离最小的部分干扰因子，按以下进行：

首先计算模拟接收天线接收到的第l簇散射波形成的响应，

s(t；θ_l)＝α_lexp(j2πv_lt)q₂(t)^Tc₂(Ω_2,l)c₁(Ω_1,l)^Tq₁(t-τ_l)u(t-τ_l)

＝α_lexp(j2πv_lt)c₂(Ω_2,l)^TU(t；τ_l)c₁(Ω_1,l)

其中U(t；τ_l)＝q₂(t)q₁(t-τ_l)^Tu(t-τ_l)是维的信号探测矩阵，表征了开关切换情况，用以确定子信道，该探测矩阵内的元素与为：

$u_{n,m}(t;{\tau }_l)=q_{2,n}\left(t\right)q_{1,m}(t-{\tau }_l)u(t-{\tau }_l)=\underset{i=1}{\overset{I}{\Sigma }}{q}_{Tsc}(t-t_{i,n,m})u(t-{\tau }_l)$

其中n＝1,...,M₂,m＝1,...,M₁,i＝1,...,I。

然后根据模拟接收到的第l簇散射波形成的响应计算模拟接收信号，接收端的接收信号即为L条多径信号的加和，并且叠加高斯噪声：

$y{\left(t\right)}_{sim}=\underset{l=1}{\overset{L}{\Sigma }}s(t;{\theta }_l)+\sqrt{\frac{N_0}{2}}W\left(t\right)$

其中，为随机白噪声，满足E[W_m(t)·W_m(t)′]＝2δ(t-t′)。y(t)_sim表示由初始化结果模拟出的接收信号。

接着定义一个代价函数，它表示原始接收信号和由初始化结果模拟估计出的接收信号间的欧氏距离：

c_k＝||y(t)-β*y(t)_sim||

y(t)_sim表示由初始化结果模拟出的接收信号。β表示可以控制欧氏距离大小的部分干扰因子。计算使得欧氏距离最小c_k,min的β，记为β′。

需要说明的是，步骤S4中代入令欧氏距离最小的部分干扰因子，利用部分干扰消除的方式对接收信号进行处理后，使用非相干最大似然估计法得到最终初始化结果，按以下进行：

将部分干扰因子代入到串行干扰消除公式当中，得到新的路径信号：

$y^{\left(l\right)}\left(t\right)=y\left(t\right)-\underset{l^{\prime }=1}{\overset{l-1}{\Sigma }}{\beta }_{\min }·s(t;\widehat{{\theta }_{l^{\prime }}})$

接着利用部分干扰消除方式处理接收信号得到新的第l条路径，过程如下：

$y_{sim}^{\left(l\right)}\left(t\right)=y\left(t\right)-\underset{l^{\prime }=1}{\overset{l-1}{\Sigma }}{\beta }^{\prime }·s(t;\widehat{{\theta }_{l^{\prime }}})$

接着利用非相干最大似然估计(NC-ML)，得到最终初始化结果，过程如下：

对τ_l进行最大似然估计，τ_l的初始化值为：

$\widehat{{\tau }_l}\left(0\right)\underset{{\tau }_l}{\arg \max }\left\{\underset{i=1}{\overset{I}{\Sigma }}\underset{n=1}{\overset{Nr}{\Sigma }}\underset{m=1}{\overset{Nt}{\Sigma }}{\left|{\int }_0^{Tsc}{y}_{sim}^{\left(l\right)}(t+t_{i,n,m})u^{*}(t-{\tau }_l)dt\right|}^2\right\}$

得到Ω_2,l的非相干最大似然估计过程如下：

${\hat{\phi }}_{2,l}\left(0\right)=\underset{{\phi }_{2,l}\left(0\right)}{\arg \max }\left\{\underset{i=1}{\overset{I}{\Sigma }}\underset{m=1}{\overset{N_t}{\Sigma }}{|\underset{n=1}{\overset{N_r}{\Sigma }}{\tilde{c}}_{2,n}{\left({\phi }_{2,l}\right)}^{*}·{\int }_0^{Tsc}{y}_{sim}^{\left(l\right)}(t+t_{i,n,m})u(t-\hat{\tau }(0\left)\right)dt|}^2\right\}$

得到后，对Ω_1,l进行非相关最大似然估计过程如下：

得到后，对υ_l进行非相干最大似然估计步骤为：

$\widehat{{\upsilon }_l}\left(0\right)=\underset{{\upsilon }_l}{\arg \max }\left\{{\left|\underset{i=1}{\overset{I}{\Sigma }}\underset{m=1}{\overset{N_t}{\Sigma }}{\tilde{c}}_{1,m}{\left({\hat{\phi }}_{1,l}\right)}^{*}\left(\underset{n=1}{\overset{N_r}{\Sigma }}\begin{array}{c}\exp (-j2{\mathrm{\pi \upsilon }}_l{t}_{i,n,m})·{\tilde{c}}_{2,n}{\left({\hat{\phi }}_{2,l}\right(0\left)\right)}^{*}\\ ·{\int }_0^{Tsc}{y}_{sim}^{\left(l\right)}(t+t_{i,n,m})u^{*}(t-{\hat{\tau }}_l(0\left)\right)dt\end{array}\right)\right|}^2\right\}$

得到后，复幅度α_l的初始化值可以由下式计算得到：

${\hat{\alpha }}_l\left(0\right)={\left[\right|c_2\left({\hat{\phi }}_{2,l}\right(0\left)\right)\left|\right|c_1\left({\hat{\phi }}_{1,l}\right(0\left)\right)\left|IPTsc\right]}^{-1}.$

$\{\underset{i=1}{\overset{I}{\Sigma }}\underset{m=1}{\overset{Nt}{\Sigma }}{\tilde{c}}_{1,m}{\left({\hat{\phi }}_{1,l}\right)}^{*}{\left(\underset{n=1}{\overset{Nr}{\Sigma }}\exp \right(-j2\pi \widehat{{\upsilon }_l}\left(0\right)t_{i,n,m})·{\tilde{c}}_{2,n}({\hat{\phi }}_{2,l}\left)\right)}^{*}·{\int }_0^{Tsc}{y}_{sim}^{\left(l\right)}(t+t_{i,n,m})u^{*}(t-\widehat{{\tau }_l}(0\left)\right)dt\}$

最终估计得到的作为基于改进的SAGE算法的最终初始化结果。

需要说明的是，步骤S4中计算Z函数，代入最终初始化结果进行迭代更新，找到使Z函数取最大值的点作为参数估计结果，按以下进行：

首先，我们计算目标函数Z函数为：

$z({\bar{\theta }}_l;x_l)=c_2{\left({\Omega }_{2,l}\right)}^H{X}_l(t;{\tau }_l,{\upsilon }_l)c_1{\left({\Omega }_{1,l}\right)}^{*}$

其中，矩阵X_l(t；τ_l,υ_l)内的元素为：

$x_l(t;{\tau }_l,{\upsilon }_l)=\underset{i=1}{\overset{I}{\Sigma }}\exp \{-j2{\mathrm{\pi \upsilon }}_l{t}_{i,n,m}\}·{\int }_0^{Tsc}u{(t-{\tau }_l)}^{*}\exp \{-j2{\mathrm{\pi \upsilon }}_{l}t\}{x}_l(t+t_{i,n,m})dt$

接着进行E步骤，计算条件期望。X_l(t)是完备数据集，不能被观测到，所以必须用可以观测到的不完备数据集的接收信号Y(t)和先前估计到的估计值来计算X_l(t)。将对X_l(t)的自然估计定义为在假设下，X_l(t)相对于观测到的数据Y(t)＝y(t)的条件期望。在SAGE算法对某一条径的参数进行估计时，需要计算该条径当前迭代次序下的完备数据集的条件期望，过程如下：

$\widehat{x_l}\left(t\right)=E_{\widehat{{\theta }_k^{\prime }}}\left[X_l\left(t\right)\right|y\left(t\right)]=y\left(t\right)-\underset{l^{\prime }=1,\ne l}{\overset{L}{\Sigma }}s(t;\widehat{{\theta }_{l^{\prime }}^{\prime }})$

式中表示上一次迭代中所估计出的最新参数集。上式的含义为：第l径信号的条件期望等于从实际的接收机接收信号y(t)的迭代估计参数值多重构出的其他L-1条径信号的加和。

E步骤进行完，接着进行M步骤，即求使得目标函数最大化的参数值。将计算得到后，将其带入目标函数的表达式中，分别估计出第l条径的时延、到达角、离开角、多普勒频移以及复振幅五个参数。通过使目标函数Z函数最大化得到待估参数的SAGE算法执行过程中的迭代顺序如下：

${\hat{\tau }}_l^{\prime \prime }=\underset{{\tau }_l}{\arg \max }z({\hat{\phi }}_{1,l}^{\prime },{\hat{\theta }}_{1,l}^{\prime },{\hat{\phi }}_{2,l}^{\prime },{\hat{\theta }}_{2,l}^{\prime },{\tau }_l,{\hat{\upsilon }}_l^{\prime };{\hat{x}}_l)$

${\hat{\upsilon }}_l^{\prime \prime }=\underset{{\upsilon }_l}{\arg \max }z({\hat{\phi }}_{1,l}^{\prime },{\hat{\theta }}_{1,l}^{\prime },{\hat{\phi }}_{2,l}^{\prime },{\hat{\theta }}_{2,l}^{\prime },{\hat{\tau }}_l^{\prime \prime },{\upsilon }_l;{\hat{x}}_l)$

${\hat{\phi }}_{1,l}^{\prime \prime }=\underset{{\phi }_{1,l}}{\arg \max }z({\phi }_{1,l},{\hat{\theta }}_{1,l}^{\prime },{\hat{\phi }}_{2,l}^{\prime },{\hat{\theta }}_{2,l}^{\prime },{\hat{\tau }}_l^{\prime \prime },{\hat{\upsilon }}_l^{\prime \prime };{\hat{x}}_l)$

${\hat{\theta }}_{1,l}^{\prime \prime }=\underset{{\theta }_{1,l}}{\arg \max }z({\hat{\phi }}_{1,l}^{\prime \prime },{\theta }_{1,l},{\hat{\phi }}_{2,l}^{\prime },{\hat{\theta }}_{2,l}^{\prime },{\hat{\tau }}_l^{\prime \prime },{\hat{\upsilon }}_l^{\prime \prime };{\hat{x}}_l)$

${\hat{\phi }}_{2,l}^{\prime \prime }=\underset{{\phi }_{2,l}}{\arg \max }z({\hat{\phi }}_{1,l}^{\prime \prime },{\hat{\theta }}_{1,l}^{\prime \prime },{\phi }_{2,l},{\hat{\theta }}_{2,l}^{\prime },{\hat{\tau }}_l^{\prime \prime },{\hat{\upsilon }}_l^{\prime \prime };{\hat{x}}_l)$

${\hat{\theta }}_{2,l}^{\prime \prime }=\underset{{\theta }_{2,l}}{\arg \max }z({\hat{\phi }}_{1,l}^{\prime \prime },{\hat{\theta }}_{1,l}^{\prime \prime },{\hat{\phi }}_{2,l}^{\prime \prime },{\theta }_{2,l},{\hat{\tau }}_l^{\prime \prime },{\hat{\upsilon }}_l^{\prime \prime };{\hat{x}}_l)$

${\hat{\alpha }}_l=\frac{z({\bar{\theta }}_l;x_l){|}_{\bar{\theta }={\left(\hat{\bar{\theta }}\right)}_l\left(x_l\right)}}{\left|c_2({\hat{\phi }}_{2,l},{\hat{\theta }}_{2,l})\right|\left|c_1({\hat{\phi }}_{1,l},{\hat{\theta }}_{1,l})\right|IPTsc}$

本发明的有益效果是：

本发明的方法计算精度高，在提高整体收敛性能的条件下对实际信道信道特征作以较为准确的估计，计算效果好，有效解决了现有技术的不足。

附图说明

图1为本发明的流程示意图；

图2为本发明在基于改进SAGE算法下的估计性能

具体实施方式

本发明的具体实现步骤如下：

如图1所示，本发明为一种基于改进后的SAGE算法的信道参数估计方法，其特征在于：所述方法包括以下步骤：

S1、利用连续干扰消除原理对接收信号进行处理后，使用非相干最大似然估计法进行初始化估计；

S2、将初始化的参数结果合成模拟的接收信号；

S3、将原接收信号与S2得到的模拟接收信号之差定义为欧式距离，计算令欧式距离最小的部分干扰因子；

S4、代入部分干扰因子，利用部分干扰消除的方式对接收信号进行处理后，使用非相干最大似然估计法得到最终初始化结果；

S5、计算Z函数，代入最终初始化结果进行迭代更新，找到使Z函数取最大值的点作为参数估计结果。

需要说明的是，步骤S1中，利用连续干扰消除原理对接收信号进行处理后，使用非相干最大似然估计法进行初始化估计按以下进行：

首先，我们将发射天线阵的发射信号定义为由此得到接收天线阵接收到的第l簇响应为：

$s(t;{\theta }_l)=\exp \left\{j2{\mathrm{\pi \upsilon }}_{l}t\right\}{C}_2\left({\Omega }_{2,l}\right)A_1{C}_1^H\left({\Omega }_{1,l}\right)u(t-{\tau }_l)$

其中包含参数集θ_l＝[Ω_1,l,Ω_2,l,τ_l,ν_l,A_l]表示第l条簇散射波的参数矢量集合。分别为离开角、到达角、传播时延、多普勒频率和极化幅度矩阵。矢量矩阵C_k(Ω)表示发/收天线的导向矢量，其计算表达式如下：

$C_k\left(\Omega \right)={\left[f_{k,m}\exp \right\{j2{\mathrm{\pi \lambda }}_0^{-1}(\Omega ·r_{m,k})\};m=1,...,M_k]}^T(k=1,2)$

Ω代表波发或波达角，r_m,k代表天线阵元的位置向量，Ω·r_m,k体现了波程差 的计算。为了在实测中使用ULA天线，必须采用特殊的天线摆放方式，这样可得波程差以及收发端导向矢量计算式为：

C_k＝[1,fm*exp(j2π*Δ/λ),fm*exp(j2π*2Δ/λ),...,fm*exp(j2π*(N_k-1)Δ/λ)]其中，fm是天线厂方向图。通过上式的计算，可以将得到的导向矢量带入SAGE，计算流程中进行参数提取。仿真时使用全向天线，即fm＝1。

接收端的接收信号即为L条多径信号的加和，并且叠加高斯噪声：

$Y\left(t\right)={[Y_1\left(t\right),...,Y_{M_2}\left(t\right)]}^T=\underset{l=1}{\overset{L}{\Sigma }}s(t;{\theta }_l)+\sqrt{\frac{N_0}{2}}W\left(t\right)$

其中，为随机白噪声，满足E[W_m(t)·W_m(t)′]＝2δ(t-t′)。定义时间窗函数为：

$q_T\left(t\right)=\left(\begin{array}{c}1,t\in [0,T)\\ 0,otherwise\end{array}\right)$

I个循环周期为一个快拍，在一个快拍内，令q_k,m(t)(k＝1代表发射端，k＝2代表接收端)表示发/收端的时间窗函数，由第m个发射阵元发射信号的时间窗为：

$q_{1,m}\left(t\right)=\underset{i=1}{\overset{I}{\Sigma }}{q}_{T_t}(t-t_{i,m})$

其中，t_i,m＝(i-1)T_cy+(m-1)T_t,m＝1,...,M₁,i＝1,...,I,发射天线阵发射的信号则为u(t)＝q₁u(t),其中

同理可以计算接收端第n个阵元接收信号的时间窗函数为：

$q_{2,n}\left(t\right)=\underset{i=1}{\overset{I}{\Sigma }}\underset{m=1}{\overset{M_1}{\Sigma }}{q}_{T_{sc}}(t-t_{i,n,m})$

其中，t_i,n,m＝(i-1)T_cy+(m-1)T_t,m＝1,...,M₁,i＝1,...I,那么它 实际上是一个值域为{0,1}的时间开关函数。则接收天线接收到的第l簇散射波形成的响应可以进一步计算为：

s(t；θ_l)＝α_lexp(j2πv_lt)q₂(t)^Tc₂(Ω_2,l)c₁(Ω_1,l)^Tq₁(t-τ_l)u(t-τ_l)

＝α_lexp(j2πv_lt)c₂(Ω_2,l)^TU(t；τ_l)c₁(Ω_1,l)

其中U(t；τ_l)＝q₂(t)q₁(t-τ_l)^Tu(t-τ_l)是维的信号探测矩阵，表征了开关切换情况，用以确定子信道，该探测矩阵内的元素与为：

$u_{n,m}(t;{\tau }_l)=q_{2,n}\left(t\right)q_{1,m}(t-{\tau }_l)u(t-{\tau }_l)=\underset{i=1}{\overset{I}{\Sigma }}{q}_{Tsc}(t-t_{i,n,m})u(t-{\tau }_l)$

其中n＝1,...,M₂,m＝1,...,M₁,i＝1,...,I。

利用连续干扰消除(Successive Interference Cancellation,SIC)方法对第l条径参数进行初始化时，首先将除第l条径的其它径的信号进行滤除。其步骤为：

$y^{\left(l\right)}\left(t\right)=y\left(t\right)-\underset{l^{\prime }=1}{\overset{l-1}{\Sigma }}s(t;\widehat{{\theta }_{l^{\prime }}\left(0\right)})$

对τ_l进行最大似然估计，τ_l的初始化值为：

$\widehat{{\tau }_l}\left(0\right)\underset{{\tau }_l}{\arg \max }\left\{\underset{i=1}{\overset{I}{\Sigma }}\underset{n=1}{\overset{Nr}{\Sigma }}\underset{m=1}{\overset{Nt}{\Sigma }}{\left|{\int }_0^{Tsc}{y}^{\left(l\right)}(t+t_{i,n,m})u^{*}(t-{\tau }_l)dt\right|}^2\right\}$

得到Ω_2,l的非相干最大似然估计过程如下：

${\hat{\phi }}_{2,l}\left(0\right)=\underset{{\phi }_{2,l}\left(0\right)}{\arg \max }\left\{\underset{i=1}{\overset{I}{\Sigma }}\underset{m=1}{\overset{N_t}{\Sigma }}{|\underset{n=1}{\overset{N_r}{\Sigma }}{\tilde{c}}_{2,n}{\left({\phi }_{2,l}\right)}^{*}·{\int }_0^{Tsc}{y}^{\left(l\right)}(t+t_{i,n,m})u(t-\hat{\tau }(0\left)\right)dt|}^2\right\}$

得到后，对Ω_1,l进行非相关最大似然估计过程如下：

得到后，对υ_l进行非相干最大似然估计步骤为：

$\widehat{{\upsilon }_l}\left(0\right)=\underset{{\upsilon }_l}{\arg \max }\left\{{\left|\underset{i=1}{\overset{I}{\Sigma }}\underset{m=1}{\overset{N_t}{\Sigma }}{\tilde{c}}_{1,m}{\left({\hat{\phi }}_{1,l}\right)}^{*}\left(\underset{n=1}{\overset{N_r}{\Sigma }}\begin{array}{c}\exp (-j2{\mathrm{\pi \upsilon }}_l{t}_{i,n,m})·{\tilde{c}}_{2,n}{\left({\hat{\phi }}_{2,l}\right(0\left)\right)}^{*}\\ ·{\int }_0^{Tsc}{y}^{\left(l\right)}(t+t_{i,n,m})u^{*}(t-{\hat{\tau }}_l(0\left)\right)dt\end{array}\right)\right|}^2\right\}$

得到后，复幅度α_l的初始化值可以由下式计算得到：

${\hat{\alpha }}_l\left(0\right)={\left[\right|c_2\left({\hat{\phi }}_{2,l}\right(0\left)\right)\left|\right|c_1\left({\hat{\phi }}_{1,l}\right(0\left)\right)\left|IPTsc\right]}^{-1}.$

$\{\underset{i=1}{\overset{I}{\Sigma }}\underset{m=1}{\overset{Nt}{\Sigma }}{\tilde{c}}_{1,m}{\left({\hat{\phi }}_{1,l}\right)}^{*}{\left(\underset{n=1}{\overset{Nr}{\Sigma }}\exp \right(-j2\pi \widehat{{\upsilon }_l}\left(0\right)t_{i,n,m})·{\tilde{c}}_{2,n}({\hat{\phi }}_{2,l}\left)\right)}^{*}·{\int }_0^{Tsc}{y}^{\left(l\right)}(t+t_{i,n,m})u^{*}(t-\widehat{{\tau }_l}(0\left)\right)dt\}$

这样依次将上一次得到的初始化值代入下一个参数的估计计算中，得到一组五维的参数即为信道参数初始化估计值。

需要说明的是，步骤S2中，将初始化的参数结果合成模拟的接收信号，按以下进行：

首先计算模拟接收天线接收到的第l簇散射波形成的响应，

s(t；θ_l)＝α_lexp(j2πv_lt)q₂(t)^Tc₂(Ω_2,l)c₁(Ω_1,l)^Tq₁(t-τ_l)u(t-τ_l)

＝α_lexp(j2πv_lt)c₂(Ω_2,l)^TU(t；τ_l)c₁(Ω_1,l)

其中U(t；τ_l)＝q₂(t)q₁(t-τ_l)^Tu(t-τ_l)是维的信号探测矩阵，表征了开关切换情况，用以确定子信道，该探测矩阵内的元素与为：

$u_{n,m}(t;{\tau }_l)=q_{2,n}\left(t\right)q_{1,m}(t-{\tau }_l)u(t-{\tau }_l)=\underset{i=1}{\overset{I}{\Sigma }}{q}_{Tsc}(t-t_{i,n,m})u(t-{\tau }_l)$

其中n＝1,...,M₂,m＝1,...,M₁,i＝1,...,I。

然后根据模拟接收到的第l簇散射波形成的响应计算模拟接收信号， 接收端的接收信号即为L条多径信号的加和，并且叠加高斯噪声：

$y{\left(t\right)}_{sim}=\underset{l=1}{\overset{L}{\Sigma }}s(t;{\theta }_l)+\sqrt{\frac{N_0}{2}}W\left(t\right)$

其中，为随机白噪声，满足E[W_m(t)·W_m(t)′]＝2δ(t-t′)。y(t)_sim表示由初始化结果模拟出的接收信号。

需要说明的是，步骤S3中，将原接收信号与前面得到的模拟接收信号之差定为欧式距离，计算令欧式距离最小的部分干扰因子，按以下进行：

首先定义一个代价函数，它表示原始接收信号和由初始化结果模拟估计出的接收信号间的欧氏距离：

c_k＝||y(t)-β*y(t)_sim||

y(t)_sim表示由初始化结果模拟出的接收信号。β表示可以控制欧氏距离大小的部分干扰因子。计算使得欧氏距离最小c_k,min的β，记为β′。

需要说明的是，步骤S4中计算Z函数，代入最终初始化结果进行迭代更新，找到使Z函数取最大值的点作为参数估计结果，按以下进行：

首先，我们计算目标函数Z函数为：

$z({\bar{\theta }}_l;x_l)=c_2{\left({\Omega }_{2,l}\right)}^H{X}_l(t;{\tau }_l,{\upsilon }_l)c_1{\left({\Omega }_{1,l}\right)}^{*}$

其中，矩阵X_l(t；τ_l,υ_l)内的元素为：

$x_l(t;{\tau }_l,{\upsilon }_l)=\underset{i=1}{\overset{I}{\Sigma }}\exp \{-j2{\mathrm{\pi \upsilon }}_l{t}_{i,n,m}\}·{\int }_0^{Tsc}u{(t-{\tau }_l)}^{*}\exp \{-j2{\mathrm{\pi \upsilon }}_{l}t\}{x}_l(t+t_{i,n,m})dt$

接着进行E步骤，计算条件期望。X_l(t)是完备数据集，不能被观测到，所以必须用可以观测到的不完备数据集的接收信号Y(t)和先前估计到的估计值来计算X_l(t)。将对X_l(t)的自然估计定义为在假设下，X_l(t)相对于观测到的数据Y(t)＝y(t)的条件期望。在SAGE算法对某一条径的参数进 行估计时，需要计算该条径当前迭代次序下的完备数据集的条件期望，过程如下：

$\widehat{x_l}\left(t\right)=E_{\widehat{{\theta }_k^{\prime }}}\left[X_l\left(t\right)\right|y\left(t\right)]=y\left(t\right)-\underset{l^{\prime }=1,\ne l}{\overset{L}{\Sigma }}s(t;\widehat{{\theta }_{l^{\prime }}^{\prime }})$

式中表示上一次迭代中所估计出的最新参数集。上式的含义为：第l径信号的条件期望等于从实际的接收机接收信号y(t)的迭代估计参数值多重构出的其他L-1条径信号的加和。

E步骤进行完，接着进行M步骤，即求使得目标函数最大化的参数值。将计算得到后，将其带入目标函数的表达式中，分别估计出第l条径的时延、到达角、离开角、多普勒频移以及复振幅五个参数。通过使目标函数Z函数最大化得到待估参数的SAGE算法执行过程中的迭代顺序如下：

${\hat{\tau }}_l^{\prime \prime }=\underset{{\tau }_l}{\arg \max }z({\hat{\phi }}_{1,l}^{\prime },{\hat{\theta }}_{1,l}^{\prime },{\hat{\phi }}_{2,l}^{\prime },{\hat{\theta }}_{2,l}^{\prime },{\tau }_l,{\hat{\upsilon }}_l^{\prime };{\hat{x}}_l)$

${\hat{\upsilon }}_l^{\prime \prime }=\underset{{\upsilon }_l}{\arg \max }z({\hat{\phi }}_{1,l}^{\prime },{\hat{\theta }}_{1,l}^{\prime },{\hat{\phi }}_{2,l}^{\prime },{\hat{\theta }}_{2,l}^{\prime },{\hat{\tau }}_l^{\prime \prime },{\upsilon }_l;{\hat{x}}_l)$

${\hat{\theta }}_{1,l}^{\prime \prime }=\underset{{\theta }_{1,l}}{\arg \max }z({\phi }_{1,l},{\hat{\theta }}_{1,l}^{\prime },{\hat{\phi }}_{2,l}^{\prime },{\hat{\theta }}_{2,l}^{\prime },{\hat{\tau }}_l^{\prime \prime },{\hat{\upsilon }}_l^{\prime \prime };{\hat{x}}_l)$

${\hat{\theta }}_{1,l}^{\prime \prime }=\underset{{\theta }_{1,l}}{\arg \max }z({\hat{\phi }}_{1,l}^{\prime \prime },{\theta }_{1,l},{\hat{\phi }}_{2,l}^{\prime },{\hat{\theta }}_{2,l}^{\prime },{\hat{\tau }}_l^{\prime \prime },{\hat{\upsilon }}_l^{\prime \prime };{\hat{x}}_l)$

${\hat{\phi }}_{2,l}^{\prime \prime }=\underset{{\phi }_{2,l}}{\arg \max }z({\hat{\phi }}_{1,l}^{\prime \prime },{\hat{\theta }}_{1,l}^{\prime \prime },{\phi }_{2,l},{\hat{\theta }}_{2,l}^{\prime },{\hat{\tau }}_l^{\prime \prime },{\hat{\upsilon }}_l^{\prime \prime };{\hat{x}}_l)$

${\hat{\theta }}_{2,l}^{\prime \prime }=\underset{{\theta }_{2,l}}{\arg \max }z({\hat{\phi }}_{1,l}^{\prime \prime },{\hat{\theta }}_{1,l}^{\prime \prime },{\hat{\phi }}_{2,l}^{\prime \prime },{\theta }_{2,l},{\hat{\tau }}_l^{\prime \prime },{\hat{\upsilon }}_l^{\prime \prime };{\hat{x}}_l)$

${\hat{\alpha }}_l=\frac{z({\bar{\theta }}_l;x_l){|}_{\bar{\theta }={\left(\hat{\bar{\theta }}\right)}_l\left(x_l\right)}}{\left|c_2({\hat{\phi }}_{2,l},{\hat{\theta }}_{2,l})\right|\left|c_1({\hat{\phi }}_{1,l},{\hat{\theta }}_{1,l})\right|IPTsc}$

对于本领域的技术人员来说，可根据以上描述的技术方案以及构思， 做出其它各种相应的改变以及变形，而所有的这些改变以及变形都应该属于本发明权利要求的保护范围之内。