Ionospheric electron density profile estimation using commercial AM broadcast signals

Journal of Atmospheric and Solar-Terrestrial Physics 130-131 (2015) 1–13

Contents lists available at ScienceDirect

Journal of Atmospheric and Solar-Terrestrial Physics journal homepage: www.elsevier.com/locate/jastp

Research Paper

Ionospheric electron density profile estimation using commercial AM broadcast signals De Yu a,n, Hong Ma b, Li Cheng b,c, Yang Li b, Yufeng Zhang b, Wenjun Chen a a

School of Physics, Huazhong University of Science and Technology, Wuhan 430074, China Department of Electronics and Information Engineering, Huazhong University of Science and Technology, Wuhan 430074, China c School of Electrical and Electronic Engineering, Wuhan Institute of Technology, Wuhan 430074, China b

art ic l e i nf o

a b s t r a c t

Article history: Received 24 October 2014 Received in revised form 8 May 2015 Accepted 9 May 2015 Available online 13 May 2015

A new method for estimating the bottom electron density profile by using commercial AM broadcast signals as non-cooperative signals is presented in this paper. Without requiring any dedicated transmitters, the required input data are the measured elevation angles of signals transmitted from the known locations of broadcast stations. The input data are inverted for the QPS model parameters depicting the electron density profile of the signal’s reflection area by using a probabilistic inversion technique. This method has been validated on synthesized data and used with the real data provided by an HF directionfinding system situated near the city of Wuhan. The estimated parameters obtained by the proposed method have been compared with vertical ionosonde data and have been used to locate the Shijiazhuang broadcast station. The simulation and experimental results indicate that the proposed ionospheric sounding method is feasible for obtaining useful electron density profiles. & 2015 Elsevier Ltd. All rights reserved.

Keywords: Ionosphere HF AM broadcast QPS model a posteriori information

1. Introduction The ionosphere is the part of the atmosphere that is ionized by solar radiation, and this ionization is responsible for the refraction of HF radio waves. The ionosphere is used as a propagation channel by Skywave OTH radars, HF AM broadcasts and other HF communications. The real-time ionospheric electron profile must be understood in order to characterize HF radio channels. Groundbased HF radio ionospheric sounding techniques, including vertical sounding, oblique sounding, and oblique backscatter sounding, are still important methods for studying ionospheric electron density profiles. There are hundreds of ionosonde stations around the world that record vertical incidence ionograms day and night, and a significant amount of work has been invested in the automatic calculation of electron density profiles using these ionograms (Galkin et al., 1996, Jiang et al., 2014, Reinisch and Huang, 1982, Reinisch and Huang, 1983, Titheridge, 1985, Zheng et al., 2013). The oblique sounding method obtains echo traces in the form of a group path as a function of frequency. Inversion algorithms are available for deriving the vertical electron density profile at the midpoint between the transmitter and the receiver of an oblique sounder path from the resulting ionograms (Chuang and Yeh, 1977, n

Corresponding author. E-mail address: [email protected] (D. Yu).

http://dx.doi.org/10.1016/j.jastp.2015.05.008 1364-6826/& 2015 Elsevier Ltd. All rights reserved.

Heaton et al., 2001, Reilly, 1985). The oblique backscatter sounding method records the backscatter signals received from distant locations on the earth’s surface via reflection by the ionosphere. Both the transmitted and received signals traverse the ionosphere, and these signals contain useful information regarding the state of the intervening ionosphere over the range of the returned signal, which can be up to a few thousand kilometers from the transmitter/receiver location. There are two types of backscatter soundings. The first type consists of transmitting a signal at a fixed azimuth by scanning in frequency and obtaining backscatter ionograms in the form of a 3-D figure with the received signal amplitude, frequency, and group path. The signal strength at the minimum group path will be strong due to the time delay focusing and the skip distance focusing. The minimum group path as a function of operating frequency is called the leading edge of the backscatter ionogram. Many researchers have been devoted to the leading edge inversion technique as a method for obtaining ionospheric electron density profiles (Dyson, 1991, Fridman and Fridman, 1994, Fridman, 1998, Norman, 2003, Rao, 1974). The second type of backscatter sounding consists of transmitting a signal at a fixed azimuth and operating frequency by scanning in elevation and obtaining backscatter ionograms in the form of a 3-D figure with the received signal amplitude, elevation angle, and group path. By inverting this type of ionograms, researchers can also obtain information about the ionospheric electron density distribution (Benito et al., 2008, Caratori and Goutelard, 1997, Jacquet et al., 2001, Norman and

2

D. Yu et al. / Journal of Atmospheric and Solar-Terrestrial Physics 130-131 (2015) 1–13

Dyson, 2006, Ruelle and Landeau, 1994). The active ionospheric sounding techniques mentioned above all need high-power transmitters, and the sounding system is complex and costly. Therefore, some researchers have started to detect the ionosphere in a passive way, by collecting and analyzing the signals emitted by non-cooperative transmitters to obtain information about the ionosphere. Earlier passive ionospheric sounding techniques used commercial AM broadcast signals to evaluate ionospheric channel characteristics, perform short-term ionospheric parameter forecasts (Jowett et al., 1989, Masrani and Riley, 1991, Piggin et al., 1996), and detect traveling ionospheric disturbances (Beley et al., 1995). Next, the Manastash Ridge Radar (MRR) was developed, which collected commercial FM broadcast signals to detect the field-aligned irregularity (FAI) of the E layer (Lind et al., 1999, Meyer et al., 2004). However, none of these passive detection techniques mentioned above can provide the ionospheric electron density distribution. Finally, the appearance of the computerized ionospheric tomography (CIT) imaging technique provided a passive way to obtain the electron density distribution. With the aid of the tomographic technique and the total electron content (TEC) along the ray paths between the ground receivers and satellites, CIT can reconstruct the electron density distribution on a global scale (Debao et al., 2008, Fremouw et al., 1992, Raymund et al., 1990). However, due to the lack of horizontal ray paths the vertical precision obtained by CIT is not good. Researchers have proposed many methods to solve this problem, such as bringing in radio occultation data by carrying receivers on low Earth orbit (LEO) ́ satellites (Tsai et al., 2001); employing ionosonde data (GarcıaFernández et al., 2003, Ma et al., 2005); absorbing incoherent scatter radar observations (Chartier et al., 2012); and using the vertical ionograms, backscatter ionograms and CHAMP satellite data together (Fridman and Nickisch, 2001, Fridman et al., 2009, Zhao et al., 2010). These methods can improve the vertical precision to a certain degree; however, these methods cannot obtain high-precision imaging because of the sparseness of ionosondes. Employing more measurement sensors in CIT is an important goal for improving the vertical precision (Bust and Mitchell, 2008), but the cost of getting more data from a dense ionosonde array is high. Therefore, it is necessary to obtain as many ionospheric measurements as possible from the pre-existing transmitters and receivers. By using commercial AM broadcast signals as non-cooperative signals, a new method for obtaining the bottom side electron density profile of the ionosphere is presented in this paper. In this method, the quasi-parabolic segment (QPS) model derived by Dyson and Bennett (1988) is chosen to depict the electron density profile, and the required input data are the elevation angles as a function of the signal frequency, which are obtained from the known locations of the broadcast stations. The input data are inverted for the QPS model parameters by using a probabilistic inversion method. This new method has been validated using synthesized data and has been applied to real data to produce a realistic ionospheric electron density profile. This method does not need any dedicated transmitters. Instead, the method just uses an HF direction finding (HFDF) system to receive and analyze the commercial AM HF broadcast signals. The profiles provided by our method can be used alone (e.g., used at a single-site location (SSL)), and they can also be absorbed into CIT to get more accurate density distributions.

2. Extract the elevation angle data by using an HFDF system Commercial AM HF broadcast stations usually transmit signals between 5 MHz and 20 MHz. These signals can be reflected by the

ionosphere and received thousands of kilometers away from the transmitter station. Currently, there are numerous commercial AM broadcast stations all over the world transmitting signals in the HF band day and night. By using an HFDF system, we can receive these signals that are transmitted by the stations located near the HFDF system at various azimuth and ground ranges. Each of these signals that is reflected from the ionosphere contains useful information regarding the state of the ionospheric electron density profile in the area where the signals have been reflected. This area is usually the midpoint of the radio signal path from the transmitter to the receiver in a one-hop propagation model. The HFDF system using an array of monopoles to receive the signals and all of the received HF band signals are filtered to identify the commercial AM signals. The transmitter location of each AM signal can then be determined and the conventional multiple signal classification (MUSIC) method is applied to do the direction of arrival (DOA) estimation: The signal subspace and noise subspace are obtained by eigen decomposition of the covariance matrix of the array signal. By using the orthogonal characteristics of eigenvectors in the signal and noise subspaces, the two-dimensional (azimuth and elevation) MUSIC spectrum is obtained, a two-dimensional searching is employed to find the direction that maximizes the MUSIC spectrum, which is the estimated arriving azimuth and elevation angle. As a result, we can get an “oblique sounding ionogram” in the form of elevation angle as a function of frequency from every broadcast station, but the ionogram only has discrete data points at the frequencies on which the broadcast station transmits its signals. A flowchart depicting the signal collecting and processing method in the HFDF system is shown in Fig. 1. As shown in the flowchart, the HFDF system samples and records a section of HF full-band signals and then performs a power spectrum analysis and selects the signals with strong power (i.e., higher than
80 dBm). The ITU-HF broadcasting schedule is referenced to determine which signals could be reasonably received at the HFDF system location at that time and from which broadcast stations they emanate. We then channelize the wide-band signals to get each of the broadcast signals we identified above and apply a two-dimensional DOA estimation to each signal to find the arriving azimuth and elevation angles. We can then filter these signals again using the information from the broadcast content (e.g., the specific broadcasting language) and the deviation (less than 3°) between the measured and geographic azimuth angles to ensure that all of the selected signals have been matched to the correct broadcast stations. Finally, we obtain a set of elevation angles that corresponds to a set of ground ranges and frequencies. This set of information is the input data for the inversion method.

3. Quasi-parabolic segment (QPS) model and the forward problem The QP layer defined by Croft and Hoogasian (1968) is given by

⎧ ⎡ 2 2⎤ ⎛ r ⎞ ⎪ ⎢ ⎛ r−rm ⎞ ⎛⎜ rb ⎞⎟ ⎥ ⎟⎟ ⎪ Nm 1−⎜⎜ , rb < r < rm⎜ b ⎟ ⎥ ⎢ ⎝ ⎠ ⎨ N (r ) = ⎝ rb−rm ⎠ ⎝ ym ⎠ r ⎦ ⎪ ⎣ ⎪ 0 , otherwise ⎩

(1)

where N (r ) represents the electron density at a radial distance r from the earth’s center, Nm represents the maximum electron density at the radial height rm , rb is the radial base height of the ionospheric layer, and ym = rm − rb is the layer semi-layer thickness. Based on the QP layer, Dyson and Bennett (1988) proposed the QPS model, which consists of multiple smoothly joined QP

D. Yu et al. / Journal of Atmospheric and Solar-Terrestrial Physics 130-131 (2015) 1–13

3

Fig. 1. Flowchart for estimating the ionospheric electron density profile by using commercial AM broadcast signals.

segments. The equations describing the QPS ionospheric model, which include the E, F2 and joining layers, may be written as shown below.

⎛ r E ⎞2 NE(r ) = aE−bE ⎜1− m ⎟ , ⎝ r ⎠

(2)

⎛ r E ⎞2 Nj(r ) = aE + bj ⎜1− m ⎟ , ⎝ r ⎠

(3)

⎛ r F ⎞2 NF2(r ) = a F2−b F2⎜1− m 2 ⎟ . ⎝ r ⎠

(4)

rj =

rmF2b F2 a F2−aE +

(

rmF2 −1 rmE

)

(

r F b F2 rm E2 −1 m

)

, (7)

⎛ r F ⎞ −rmF2b F2⎜1− mr 2 ⎟ ⎝ j ⎠ bj = . ⎛ rmE ⎞ rmE⎜1− r ⎟ ⎝ j ⎠

(8)

The electron density profile is then determined via the QPS model using six parameters: fo E, rbE, rmE, fo F2, rbF2, and rmF2 Ignoring the geomagnetic field, the ground range of a radio signal propagating through the ionosphere is given by

Eqs. (2), (3) and (4) represent the electron density distribution of the E layer, joining layer and F2 layer, respectively. For the E and F2 layers

D = 2re(γ −β ) + 2re2 cos β

a = Nm = 12407fo2 ,

where μ2 = 1−fp2 /f 2 is the square of the refractive index,

(

(5)

2

)

b = Nm rb/ym .

(6)

where fo is the critical frequency of the layer. By matching the electron density and its gradient of the joining layer with that of the E and F2 layers at the joining points, we obtain the joining point height r j of the joining and F2 layers and the joining layer parameter bj .

∫r

rt b

dr 2 2

r r μ −re2cos2β

, (9)

fp = N/12407 is the ionospheric plasma frequency, f is the radio frequency, β is the initial elevation angle of the ray, rt is the radial height at which the radio signal is reflected, re is the radius of the earth, and γ is the elevation angle of the ray at the base of the ionosphere. For a ray that propagates through n segments in a spherically stratified ionosphere described by the QPS model, there is an explicit equation in the model for finding the ground range (Dyson and Bennett, 1988)

D. Yu et al. / Journal of Atmospheric and Solar-Terrestrial Physics 130-131 (2015) 1–13

2000

45

1800

40

1600

35 Elevation Angle(Degree)

Ground Range(Km)

4

1400 1200 1000 800

20

10 5

10

15

20

25

30

35

40

45

50

Elevation Angle(Degree) Fig. 2. Ground range as a function of elevation angle under the following conditions: f ¼11.76 MHz, foE ¼3.27 MHz, rbE ¼90 km, rmE ¼110 km, foF2 ¼ 10.46 MHz, rbF2 ¼ 179.02 km, and rmF2 ¼ 268.53 km. The ionospheric parameter used here is the IRI2012 model value at the middle point between Beijing and the receiver location at 0445 UT, 31 October 2012. n ⎧ ⎫ γ −β + re cos β ∑ (Ui−L i )⎬ D = 2re⎨ , ⎪ ⎪ ⎩ ⎭ i=1 ⎪

⎪

(10)

where Ui and L i represent the values of the integral in Eq. (9) at the upper and lower bounds, respectively. Hereinafter, the following parameters will be noted as described: the ionospheric QPS model parameters will be noted as T m = ⎡⎣f E, r E, r E, f F , r F , r F ⎤⎦ , the frequencies of the received o

25

15

600 400

30

b

m

o 2

b 2

m

T signals as f = ⎡⎣f1 , f2 ⋯fn ⎤⎦ , the corresponding ground ranges from T

broadcast stations to the receiver station as D = [D1, D2⋯Dn] , and T

the arriving elevation angles as β = ⎡⎣β1, β2⋯βn⎤⎦ When the ionospheric parameter m , signal frequencies f and radio launching elevation angle β are known, the ground range D of the ray propagation can be determined by Eq. (9) to obtain a theoretical relationship between m , β and D as

D = D(m, β).

(11)

Fig. 2 shows an example of the propagating ground range D as a function of the launch elevation angle β of a signal transmitting from Beijing when the ionospheric parameters are known. Similarly, we can obtain another theoretical relationship between m , β and D (Eq. (12)) when the ionospheric parameter m , the signal frequencies f and the ground range D from the broadcast station to the receiver location are known. The β are determined by means of scanning β at a certain interval to calculate the ground range at each fi and by finding the βi for which the ground range equals Di as the arriving elevation angle (which equals the launching elevation angle under the condition of a spherically stratified ionosphere).

β = β(m, D).

(12)

Fig. 3 shows an example of the different launch elevation angles β of signals transmitting from Beijing on different frequencies when the ionospheric parameters are known and the ground range D is fixed. 4. Inversion method Eqs. (11) and (12) solve the forward problem: determining the ground range D with the ionospheric parameter m when signal

5

4

6

8

10 12 Frequency(MHz)

14

16

18

Fig. 3. Elevation angle as a function of signal frequency under the following conditions: D ¼1059.6 km, f¼ 11.76 MHz, foE¼ 3.27 MHz, rbE ¼90 km, rmE¼ 110 km, foF2 ¼10.46 MHz, rbF2 ¼179.02 km, and rmF2 ¼ 268.53 km. The ionospheric parameter used here is the IRI2012 model value at the middle point between Beijing and the receiver location at 0445 UT, 31 October 2012.

frequencies f and launching elevation angles β are known or determining the arriving elevation angles β with the ionospheric parameter m when signal frequencies f and the ground ranges D are known. Here we present a probabilistic approach to solving the inverse problem: determining the ionospheric parameter m with the arriving elevation angles β when signal frequencies f and the ground ranges D are known. A consistent formulation of inverse problems based on probability theory was proposed by Tarantola, a French mathematician, in the 1980s. Central to this method is the concept of the “state of information” over a parameter set. The measurements of the observable parameters (data), the a priori information concerning model parameters, and the information on the theoretical correlations between the observable parameters and model parameters can all be described by probability densities. Therefore, the general inverse problem can be set as a problem of “combining” all of this information to obtain a posteriori information concerning the model parameters (Tarantola, 2005). The set of model parameters that maximizes the a posteriori probability density is then found by an optimization algorithm and is expected to be the result of inversion. 4.1. The a priori information The a priori information on model parameters is the information obtained independently without the observed values. It describes all of the knowledge we know before the measurement. The a priori probability density on ionospheric model parameters can be obtained by analyzing the long-term ionosonde data. Here we did not bring in that information, and the a priori probability density on model parameters should be constant

ρm (m) = const .

(13)

The a priori information on data is the information obtained by measurement. All physical measurements are subject to uncertainties. Therefore, the result of a measurement act is not simply a set of “observed values” but a “state of information” acquired on some observable parameters. We assume that our measurements are subject to independent Gaussian uncertainties; therefore, the a priori probability densities on the ground range data D between the transmitter and receiver and the a priori

D. Yu et al. / Journal of Atmospheric and Solar-Terrestrial Physics 130-131 (2015) 1–13

probability densities on the arriving elevation angle data β can be written as

35

⎡ 1 ⎤ T ρ D (D) = kexp⎢− (D−Dobs) CD−1(D−Dobs)⎥, ⎣ 2 ⎦

30

(

)

(

)

(15)

respectively. In Eqs. (14) and (15), Dobs (βobs ) are observed values of the ground ranges (arriving elevation angles), k is the normalization constant of the probability densities, CD and Cβ are the covariance matrices describing the Gaussian type uncertainties

⎛σ 2 0 ⎜ Dobs1 ⎜ 0 2 σ Dobs2 CD = ⎜ ⎜ 0 0 ⎜ ⎜ 0 0 ⎝ ⎛σ 2 0 ⎜ βobs1 ⎜ 0 2 σ βobs2 Cβ = ⎜ ⎜ 0 0 ⎜ ⎜ 0 0 ⎝

0 0 ⋱ 0

0 ⎞ ⎟ 0 ⎟⎟ , 0 ⎟ ⎟ σ D2obsn ⎟⎠

20

15

10

0

(16)

0 ⎞ ⎟ 0 0 ⎟⎟ ⋱ 0 ⎟ ⎟ 0 σ β2obsn ⎟⎠

9

10

11

12

(17)

15

16

17

18

(18)

Using Eqs. (13), (14) and (15), we obtain

250

200

150

ρ(m, D, β) = ⎧ 1⎡ T kexp⎨− ⎢(D−Dobs) CD−1(D−Dobs) ⎩ 2⎣ T ⎤⎫ + β−βobs Cβ−1 β−βobs ⎥⎬ ⎦⎭

(

"True" profile Inverted profile

300

Height(Km)

ρ(m, D, β) = ρm (m)ρ D (D)ρβ (β)

)

100

(19) 50

0

2

4

4.2. The theoretical information From Section 3, we know there are some physical theoretical correlations among the ionospheric model parameter m , the signal arriving elevation angle data β and the ground range data D . The homogeneous probability densities of these three parameter sets are constant ( μm (m) = const ., μβ (β) = const. and μD (D) = const .) and will be accounted for in the normalization constant k later in this paper. The correlation of Eq. (11) can be described by a probability density

θ (m, D, β) = kδ(D−Dcal(m, β)).

)

θ (m, D, β) = kδ β−βcal(m, D) .

6 8 Plasma Frequency(MHz)

10

12

Fig. 5. The ‘true’ electron density profile and the inverted profile calculated with the set of synthesized measured elevation angle data. Table 1 Comparisons of the true and inverted QPS model parameters.

True Inverted

foE (MHz)

rbE (Km)

rmE (Km)

foF2 (MHz)

rbF2 (Km)

rmF2 (Km)

3.27 3.40

90.00 85.27

110.00 113.58

10.46 10.40

179.02 172.74

268.54 265.24

(20)

Similarly, the correlation of Eq. (12) can be described by another probability density

(

13 14 Frequency(MHz)

Fig. 4. A set of synthesized true elevation angles (red circles), measured elevation angles (blue stars) and the elevation angles calculated using the inverted ionospheric electron density profile (black squares) (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.).

0

)

25

5

where σDobsi and σ β are the standard deviations of the ith ground obsi range and arriving elevation angle, respectively. The a prior information we have on both the model parameters and observable parameters (data) can then be described by a joint probability density

(

Ture Elevations Measured Elevations Inverted Elevations

(14) Elevation Angle(Degree)

⎡ 1 ⎤ T ρβ (β) = kexp⎢− β−βobs Cβ−1 β−βobs ⎥, ⎣ 2 ⎦

5

(21)

where Dcal(m, β) and βcal(m, D) are the ground range and arriving elevation angle calculated by Eqs. (11) and (12), respectively. 4.3. The a posteriori information on model parameters The a posteriori probability density on the ionospheric model parameters, the arriving elevation angles and the ground ranges

can be produced by the conjunction of the a priori information and the theoretical information (Tarantola, 2005)

σ (m, D, β) = kρ(m, D, β)θ (m, D, β).

(22)

The a posteriori information on the ionospheric model parameter m is then given by the marginal probability density

σm(m) =

∬D, β σ (m, D, β)dDdβ.

By substituting Eqs. (19) and (22) into (23), we have

(23)

6

D. Yu et al. / Journal of Atmospheric and Solar-Terrestrial Physics 130-131 (2015) 1–13

35

35

True Elevations

Measured Elevations

30 Elevation Angle(Degree)

Elevation Angle(Degree)

30 25 20 15 10 5

25 20 15 10 5

0

10

12

14

16

0

18

10

12

Frequency(MHz)

35

16

18

"True" profile Inverted profile

250 25 Height(Km)

Elevation Angle(Degree)

300

Elevations calculated with estimated parameters

30

14

Frequency(MHz)

20 15

200

150

10 100 5 0

10

12

14

16

18

50

0

2

4

Frequency(MHz)

6

8

10

12

Plasma Frequency(MHz)

Fig. 6. Simulation test results: (a) the set of true elevation angles; (b) 100 sets of simulated measured elevation angles; (c) 100 sets of inverted elevation angles; and (d) 100 sets of inverted ionospheric profiles and the true profile.

Table 2 Statistical results of 100 times independently simulation tests.

σm(m) =

T ⎤⎫ + β−βobs Cβ−1 β−βobs ⎥⎬dβ. ⎦⎭

foE (MHz) rbE (km) rmE (km) foF2 (MHz) rbF2 (km) rmF2 (km) True Inverted mean Inverted RMSE

3.27 3.41

90.00 89.94

110.00 110.43

10.46 10.48

179.02 181.82

268.54 266.11

0.22

4.08

3.76

0.41

10.05

16.86

(

(

)

(26)

(

∬

(24)

In addition, by substituting Eqs. (20) and (22) into (23) we have ⎧

⎡

⎤⎫

∬D, β kexp⎨⎩− 21 ⎣⎢(D−Dobs )T CD−1(D−Dobs ) + (β−βobs ) Cβ−1(β−βobs )⎦⎥⎬⎭

δ(β−βcal (m, D))dDdβ .

By integrating D and β , we obtain

T

⎧

∫D kexp⎨⎩− 21 ⎡⎣⎢(D−Dobs)T CD−1(D−Dobs) T ⎤⎫ + βcal(m, D)−βobs Cβ−1 βcal(m, D)−βobs ⎥⎬dD. ⎦⎭

⎧ 1⎡ T ⎤⎫ T k exp⎨− ⎢(D−Dobs ) CD−1(D−Dobs ) + (β−βobs ) C β−1(β−βobs )⎥⎬ ⎦⎭ ⎩ 2⎣ D, β

δ(D−Dcal (m,β))dDdβ .

σm(m) =

)

and

σm(m) =

σm(m) =

⎧

∫β kexp⎨⎩− 21 ⎡⎣⎢(Dcal(m,β)−Dobs)T CD−1(Dcal(m,β)−Dobs)

)

(

)

(27)

In this real situation, the ground ranges D are calculated using the longitude and latitude data of the transmitter and receiver locations. The longitude and latitude data have very high accuracy. Therefore, the ground ranges D have very small measurement errors, which can be neglected with respect to the value of D. Eq. (14) can be written as

(25) ⎡ 1 ⎤ T ρ D (D) = kexp⎢− (D−Dobs ) CD−1(D−Dobs )⎥ ≈ δ(D−Dobs ). ⎣ 2 ⎦ Substituting Eq. (28) into (27), we have

(28)

D. Yu et al. / Journal of Atmospheric and Solar-Terrestrial Physics 130-131 (2015) 1–13

45 40

1100

35

Ground Range(Km)

Elevation Angle(Degree)

1200

True Elevations Measured Elevations

30 25 20

True Range Calculated Range with True Profile and Measured Elevation

1000 900 800 700

15 10 11.7

11.72

11.74

11.76

11.78

600 11.7

11.8

11.72

Frequency(MHz)

1200

11.74

11.76

11.78

11.8

Frequency(MHz)

1200

True Range Calculated Range with Inverted Profile and Measured Elevation

1100 Ground Range(Km)

1100 Ground Range(Km)

7

1000 900 800 700

True Range Calculated Range with Inverted Profile and True Elevation

1000 900 800 700

600 11.7

11.72

11.74

11.76

11.78

11.8

Frequency(MHz)

600 11.7

11.72

11.74

11.76

11.78

11.8

Frequency(MHz)

Fig. 7. The location results for the Shijiazhuang station on 3 different frequencies: (a) the true and measured elevations of the signals; (b) the ground ranges calculated with the true ionospheric profile and the measured elevations; (c) the ground ranges calculated with the inverted ionospheric profile and the measured elevations; and (d) the ground ranges calculated with the inverted profile and the true elevations. Table 3 Statistical results of ground ranges calculated under different conditions. Signal frequency (MHz)

True range (km)

Calculated range in Fig. 7(b) (km) Mean RMSE

Calculated range in Fig. 7(c) (km) Mean RMSE

Calculated range in Fig. 7(d) (km) Mean RMSE

11.72 11.75 11.76

845.10 845.10 845.10

859.92 848.36 847.92

859.50 847.54 846.99

844.17 844.13 844.12

σm(m) =

50.12 50.35 46.24

⎡

∫D kexp⎢⎣− 21 (βcal(m, D)−βobs)

T

52.89 50.95 47.23

11.21 11.18 11.17

⎤ Cβ−1 βcal(m, D)−βobs ⎥δ(D−Dobs )d ⎦

(

)

⎤ ⎡ 1 T D= kexp⎢− βcal(m,Dobs )−βobs Cβ−1 βcal(m,Dobs )−βobs ⎥ = kρβ βcal(m,Dobs ) . ⎦ ⎣ 2

(

)

(

)

(

)

(29)

As Eqs. (26) and (27) show, the a posteriori probability density on model parameters should be obtained by calculating an integral of D or β when the measurement errors are non-negligible. In the case that the measurement of D is precise, the a priori probability density could be written as a Dirac delta function, and the integration could be simplified to Eq. (29). We can see that the probability density function used in early research (Benito et al., 2008) was improper due to the presence of non-negligible uncertainties in both the measured group path data and the transmitting elevation angle data. In the measured data, there are usually a small number of outliers that are difficult to

Fig. 8. The geographical distribution of the ionosonde stations (triangles), the broadcast stations (squares), the receiver (star) and the signals’ reflection areas (circles).

8

D. Yu et al. / Journal of Atmospheric and Solar-Terrestrial Physics 130-131 (2015) 1–13

and the genetic algorithm (Goldberg, 1989) are the two global optimization algorithms most widely used in the geophysical inversion. Warrington, et al. (2009) have tested these two algorithms and found that the genetic algorithm converges faster and provides better results in an ionospheric inversion than the simulated annealing algorithm. Consequently, the genetic algorithm is used to find the optimal ionospheric model parameters in this paper. Previous research (LIU et al., 1994) has shown that the parameters of the E layer in the IRI model are similar to the observed values in China. However, the difference between the parameters of the F2 layer in the IRI model and the observational value is large. Referring to Jiang et al. (2014), the search range for the ionospheric model parameters is set as

Table 4 The location, ground range and azimuth of the broadcast stations. Station code

Location

Latitude (deg)

Longitude (deg)

Ground range (km)

Azimuth (deg)

BJI KUN TWN BEI

Baoji Kunming Taiwan Beijing

34.50 25.17 23.55 39.88

107.17 102.83 120.42 116.57

757.86 1237.30 1021.58 1059.64

47.08 110.18 214.80 340.69

45

⎧ f E −0.5 MHz ≤ f E ≤ f E + 0.5MHz o o IRI ⎪ o IRI ⎪ rbEIRI−5 km ≤ rbE ≤ rbEIRI + 5 km ⎪ ⎪ rmEIRI−5 km ≤ rmE ≤ rmEIRI + 5 km ⎨ ⎪ fo F2IRI−3 MHz ≤ fo F2 ≤ fo F2IRI + 3 MHz ⎪ ⎪ rbF2IRI−50 km ≤ rbF2 ≤ rbF2IRI + 50 km ⎪ r F −50 km ≤ r F ≤ r F + 50 km ⎩ m 2IRI m 2 m 2IRI

Elevation Angle(Degree)

40

35

30

where fo EIRI , rbEIRI , rmEIRI , fo F2IRI , rbF2IRI , and rmF2IRI are the predicated values of IRI2012.

25

20

5. Validation on synthesized data 15 9

10

11

12

13

14

15

16

17

18

Signal Frequency(MHz) Fig. 9. One set of arriving elevation angles for the broadcast signals transmitting from Beijing, extracted from the signals recorded at 0445UT on 31 October 2012.

eliminate but that could lead to an unacceptable inverse result if the Gaussian uncertainty assumption is used. However, an effective method is proposed (Benito et al., 2008) to reduce the impact of the outliers: using a long-tailed probability density (Cauchy distribution) instead of the Gaussian distribution. This approach is taken into account in our method, and the a priori probability densities on the arriving elevation angle data β can be written as

⎛ ⎜ ρβ (β) = k ∏ ⎜ i=1 ⎜ 1 + ⎝ n

1

⎞ ⎟ . 2⎟ ⎟ ⎠

((β −β )/σ ) i

obsi

βobsi

(30)

These uncertainties are presumed to be Cauchy distributed and independent of each other. Substituting Eq. (30) into Eq. (29), we obtain a new expression for σm(m)

⎛ n ⎜ σm(m) = k ∏ ⎜ i=1 ⎜ 1 + ⎝

1

(( β

cali

⎞ ⎟ . 2⎟ ⎟ ⎠

(m,Dobsi)−βobsi)/σ βobsi)

(31)

Eq. (31) is the final expression that we used in our inversion method. As there are no a priori information under considering in this paper, at least 6 elevation data should be used to produce the σm(m) in order to obtain an unique solution. 4.4. Obtain the optimal model parameters By using an optimization algorithm to search in a possible range of the ionospheric model parameters, we can obtain the maximum of the a posteriori probability density function σm(m). The corresponding model parameters m are the optimal parameters and will be output as the inversion result. The simulated annealing algorithm (Kirkpatrick, et al., 1983)

A simulation test was performed to verify the feasibility of this inversion method. The ionospheric model parameters chosen for the simulation were the predicated values from IRI2012 at 0445 UT, 31 October 2012, and the signals were transmitted from the broadcast station in Beijing. The ground range between the transmitter and receiver is 1059.6 km. According to the ITU-HF broadcasting schedule at that time, the HFDF system could receive 21 signals from the Beijing broadcast station on the following frequencies (in MHz): 9.440, 9.610, 9.645, 9.675, 9.685, 11.620, 11.800, 11.905, 11.935, 11.960, 12.045, 15.130, 15.270, 15.380, 15.480, 15.500, 15.710, 17.550, 17.565, 17.605, and 17.625. The radius of the earth is set as re ¼ 6371 km. The ionospheric model parameters were given by IRI2012: fo E = 3.27 MHz , rbE = re + 90 km , rmE = re + 110 km , fo F2 = 10.46 MHz , rbF2 = re + 179.02 km , and rmF2 = re + 268.54 km . After calculating a set of “true” arriving elevation angles by using the method of solving the forward problem described in Section 4.3, we added zero mean Gaussian noise (3° of standard deviation between 7 and 11 MHz, 2° between 11 and 15 MHz and 1° between 15 and 20 MHz, according to the real DOA estimation accuracy of the HFDF system) to simulate a set of measured elevation angle data to test the inversion method, as shown in Fig. 4. Substituting this set of synthesized measured elevation angles as the input data of Eq. (31) and using the GA optimization algorithm, the optimal model parameters can be determined from the range described in Section 4.4 through the use of inversion. The inverted profile and the true profiles are shown in Fig. 5, and the QPS model parameters are listed in Table 1. Under the same conditions, the simulation test was repeated 100 times independently, and the results are shown in Fig. 6. Define the root-mean-square error (RMSE) between the inverted and the true value as

mrmse =

1 N

N

∑ (minverted _ i − mtrue)2 , i=1

(32)

where N is the number of times the test was repeated. The statistical analyses for the 100 simulation tests are shown

D. Yu et al. / Journal of Atmospheric and Solar-Terrestrial Physics 130-131 (2015) 1–13

BJI

9

KUN

13 IRI2012 Estimated Ionostation

12 11

12 foF2(MHz)

10 foF2(MHz)

IRI2012 Estimated Ionostation

14

9 8 7

10

8

6 5

6

4 3

4

6

8

10

12

4

14

4

6

UTC−Time(hour)

8

10

12

14

UTC−Time(hour)

TWN

BEI 12 IRI2012 Estimated Ionostation

14

IRI2012 Estimated Ionostation

11 10

12 foF2(MHz)

foF2(MHz)

9 10

8

8 7 6 5

6

4 4

4

6

8

10

12

14

UTC−Time(hour)

3

4

6

8

10

12

14

UTC−Time(hour)

Fig. 10. The model parameter foF2 for the 4 signal reflection areas estimated by the proposed inversion method (red squares), predicted by IRI2012 (blue line) and measured by ionosonde (black plus line) for the period between 0445UT and 1430UT on 31 October 2012 (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.).

in Table 2. Fig. 5 demonstrates that the electron density profile estimated by the inversion method agrees quite well with the true profile. Similarly, the 2 sets of model parameters listed in Table 1 are very close. These results indicate that the inversion method of using elevation angle data from non-cooperative broadcast stations to infer the ionospheric profile is feasible, and the results of multiple independent tests shown in Fig. 6 and Table 2 confirm this conclusion. In contrast, Table 2 shows that the parameters for the E layer agree better than those of the F2 layer with the true parameters. This result is because the elevation angle of a signal only contains the electron density information of the ionospheric area where it propagates. The higher the area is the less signals can propagate through it, resulting in less electron density information being contained in the elevation angle data. Despite this issue, the profiles of the E layer and F2 layer conform well to the true profile, except for the area very close to the height of rmF2. Therefore, we believe that the real-time ionospheric electron density profile is obtained and can be used to reproduce the propagation characteristics of the true profile. The simulation test to locate an unknown broadcast transmitter by using the estimated parameters is described below. The “unknown” station is set to be Shijiazhuang, located 845.08 km away from the receiver at a similar azimuth (differs by

9°) to the known station, Beijing, where the signals used as input data for the inversion method are launched. The signals launched from the two stations propagate through the same general area of the ionosphere, thus the ionospheric parameters estimated from signals from Beijing can be used to calculate the ground ranges of Shijiazhuang when the elevation angles have been measured. By combining the ground range with the measured azimuth angle, the location of the unknown station can be obtained. According to the ITU-HF broadcasting schedule at test time, the HFDF system could receive 3 signals on different frequencies from Shijiazhuang (in MHz): 11.720, 11.750, and 11.760. A hundred sets of measured elevation angles of these signals are synthesized in the same way as in Beijing. The calculations of the ground ranges using the measured elevations have been performed for the conditions of the true ionospheric profile and the inverted profile, respectively. For comparison, ground ranges have been determined using the true elevations for the condition of the inverted profile. The results are shown in Fig. 7 and Table 3. There are two components contributing to the total location error, one is the uncertainty of the measured elevations and the other is the deviation between the true and the inverted ionospheric profile. It is observed that the location error caused by the inaccurate elevations (as shown in Fig. 7(b)) is much larger than that caused by the inaccurate electron profile (as shown in Fig. 7

10

D. Yu et al. / Journal of Atmospheric and Solar-Terrestrial Physics 130-131 (2015) 1–13

KUN 12

11

11

10

10

9

9 foF2(MHz)

foF2(MHz)

BJI 12

8 7 6

7 6

IRI2012 Estimated Ionostation

5 4 3

8

23

24

25 26 UTC−Time(hour)

27

IRI2012 Estimated Ionostation

5 4 3

28

23

24

12

11

11

10

10

9

9

8 7 6

28

8 7 6

IRI2012 Estimated Ionostation

5 4 3

27

BEI

12

foF2(MHz)

foF2(MHz)

TWN

25 26 UTC−Time(hour)

23

24

25

26

27

IRI2012 Estimated Ionostation

5 4 28

UTC−Time(hour)

3

23

24

25

26

27

28

UTC−Time(hour)

Fig. 11. The model parameter foF2 for the 4 signal reflection areas estimated by the proposed inversion method (red squares), predicted by IRI2012 (blue line) and measured by ionosonde (black plus line) for the period between 31 October 2012, 2300UT and 1 November 2012, 0330UT. (On the x-axis, 25–28 represent 0100UT–0400UT on 1 November.) (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.).

(d)). Consequently, the total location error (as shown in Fig. 7(c)) can be approximated as just the error caused by the inaccurate elevations. That is, although there are obvious differences between the estimated and true profiles in the area that is very close to the height of rmF2, the estimated profiles can approximately reproduce the propagation characteristics of the true profile and could be used in the single-site location (SSL).

6. Results from real data This method has been applied to real data obtained by the HFDF system. The resulting ionospheric model parameters have been compared with the values from IRI2012 and the ionosonde data. To further verify these results, the single-site location of a known broadcast is implemented based on both the estimated parameters and the IRI2012 model parameters. 6.1. The estimated model parameters The elevation angle data used here are derived from 2 sets of signals recorded by the HFDF system near Wuhan. One set is recorded between 31 October 2012, 0445UT and 31 October 2012,

1430UT, and the other set is recorded between 31 October 2012, 2300UT and 1 November 2012, 0330UT. The length of every section of signal is 300 ms and is recorded approximately every half an hour during the 2 recording intervals. In the end, a total of 30 sections of signal are processed in accordance with the flowchart in Fig. 1. Fig. 8 shows the geographical distribution of the ionosonde stations, the broadcast stations, the receiver location and the locations of the signal reflection areas. The specific parameters of the 4 broadcasts are listed in Table 4. Fig. 9 shows an example of the observed elevation angle data, which contains 23 arriving elevation angles from different signals all transmitted from Beijing on 23 different frequencies. The elevation data of each broadcast have been inverted to obtain the ionospheric model parameters of the corresponding reflection area. Meanwhile, the true value of an ionospheric model parameter (foF2) of this area is obtained from the ionosonde stations by using the geostatistics kriging interpolation (LIU et al., 2008). The result allows for the foF2 values estimated by our method, predicted by IRI2012 and measured by ionosonde to be put together, as shown in Figs. 10 and 11. The 3 values of the ionospheric model parameter foF2 obtained from the 3 different data sources are shown in Fig. 10 for times

D. Yu et al. / Journal of Atmospheric and Solar-Terrestrial Physics 130-131 (2015) 1–13

11

1500 True Range Range Located with IRI2012 parameters Range Located with Inverted Parameters

1400 1300

Ground Range(Km)

1200 1100 1000 900 800 700 600 500

0

5

10

15

20

25

30

UTC−Time(hour) Fig. 13. The location results obtained by using the IRI2012 parameters and the estimated parameters.

except that the recording period is from 0700LT to 1200LT, leading to a rising trend. These results verify that this inversion method is feasible and that the estimated parameters are closer to the “true” values than the parameters predicted by the IRI2012 model. 6.2. Locating the SZG station with the estimated parameters

Fig. 12. The geographical locations of the Beijing station (BEI), the Shijiazhuang station (SZG) and the receiver (RT).

varying from 0400UT (1200LT) to 1500UT (2300LT). The figure shows that the values of foF2 at 4 areas from the IRI2012 model and ionosonde are both decreasing, however, with obvious deviations. The estimated values from the proposed method mainly fall on the curve obtained from the ionosonde, which could be used here as the “true” value. Some of the recording time points do not have an estimated value because the station was either not broadcasting at that time or the receiver was not in the coverage area due to a change in the ionospheric environment. Fig. 11 is similar to Fig. 10

The geographical locations of the Beijing station (BEI), the Shijiazhuang station (SZG) and the receiver (RT) are shown in Fig. 12. The signals transmitted from Beijing and Shijiazhuang propagate through approximately the same area of the ionosphere before being received by the receiver near Wuhan. Therefore, the parameters estimated with the elevation data from the Beijing station can be used together with the elevation angle data from Shijiazhuang to locate the Shijiazhuang station. Meanwhile, another location result can be obtained by using the IRI2012 model parameters. Table 5 and Fig. 13 show the location results using the parameters predicted by the IRI2012 model and the parameters estimated by the inversion method. There are 37 effective elevation angles of signals transmitted from SZG on three different frequencies (11.720 MHz, 11.750 MHz and 11.760 MHz), each corresponding to 2 location results obtained by using the IRI2012 parameters and the estimated parameters. Both the mean deviations and the RMSE of the result obtained by using the estimated parameters are better than when using the IRI2012 parameters. The mean error of using the estimated parameters is only 0.40%, which is much better than when using the IRI2012 parameters (
9.80%). However, the difference between the two RMSE values is not so obvious (87.66 km and 105.90 km). This smaller difference is because the RMSE comes mainly from the uncertainties of the elevation angles and not the electron density profile, which is demonstrated in Section 5.

Table 5 The location results obtained on three frequencies by using the IRI2012 parameters and the estimated parameters (The actual ground distance range is 845.08 km). Frequency (MHz)

11.720 11.750 11.760 Total

Occurrence number

14 9 14 37

Using IRI2012 parameters

Using inverted parameters

Mean (km)

Relative error (%)

RMSE (km)

Mean (km)

Relative error (%)

RMSE (km)

735.15 777.90 779.26 762.24


13.01
7.95
7.79
9.80

120.35 87.86 100.95 105.90

822.53 863.58 864.71 848.48


2.67 2.18 2.32 0.40

84.40 92.32 92.32 87.66

12

D. Yu et al. / Journal of Atmospheric and Solar-Terrestrial Physics 130-131 (2015) 1–13

All of the results described above indicate that the detection of the ionosphere using this inversion method is effective and that the resulting ionospheric parameters can be used in SSL to obtain better location results than simply using the IRI model parameters.

7. Conclusion A new method to estimate the bottom side electron density profile of the ionosphere is presented in this paper. By using the commercial HF broadcast transmitters as non-cooperative illuminators, the proposed method does not need any dedicated transmitters. After screening and estimating the DOA of the signals, each group of ground ranges and elevation angles is inverted to obtain a set of QPS model parameters. A probabilistic approach to inverse problems is employed in this inverted method, the following information of which can all be described using probability densities: the information from the measurements of the observable parameters (data), the a priori information on model parameters, and the information on the theoretical correlations between the observable data and model parameters. The inverse problem can then be set as a problem of combining all of this information to obtain a posteriori information on model parameters and find the maximum of the a posteriori probability density function by using the GA algorithm. The proposed method is first validated on synthesized data and then used with real data. The estimated model parameters are much closer to the true values than the values provided by IRI2012. The estimated parameters are then used to locate an HF broadcast. As expected, the method’s performance, which has a mean deviation of 0.40% and an RMSE of 87.66 km, is also much better than using the parameters predicted by IRI2012. The results of simulated and real data verify that the proposed method of obtaining electron density profiles from the elevation angle data of commercial HF broadcast signals is feasible. Although there is an obvious difference between the estimated and true profiles at the height that is very close to rmF2, the estimated profile can be used to reproduce the propagation characteristics of the true profile.

Appendix A. Supplementary material Supplementary data associated with this article can be found in the online version at http://dx.doi.org/10.1016/j.jastp.2015.05.008

References Beley, V.S., Galushko, V.G., Yampolski, Y.M., 1995. Traveling ionospheric disturbance diagnostics using HF signal trajectory parameter variations. Radio Sci. 30 (6), 1739–1752. Benito, E., Bourdillon, A., Saillant, S., Rannou, V., Molinié, J.P., 2008. Inversion of HF backscatter ionograms using elevation scans. J. Atmos. Sol.-Terr. Phys. 70 (15), 1935–1948. Bust, G.S., Mitchell, C.N., 2008. History, current state, and future directions of ionospheric imaging. Rev. Geophys. 46 (1), RG1003. Caratori, J., Goutelard, C., 1997. Derivation of horizontal ionospheric gradients from variable azimuth and elevation backscatter ionograms. Radio Sci. 32 (1), 181–190. Chartier, A.T., Smith, N.D., Mitchell, C.N., Jackson, D.R., Patilongo, P.J.C., 2012. The use of ionosondes in GPS ionospheric tomography at low latitudes. J. Geophys. Res.: Space Phys. 117 (A10), A10326. Chuang, S., Yeh, K., 1977. A method for inverting oblique sounding data in the ionosphere. Radio Sci. 12 (1), 135–140. Croft, T.A., Hoogasian, H., 1968. Exact ray calculations in a quasi-parabolic ionosphere with no magnetic field. Radio Sci. 3, 69. Debao, W., Yunbin, Y., Jikun, O., Kefei, Z., Kai, L., 2008. A hybrid reconstruction algorithm for 3-D ionospheric tomography. IEEE Trans. Geosci. Remote Sens. 46 (6), 1733–1739. Dyson, P., 1991. A simple method of backscatter ionogram analysis. J. Atmos. Terr. Phys. 53 (1–2), 75–88.

Dyson, P., Bennett, J., 1988. A model of the vertical distribution of the electron concentration in the ionosphere and its application to oblique propagation studies. J. Atmos. Terr. Phys. 50 (3), 251–262. Fremouw, E.J., Secan, J.A., Howe, B.M., 1992. Application of stochastic inverse theory to ionospheric tomography. Radio Sci. 27 (5), 721–732. Fridman, O.V., Fridman, S.V., 1994. A method of determining horizontal structure of the ionosphere from backscatter ionograms. J. Atmos. Terr. Phys. 56 (1), 115–131. Fridman, S.V., 1998. Reconstruction of a three-dimensional ionosphere from backscatter and vertical ionograms measured by over-the-horizon radar. Radio Sci. 33 (4), 1159–1171. Fridman, S.V., Nickisch, L.J., 2001. Generalization of ionospheric tomography on diverse data sources: Reconstruction of the three-dimensional ionosphere from simultaneous vertical ionograms, backscatter ionograms, and total electron content data. Radio Sci. 36 (5), 1129–1139. Fridman, S.V., Nickisch, L.J., Hausman, M., 2009. Personal-computer-based system for real-time reconstruction of the three-dimensional ionosphere using data from diverse sources. Radio Sci. 44 (3), RS3008. Galkin, I.A., Reinisch, B.W., Ososkov, G.A., Zaznobina, E.G., Neshyba, S.P., 1996. Feedback neural networks for ARTIST ionogram processing. Radio Sci. 31 (5), 1119–1128. ́ Garcıa-Fernández, M., Hernández-Pajares, M., Juan, J.M., Sanz, J., Orús, R., Coisson, P., Nava, B., Radicella, S.M., 2003. Combining ionosonde with ground GPS data for electron density estimation. J. Atmos. Sol.-Terr. Phys. 65 (6), 683–691. Goldberg, D., 1989. Algorithms in Search, Optimization, and Machine Learning. Addison-Wesley, Reading, MA. Heaton, J.A.T., Cannon, P.S., Rogers, N.C., Mitchell, C.N., Kersley, L., 2001. Validation of electron density profiles derived from oblique ionograms over the United Kingdom. Radio Sci. 36 (5), 1149–1156. Jacquet, F., Caratori, J., Dorizzi, B., 2001. Derivation of the ionospheric structure over extended geographical regions from azimuth and elevation scan backscatter ionograms: a modular neural network approach. Radio Sci. 36 (5), 1043–1052. Jiang, C., Yang, G., Zhao, Z., Zhang, Y., Zhu, P., Sun, H., Zhou, C., 2014. A method for the automatic calculation of electron density profiles from vertical incidence ionograms. J. Atmos. Sol.-Terr. Phys. 107 (0), 20–29. Jowett, A.P., Darnell, M., Riley, N.G., 1989. Passive monitoring of chirpsounder transmitters as an aid to HF frequency management, IEE Colloquium on Adaptive HF Management, pp. 111–115. Kirkpatrick, S., Gelatt, C., Vecchi, M., 1983. Optimization by simulated annealing. Science 220, 671–680. Lind, F.D., Sahr, J.D., Gidner, D.M., 1999. First passive radar observations of auroral E-region irregularities. Geophys. Res. Lett. 26 (14), 2155–2158. LIU, R.-Y., LIU, G.-H., WU, J., ZHANG, B.-C., HUANG, J.-Y., HU, H.-Q., XU, Z.-H., 2008. Ionospheric foF2 reconstruction and its application to the short‐term forecasting in china region. Chin. J. Geophys. 51 (2), 206–213. LIU, R.-Y., QUAN, K.-H., DAI, K.-L., LUO, F.-G., SUN, X.-R., LI, Z.-Q., 1994. A corrected method of the international reference ionosphere to be used in Chinese region. Chin. J. Geophys. 37 (4), 422–432. Ma, X.F., Maruyama, T., Ma, G., Takeda, T., 2005. Three-dimensional ionospheric tomography using observation data of GPS ground receivers and ionosonde by neural network. J. Geophys. Res.: Space Phys. 110 (A5), A05308. Masrani, K., Riley, N.G., 1991. Passive sensing of radio signals in automated HF systems and the updating of off-line prediction models, Seventh International Conference on (IEE) Antennas and Propagation, ICAP 91, 892, pp. 893–897. Meyer, M.G., Sahr, J.D., Morabito, A., 2004. A statistical study of subauroral E-region coherent backscatter observed near 100 MHz with passive radar. J. Geophys. Res.: Space Phys. 109 (A7), A07308. Norman, R., 2003. Backscatter ionogram inversion, Proceedings of the International Radar Conference, 2003, pp. 368–374. Norman, R., Dyson, P., 2006. HF radar backscatter inversion technique. Radio Sci. 41 (4), RS4010. Piggin, P.W., Darnell, M., Gallagher, M., 1996. Passive monitoring for improved HF frequency management, IEE Colloquium on Frequency Selection and Management Techniques for HF Communications, pp. 16/11–16/16. Rao, N.N., 1974. Inversion of sweep-frequency sky-wave backscatter leading edge for quasiparabolic ionospheric layer parameters. Radio Sci. 9 (10), 845–847. Raymund, T.D., Austen, J.R., Franke, S.J., Liu, C.H., Klobuchar, J.A., Stalker, J., 1990. Application of computerized tomography to the investigation of ionospheric structures. Radio Sci. 25 (5), 771–789. Reilly, M.H., 1985. Ionospheric true height profiles from oblique ionograms. Radio Science. 20 (3), 280–286. Reinisch, B.W., Huang, X., 1982. Automatic calculation of electron density profiles from digital ionograms 1. Automatic O and X trace identification for topside ionograms. Radio Sci. 17 (2), 421–434. Reinisch, B.W., Huang, X., 1983. Automatic calculation of electron density profiles from digital ionograms 3. Processing of bottomside ionograms. Radio Sci. 17 (2), 421–434. Ruelle, N., Landeau, T., 1994. Interpretation of elevation-scan HF backscatter data from Losquet Island radar. J. Atmos. Terr. Phys. 56 (1), 103–114. Tarantola, A., 2005. Inverse problem theory and methods for model parameter estimation. Society for Industrial Mathematics. Titheridge, J., 1985. Ionogram analysis with the generalised program POLAN. World Data Center A for Solar-Terrestrial Physics, Boulder, CO (USA), Rep. UAG-93. Tsai, L., Tsai, W., Schreiner, W., Berkey, F., Liu, J., 2001. Comparisons of GPS/MET retrieved ionospheric electron density and ground based ionosonde data. Earth Planets Space 53 (3), 193–205.

D. Yu et al. / Journal of Atmospheric and Solar-Terrestrial Physics 130-131 (2015) 1–13

Warrington, E., Bourdillon, A., Benito, E., Bianchi, C., Monilie, J., Muriuki, M., Pietrella, M., Rannou, M., Rothkaehl, H., Saillant, S., 2009. Aspects of HF radio propagation. Ann. Geophys. 52, 301–321. Zhao, H.S., Xu, Z.W., Wu, J., Wang, Z.G., 2010. Ionospheric tomography by combining vertical and oblique sounding data with TEC retrieved from a tri-band b eacon. J. Geophys. Res.: Space Phys. 115 (A10), A10303.

13

Zheng, H., Ji, G., Wang, G., Zhao, Z., He, S., 2013. Automatic scaling of F layer from ionograms based on image processing and analysis. J. Atmos. Sol.-Terr. Phys. 105–106 (0), 110–118.