CN114362852B - Doppler parameter estimation method based on improved SAGE - Google Patents
Doppler parameter estimation method based on improved SAGE Download PDFInfo
- Publication number
- CN114362852B CN114362852B CN202111474809.4A CN202111474809A CN114362852B CN 114362852 B CN114362852 B CN 114362852B CN 202111474809 A CN202111474809 A CN 202111474809A CN 114362852 B CN114362852 B CN 114362852B
- Authority
- CN
- China
- Prior art keywords
- doppler
- multipath
- doppler frequency
- estimation
- range
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Abstract
The invention relates to a Doppler parameter estimation method based on improved SAGE, which reasonably assumes a probability density function of a Doppler power spectrum by utilizing prior information of scatterer distribution in a channel, sets an initial value and carries out layered estimation on a multipath Doppler parameter in an SAGE algorithm based on prior probability distribution of Doppler frequency, and preferentially estimates multipath components in an interval with high prior probability of the Doppler frequency. Compared with the prior art, the method reduces the influence of the period continuation part in the likelihood spectrum on the estimation process by exchanging the estimation sequence of the multipath, and has the advantage of high reliability.
Description
Technical Field
The invention relates to the field of wireless channel parameter estimation, in particular to a Doppler parameter estimation method based on improved SAGE.
Background
With the intensive research of wireless channels, various channel parameter estimation algorithms are developed vigorously, including a space alternating generalized expectation-maximization algorithm, namely a SAGE algorithm. However, the estimation algorithm does not take into account the case where the doppler estimation range is limited. When the transceiving terminals move at a high speed relatively or objects in the environment move at a high speed, the doppler frequency of the multipath may exceed the maximum doppler estimation range which can be reached by snapshot measurement data (which is determined by the time interval between adjacent frames in a snapshot, and the time interval may be limited by a certain length no matter the measurement is active or passive), so that the estimation process is affected by a periodically extended likelihood peak in a likelihood spectrum, and the estimation result is wrong. Although the literature indicates that for MIMO channel measurement in TDM mode, the likelihood spectrum is measured by a measurement period T cy Switching period T of transmitting terminal t Receiving end switching period T r And the like. According to the document "Doppler frequency estimation for channel sounding using switched multi-element transmission and received antennas" (Yin X, flow B H, journal P, et al. Ieee Global Telecommunications conference. Ieee, 2004), when the parameters of the multipath components other than the l-th multipath component are perfectly estimated, the likelihood function of the multipath component l can be written as:
wherein, the first and the second end of the pipe are connected with each other,
is the l-th multipath component with respect to the Doppler frequency v l The likelihood function of (a) is, device for selecting or keeping>
Is the hidden data for the lth multipath component>
y (t) is the received signal>
For the multi-path signal component reconstructed on the basis of the parameter estimate for the currently obtained l-th multi-path component, ->
I is the number of measurement cycles, M 1 And M 2 Respectively representing the number of arrays at the transmitting end and the receiving end, P representing the power of the transmitted signal, T s For the transmission time, α, of a single transceiving antenna pair l Is the amplitude of the ith multipath component, c 1 (Ω 1,l ) And c 2 (Ω 2,l ) With respect to the departure angle omega of the transmitting and receiving arrays, respectively 1,l Angle of sum of arrival omega 2,l Guide vector of v l ' is the true Doppler shift, V (V), of the l-th multipath component l ) Is the noise term.
Wherein, it is worth noting
Thus, the value of the period extension of the likelihood spectrum is taken from T t And T r The determined influence of the envelope of the skin. Therefore, theoretically, the value of the likelihood spectrum of the periodic extension is lower than that of the main peak, and meanwhile, in the expectation step (E-step), the periodic extension of the likelihood spectrum is reconstructed simultaneously in the signal reconstruction process, and the periodic extension is eliminated simultaneously when hidden data is solved. Thus, the theoretical maximum Doppler estimation range is greater than T cy Is determined from the reciprocal ofIn (c) is used.
However, in practical cases, according to different T t And T r And (3) value taking, the outer envelope can be very flat, and meanwhile, the main peak and the period continuation peak in the likelihood spectrum fluctuate due to the existence of channel noise, and the estimation error of the SAGE algorithm is added. The comprehensive action of the factors can make the main peak and the periodic continuation peak in the likelihood function very close, which affects the sequence interference elimination process and makes the estimation wrong.
Disclosure of Invention
The invention aims to overcome the defect that the cycle extension in the likelihood spectrum in the prior art has adverse effect on the estimation process, and provides a Doppler parameter estimation method based on improved SAGE with high reliability.
The purpose of the invention can be realized by the following technical scheme:
the invention provides a Doppler parameter estimation method based on improved SAGE, which reasonably assumes a probability density function of a Doppler power spectrum by utilizing prior information of scatterer distribution in a channel, sets an initial value and carries out layered estimation on a multipath Doppler parameter in a space alternative generalized expectation maximization SAGE algorithm based on prior probability distribution of Doppler frequency, and preferentially estimates multipath components in an interval with high prior probability of the Doppler frequency.
Preferably, the method comprises the steps of:
s1, determining initial values and estimation ranges of relevant parameters except Doppler frequency, and giving prior distribution of multipath Doppler frequency; setting the number of layers, and determining the Doppler frequency estimation range and the initial value of each layer;
s2, utilizing an initial Doppler frequency estimation range, adopting expectation maximization on a current estimation path, and updating and acquiring parameter estimation of a current multipath component based on a coordinate state;
s3, judging whether the Doppler frequency estimation result of the current multipath is smaller than a set threshold value from the boundary of the Doppler frequency estimation range of the current layer: if yes, and the current Doppler frequency estimation range does not reach the maximum Doppler estimation range, expanding a layer of Doppler frequency estimation range, and repeating the step S2; otherwise, finishing the parameter estimation of the current path, and turning to the step S4;
s4, estimating parameters of the next multipath component, and repeating the steps S2 to S3 until all paths are estimated;
s5, comparing a difference value of likelihood function values between two iterations with a set threshold value, and judging whether SAGE iteration estimation is converged; if not, repeating the steps to perform next iterative update, otherwise, finishing the estimation process to judge whether SAGE iterative estimation is converged, if not, repeating S2-S4 to perform next iterative update, otherwise, finishing the estimation process.
Preferably, the initial value and the estimation range of the relevant parameter in step S1 include a multipath time delay, a complex amplitude value and a multipath number outside the multipath doppler frequency.
Preferably, said step S1 comprises the following sub-steps:
step S11, obtaining an initial Doppler frequency estimation range:
wherein c is the speed of light, f c In order to measure the center frequency of the system,
for the relative radial movement speed of the transmitting and receiving ends, T f For measuring adjacent frame time intervals, v, of signals max The maximum relative motion speed of the receiving and transmitting end in the current snapshot measurement process is obtained;
step S12, determining the range of the theoretical Doppler frequency:
wherein v is max For the maximum relative of the receiving and transmitting ends in the current snapshot measurement processMagnitude of movement velocity, c is speed of light, f c Is the center frequency of the measuring system;
step S13, determining the Doppler frequency estimation range of each layer:
the widened doppler estimation range in each layer estimation is as follows:
wherein the content of the first and second substances,
and &>
Respectively an upper bound and a lower bound of an initial Doppler frequency estimation range; inf (F) D,max ) And sup (F) D,max ) An upper bound and a lower bound of a maximum doppler estimation range, respectively; k is a set layering number;
the k-th layer doppler frequency estimation range is:
wherein v is max The maximum relative movement speed of the receiving and transmitting end in the current snapshot measurement process, c is the light speed, f c Is the center frequency of the measurement system;
and &>
Upper and lower bounds, f, respectively, of the initial Doppler frequency estimation range D,layer A widened Doppler estimation range is estimated in each layer;
step S14, setting the initial value grid point interval of the multipath on the Doppler domain as follows:
wherein the content of the first and second substances,
and &>
Upper and lower bounds, f, of the k-th layer Doppler frequency estimation range, respectively D,sidelobe Is the side lobe width of the doppler domain likelihood function.
Preferably, the doppler domain likelihood function in step S14 adopts Bartlett likelihood spectrum, and its side lobe width f D,sidelobe The expression is as follows:
wherein, T s Is the duration of a snapshot.
Preferably, the step S1 is:
s11, constructing a geodetic coordinate system, and determining the positions of a ground end and a receiving end in the coordinate system:
constructing a geodetic coordinate system by taking one point on a map as a coordinate O point, taking a ground plane as an XY plane and taking a ground perpendicular line as a Z axis, wherein
A set of orthonormal bases under the geodetic coordinate system;
determining the geometric position of the ground end at the current moment in the geodetic coordinate system
The geometric position of the receiving terminal at the current moment is->
S12, acquiring the maximum multipath time delay based on the channel impact response h (t, tau) and the power time delay spectrum P (t, tau) of the current timeExtension tau max So as to obtain the maximum multipath propagation distance c tau in the channel max Wherein c is the speed of light;
the expression of the channel impulse response h (t, τ) and the power delay spectrum P (t, τ) at the current moment is as follows:
h(t,τ)＝IFFT(Y(f)*conj(U(f))/P(U(f)))
P(t,τ)＝E[|h(t,τ)| 2 ]
wherein, IFFT (·) is inverse discrete fourier transform, Y (f) is frequency domain representation of received signal snapshot Y (t) at the current time, and U (f) is frequency domain representation of transmitted signal snapshot U (t) at the current time;
step S13, focusing on the transmitting and receiving ends, the multipath maximum propagation distance c tau in the channel max Determining an ellipsoid corresponding to the multipath distribution range for the length of the long axis; the ellipsoid equation is:
[x',y',z'] T ＝R([x,y,z] T -[x 0 ,y 0 ,z 0 ] T ),
wherein a, b and c are respectively a long half shaft, a middle half shaft and a short half shaft of an ellipsoid, and the expression is
[x 0 ,y 0 ,z 0 ]Is the coordinate of the central point of the ellipsoid in the geodetic coordinate system, [ x, y, z [ ]]Is the coordinate in the earth coordinate system, [ x ', y ', z ']Is a coordinate under a Cartesian coordinate system which is constructed by taking the direction of the long axis, the middle axis and the short axis of the ellipsoid as the direction of a coordinate axis, and is used for judging whether the ellipsoid is normal or normal>
A set of standard orthogonal bases under the ellipsoid coordinate system is adopted, and R is a rotation matrix transformed from the geodetic coordinate system to the ellipsoid coordinate system;
step S14, obtaining a scatterer distribution three-dimensional region V = V 1 ∪V 2 ；
Obtaining a three-dimensional region V corresponding to a scatterer distribution prior hypothesis 1 ；
The three-dimensional region V 1 A three-dimensional area which takes a coordinate range sigma on an XY plane as a bottom surface and takes the average height h of the building in the coordinate range sigma on the XY plane as a height; wherein, the coordinate range σ of the scatterer distributed on the XY plane is:
the V is 2 A cylindrical region defined by the size of the fuselage and distributed for the fuselage diffuser;
step S15, based on the moving speed of the transmitting and receiving end, the geometric position of the transmitting and receiving end and the scatterer, obtaining the prior probability distribution F (F) of the multipath Doppler frequency D )；
And S16, setting primary layering number, determining the Doppler frequency estimation range and the initial value of each layer, and then performing secondary layering to obtain the Doppler frequency estimation range and the initial value of each layer after secondary layering.
Preferably, the step S15 specifically includes: the prior probability distribution of the multipath Doppler frequency is obtained by adopting Monte Carlo simulation, and the specific process is as follows:
by randomly scattering points in a three-dimensional region V of a scatterer, only one scattering is considered in addition to a direct view component, the Doppler frequency of each scattering point is calculated, and a radio wave diffusion factor is used
And a value of [0,1]The attenuation factors are weighted and normalized to obtain the prior summary of the Doppler power spectrumRate distribution F (F) D,LOS )。
Preferably, the step S16 specifically includes:
determining a prior probability distribution F (F) of a multipath Doppler frequency D ) Probability distribution interval of [ f ] D,min ,f D,max ]Setting the first number of layers to be K 1 Then, the k-th layer of the first time layer 1 The range of the layer doppler frequency estimation is:
wherein, F -1 (. Is a prior probability distribution F (F) of the Doppler power spectrum D ) The inverse function of (c); f (F) D,LOS ) Is prior probability distribution of the normalized Doppler power spectrum;
considering the jitter factor of the airplane body, carrying out secondary layering, wherein the kth layering 2 The layer doppler frequency estimation range is:
wherein, K 2 Setting a secondary layering number; inf (F) D,max ) And sup (F) D,max ) Upper and lower bounds, respectively, of the maximum doppler estimation range;
and &>
Respectively being the kth in one delamination 1 The upper and lower bounds of the range of the layer doppler frequency estimates.
Preferably, the step S2 specifically includes:
the signal model adopted is
Wherein alpha is l ,f D,l ,τ l Respectively representing the complex amplitude, doppler frequency and time delay of the l-th multipath, N (t) and N l (t) is Gaussian white noise Process with 0-means, N 0 Is a normal number, beta, related to the noise power l L =1, 2.., L is a decomposition of the noise part;
the expectation maximization and coordinate state updating process specifically comprises the following steps:
first, the hidden data expectation of the I-th multipath signal is calculated
Wherein the content of the first and second substances,
is the l-th multipath signal reconstructed according to the initial value of the parameter>
The parameter value of the ith multipath signal in the current iteration process, L is the number of paths, and y (t) always receives the signal, so that the time delay, doppler frequency and complex amplitude of the ith multipath in the current iteration are respectively estimated as follows:
wherein, Λ (theta) l ；x l ) Is the likelihood function of the ith multipath component, θ l ＝[τ l ,f D,l ,α l ]Is the l-th multipath component x l Including the time delay tau l Doppler frequency f D,l And a complex amplitude a l (ii) a The superscripts i-1 and i represent the last and current iterations, respectively, H is the Hamiltonian, c (τ, f) D ) Is a steering vector in the combined domain of delay and doppler.
Preferably, the threshold value in step S3 is 2f D,sidelobe Wherein f is D,sidelobe Is the side lobe width of the doppler domain likelihood function.
Compared with the prior art, the invention has the following advantages:
the invention improves the traditional SAGE algorithm based on the prior information of the scatterer distribution in the air-to-ground channel, exchanges the estimation sequence of multipath by adopting a layered processing method, effectively reduces the influence of the periodic extension of the Bartlett likelihood function spectrum on the estimation process, and can more reliably estimate the Doppler value of the multipath component under the special scene of the air-to-ground channel when the length of a signal frame is limited and the Doppler estimation range is limited.
Drawings
FIG. 1 is a schematic flow diagram of the process of the present invention;
FIG. 2 is a schematic diagram of an air-to-ground channel propagation scenario after a directional antenna is used in the embodiment;
FIG. 3 is a schematic diagram of a second air-to-ground channel propagation scenario in which a directional antenna is employed in the embodiment;
FIG. 4 is a power delay profile in an embodiment;
FIG. 5 is a diagram illustrating a cumulative distribution function of prior Doppler spectra obtained by the Monte Carlo method in the embodiment;
FIG. 6 is a result of multipath Doppler estimation obtained by using SAGE algorithm in the present invention in the embodiment;
FIG. 7 is a multipath Doppler estimation result obtained by using a conventional SAGE algorithm in the embodiment;
fig. 8 is a graph of the doppler frequency theoretical value of the LOS path over time in an embodiment.
Detailed Description
The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the drawings in the embodiments of the present invention, and it is obvious that the described embodiments are some, not all, embodiments of the present invention. All other embodiments, which can be obtained by a person skilled in the art without any inventive step based on the embodiments of the present invention, shall fall within the scope of protection of the present invention.
The invention is based on that in an air-to-ground channel, when the height of an airborne end from the ground is larger, multipath components are mainly introduced by a scatterer near a transceiver, doppler frequency components in the channel are centrally distributed near main path Doppler frequency, which is prior information, the Doppler frequency estimation range is reasonably layered, multipath components closer to the main path Doppler frequency are preferentially estimated, and then the estimation range is continuously expanded, so that the influence of periodic extension in a likelihood spectrum on the estimation process is reduced, and the estimation result is closer to the actual condition.
Example 1
As shown in the second air-to-ground channel propagation scenario of fig. 3, based on the geometric position relationship between the transceiver end and scatterer distribution in this scenario, the doppler frequency of the multipath component due to the scatterer is distributed near the LOS (Line of Sight) path, and a probability distribution with a certain bias is shown approximately around the doppler frequency of the LOS path. In addition, in air-to-ground scenarios, directional tracking antennas are often employed to improve channel gain, in this case due to directional antenna radiationThe effect of the radio mode and the propagation distance of the multipath will be further contracted, taking the channel measurement in this experimental example as an example, the center frequency is 27.5GHz, the bandwidth is 420MHz, the sampling rate is 2.24GHz, and the interval duration of the adjacent frames in the snapshot is 1.4629 × 10 -4 s, the doppler estimate range determined thereby is within ± 3418Hz, less than the maximum doppler frequency possible due to the maximum flight speed ± 5041.7Hz. The ground end adopts a directional tracking antenna with the beam width of only 0.4 degree, the airborne end is an omnidirectional antenna, the airplane horizontally cruises, the maximum flying speed is 55m/s, and the flying height is 2000 m. As can be seen from the power delay profile shown in fig. 4 (all snap-shot LOS paths have been time-delayed aligned), the propagation distance of multipath compared to LOS paths is within 15 meters, the scatterers will be distributed within the ellipsoid as shown in fig. 3, and the doppler frequency of multipath will also be very concentrated near LOS paths.
Aiming at the above situation, the specific implementation process of the algorithm in the invention is shown in fig. 1, the method reasonably assumes the probability density function of the Doppler power spectrum by using the prior information of the scattering body distribution in the channel, and sets the initial value and carries out layered estimation on the Doppler parameters of the multi-path in the SAGE algorithm based on the prior probability distribution, preferentially estimates the multi-path component in the interval in which the Doppler frequency is most likely to appear, and comprises the following steps:
s1, according to the values of multipath time delay in the power time delay spectrum, the ground environment characteristics during measurement and the size of an airplane, the size of an ellipsoid is set to be that the long axis is 15 meters longer than the focal length, the height of a ground scatterer is below 20 meters, and scattering points of a fuselage are distributed in a cylinder which takes a transmitting antenna as a center, is 3 meters high and has a radius of 15 meters. The cumulative distribution function of the prior doppler power spectrum of the partial snapshots shown in fig. 5 can be obtained by monte carlo simulation. It can be known that the prior doppler distribution range in the cruise scene is below 300Hz and is much smaller than ± 3418Hz determined by the reciprocal of the frame interval, so that a second expansion method is not needed;
after analysis, a first method is adopted to determine the initial estimation range of the doppler frequency according to the relative motion speed (recorded in a flight log) of the transmitting and receiving end in the current snapshot and the time interval information of the adjacent frames, and the method specifically comprises the following steps:
step S11, obtaining an initial Doppler frequency estimation range:
where c is the speed of light, f c In order to measure the center frequency of the system,
for the relative radial movement speed of the transmitting and receiving ends, T f For measuring adjacent frame time intervals, v, of signals max The maximum relative movement speed of the transceiving end in the current snapshot measurement process is obtained;
step S12, determining the range of the theoretical Doppler frequency:
wherein v is max The maximum relative movement speed of the receiving and transmitting end in the current snapshot measurement process, c is the light speed, f c Is the center frequency of the measuring system;
step S13, determining the Doppler frequency estimation range of each layer:
the widened doppler estimation range in each layer estimation is as follows:
wherein, the first and the second end of the pipe are connected with each other,
and &>
The difference between the upper and lower limits of the Doppler frequency initial estimation range and the maximum Doppler estimation range is respectively;
the k-th layer doppler frequency estimation range is:
wherein K is a set layering number;
step S14, setting the initial value grid point interval of the multipath on the Doppler domain as
Wherein, f D,sidelobe Is the side lobe width of the Doppler domain Bartlett likelihood spectrum and is expressed as->
T s Duration of a snapshot; setting the initial possible value of the current multipath Doppler frequency, wherein the initial value of the current multipath Doppler is determined according to the position of the highest peak in the Doppler spectrum;
the initial possible value of the multipath delay can be determined according to the highest peak position I in the Power Delay Profile (PDP) 0 Determining that the cell interval [ I ] is divided around the peak time delay 0 -1/f s ,I 0 +1/f s ]And the grid point is used as the possible values of the initial estimation range and the current multipath time delay, wherein f s The sampling rate of the receiving end;
the method of step S1 is suitable for the case where the multipath doppler frequency spread is very limited, including the scenario where a narrow beam directional antenna is used;
s2, utilizing the initial Doppler frequency estimation range, adopting expectation maximization to the current estimation path, and updating and acquiring the parameter estimation of the current multipath component based on the coordinate state, wherein the specific process is as follows:
the signal model used was:
wherein alpha is l ,f D,l ,τ l Respectively representing the complex amplitude, doppler frequency and time delay of the first multipath, N (t) and N l (t) is 0-mean Gaussian white noise Process, N 0 Is a normal number, beta, related to the noise power l L =1, 2.., L is a decomposition of the noise part;
the expectation maximization and coordinate state updating process specifically comprises the following steps:
first, the hidden data expectation of the I-th multipath signal is calculated
Wherein the content of the first and second substances,
is the l-th multipath signal reconstructed according to the initial value of the parameter>
The parameter value of the ith multipath signal in the current iteration process, L is the number of paths, and y (t) always receives the signal, so that the time delay, doppler frequency and complex amplitude of the ith multipath in the current iteration are respectively estimated as follows:
wherein, Λ (theta) l ；x l ) Is the likelihood function of the ith multipath component, θ l ＝[τ l ,f D,l ,α l ]Is the l-th multipath component x l Including the time delay tau l Doppler frequency f D,l And a complex amplitude a l (ii) a The superscripts i-1 and i denote the last and current iteration, respectively, H is the Hamiltonian, c (τ, f) D ) A steering vector on a time delay and Doppler joint domain;
in the above process, the Bartlett likelihood spectrum used is:
wherein, c (τ, f) D ) As steering vectors in the combined delay and doppler domain,
expectation of hidden data for the ith multipath signal>
The covariance matrix of (a);
s3, judging whether the boundary of the Doppler estimation result of the current multipath and the Doppler estimation range of the current layer is smaller than a set threshold value 2f D,sidelobe (ii) a If yes, and the current Doppler estimation range does not reach the maximum Doppler estimation range, expanding a layer of Doppler estimation range, and repeating the step S2; otherwise, ending the parameter estimation of the current path;
s4, estimating parameters of the next multipath component, and repeating the steps S2 to S3 until all paths are estimated, wherein the number L of the paths is set to be 15;
s5, comparing the difference value of the likelihood function values between two iterations with a set threshold (set to be 0.005 here), judging whether SAGE iterative estimation is converged, if not, repeating the steps to perform next iterative updating, otherwise, finishing the estimation process;
the likelihood function Λ (θ; y) is expressed as:
wherein, Λ (theta; y) is a likelihood function related to the parameter set theta to be estimated and the value of the received signal y (t), and N is 0 Is the white noise power in the current environment, D 0 For the acquisition time window of the current snapshot of data,
for the operation of the real part, s (t; theta) is a hidden data part excluding white Gaussian noise in the received signal at time t, and theta l The parameter to be estimated for the ith multipath component includes the complex amplitude α l Time delay tau l And a Doppler frequency v l 。
The multipath doppler frequency estimation result obtained after the above process is implemented is shown in fig. 6. In contrast, FIG. 7 shows the multipath Doppler estimates obtained using the traditional SAGE method, and FIG. 8 shows the LOS path Doppler frequency calculated from the GPS position, velocity, heading of the aircraft and the GPS position of the ground station recorded in the flight log. It can be seen that, when the conventional SAGE algorithm is used for estimation, the doppler estimation result is erroneous due to the influence of the doppler domain period prolongation, and the improved SAGE algorithm described in this embodiment can be used for more reliably estimating the doppler frequency of the multipath in the channel.
The embodiment provides an SAGE method when the Doppler estimation range is limited under the condition of high speed in an air-to-ground channel, reasonable assumption is carried out on a probability density function of a Doppler power spectrum by utilizing prior information of scatterer distribution in the channel, initial value setting and layered estimation are carried out on Doppler parameters of multiple paths in an SAGE algorithm based on the prior probability distribution, and the multiple path components in an interval where Doppler frequency is most likely to appear are preferentially estimated. And exchanging the estimation sequence of the multipath, and reducing the influence of the period extension part in the likelihood spectrum on the estimation process so as to obtain a more real and reliable estimation result.
The method of this embodiment 1 is suitable for the case where the multipath component doppler is very concentrated, and in this case, most of the multipath component doppler frequency falls within the range where aliasing does not occur, and can be centered on the direct-view doppler frequency.
Example 2
As shown in fig. 2, when the airborne terminal is higher than the ground, except for the direct-view component, the multipath component in the received signal is mainly generated by scatterers near the ground terminal and corresponding to the airframe structure, and the scatterers are distributed in an ellipsoid whose size is determined by the longest transmission distance with the ground terminal and the airborne terminal as the focus.
In this embodiment, for the above situation, the step S1 adopts an expanding method, which includes the following sub-steps:
s11, constructing a geodetic coordinate system, and determining the positions of a ground end and a receiving end in the coordinate system:
constructing a geodetic coordinate system by taking one point on a map as a coordinate O point, taking a ground plane as an XY plane and taking a ground perpendicular line as a Z axis, wherein
A set of orthonormal bases under the geodetic coordinate system;
determining the geometric position of the ground end at the current moment in the geodetic coordinate system
The geometric position of the receiving end at the current moment is->
S12, acquiring the maximum multipath delay spread tau based on the channel impact response h (t, tau) and the power delay spectrum P (t, tau) of the current time max So as to obtain the maximum multipath propagation distance c tau in the channel max Wherein c is the speed of light;
the expression of the channel impulse response h (t, τ) and the power delay spectrum P (t, τ) at the current moment is as follows:
h(t,τ)＝IFFT(Y(f)*conj(U(f))/P(U(f)))
P(t,τ)＝E[|h(t,τ)| 2 ]
wherein, IFFT (·) is inverse discrete fourier transform, Y (f) is frequency domain representation of received signal snapshot Y (t) at the current time, and U (f) is frequency domain representation of transmitted signal snapshot U (t) at the current time;
step S13, focusing on the transmitting and receiving end, the multipath maximum propagation distance c tau in the channel max Determining an ellipsoid corresponding to the multi-path distribution for the length of the long axis; the ellipsoid equation is:
[x',y',z'] T ＝R([x,y,z] T -[x 0 ,y 0 ,z 0 ] T ),
wherein a, b and c are respectively a long half shaft, a middle half shaft and a short half shaft of an ellipsoid, and the expression is
[x 0 ,y 0 ,z 0 ]Is the coordinate of the central point of the ellipsoid in the geodetic coordinate system, [ x, y, z [ ]]In the geodetic coordinate systemOf [ x ', y ', z ']A coordinate system which is constructed by taking the direction of the long axis, the middle axis and the short axis of the ellipsoid as the coordinate axis direction, and a device for collecting and combining the coordinates under the Cartesian coordinate system>
A set of standard orthogonal bases under the ellipsoid coordinate system is adopted, and R is a rotation matrix transformed from the geodetic coordinate system to the ellipsoid coordinate system;
step S14, obtaining a scatterer distribution three-dimensional region V = V 1 ∪V 2 ；
It is assumed that scatterers are uniformly distributed within the three-dimensional region V and that the ground scatterer velocity is 0.
Obtaining a three-dimensional region V corresponding to a scatterer distribution prior hypothesis 1 ；
The three-dimensional region V 1 A three-dimensional area with the coordinate range sigma on the XY plane as the bottom surface and the average height h of the building in the coordinate range sigma on the XY plane as the height; regardless of the ground having large undulations, it can be assumed that the ground plane is z =0, and the coordinate range σ of the scatterer distributed on the XY plane is:
the V is 2 A cylindrical region defined by the size of the fuselage and distributed for the fuselage diffuser;
step S15, based on the mobile speed of the transmitting and receiving end, the geometric position of the transmitting and receiving end and the scatterer, the prior probability distribution F (F) of the multipath Doppler frequency is obtained D ) The specific process is as follows: by randomly scattering points in a three-dimensional region V distributed on a scatterer, only one scattering is considered in addition to a direct view component, the Doppler frequency of each scattering point is calculated, and a radio wave diffusion factor is used
And a value of [0,1]Weighting attenuation factors (considering the attenuation of the radio wave by the dielectric characteristics of scatterers in the propagation of reflection, scattering and the like), and normalizing to obtain a Doppler power spectrumPrior probability distribution F (F) D,LOS )。
Step S16, setting the first-level layering number, determining the Doppler frequency estimation range and the initial value of each layer, and then performing secondary layering to obtain the Doppler frequency estimation range and the initial value of each layer after secondary layering, wherein the specific process comprises the following steps:
determining a prior probability distribution F (F) of a multipath Doppler frequency D ) Probability distribution interval of [ f ] D,min ,f D,max ]Setting the first number of layers to be K 1 Then, the k-th layer of the first time layer 1 The range of the layer doppler frequency estimation is:
wherein, F -1 (. Cndot.) is the prior probability distribution F (F) of the Doppler power spectrum D ) The inverse function of (c); f (F) D,LOS ) Is prior probability distribution of the normalized Doppler power spectrum;
considering factors such as fuselage jitter, there may be multipath components with Doppler values greater than the prior Doppler frequency range, so that secondary layering is performed, wherein the kth layer in the secondary layering 2 The range of the layer doppler frequency estimation is:
wherein, K 2 Setting a secondary layering number; inf (F) D,max ) And sup (F) D,max ) Upper and lower bounds, respectively, of the maximum doppler estimation range;
and &>
Respectively being the kth in one delamination 1 The upper and lower bounds of the layer doppler frequency estimation range.
The parameter setting is subject to the actual application scenario, and the rest steps are the same as those in embodiment 1.
The method of this embodiment 2 is suitable for the case that the multipath doppler frequency distribution range is larger than the estimable range determined by the reciprocal of the frame interval duration, and needs to use the prior distribution as guidance.
While the invention has been described with reference to specific embodiments, it will be understood by those skilled in the art that various changes in form and details may be made therein without departing from the spirit and scope of the invention. Therefore, the protection scope of the present invention shall be subject to the protection scope of the claims.
Claims (9)
1. A Doppler parameter estimation method based on improved SAGE is characterized in that the method reasonably assumes a probability density function of a Doppler power spectrum by using prior information of scatterer distribution in a channel, performs initial value setting and layered estimation on a multipath Doppler parameter in a space alternative generalized expectation maximization SAGE algorithm based on prior probability distribution of Doppler frequency, and preferentially estimates multipath components in an interval with high prior probability of Doppler frequency;
s1, determining initial values and estimation ranges of relevant parameters except Doppler frequency, and giving prior distribution of multipath Doppler frequency; setting the number of layers, and determining the Doppler frequency estimation range and the initial value of each layer;
s2, utilizing an initial Doppler frequency estimation range, adopting expectation maximization on a current estimation path, and updating and acquiring parameter estimation of a current multipath component based on a coordinate state;
s3, judging whether the Doppler frequency estimation result of the current multipath is smaller than a set threshold value from the boundary of the Doppler frequency estimation range of the current layer: if yes, and the current Doppler frequency estimation range does not reach the maximum Doppler estimation range, expanding a layer of Doppler frequency estimation range, and repeating the step S2; otherwise, finishing the parameter estimation of the current path, and turning to the step S4;
s4, estimating parameters of the next multipath component, and repeating the steps S2 to S3 until all paths are estimated;
s5, comparing a difference value of likelihood function values between two iterations with a set threshold value, and judging whether SAGE iteration estimation is converged; if not, repeating S2-S4 to carry out the next iteration updating, otherwise, ending the estimation process.
2. The improved SAGE-based Doppler parameter estimation method according to claim 1, wherein the initial values and estimation ranges of the relevant parameters in the step S1 comprise multipath time delay, complex amplitude and multipath quantity outside multipath Doppler frequency.
3. The improved SAGE based Doppler parameter estimation method according to claim 2, wherein the step S1 comprises the following sub-steps:
step S11, obtaining an initial Doppler frequency estimation range:
wherein c is the speed of light, f c In order to measure the center frequency of the system,
for the relative radial movement speed of the transmitting and receiving ends, T f For measuring adjacent frame time intervals, v, of signals max The maximum relative movement speed of the transceiving end in the current snapshot measurement process is obtained;
step S12, determining a range of theoretical doppler frequencies:
wherein v is max The maximum relative movement speed of the receiving and transmitting end in the current snapshot measurement process, c is the speed of light, f c Is the center frequency of the measurement system;
step S13, determining the Doppler frequency estimation range of each layer:
the widened doppler estimation range in each layer estimation is as follows:
wherein, the first and the second end of the pipe are connected with each other,
and &>
Respectively an upper bound and a lower bound of an initial Doppler frequency estimation range; inf (F) D,max ) And sup (F) D,max ) Upper and lower bounds, respectively, of a theoretical doppler frequency range; k is a set layering number;
the k-th layer doppler frequency estimation range is:
wherein v is max The maximum relative movement speed of the receiving and transmitting end in the current snapshot measurement process, c is the speed of light, f c Is the center frequency of the measurement system;
and &>
Upper and lower bounds, f, respectively, of the initial Doppler frequency estimation range D,layer Multiple broadening in estimation for each layerA Doppler estimation range;
step S14, setting the initial value grid point interval of the multipath on the Doppler domain as:
wherein the content of the first and second substances,
and &>
Upper and lower bounds, f, of the k-th layer Doppler frequency estimation range, respectively D,sidelobe Is the side lobe width of the doppler domain likelihood function.
4. The improved SAGE-based Doppler parameter estimation method according to claim 3, wherein the Doppler domain likelihood function in the step S14 adopts Bartlett likelihood spectrum, and the side lobe width f of the Bartlett likelihood spectrum is D,sidelobe The expression is as follows:
wherein, T s Is the duration of a snapshot.
5. The improved SAGE based Doppler parameter estimation method as claimed in claim 1, wherein the step S1 is as follows:
s11, constructing a geodetic coordinate system, and determining the positions of a ground end and a receiving end in the coordinate system:
constructing a geodetic coordinate system by taking one point on a map as a coordinate O point, taking a ground plane as an XY plane and taking a ground perpendicular line as a Z axis, wherein
Is the earth groundA set of orthonormal bases under a coordinate system;
determining the geometric position of the ground end at the current moment in the geodetic coordinate system
The geometric position of the receiving terminal at the current moment is->
S12, acquiring the maximum multipath delay spread tau based on the channel impact response h (t, tau) and the power delay spectrum P (t, tau) of the current time max So as to obtain the maximum propagation distance of multipath in the channel as c * τ max Wherein c is the speed of light;
the expression of the channel impulse response h (t, τ) and the power delay spectrum P (t, τ) at the current time is as follows:
h(t,τ)＝IFFT(Y(f)*conj(U(f))/P(U(f)))
P(t,τ)＝E[|h(t,τ)| 2 ]
wherein, IFFT (·) is inverse discrete fourier transform, Y (f) is frequency domain representation of received signal snapshot Y (t) at the current time, and U (f) is frequency domain representation of transmitted signal snapshot U (t) at the current time;
step S13, taking the transmitting and receiving end as the focus, the multipath maximum propagation distance c in the channel * τ max Determining an ellipsoid corresponding to the multipath distribution range for the length of the long axis; the ellipsoid equation is:
[x',y',z'] T ＝R([x,y,z] T -[x 0 ,y 0 ,z 0 ] T ),
wherein a, b and c are respectively a long half shaft, a middle half shaft and a short half shaft of an ellipsoid, and the expression is
c = b, | | | is norm; [ x ] of 0 ,y 0 ,z 0 ]Is the coordinate of the central point of the ellipsoid under the geodetic coordinate system, [ x, y, z ]]Is the coordinate in the earth coordinate system, [ x ', y ', z ']A coordinate system which is constructed by taking the direction of the long axis, the middle axis and the short axis of the ellipsoid as the coordinate axis direction, and a device for collecting and combining the coordinates under the Cartesian coordinate system>
A set of standard orthogonal bases under an ellipsoid coordinate system, wherein R is a rotation matrix transformed from the geodetic coordinate system to the ellipsoid coordinate system;
step S14, obtaining a scatterer distribution three-dimensional region V = V 1 ∪V 2 ；
Obtaining a three-dimensional region V corresponding to a scatterer distribution prior hypothesis 1 ；
The three-dimensional region V 1 A three-dimensional area with the coordinate range sigma on the XY plane as the bottom surface and the average height h of the building in the coordinate range sigma on the XY plane as the height; wherein, the coordinate range σ of the scatterer distributed on the XY plane is:
the V is 2 A cylindrical region distributed for the fuselage scatterers and determined by the size of the fuselage;
step S15, based on the moving speed of the transmitting and receiving end, the geometric position of the transmitting and receiving end and the scatterer, obtaining the prior probability distribution F (F) of the multipath Doppler frequency D )；
And S16, setting primary layering number, determining the Doppler frequency estimation range and the initial value of each layer, and then performing secondary layering to obtain the Doppler frequency estimation range and the initial value of each layer after secondary layering.
6. The improved SAGE-based Doppler parameter estimation method according to claim 5, wherein the step S15 is specifically as follows: the prior probability distribution of the multipath Doppler frequency is obtained by adopting Monte Carlo simulation, and the specific process is as follows:
by randomly scattering points in a three-dimensional region V distributed on a scatterer, only one scattering is considered in addition to a direct view component, the Doppler frequency of each scattering point is calculated, and a radio wave diffusion factor is used
And a value of [0,1]Weighting the attenuation factors in the Doppler power spectrum, and obtaining prior probability distribution F (F) of the Doppler power spectrum after normalization D )。
7. The method of claim 6, wherein the step S16 is specifically as follows:
determining a prior probability distribution F (F) of multipath Doppler frequencies D ) Probability distribution interval of [ f ] D,min ,f D,max ]Setting the first number of layers to be K 1 Then the k-th layer of the first layer 1 The layer doppler frequency estimation range is:
wherein, F -1 (. Is a prior probability distribution F (F) of the Doppler power spectrum D ) The inverse function of (a);
considering the jitter factor of the airplane body, carrying out secondary layering, wherein the kth layering 2 The range of the layer doppler frequency estimation is:
wherein, K 2 Setting a secondary layering number; inf (F) D,max ) And sup (F) D,max ) Upper and lower bounds, respectively, of a theoretical doppler frequency range;
and &>
Respectively being the kth in one delamination 1 The upper and lower bounds of the layer doppler frequency estimation range.
8. The improved SAGE-based Doppler parameter estimation method according to claim 1, wherein the step S2 specifically comprises:
the signal model adopted is
Wherein alpha is l ,f D,l ,τ l Respectively representing the complex amplitude, doppler frequency and time delay of the l-th multipath, N (t) and N l (t) is 0-mean Gaussian white noise Process, N 0 Is a normal number, beta, related to the noise power l L =1,2, L is the noise partA decomposition of (2);
the expectation maximization and coordinate state updating process specifically comprises the following steps:
first, the hidden data expectation of the I-th multipath signal is calculated
Wherein the content of the first and second substances,
is the l-th multipath signal reconstructed according to the initial value of the parameter>
The parameter value of the ith multipath signal in the current iteration process is obtained, L is the number of paths, and under the condition of always receiving the signal y (t), the estimation of the delay, the Doppler frequency and the complex amplitude of the ith multipath in the current iteration is respectively as follows:
wherein, Λ (theta) l ；x l ) Is the likelihood function of the ith multipath component, θ l ＝[τ l ,f D,l ,α l ]Is the l-th multipath component x l Including the time delay tau l Doppler frequency f D,l And a complex amplitude a l (ii) a The superscripts i-1 and i denote the last and current iteration, respectively, H is the Hamiltonian, c (τ, f) D ) Is a steering vector in the joint domain of delay and doppler.
9. The improved SAGE-based Doppler parameter estimation method as claimed in claim 8, wherein the threshold value in the step S3 is 2f D,sidelobe Wherein f is D,sidelobe Is the side lobe width of the doppler domain likelihood function.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111474809.4A CN114362852B (en) | 2021-12-03 | 2021-12-03 | Doppler parameter estimation method based on improved SAGE |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111474809.4A CN114362852B (en) | 2021-12-03 | 2021-12-03 | Doppler parameter estimation method based on improved SAGE |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114362852A CN114362852A (en) | 2022-04-15 |
| CN114362852B true CN114362852B (en) | 2023-03-28 |
Family
ID=81097048
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202111474809.4A Active CN114362852B (en) | 2021-12-03 | 2021-12-03 | Doppler parameter estimation method based on improved SAGE |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114362852B (en) |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105978833A (en) * | 2016-06-29 | 2016-09-28 | 西安电子科技大学 | Improved SAGE channel parameter estimation method |
| CN106713191A (en) * | 2017-02-28 | 2017-05-24 | 西安电子科技大学 | Multistage searching SAGE method |
| WO2019138156A1 (en) * | 2018-01-12 | 2019-07-18 | Nokia Technologies Oy | Profiled channel impulse response for accurate multipath parameter estimation |
| CN113225274A (en) * | 2021-04-14 | 2021-08-06 | 西安宇飞电子技术有限公司 | Multi-path channel model measuring method for fast moving |
-
2021
- 2021-12-03 CN CN202111474809.4A patent/CN114362852B/en active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105978833A (en) * | 2016-06-29 | 2016-09-28 | 西安电子科技大学 | Improved SAGE channel parameter estimation method |
| CN106713191A (en) * | 2017-02-28 | 2017-05-24 | 西安电子科技大学 | Multistage searching SAGE method |
| WO2019138156A1 (en) * | 2018-01-12 | 2019-07-18 | Nokia Technologies Oy | Profiled channel impulse response for accurate multipath parameter estimation |
| CN113225274A (en) * | 2021-04-14 | 2021-08-06 | 西安宇飞电子技术有限公司 | Multi-path channel model measuring method for fast moving |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114362852A (en) | 2022-04-15 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US9084217B2 (en) | Single-site localization via multipath fingerprinting | |
| CN110412559B (en) | Non-coherent fusion target detection method for MIMO radar of distributed unmanned aerial vehicle | |
| Honkavirta et al. | A comparative survey of WLAN location fingerprinting methods | |
| US8478294B2 (en) | Method and system for mobile station location | |
| CN103744076B (en) | MIMO radar moving target detection method based on non-convex optimization | |
| CN109765521A (en) | A kind of Beam Domain imaging method based on Subarray partition | |
| Khan et al. | Enhanced hybrid positioning in wireless networks I: AoA-ToA | |
| CN114371445A (en) | Multi-radiation source direct positioning method based on single unmanned aerial vehicle | |
| Haneda et al. | Radio propagation modeling methods and tools | |
| CN115150744A (en) | Indoor signal interference source positioning method for large conference venue | |
| Foroozan et al. | Direction finding algorithms for time reversal MIMO radars | |
| Yajnanarayana et al. | Multistatic Sensing of Passive Targets Using 6G Cellular Infrastructure | |
| Ding et al. | A time-varying transition channel model for air-ground communication | |
| CN114362852B (en) | Doppler parameter estimation method based on improved SAGE | |
| WO2022166477A1 (en) | Positioning method and apparatus, base station, computer device, and storage medium | |
| CN115052246A (en) | Broadband signal direct positioning method based on multi-frequency cost function fusion under unknown attenuation coefficient | |
| Tran et al. | Complexity reduction for hybrid toa/aoa localization in uav-assisted wsns | |
| Wan et al. | A range-Doppler-angle estimation method for passive bistatic radar | |
| US10085120B1 (en) | Copy aided geolocation | |
| Wang et al. | A comparison of outdoor-to-indoor wideband propagation at S-band and C-band for ranging | |
| Yassin et al. | 3D localization and mapping using mm-Wave: What are the opportunities in vehicular and indoor environments? | |
| Kuchar et al. | Azimuth, elevation, and delay of signals at mobile station site | |
| Monnoyer et al. | Grid Hopping: Accelerating Direct Estimation Algorithms for Multistatic FMCW Radar | |
| Riback et al. | Statistical analysis of measured radio channels for future generation mobile communication systems | |
| Pajusco | Double characterisations of power angular spectrum in macro-cell environment |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |