Wavefront Correction of Ionospherically Propagated HF Radio Waves Using Covariance Matching Techniques

High Frequency (HF) radio waves propagating in the ionospheric random inhomogeneous media exhibit a spatial nonlinearity wavefront which may limit the performance of conventional high-resolution methods for HF sky wave radar systems. In this paper, the spatial correlation function of wavefront is theoretically derived on condition that the radio waves propagate through the ionospheric structure containing irregularities. With this function, the influence of wavefront distortions on the array covariance matrix can be quantitatively described with the spatial coherence matrix, which is characterized with the coherence loss parameter. Therefore, the problem of wavefront correction is recast as the determination of coherence loss parameter and this is solved by the covariance matching (CM) technique. The effectiveness of the proposed method is evaluated both by the simulated and real radar data. It is shown numerically that an improved direction of arrival (DOA) estimation performance can be achieved with the corrected array covariance matrix.


Introduction
The use of High Frequency (HF) sky wave propagation as a means of target detection and tracking over the horizon has been investigated for many years [1], [2].It has showed powerful strength in large-scale surveillance capacities at an economic cost.In addition to this, it can be also adopted in the applications of ocean state estimation such as wind, ocean waves and ocean currents [3].
However, as for the ionospherically propagated HF radio wave, the contamination suffered during the ionospheric transit usually influences the accurate determination of target as well as the ocean state information extraction [4], [5].One type of the contamination results from the fine-scale irregular ionospheric structure which consists of irregularities of various electron density fluctuations.Such irregular structure can manifest in its distortions to the radio wavefronts of signals.Borne out by experiments [6], the distorted wavefronts exhibit progressive decorrelation between the receiving array sensors as the spatial separation distance increases.Coherence loss may directly result in the degradation of the performance of eigenstructure based high resolution techniques [7].This is because, as for the perfect wavefront situation, each wavefront contributes a rank one component to the array covariance matrix.When taking the wavefront distortions into consideration, the array covariance matrix will be tapered by another random matrix in the manipulation of Schur-Hadamard product (element-wise multiplication operator) [8].Such modulation can lead to the increase of the effective rank and cause the eigensubspace leakage problem [9].As for the conventional MUltiple SIgnal Classification (MUSIC) algorithm, the leaked signal subspace can give rise to the smearing in the direction of arrival (DOA), i.e. the wavefront appears to arrive from a spread of angles centered at the true direction.
To enhance the performance of MUSIC algorithm in DOA estimation for HF sky wave radar systems, special attention has to be paid to the role of spatial coherence in array processing.In this paper, an experimental investigation is first undertaken to demonstrate the breakdown of MUSIC algorithm when wavefronts have reduced spatial coherence.Afterwards, a model characterizing the spatial coherence loss after propagating through the ionosphere containing irregularities is theoretically constructed.Combined with the CM technique [10], [11], the restored array covariance matrix can be obtained and thus results in a noticeable performance improvement of the conventional MUSIC algorithm.
The remainder of the paper is organized as follows.In Sec. 2, the spatial correlation characteristics are studied with a large aperture antenna array.Afterwards, the distortions on the wavefront are modelled and corrected in Sec. 3. To verify the accuracy of the proposed model and the effectiveness of wavefront correction method, experiments on both simulated and real radar data are undertaken in Sec. 4. Eventually, a brief summary is given in Sec. 5.

Notation: (•)
T and (•) H are the transpose and Hermitian transpose operation, * is the complex conjugate, • indicates the Schur-Hadamard product, ∅ is the inverse of Schur-Hadamard product, E {•} is the expectation operator, • F is the Frobenius norm, tr {•} is the sum of diagonal elements.

Problem Formulation from an Experimental Study on the Signal Wavefront
In order to further study the wavefront signatures of ionospherically propagated radio waves, measurements of the beacon signal transmission received on a large aperture antenna array have been made.A brief description of the experimental arrangement is given below.
The beacon transmitter is located on an island with the baseline distance of 890 km far from the receiving station.The transmitted signal impinges on a five element array from the array normal direction where the elements are 32.8 m apart.The experiment was conducted at 11:00 on 26 December 2014 with the operating frequency of 15.8 MHz.The geometric layout of the experimental system is illustrated in Fig. 1.To track the specific propagation path of each ray from a given radar and beam direction, a 2D ray tracing technique is used based on work in [12] and integrated using an adaptive Runge-Kutta numerical method.Ionospheric profiles are generated by the latest International Reference Ionosphere (IRI-2011) [13].Figure 2 illustrates the ray-tracing result where rays are assumed to be launched over an elevation range from 5 • to 55 • at the selected operating frequency.As observed, the reflection happens at the height of 192 km in the ionospheric F region.The irregularities in this area can be regarded as density structures within the ionospheric plasma that are small compared to the scale size of the overall ionosphere.Propagating through such fine irregular structure, the correlation between closely-spaced ray paths will be degraded.To illustrate this, Fig. 3(a) gives the measured results of correlation coefficient for the signal amplitude, phase and combined amplitude-phase with respect to the first antenna element.Herein, three typical correlation associations are concluded from a number of experimental observations.According to Fig. 3(a), types 1 through 3 illustrate the intrinsic dependence of signal correlation against the amplitude and phase.With regard to type 1 and 2, it is seen that the correlation degradation of signal amplitude has a little influence on the array sampled complex signal (with both amplitude and phase).However, as observed in type 3, significant spatial coherence loss between each element can be easily spotted especially in the presence of phase correlation reduction.Therefore, we could conclude that the signal correlation has a stronger dependence on the phase correlation than the amplitude correlation.Combined with the corresponding MUSIC based DOA estimation result (128 snapshots are used) illustrated in Fig. 3(b), it is straightforward to notice that the spatial correlation reduction directly results in the DOA performance degradation.To explain the result in Fig. 3(b), we consider q number of narrowband stationary zero mean mutually uncorrelated plane waves impinging on a uniform linear array (ULA) of M elements and the array M × 1 vector observations can be modeled as where is the wavefront vector corresponding to the source direction θ l , l = 1, • • • , q, s l is the amplitude of the l-th source, n is the zero mean white noise process, d is the element interval and λ is the radar wavelength.
To involve the influence of wavefront distortions, a l is recast to consider the time varying wavefront fluctuation by introducing the spatial perturbation function and thus where the random amplitude and phase perturbation lumped into the wavefront throughout each channel.Consider that each source has the same spatial perturbation function and, in this case, the array matrix is given by where σ 2 is the noise variance, I is the identity matrix and F is defined as the spatial coherence matrix characterizing the coherence loss extent due to the properties of propagation medium which is primarily determined by the spatial correlation function.In fact, (3) is a relatively simply model derived on condition that each source is modulated by one same random process and much more complicated circumstances are discussed in the Appendix.
According to (3), the issue of angle spreading around the true direction illustrated in Fig. 3(b) can find its explanation through the spatial coherence matrix modulation.Therefore, to yield the method for wavefront distortion correction, a model is needed to explain the intrinsic perturbation mechanism.

Modeling and Correction of Wavefront Distortion
As alluded to earlier, the correlation of signals sampled at the wavefront has a much stronger dependency on its phase correlation.Therefore, in what follows, the spatial correlation function of perturbed phases for ionospherically propagated HF radio wavefronts will be theoretically constructed and then a CM technique based method is derived to eliminate the influence on the array covariance matrix induced by the wavefront distortion.

A Spatial Coherence Model for Imperfect Wavefronts
In general, the ionosphere is regarded as the inhomogeneous and non-stationary medium.The ionospheric irregularity component captured as the density structure within the ionospheric plasma will impart a random fluctuation on the wavefront and thus reduce the correlation between closelyspaced ray paths.To study this effect, the spatial correlation function of ionospheric perturbed phases is derived in the following part.
To facilitate the analysis and calculation, we ignore the influence of geomagnetic field and collisions and thus yield the basic dispersion relationship of electromagnetic waves expressed by ( 4) combined with the Appleton-Hartree formula [14] where N is the index of refraction.ω p , defined as ω p = e 2 n/(ε 0 m 0 ), is electron plasma frequency, e is the charge on an electron, m 0 is the electron mass, n is the total ionosphere plasma density profile, ε 0 is the permittivity of free space, c is the velocity of light and k is radar wave number.Afterwards, we consider a Cartesian system of coordinates, with x east, y north and z vertical.In the ionosphere, n is a function of space which is assumed as linearly increased with altitude z independent of the horizontal coordinate (x, y) and thus Further, the phase accumulated along a path to a single scattering event on the ground can be estimated by Following [15], n is reckoned as the combination of the quiescent background part n 0 and an irregular part n 1 , i.e. n = n 0 + n 1 .Combined with such relationship, we use the Taylor series expansion of k z (n) at n = n 0 up to the first order term and we have Substituting ( 7) into (6), the accumulated phase involving the irregularity influence can be calculated by where r e = e 2 / 4πε 0 m 0 c 2 is the classical electron radius (≈ 2.8 × 10 −15 m), z t is the height corresponding to the reflection and z t = 0 indicates the bottom of the ionosphere.
Based upon (8), the spatial autocorrelation function over the horizon plane (x, y) is estimated by the ensemble average of To estimate R ϕ 1 , we let u = z − z /z t and v = z + z /z t so that the Jacobian determinant is Therefore, ( 9) can be reduced into (10) with where µ (•) is the step function.Using the Fourier transform implementation of (10), the perturbed phase spectrum is yielded by Invoking ( 10) and (11), it is seen that R ϕ 1 can be determined when S n 1 is specified.Herein, the fourth order power law spectrum in [16] is exploited so that R ϕ 1 is finally calculated by taking the inverse Fourier transform implementation of (11).In lieu of a direct measurement of correlation, the complex signal form containing the perturbed phases is considered to measure the correlation strength  13), (14) and the average measurement of 500 independent statistical tests of the real radar data.The blue error bars indicate the deviation range of the measured correlation.
Substituting the calculation result of R ϕ 1 into (12), we have where ϕ 2 1 represents the mean-square phase fluctuation, κ 0 is the ionosphere outer scale parameter and ρ c indicates the correlation distance in space.
Resorting to (13), it is a physics-based model where the parameters primarily depend on the mean square phase fluctuation magnitude and the outer scale length of the plasma density irregularities.As for the de facto condition, it is difficult to get these parameters because they are generally obtained from the experimental observations from radar, in situ, and satellite.Therefore, a simplified model of ( 13) is needed.To this end, we consider a wavefront model on the basis of two assumptions: 1) the amplitude remains constant and the phase follows the independent and identically distributed zero mean Gaussian process with the phase variance σ 2 ϕ and 2) the correlation function subjects to a strictly decreasing of the separation distance.By doing this, the correlation between the i th and m th sensor is determined by and thus the spatial coherence matrix in (3) is reasonable to be calculated by where F is a real-valued symmetric Toeplitz positive definite matrix.To verify the rationality of the simplified model ( 14), Fig. 4 illustrates the model fitting result based on least square method between ( 13), ( 14) and the average measurement of 500 independent statistical tests of the real radar data.As noticed, the model of ( 14) is a decent tradeoff between the complexity and accuracy and can also agree well with the theoretical and measured results.

Wavefront Correction Using Covariance Matching Approach
Resorting to (3), it follows that the problem of wavefront correction is recast as demodulating the influence of the spatial coherence matrix.In order to address this problem, the CM technique [10], [11] is used.To explain this, we denote the corrected array covariance matrix as Rx = Rx ∅F.
The cost function is given by where Θ is the vector space of unknown parameters.The required estimates of parameters will be chosen by making the cost function (17) to the minimum in the least square fit sense and ( 17) is equivalently reduced into From (18), it is seen that the parameters are estimated through the multidimensional search over Θ.To reduce the search dimension, a dimension-reduced optimization target function of ( 18) is given by (19) (at the bottom of this page).

In (19), P
is the matrix of wavefront vectors corresponding to each source.By doing this, the finally obtained set of estimated parameters can be reduced from q 2 + q + 2 × 1 to q + 1 × 1.To determine the q + 1 × 1 parameters, particle swarm optimization (PSO) algorithm [17] is used to solve this multidimensional parameter estimation problem for its low computational effort and rapid global convergence.

Simulation and Experiment
In simulations, a ULA of 16 elements with half wavelength interelement spacing is used.Two conditions corresponding to a single source and two uncorrelated equipowered sources are considered with signal-to-noise SNR = 20 dB.The array covariance matrix is estimated with the number of 32 snapshots.The phase variance of the distorted wavefront is σ 2 ϕ = 1 ( σ 2 ϕ is also renamed as the coherence parameter whose magnitude indicates the extent of coherence loss) and the spatial coherence matrix F is constructed according to (15).According to the simulated results in Fig. 5, it is straightforward to notice that angle spreading and deviation around the true direction are the most significant effects triggered by the distorted wavefront.With the proposed method, the imperfect wavefront is well restored which can be verified by the noticeable performance improvement of MUSIC algorithm.Inspiringly, the present method can be further developed to suppress the covariance matrix tapering influence in 2D space time adaptive processing [9].
The performance of PSO method is illustrated in Fig. 6.For comparison, the genetic algorithm (GA) in [10] is considered.Table 1 and Tab. 2 give the parameters applied for these two algorithms.Figure 6(a) shows the distribution of target function (19) over the parameter searching range.It is seen that the target function has a global minimization value at the true DOA direction and coherence parameter.Different from [7], any priori knowledge of F is not required to be given and all parameters will be yielded simultaneously once (19) converges to the global minimum.According to Fig. 6(b), it illustrates the convergence performance between PSO and GA algorithms against the number of iterations.Comparatively, PSO algorithm can converge to the optimal Θ0 = arg min solution in much smaller iterative steps and thus appears to be more appropriate for the requirement of real-time processing.In addition, the fast convergence performance makes the proposed method reasonable to be applied for dealing with more complicated scenarios discussed in the Appendix.As shown in Fig. 7, the validity of the wavefront correction method is verified by the real radar data.Here, we assume that the measured result of correlation coefficient corresponds to Type 3 in Fig. 3(b).As observed from this case, the received wavefront signal has an evidently degraded correlation between each antenna element.With the proposed wavefront correction method, it is seen that the output of conventional MUSIC spectrum has been improved to the condition where only the amplitude distortion is left to be corrected on the wavefront.

Conclusions
In the present paper, special attention has been paid to correct the distorted wavefronts of ionospherically propagated HF radio signals so as to enhance the performance of conventional MUSIC algorithm.To study the influence on the array covariance matrix due to the imperfect wavefront, the spatial coherence matrix is constructed with a reasonably simplified spatial correlation function.On the basis of CM technique, the distorted array covariance matrix can be restored by eliminating the influence induced by the spatial decorrelation of the received wavefront signals.In the end, simulation and experimental radar data have demonstrated the effectiveness of the proposed method.
Future work will study the amplitude fluctuation mechanism on the wavefront as well as the corresponding technique to restore the array covariance matrix.

Fig. 2 .
Fig. 2. Ray tracing output of the experiment.The blue solid line indicates the ray trajectory travelling from the transmitter to the receiver.

Fig. 3 .
Fig. 3. Illustrations of typical correlation measurement result and the corresponding DOA estimation based on MUSIC algorithm.A correlation coefficient of 0.7 is assumed as the threshold indicating the correlation strength.

Fig. 4 .
Fig. 4. Spatial correlation model comparison between (13),(14) and the average measurement of 500 independent statistical tests of the real radar data.The blue error bars indicate the deviation range of the measured correlation.
A single source case with DOA of 0 • .Two uncorrelated sources with DOAs of −5 • and 5 • .

Fig. 5 .
Fig. 5. Simulated results of DOA estimation with corrected wavefronts based on the proposed method.

Fig. 7 .
Fig. 7. Experimental result of DOA estimation with corrected wavefronts based on the proposed method.