HF Radio Angle-of-Arrival Measurements and Ionosonde Positioning

Since 2010 a 2nd generation NOAA MF/HF radar, also referred to as the VIPIR ionosonde, has been operated at Hualien, Taiwan (23.8973°N, 121.5503°E). The Hualien VIPIR ionosonde is a modern ionospheric radar, fully digitizing complex signal records and using multiple parallel receiver channels for simultaneous signal measurements from multiple spaced receiving antennas. This paper considers radio direction finding based on interferometric phase measurements from a horizontal antenna array in the Hualien VIPIR ionosonde system. We applied the Hermite normal form method to solve the phasemeasurement aliasing and least squares problems and improve the radio angle-of-arrival (AOA) measurements. Backward ray-tracing simulation has been proposed to determine radio transmitter position. This paper presents a numerical, step by step ray-tracing method based on the IGRF superimposed onto a phenomenological ionospheric electron density model, the TaiWan Ionospheric Model (TWIM). The proposed methodology is successfully applied to locate two experimental HF radio transmitters at Longquan and Chungli with distance errors within 5 km and less than 5% of the great circle distances.


INTRODUCTION
Ground-based and vertical incidence ionospheric sounding radar (also termed ionospheric sounder or ionosonde) has been used for providing viable and inexpensive ionospheric measurements using the radio total reflection method.Ionospheric sounders usually transmit pulsed (and/ or coded) waveforms in the medium or high-frequency (MF/ HF) band and measure the virtual ranges, i.e., envelop group delay of ionospheric echoes.The virtual range of the ionospheric echoes as a function of the radio carrier frequency is the fundamental ionosonde data product and is named the ionogram.Many earlier investigations (Wilkinson 1998;Wright and Pitteway 1999;Tsai and Berkey 2000;Reinisch et al. 2005;Bullett et al. 2010) provided a means for autonomously characterizing ionograms into ionospheric parameters and reducing ionograms to their electron density (N e ) profiles.Synoptic ionospheric measurements derived from ionograms have provided valuable information on the physical characteristics of the ionosphere and for MF and HF radio communications and wave propagation predictions (e.g., URSI-88 coefficients) (Piggott and Rawer 1972).In 2010 the National Central University installed a new digital ionosonde system at Hualien, Taiwan (23.8973°N, 121.5503°E).The digital ionosonde was developed by Scion Associates, supported by the US Air Force Research Laboratory and their Small Business Innovative Research Program (Grubb et al. 2008;Bullett et al. 2010).This project was inspired by the capabilities of the original National Oceanic and Atmospheric Administration (NOAA) HF radar/dynasonde (Grubb 1979;Wright and Pitteway 1982).Bullett et al. (2010) named the digital ionosonde system Vertical Incidence Pulsed Ionospheric Radar (VIPIR).We name our system the Hualien VIPIR ionosonde in the following study.
As the original NOAA HF radar, the VIPIR ionosonde system was an evolution of the traditional pulsed ionosondes and includes new capabilities related to phase coherent measurements on multiple parallel receiver channels and spaced antennas.In ionospheric sounding applications the system can transmit sets of radio pulses at designated frequencies and measure the complex amplitudes of ionospheric echoes.This mode of operation is referred to as active.The VIPIR ionosonde system can also be operated in a passive model, where it does not transmit radio pulses, but the receiving antennas and channels still measure signal amplitudes and phases from any background signals.
The background signals can include those from stationary radio transmitters.These can be propagated either as ground waves or as single or multi-hop reflections from the ionosphere.It is well known that sky wave reception can be obtained over distances of a few thousand kilometers or more.The passive-mode observations can be used for solving background HF radio surveillance problems, notably transmitter localization based on radio angle-of-arrival (AOA) determination of sky waves and further ray tracing in the ionosphere.The general AOA determination problem considers receiving antennas with arbitrary locations and arbitrary directional characteristics in a noise/interference arbitrary covariance matrix environment.It has been described and solved using a direction finding system for example the minimum variance estimator (MVE) algorithm (Capon 1969;Carlson 1988) or the multiple signal classification (MUSIC) algorithm (Schmidt 1986).An excellent treatise by Gething (1991) reviewed the physical principles involved in radio-direction finding in the HF band.In this paper we consider only direction finding from a groundbased ionosonde system, especially the VIPIR ionosonde, and on signals propagated via the ionosphere.The problem can be simplified with a general plane wave approach at specified frequencies and radio AOA can be determined using interferometric phase measurements from a horizontal antenna array.Various authors have addressed the phase aliasing problem from interferometric measurements in the NOAA HF radar (Pitteway and Wright 1992;Tsai et al. 1997;Harris 2003).The aliasing problem could become thornier when the antenna array configuration and/or the sounding pulse pattern are different.HF radio ray tracing in the ionosphere mainly involves refraction and reflection due to the ionospheric structure on a scale that is large compared with the radio wavelength.Yeoman et al. (2001) used a modified version of the numerical ray-tracing code developed by Jones and Stephenson (1975) and evaluated ground range accuracy in the Super Dual Auroral Radar Network over-the-horizon HF radar systems.They found that the standard range-finding algorithm was accurate to within 16 km for direct ionospheric scattering and 60 km for 1.5-hop backscattering.More ray tracing investigations (Jones 1966;Anderson et al. 1989;Norman and Cannon 1997) combined scattering and time varying nature of the ionosphere, providing additional field strength, absorption and Doppler shift parameters.A detailed discussion and comparison of these parameters is beyond the scope of this study.Here we attempt to locate a ground transmitting station after radio AOA determination by applying a numerical and stepped ray-tracing method developed by Tsai et al. (2009).This paper is organized as follows; section 2 introduces the Hualien VIPIR ionosonde system and its measurements in ionospheric sounding and background radio surveillance.The phase coherent measurements and analysis for least squares and aliases through passive-mode experiments are described in section 3. Section 4 describes radio AOA estimates of background interferences and presents sky-wave ray-tracing results from the radio source location system.Section 5 summarizes the Hualien VIPIR ionosonde observation results and their uses for radio transmitter position localization.

THE HUALIEN VIPIR IONOSONDE MEASUREMENTS
The Hualien VIPIR ionosonde system consists of two parallel and digital radio-frequency (RF) exciter channels generating RF pulses at milliwatt level and 0.5 -30 MHz frequency range for both transmission and calibration and eight parallel receiver channels connected to wide-band active dipole antennas with 10 dB preamplifiers.A wide-band power amplifier can raise the transmitted power level to 4 kilowatts.A half-delta wire antenna of 30 m length at 8 m height is used for transmission.The parallel receiver channels can simultaneously measure complex signal return amplitudes versus pulse envelop group delay and radio frequency to permit deriving the amplitude, phase and virtual range (to 5 microsecond resolution, i.e., 1.5 km).The system usually provides the so called 'P-mode ionograms' which record signal amplitude and phase at all range bins.An example ionogram is given in Fig. 1, which shows the relative amplitudes in color for signals including ionospheric echoes, background interferences, calibration signals at ~60 km, the strong signals at ~0 km from the transmitting antenna directly, and other background noises.The receiving dipole array has a scissorshaped layout with 34 m element spacing, as shown in Fig. 2. From the system exciter an overall calibration pulse is used to verify coherence and provide a mechanism for calibrating amplitude and phase differences external to the antenna preamplifier circuits, transmission cables and receivers.The upper panel in Fig. 3 shows an example of the observed, i.e., uncalibrated and 2π aliased, phases of the calibration pulse into eight receivers in the MF and HF bands.Because 2π aliasing is inherent to phase measurements, the phase profiles lie in the half open interval (-π, π] and are mostly sawtoothshaped functions of the calibration signal frequency.Some phase variations are induced by strong background radio interferences which also produce stronger and/or disturbed signal amplitudes shown in the bottom panel of Fig. 3.Because there are different lengths in the cables connecting dipole antennas and receivers, the phase profiles are not consistent from the eight receivers.After applying phase versus frequency linear fitting and calibration, the calibrated phases are shown in the center panel in Fig. 3 and are consistent with differences less than 2° for the eight receivers except for interfered frequencies. As in dynasondes, the Hualien VIPIR ionosonde utilizes an interferometric receiving dipole array to receive ionospheric echoes of transmitted pulse set and/or background radio interferences in the MF and HF bands.The system also transmits sets of radio pulses composed of four pulses with an interpulse interval of 10 msec and a radio carrier frequency configuration of f, f + δf, f + δf, and f, where f is the 1 st pulse carrier frequency within 0.5 and 30 MHz and δf is a radio frequency difference of 10 kHz.For each pulse set the measured parameters have one time-of-flight and thirty-two complex amplitudes from eight parallel receivers.The wave phase measurement methods are treated from the radio plane wave point of view.As in previous investigations reported by Pitteway and Wright (1992) and Tsai et al. (1997), an interferometric receiving array and sounding pulse set are used to introduce the six phase parameters ( o z , x z , y z , p z , t z , and f z ).Here o z is the reference phase obtained from the first pulse of each pulse set and at a fixed reference position, e.g., the center of Dipole "o" or "&" in Fig. 2. x z and y z are the phase changes due to antenna separation along the x-and y-axis directions separately.p z is the phase change from a north-south oriented dipole to an east-west oriented  where eight dipole antennas are labeled as "E", "e", "o", "w", "N", "n", "&", and "s" antenna separately, and the crossed-dipole center is a reference point for phase parameter analysis.dipole and is used to identify an ionospheric echo polarization in the ordinary mode (o-mode) or extraordinary mode (x-mode).t z is the phase change due to the time interval between pulses.f z is the phase change resulting from the small radio frequency increment δf (10 kHz in this study) within a pulse set and is used to estimate a group delay approximation using the stationary phase principle.Thus, ionospheric echo polarization, AOA, line-of-sight Doppler shift and phase path shall be derived using the six phase parameter determination.Ignoring the phase measurement residuals the phase path obtained for an individual echo is a linear function of these six phase parameters, which depend on the receiving array layout and the transmitted pulse-set pattern configuration.The configuration "quality" is quantified using the derived phase parameter ambiguities and the associated relative confidence limit factors.For details concerning the quality and receiving array and transmitted pulse-set pattern configuration evaluation the interested reader is referred to Tsai et al. (1997).
In addition to ionospheric sounding, the Hualien VIPIR ionosonde can also be operated in a passive mode for background radio surveillance in the MF and HF bands.The resulting signal amplitude image, as a function of range and frequency, is similar to Fig. 1, but there are no direct transmitted signals or calibration signals from the system and no echoes from the ionosphere or other reflectors.Figure 4 presents the observed 24-hour background radio signals received at Hualien, Taiwan, in relative amplitude as a function of the radio frequency and local time (LT) on the 121 th day in 2012.The signals within 1 to 1.6 MHz frequencies are ground-wave radios transmitted from local amplitude modulation stations.Most of the other strong interferences should be sky-wave radios from ground-based transmitters and after one-or multi-hop propagation through the ionosphere.As Fig. 4 shows the strong interferences are within 3 -12 MHz at night and 9 -18 MHz in the daytime due to the different night and day ionospheric characteristics.In the following study we focus on x z and y z determination and then sky-wave radio AOA from phase coherent measurements in the presence of measurement errors.

LEAST SQUARES AND ALIASES ESTIMATES OF ERROR AND PHASE PARAMETERS
To simplify phase parameter analyses for AOA, we used six antennas (labeled "e", "w", "o", "s", "&", and "n" dipoles in Fig. 2) and one pulse only in the sounding configuration and excluded t z and f z phase parameters.We analyzed the experimental phase data recorded in a passive and fixed-frequency mode.The signals were transmitted in a continuous wave mode at two frequencies of 7.7 and 8.6 MHz from different sites at Longquan (22.6702°N, 120.5971°E) and Chungli (24.9679°N, 121.1875°E), respectively, received by the Hualien VIPIR ionosonde.The experiment took place from 6 to 8 UT on the 121th day in 2012.It was a quiet geomagnetic period when Dst, Kp, and Ap indices were -6 nT, 1-, and 3, respectively (http://www.swpc.noaa.gov,http://swdcwww.kugi.kyoto-u.ac.jp).The radio AOA analysis problem is to determine x z and y z first from aliased phase measurements in the presence of measurement residuals.In a plane wave approximation the measured phases at Dipole "e", "w", "o", "s", "&", and "n" are sums of the phase parameters, aliasing items, residuals and could be represented as the following equations and in matrix form as: The phase parameter determination problem solution consists of finding the decomposition in Eq. ( 1), which contains a four-component state vector r 2 , r 3 , r 4 , r 5 , r 6 ] T ).The matrix A could be defined as the linear transformation matrix for the phase parameters in this study.The aliasing items are integers in phase cycle units.The residuals may include measurement errors and also background noise, systematic errors, antenna altitude and orientation bias, multi-reflection, etc.As usual we could assume that the residuals are small and approach zero.Meanwhile, if the number of state vector components is less than the number of measurements, i.e., the phase parameters are over determined and the matrix (A T A) is nonsingular and a least square estimate exists such that the sum of the residual squares is minimized.The state vector and residual vector can be estimated using  If it is the case that matrix M is an integer matrix then the least squares estimates from Eqs. ( 2) and (3) can be used to separate the phase-measurement aliasing items with a small residual assumption and then determine the phase parameters and residuals.However, this is not the case in the present analyses.We further note that any nonsingular realization of A T A results in the property that M is equivalent to M k for any integer k.This property of matrix M makes this matrix an idempotent matrix (Gantmakher 1959).We could define s as the expected sum of the residual squares and obtain where tr( ) is the trace operator, which is the sum of the diagonal elements of a square matrix.From the matrices theory, because the matrix M is idempotent and has a trace value of 2, we could obtain a four-order eigenvalue of 0 and a two-order eigenvalue of 1 for M (Gantmakher 1959).
Actually, the least squares residual vector could be linear to an error vector E, which has two elements assumed to have an independent Gaussian distribution with a bias and a standard deviation of σ.Thus, we obtain where B is a 6 × 2 matrix whose columns are a basis for the subspace of all vectors orthogonal to the columns of A. The sum of measured phases and aliasing items could be thus approximated using a least squares estimate linear function of the phase parameters and two independent error elements as in the following equation.
[ , ] Equation ( 6) is difficult to use as it stands, so we seek to modify it into a simple form by changing the matrix bases.The problem is also a basis reduction in the lattice theory.Harris (2003) applied the Hermite normal form of an integer matrix together with the Lenstra, Lenstra, and Lovász (LLL) basis reduction to determine the phase parameters and independent error elements of the dynasonde dipole arrays and pulse patterns including the Skiymet array (Jones et al. 1998), the L WHERPOL array (Jarvis and Dudeney 1986), and the X WHERPOL array (Wright and Pitteway 1982).WERPOL is derived from "where polarization", since it enables the system to determine the spatial location of an ionospheric echo and its polarization.In this study we do not follow the LLL basis reduction but apply an alternative method from Grötschel et al. (1988).If n is larger than or equal to m, every rational n × m matrix A of full column rank can be brought into a Hermite normal form using a series of elementary row operations.The series of elementary row operations produce a unimodular matrix U such that, UA is in a Hermite normal form.The unimodular matrix U could be not unique, but the Hermite normal form of a n × m matrix A is uniquely determined by A and has a form [ Al , 0] T , where Al is a nonsingular, upper triangular, nonnegative matrix in which each column has a unique maximum entry.Meanwhile, if the matrix A as a whole is integral, its Hermite normal form will also be an integer matrix.Thus, applying basis reduction on Eqs. ( 1) and ( 6), we obtain (7) Because the matrix U is now an integer matrix, two least squares error parameters e 1 and e 2 can be unaliased and estimated directly under the assumption that the two error parameters are small and have zero means.Figure 5a shows the scatter plot and percentage profiles of the two aliased errors, which were derived using Eq. ( 7) with an initial zero aliasing vector N for the measured experimental radio signal phases at 7.7 MHz.In the scatter plot the aliased errors are clustered at the points with integer coordinates, which mean the phase aliases in the cycle represent two linear combinations of the aliasing items.The results conform to expectations and the unaliased errors can be obtained directly by removing their integral phase cycles to approach zero as shown in Fig. 5b.The aliased and unaliased errors of the other 8.6 MHz experimental radio phase measurement are similar to that shown in Fig. 5 and are thus not shown in this paper.From Eq. ( 7) we also note that the matrix Al is a unit matrix in the present analyses.Because the reference phase o z will lie in the half open interval (-π, π], i.e., (-0.5, 0.5] in cycle, the phase parameter p z due to receiving dipole orientation change is known to be about ± 90° for o-and x-mode sky waves, respectively, and lies in (-π, π] for a combination wave.The phase parameters x z and y z are constrained by the antenna separations along the x-and y-axis directions separately.In general we could impose the constraint: ; ; in cycle for the least aliases.It follows that the aliasing items could be obtained by where we define $ ^h as a vector operator to return the least integer value greater than or equal to each vector element, known as the 'Ceil' function in many programming languages.The phase parameters can be solved using least squares and aliases by substituting Eq. ( 9) for N into Eq.( 7).
If any diagonal element of Al is not equal to one, the aliasing range of the corresponding phase parameter will be divided by the diagonal integer.The same as the aliased error parameter determination, aliased x z and y z were derived using Eq. ( 7) with an initial zero aliasing vector N. shows the aliased x z and y z values plotted against each other and their percentage profiles for both experimental radio measurements at 7.7 MHz (in red) and 8.6 MHz (in blue).The least squares and aliases x z and y z values are successfully determined from Eqs. ( 7) and ( 9) and shown in Fig. 7b after the aliasing items have been included.The resulting x z and y z values have statistical peaks of 0.122 (0.044) and 0.156 (-0.244) cycles, respectively, for the experimental 7.7 (8.6) MHz radio measurements.

POSITION DETERMINATION OF HF RADIO TRANSMITTING STATION
In this paper we consider radio direction finding based on interferometric phase measurements from a horizontal antenna array in the Hualien VIPIR ionosonde system.When a plane radio wave of wavelength λ arrives at identical and spaced receiving antennas, located similarly with respect to the ground and separated by known distances d x and d y along the x-and y-axis directions, respectively, the phase changes due to antenna separation can be defined using the following equations.where Δ is the elevation angle measured above the horizontal and y a is the azimuth angle from the y-axis direction.After finding the phase decomposition in Eq. ( 7), the least squares and aliases x z and y z values can be successfully estimated and the elevation and azimuth angles determined.For a statistical peak x z value of 0.122 (0.044) and a peak y z value of 0.156 (-0.244) in cycle, we obtain the corresponding elevation and azimuth angles to be 75.8°(75.3°) and 214.8° (346.9°),respectively, for the experimental 7.7 (8.6) MHz radio measurements.The results are listed in Table 1.Since the early days of radio communications, HF sky waves have been used extensively for long-distance radio service.HF radio ray tracing in the ionosphere mainly involves refraction and reflection due to the ionospheric structure on a scale that is largely compared with the radio wavelength.This study presents HF sky-wave transmitting station determination using a numerical and stepped ray-tracing method on a phenomenological ionospheric N e model, the TaiWan-Ionospheric Model (TWIM) (Tsai et al. 2009).The TWIM is constructed from a data base including the FormoSat3/Constellation Observing System for Meteorology, Ionosphere and Climate (FS3/COSMIC) ionospheric radio occultation (RO) data, Hualien VIPIR ionosonde data and global ionosonde foF2 and foE data.Generally, FS3/COSMIC can perform over 2500 ionospheric RO measurements per day with 80% of the RO measurements, roughly 2000 measurements on average, successfully retrieved into vertical N e profiles in 2006 and 1 st half of 2007.The number of profiles retrieved in 2013 declined to about 1000 vertical N e profiles due to FS3 spacecraft bus and GPS receiver/antenna payload degradation.The global ionosonde foF2 and foE values were automatically scaled from the ground-based ionosonde data from a network of 81 stations provided by the National Geophysical Data Center (NGDC) of NOAA and from http:// ngdc.noaa.gov/ionosonde/data/,which aggregates measurements and analyses from various international collaborators and augments these data with additional quality control and ionogram reduction.The three-dimensional TWIM consists of vertically-fitted α-Chapman-type layers with distinct F2, F1, E, and D layers, for which the layer parameters such as peak density, peak density height and scale height are represented by surface spherical harmonics.After noise separation optimization analyses we determined the cutoffs which are an order of 3 and a degree of 20 for each order of surface spherical harmonics.For details concerning TWIM performance and evaluation, interested readers are referred to Tsai et al. (2009) and Macalalad et al. (2013).It should be noted that in order to ensure good global and dense data coverage, every half-hour TWIM coefficients are estimated using a 30-day FS3/COSMIC data set and one-day ionosonde foF2 and foE data.We also instituted a one-hour universal time window and more weight (using Gaussian weighting) to days closer to the day of interest.For example, to generate the TWIM model at 7:30 UT on the 121 th day of 2012, we collected all of the RO observations at this universal time from the 92 th day to the 121 th day in 2012.The FS3/ COSMIC cannot provide dense enough RO observations to observe disturbed N e profiles when magnetic storms occur.Therefore, the TWIM could be suggested for use when the ionosphere is quiet for geomagnetic fields as the choice for the 121 th day of 2012 in this study.Note also that the TWIM is a three-dimensional and continuous model and the continuity of the first and second derivatives for the practical schemes can be maintained.This is important for providing reliable radio propagation predictions described as follows.
A three-dimensional TWIM ray-tracing algorithm developed by Tsai et al. (2010) was applied to determine the experimental HF radio transmitter positions.The TWIM ray tracing was simulated upon independent vectors of initial ray direction from the AOA determination, ionospheric N e model and the magnetic field.Two different geomagnetic field models were used and superimposed onto the TWIM N e maps to show the characteristic features of the geomagnetic configuration in this study.The first is an Earth-centered dipole field assumed to be inversely proportional to the cube of ray altitude from the Earth center.A better one is the Version 11 of the International Geomagnetic Reference Field (IGRF) (Finlay et al. 2010; http://www.ngdc.noaa.gov/geomagmodels/IGRFWMM.jsp), which has a carefully managed spherical-harmonic representation.The  Using a large time-of-flight interval is clearly more efficient, however, one needs to limit the truncation error.We considered a ray-tracing path step about one order less than the TWIM vertical scale and a time-of-flight interval value of 10 -5 sec has been used in this study.From the specified ray-tracing universal time the corresponding surface spherical harmonic coefficients of the TWIM are read to reconstruct a three-dimensional ionosphere.Once the arriving elevation and azimuth angles of a ray have been derived ray tracing begins from the Hualien VIPIR ionosonde site as the initial aiming point for determining the wave vector of the beginning path step.For convenience, assume that the ray is traveling straight along each short step distance within a constant time-of-flight and then enters a new path step.Onto the new step a TWIM N e gradient is calculated as the interface normal vector and Snell's law is applied to derive a new wave vector.The ray tracing process is then repeated step by step as the ray travels up through the ionosphere and is reflected back down to the ground.
Figure 7 shows the simulated ray-tracing paths of the experimental 7.7 and 8.6 MHz radios based on the IGRF superimposed on the TWIM at 7:40 UT (15:40 LT) on the 121 th day of 2012.The simulation resulted in two traces including the o-and x-mode rays in white (gray) and yellow (blue) colors, respectively, for the experimental 7.7 (8.6) MHz radio with the same AOA, which were determined using the statistical least squares and aliases peak x z and y z values, as shown in Fig. 6 and Table 1.The simulated ray traces are shown and projected onto the Earth's surface (Fig. 7a), the meridional and vertical plane (Fig. 7b) and the longitudinal and vertical plane (Fig. 7c), where the background imaging colors represent the modeled plasma frequencies from the TWIM.The simulated transmitter coordinates are also listed in Table 1.The positioning results have distance errors of 0.24 (4.8) and 5.0 (4.7) km based on o-and x-mode ray tracing, respectively, for the experimental 7.7 (8.6) MHz radio transmitted from Longquan (Chungli), Taiwan.In relation to a great circle distance of 124 (176) km the percentage errors are 0.2% (2.9%) and 4.0% (2.8%) based on o-and xmode rays, respectively, for the experimental 7.7 (8.6) MHz radio.We further note that though the 8.6 MHz o-and xmode error (0.24 and 5.0 km) difference is larger than the 7.7 MHz error (4.8 and 4.7 km) difference, the location distance is 5.1 km, which is less than 9.2 km for 8.6 MHz.
Other ray-tracing results based on the Earth-centered geomagnetic dipole field were conducted in a similar manner as shown in Fig. 7 and are not shown in this paper.However, the positioning results are also listed in Table 1.Compared with the ray-tracing results based on the IGRF the positioning errors using the dipole field have greater distance errors of 9.2 (8.6) and 10.3 (8.4) km based on o-and x-mode rays, respectively, for the experimental 7.7 (8.6) MHz radio.The results confirm the expectations from more accurate geomagnetic fields from the IGRF.

SUMMARY OF RESULTS
The Hualien VIPIR ionosonde has been operated in active mode for ionospheric sounding and in passive mode for detecting background radio interferences.The passivemode observations are useful for background HF radio surveillance and solving transmitter localization problems.The problem of interest is the position determination of a licensed or unlicensed source of radio transmission.We have proposed a HF sky-wave radio positioning technique using interferometric phase measurements from the Hualien VIPIR ionosonde and further ray tracing on the IGRF or an Earth-centered geomagnetic dipole field superimposed on the TWIM.The phase parameters due to antenna separation along the x-and y-axis directions and different dipole orientation were determined and used to obtain radio AOA.An estimate of least square error phase parameters exists when the phase parameters are over determined from interferometric phase measurements and the matrix (A T A) is nonsingular, where A was defined as the linear transformation matrix for the phase parameters.We also applied the Hermite normal form and lattice methods to solve the phase parameter and error determination lease alias problem.The derived AOA value can be used as an initial aiming direction for a backward numerical ray-tracing simulation on the TWIM.Two different geomagnetic field models of the IGRF and an Earth-centered dipole field were used and superimposed onto the TWIM for ray tracing.Based on the IGRF the positioning results have distance errors of 0.24 (4.8) and 5.0 (4.7) km for o-and x-mode rays, respectively, of the experimental 7.7 (8.6) MHz radio transmitted from Longquan (Chungli).Related to a great circle distance of 124 (176) km the experimental 7.7 MHz (8.6 MHz) radio percentage errors are 0.2% (2.9%) and 4.0% (2.8%) based on o-and x-mode ray tracing, respectively.The positioning results using the Earth-centered geomagnetic dipole field have greater distance errors as expected.We believe that the dominant error in HF radio transmitter position is due to the idealistic TWIM ionospheric model.In future investigations a follow-on mission evaluation is planned to extend the current FS3/COSMIC mission.It will place twelve satellites and enhance the GPS receiver capability to accommodate signals from the Tri-G system (GPS, GLONASS, and Galileo).The RO profiles in the joint follow-on mission should be increased by a factor of ~5 and be no less than 8000 on average per day.It is expected that more dense GPS RO measurements will result in a more realistic ionospheric N e model and subsequently improved HF sky-wave radio positioning.

Fig. 1 .
Fig. 1.An example of a typical ionogram; relative amplitudes are shown in color.

Fig. 3 .
Fig. 3. Examples of measured calibration phase and amplitude profiles onto eight receivers in different colors.The upper panel shows observed, i.e., uncalibrated and aliased, phases of the calibration pulses into eight receivers in MF and HF bands, and the corresponding calibrated phases and measured signal relative amplitudes are shown in the middle and bottom panels, respectively.

Fig. 4 .
Fig. 4. Surveillance of MF and HF background radio interferences at Hualien, Taiwan, on the 121 th day in 2012.

Fig. 5 .
Fig. 5. Scatter plots of two error parameters against each other and their percentage profiles for the measured experimental radios at 7.7 MHz.Aliased and unaliased errors are shown in the upper panel (a) and the lower panel (b), respectively.The relative signal amplitudes are shown in color.

Fig. 6 .
Fig. 6.Scatter plots of phase parameter x z against y z and their percentage profiles shown in blue and red colors for the measured experimental radios at 7.7 and 8.6 MHz, respectively.Aliased x z and y z values are shown in the upper panel (a), and unaliased x z and y z values are shown in the lower panel (b).

Fig. 7 .
Fig. 7. Simulated ray-tracing paths of the experimental 7.7 and 8.6 MHz radios, which are received by the Hualien VIPIR ionosonde and have the AOA values listed in Table 1.The upper (a), middle (b), and bottom (c) panels show the projected ray-tracing paths onto the Earth's surface in the meridional and vertical planes and the longitudinal and vertical planes, respectively.The simulated o-and x-mode traces are shown in white (gray) and yellow (blue) colors, respectively, for the experimental 7.7 (8.6) MHz radio, and the background imaging colors represent the modeled plasma frequencies from the TWIM.
procedure starts by specifying a radio frequency, type of o or x-mode wave characteristic, receiver position, arriving elevation and azimuth angles, universal ray tracing time and time-of-flight interval for backward ray tracing.

Table 1 .
The positioning results for two experimental radio transmitters at Longquan and Chungli.