CN103747456A - Modeling method based on three-dimensional spatial domain multi-antenna MIMO (Multiple input and multiple output) statistical channel - Google Patents

Abstract

The invention discloses a modeling method specific to an indoor 3D (Three-Dimensional) spatial statistical channel. The modeling method can be applied to a distributed multiple input and multiple output (MIMO) system. When a base station (BS) is designed with a directional antenna, important space-time channel parameters such as Doppler Power Spectrum (DS), space-time correlation and channel capacity are concluded by taking a Poda signal at an MS (Mobile Station) end as a basis; space-time model parameters and the mechanism relationship among included angles of directional antenna main lobes are resolved by using channel space-time feature parameters; meanwhile, the system transmission performance of an MIMO multi-antenna linear array and a circular array under a 3D spatial statistical channel model is researched, so that accurate, flexible diverse channel models are provided for an MIMO multi-antenna space-time processing algorithm and a simulation wireless communication system, and the research of the spatial statistical channel model is expanded.

Description

Modeling method based on three dimensions territory many antennas MIMO statistical channel

Technical field

The invention belongs to wireless communication technology field, especially relate to a kind of modeling method of the many antennas mimo channel based on three dimensions territory.

Background technology

Mobile communication is a part with fastest developing speed in communications industry, and the essence of mobile communication is to utilize wireless channel to carry out effective transmission of information.Performance of mobile communication system is mainly subject to the restriction of characteristics of radio channels.Multiple-input and multiple-output (MIMO:Multiple Input Multiple Output) technology can greatly increase power system capacity and improve the advantage of radio link quality and received increasing attention and concern with it, and also therefore by many people, being predicted as the 4th generation of moving communicating field is one of 4G core technology.MIMO technology is a very brand-new technology, and MIMO technology has many outstanding performances, such as the higher availability of frequency spectrum, and lower transmission error rates etc.Along with the development of MIMO technology, to being applied to the research of mimo system antenna, become very urgent, also there is theory significance and practical value simultaneously.

Propagation path between base station (BS:base station) and travelling carriage (MS:mobile station) is generally distributed with complicated landform, and its channel is on-fixed and unpredictalbe often.Multipath effect is the multipath fading in mobile telecommunication channel, is one of main contents of wireless channel research, therefore sets up one accurately and the important step that effective channel model is structure mobile communication system.Ertel and Petrus have proposed scattering object spatial distribution circle model (GBSBM:geometrically based single bounce model) and elliptical model (EBSBM:Ellipse based single bounce model).Simulation result shows that GBSBM model can estimate important parameter under macrocell (Macrocell) mobile communication environment, and EBSBM model can be estimated important channel parameter under the mobile communication environment of Microcell (Microcell).Because the estimated result of GBSBM and EBSBM model is not accurate enough, Jiang provides based on scattering object Rayleigh and index (Exponential) and divides the round model planting, and Janaswamy and Zhou have proposed respectively round model and hollow ring model (HSDM:hollow-disc scatter density model) that scattering object Gauss (Gaussian) distributes.Research about channel model is all based in 2D plane above, derive and calculate also just rest on empty in horizontal plane time channel parameter, can not well reflect objective reality, for the complexity of space angle, can not have one to describe accurately.

On the basis of studying at 2D channel model based on forefathers, Janaswamy etc. have proposed 3d space statistical channel model, the angle of base station BS and mobile station MS end is refined as to the angle on horizontal plane and perpendicular, particularly Nawaz has proposed outdoor directional antenna and has divided the 3D planting model, has derived AOA and TOA probability density function on horizontal plane and vertical plane.But above research is mostly the research about outdoor macro cell 3D model, and for the exploration under the environment of indoor Microcell, still has certain blank.

Summary of the invention

For still lacking indoor Microcell this defect of 3D channel model in prior art, the invention discloses a kind of for indoor 3d space statistical channel modeling method, during for MIMO multi-antenna space, Processing Algorithm and emulation wireless communication system provide accurate and flexible and changeable channel model, have expanded the research of spatial statistics channel model.

In the mobile communication environment of the indoor Microcell of 3D, because complex environment can affect the variation of spatial channel parameters, and signal propagation can be subject to stopping of wall, ground, building and other object, increased the weight of multipath fading effect, thereby make system be subject to great uncertainty, affect us for the research and analysis of wireless channel.Therefore, the present invention is when base station BS is designed with oriented antenna, be directed to the complexity of 3d space angle, the angle of BS and MS end is refined as to the angle on horizontal plane and perpendicular, channel parameter during important empty comprehensively and under the indoor Microcell of the description 3D environment of accurate and flexible, is embodied in following technical scheme:

Step 1, the base station BS of take is set up three-dimensional system of coordinate as initial point, and wherein the main lobe angle of the directional antenna based on BS end is 2 α, and the angle of BS on horizontal plane is φ _b, the angle on perpendicular is β _b; The angle of MS on horizontal plane and perpendicular is respectively φ _mand β _m, the scattering object in scattering region is respectively r to the distance of BS and MS _band r _m, the distance between BS and MS is D, the major axis of 3D channel model and minor axis length are respectively a and b, and D<a, b≤a, all scattering objects are evenly distributed in scattering region, and the volume of scattering region is V=2a ²b α/3;

Step 2, is divided into P by scattering region ₁and P ₂two parts, wherein, P ₂refer to that the scattering object of electromagnetic signal on the vertical scattering boundary of grey reflexes to the propagation path region on MS, and P ₁scattering region I _regionremaining part;

Scattering region P ₁and P ₂be expressed as

${\mathrm{<figure-callout\; id="1"\; label="P"\; filenames="140120170431.png,140120170435.png"\; state="\{\{state\}\}">P</figure-callout>}}_1\&RightArrow;\left\{\begin{array}{c}0\&le;{\mathrm{\&beta;}}_M\&le;{\mathrm{\&beta;}}_t\\ \mathrm{or}\\ {\mathrm{\&phi;}}_{\mathrm{<figure-callout\; id="1"\; label="t"\; filenames="140120170431.png,140120170435.png"\; state="\{\{state\}\}">t</figure-callout>}1}\&le;\left|{\mathrm{\&phi;}}_m\right|\&le;{\mathrm{\&phi;}}_{t2}\end{array}\right.{\mathrm{<figure-callout\; id="2"\; label="P"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">P</figure-callout>}}_2\&RightArrow;\left\{\begin{array}{c}{\mathrm{\&beta;}}_t\&le;{\mathrm{\&beta;}}_M\&le;\mathrm{\&pi;}/2\\ \mathrm{or}\\ {\mathrm{\&phi;}}_{\mathrm{<figure-callout\; id="2"\; label="t"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">t</figure-callout>}2}\&le;\left|{\mathrm{\&phi;}}_m\right|\&le;2\mathrm{\&pi;}-{\mathrm{\&phi;}}_{\mathrm{<figure-callout\; id="2"\; label="t"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">t</figure-callout>}2}\end{array}\right.$

Wherein, β _tfor the vertical angle of plane P MQ,

${\mathrm{\&beta;}}_t=\left\{\begin{array}{cc}{\cot }^{-1}\left(\frac{\mathrm{\&alpha;}D\csc (\mathrm{\&alpha;}+{\mathrm{\&phi;}}_m)\sin \mathrm{\&alpha;}}{b\sqrt{a^2-{\mathrm{<figure-callout\; id="2"\; label="D"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">D</figure-callout>}}^2{\csc }^2(\mathrm{\&alpha;}+{\mathrm{\&phi;}}_m){\sin }^2{\mathrm{\&phi;}}_m}}\right),& {\mathrm{\&phi;}}_1\&le;\left|{\mathrm{\&phi;}}_m\right|\&le;{\mathrm{\&phi;}}_2\\ 0,& \mathrm{otherwise}\end{array}\right.,$

Wherein, angle φ ₁and φ ₂for β _m=0 o'clock, the angle of MS end on xoy horizontal plane,

φ _t1for the start angle of the section PMQ on vertical direction on level angle,

φ _t1=0，0≤β _M≤π／2，

φ _t2for the termination point of the section PMQ on vertical direction on level angle,

${\mathrm{\&phi;}}_{\mathrm{<figure-callout\; id="2"\; label="t"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">t</figure-callout>}2}=\left\{\begin{array}{cc}\arccos \left\{\frac{{\mathrm{<figure-callout\; id="2"\; label="PM"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">PM</figure-callout>}}^2+{\mathrm{QM}}^2-{\mathrm{<figure-callout\; id="2"\; label="PQ"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">PQ</figure-callout>}}^2}{2\mathrm{PM}\&times;\mathrm{QM}}\right\},& 0\&le;{\mathrm{\&beta;}}_M\&le;\arctan \left(\frac{b}{D\sin \mathrm{\&alpha;}}\right)\\ 0,& \arctan \left(\frac{b}{D\sin \mathrm{\&alpha;}}\right)\&le;{\mathrm{\&beta;}}_M\&le;\mathrm{\&pi;}/2\end{array}\right.,$

Wherein, P and Q are section on perpendicular and the intersection point of scattering boundary, and the dihedral angle of PMQ plane and xoy plane is β _m, and

$\begin{array}{c}\mathrm{PQ}=\sqrt{a^2-\frac{a^2}{b^2}{\mathrm{<figure-callout\; id="2"\; label="D"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">D</figure-callout>}}^2{\mathrm{<figure-callout\; id="2"\; label="sin"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">sin</figure-callout>}}^2{\mathrm{\&alpha;}\tan }^2{\mathrm{\&beta;}}_M}\\ \mathrm{QM}=\sqrt{{\mathrm{<figure-callout\; id="2"\; label="d"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">d</figure-callout>}}^2+{\mathrm{<figure-callout\; id="2"\; label="D"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">D</figure-callout>}}^2{\mathrm{<figure-callout\; id="2"\; label="sin"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">sin</figure-callout>}}^2\mathrm{\&alpha;}{\tan }^2{\mathrm{\&beta;}}_M}\\ \mathrm{PM}=\sqrt{{\mathrm{<figure-callout\; id="2"\; label="D"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">D</figure-callout>}}^2+{\mathrm{<figure-callout\; id="2"\; label="D"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">D</figure-callout>}}^2{\mathrm{<figure-callout\; id="2"\; label="sin"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">sin</figure-callout>}}^2\mathrm{\&alpha;}{\tan }^2{\mathrm{\&beta;}}_M+{\mathrm{PQ}}^2-2D\&times;\mathrm{PQ}\cos \mathrm{\&alpha;}}\end{array},$

Step 3: by the distribution function f (x of scattering object _m, y _m, z _m)=1/V obtains by Jacobi's change type:

$p(r_m,{\mathrm{\&phi;}}_m,{\mathrm{\&beta;}}_m)=\frac{f(x_m,y_m,z_m)}{\left|J(x_m,y_m,z_m)\right|}{|}_{\begin{array}{c}{x}_m=r_m\cos {\mathrm{\&beta;}}_m\cos {\mathrm{\&phi;}}_m\\ y_m=r_m\cos {\mathrm{\&beta;}}_m\sin {\mathrm{\&phi;}}_m\\ z_m=r_m\sin {\mathrm{\&beta;}}_m\end{array}}=\frac{{\mathrm{<figure-callout\; id="2"\; label="r\; m"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">r</figure-callout>}}_m^{}\cos {\mathrm{\&beta;}}_m}{V}$

Wherein, r _mfor the distance of MS to scattering boundary, by following formula, calculate:

$r_m({\mathrm{\&phi;}}_m,{\mathrm{\&beta;}}_m)=\left\{\begin{array}{cc}D\sin \mathrm{\&alpha;}\csc (\mathrm{\&alpha;}+{\mathrm{\&phi;}}_m)\sec {\mathrm{\&beta;}}_m,& {\mathrm{<figure-callout\; id="1"\; label="P"\; filenames="140120170431.png,140120170435.png"\; state="\{\{state\}\}">P</figure-callout>}}_1\\ \frac{1}{{\mathrm{<figure-callout\; id="2"\; label="b"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">b</figure-callout>}}^2{\mathrm{<figure-callout\; id="2"\; label="cos"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">cos</figure-callout>}}^2{\mathrm{\&beta;}}_m+{\mathrm{<figure-callout\; id="2"\; label="a"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">a</figure-callout>}}^2{\mathrm{<figure-callout\; id="2"\; label="sin"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">sin</figure-callout>}}^2{\mathrm{\&beta;}}_m}& \\ \&times;\{{\mathrm{<figure-callout\; id="2"\; label="Db"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">Db</figure-callout>}}^2\cos {\mathrm{\&beta;}}_m\cos {\mathrm{\&phi;}}_m+\sqrt{{\left({\mathrm{<figure-callout\; id="2"\; label="Db"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">Db</figure-callout>}}^2\cos {\mathrm{\&beta;}}_m\cos {\mathrm{\&phi;}}_m\right)}^2}& \\ \stackrel{\&OverBar;}{-(b^2{\mathrm{<figure-callout\; id="2"\; label="cos"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">cos</figure-callout>}}^2{\mathrm{\&beta;}}_m+a^2{\mathrm{<figure-callout\; id="2"\; label="sin"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">sin</figure-callout>}}^2{\mathrm{\&beta;}}_m)(b^2{D}^2-a^2{b}^2)}\},& {\mathrm{<figure-callout\; id="2"\; label="P"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">P</figure-callout>}}_2\end{array}\right.$

By r _mjacobi's change type in this step of substitution, obtains the AOA joint probability density function of MS:

$p({\mathrm{\&phi;}}_m,{\mathrm{\&beta;}}_m)=\left\{\begin{array}{cc}\frac{\cos {\mathrm{\&beta;}}_m}{3V}{\left\{D\sin \mathrm{\&alpha;}\csc (\mathrm{\&alpha;}+{\mathrm{\&phi;}}_m)\sec {\mathrm{\&beta;}}_m\right\}}^3,& {\mathrm{<figure-callout\; id="1"\; label="P"\; filenames="140120170431.png,140120170435.png"\; state="\{\{state\}\}">P</figure-callout>}}_1\\ \frac{\cos {\mathrm{\&beta;}}_m}{3V}\{\frac{1}{{\mathrm{<figure-callout\; id="2"\; label="b"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">b</figure-callout>}}^2{\mathrm{<figure-callout\; id="2"\; label="cos"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">cos</figure-callout>}}^2{\mathrm{\&beta;}}_m+a^2{\mathrm{<figure-callout\; id="2"\; label="sin"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">sin</figure-callout>}}^2{\mathrm{\&beta;}}_m}& \\ \&times;({\mathrm{<figure-callout\; id="2"\; label="Db"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">Db</figure-callout>}}^2\cos {\mathrm{\&beta;}}_m\cos {\mathrm{\&phi;}}_m+\sqrt{{\left({\mathrm{<figure-callout\; id="2"\; label="Db"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">Db</figure-callout>}}^2\cos {\mathrm{\&beta;}}_m\cos {\mathrm{\&phi;}}_m\right)}^2}& \\ \stackrel{\&OverBar;}{-({\mathrm{<figure-callout\; id="2"\; label="b"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">b</figure-callout>}}^2{\mathrm{<figure-callout\; id="2"\; label="cos"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">cos</figure-callout>}}^2{\mathrm{\&beta;}}_m+a^2{\mathrm{<figure-callout\; id="2"\; label="sin"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">sin</figure-callout>}}^2{\mathrm{\&beta;}}_m(b^2{D}^2-a^2{b}^2)}){\}}^3,& {\mathrm{<figure-callout\; id="2"\; label="P"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">P</figure-callout>}}_2\end{array}\right.;$

Step 4: Jacobi's change type in step 3 is obtained by Jacobi's change type again:

$\begin{array}{c}p(\mathrm{\&tau;},{\mathrm{\&phi;}}_m,{\mathrm{\&beta;}}_m)=\frac{p(r_m,{\mathrm{\&phi;}}_m,{\mathrm{\&beta;}}_m)}{\left|J(r_m,{\mathrm{\&phi;}}_m,{\mathrm{\&beta;}}_m)\right|}{|}_{{\mathrm{<figure-callout\; id="1"\; label="P"\; filenames="140120170431.png,140120170435.png"\; state="\{\{state\}\}">P</figure-callout>}}_1,{\mathrm{<figure-callout\; id="2"\; label="P"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">P</figure-callout>}}_2}\\ =\frac{c{({\mathrm{<figure-callout\; id="2"\; label="c"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">c</figure-callout>}}^2{\mathrm{\&tau;}}^2-D^2)}^2({\mathrm{<figure-callout\; id="2"\; label="D"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">D</figure-callout>}}^2+{\mathrm{<figure-callout\; id="2"\; label="c"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">c</figure-callout>}}^2{\mathrm{\&tau;}}^2-2\mathrm{c\&tau;}D\cos {\mathrm{\&beta;}}_m\cos {\mathrm{\&phi;}}_m)\cos {\mathrm{\&beta;}}_m}{8V{(D\cos {\mathrm{\&beta;}}_m\cos {\mathrm{\&phi;}}_m-\mathrm{c\&tau;})}^4}\end{array}$

Wherein, τ is for arriving time delay, minimum value τ _min=D/c, maximum delay is

By above formula to β _mcarry out integration, obtain the delay character TOA of mobile station MS and the joint probability density function of horizontal plane angle is:

$\begin{array}{c}p(\mathrm{\&tau;},{\mathrm{\&phi;}}_m)={\&Integral;}_0^{\mathrm{\&pi;}/2}p(\mathrm{\&tau;},{\mathrm{\&phi;}}_m,{\mathrm{\&beta;}}_m)d{\mathrm{\&beta;}}_m\\ =\frac{{\mathrm{<figure-callout\; id="1"\; label="k"\; filenames="140120170431.png,140120170435.png"\; state="\{\{state\}\}">k</figure-callout>}}_1}{48\mathrm{<figure-callout\; id="5"\; label="V\; k"\; filenames="140120170435.png"\; state="\{\{state\}\}">V</figure-callout>}{k}_{}^{}{}_5^2{(k_4^2-{\mathrm{<figure-callout\; id="5"\; label="k"\; filenames="140120170435.png"\; state="\{\{state\}\}">k</figure-callout>}}_5^2)}^{7/2}}\&times;\{\sqrt{k_4^2-{\mathrm{<figure-callout\; id="5"\; label="k"\; filenames="140120170435.png"\; state="\{\{state\}\}">k</figure-callout>}}_5^2}(k_3{k}_5(2k_4^3+13k_4{\mathrm{<figure-callout\; id="5"\; label="k"\; filenames="140120170435.png"\; state="\{\{state\}\}">k</figure-callout>}}_5^2)+k_2(k_4^4-10k_4^2{\mathrm{<figure-callout\; id="5"\; label="k"\; filenames="140120170435.png"\; state="\{\{state\}\}">k</figure-callout>}}_5^2-6{\mathrm{<figure-callout\; id="5"\; label="k"\; filenames="140120170435.png"\; state="\{\{state\}\}">k</figure-callout>}}_5^4))\\ -6{\mathrm{<figure-callout\; id="5"\; label="k"\; filenames="140120170435.png"\; state="\{\{state\}\}">k</figure-callout>}}_5^2(k_3{k}_5(4k_4^2+{\mathrm{<figure-callout\; id="5"\; label="k"\; filenames="140120170435.png"\; state="\{\{state\}\}">k</figure-callout>}}_5^2)-k_2(k_4^3+4k_4{\mathrm{<figure-callout\; id="5"\; label="k"\; filenames="140120170435.png"\; state="\{\{state\}\}">k</figure-callout>}}_5^2))a\tanh \left(\frac{k_4+k_5}{\sqrt{k_4^2-{\mathrm{<figure-callout\; id="5"\; label="k"\; filenames="140120170435.png"\; state="\{\{state\}\}">k</figure-callout>}}_5^2}}\right)\}\end{array}$

To φ _mcarry out integration, obtain the delay character TOA of MS and the joint probability density function of perpendicular angle is:

$\begin{array}{c}p(\mathrm{\&tau;},{\mathrm{\&beta;}}_m)={\&Integral;}_0^{2\mathrm{\&pi;}}p(\mathrm{\&tau;},{\mathrm{\&phi;}}_b,{\mathrm{\&beta;}}_b)d{\mathrm{\&phi;}}_m\\ =\frac{k_6(k_7{\mathrm{<figure-callout\; id="8"\; label="k"\; filenames="140120170423.png,140120170427.png"\; state="\{\{state\}\}">k</figure-callout>}}_8^3-3k_2{\mathrm{<figure-callout\; id="8"\; label="k"\; filenames="140120170423.png,140120170427.png"\; state="\{\{state\}\}">k</figure-callout>}}_8^2{\mathrm{<figure-callout\; id="5"\; label="k"\; filenames="140120170435.png"\; state="\{\{state\}\}">k</figure-callout>}}_5+4k_7{k}_8{\mathrm{<figure-callout\; id="5"\; label="k"\; filenames="140120170435.png"\; state="\{\{state\}\}">k</figure-callout>}}_5^2-2k_2{\mathrm{<figure-callout\; id="5"\; label="k"\; filenames="140120170435.png"\; state="\{\{state\}\}">k</figure-callout>}}_5^3)}{8V{(k_8-k_5)}^3{({\mathrm{<figure-callout\; id="8"\; label="k"\; filenames="140120170423.png,140120170427.png"\; state="\{\{state\}\}">k</figure-callout>}}_8+k_5)}^4}\&times;\sqrt{1+\frac{2k_8}{-{\mathrm{<figure-callout\; id="8"\; label="k"\; filenames="140120170423.png,140120170427.png"\; state="\{\{state\}\}">k</figure-callout>}}_8+{\mathrm{<figure-callout\; id="5"\; label="k"\; filenames="140120170435.png"\; state="\{\{state\}\}">k</figure-callout>}}_5}}\mathrm{\&pi;}\end{array}$

Wherein,

$\left\{\begin{array}{c}{\mathrm{<figure-callout\; id="1"\; label="k"\; filenames="140120170431.png,140120170435.png"\; state="\{\{state\}\}">k</figure-callout>}}_1=c{({\mathrm{<figure-callout\; id="2"\; label="c"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">c</figure-callout>}}^2{\mathrm{\&tau;}}^2-D^2)}^2\\ {\mathrm{<figure-callout\; id="2"\; label="k"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">k</figure-callout>}}_2={\mathrm{<figure-callout\; id="2"\; label="c"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">c</figure-callout>}}^2{\mathrm{\&tau;}}^2+D^2\\ {\mathrm{<figure-callout\; id="3"\; label="k"\; filenames="140120170442.png,140120170446.png"\; state="\{\{state\}\}">k</figure-callout>}}_3=2\mathrm{c\&tau;}D\cos {\mathrm{\&phi;}}_m\\ k_4=D\cos {\mathrm{\&phi;}}_{\mathrm{<figure-callout\; id="5"\; label="m\; k"\; filenames="140120170435.png"\; state="\{\{state\}\}">m</figure-callout>}}& \\ k_{}{}_5=\mathrm{c\&tau;}\end{array}\right.\mathrm{and}\left\{\begin{array}{c}{\mathrm{<figure-callout\; id="6"\; label="k"\; filenames="140120170458.png,140120170502.png"\; state="\{\{state\}\}">k</figure-callout>}}_6=c{({\mathrm{<figure-callout\; id="2"\; label="c"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">c</figure-callout>}}^2{\mathrm{\&tau;}}^2-D^2)}^2\cos {\mathrm{\&beta;}}_m\\ k_7={\mathrm{<figure-callout\; id="2"\; label="c"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">c</figure-callout>}}^2{\mathrm{\&tau;}}^2+D^2\\ {\mathrm{<figure-callout\; id="8"\; label="k"\; filenames="140120170423.png,140120170427.png"\; state="\{\{state\}\}">k</figure-callout>}}_8=2\mathrm{c\&tau;}D\cos {\mathrm{\&beta;}}_{\mathrm{<figure-callout\; id="9"\; label="m\; k"\; filenames="140120170427.png"\; state="\{\{state\}\}">m</figure-callout>}}& \\ k_{}{}_9=D\cos {\mathrm{\&beta;}}_{\mathrm{<figure-callout\; id="10"\; label="m\; k"\; filenames="140120170423.png,140120170427.png"\; state="\{\{state\}\}">m</figure-callout>}}& \\ k_{}{}_{10}=\mathrm{c\&tau;}\end{array}\right.;$

Step 5: by classical Clarke model, the probability density function of Doppler frequency shift is as follows:

$f_{\mathrm{\&upsi;}}\left(f\right)=\frac{p({\mathrm{\&phi;}}_{\mathrm{\&upsi;}}+|{\cos }^{-1}(f/f_m\cos \mathrm{\&beta;})|,{\mathrm{\&beta;}}_m)}{f_m\sqrt{{1-(f/f_m)}^2}}+\frac{p({\mathrm{\&phi;}}_{\mathrm{\&upsi;}}-|{\cos }^{-1}(f/f_m\cos \mathrm{\&beta;})|,{\mathrm{\&beta;}}_m)}{f_m\sqrt{{1-(f/f_m)}^2}},\left|f\right|\&le;f_m,$

Wherein, f _mthe maximum doppler frequency of mobile station MS, φ _υbe the moving direction of MS, υ is translational speed, β _mangle for plane P MQ and horizontal plane xoy;

Step 6: calculate the Spatial fading correlation function between array element m and n by following formula:

Step 7: calculate channel capacity by following formula:

$C={\log }_2\det (I_{\mathrm{Nr}}+\frac{P}{N_t{\mathrm{\&sigma;}}^2}{\mathrm{HH}}^H)$

Shi Zhong C unit is (bits/s/Hz), I _nrfor N _rdimension unit matrix, P/ σ ²for channel signal to noise ratio snr, N _tfor transmitting terminal antenna amount and N _rfor receiving terminal antenna amount.

Further, for the mimo channel of space correlation, channel matrix H is expressed as

$H={\mathrm{<figure-callout\; id="1"\; label="R\; r"\; filenames="140120170431.png,140120170435.png"\; state="\{\{state\}\}">R</figure-callout>}}_r^{}{H}_w{\left({\mathrm{<figure-callout\; id="1"\; label="R\; t"\; filenames="140120170431.png,140120170435.png"\; state="\{\{state\}\}">R</figure-callout>}}_t^{}\right)}^T$

Wherein, R _rfor correlation matrix between the array element of receiving terminal, R _tfor correlation matrix between transmitting terminal array element, H _wfor the multiple gaussian random matrix with distributing, subscript ^tand ^htransposition and the conjugate transpose of difference representing matrix.

Further, in described step 6 by following two kinds of modes, obtain:

Step 6-1, when MIMO array is linear ULA, the incoming signal steric direction vector that receives MIMO ULA is

Wherein,

the angle of incoming signal on horizontal plane; θ is the angle of incoming signal on perpendicular, k _w=2 π/λ, d is antenna spacing, and λ is incoming signal wavelength, and L is receiving terminal number of antennas, [] ^trepresenting matrix transposition;

Step 6-2, the UCA that is r for radius, its steering vector is

Wherein, ζ=k _wrsin θ,

Beneficial effect: the present invention provides a kind of modeling method for improving the indoor distributed mimo system channel capacity of 3D, the 3D communication environment of Microcell in simulating chamber well, its channel parameter estimation result meets theory and experience, has expanded research and the application of 3d space statistical channel model.The present invention can be applied in distributed MIMO system, when base station BS is designed with oriented antenna, the ripple of the MS of take end reaches signal research as basis, channel parameter while deriving important sky, as Doppler's power spectrum (DS), space-time correlation and channel capacity etc., while utilizing channel empty, characterisitic parameter has been resolved spatial model parameter b/a, mechanism relation between D/a and oriented antenna main lobe angle 2 α, and studied many antennas of MIMO linear array (ULA:uniform linear array) and circular array (UCA:uniform circular array) the system transmission performance under 3d space statistical channel model simultaneously, thereby during for assessment multiple-antenna MIMO system sky, Processing Algorithm and emulation wireless communication system provide strong instrument.

Accompanying drawing explanation

Fig. 1 is for realizing the 3d space statistical channel model that the invention provides modeling method;

Fig. 2 is the scattering region distribution map of 3D model;

Fig. 3 is the profile of 3D statistical channel model on perpendicular;

Fig. 4 is the vertical view of 3D statistical channel model;

Fig. 5 is the Doppler frequency shift schematic diagram that the mobility of mobile station MS produces;

Fig. 6 (a) is four unit MIMO ULA Array Model (b) four unit UCA Array Model.

Fig. 7 is the AOA probability density distribution schematic diagram (a=100m, D=50m) of mobile station MS on horizontal plane;

Fig. 8 is the AOA probability density distribution schematic diagram (a=100m, D=50m) of MS on perpendicular;

Fig. 9 is the TOA/AOA joint probability density distribution schematic diagram (a=100m, b=50m, α=60 °) of mobile station MS on horizontal plane;

Figure 10 is the TOA/AOA joint probability density distribution schematic diagram (a=100m, b=50m, α=60 °) of MS on perpendicular;

To be perpendicular angle β M distribute on Doppler's power of mobile station MS Figure 11 affects schematic diagram (a=100m, b=50m, D=50m, α=60 °, φ _υ=90 °);

Figure 12 be main lobe angle α and spatial parameter b/a Doppler's power of mobile station MS is distributed affect schematic diagram (a=100m, D=50m, φ _υ=90 °);

Figure 13 is the direction of motion φ of mobile station MS _υand spatial parameter D/a Doppler's power is distributed affect schematic diagram (a=100m, b=50m, α=60 °);

Figure 14 is perpendicular angle β _mwith main lobe angle α on MIMO ULA array manifold correlation, affect schematic diagram (a=100m, b=50m, D=50m);

Figure 15 is perpendicular angle β _mwith main lobe angle α on UCA array manifold correlation, affect schematic diagram (a=100m, b=50m, D=50m);

Figure 16 is perpendicular angle β _mwith spatial parameter D/a on MIMO ULA array manifold correlation affect schematic diagram (a=100m, b=50m, D=50m);

Figure 17 is perpendicular angle β _mwith spatial parameter D/a on UCA array manifold correlation, affect schematic diagram (a=100m, b=50m, D=50m);

Figure 18 is perpendicular angle β _mwith the affect schematic diagram of main lobe angle α on MIMO ULA array channel capacity;

Figure 19 is perpendicular angle β _mwith main lobe angle α on UCA array channel capacity, affect schematic diagram (a=100m, b=50m, D=50m, N _r=4, SNR=20dB);

Figure 20 is perpendicular angle β _mwith spatial parameter D/a on MIMO ULA array channel capacity, affect schematic diagram (a=100m, b=50m, D=50m, N _r=4, SNR=20dB);

Figure 21 is perpendicular angle β _mwith spatial parameter D/a on UCA array channel capacity, affect schematic diagram (a=100m, b=50m, D=50m, N _r=4, SNR=20dB).

Embodiment

Below with reference to specific embodiment, technical scheme provided by the invention is elaborated, should understands following embodiment and only for the present invention is described, is not used in and limits the scope of the invention.

Modeling method about 3d space territory many antennas MIMO statistical channel provided by the invention, as shown in Figure 1 and Figure 2, comprising base station BS, mobile station MS and scattering object point s, specifically, comprise the steps: its 3d space statistical channel model

Step 1, we take BS and set up three-dimensional system of coordinate as initial point, and the main lobe angle of the directional antenna of BS end is 2 α, and the angle of BS on horizontal plane is φ _b, the angle on perpendicular is β _b; In like manner, the angle of MS on horizontal plane and perpendicular is respectively φ _mand β _m.Scattering object in scattering region is respectively r to the distance of BS and MS _band r _m.Suppose that the distance between BS and MS is D, the major axis of 3D channel model and minor axis length are respectively a and b, suppose D<a and b≤a.Under this model, continue to use uniform diffuser distributional assumption, all scattering objects are evenly distributed on scattering region I _regioninterior (grey vertical plane is the boundary of scattering region and non-scattering region), by numerical computations, the volume of trying to achieve scattering region is V=2a ²b α/3.

Step 2, has directional antenna based on base station BS end, makes scattering region I _regionpresentation space territory shape, for the needs of studying, by scattering region I _regionbe divided into P ₁and P ₂two parts.Wherein, P ₂refer to that the scattering object of electromagnetic signal on the vertical scattering boundary of grey reflexes to the propagation path region on MS, and P ₁scattering region I _regionremaining part.

As shown in Figure 3, during section PMQ in definition on vertical direction, from level angle and vertical angle, on level angle, start angle and termination point are respectively φ _t1and φ _t2, can be solved to

φ _t1=0，0≤β _M≤π／2 (1)

${\mathrm{\&phi;}}_{\mathrm{<figure-callout\; id="2"\; label="t"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">t</figure-callout>}2}=\left\{\begin{array}{cc}\arccos \left\{\frac{{\mathrm{<figure-callout\; id="2"\; label="PM"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">PM</figure-callout>}}^2+{\mathrm{QM}}^2-{\mathrm{<figure-callout\; id="2"\; label="PQ"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">PQ</figure-callout>}}^2}{2\mathrm{PM}\&times;\mathrm{QM}}\right\},& 0\&le;{\mathrm{\&beta;}}_M\&le;\arctan \left(\frac{b}{D\sin \mathrm{\&alpha;}}\right)\\ 0,& \arctan \left(\frac{b}{D\sin \mathrm{\&alpha;}}\right)\&le;{\mathrm{\&beta;}}_M\&le;\mathrm{\&pi;}/2\end{array}\right.---\left(2\right)$

Wherein, P and Q are section on perpendicular and the intersection point of scattering boundary, and the dihedral angle of PMQ plane and xoy plane is β _m, and

$\begin{array}{c}\mathrm{PQ}=\sqrt{a^2-\frac{a^2}{b^2}{\mathrm{<figure-callout\; id="2"\; label="D"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">D</figure-callout>}}^2{\mathrm{<figure-callout\; id="2"\; label="sin"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">sin</figure-callout>}}^2{\mathrm{\&alpha;}\tan }^2{\mathrm{\&beta;}}_M}\\ \mathrm{QM}=\sqrt{{\mathrm{<figure-callout\; id="2"\; label="d"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">d</figure-callout>}}^2+{\mathrm{<figure-callout\; id="2"\; label="D"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">D</figure-callout>}}^2{\mathrm{<figure-callout\; id="2"\; label="sin"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">sin</figure-callout>}}^2\mathrm{\&alpha;}{\tan }^2{\mathrm{\&beta;}}_M}\\ \mathrm{PM}=\sqrt{{\mathrm{<figure-callout\; id="2"\; label="D"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">D</figure-callout>}}^2+{\mathrm{<figure-callout\; id="2"\; label="D"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">D</figure-callout>}}^2{\mathrm{<figure-callout\; id="2"\; label="sin"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">sin</figure-callout>}}^2\mathrm{\&alpha;}{\tan }^2{\mathrm{\&beta;}}_M+{\mathrm{PQ}}^2-2D\&times;\mathrm{PQ}\cos \mathrm{\&alpha;}}\end{array}---\left(3\right)$

And the vertical angle beta t of plane P MQ is about parameter phi _mfunctional relation by following formula, calculate:

${\mathrm{\&beta;}}_t=\left\{\begin{array}{cc}{\cot }^{-1}\left(\frac{\mathrm{\&alpha;}D\csc (\mathrm{\&alpha;}+{\mathrm{\&phi;}}_m)\sin \mathrm{\&alpha;}}{b\sqrt{a^2-{\mathrm{<figure-callout\; id="2"\; label="D"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">D</figure-callout>}}^2{\csc }^2(\mathrm{\&alpha;}+{\mathrm{\&phi;}}_m){\mathrm{<figure-callout\; id="2"\; label="sin"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">sin</figure-callout>}}^2{\mathrm{\&phi;}}_m}}\right),& {\mathrm{\&phi;}}_1\&le;\left|{\mathrm{\&phi;}}_m\right|\&le;{\mathrm{\&phi;}}_2\\ 0,& \mathrm{otherwise}\end{array}\right.---\left(4\right)$

Wherein, angle φ ₁and φ ₂for β _m=0 o'clock, when 3D statistical channel model is converted to 2D model, the angle of MS end on xoy horizontal plane,

therefore, scattering region P ₁and P ₂can be expressed as

${\mathrm{<figure-callout\; id="1"\; label="P"\; filenames="140120170431.png,140120170435.png"\; state="\{\{state\}\}">P</figure-callout>}}_1\&RightArrow;\left\{\begin{array}{c}0\&le;{\mathrm{\&beta;}}_M\&le;{\mathrm{\&beta;}}_t\\ \mathrm{or}\\ {\mathrm{\&phi;}}_{\mathrm{<figure-callout\; id="1"\; label="t"\; filenames="140120170431.png,140120170435.png"\; state="\{\{state\}\}">t</figure-callout>}1}\&le;\left|{\mathrm{\&phi;}}_m\right|\&le;{\mathrm{\&phi;}}_{t2}\end{array}\right.{\mathrm{<figure-callout\; id="2"\; label="P"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">P</figure-callout>}}_2\&RightArrow;\left\{\begin{array}{c}{\mathrm{\&beta;}}_t\&le;{\mathrm{\&beta;}}_M\&le;\mathrm{\&pi;}/2\\ \mathrm{or}\\ {\mathrm{\&phi;}}_{\mathrm{<figure-callout\; id="2"\; label="t"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">t</figure-callout>}2}\&le;\left|{\mathrm{\&phi;}}_m\right|\&le;2\mathrm{\&pi;}-{\mathrm{\&phi;}}_{t2}\end{array}\right.---\left(5\right)$

Step 3: channel parameter in order to describe accurately this model important empty time, we start with from MS end, suppose that scattering object is evenly distributed in scattering region, and its distribution function is f (x _m, y _m, z _m)=1/V, therefore can obtain by Jacobi's change type

$p(r_m,{\mathrm{\&phi;}}_m,{\mathrm{\&beta;}}_m)=\frac{f(x_m,y_m,z_m)}{\left|J(x_m,y_m,z_m)\right|}{|}_{\begin{array}{c}{x}_m=r_m\cos {\mathrm{\&beta;}}_m\cos {\mathrm{\&phi;}}_m\\ y_m=r_m\cos {\mathrm{\&beta;}}_m\sin {\mathrm{\&phi;}}_m\\ z_m=r_m\sin {\mathrm{\&beta;}}_m\end{array}}=\frac{{\mathrm{<figure-callout\; id="2"\; label="r\; m"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">r</figure-callout>}}_m^{}\cos {\mathrm{\&beta;}}_m}{V}---\left(6\right)$

Wherein, MS is to the distance r of scattering boundary _mby following formula, solve

$r_m({\mathrm{\&phi;}}_m,{\mathrm{\&beta;}}_m)=\left\{\begin{array}{cc}D\sin \mathrm{\&alpha;}\csc (\mathrm{\&alpha;}+{\mathrm{\&phi;}}_m)\sec {\mathrm{\&beta;}}_m,& {\mathrm{<figure-callout\; id="1"\; label="P"\; filenames="140120170431.png,140120170435.png"\; state="\{\{state\}\}">P</figure-callout>}}_1\\ \frac{1}{{\mathrm{<figure-callout\; id="2"\; label="b"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">b</figure-callout>}}^2{\mathrm{<figure-callout\; id="2"\; label="cos"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">cos</figure-callout>}}^2{\mathrm{\&beta;}}_m+a^2{\mathrm{<figure-callout\; id="2"\; label="sin"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">sin</figure-callout>}}^2{\mathrm{\&beta;}}_m}& \\ \&times;\{{\mathrm{<figure-callout\; id="2"\; label="Db"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">Db</figure-callout>}}^2\cos {\mathrm{\&beta;}}_m\cos {\mathrm{\&phi;}}_m+\sqrt{{\left({\mathrm{<figure-callout\; id="2"\; label="Db"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">Db</figure-callout>}}^2\cos {\mathrm{\&beta;}}_m\cos {\mathrm{\&phi;}}_m\right)}^2}& \\ \stackrel{\&OverBar;}{-({\mathrm{<figure-callout\; id="2"\; label="b"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">b</figure-callout>}}^2{\mathrm{<figure-callout\; id="2"\; label="cos"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">cos</figure-callout>}}^2{\mathrm{\&beta;}}_m+a^2{\mathrm{<figure-callout\; id="2"\; label="sin"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">sin</figure-callout>}}^2{\mathrm{\&beta;}}_m)(b^2{D}^2-a^2{b}^2)}\},& P_2\end{array}\right.---\left(7\right)$

Formula (7) is updated in formula (6), and the AOA joint probability density function that can obtain MS is

$p({\mathrm{\&phi;}}_m,{\mathrm{\&beta;}}_m)=\left\{\begin{array}{cc}\frac{\cos {\mathrm{\&beta;}}_m}{3V}{\left\{D\sin \mathrm{\&alpha;}\csc (\mathrm{\&alpha;}+{\mathrm{\&phi;}}_m)\sec {\mathrm{\&beta;}}_m\right\}}^3,& {\mathrm{<figure-callout\; id="1"\; label="P"\; filenames="140120170431.png,140120170435.png"\; state="\{\{state\}\}">P</figure-callout>}}_1\\ \frac{\cos {\mathrm{\&beta;}}_m}{3V}\{\frac{1}{{\mathrm{<figure-callout\; id="2"\; label="b"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">b</figure-callout>}}^2{\mathrm{<figure-callout\; id="2"\; label="cos"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">cos</figure-callout>}}^2{\mathrm{\&beta;}}_m+a^2{\mathrm{<figure-callout\; id="2"\; label="sin"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">sin</figure-callout>}}^2{\mathrm{\&beta;}}_m}& \\ \&times;({\mathrm{<figure-callout\; id="2"\; label="Db"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">Db</figure-callout>}}^2\cos {\mathrm{\&beta;}}_m\cos {\mathrm{\&phi;}}_m+\sqrt{{\left({\mathrm{<figure-callout\; id="2"\; label="Db"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">Db</figure-callout>}}^2\cos {\mathrm{\&beta;}}_m\cos {\mathrm{\&phi;}}_m\right)}^2}& \\ \stackrel{\&OverBar;}{-({\mathrm{<figure-callout\; id="2"\; label="b"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">b</figure-callout>}}^2{\mathrm{<figure-callout\; id="2"\; label="cos"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">cos</figure-callout>}}^2{\mathrm{\&beta;}}_m+a^2{\mathrm{<figure-callout\; id="2"\; label="sin"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">sin</figure-callout>}}^2{\mathrm{\&beta;}}_m(b^2{D}^2-a^2{b}^2)}){\}}^3,& P_2\end{array}\right.---\left(8\right)$

Step 4: the electromagnetic signal to the reflection of any scattering object, from mobile station MS to base station BS between as there is wave path, arrive time delay τ=(r _m+ r _b)/c, wherein c is the light velocity, if there is direct path (LoS:Line of Sight), the minimum value that arrives time delay is τ _min=D/c, maximum delay is

${\mathrm{\&tau;}}_{\max }=\frac{1}{c}(a+\sqrt{a^2-2\mathrm{aD}\cos \mathrm{\&alpha;}+D^2})---\left(9\right)$

For the research of the delay character TOA of MS, similar with the research method of above-mentioned arrival angle AOA, formula (6) can be obtained by Jacobi's change type

$\begin{array}{c}p(\mathrm{\&tau;},{\mathrm{\&phi;}}_m,{\mathrm{\&beta;}}_m)=\frac{p(r_m,{\mathrm{\&phi;}}_m,{\mathrm{\&beta;}}_m)}{\left|J(r_m,{\mathrm{\&phi;}}_m,{\mathrm{\&beta;}}_m)\right|}{|}_{{\mathrm{<figure-callout\; id="1"\; label="P"\; filenames="140120170431.png,140120170435.png"\; state="\{\{state\}\}">P</figure-callout>}}_1,{\mathrm{<figure-callout\; id="2"\; label="P"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">P</figure-callout>}}_2}\\ =\frac{c{({\mathrm{<figure-callout\; id="2"\; label="c"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">c</figure-callout>}}^2{\mathrm{\&tau;}}^2-D^2)}^2({\mathrm{<figure-callout\; id="2"\; label="D"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">D</figure-callout>}}^2+{\mathrm{<figure-callout\; id="2"\; label="c"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">c</figure-callout>}}^2{\mathrm{\&tau;}}^2-2\mathrm{c\&tau;}D\cos {\mathrm{\&beta;}}_m\cos {\mathrm{\&phi;}}_m)\cos {\mathrm{\&beta;}}_m}{8V{(D\cos {\mathrm{\&beta;}}_m\cos {\mathrm{\&phi;}}_m-\mathrm{c\&tau;})}^4}\end{array}---\left(10\right)$

By formula (10) to β _mcarry out integration, the delay character TOA of mobile station MS and the joint probability density function of horizontal plane angle are

$\begin{array}{c}p(\mathrm{\&tau;},{\mathrm{\&phi;}}_m)={\&Integral;}_0^{\mathrm{\&pi;}/2}p(\mathrm{\&tau;},{\mathrm{\&phi;}}_m,{\mathrm{\&beta;}}_m)d{\mathrm{\&beta;}}_m\\ =\frac{{\mathrm{<figure-callout\; id="1"\; label="k"\; filenames="140120170431.png,140120170435.png"\; state="\{\{state\}\}">k</figure-callout>}}_1}{48\mathrm{<figure-callout\; id="5"\; label="V\; k"\; filenames="140120170435.png"\; state="\{\{state\}\}">V</figure-callout>}{k}_{}^{}{}_5^2{(k_4^2-{\mathrm{<figure-callout\; id="5"\; label="k"\; filenames="140120170435.png"\; state="\{\{state\}\}">k</figure-callout>}}_5^2)}^{7/2}}\&times;\{\sqrt{k_4^2-{\mathrm{<figure-callout\; id="5"\; label="k"\; filenames="140120170435.png"\; state="\{\{state\}\}">k</figure-callout>}}_5^2}(k_3{k}_5(2k_4^3+13k_4{\mathrm{<figure-callout\; id="5"\; label="k"\; filenames="140120170435.png"\; state="\{\{state\}\}">k</figure-callout>}}_5^2)+k_2(k_4^4-10k_4^2{\mathrm{<figure-callout\; id="5"\; label="k"\; filenames="140120170435.png"\; state="\{\{state\}\}">k</figure-callout>}}_5^2-6{\mathrm{<figure-callout\; id="5"\; label="k"\; filenames="140120170435.png"\; state="\{\{state\}\}">k</figure-callout>}}_5^4))\\ -6{\mathrm{<figure-callout\; id="5"\; label="k"\; filenames="140120170435.png"\; state="\{\{state\}\}">k</figure-callout>}}_5^2(k_3{k}_5(4k_4^2+{\mathrm{<figure-callout\; id="5"\; label="k"\; filenames="140120170435.png"\; state="\{\{state\}\}">k</figure-callout>}}_5^2)-k_2(k_4^3+4k_4{\mathrm{<figure-callout\; id="5"\; label="k"\; filenames="140120170435.png"\; state="\{\{state\}\}">k</figure-callout>}}_5^2))a\tanh \left(\frac{k_4+k_5}{\sqrt{k_4^2-{\mathrm{<figure-callout\; id="5"\; label="k"\; filenames="140120170435.png"\; state="\{\{state\}\}">k</figure-callout>}}_5^2}}\right)\}\end{array}---\left(11\right)$

By formula (10) to φ _mcarry out integration, the delay character TOA of mobile station MS and the joint probability density function of perpendicular angle are

$\begin{array}{c}p(\mathrm{\&tau;},{\mathrm{\&beta;}}_m)={\&Integral;}_0^{2\mathrm{\&pi;}}p(\mathrm{\&tau;},{\mathrm{\&phi;}}_b,{\mathrm{\&beta;}}_b)d{\mathrm{\&phi;}}_m\\ =\frac{k_6(k_7{\mathrm{<figure-callout\; id="8"\; label="k"\; filenames="140120170423.png,140120170427.png"\; state="\{\{state\}\}">k</figure-callout>}}_8^3-3k_2{\mathrm{<figure-callout\; id="8"\; label="k"\; filenames="140120170423.png,140120170427.png"\; state="\{\{state\}\}">k</figure-callout>}}_8^2{\mathrm{<figure-callout\; id="5"\; label="k"\; filenames="140120170435.png"\; state="\{\{state\}\}">k</figure-callout>}}_5+4k_7{k}_8{\mathrm{<figure-callout\; id="5"\; label="k"\; filenames="140120170435.png"\; state="\{\{state\}\}">k</figure-callout>}}_5^2-2k_2{\mathrm{<figure-callout\; id="5"\; label="k"\; filenames="140120170435.png"\; state="\{\{state\}\}">k</figure-callout>}}_5^3)}{8V{(k_8-k_5)}^3{({\mathrm{<figure-callout\; id="8"\; label="k"\; filenames="140120170423.png,140120170427.png"\; state="\{\{state\}\}">k</figure-callout>}}_8+k_5)}^4}\&times;\sqrt{1+\frac{2k_8}{-{\mathrm{<figure-callout\; id="8"\; label="k"\; filenames="140120170423.png,140120170427.png"\; state="\{\{state\}\}">k</figure-callout>}}_8+{\mathrm{<figure-callout\; id="5"\; label="k"\; filenames="140120170435.png"\; state="\{\{state\}\}">k</figure-callout>}}_5}}\mathrm{\&pi;}\end{array}---\left(12\right)$

Wherein,

$\left\{\begin{array}{c}{\mathrm{<figure-callout\; id="1"\; label="k"\; filenames="140120170431.png,140120170435.png"\; state="\{\{state\}\}">k</figure-callout>}}_1=c{({\mathrm{<figure-callout\; id="2"\; label="c"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">c</figure-callout>}}^2{\mathrm{\&tau;}}^2-D^2)}^2\\ {\mathrm{<figure-callout\; id="2"\; label="k"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">k</figure-callout>}}_2={\mathrm{<figure-callout\; id="2"\; label="c"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">c</figure-callout>}}^2{\mathrm{\&tau;}}^2+D^2\\ {\mathrm{<figure-callout\; id="3"\; label="k"\; filenames="140120170442.png,140120170446.png"\; state="\{\{state\}\}">k</figure-callout>}}_3=2\mathrm{c\&tau;}D\cos {\mathrm{\&phi;}}_m\\ k_4=D\cos {\mathrm{\&phi;}}_m\\ {\mathrm{<figure-callout\; id="5"\; label="k"\; filenames="140120170435.png"\; state="\{\{state\}\}">k</figure-callout>}}_5=\mathrm{c\&tau;}\end{array}\right.\mathrm{and}\left\{\begin{array}{c}{\mathrm{<figure-callout\; id="6"\; label="k"\; filenames="140120170458.png,140120170502.png"\; state="\{\{state\}\}">k</figure-callout>}}_6=c{({\mathrm{<figure-callout\; id="2"\; label="c"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">c</figure-callout>}}^2{\mathrm{\&tau;}}^2-D^2)}^2\cos {\mathrm{\&beta;}}_m\\ k_7={\mathrm{<figure-callout\; id="2"\; label="c"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">c</figure-callout>}}^2{\mathrm{\&tau;}}^2+D^2\\ {\mathrm{<figure-callout\; id="8"\; label="k"\; filenames="140120170423.png,140120170427.png"\; state="\{\{state\}\}">k</figure-callout>}}_8=2\mathrm{c\&tau;}D\cos {\mathrm{\&beta;}}_m\\ {\mathrm{<figure-callout\; id="9"\; label="k"\; filenames="140120170427.png"\; state="\{\{state\}\}">k</figure-callout>}}_9=D\cos {\mathrm{\&beta;}}_m\\ {\mathrm{<figure-callout\; id="10"\; label="k"\; filenames="140120170423.png,140120170427.png"\; state="\{\{state\}\}">k</figure-callout>}}_{10}=\mathrm{c\&tau;}\end{array}\right.;$

Step 5: in 3d space channel model, the mobility of mobile station MS can make receiving end signal produce Doppler frequency shift (as shown in Figure 5), wherein, the translational speed υ of MS, the MS included angle on horizontal plane _m, the angle β of plane P MQ and horizontal plane xoy _mthree has determined the motion direction of MS.By classical Clarke model, the probability density function of Doppler frequency shift can be derived as

$f_{\mathrm{\&upsi;}}\left(f\right)=\frac{p({\mathrm{\&phi;}}_{\mathrm{\&upsi;}}+|{\cos }^{-1}(f/f_m\cos \mathrm{\&beta;})|,{\mathrm{\&beta;}}_m)}{f_m\sqrt{{1-(f/f_m)}^2}}+\frac{p({\mathrm{\&phi;}}_{\mathrm{\&upsi;}}-|{\cos }^{-1}(f/f_m\cos \mathrm{\&beta;})|,{\mathrm{\&beta;}}_m)}{f_m\sqrt{{1-(f/f_m)}^2}},\left|f\right|\&le;f_m,---\left(13\right)$

Wherein, f _mfor the maximum doppler frequency of mobile station MS, φ _υmoving direction for MS.

Step 6: the geometry of MIMO multi-antenna array can be arbitrarily, the difference according to antenna element array element in spatial distribution form, can be divided into linear ULA and circle UCA etc.Linear array has simple in structure, and circular array has the feature of isotropic directivity, and that at moving communicating field, applies is more.At the incoming signal steric direction vector that receives MIMO ULA, be

Wherein,

the angle of incoming signal on horizontal plane; θ is the angle of incoming signal on perpendicular.And in formula (14), k _w=2 π/λ, d is antenna spacing, and λ is incoming signal wavelength, and L is receiving terminal number of antennas, [] ^trepresenting matrix transposition.The UCA that is r for radius, its steering vector is

In formula, ζ=k _wr sin θ,

Therefore,, for MIMO ULA and the UCA array antenna of space structure shown in Fig. 6, the Spatial fading correlation function between array element m and n can be expressed as

Step 7: in performance of mobile communication system is analyzed, channel capacity, fundamentally having determined the performance of wireless system, has great importance for system.Suppose that transmitting terminal does not have any channel information, and receiving terminal is known state, channel capacity can be expressed as

$C={\mathrm{<figure-callout\; id="2"\; label="log"\; filenames="140120170338.png,140120170409.png"\; state="\{\{state\}\}">log</figure-callout>}}_2\det (I_{\mathrm{Nr}}+\frac{P}{N_t{\mathrm{\&sigma;}}^2}{\mathrm{HH}}^H)---\left(17\right)$

Shi Zhong C unit is (bits/s/Hz), I _nrfor N _rdimension unit matrix, P/ σ ²for channel signal to noise ratio snr, N _tfor transmitting terminal antenna amount and N _rfor receiving terminal antenna amount.For the mimo channel of space correlation, channel matrix H can be expressed as

$H={\mathrm{<figure-callout\; id="1"\; label="R\; r"\; filenames="140120170431.png,140120170435.png"\; state="\{\{state\}\}">R</figure-callout>}}_r^{}{H}_w{\left({\mathrm{<figure-callout\; id="1"\; label="R\; t"\; filenames="140120170431.png,140120170435.png"\; state="\{\{state\}\}">R</figure-callout>}}_t^{}\right)}^T---\left(18\right)$

Wherein, R _rfor correlation matrix between the array element of receiving terminal, R _tfor correlation matrix between transmitting terminal array element, H _wfor the multiple gaussian random matrix with distributing.Subscript ^tand ^htransposition and the conjugate transpose of difference representing matrix.

Based on above-mentioned steps, we define several groups of different model parameters, will in each important parameter expression formula substitution of the channel obtaining by above-mentioned modeling process matlab, carry out numerical simulation and calculating, obtain result as shown in Fig. 7 to 21:

Fig. 7 is the probability density distribution figure schematic diagram of MS on horizontal plane (a=100m wherein, D=50m), based on BS end, have directional antenna, and due to the symmetry (seeing Fig. 4) of space structure, make AOA probability density distribution figure be cut out angle non-existent two parts region, and left and right is in symmetry status.From figure, can also find, ripple reaches the φ of signal AOA distribution function in spatial model _m=φ ₂there is noncontinuity feature in point, and at (0, φ ₂) scope Nei Boda signal AOA probability density presents the trend of first increases and then decreases, and at (φ ₂, π) in, probability density presents the trend of monotone decreasing.The probability density distribution schematic diagram of MS on perpendicular be (wherein a=100m, D=50m) as shown in Figure 8, from this figure, can find, works as β _m=0 o'clock, under omnidirectional antenna away from the scattering object of MS many than in directional antenna situation, make the reflection of electromagnetic signal and refraction probability larger, it is relatively little that ripple reaches signal AOA probability density.And be accompanied by β _mbe increased to tan ^-1(b/D), time, the probability density under omni-directional antenna case and the probability density under directional antenna are tending towards overlapping.From figure, can also find, ripple reaches the β of signal AOA distribution function in spatial model _m=tan ^-1(b/D) there is noncontinuity feature, and be accompanied by β _mincrease, in the scope Nei Boda signal AOA of (0, pi/2) probability density monotone decreasing to zero.

Fig. 9 is the TOA/AOA probability density distribution schematic diagram (a=100m of MS on horizontal plane, b=50m, α=60 °), because scattering region is about xoz plane symmetry, make joint probability distribution figure be dug up non-existent two panel region of delay character τ, and left and right is in symmetry status.From figure, can also find, it is at low-angle or minimal time delay τ substantially that mobile station MS receives signal _minplace, and at wide-angle (φ _m=180 °) and its probability density distribution of maximum delay place less, and its probability density function values exists

point reaches peak value, wherein, and τ _p=D/c and the TOA/AOA probability density distribution schematic diagram of MS on perpendicular be (a=100m, b=50m, α=60 °) as shown in figure 10, when propagation delay time τ surpasses maximum delay τ _max, while having surpassed scattering region, its probability density is zero.From figure, can also find, as perpendicular angle β _mor propagation delay time τ is when increase, away from the scattering object of MS end, increase gradually, the reflection of electromagnetic signal and refraction probability are larger, make TOA probability density relatively low; In this figure, also shown that in addition joint probability function value exists

reach peak value, wherein, τ _p=D/c and

Perpendicular angle β _mwhat Doppler's power of mobile station MS was distributed affects schematic diagram (a=100m, b=50m, D=50m, α=60 °, φ as shown in figure 11 _υ=90 °).From figure, can find, work as φ _υ=90 °, when the moving direction of MS is perpendicular to direct path LoS, scattering object can produce forward and reverse frequency displacement, and at β _m≤ β _ttime, perpendicular angle β no matter _mvalue how, its energy all can concentrate on Doppler frequency near zero point.And at β _min the time of=30 °, be accompanied by the increase (being the distance increase between BS and MS) of parameter D/a, Doppler's performance number increases gradually.From this figure, can also find, be accompanied by β _mincrease, the energy that Doppler concentrates on zero point reduces gradually, this is because the section that ripple reaches signal in the vertical direction reduces gradually, and electromagnetic signal reflection probability is reduced, therefore cause Doppler's power spectrum energy to reduce, and works as β _m>=β _ttime, Doppler's power distributes and is tending towards gradually classical U Clarke model.Above-mentioned simulation result meets the research of (Nawaz) in the past, and this shows that with regard to surface modeling method provided by the invention is applicable to describe 3d space statistical channel model parameter.

Main lobe angle α and spatial parameter b/a distribute on Doppler's power of mobile station MS affects schematic diagram (a=100m, D=50m, φ as shown in figure 12 _υ=90 °), when mobile station MS moving direction is vertical with LoS path, because scattering region is about xoz plane symmetry, so Doppler's power spectrum is symmetrical about frequency zero.While increasing to α=180 ° along with main lobe angle, Doppler's power spectrum is tending towards Clarke U-shaped classical model gradually.And when main lobe angle α=60 °≤180 °, when base station BS is designed with directional antenna, signal is relatively less at the multipath component in space, now the dimensional energy of Doppler's power spectrum mainly concentrates near frequency zero.In figure, also shown the impact that spatial parameter b/a distributes on Doppler's power, from figure, can also find, when b/a is larger, Doppler's power is relatively large, this is because be accompanied by the increase of parameter b/a, it is many that scattering object in scattering region becomes gradually, the reflection of electromagnetic signal and refraction probability are larger, thereby cause the energy of Doppler's power spectrum relatively large, this is consistent with result of study in the past, the channel parameter estimation result that shows this model meets theory and experience, has expanded research and the application of 3d space statistical channel model.

Figure 13 is the direction of motion φ of mobile station MS _υand spatial parameter D/a Doppler's power is distributed affect schematic diagram (a=100m, b=50m, α=60 °).Work as φ _υ=5 °, during the mobile arrival bearing of MS, Doppler's power spectrum that ripple reaches signal tilts and presents asymmetrical shape, and in the value of parameter D/a hour, Doppler's power spectrum is tilted to the right, and illustrate that Doppler frequency is born component proportion larger; And be accompanied by the increase of parameter D/a, Doppler frequency positive component proportion increases gradually, the power spectrum that makes Doppler gradually left direction tilt.Otherwise at φ _υ=175 °, when the moving direction of MS deviates from arrival bearing, with above-mentioned φ _υ=5 ° produce antipodal effect.From figure, can also find (φ when MS moves perpendicular to direct path _υ=90 °), regardless of the distance between MS and BS, the dimensional energy of Doppler's power spectrum all concentrates near frequency zero.

Figure 14, Figure 15, Figure 16, Figure 17 are respectively scattering object and are uniformly distributed down, the spatial coherence distribution schematic diagram (a=100m, b=50m, D=50m) between MIMO ULA and UCA array element (1,2).From figure, can find, spatial coherence reduces along with the increase of d/ λ or r/ λ and levels off to 0 value, and is accompanied by vertical angle beta _mincrease, the correlation between array dies down gradually.In figure, also shown main lobe angle α and the impact of spatial parameter D/a on correlation, from figure, can also find, be accompanied by reducing of main lobe angle α, the correlation between array weakens gradually; And being accompanied by the increase of distance between base station BS and MS (D/a), correlation is grow gradually.

Figure 18, Figure 19, Figure 20, Figure 21 are respectively distribution situation schematic diagram (a=100m, b=50m, D=50m, the N of MIMO ULA and UCA aerial array channel capacity _r=4, SNR=20dB).From figure, can find, at two array element distance d/ λ or r/ λ hour, the correlation between array element is larger, and channel capacity is less.Main lobe angle α and the impact of spatial parameter D/a on channel capacity in figure, have also been shown, from figure, can also find, be accompanied by reducing of main lobe angle α, the correlation between array weakens gradually, channel capacity is increased, thereby embody the superiority of directional antenna; In addition, be accompanied by the increase of distance between base station BS and MS, correlation is grow gradually, and channel capacity reduces gradually.

From illustrating above, emulation experiment numerical result and 2D and the contrast of 3D multidiameter fading channel show that the channel parameter estimation result of this model meets theory and experience, by the step provided by the invention 3D communication environment of Microcell in simulating chamber well, research and the application of 3d space statistical channel model have been expanded.The present invention can computing 3d space statistical channel model important empty time channel parameter, as ripple reaches the direction of arrival degree (AOA:Angle of Arrival) of signal on horizontal plane and perpendicular, the time of advent (TOA:Time Of Arrival), the space-time correlation of Doppler effect (DS:Doppler Shift) and arriving signal, and studied many antennas of MIMO linear array (ULA:uniform linear array) and circular array (UCA:uniform circular array) the system transmission performance under 3d space statistical channel model simultaneously.In addition, further every input parameter of research model and important relation between channel parameter when empty, from the analogous diagram of implementing to provide, can obviously find out, 3d space is provided with the capacity that directional antenna can improve space channel effectively, and, reduce the capacity that space D between BS and MS can further improve channel.The present invention provides a kind of channel modeling method comparatively exactly for 3D interior space statistical channel, thus for assessment multiple-antenna MIMO system when empty Processing Algorithm and emulation wireless communication system strong instrument is provided, there is far-reaching Research Significance.

The disclosed technological means of the present invention program is not limited only to the disclosed technological means of above-mentioned execution mode, also comprises the technical scheme being comprised of above technical characterictic combination in any.It should be pointed out that for those skilled in the art, under the premise without departing from the principles of the invention, can also make some improvements and modifications, these improvements and modifications are also considered as protection scope of the present invention.

Claims (3)

1. the modeling method based on three dimensions territory many antennas MIMO statistical channel, is characterized in that, comprises the steps:

Step 1, the base station BS of take is set up three-dimensional system of coordinate as initial point, and wherein the main lobe angle of the directional antenna based on BS end is 2 α, and the angle of BS on horizontal plane is φ _b, the angle on perpendicular is β _b; The angle of MS on horizontal plane and perpendicular is respectively φ _mand β _m, the scattering object in scattering region is respectively r to the distance of BS and MS _band r _m, the distance between BS and MS is D, the major axis of 3D channel model and minor axis length are respectively a and b, and D<a, b≤a, all scattering objects are evenly distributed in scattering region, and the volume of scattering region is V=2a ²b α/3;

Step 2, is divided into P by scattering region ₁and P ₂two parts, wherein, P ₂refer to that the scattering object of electromagnetic signal on the vertical scattering boundary of grey reflexes to the propagation path region on MS, and P ₁scattering region I _regionremaining part;

Scattering region P ₁and P ₂be expressed as

$P_1\&RightArrow;\left\{\begin{array}{c}0\&le;{\mathrm{\&beta;}}_M\&le;{\mathrm{\&beta;}}_t\\ \mathrm{or}\\ {\mathrm{\&phi;}}_{t1}\&le;\left|{\mathrm{\&phi;}}_m\right|\&le;{\mathrm{\&phi;}}_{t2}\end{array}\right.P_2\&RightArrow;\left\{\begin{array}{c}{\mathrm{\&beta;}}_t\&le;{\mathrm{\&beta;}}_M\&le;\mathrm{\&pi;}/2\\ \mathrm{or}\\ {\mathrm{\&phi;}}_{t2}\&le;\left|{\mathrm{\&phi;}}_m\right|\&le;2\mathrm{\&pi;}-{\mathrm{\&phi;}}_{t2}\end{array}\right.$

Wherein, β _tfor the vertical angle of plane P MQ,

${\mathrm{\&beta;}}_t=\left\{\begin{array}{cc}{\cot }^{-1}\left(\frac{\mathrm{\&alpha;}D\csc (\mathrm{\&alpha;}+{\mathrm{\&phi;}}_m)\sin \mathrm{\&alpha;}}{b\sqrt{a^2-D^2{\csc }^2(\mathrm{\&alpha;}+{\mathrm{\&phi;}}_m){\sin }^2{\mathrm{\&phi;}}_m}}\right),& {\mathrm{\&phi;}}_1\&le;\left|{\mathrm{\&phi;}}_m\right|\&le;{\mathrm{\&phi;}}_2\\ 0,& \mathrm{otherwise}\end{array}\right.,$

Wherein, angle φ ₁and φ ₂for β _m=0 o'clock, the angle of MS end on xoy horizontal plane,

φ _t1for the start angle of the section PMQ on vertical direction on level angle,

φ _t1=0，0≤β _M≤π／2，

φ _t2for the termination point of the section PMQ on vertical direction on level angle,

${\mathrm{\&phi;}}_{t2}=\left\{\begin{array}{cc}\arccos \left\{\frac{{\mathrm{PM}}^2+{\mathrm{QM}}^2-{\mathrm{PQ}}^2}{2\mathrm{PM}\&times;\mathrm{QM}}\right\},& 0\&le;{\mathrm{\&beta;}}_M\&le;\arctan \left(\frac{b}{D\sin \mathrm{\&alpha;}}\right)\\ 0,& \arctan \left(\frac{b}{D\sin \mathrm{\&alpha;}}\right)\&le;{\mathrm{\&beta;}}_M\&le;\mathrm{\&pi;}/2\end{array}\right.,$

Wherein, P and Q are section on perpendicular and the intersection point of scattering boundary, and the dihedral angle of PMQ plane and xoy plane is β _m, and

$\begin{array}{c}\mathrm{PQ}=\sqrt{a^2-\frac{a^2}{b^2}{D}^2{\sin }^2{\mathrm{\&alpha;}\tan }^2{\mathrm{\&beta;}}_M}\\ \mathrm{QM}=\sqrt{d^2+D^2{\sin }^2\mathrm{\&alpha;}{\tan }^2{\mathrm{\&beta;}}_M}\\ \mathrm{PM}=\sqrt{D^2+D^2{\sin }^2\mathrm{\&alpha;}{\tan }^2{\mathrm{\&beta;}}_M+{\mathrm{PQ}}^2-2D\&times;\mathrm{PQ}\cos \mathrm{\&alpha;}}\end{array},$

Step 3: by the distribution function f (x of scattering object _m, y _m, z _m)=1/V obtains by Jacobi's change type:

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

Wherein, r _mfor the distance of MS to scattering boundary, by following formula, calculate:

$r_m({\mathrm{\&phi;}}_m,{\mathrm{\&beta;}}_m)=\left\{\begin{array}{cc}D\sin \mathrm{\&alpha;}\csc (\mathrm{\&alpha;}+{\mathrm{\&phi;}}_m)\sec {\mathrm{\&beta;}}_m,& P_1\\ \frac{1}{b^2{\cos }^2{\mathrm{\&beta;}}_m+a^2{\sin }^2{\mathrm{\&beta;}}_m}& \\ \&times;\{{\mathrm{Db}}^2\cos {\mathrm{\&beta;}}_m\cos {\mathrm{\&phi;}}_m+\sqrt{{\left({\mathrm{Db}}^2\cos {\mathrm{\&beta;}}_m\cos {\mathrm{\&phi;}}_m\right)}^2}& \\ \stackrel{\&OverBar;}{-(b^2{\cos }^2{\mathrm{\&beta;}}_m+a^2{\sin }^2{\mathrm{\&beta;}}_m)(b^2{D}^2-a^2{b}^2)}\},& P_2\end{array}\right.$

By r _mjacobi's change type in this step of substitution, obtains the AOA joint probability density function of MS:

$p({\mathrm{\&phi;}}_m,{\mathrm{\&beta;}}_m)=\left\{\begin{array}{cc}\frac{\cos {\mathrm{\&beta;}}_m}{3V}{\left\{D\sin \mathrm{\&alpha;}\csc (\mathrm{\&alpha;}+{\mathrm{\&phi;}}_m)\sec {\mathrm{\&beta;}}_m\right\}}^3,& P_1\\ \frac{\cos {\mathrm{\&beta;}}_m}{3V}\{\frac{1}{b^2{\cos }^2{\mathrm{\&beta;}}_m+a^2{\sin }^2{\mathrm{\&beta;}}_m}& \\ \&times;({\mathrm{Db}}^2\cos {\mathrm{\&beta;}}_m\cos {\mathrm{\&phi;}}_m+\sqrt{{\left({\mathrm{Db}}^2\cos {\mathrm{\&beta;}}_m\cos {\mathrm{\&phi;}}_m\right)}^2}& \\ \stackrel{\&OverBar;}{-(b^2{\cos }^2{\mathrm{\&beta;}}_m+a^2{\sin }^2{\mathrm{\&beta;}}_m(b^2{D}^2-a^2{b}^2)}){\}}^3,& P_2\end{array}\right.;$

Step 4: Jacobi's change type in step 3 is obtained by Jacobi's change type again:

$\begin{array}{c}p(\mathrm{\&tau;},{\mathrm{\&phi;}}_m,{\mathrm{\&beta;}}_m)=\frac{p(r_m,{\mathrm{\&phi;}}_m,{\mathrm{\&beta;}}_m)}{\left|J(r_m,{\mathrm{\&phi;}}_m,{\mathrm{\&beta;}}_m)\right|}{|}_{P_1,P_2}\\ =\frac{c{(c^2{\mathrm{\&tau;}}^2-D^2)}^2(D^2+c^2{\mathrm{\&tau;}}^2-2\mathrm{c\&tau;}D\cos {\mathrm{\&beta;}}_m\cos {\mathrm{\&phi;}}_m)\cos {\mathrm{\&beta;}}_m}{8V{(D\cos {\mathrm{\&beta;}}_m\cos {\mathrm{\&phi;}}_m-\mathrm{c\&tau;})}^4}\end{array}$

Wherein, τ is for arriving time delay, minimum value τ _min=D/c, maximum delay is

By above formula to β _mcarry out integration, obtain the delay character TOA of mobile station MS and the joint probability density function of horizontal plane angle is:

$\begin{array}{c}p(\mathrm{\&tau;},{\mathrm{\&phi;}}_m)={\&Integral;}_0^{\mathrm{\&pi;}/2}p(\mathrm{\&tau;},{\mathrm{\&phi;}}_m,{\mathrm{\&beta;}}_m)d{\mathrm{\&beta;}}_m\\ =\frac{k_1}{48V{k}_5^2{(k_4^2-k_5^2)}^{7/2}}\&times;\{\sqrt{k_4^2-k_5^2}(k_3{k}_5(2k_4^3+13k_4{k}_5^2)+k_2(k_4^4-10k_4^2{k}_5^2-6k_5^4))\\ -6k_5^2(k_3{k}_5(4k_4^2+k_5^2)-k_2(k_4^3+4k_4{k}_5^2))a\tanh \left(\frac{k_4+k_5}{\sqrt{k_4^2-k_5^2}}\right)\}\end{array}$

To φ _mcarry out integration, obtain the delay character TOA of MS and the joint probability density function of perpendicular angle is:

$\begin{array}{c}p(\mathrm{\&tau;},{\mathrm{\&beta;}}_m)={\&Integral;}_0^{2\mathrm{\&pi;}}p(\mathrm{\&tau;},{\mathrm{\&phi;}}_b,{\mathrm{\&beta;}}_b)d{\mathrm{\&phi;}}_m\\ =\frac{k_6(k_7{k}_8^3-3k_2{k}_8^2{k}_5+4k_7{k}_8{k}_5^2-2k_2{k}_5^3)}{8V{(k_8-k_5)}^3{(k_8+k_5)}^4}\&times;\sqrt{1+\frac{2k_8}{-k_8+k_5}}\mathrm{\&pi;}\end{array}$

Wherein,

$\left\{\begin{array}{c}{k}_1=c{(c^2{\mathrm{\&tau;}}^2-D^2)}^2\\ k_2=c^2{\mathrm{\&tau;}}^2+D^2\\ k_3=2\mathrm{c\&tau;}D\cos {\mathrm{\&phi;}}_m\\ k_4=D\cos {\mathrm{\&phi;}}_m\\ k_5=\mathrm{c\&tau;}\end{array}\right.\mathrm{and}\left\{\begin{array}{c}{k}_6=c{(c^2{\mathrm{\&tau;}}^2-D^2)}^2\cos {\mathrm{\&beta;}}_m\\ k_7=c^2{\mathrm{\&tau;}}^2+D^2\\ k_8=2\mathrm{c\&tau;}D\cos {\mathrm{\&beta;}}_m\\ k_9=D\cos {\mathrm{\&beta;}}_m\\ k_{10}=\mathrm{c\&tau;}\end{array}\right.;$

Step 5: by classical Clarke model, the probability density function of Doppler frequency shift is as follows:

$f_{\mathrm{\&upsi;}}\left(f\right)=\frac{p({\mathrm{\&phi;}}_{\mathrm{\&upsi;}}+|{\cos }^{-1}(f/f_m\cos \mathrm{\&beta;})|,{\mathrm{\&beta;}}_m)}{f_m\sqrt{{1-(f/f_m)}^2}}+\frac{p({\mathrm{\&phi;}}_{\mathrm{\&upsi;}}-|{\cos }^{-1}(f/f_m\cos \mathrm{\&beta;})|,{\mathrm{\&beta;}}_m)}{f_m\sqrt{{1-(f/f_m)}^2}},\left|f\right|\&le;f_m,$

Wherein, f _mthe maximum doppler frequency of mobile station MS, φ _υbe the moving direction of MS, υ is translational speed, β _mangle for plane P MQ and horizontal plane xoy;

Step 6: calculate the Spatial fading correlation function between array element m and n by following formula:

Step 7: calculate channel capacity by following formula:

$C={\log }_2\det (I_{\mathrm{Nr}}+\frac{P}{N_t{\mathrm{\&sigma;}}^2}{\mathrm{HH}}^H)$

Shi Zhong C unit is (bits/s/Hz), I _nrfor N _rdimension unit matrix, P/ σ ²for channel signal to noise ratio snr, N _tfor transmitting terminal antenna amount and N _rfor receiving terminal antenna amount.

2. the modeling method based on three dimensions territory many antennas MIMO statistical channel according to claim 1, is characterized in that: in described step 7, for the mimo channel of space correlation, channel matrix H is expressed as

$H=R_r^{1/2}{H}_w{\left(R_t^{1/2}\right)}^T$

Wherein, R _rfor correlation matrix between the array element of receiving terminal, R _tfor correlation matrix between transmitting terminal array element, H _wfor the multiple gaussian random matrix with distributing, subscript ^tand ^htransposition and the conjugate transpose of difference representing matrix.

3. the modeling method based on three dimensions territory many antennas MIMO statistical channel according to claim 1 and 2, is characterized in that, in described step 6

by following two kinds of modes, obtain:

Step 6-1, when MIMO array is linear ULA, the incoming signal steric direction vector that receives MIMO ULA is

Wherein,

the angle of incoming signal on horizontal plane; θ is the angle of incoming signal on perpendicular,

d is antenna spacing, and λ is incoming signal wavelength, and L is receiving terminal number of antennas, [] ^trepresenting matrix transposition;

Step 6-2, when MIMO array is radius while being the UCA of r, its incoming signal steric direction vector is:

Wherein, ζ=k _wrsin θ,

Images