CN103763715B - 基于多天线mimo 3d空心椭球的统计信道建模方法 - Google Patents

Abstract

本发明公开了一种合理的综合改进的3D空心椭球的统计信道建模方法，包括建立空心椭球统计信道模型，计算3D空心椭球移动台端方位角和俯仰角的AOA概率密度函数，计算基站端方位角和俯仰角的AOA概率密度函数，计算多普勒频移DS的概率密度函数的步骤，建立了多天线MIMO 3D空间域信道模型，能更加准确灵活方便地估计宏小区和微小区等移动通信环境，有效的提高电磁信号波达角度（AOA：angle of arrival）、波达时间（TOA:time of arrival）以及多普勒频移(DS:Doppler Shift)等信道参数估计的准确性，拓展了空间统计信道模型的研究和应用。

Description

基于多天线MIMO 3D空心椭球的统计信道建模方法

技术领域

本发明属于无线通信技术领域，尤其是涉及一种基于多天线MIMO 3D空心椭球的统计信道建模方法。

背景技术

近年来，在无线通信中对MIMO(multiple-input multiple-output)系统的研究越来越多，建立能够准确描述信道多径效应的无线信道模型，对于研究无线移动通信系统有着重要的作用。

鉴于3D空间域信道模型描述波达信号的精确性，有一些研究员提出了3D空间模型，这些信道模型更加符合实际移动通信环境。K.B.Baltzis提出了一种3D信道模型，基站BS的位置在离地面有一定的高度，而散射体分布空间则是一个2-D平面空间的圆形区域，这种模型更精确地符合宏小区(Macrocell)环境的时空特征。S.J.Nawaz提出一种在基站BS有指向性天线的3D空间信道模型，引入了主瓣宽度为2α的指向性天线。但是在某些特定的市区蜂窝小区移动通信环境下，比如在某个足球体育场上或者在某个大型购物商场里，而发射信号台在某幢高楼上，这使得移动台周围散射体可能很少，甚至于为零。目前尚缺乏针对微小区的较为精确的信道模型。

此外，在移动通信系统性能中，信道容量能够全面的显示MIMO系统的性能特点。在研究MIMO信道容量之前，首先需要研究的是MIMO阵列天线的各阵元间的空间衰落相关性(SFC:spatial fading correlation)。而阵元间相关性主要取决于入射信号的角度扩展，但是现有信道模型中入射信号的到达角度(AOA：angle of arrival)分布都是均匀分布、高斯分布以及拉普拉斯分布，缺乏真实的空间统计信道仿真。

发明内容

为解决上述问题，本发明公开了一种合理的综合改进的3D空心椭球的统计信道建模方法，建立多天线MIMO 3D空间域信道模型，并推导出实际的方位角(AA，azimuth ofarrival)和俯仰角(EA，elevation of arrival)的AOA概率密度函数。

为了达到上述目的，本发明提供如下技术方案：

一种基于多天线MIMO 3D空心椭球的统计信道建模方法，包括如下步骤：

步骤一：建立多天线MIMO 3D空心椭球统计信道的模型，包括基站和移动台，以移动台MS为原点，建立三维坐标系和空心椭球模型，所述椭圆球体内分布有散射体，以移动台MS为中心的椭圆球体内部空心空间内不存在散射体，模型中r_b为基站到某个散射体的距离，φ_b为基站到某个散射体的方位角，β_b为基站到某个散射体的俯仰角，r_m为移动台到某个散射体的距离，φ_m为移动台到某个散射体的方位角，β_m为移动台到某个散射体的俯仰角，D_Los为基站和移动台之间的直达距离，D为基站和移动台之间的水平距离，H为基站离地面的垂直高度，a₁为散射体3D椭圆球体空间的长轴，a₂为散射体3D椭圆球体空间的短轴，ρ_b为r_b在水平面内的投影，b₁为散射体3D椭圆球体内部空心空间的长轴，b₂为散射体3D椭圆球体内部空心空间的短轴；

步骤二：计算到达角度AOA的概率密度函数：

步骤二-1：计算3D空心椭球移动台端方位角和俯仰角的AOA概率密度函数

(1)定义散射体分布函数：

$f(x_m,y_m,z_m)=\left\{\begin{array}{cc}\frac{1}{V}& (x_m,y_m,x_m)\∈I\\ 0& otherwise\end{array}\right.$

式中，I是空心椭圆球体的散射体分布空间，空间体积V为：

$V=V_{he1}-V_{he2}=\frac{2}{3}\mathrm{\π}(a_1^2{b}_1-a_2^2{b}_2)$

(2)通过雅可比式将坐标(x_m，y_m，z_m)转换为(r_m，φ_m，β_m)：

$\begin{array}{c}p(r_m,{\mathrm{\β}}_m,{\mathrm{\φ}}_m)=\frac{f(x_m,y_m,z_m)}{\left|J(x_m,y_m,z_m)\right|}{|}_{\underset{z_m=r_m{\mathrm{sin\β}}_m}{\underset{y_m=r_m{\mathrm{cos\β}}_m{\mathrm{cos\φ}}_m}{x_m=r_m{\mathrm{cos\β}}_m{\mathrm{sin\φ}}_m}}}\\ =\frac{r_m^2{\mathrm{cos\β}}_m}{V}\end{array}$

(3)由上式对r_m进行积分得到p(β_m，φ_m)概率密度函数：

$\begin{array}{c}p({\mathrm{\β}}_m,{\mathrm{\φ}}_m)={\∫}_{r_{\min }}^{r_{\max }}p(r_m,{\mathrm{\β}}_m,{\mathrm{\φ}}_m){\mathrm{dr}}_m\\ =\frac{{\mathrm{cos\β}}_m}{2\mathrm{\π}(a_1^2{b}_1-a_2^2{b}_2)}\[\frac{a_1^3{b}_1^3}{{(b_1^2{\cos }^2{\mathrm{\β}}_m+a_1^2{\sin }^2{\mathrm{\β}}_m)}^{3/2}}\\ -\frac{a_2^3{b}_2^3}{{(b_2^2{\cos }^2{\mathrm{\β}}_m+a_2^2{\sin }^2{\mathrm{\β}}_m)}^{3/2}}\]\\ 0\≤{\mathrm{\φ}}_m\≤2\mathrm{\π},0\≤{\mathrm{\β}}_m\≤\mathrm{\π}/2\end{array}$

其中，

(4)由p(β_m，φ_m)概率密度函数对β_m直接积分得到方位角的边缘密度函数p(φ_m)：

$p\left({\mathrm{\φ}}_m\right)={\∫}_0^{\mathrm{\π}/2}p({\mathrm{\β}}_m,{\mathrm{\φ}}_m){\mathrm{d\β}}_m=\frac{1}{2\mathrm{\π}}$

(5)由p(β_m，φ_m)概率密度函数对φ_m积分得到俯仰角的边缘密度函数p(β_m)：

$\begin{array}{c}p\left({\mathrm{\β}}_m\right)={\∫}_0^{2\mathrm{\π}}p({\mathrm{\β}}_m,{\mathrm{\φ}}_m){\mathrm{d\φ}}_m\\ =\frac{{\mathrm{cos\β}}_m}{(a_1^2{b}_1-a_2^2{b}_2)}\[\frac{a_1^3{b}_1^3}{{(b_1^2{\cos }^2{\mathrm{\β}}_m+a_1^2{\sin }^2{\mathrm{\β}}_m)}^{3/2}}\\ -\frac{a_2^3{b}_2^3}{{(b_2^2{\cos }^2{\mathrm{\β}}_m+a_2^2{\sin }^2{\mathrm{\β}}_m)}^{3/2}}\]\end{array}$

步骤二-2：计算基站端方位角和俯仰角的AOA概率密度函数

(1)计算基站端的方位角和俯仰角的联合分布函数p(β_b，φ_b)：

$p({\mathrm{\β}}_b,{\mathrm{\φ}}_b)=\frac{({\mathrm{\ρ}}_{b_2}^3-{\mathrm{\ρ}}_{b_1}^3)}{2{\mathrm{\πcos}}^2{\mathrm{\β}}_b(a_1^2{b}_1-a_2^2{b}_2)}$

式中，ρ_b1、ρ_b2分别是和在水平面内的投影，而和是在BS端一定角度(β_b，φ_b)情况下与散射体空间相交的距离，其中，

(2)计算基站方位角的边缘密度函数p(φ_b)：

$p\left({\mathrm{\φ}}_b\right)=\∫{\∫}_{S_{\mathrm{\φ}}}\frac{1}{V}{\mathrm{\ρ}}_b{\mathrm{d\ρ}}_b{\mathrm{dz}}_b=\frac{D{\mathrm{cos\φ}}_b}{V}{A}_{\mathrm{\φ}}$

(3)根据散射体区域投影的面积A_φ，计算AOA方位角的边缘密度函数：

$\begin{array}{c}p\left({\mathrm{\φ}}_b\right)=\\ \left\{\begin{array}{c}\frac{3{\mathrm{cos\φ}}_b\[a_1{b}_1(a_1^2-D^2{\sin }^2{\mathrm{\φ}}_b)-a_2{b}_2(a_2^2-D^2{\sin }^2{\mathrm{\φ}}_b)\]}{4a_1{a}_2(a_1^2{b}_1-a_2^2{b}_2)}\\ -{\mathrm{\φ}}_m^{\′}\≤{\mathrm{\φ}}_b\≤{\mathrm{\φ}}_m^{\′}\\ \frac{3b_1{\mathrm{cos\φ}}_b(a_1^2-D^2{\sin }^2{\mathrm{\φ}}_b)}{4a_1(a_1^2{b}_1-a_2^2{b}_2)}\\ {\mathrm{\φ}}_m^{\′}\≤{\mathrm{\φ}}_b\≤{\mathrm{\φ}}_{m}and\\ -{\mathrm{\φ}}_m\≤{\mathrm{\φ}}_b\≤-{\mathrm{\φ}}_m^{\′}\\ \begin{array}{cc}0& otherwise\end{array}\end{array}\right.\end{array}$

其中，sinφ_m＝(a₁/D)和sinφ′_m＝(a₂/D)；

(4)计算垂直面内边缘密度函数：

$p\left({\mathrm{\β}}_b\right)=\left\{\begin{array}{c}\frac{1}{2{\mathrm{\πcos}}^2{\mathrm{\β}}_b(a{1}^{2}b1-a{2}^{2}b2)}\underset{-{\mathrm{\φ}}_{\max }}{\overset{{\mathrm{\φ}}_{\max }}{\∫}}({\mathrm{\ρ}}_{b_2}^3-{\mathrm{\ρ}}_{b_1}^3){\mathrm{d\φ}}_b\\ {\mathrm{\β}}_{\min }\≤{\mathrm{\β}}_b\≤{\mathrm{\β}}_2,\\ {\mathrm{\β}}_4\≤{\mathrm{\β}}_b\≤{\mathrm{\β}}_{\max }\\ \frac{1}{2{\mathrm{\πcos}}^2{\mathrm{\β}}_b(a{1}^{2}b1-a{2}^{2}b2)}\[\underset{-{\mathrm{\φ}}_{\max }}{\overset{{\mathrm{\φ}}_{\max }}{\∫}}({\mathrm{\ρ}}_{b_2}^3-{\mathrm{\ρ}}_{b_1}^3){\mathrm{d\φ}}_b-\\ \underset{-{\mathrm{\φ}}_{\max }^{\′}}{\overset{{\mathrm{\φ}}_{\max }^{\′}}{\∫}}({\mathrm{\ρ}}_{b_2}^{\′3}-{\mathrm{\ρ}}_{b_1}^{\′3}){\mathrm{d\φ}}_b\]\\ {\mathrm{\β}}_2\≤{\mathrm{\β}}_b\≤{\mathrm{\β}}_4\end{array}\right.$

其中，β_min≤β_b≤β_max，且

${\mathrm{\β}}_{\min }=arctan\left(\frac{HD-\sqrt{H^2{a}_1^2+b_1^2(D^2-a_1^2)}}{D^2-a_1^2}\right)$

${\mathrm{\β}}_{max}=arctan\left(\frac{H}{D-a_1}\right)$

${\mathrm{\β}}_1=arctan\left(\frac{H}{D+a_1}\right)$

${\mathrm{\β}}_2=arctan\left(\frac{HD-\sqrt{H^2{a}_2^2+b_2^2(D^2-a_2^2)}}{D^2-a_2^2}\right)$

${\mathrm{\β}}_3=arctan\left(\frac{H}{D+a_2}\right)$

${\mathrm{\β}}_4=arctan\left(\frac{H}{D-a_2}\right);$

步骤三：计算到达时间TOA的联合概率密度函数

(1)计算来波信号的AOA/TOA联合密度函数：

$p(\mathrm{\τ},{\mathrm{\β}}_m,{\mathrm{\φ}}_m)=\frac{p(r_m,{\mathrm{\β}}_m,{\mathrm{\φ}}_m)}{\left|J(r_m,{\mathrm{\β}}_m,{\mathrm{\φ}}_m)\right|}$

其中，

其中J(r_m，β_m，φ_m)为雅可比转换式，且

$\begin{array}{c}J(r_m,{\mathrm{\β}}_m,{\mathrm{\φ}}_m)=|\frac{\∂r_m}{\∂\mathrm{\τ}}{|}^{-1}\\ =\frac{2{(D{\mathrm{cos\β}}_m{\mathrm{cos\φ}}_m-c\mathrm{\τ}+H{\mathrm{sin\β}}_m)}^2}{c(D_{Los}^2-c^2{\mathrm{\τ}}^2-2c\mathrm{\τ}(D{\mathrm{cos\β}}_m{\mathrm{cos\φ}}_m+H{\mathrm{sin\β}}_m\left)\right)}\end{array}$

(2)通过以上式子，计算AOA/TOA的联合密度函数：

$\begin{array}{c}p(\mathrm{\τ},{\mathrm{\β}}_m,{\mathrm{\φ}}_m)\\ =\frac{(c^2{\mathrm{\τ}}^2+D_{Los}^2-2c\mathrm{\τ}(D{\mathrm{cos\β}}_m{\mathrm{cos\φ}}_m+H{\mathrm{sin\β}}_m\left)\right)}{8V{(D{\mathrm{cos\β}}_m{\mathrm{cos\φ}}_m-c\mathrm{\τ}+H{\mathrm{sin\β}}_m)}^4}\\ \×\[c{(D_{Los}^2-c^2{\mathrm{\τ}}^2)}^2{\mathrm{cos\β}}_m\]\end{array}$

其中到达时延c为光速；

对上式进行β_m积分计算移动台MS在方位角的TOA联合分布函数：

$p(\mathrm{\τ},{\mathrm{\φ}}_m)={\∫}_0^{\mathrm{\π}/2}p(\mathrm{\τ},{\mathrm{\β}}_m,{\mathrm{\φ}}_m){\mathrm{d\β}}_m,$

其中，p(τ，φ_m)的闭式表达式如下：

$\begin{array}{c}p(\mathrm{\τ},{\mathrm{\φ}}_m)=\\ \frac{1}{6{(k_1^2+k_2^2-k_3^2)}^{7/2}}\{\frac{-1}{{(k_1-k_3)}^3}\[k_2\sqrt{k_1^2+k_2^2-k_3^2}\[6k_1^5{k}_3+14k_3^6+k_1^4(22k_3^2-5k_4)-19k_3^4{k}_4+k_1^3{k}_3(12k_2^2-50k_3^2+k_4)-\\ 2k_2^4(k_3^2+k_4)+6k_2^2{k}_3^2(3k_3^2+k_4)+k_1^2(24k_3^4+k_2^2(20k_3^2-7k_4)-6k_3^2{k}_4)+k_1{k}_3(6k_2^4-16k_3+29k_3^2{k}_4+k_2^2(-50k_3^2+k))\]+\\ 6k_1{(k_1-k_3)}^3\[2k_3^2(k_3^2-2k_4)+k_1^2(8k_3^2-k_4)+k_2^2(8k_3^2-k_4)\]\arctan h\left(\frac{k_2}{k_1^2+k_2^2-k_3^2}\right)\]-\\ \frac{1}{{(k_2-k_3)}^2}\[\sqrt{k_1^2+k_2^2-k_3^2}\[k_1^4(6k_2{k}_3-4k_2^2-k_4))+\\ (k_2-2k_2){(k_2-k_3)}^2(6k_2^2{k}_3-4k_2{k}_3^2+2k_2{k}_3-3k_3{k}_4)+k_1^2(k_2-k_3)(12k_2^2{k}_3+26k_3^3-10k_3{k}_4+k_2(-20k_3^2+k^4\left)\right)\]+\\ 6k_1{(k_2-k_3)}^2\[2k_3^2(k_3^2-2k_4)+k_1^2(8k_3^2-k_4)+k_2^2(8k_3^2-k_4)\]\arctan h\left(\frac{k_1-k_2+k_3}{k_1^2+k_2^2-k_3^2}\right)\}\end{array}$

其中

k₁＝Dcosφ_m

k₂＝H

k₃＝cτ

$k_4=c^2{\mathrm{\τ}}^2+D_{Los}^2$

对AOA/TOA的联合密度函数进行φ_m积分计算移动台MS在俯仰角的TOA联合分布函数：

$\begin{array}{c}p(\mathrm{\τ},{\mathrm{\β}}_m)={\∫}_0^{2\mathrm{\π}}p(\mathrm{\τ},{\mathrm{\β}}_m,{\mathrm{\φ}}_m){\mathrm{d\φ}}_m\\ =\sqrt{\frac{-k_1+k_2}{k_1+k_2}}\frac{(3k_1^2{k}_2{k}_3+2k_2^3{k}_3+k_1^3{k}_4+k_1^3{k}_4+4k_1{k}_2^2{k}_4)\mathrm{\π}}{{(k_1-k_2)}^4{(k_1+k_2)}^3}\end{array}$

其中

k₁＝Dcosβ_m

k₂＝Hsinβ_m-cτ

$k_3=c^2{\mathrm{\τ}}^2+D_{Los}^2-2c\mathrm{\τ}H{\mathrm{sin\β}}_m$

k₄＝2cτDcosβ_m

步骤四：计算多普勒频移DS的概率密度函数

计算多普勒频移的累积分布函数：

$\begin{array}{c}F\left(\mathrm{\γ}\right)=P_{rob}\left\{{\mathrm{cos\β}}_m{\mathrm{cos\θ}}_m\right\}\\ =P_{rob}\left\{{\mathrm{cos\β}}_m{\mathrm{cos\φ}}_m\right\}\left(p\right({\mathrm{\φ}}_m)=p({\mathrm{\θ}}_m)=\frac{1}{2\mathrm{\π}})\\ =2{\∫}_0^{\mathrm{\π}/2}p\left({\mathrm{\β}}_m\right)\[{\∫}_{\arccos \[\mathrm{\γ}/{\mathrm{cos\β}}_m\]}^{\mathrm{\π}}p\left({\mathrm{\φ}}_m\right){\mathrm{d\φ}}_m\]{\mathrm{d\β}}_m\\ =\frac{1}{\mathrm{\π}}{\∫}_0^{\mathrm{\π}/2}p\left({\mathrm{\β}}_m\right)\[{\∫}_{\arccos \[\mathrm{\γ}/{\mathrm{cos\β}}_m\]}^{\mathrm{\π}}{\mathrm{d\φ}}_m\]{\mathrm{d\β}}_m\end{array}$

其中，P_rob{A}表示A事件发生的概率，γ≡f_DS/f_m，所述其中v为移动台的移动速度，f_c为信号的载波频率，f_m＝vf_c/c是最大多普勒频移，θ_m和β_m分别为移动台MS到达角度的方位角和俯仰角，

对F(γ)求导计算多普勒频移的概率密度函数：

$p\left(\mathrm{\γ}\right)=\frac{\∂F_{\mathrm{\γ}}\left(\mathrm{\γ}\right)}{\∂\mathrm{\γ}}=\frac{1}{\mathrm{\π}}{\∫}_0^{\arccos \left|\mathrm{\γ}\right|}\frac{p\left({\mathrm{\β}}_m\right)}{\sqrt{{\cos }^2{\mathrm{\β}}_m-{\mathrm{\γ}}^2}}{\mathrm{d\β}}_m$

其中，|γ|≤cosβ_m。

具体的，所述的散射体均匀分布在空心椭球体区域内。

与现有多天线MIMO系统建模方法相比，本发明提供的3D空心椭球体信道建模方法，能更加准确灵活方便地估计宏小区和微小区等移动通信环境，有效的提高电磁信号波达角度(AOA：angle of arrival)、波达时间(TOA:time of arrival)以及多普勒频移(DS:Doppler Shift)等信道参数估计的准确性。本发明拓展了空间统计信道模型的研究和应用，为评估多天线MIMO系统空时处理算法和仿真无线通信系统提供了有力的工具。

附图说明

图1为本发明提出的3D空心椭球信道理论模型图；

图2为散射体相交面的3D图；

图3为ρ_b穿过内部空心椭圆时的垂直平面图；

图4为ρ_b没有穿过内部空心椭圆时的垂直平面图；

图5为散射体相交面的投影图；

图6为参数a₁/b₁对MS端EA边缘密度函数的影响示意图(D＝500m，H＝100m，a1＝100m，a₁/b₁＝a₂/b₂)；

图7为不同空心形状对MS端EA边缘密度函数的影响示意图(D＝500m，H＝100m，a1＝b1＝100m)；

图8为参数2a/a1对AA边缘密度函数的影响示意图(D＝500m，H＝100m，a1＝b1＝100m，a2＝b2)；

图9为垂直高度H对BS端EA边缘密度函数的影响示意图(D＝500m，a1＝b1＝100m，a2＝b2＝50m)；

图10为参数a2/a1对BS端EA边缘密度函数的影响示意图(D＝500m，H＝100m，a1＝b1＝100m，a2/b2)；

图11为MS端AA的TOA联合概率密度分布示意图(D＝500m，H＝100m，a₁＝b₁＝100m)；

图12为MS端EA的TOA联合概率密度分布示意图(D＝500m，H＝100m，a₁＝b₁＝100m)；

图13为参数a₁/b₁对多普勒频移的概率密度函数的影响示意图(D＝500m，H＝100m，a₁＝100m，a₁/b₁＝a₂/b₂)；

图14为不同空心形状对多普勒频移的概率密度函数的影响示意图(D＝500m，H＝100m，a1＝b1＝100m)；

具体实施方式

以下将结合具体实施例对本发明提供的技术方案进行详细说明，应理解下述具体实施方式仅用于说明本发明而不用于限制本发明的范围。

本发明提供的统计信道建模方法，需要建立3D空心椭球体信道模型，具体包括如下步骤：

步骤一：建立如图1所示的多天线MIMO 3D空心椭球统计信道的模型，包括基站和移动台，以接收信号移动台MS为中心，以移动台MS为原点，建立三维坐标系和空心椭球模型，假设所有散射体集中分布在整个椭圆球体内，假设在移动台MS周围不存在散射体，这个区域我们假设为一个半椭圆球体(称之为椭圆球体内部空心空间)，图1中各参数符号定义如下：

r_b：基站到某个散射体的距离

φ_b：基站到某个散射体的方位角

β_b：基站到某个散射体的俯仰角

r_m：移动台到某个散射体的距离

φ_m：移动台到某个散射体的方位角

β_m：移动台到某个散射体的俯仰角

D_Los：基站和移动台之间的直达距离

D：基站和移动台之间的水平距离

H：基站离地面的垂直高度

ρ_b：r_b在水平面内的投影

a₁和a₂：散射体3D椭圆球体空间的长轴和短轴

b₁和b₂：散射体3D椭圆球体内部空心空间的长轴和短轴

步骤二：计算到达角度(AOA)的概率密度函数：

步骤二-1：计算3D空心椭球移动台端方位角和俯仰角的AOA概率密度函数

假设散射体分布是均匀分布，其散射体分布函数为：

$f(x_m,y_m,z_m)=\left\{\begin{array}{cc}\frac{1}{V}& (x_m,y_m,x_m)\∈I\\ 0& otherwise\end{array}\right.---\left(1\right)$

式中，I是空心椭圆球体的散射体分布空间，空间体积可以写为：

$V=V_{he1}-V_{he2}=\frac{2}{3}\mathrm{\π}(a_1^2{b}_1-a_2^2{b}_2)---\left(2\right)$

由图1可以看出，任意散射体点反射产生的来波信号AOA的方位角和俯仰角分别为φ_m和β_m，散射体点到移动台MS的距离为r_m，我们通过雅可比式将坐标(x_m，y_m，z_m)转换为(r_m，φ_m，β_m)：

$\begin{array}{c}p(r_m,{\mathrm{\β}}_m,{\mathrm{\φ}}_m)=\frac{f(x_m,y_m,z_m)}{\left|J(x_m,y_m,z_m)\right|}{|}_{\underset{z_m=r_m{\mathrm{sin\β}}_m}{\underset{y_m=r_m{\mathrm{cos\β}}_m{\mathrm{cos\φ}}_m}{x_m=r_m{\mathrm{cos\β}}_m{\mathrm{sin\φ}}_m}}}\\ =\frac{r_m^2{\mathrm{cos\β}}_m}{V}\end{array}---\left(3\right)$

上式(3)对r_m进行积分可以得到p(β_m，φ_m)概率密度函数：

$\begin{array}{c}p({\mathrm{\β}}_m,{\mathrm{\φ}}_m)={\∫}_{r_{\min }}^{r_{\max }}p(r_m,{\mathrm{\β}}_m,{\mathrm{\φ}}_m){\mathrm{dr}}_m\\ =\frac{{\mathrm{cos\β}}_m}{2\mathrm{\π}(a_1^2{b}_1-a_2^2{b}_2)}\[\frac{a_1^3{b}_1^3}{{(b_1^2{\cos }^2{\mathrm{\β}}_m+a_1^2{\sin }^2{\mathrm{\β}}_m)}^{3/2}}\\ -\frac{a_2^3{b}_2^3}{{(b_2^2{\cos }^2{\mathrm{\β}}_m+a_2^2{\sin }^2{\mathrm{\β}}_m)}^{3/2}}\]\\ 0\≤{\mathrm{\φ}}_m\≤2\mathrm{\π},0\≤{\mathrm{\β}}_m\≤\mathrm{\π}/2\end{array}---\left(4\right)$

其中r_max和r_min表达为下式：

$r_{max}=\frac{a_1{b}_1}{\sqrt{(b_1^2{\cos }^2{\mathrm{\β}}_m+a_1^2{\sin }^2{\mathrm{\β}}_m)}}---\left(5\right)$

$r_{\min }=\frac{a_2{b}_2}{\sqrt{(b_2^2{\cos }^2{\mathrm{\β}}_m+a_2^2{\sin }^2{\mathrm{\β}}_m)}}---\left(6\right)$

由式(4)对β_m直接积分可得到方位角的边缘密度函数p(φ_m)：

$p\left({\mathrm{\φ}}_m\right)={\∫}_0^{\mathrm{\π}/2}p({\mathrm{\β}}_m,{\mathrm{\φ}}_m){\mathrm{d\β}}_m=\frac{1}{2\mathrm{\π}}---\left(7\right)$

由式(4)对φ_m积分可得到俯仰角的边缘密度函数p(β_m)：

$\begin{array}{c}p\left({\mathrm{\β}}_m\right)={\∫}_0^{2\mathrm{\π}}p({\mathrm{\β}}_m,{\mathrm{\φ}}_m){\mathrm{d\φ}}_m\\ =\frac{{\mathrm{cos\β}}_m}{(a_1^2{b}_1-a_2^2{b}_2)}\[\frac{a_1^3{b}_1^3}{{(b_1^2{\cos }^2{\mathrm{\β}}_m+a_1^2{\sin }^2{\mathrm{\β}}_m)}^{3/2}}\\ -\frac{a_2^3{b}_2^3}{{(b_2^2{\cos }^2{\mathrm{\β}}_m+a_2^2{\sin }^2{\mathrm{\β}}_m)}^{3/2}}\]\end{array}---\left(8\right)$

步骤二-2：计算基站端方位角和俯仰角的AOA概率密度函数

基站端的方位角和俯仰角的联合分布函数p(β_b，φ_b)为：

$p({\mathrm{\β}}_b,{\mathrm{\φ}}_b)=\frac{({\mathrm{\ρ}}_{b_2}^3-{\mathrm{\ρ}}_{b_1}^3)}{2{\mathrm{\πcos}}^2{\mathrm{\β}}_b(a_1^2{b}_1-a_2^2{b}_2)}---\left(9\right)$

式中，ρ_b1、ρ_b2分别是和在水平面内的投影，而和是在BS端一定角度(β_b，φ_b)情况下与散射体空间相交的距离(如图2所示)。如图5所示，对于给定的β_b，散射体空间相交面在水平内的投影是个椭圆，此时ρ_b1、ρ_b2可以得到：

其中，

$P=1+\frac{a_1^2}{b_1^2}{\tan }^2{\mathrm{\β}}_b---\left(11\right)$

$Q=D{\mathrm{cos\φ}}_b+\frac{a_1^2}{b_1^2}H{\mathrm{tan\β}}_b---\left(12\right)$

$R=D^2-a_1^2+\frac{a_1^2}{b_1^2}{H}^2---\left(13\right)$

而基站方位角的边缘密度函数p(φ_b)表达式为：

$p\left({\mathrm{\φ}}_b\right)=\∫{\∫}_{S_{\mathrm{\φ}}}\frac{1}{V}{\mathrm{\ρ}}_b{\mathrm{d\ρ}}_b{\mathrm{dz}}_b=\frac{D{\mathrm{cos\φ}}_b}{V}{A}_{\mathrm{\φ}}---\left(14\right)$

图3、图4是信道模型在ρ_b平面内的投影，投影椭圆的中心在ρ_b＝Dcosφ_b处，其长轴和短轴分别如图3所示。另外定义sinφ_m＝(a₁/D)和sinφ′_m＝(a₂/D)。如图3示，分别计算2种情况下的散射体区域投影的面积A_φ，

将A_φ的计算公式代入式(14)，AOA方位角的边缘密度函数表达式可以归纳为：

$\begin{array}{c}p\left({\mathrm{\φ}}_b\right)=\\ \left\{\begin{array}{c}\frac{3{\mathrm{cos\φ}}_b\[a_1{b}_1(a_1^2-D^2{\sin }^2{\mathrm{\φ}}_b)-a_2{b}_2(a_2^2-D^2{\sin }^2{\mathrm{\φ}}_b)\]}{4a_1{a}_2(a_1^2{b}_1-a_2^2{b}_2)}\\ -{\mathrm{\φ}}_m^{\′}\≤{\mathrm{\φ}}_b\≤{\mathrm{\φ}}_m^{\′}\\ \frac{3b_1{\mathrm{cos\φ}}_b(a_1^2-D^2{\sin }^2{\mathrm{\φ}}_b)}{4a_1(a_1^2{b}_1-a_2^2{b}_2)}\\ {\mathrm{\φ}}_m^{\′}\≤{\mathrm{\φ}}_b\≤{\mathrm{\φ}}_{m}and\\ -{\mathrm{\φ}}_m\≤{\mathrm{\φ}}_b\≤-{\mathrm{\φ}}_m^{\′}\\ \begin{array}{cc}0& otherwise\end{array}\end{array}\right.---\left(15\right)\end{array}$

如图3所示，基站BS垂直面角β_b的范围为β_min≤β_b≤β_max，其中

${\mathrm{\β}}_{\min }=arctan\left(\frac{HD-\sqrt{H^2{a}_1^2+b_1^2(D^2-a_1^2)}}{D^2-a_1^2}\right)---\left(16\right)$

${\mathrm{\β}}_{max}=arctan\left(\frac{H}{D-a_1}\right)---\left(17\right)$

在计算俯仰角的边缘密度函数时需要用到β₁、β₂、β₃和β₄，其表达式如下：

${\mathrm{\β}}_1=arctan\left(\frac{H}{D+a_1}\right)---\left(18\right)$

${\mathrm{\β}}_2=arctan\left(\frac{HD-\sqrt{H^2{a}_2^2+b_2^2(D^2-a_2^2)}}{D^2-a_2^2}\right)---\left(19\right)$

${\mathrm{\β}}_3=arctan\left(\frac{H}{D+a_2}\right)---\left(20\right)$

${\mathrm{\β}}_4=arctan\left(\frac{H}{D-a_2}\right)---\left(21\right)$

如图3、图4和图5所示，垂直面角β_b可以分为5种情况(β_min≤β_b≤β₁，β₂≤β_b≤β₃，β₃≤β_b≤β₄，β₄≤β_b≤β₅，β₅≤β_b≤β_max)计算出垂直面内的边缘密度函数。当β₂≤β_b≤β₃时，与散射体空间相交面在水平内的投影内部出现空心圆，此时计算边缘密度函数需要减去空心圆的部分。另外当β_b≥β₁，投影部分不在是个完整的圆，出现曲线OO′，此时ρ_b2＝H/tanβ_b(如图4)。综合考虑这5种情况，可推导出垂直面内边缘密度函数：

$p\left({\mathrm{\β}}_b\right)=\left\{\begin{array}{c}\frac{1}{2{\mathrm{\πcos}}^2{\mathrm{\β}}_b(a{1}^{2}b1-a{2}^{2}b2)}\underset{-{\mathrm{\φ}}_{\max }}{\overset{{\mathrm{\φ}}_{\max }}{\∫}}({\mathrm{\ρ}}_{b_2}^3-{\mathrm{\ρ}}_{b_1}^3){\mathrm{d\φ}}_b\\ {\mathrm{\β}}_{\min }\≤{\mathrm{\β}}_b\≤{\mathrm{\β}}_2,\\ {\mathrm{\β}}_4\≤{\mathrm{\β}}_b\≤{\mathrm{\β}}_{\max }\\ \frac{1}{2{\mathrm{\πcos}}^2{\mathrm{\β}}_b(a{1}^{2}b1-a{2}^{2}b2)}\[\underset{-{\mathrm{\φ}}_{\max }}{\overset{{\mathrm{\φ}}_{\max }}{\∫}}({\mathrm{\ρ}}_{b_2}^3-{\mathrm{\ρ}}_{b_1}^3){\mathrm{d\φ}}_b-\\ \underset{-{\mathrm{\φ}}_{\max }^{\′}}{\overset{{\mathrm{\φ}}_{\max }^{\′}}{\∫}}({\mathrm{\ρ}}_{b_2}^{\′3}-{\mathrm{\ρ}}_{b_1}^{\′3}){\mathrm{d\φ}}_b\]\\ {\mathrm{\β}}_2\≤{\mathrm{\β}}_b\≤{\mathrm{\β}}_4\end{array}\right.---\left(22\right)$

其中ρ_b1、ρ_b2、ρ′_b1和ρ′_b2分别是β_b在某一取值范围内相对应的值。

步骤三：计算到达时间(TOA)的联合概率密度函数

对任何散射体反射的信号，从基站BS到移动台MS存在传播路径，其到达时延τ为：

$\mathrm{\τ}=\frac{r_m+r_b}{c}---\left(23\right)$

其中c为光速。通过计算可以到r_b关于(r_m，β_m，φ_m)的表达式：

$\begin{array}{c}{r}_b(r_m,{\mathrm{\β}}_m,{\mathrm{\φ}}_m)\\ =\sqrt{r_m^2+D_{Los}^2-2r_m(D{\mathrm{cos\β}}_m{\mathrm{cos\φ}}_m+H{\mathrm{sin\β}}_m)}\end{array}---\left(24\right)$

将(24)代入(23)中，然后两边平方化简得到r_m(τ，β_m，φ_m)的表达式：

$r_m(\mathrm{\τ},{\mathrm{\β}}_m,{\mathrm{\φ}}_m)=\frac{c^2{\mathrm{\τ}}^2-D_{Los}^2}{2(c\mathrm{\τ}-D{\mathrm{cos\β}}_m{\mathrm{cos\φ}}_m-H{\mathrm{sin\β}}_m)}---\left(25\right)$

来波信号的AOA/TOA联合密度函数表达式为：

$p(\mathrm{\τ},{\mathrm{\β}}_m,{\mathrm{\φ}}_m)=\frac{p(r_m,{\mathrm{\β}}_m,{\mathrm{\φ}}_m)}{\left|J(r_m,{\mathrm{\β}}_m,{\mathrm{\φ}}_m)\right|}---\left(26\right)$

其中J(r_m，β_m，φ_m)为雅可比转换式，有

$\begin{array}{c}J(r_m,{\mathrm{\β}}_m,{\mathrm{\φ}}_m)=|\frac{\∂r_m}{\∂\mathrm{\τ}}{|}^{-1}\\ =\frac{2{(D{\mathrm{cos\β}}_m{\mathrm{cos\φ}}_m-c\mathrm{\τ}+H{\mathrm{sin\β}}_m)}^2}{c(D_{Los}^2-c^2{\mathrm{\τ}}^2-2c\mathrm{\τ}(D{\mathrm{cos\β}}_m{\mathrm{cos\φ}}_m+H{\mathrm{sin\β}}_m))}\end{array}---\left(27\right)$

将(25)和(27)代入(26)中，可以得到AOA/TOA的联合密度函数：

$\begin{array}{c}p(\mathrm{\τ},{\mathrm{\β}}_m,{\mathrm{\φ}}_m)\\ =\frac{(c^2{\mathrm{\τ}}^2+D_{Los}^2-2c\mathrm{\τ}(D{\mathrm{cos\β}}_m{\mathrm{cos\φ}}_m+H{\mathrm{sin\β}}_m\left)\right)}{8V{(D{\mathrm{cos\β}}_m{\mathrm{cos\φ}}_m-c\mathrm{\τ}+H{\mathrm{sin\β}}_m)}^4}\\ \×\[c{(D_{Los}^2-c^2{\mathrm{\τ}}^2)}^2{\mathrm{cos\β}}_m\]\end{array}---\left(28\right)$

对上式进行β_m积分计算可以得到移动台MS在方位角的TOA联合分布函数：

$p(\mathrm{\τ},{\mathrm{\φ}}_m)={\∫}_0^{\mathrm{\π}/2}p(\mathrm{\τ},{\mathrm{\β}}_m,{\mathrm{\φ}}_m){\mathrm{d\β}}_m---\left(29\right)$

其中，p(τ，φ_m)的闭式表达式如下：

$\begin{array}{c}p(\mathrm{\τ},{\mathrm{\φ}}_m)=\\ \frac{1}{6{(k_1^2+k_2^2-k_3^2)}^{7/2}}\{\frac{-1}{{(k_1-k_3)}^3}\[k_2\sqrt{k_1^2+k_2^2-k_3^2}\[6k_1^5{k}_3+14k_3^6+k_1^4(22k_3^2-5k_4)-19k_3^4{k}_4+k_1^3{k}_3(12k_2^2-50k_3^2+k_4)-\\ 2k_2^4(k_3^2+k_4)+6k_2^2{k}_3^2(3k_3^2+k_4)+k_1^2(24k_3^4+k_2^2(20k_3^2-7k_4)-6k_3^2{k}_4)+k_1{k}_3(6k_2^4-16k_3+29k_3^2{k}_4+k_2^2(-50k_3^2+k))\]+\\ 6k_1{(k_1-k_3)}^3\[2k_3^2(k_3^2-2k_4)+k_1^2(8k_3^2-k_4)+k_2^2(8k_3^2-k_4)\]\arctan h\left(\frac{k_2}{k_1^2+k_2^2-k_3^2}\right)\]-\\ \frac{1}{{(k_2-k_3)}^2}\[\sqrt{k_1^2+k_2^2-k_3^2}\[k_1^4(6k_2{k}_3-4k_2^2-k_4))+\\ (k_2-2k_2){(k_2-k_3)}^2(6k_2^2{k}_3-4k_2{k}_3^2+2k_2{k}_3-3k_3{k}_4)+k_1^2(k_2-k_3)(12k_2^2{k}_3+26k_3^3-10k_3{k}_4+k_2(-20k_3^2+k^4\left)\right)\]+\\ 6k_1{(k_2-k_3)}^2\[2k_3^2(k_3^2-2k_4)+k_1^2(8k_3^2-k_4)+k_2^2(8k_3^2-k_4)\]\arctan h\left(\frac{k_1-k_2+k_3}{k_1^2+k_2^2-k_3^2}\right)\}\end{array}---\left(30\right)$

其中

k₁＝Dcosφ_m

k₂＝H

k₃＝cτ

$k_4=c^2{\mathrm{\τ}}^2+D_{Los}^2---\left(31\right)$

对式(28)进行φ_m积分计算可以得到移动台MS在俯仰角的TOA联合分布函数：

$\begin{array}{c}p(\mathrm{\τ},{\mathrm{\β}}_m)={\∫}_0^{2\mathrm{\π}}p(\mathrm{\τ},{\mathrm{\β}}_m,{\mathrm{\φ}}_m){\mathrm{d\φ}}_m\\ =\sqrt{\frac{-k_1+k_2}{k_1+k_2}}\frac{(3k_1^2{k}_2{k}_3+2k_2^3{k}_3+k_1^3{k}_4+4k_1{k}_2^2{k}_4)\mathrm{\π}}{{(k_1-k_2)}^4{(k_1+k_2)}^3}\end{array}---\left(32\right)$

其中

k₁＝Dcosβ_m

k₂＝Hsinβ_m-cτ

$k_3=c^2{\mathrm{\τ}}^2+D_{Los}^2-2c\mathrm{\τ}H{\mathrm{sin\β}}_m$

k₄=2cτDcosβm (33)

步骤四：计算多普勒频移(DS)的概率密度函数

在3D空间信道模型中，移动台MS的移动特性会使信号产生多普勒频移。因此，多普勒频移与传播路径中的角度的关系为：

$f_{DS}=\frac{\mathrm{\υ}}{c}{f}_c{\mathrm{cos\θ}}_m{\mathrm{cos\β}}_m=f_m{\mathrm{cos\θ}}_m{\mathrm{cos\β}}_m---\left(34\right)$

其中，f_c为信号的载波频率，r_m＝vf_c/_c是最大多普勒频移，θ_m和β_m分别为移动台MS到达角度的方位角和俯仰角。我们定义γ≡f_DS/f_m，则有

γ＝cosθ_mcosβ_m (35)

这里我们通过先计算多普勒频移的累积分布函数(CDFs：cumulative densityfunctions):

$\begin{array}{c}F\left(\mathrm{\γ}\right)=P_{rob}\left\{{\mathrm{cos\β}}_m{\mathrm{cos\θ}}_m\right\}\\ =P_{rob}\left\{{\mathrm{cos\β}}_m{\mathrm{cos\φ}}_m\right\}\left(p\right({\mathrm{\φ}}_m)=p({\mathrm{\θ}}_m)=\frac{1}{2\mathrm{\π}})\\ =2{\∫}_0^{\mathrm{\π}/2}p\left({\mathrm{\β}}_m\right)\[{\∫}_{\arccos \[\mathrm{\γ}/{\mathrm{cos\β}}_m\]}^{\mathrm{\π}}p\left({\mathrm{\φ}}_m\right){\mathrm{d\φ}}_m\]{\mathrm{d\β}}_m\\ =\frac{1}{\mathrm{\π}}{\∫}_0^{\mathrm{\π}/2}p\left({\mathrm{\β}}_m\right)\[{\∫}_{\arccos \[\mathrm{\γ}/{\mathrm{cos\β}}_m\]}^{\mathrm{\π}}{\mathrm{d\φ}}_m\]{\mathrm{d\β}}_m\end{array}---\left(36\right)$

其中P_rob{A}表示A事件发生的概率。对F(γ)求导就能得到多普勒频移的概率密度函数：

$p\left(\mathrm{\γ}\right)=\frac{\∂F_{\mathrm{\γ}}\left(\mathrm{\γ}\right)}{\∂\mathrm{\γ}}=\frac{1}{\mathrm{\π}}{\∫}_0^{\mathrm{\π}/2}\frac{p\left({\mathrm{\β}}_m\right)}{\sqrt{{\cos }^2{\mathrm{\β}}_m-{\mathrm{\γ}}^2}}{\mathrm{d\β}}_m---\left(37\right)$

由于|γ|≤cosβ_m，则上式的上限被限制成β_m≤arccos|γ|，

$p\left(\mathrm{\γ}\right)=\frac{\∂F_{\mathrm{\γ}}\left(\mathrm{\γ}\right)}{\∂\mathrm{\γ}}=\frac{1}{\mathrm{\π}}{\∫}_0^{\arccos \left|\mathrm{\γ}\right|}\frac{p\left({\mathrm{\β}}_m\right)}{\sqrt{{\cos }^2{\mathrm{\β}}_m-{\mathrm{\γ}}^2}}{\mathrm{d\β}}_m---\left(38\right)$

与2-D空间模型比较，本发明在移动台(MS:mobile station)端和基站(BS:basestation)端分别推导出了方位角和俯仰角的AOA边缘密度函数，拓展了波达信号在垂直面内的参数估计。在此基础上，可以利用此信道的多径分量能够有效地计算出各阵元之间的SFC，为以后研究多天线MIMO阵列系统的信道容量性能提供了有力的支撑。与S.J.Nawaz的3D椭圆模型相比较，本发明给出了波达信号TOA联合密度函数的解析表达式，更好地阐明和刻画出了TOA的变化规律。在3D空间信道模型中，移动台MS的移动特性会使信号产生多普勒频移，利用方位角和俯仰角的概率密度函数之间的相互独立性，推导出多普勒频移的表达式。本模型能更加准确灵活方便地估计宏小区和微小区等移动通信环境。

实验结果与分析：

基于本发明建模方法推导出的AOA、TOA和多普勒频移表达式，我们进行仿真实验并具体分析如下：

1.AOA结果分析

图6中显示出了参数a₁/b₁对MS端EA边缘密度函数的影响。从图中可以看出，随着a₁/b₁不断增大，AOA俯仰角的概率密度主要集中在β_m＝0°处，而在大角度处β_m＞50°其概率分布非常小。这是因为随着a₁/b₁不断增打，也就是椭圆体短轴的不断变小，散射体在EA小角度的数量比在大角度的数量要多很多，这就导致了其信号的反射概率比较大。特别地，当a₁/b₁＝1时，EA的边缘密度函数曲线正好是一条余弦函数曲线。图7显示了在空心椭圆球体取不同形状下，MS端AOA俯仰角边缘密度函数的变化规律。当a₂＜b₂，与图6相比较，随着短轴的减小，图7显示的密度函数曲线变化趋势波动比较小，但是还是一条单调变化的曲线。而当a₂＞b₂，此时的曲线出现了极大值，极大值出现的地方随着长轴a₂的增大而向EA大角度靠拢。特别地，当a2增大到100m极限长度时，其概率密度函数分布在β_m＝0°的取值为零，其最大概率在β_m＝40°左右。

在图8中显示了参数a2/a1对AA边缘密度函数的影响。从图中可以发现，AOA方位角的边缘密度概率分布呈现左右对称的特征。图8清晰的显示了随着散射体空心率(a2/a1的比值)的变化，决定了AA概率分布的变化规律。随着空心率趋近于零，本文提出的模型就渐变成了Janaswamy模型。图中显示，在φ_b＝0°处AA的概率密度存在极大值，而在|φ_b|＝φ_m，其概率密度为零值，与Janaswamy模型的结果是一致的。这是因为在靠近BS端的散射体数量比较多从而导致其反射概率比较大，而在MS两侧散射体反射概率就比较小。

图9是在散射体均匀分布的空心3D空间模型下，垂直高度H对EA边缘密度函数的影响。从图中可以发现，随着垂直高度的不断变大，其概率密度曲线随着β_b的范围向右移动。从式(22)和式(23)可以看出，天线的垂直高度决定着BS端EA的范围。和Janaswamy模型相比较，显示在case1的情况下其概率密度值都是逐渐增大的。而由于本模型存在空心散射体区域，在case3情况下其概率密度值是非常缓慢地增加，基本是接近于一条直线。图10给出了空心率对BS端EA边缘密度函数的影响。随着空心率的变大，也就是说当散射体空心区域的范围逐渐变大时，EA的边缘密度函数值是在逐渐降低的，这符合信道的统计特征。

2.TOA结果分析

Janaswamy没有研究3D模型的TOA概率密度函数，本文基于改进的三维域空间模型，分别给出了MS端EA和AA的TOA联合概率密度函数。如图11所示在散射体均匀分布下，MS端AA的时延特性TOA联合概率密度。从图中可以发现，MS端的TOA最大概率主要集中在最小时延τ_min和水平面角φ_m＝0°的区域。作为对比图，同时给出了Janaswamy均匀分布3D模型的TOA估计结果。由图11可以看出，随着空心率(a2/a1的比值)的不断变大，TOA联合密度函数分布主要集中在大角度和长时延处。尤其a₂→a₁时，可以直观的看出3D椭圆模型趋近于3D椭圆环模型时，其TOA联合分布函数呈现“人型”分布。如果假设a₂＝b₂→0、H→0和b₁→0，本文提出的模型就转化为Ertel.R圆模型，两者的分析结果趋于一致。

如图12所示在散射体均匀分布下，MS端EA的时延特性TOA联合概率密度。移动台MS接收信号表现集中在全方位和短时延。随着空心率的增加，受MS周围远处散射体的影响，BS所接收信号经历角度到达和时延也不断增大。但是不管空心率怎么变化，TOA联合概率密度峰值还是在点达到峰值，其中和φ_m＝0⁺。这是因为信号在直达视距LOS上的散射体反射概率是最大的，数值结果符合信道定性分析。

3.DS结果分析

图13给出了参数a₁/b₁(ε₁)对多普勒频移的概率密度函数的影响。图14显示了不同形状的3D散射体空间对多普勒频移概率密度的影响。在图13中，随着ε₁逐渐增大，多普勒频移的概率密度值|γ|＝1不断变大。从图14中可以发现，当b₂＜a₂＝50m时，与图13相比较，其图形是倒U型，其概率密度的最大值在|γ|＝0处。而当b₂＞a₂＝50m时，最大值在|γ|＝1处。这是因为随着移动台周围空心散射体区域的逐渐变大，EA的密度函数就有可能逐渐变小，从而导致多普勒频移密度函数的改变。当ε₁→∞，即b₁＝0时，其3D空间模型就变成了经典Clarke 2-D模型，EA的密度函数可以看作p(β_m)＝δ(β_m)，其中δ(·)是狄拉克δ函数。则经典Clarke 2-D模型多普勒频移概率密度函数为

通过上述分析可知，本发明提供的建模方法，其仿真结果符合理论和经验，能够很好地模拟宏小区和微小区等移动通信环境，此外，还能进一步研究模型的各项输入参数与重要空时信道参数之间的关系，拓展了空间统计信道模型的研究和应用。

本发明方案所公开的技术手段不仅限于上述实施方式所公开的技术手段，还包括由以上技术特征任意组合所组成的技术方案。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也视为本发明的保护范围。

Claims (2)

1.一种基于多天线MIMO 3D空心椭球的统计信道建模方法，其特征在于，包括如下步骤：

步骤一：建立多天线MIMO 3D空心椭球统计信道的模型，包括基站和移动台，以移动台MS为原点，建立三维坐标系和空心椭球模型，所述空心椭球模型内分布有散射体，以移动台MS为中心的椭圆球体内部空心空间内不存在散射体，模型中r_b为基站到某个散射体的距离，φ_b为基站到某个散射体的方位角，β_b为基站到某个散射体的俯仰角，r_m为移动台到某个散射体的距离，φ_m为移动台到某个散射体的方位角，β_m为移动台到某个散射体的俯仰角，D_Los为基站和移动台之间的直达距离，D为基站和移动台之间的水平距离，H为基站离地面的垂直高度，a₁为散射体3D椭圆球体空间的长轴，a₂为散射体3D椭圆球体空间的短轴，ρ_b为r_b在水平面内的投影，b₁为散射体3D椭圆球体内部空心空间的长轴，b₂为散射体3D椭圆球体内部空心空间的短轴；

步骤二：计算到达角度AOA的概率密度函数：

步骤二-1：计算3D空心椭球移动台端方位角和俯仰角的AOA概率密度函数

(1)定义散射体分布函数：

$f(x_m,y_m,z_m)=\left\{\begin{array}{cc}\frac{1}{V}& (x_m,y_m,x_m)\∈I\\ 0& otherwise\end{array}\right.$

式中，I是空心椭圆球体的散射体分布空间，空间体积V为：

$V=\frac{2}{3}\mathrm{\π}\left(a_1^2{b}_1-a_2^2{b}_2\right)$

(2)通过雅可比式将坐标(x_m，y_m，z_m)转换为(r_m，φ_m，β_m)：

$\begin{array}{c}p(r_m,{\mathrm{\β}}_m,{\mathrm{\φ}}_m)=\frac{f\left(x_m,y_m,z_m\right)}{\left|J\left(x_m,y_m,z_m\right)\right|}{|}_{\begin{array}{c}{x}_m=r_m{\mathrm{cos\β}}_m{\mathrm{sin\φ}}_m\\ y_m=r_m{\mathrm{cos\β}}_m{\mathrm{cos\φ}}_m\\ z_m=r_m{\mathrm{sin\β}}_m\end{array}}\\ =\frac{r_m^2{\mathrm{cos\β}}_m}{V}\end{array}$

(3)由上式对r_m进行积分得到p(β_m，φ_m)概率密度函数：

$\begin{array}{c}p\left({\mathrm{\β}}_m,{\mathrm{\φ}}_m\right)={\∫}_{r_{\min }}^{r_{\max }}p(r_m,{\mathrm{\β}}_m,{\mathrm{\φ}}_m){\mathrm{dr}}_m\\ =\frac{{\mathrm{cos\β}}_m}{2\mathrm{\π}\left(a_1^2{b}_1-a_2^2{b}_2\right)}\[\frac{a_1^3{b}_1^3}{{\left(b_1^2{\cos }^2{\mathrm{\β}}_m+a_1^2{\sin }^2{\mathrm{\β}}_m\right)}^{3/2}}\\ -\frac{a_2^3{b}_2^3}{{\left(b_2^2{\cos }^2{\mathrm{\β}}_m+a_2^2{\sin }^2{\mathrm{\β}}_m\right)}^{3/2}}\]\\ 0\≤{\mathrm{\φ}}_m\≤2\mathrm{\π},0\≤{\mathrm{\β}}_m\≤/2\end{array}$

其中，

(4)由p(β_m，φ_m)概率密度函数对β_m直接积分得到方位角的边缘密度函数p(φ_m)：

$p\left({\mathrm{\φ}}_m\right)={\∫}_0^{\mathrm{\π}/2}p({\mathrm{\β}}_m,{\mathrm{\φ}}_m){\mathrm{d\β}}_m=\frac{1}{2\mathrm{\π}}$

(5)由p(β_m，φ_m)概率密度函数对φ_m积分得到俯仰角的边缘密度函数p(β_m)：

$\begin{array}{c}p\left({\mathrm{\β}}_m\right)={\∫}_0^{2\mathrm{\π}}p({\mathrm{\β}}_m,{\mathrm{\φ}}_m){\mathrm{d\φ}}_m\\ =\frac{{\mathrm{cos\β}}_m}{\left(a_1^2{b}_1-a_2^2{b}_2\right)}\[\frac{a_1^3{b}_1^3}{{\left(b_1^2{\cos }^2{\mathrm{\β}}_m+a_1^2{\sin }^2{\mathrm{\β}}_m\right)}^{3/2}}\\ -\frac{a_2^3{b}_2^3}{{\left(b_2^2{\cos }^2{\mathrm{\β}}_m+a_2^2{\sin }^2{\mathrm{\β}}_m\right)}^{3/2}}\]\end{array}$

步骤二-2：计算基站端方位角和俯仰角的AOA概率密度函数

(1)计算基站端的方位角和俯仰角的联合分布函数p(β_b，φ_b)：

$p({\mathrm{\β}}_b,{\mathrm{\φ}}_b)=\frac{({\mathrm{\ρ}}_{b_2}^3-{\mathrm{\ρ}}_{b_1}^3)}{2{\mathrm{\πcos}}^2{\mathrm{\β}}_b(a_1^2{b}_1-a_2^2{b}_2)}$

式中，ρ_b1、ρ_b2分别是和在水平面内的投影，而和是在BS端一定角度(β_b，φ_b)情况下与散射体空间相交的距离，其中，

(2)计算基站方位角的边缘密度函数p(φ_b)：

$p\left({\mathrm{\φ}}_b\right)=\∫{\∫}_{S_{\mathrm{\φ}}}\frac{1}{V}{\mathrm{\ρ}}_b{\mathrm{d\ρ}}_b{\mathrm{dz}}_b=\frac{D{\mathrm{cos\φ}}_b}{V}{A}_{\mathrm{\φ}}$

(3)根据散射体区域投影的面积A_φ，计算AOA方位角的边缘密度函数：

$\begin{array}{c}p\left({\mathrm{\φ}}_b\right)=\\ \left\{\begin{array}{c}\frac{3{\mathrm{cos\φ}}_b\[a_1{b}_1\left(a_1^2-D^2{\sin }^2{\mathrm{\φ}}_b\right)-a_2{b}_2\left(a_2^2-D^2{\sin }^2{\mathrm{\φ}}_b\right)\]}{4a_1{a}_2\left(a_1^2{b}_1-a_2^2{b}_2\right)}\\ -{\mathrm{\φ}}_m^{\′}\≤{\mathrm{\φ}}_b\≤{\mathrm{\φ}}_m^{\′}\\ \frac{3b_1{\mathrm{cos\φ}}_b\left(a_1^2-D^2{\sin }^2{\mathrm{\φ}}_b\right)}{4a_1\left(a_1^2{b}_1-a_2^2{b}_2\right)}\\ {\mathrm{\φ}}_m^{\′}\≤{\mathrm{\φ}}_b\≤{\mathrm{\φ}}_{m}and\\ -{\mathrm{\φ}}_m\≤{\mathrm{\φ}}_b\≤-{\mathrm{\φ}}_m^{\′}\\ \begin{array}{cc}0& otherwise\end{array}\end{array}\right.\end{array}$

sinφ_m＝(a₁/D)和sinφ′_m＝(a₂/D)；

(4)计算垂直面内边缘密度函数：

$p\left({\mathrm{\β}}_b\right)=\left\{\begin{array}{c}\frac{1}{2{\mathrm{\πcos}}^2{\mathrm{\β}}_b\left(a{1}^{2}b1-a{2}^{2}b2\right)}\underset{-{\mathrm{\φ}}_{\max }}{\overset{{\mathrm{\φ}}_{\max }}{\∫}}\left({\mathrm{\ρ}}_{b_2}^3-{\mathrm{\ρ}}_{b_1}^3\right){\mathrm{d\φ}}_b\\ {\mathrm{\β}}_{\min }\≤{\mathrm{\β}}_b\≤{\mathrm{\β}}_2,\\ {\mathrm{\β}}_4\≤{\mathrm{\β}}_b\≤{\mathrm{\β}}_{\max }\\ \frac{1}{2{\mathrm{\πcos}}^2{\mathrm{\β}}_b\left(a{1}^{2}b1-a{2}^{2}b2\right)}\[\underset{-{\mathrm{\φ}}_{\max }}{\overset{{\mathrm{\φ}}_{\max }}{\∫}}\left({\mathrm{\ρ}}_{b_2}^3-{\mathrm{\ρ}}_{b_1}^3\right){\mathrm{d\φ}}_b-\\ \underset{-{\mathrm{\φ}}_{\max }^{\′}}{\overset{{\mathrm{\φ}}_{\max }^{\′}}{\∫}}\left({\mathrm{\ρ}}_{b_2}^{\′3}-{\mathrm{\ρ}}_{b_1}^{\′3}\right){\mathrm{d\φ}}_b\]\\ {\mathrm{\β}}_2\≤{\mathrm{\β}}_b\≤{\mathrm{\β}}_4\end{array}\right.$

其中，其中ρ_b1、ρ_b2、ρ′_b1和ρ′_b2分别是β_b在某一取值范围内相对应的值，β_min≤β_b≤β_max，且

${\mathrm{\β}}_{\min }=arctan\left(\frac{HD-\sqrt{H^2{a}_1^2+b_1^2(D^2-a_1^2)}}{D^2-a_1^2}\right)$

${\mathrm{\β}}_{max}=arctan\left(\frac{H}{D-a_1}\right)$

${\mathrm{\β}}_1=arctan\left(\frac{H}{D+a_1}\right)$

${\mathrm{\β}}_2=arctan\left(\frac{HD-\sqrt{H^2{a}_2^2+b_2^2(D^2-a_2^2)}}{D^2-a_2^2}\right)$

${\mathrm{\β}}_3=arctan\left(\frac{H}{D+a_2}\right)$

${\mathrm{\β}}_4=arctan\left(\frac{H}{D-a_2}\right);$

步骤三：计算到达时间TOA的联合概率密度函数

(1)计算来波信号的AOA/TOA联合密度函数：

$p(\mathrm{\τ},{\mathrm{\β}}_m,{\mathrm{\φ}}_m)=\frac{p(r_m,{\mathrm{\β}}_m,{\mathrm{\φ}}_m)}{\left|J(r_m,{\mathrm{\β}}_m,{\mathrm{\φ}}_m)\right|}$

其中，

其中J(r_m，β_m，φ_m)为雅可比转换式，且

$\begin{array}{c}J\left(r_m,{\mathrm{\β}}_m,{\mathrm{\φ}}_m\right)=|\frac{\∂r_m}{\∂\mathrm{\τ}}{|}^{-1}\\ =\frac{2{\left(D{\mathrm{cos\β}}_m{\mathrm{cos\φ}}_m-c\mathrm{\τ}+H{\mathrm{sin\β}}_m\right)}^2}{c\left(S_{Los}^2-c^2{\mathrm{\τ}}^2-2c\mathrm{\τ}\left(D{\mathrm{cos\β}}_m{\mathrm{cos\φ}}_m+H{\mathrm{sin\β}}_m\right)\right)}\end{array}$

(2)通过以上式子，计算AOA/TOA的联合密度函数：

$\begin{array}{c}p(\mathrm{\τ},{\mathrm{\β}}_m,{\mathrm{\φ}}_m)\\ =\frac{(c^2{\mathrm{\τ}}^2+D_{Los}^2-2c\mathrm{\τ}(D{\mathrm{cos\β}}_m{\mathrm{cos\φ}}_m+H{\mathrm{sin\β}}_m\left)\right)}{8V{(D{\mathrm{cos\β}}_m{\mathrm{cos\φ}}_m-c\mathrm{\τ}+H{\mathrm{sin\β}}_m)}^4}\\ \×\[c{\left(D_{Los}^2-c^2{\mathrm{\τ}}^2\right)}^2{\mathrm{cos\β}}_m\]\end{array}$

其中到达时延c为光速；

对上式进行β_m积分计算移动台MS在方位角的TOA联合分布函数：

$p\left(\mathrm{\τ},{\mathrm{\φ}}_m\right)={\∫}_0^{\mathrm{\π}/2}p(\mathrm{\τ},{\mathrm{\β}}_m,{\mathrm{\φ}}_m){\mathrm{d\β}}_m,$

其中，p(τ，φ_m)的闭式表达式如下：

$\begin{array}{c}p\left(\mathrm{\τ},{\mathrm{\φ}}_m\right)=\\ \frac{1}{6{\left(k_1^2+k_2^2-k_3^2\right)}^{7/2}}\{\frac{-1}{{\left(k_1-k_3\right)}^3}\[k_2\sqrt{k_1^2+k_2^2-k_3^2}\[6k_1^5{k}_3+14k_3^6+k_1^4\left(22k_3^2-5k_4\right)-19k_3^4{k}_4+k_1^3{k}_3\left(12k_2^2-50k_3^2+k_4\right)-\\ 2k_2^4\left(k_3^2+k_4\right)+6k_2^2{k}_3^2\left(3k_3^2+k_4\right)+k_1^2\left(24k_3^4+k_2^2\left(20k_3^2-7k_4\right)-6k_3^2{k}_4\right)+k_1{k}_3\left(6k_2^4-16k_3+29k_3^2{k}_4+k_2^2\left(-50k_3^2+k_4\right)\right)\]+\\ 6k_1{\left(k_1-k_3\right)}^3\[2k_3^2\left(k_3^2-2k_4\right)+k_1^2\left(8k_3^2-k_4\right)+k_2^2\left(8k_3^2-k_4\right)\]\arctan h\left(\frac{k_2}{k_1^2+k_2^2-k_3^2}\right)\]-\\ \frac{1}{{\left(k_2-k_3\right)}^2}\[\sqrt{k_1^2+k_2^2-k_3^2}\[k_1^4(6k_2{k}_3-4k_2^2-k_4))+\\ \left(k_2-2k_2\right){\left(k_2-k_3\right)}^2\left(6k_2^2{k}_3-4k_2{k}_3^2+2k_2{k}_3-3k_3{k}_4\right)+k_1^2\left(k_2-k_3\right)\left(12k_2^2{k}_3+26k_3^3-10k_3{k}_4+k_2\left(-20k_3^2+k^4\right)\right)\]+\\ 6k_1{\left(k_2-k_3\right)}^2\[2k_3^2\left(k_3^2-2k_4\right)+k_1^2\left(8k_3^2-k_4\right)+k_2^2\left(8k_3^2-k_4\right)\]\arctan h\left(\frac{k_1-k_2+k_3}{k_1^2+k_2^2-k_3^2}\right)\]\}\end{array}$

其中

k₁＝D cosφ_m

k₂＝H

k₃＝cτ

$k_4=c^2{\mathrm{\τ}}^2+D_{Los}^2$

对AOA/TOA的联合密度函数进行φ_m积分计算移动台MS在俯仰角的TOA联合分布函数：

$\begin{array}{c}p\left(\mathrm{\τ},{\mathrm{\β}}_m\right)={\∫}_0^{2\mathrm{\π}}p\left(\mathrm{\τ},{\mathrm{\β}}_m,{\mathrm{\φ}}_m\right){\mathrm{d\φ}}_m\\ =\sqrt{\frac{-k_1+k_2}{k_1+k_2}}=\frac{(3k_1^2{k}_2{k}_3+2k_2^3{k}_3+k_1^3{k}_4+4k_1{k}_2^2{k}_4)\mathrm{\π}}{{(k_1-k_2)}^4{(k_1+k_2)}^3}\end{array}$

其中

k₁＝D cosβ_m

k₂＝H sinβ_m-cτ

$k_3=c^2{\mathrm{\τ}}^2+D_{\mathrm{Los}}^2-2\mathrm{c\τ}H\sin {\mathrm{\β}}_m$

k₄＝2cτD cosβ_m；

步骤四：计算多普勒频移DS的概率密度函数

计算多普勒频移的累积分布函数：

$\begin{array}{c}F\left(\mathrm{\γ}\right)=P_{rob}\left\{{\mathrm{cos\β}}_m{\mathrm{cos\θ}}_m\right\}\\ =P_{rob}\left\{{\mathrm{cos\β}}_m{\mathrm{cos\θ}}_m\right\},\left(p\left({\mathrm{\φ}}_m\right)=p\left({\mathrm{\θ}}_m\right)-\frac{1}{2\mathrm{\π}}\right)\\ =2{\∫}_0^{\mathrm{\π}/2}p\left({\mathrm{\β}}_m\right)\[{\∫}_{\arccos \[\mathrm{\γ}/{\mathrm{cos\β}}_m\]}^{\mathrm{\π}}p\left({\mathrm{\φ}}_m\right){\mathrm{d\φ}}_m\]{\mathrm{d\β}}_m\\ =\frac{1}{\mathrm{\π}}{\∫}_0^{\mathrm{\π}/2}p\left({\mathrm{\β}}_m\right)\[{\∫}_{\arccos \[\mathrm{\γ}/{\mathrm{cos\β}}_m\]}^{\mathrm{\π}}{\mathrm{d\φ}}_m\]{\mathrm{d\β}}_m\end{array}$

其中，P_rob{A}表示A事件发生的概率，γ≡f_DS/f_m，所述其中v为移动台的移动速度，f_c为信号的载波频率，f_m＝vf_c/c是最大多普勒频移，θ_m和β_m分别为移动台MS到达角度的方位角和俯仰角，

对F(γ)求导计算多普勒频移的概率密度函数：

$p\left(\mathrm{\γ}\right)=\frac{\∂F_{\mathrm{\γ}}\left(\mathrm{\γ}\right)}{\∂\mathrm{\γ}}=\frac{1}{\mathrm{\π}}{\∫}_0^{\arccos \left|\mathrm{\γ}\right|}\frac{p\left({\mathrm{\β}}_m\right)}{\sqrt{{\cos }^2{\mathrm{\β}}_m-{\mathrm{\γ}}^2}}{\mathrm{d\β}}_m$

其中，|γ|≤cosβ_m。

2.根据权利要求1所述的基于多天线MIMO 3D空心椭球的统计信道建模方法，其特征在于：所述的散射体均匀分布在空心椭球体区域内。