A method for the automatic calculation of electron density profiles from vertical incidence ionograms

Journal of Atmospheric and Solar-Terrestrial Physics 107 (2014) 20–29

Contents lists available at ScienceDirect

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

A method for the automatic calculation of electron density profiles from vertical incidence ionograms Chunhua Jiang, Guobin Yang n, Zhengyu Zhao, Yuannong Zhang, Peng Zhu, Hengqing Sun, Chen Zhou Ionospheric Laboratory, School of Electronic Information, Wuhan University, 430079 Wuhan, Hubei, China

art ic l e i nf o

a b s t r a c t

Article history: Received 1 July 2013 Received in revised form 17 October 2013 Accepted 22 October 2013 Available online 29 October 2013

Vertical incidence ionograms indicate ionospheric characteristics over the ionosonde station, from which electron density profiles can be derived. This paper describes a method for the automatic calculation of electron density profiles from vertical incidence ionograms. First, the method calculates the initial parameters of the quasi-parabolic segments (QPS) model by using the International Reference Ionosphere (IRI) model, the Nequick2 model, image processing techniques, and the Empirical Orthogonal Function (EOF). Once the initial parameters have been calculated the method then adjusts those to obtain the electron density profiles and synthesized traces that match the recorded ionograms. The algorithm then selects the best-fit synthesized trace and corresponding parameters as output from the candidate ones. Furthermore, the corresponding electron density profile of the recorded ionogram is calculated by using the best-fit parameters of the QPS model. To further test the feasibility of the proposed method, we apply it to some ionograms that were recorded at Wuhan during different seasons. As a result, our results demonstrate that the proposed method is feasible for the automatic calculation of electron density profiles. & 2013 Elsevier Ltd. All rights reserved.

Keywords: Ionogram Automatic scaling Electron density profile Quasi-parabolic segments

1. Introduction With the development of modern advanced ionospheric sounders, the automatic processing of digital ionograms, which is essential for research in ionosphere characteristics and space weather applications, has gathered increasing attention. The automatic processing technique consists of two steps: the first part is the automatic scaling of the ionograms; and the second part is the automatic calculation of electron density profiles from ionograms. Wuhan University developed a versatile advanced ionospheric sounder (Chen et al., 2009) that can carry out both backscatter and vertical incidence soundings. At the moment, the ionosonde is not equipped with an automatic scaling system that can output in near real-time both the main ionospheric characteristics and the vertical electron density profile from the vertical incidence ionograms. Therefore, a method is proposed for the automatic calculation of electron density profiles from vertical incidence ionograms in this paper. At present, there are two methods that are widely used to calculate the electron density profile from the vertical ionogram: the NhPC program, included in the ARTIST system, was developed by University of Lowell (Reinisch and Xueqin, 1983), and it can automatically calculate electron density profiles for Digisonde

n

Corresponding author. Tel.: þ 8613407195656. E-mail address: [email protected] (G. Yang).

1364-6826/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jastp.2013.10.012

ionograms; the Polan program, developed by Titheridge (1985), can calculate electron density profiles using some required input parameters. Recently, Zabotin et al. (2006) developed an iterative ray-tracing approach, is named NeXtYZ, to recover the parameters of a quite sophisticated three-dimensional model of the local electron density distribution. The INGV (National Institute of Geophysics and Vulcanology) developed an automatic processing software named Autoscala (Scotto and Pezzopane, 2002; Pezzopane and Scotto, 2005, 2008) that can automatically scale ionograms using an image processing technique characterized by several filters (Scotto and Pezzopane, 2008; Pezzopane and Scotto, 2010). In 2009, Scotto (2009) added a routine for the real-time computation of electron density profiles. The model uses the Reinisch and Huang formulation (Reinisch and Huang, 2000) for the bottom-side F2 profile and the F1 layer and uses the polynomial anchor points for the E valley and the E bottom-side profile. Autoscala adjusts the model parameters, which are 12 free parameters, to fit a recorded ionogram and then selects the best fit as the output result. Vesnin et al. (2011) presented a method for calculating electron density profiles from vertical sounding ionograms based on modifying the initial profile, which can be obtained from the IRI model prediction or from the previous sounding of the ionosphere. This paper presents a method for the automatic calculation of electron density profiles from vertical incidence ionograms. An ionosphere model named QPS (Croft and Hoogasian, 1969; Dyson and Bennett, 1988), that is widely used to calculate the electron density

C. Jiang et al. / Journal of Atmospheric and Solar-Terrestrial Physics 107 (2014) 20–29

800

700

600 Virtual height (km)

profile from ionogram, is applied to this method for the calculation of electron density profile. Firstly, we determine the initial frequencies for the ionospheric layers from vertical incidence ionograms by using the IRI model, the Nequick2 model and image processing technique. Then the proposed method determines the initial heights for the ionospheric layers by selecting the best-fit trace from the candidate traces that are synthesized by the IRI model. To reduce the calculation time, the EOF is used to reduce the size of the candidate traces. Once the initial parameters of the QPS model have been determined, the method fine-tunes those parameters to obtain the best-fit parameters of the QPS model based on the fit match result between the recorded ionogram and the synthesized trace. Furthermore, it can calculate the corresponding electron density profile using the best-fit parameters of the QPS model. In the following paragraphs, the present study describes in detail the procedures of determining the initial parameters and the corresponding electron density profile. Moreover, the present method is applied to the ionograms recorded at Wuhan for testing its feasibility.

21

500

400

300

200

100 2

4

6

8

10

12

14

Frequency (MHz)

Fig. 1. A typical synthesized trace.

2. Quasi-parabolic segments The quasi-parabolic segments model is widely used for inversion of ionograms. Much work (Rao, 1974; Benito et al., 2008; Song et al., 2011) has been devoted to the inversion of backscatter ionograms. Norman (2003) introduced the QPS model for the inversion of vertical incidence ionograms. As mentioned above, previous studies mostly calculated the differences between the simulated and the recorded traces to obtain the best-fit parameters, thus it must guarantee that the recorded trace is accurately identified. To avoid the identification of the recorded trace for the recorded ionogram, the proposed method uses the complete set of data to calculate the fit match between the synthesized traces and recorded ionogram. The QPS model consists of five segments that describe the ionosphere layers (E, F1, and F2) and the joining layers. Each of the ionosphere layers (E, F1, and F2) is characterized by three parameters: the critical frequency f c of the layer, the peak height hm of the layer, and the semi-thickness ym of the layer. The joining layer is then determined using the parameters of the two corresponding layers. In other words, nine parameters for the QPS model will be selected in this study to calculate the electron density profile. Once electron density profiles have been calculated by using the QPS model, the proposed method then synthesizes the candidate traces using the well-known formulae (1) and (2). Furthermore, we can obtain the best-fit synthesized trace and corresponding parameters from the candidate traces. Z hr h′ðf Þ ¼ u′ dh ð1Þ hb

1 u′ ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 1  f p =f

f oF2 ¼ 11:8ðMHzÞ; ym F2 ¼ 90ðkmÞ; hm F2 ¼ 280ðkmÞ

3. The calculation of the electron density profiles The further the initial parameters are away from the best-fit parameters, the longer the calculation time is. Therefore, to reduce the calculation time for the present method, it is important for determining the initial parameters. This section introduces how to determine the initial parameters of the QPS model for recorded ionograms. Once the initial parameters have been determined, the proposed method fine-tunes the initial parameters to obtain the best-fit synthesized trace and corresponding electron density profile for recorded ionograms. 3.1. The initial parameters of the QPS model 3.1.1. Initial values for the critical frequency of the ionosphere layers 3.1.1.1. E layer. Previous research (Liu et al., 1994) has shown that the critical frequency of the E layer of the IRI model is similar to the observed values in China. In this paper, we assign the value of the IRI model to the initial critical frequency of the E layer. The initial frequency of the E layer is represented by the following formula: f oEinitial ¼ f oEmodel

ð3Þ

ð2Þ

where f is the working frequency, hb is the base of the ionosphere, hr is the height of the reflecting position, u′ is the group refraction index, and f p is the plasma frequency. In the procedure of calculation of the trace for the electron density profile, the simple numerical integration can cause the instability. Both Titheridge (1985, 1988) and Scotto et al. (2012) developed algorithms to address this problem. In the present work, the proposed method uses the algorithm developed by Scotto et al. (2012) to remove the phenomenon. Fig. 1 shows a typical synthesized trace, where the parameters of the ionospheric layers are represented by the following equations: f oE ¼ 4:2ðMHzÞ; ym E ¼ 20ðkmÞ; hm E ¼ 110ðkmÞ

f oF1 ¼ 5:6ðMHzÞ; ym F1 ¼ 80ðkmÞ; hm F1 ¼ 190ðkmÞ

3.1.1.2. F1 layer. Based on the results reported by Liu et al. (1994), the predicted critical frequency of the F1 layer of the IRI model is close to the observed values in China. The mean variation of f oF1 is approximately 0–0.1 MHz, and the standard error of f oF1 is approximately 0.3–0.6 MHz. As a result, it is suitable to define the predicted values of the IRI model as the initial value of the critical frequency of the F1 layer. However, the occurrence probability of the F1 layer does not agree well with the observed values in China. Thus, this work addressed this deficiency of the IRI model by using the Nequick2 model (Nava et al., 2008), which is mainly used for GPS applications associated with the ionospheric time delay. If the F1 layer is present in the IRI model, the corresponding value in the IRI model is defined as the initial critical frequency of

22

C. Jiang et al. / Journal of Atmospheric and Solar-Terrestrial Physics 107 (2014) 20–29

the F1 layer: f oF1initial ¼ f oF1model

ð4Þ

If the F1 layer is not present in the IRI model, the critical frequency of the F1 layer is represented by formula (5), which was described in detail by Nava et al. (2008): 8 if f oEmodel Z 2 > < f oF1initial ¼ 1:4nf oEmodel if f oEmodel o 2 f oF1initial ¼ 0 > : f oF1 ¼ 0:85n 1:4nf oE if 1:4f oE 40:85nf oF2 initial

model

model

initial

Defining a processing window

Noise and dilating processing for the recorded ionograms

Calculating the projection values for the recorded iongrams

ð5Þ 3.1.1.3. F2 layer. In China, the difference between the critical frequency of the F2 layer in the IRI model and that in the observation value is large. Thus, this work determines the initial value of f oF2 by using an image processing technique named image projection (Sun et al., 2000). The image projection projects the image data, which are represented by the two-dimensional data, to the one-dimensional coordinate to get the one-dimensional image data which represent the characteristics of the image. Generally, the one-dimensional coordinate is the x axis or the y axis that depends on which one characteristics of the image we want to obtain. In the present study, the proposed method requires the characteristics of the frequency to determine the initial values of the critical frequency for the recorded ionogram. Hence, it projects the recorded ionograms to the axis of the frequency. The recorded ionogram can be represented by a two-dimensional matrix AðM; NÞ , where M is the number of the frequency points and N is the number of the virtual height points. Furthermore, the projection values of the frequency for the recorded ionogram can be calculated by formula (6). N

FðiÞ ¼ ∑ Aði; jÞ; where i o ¼ M; jo ¼ N j¼1

ð6Þ

In order to avoid the noise interference and the disconnection of the trace of the F2 layer in the recorded ionograms, the proposed method need to some pre-processing techniques including median filter and image dilating algorithms. Moreover, to eliminate the echo of the Es layer that affects the projection values of the recorded ionograms on the frequency axis, it defines a processing window to which the projection algorithm is applied. The proposed method defines the processing window as W½M; N, whose numbers are represented by the following formulae: M ¼ int ½ðΔf w Þ=Δf  þ 1

ð7Þ 0

0 0 N ¼ int ½ðhmax  hmin Þ=Δh  þ 1;

ð8Þ

where Δf w is the horizontal size of the processing window, Δf is 0 the resolution of the frequency, hmax is the maximum height of the 0 processing window, hmin is the minimum height of the processing 0 window, and Δh is the resolution of the virtual height. The present method defines Δf w as the width of the working 0 frequency (2–15 MHz). The value of hmax varies from 300 km to 0 600 km, and hmin varies from 200 km to 300 km. In addition, the 0 0 variables hmax and hmin are dynamic with the seasonal and diurnal changes. Fig. 2 describes the flowchart of the steps used to determine the initial value of the critical frequency for the F2 layer, and Nth is defined as 1 in this paper. Fig. 3 illustrates the characteristics of the projection on the frequency for the typical synthesized trace including the ordinary and extraordinary waves. From the illustrations of the Fig. 3, we find that the first rightmost peak value represents the position of f xF2 and the second rightmost peak

Filtering the projection values and finding the position (Pos_1th) of the rightmost peak value. Then calculating the number Npeaks of the projection peak values at the range (Pos_1th-1.5MHz, Pos_1th)

NO Npeaks>Nth?

YES Selecting the two rightmost peak values as the critical values of the F2 layer

The output of the electronic density profile is NA

Fig. 2. Flowchart of the steps used to establish the initial value for the F2 critical frequency. N th is the threshold value of the number of projection peak values within the range (Pos_1th–1.5 MHz, Pos_1th).

value represents the position of f oF2. In the projection values, because of the contribution of the trace of the extraordinary wave, the projection value of the position of f oF2 is mostly higher than the X model value. Fig. 4 illustrates the characteristics of the projection on the frequency for the recorded ionogram, in which case the Npeaks value is higher than the Nth value. Fig. 5 illustrates the characteristics of the projection on the frequency for the recorded ionogram, in which case the Npeaks value is no higher than the Nth value. As mentioned above, the projection algorithm can effectively determine the initial critical frequency f oF2initial for the F2 layer.

3.1.2. Initial values of the height of the ionosphere layers To determine the initial values of the heights for the ionosphere layers, the proposed method selects the initial heights for the ionospheric layers from the candidate sample parameters that are built in this section. At the beginning, the presented work builds the candidate sample parameters of the QPS model based on the IRI predicted values. In order to represent completely the characteristics of the height for the ionosphere as much as possible, this method obtains the candidate sample parameters in which the year varies from 1993 to 2003 over Wuhan. However, the size of the candidate sample parameters from the IRI model is large. To reduce the calculation time, it reduces the size of the candidate sample parameters using the EOF (Huang, 2004). The EOF is an empirical model that is used to analyze the main characteristics of the history sample data. (Dvinskikh 1988; Dvinskikh and Naidenova, 1991) firstly introduced the EOF to the ionospheric research. Recently, the EOF is widely used to research the characteristics of the ionosphere (Wang et al., 2004; Mao et al., 2005; Ding et al., 2007). When the candidate sample parameters of the QPS model have been built, the proposed method can calculate electron density

C. Jiang et al. / Journal of Atmospheric and Solar-Terrestrial Physics 107 (2014) 20–29

800

800

700

700

600

600

23

280

500

400

500

400

300

200

200

4

265

255

100 2

fxF2

270

260

300

100

275

Projection value

Virtual height (km)

Virtual height (km)

foF2

6 8 10 12 14 Frequency (MHz)

2

4

250

6 8 10 12 14 Frequency (MHz)

2

4 6 8 10 12 14 Frequency (MHz)

Fig. 3. An example of the calculation of the projection values for the synthesized trace: (a) the synthesized trace including the ordinary (the green curve) and extraordinary (the red curve) waves; (b) the synthesized trace is represented by the two-dimensional matrix and applied by the image dilating algorithm; (c) the projection values of the frequency for the synthesized trace. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

SNR(dB) 800

800 70

700 60

600

50

500

40

400

30

300

20

200 100 2

4

6

8

10

12

14

Virtual height (km)

Virtual height (km)

700

600 500 400 300

10

200

0

100 2

4

6

800

140

700

130

600

Projection value

Virtual height (km)

8

10

12

14

Frequency (MHz)

Frequency (MHz)

500 400 300

foF2 fxF2 at Pos_1th

120 110

1.5MHz 100 90

200 100

80 2

4

6

8

10

12

14

Frequency (MHz)

2

4

6

8

10

12

14

Frequency (MHz)

Fig. 4. An example of the calculation of the projection values for the recorded ionogram that the condition ‘Npeaks4Nth' is met: (a) the original ionogram, (b) the corresponding binary ionogram that is processed by the median filter, (c) the binary ionogram processed by the image dilating algorithm, (d) the projection values of the frequency for the recorded ionogram.

profiles, and then synthesize the corresponding candidate traces using formulae (1) and (2). In this section, this method uses the projection algorithm, which was used to determine the initial values

of the critical frequency for the F2 layer, to obtain the characteristics of the height for the synthesized traces and the recorded ionograms by projecting the data to the axis of the virtual height. When the

24

C. Jiang et al. / Journal of Atmospheric and Solar-Terrestrial Physics 107 (2014) 20–29

SNR(dB) 800

800 80

700

700 600

60

500

50

400

40 30

300

Virtual height (km)

Virtual height (km)

70 600 500 400 300

20 200

200

10

100 2

4

6

8

10

12

100

0

14

2

4

6

8

10

12

14

Frequency (MHz)

Frequency (MHz)

100

800

1.5MHz Pos_1th

95

600

Projection value

Virtual height (km)

700

500 400 300

90 85 80 75

200 100

70 2

4

6

8

10

12

14

2

4

6

8

10

12

14

Frequency (MHz)

Frequency (MHz)

Fig. 5. An example of the calculation of the projection values for the recorded ionogram that the condition ‘Npeaks4Nth' is not met: (a) the original ionogram; (b)the corresponding binary ionogram that is processed by the median filter; (c) the binary ionogram processed by the image dilating algorithm; (d) the projection values of the frequency for the recorded ionogram.

projection values of the synthesized traces and the recorded ionograms have been calculated, this method selects the best-fit height parameters from the candidate sample parameters as the initial heights based on the match results between the projections values for the synthesized traces and that for the recorded ionograms. In order to build the candidate sample parameters ðf o; hm ; ym Þ of the QPS model for the ionospheric layers, we use formulae (9)–(11) to represent the parameters of the ionospheric layers. The parameters of the QPS model for the E layer are defined as f oE ¼ f oEIRI hm E ¼ hm EIRI ym E ¼ hm EIRI  90

N

ZðtÞ ¼ X þ ∑ ai ðtÞ  Ei n¼1

ð9Þ

The parameters of the QPS model for the F1 layer are defined as f oF1 ¼ f oF1IRI hm F1 ¼ hm F1IRI ym F1 ¼ hm F1IRI  hm EIRI

ð10Þ

The parameters of the QPS model for the F2 layer are defined as f oF2 ¼ f oF2IRI hm F2 ¼ hm F2IRI ym F2 ¼ B0IRI ;

ð11Þ

In formulae (9)–(11), the script IRI means that the value is from the IRI model. The candidate sample parameters can be represented by the one-dimensional variable XðtÞ: XðtÞ ¼ ½f oEðtÞ; hm EðtÞ; ym EðtÞ; f oF1ðtÞ; hm F1ðtÞ; ym F1ðtÞ; f oF2ðtÞ; hm F2ðtÞ; ym F2ðtÞ;

where f oEðtÞ, hm EðtÞ, and ym EðtÞ are the parameters for the E layer; f oF1ðtÞ, hm F1ðtÞ, and ym F1ðtÞ are the parameters for the F1 layer; f oF2ðtÞ, hm F2ðtÞ, and ym F2ðtÞ are the parameters for the F2 layer; and t is the time variable, which varies from 1993 to 2003 with a 1-h step size. To reduce the size of the candidate sample parameters XðtÞ, the EOF is applied to the candidate sample parameters XðtÞ and it can build a linear expression for the candidate sample parameters XðtÞ. Thus the smaller candidate parameters can be represented by formula (13).

ð12Þ

ð13Þ

where ZðtÞ is the smaller sample parameters, X is the average of XðtÞ, Nis the number of quantities of XðtÞ, ai ðtÞ is the ith-order time index of the orthogonal function, and Ei is the ith-order orthogonal function. Based on previous studies (Wang et al., 2004; Mao et al., 2005; Ding et al., 2007), this method uses the first four orders to represent the main characteristics of the candidate sample parameters. By changing the first four orders of the time index (a1 ¼ minða1 ðtÞÞ  maxða1 ðtÞÞ, a2 ¼ minða2 ðtÞÞ  maxða2 ðtÞÞ,a3 ¼ minða3 ðtÞÞ  max ða3 ðtÞÞ, and a4 ¼ minða4 ðtÞÞ  maxða4 ðtÞÞ) of the orthogonal function, the proposed method obtains the new dataset ZðtÞ of the candidate sample parameters. Using the different steps of the time index, the method can obtain the sample parameters ZðtÞ with different sizes. In the later procedure, this method uses the candidate sample parameters ZðtÞ instead of the dataset XðtÞ. When the candidate sample parameters ZðtÞ have been built, we can calculate the

C. Jiang et al. / Journal of Atmospheric and Solar-Terrestrial Physics 107 (2014) 20–29

corresponding candidate traces. Both the synthesized traces and the recorded ionograms can be represented by a two-dimensional matrix, the proposed method then defines the recorded ionograms as AðM; NÞand the synthesized trace as A′ðM; NÞ, where M is the number of the frequency points and N is the number of the virtual height points. The projection values of the recorded ionogram and the synthesized trace can be obtained by formulae (14) and (15), respectively. M

FðjÞ ¼ ∑ Aði; jÞ; where i o ¼ M; j o ¼ N

ð14Þ

i¼1 M

F′ðjÞ ¼ ∑ A′ði; jÞ; where i o ¼ M; j o ¼ N

ð15Þ

i¼1

25

3.2. Selecting the best-fit parameters by adjusting the initial parameters When the initial parameters of the QPS model have been determined, they are ðf oEinitial ; hm Einitial ; ym Einitial ; f oF1initial ; hm F1initial ; ym F1initial ; f oF2initial ; hm F2initial ; ym F2initial Þ, and close to the values of the best-fit parameters. The primary task discussed in this section is the adjustment of the parameters within a small range. By the adjustment of the initial parameters, the proposed method can obtain the best-fit synthesized trace. Furthermore, it can calculate the corresponding electron density profile for the recorded ionogram by using the corresponding best-fit parameters. To calculate the fit quality between the synthesized trace

When the projection values have been calculated, the proposed method obtains the best-fit parameters for the recorded ionograms based on the correlation values between F and F′ using formula (16).

800

700

n

C ¼ ∑ F′ðjÞnFðjÞ

600 Virtual height (km)

j ¼ n0

n ¼ h′max =Δh′ n0 ¼ h′min =Δh′;

ð16Þ

where h′max varies from 300 km to 600 km with the seasonal and diurnal changes, h′min is defined as 90 km, and Δh′ is the resolution of the virtual height. In this section, the proposed method can obtain the best-fit sample parameters ðhm Einitial ; ym Einitial ; hm F1initial ; ym F1initial ; hm F2initial ; ym F2initial Þ from the candidate sample parameters ZðtÞ. Fig. 6 illustrates the matched result with the image projection technique and shows that the projection is able to identify satisfactory initial values for the heights of the ionosphere layers. The green curve shown in Fig. 7 is the best-fit synthesized trace, which was obtained by adjusting the initial parameters, and the red curve is the corresponding electron density profile which represents the best estimation of the electron density profile associated to the recorded ionogram.

500

400

300

200

100 2

4

6

8

10

12

14

Frequency (MHz)

Fig. 7. The best-fit synthesized trace with the best-fit parameters, which were obtained by adjusting the initial parameters. The green curve indicates the synthesized trace, and the red curve is the corresponding electronic density profile that was automatically calculated using the proposed method. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

0.9

800

0.8

700

0.7

Normalized value

Virtual height (km)

600

500

400

0.6 0.5 0.4 0.3

300 0.2 200 0.1 100 2

4

6

8

10

Frequency (MHz)

12

14

0 100

200

300

400

500

600

Virtual height (km)

Fig. 6. (a) A recorded ionogram at Wuhan with the best-fit initial sample parameters of the height for ionospheric layers (the green curve); (b) the projection values for the recorded ionogram (the blue curve) and the synthesized trace (the red curve). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

26

C. Jiang et al. / Journal of Atmospheric and Solar-Terrestrial Physics 107 (2014) 20–29

and recorded ionogram, the proposed method transforms the twodimensional matrixes AðM; NÞ and A′ðM; NÞdescribed in formulae (14) and (15) to binary matrixes AðM; NÞbin andA′ðM; NÞbin , and

calculates the fit quality using formula (17). C¼ t¼

Fine-tuning the initial parameters of the QPS model for the ionospheric layers

T¼

Calculating the electronic density profiles and the synthesized traces for the adjusted parameters

Calculating the fit quality between the synthesized traces and the recorded ionograms

Selecting the highest fit quality value Cmax

NO Cmax>Cth YES Outputting the best-fit parameters and the corresponding electron density profile

Outputting the initial parameters and the corresponding electron density profile

m ¼ M;n ¼ N

∑

m ¼ 1;n ¼ 1

m ¼ M;n ¼ N

800

700

700

600

600

400 300

∑

m ¼ 1;n ¼ 1

800

500

Aðm; nÞbin ¼ A′ðm; nÞbin ; where A′ðm; nÞbin ¼ 1 A′ðm; nÞbin ¼ 1

200

500 400 300 200

100

100 2

4

6

8

10

12

14

2

4

6

Frequency (MHz)

8

10

12

14

12

14

Frequency (MHz)

800

800

700

700

600

600

Virtual height (km)

Virtual height (km)

ð17Þ

where t is the number points of the fit values and T is the total number of the synthesized trace. Fig. 8 illustrates a flowchart of the steps used to select best-fit parameters by adjusting the parameters for the ionosphere layers, where C th is defined as 50%. Because the bottom layers (E and F1) can affect the echo traces of the radio reflection of the upper layers (F1 and F2), the proposed method firstly adjusts the parameters for the E layer, then that for the F1 layer, and lastly that for the F2 layer to obtain the best-fit parameters for the ionosphere layers. To obtain the best-fit trace of the E layer, the proposed method addresses the following two cases. If the F1 layer is present, hm E is varied from hm Einitial to hm Einitial þ 20 with a 10-km step, ym E is defined as hm E 90, and f oE is varied from f oEinitial 0:5 to f oEinitial þ 1 with a 0.1-MHz step. If the F1 layer is not present, hm E is varied from hm Einitial tohm Einitial þ100, and the other parameters are the same as in the first case. To obtain the best-fit parameters for the F1 layer, hm F1 is varied from hm F1initial  50 to hm F1initial þ50 with a 10-km step, ym F1 is varied from ym F1initial  50 to ym F1initial þ 50 with a 10-km step, and f oF1 is varied from f oF1initial  0:5 to f oF1initial þ 0:5 with a 0.1-MHz step. Lastly, to obtain the best-fit parameters for the F2 layer, hm F2 is varied from hm F2initial  50 to hm F2initial þ50 with a 10-km step, ym F2 is varied from ym F2initial  50 to ym F2initial þ 50 with a 10-km step, and f oF2 is varied from f oF2initial  0:3 to f oF2initial þ 0:2 with a 0.1-MHz step.

Virtual height (km)

Virtual height (km)

Fig. 8. Flowchart of the steps used to select the best-fit parameters by adjusting the initial parameters for the ionospheric layers, and C th is the threshold value of the fit match.

t  100% T

500 400 300 200

500 400 300 200

100

100 2

4

6

8

10

Frequency (MHz)

12

14

2

4

6

8

10

Frequency (MHz)

Fig. 9. The best-fit Ionograms without the echoes of the F1 layer recorded at Wuhan: (a) on 28 June 2012 at 19:00 LT; (b) on 14 August 2012 at 21:05 LT; (c) on 1 September 2012 at 22:40 LT; (d) on 17 December 2012 at 03:34 LT.

800

800

700

700

600

600

Virtual height (km)

Virtual height (km)

C. Jiang et al. / Journal of Atmospheric and Solar-Terrestrial Physics 107 (2014) 20–29

500 400 300

500 400 300 200

200

100

100 2

4

6

8

10

12

2

14

4

6

8

10

12

14

12

14

Frequency (MHz)

Frequency (MHz)

800

800

700

700

600

600

Virtual height (km)

Virtual height (km)

27

500 400 300

500 400 300 200

200

100

100 2

4

6

8

10

12

14

2

4

6

8

10

Frequency (MHz)

Frequency (MHz)

Fig. 10. The best-fit Ionograms with the echoes of the F1 layer recorded at Wuhan: (a) on 17 May 2012 at 11:55 LT; (b) on 27 June 2012 at 13:55 LT; (c) on 14 August 2012 at 13:10 LT; (d) on 16 December 2012 at 09:55 LT.

4. Testing the method with recorded ionograms

800

We have determined the initial parameters of the QPS model for the ionospheric layers, which effectively reduce the calculation time for obtaining the best-fit parameters. The average time required to complete the procedure for the automatic calculation of electron density profiles does not exceed 60 s on a PC equipped with a Core i3 CPU, a working frequency of 2.5 GHz and a memory capacity of 2 GB. To verify the feasibility of the proposed method that is used to automatically calculate electron density profiles, some ionograms recorded at Wuhan in winter, summer, and equinoctial months are investigated. In this test, the ionograms are classified into three categories: the first category includes the recorded ionograms without the F1 layer, the second category includes the recorded ionograms with the F1 layer, and the third category includes the fit mismatches between the recorded ionograms and the synthesized traces. Fig. 9 shows the recorded ionograms belonged to the first category, and these recorded ionograms were recorded in summer, winter, and equinoctial months. The results shown in Fig. 9a, c and d demonstrate that the present method performs well for the analysis of disconnected ionograms. Fig. 10 shows the recorded ionograms belonged to the second category and these ionograms were recorded in summer and equinoctial months. Although the echo of the F1 layer is weak so that it is difficult to be observed (Fig. 10b), the proposed method can automatically calculate the electron density profile. From the calculated electron density profile shown in Fig. 10b, the peak height of the maximum electron density is not always defined as the virtual height at the frequency about 0:834nf oF2, and Piggott and Rawer (1972) discussed the difference. Because of the model of the F1 layer, the proposed method does not agree well with the echo of the F1 layer with a ledge F1 layer (Fig. 10d). Baker (1990) and

700

Virtual height (km)

600 500 400 300 200 100 2

4

6

8

10

12

14

Frequency (MHz) Fig. 11. The fit mismatch Ionogram with the echoes of the F0.5 layer recorded at Wuhan on 30 July 2012 at 10:45LT.

Baltazart and Wilkinson (1995) discussed this problem, and the future works will focus on improving the performance of the model of the F1 layer. The present method is not designed to identify the additional layers (such as F0.5 layer, F1.5 layer or F3 layer), and then it produces the fit mismatches for that case. Fig. 11 describes the phenomenon with the F0.5 layer. Fig. 12 illustrates its application during the spread-F stage at the daytime (Fig. 12a) and the nighttime (Fig. 12b), the results indicate that it is performed not well. In addition, because the multiple reflection echo of the Es layer can affect the initial values of the critical frequency for the F2 layer and the fit match between the recorded ionogram and the

C. Jiang et al. / Journal of Atmospheric and Solar-Terrestrial Physics 107 (2014) 20–29

800

800

700

700

600

600

Virtual height (km)

Virtual height (km)

28

500 400 300

500 400 300 200

200

100

100 2

4

6

8

10

12

14

Frequency (MHz)

2

4

6

8

10

12

14

Frequency (MHz)

Fig. 12. The fit mismatch Ionograms during the spread-F stage recorded at Wuhan: (a) on 22 December 2012 at 15:30 LT; (b) on 17 December 2012 at 05:12LT.

synthesized trace, we do not show the samples in that case the proposed method is invalidated. As a result, the proposed method should be carefully applied in these cases mentioned above.

5. Conclusions This paper describes a method for the automatic calculation of electron density profiles from vertical incidence ionograms. This method determines the initial parameters using the International Reference Ionosphere model, the Nequick2 model, the image processing technique and the EOF. Furthermore, it fine-tunes the initial parameters to obtain the best-fit parameters for the recorded ionograms. In addition, this work uses the QPS model, which has been proven to be reliable and robust for the calculation of the electron density profile by many studies, to calculate the electron density profiles. Although the proposed method is performed not well for some cases (the additional layer, spread-F and multiple Es), our results indicate that it is feasible for the fit match between the recorded ionograms and the synthesized traces, and that inspires us to develop and improve it. Thus, the proposed method still requires some adjustments to improve its accuracy and performance, and future studies will focus on the application of ionograms recorded at different geographic locations and during the spread-F stage, the application of ionograms that are affected by the multiple echoes of the Es layer, and the improvement of the performance of the model of the F1 layer with a ledge F1 layer.

Acknowledgments This work was supported by the National Natural Science Foundation of China (NSFC grant no. 41304127 and NSFC no. 41204111). The authors are grateful to two anonymous reviewers for their assistance in evaluating this paper. References Baker, D.C., 1990. Extension of the multi-parabolic ionospheric model to the F1 layer L condition. Electron. Lett. Commun. 26, 1301–1303. Baltazart, V.P., Wilkinson, 1995. F1 layer modeling of ionospheric electron density distribution. Electron. Lett. Commun. 26, 1816–1817. Benito, E., Bourdillon, A., Saillant, S., Rannou, V., Molinie, J.P., 2008. Inversion of HF backscatter ionograms using elevation scans. J. Atmos. Sol.-Terr. Phys. 70 (15), 1935–1948. Chen, G., Zhao, Z., Li, S., Shi, S., 2009. WIOBSS: The Chinese low-power digital ionosonde for ionospheric backscattering detection. Adv. Space Res. 43 (9), 1343–1348. Croft., T.A., Hoogasian, H., 1969. Exact ray calculations in a quasi-parabolic ionosphere. Radio Sci. 3 (1969), 69–74.

Ding, Z.H., Ning, B.Q., Wan, W.X., 2007. Real-time automatic scaling and analysis of ionospheric ionogram parameters. Chin. J. Geophys. 50, 837–847, http://dx.doi. org/10.1002/cjg2.1101. Dvinskikh, N.I., 1988. Expansion of ionospheric characteristics fields in empirical orthogonal functions. Adv. Space Res. 8 (4), 179–187, http://dx.doi.org/10.1016/ 0273-1177(88)90238-4. Dvinskikh, N.I., Naidenova, M.Ya., 1991. An adaptable regional empirical ionospheric model. Adv. Space Res. 11 (10), 7–10, http://dx.doi.org/10.1016/02731177(91)90312-8. Dyson, P.L., Bennett, J.A., 1988. A model of the vertical distribution of the electron concentration in the ionosphere and its application to oblique propagation studies. J. Atmos. Terres. Phys. 50 (3), 251–262. Huang, J.Y., 2004. Method of Meteorological Statistical Analysis and Forecast, 3th edition Meteorological Press, Beijing, China, pp. 170–197. (In Chinese). 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. V37 (04), 422–432. Mao, T., WAN, W., Xing, L., Li, B., 2005. An EOF based empirical model of TEC over Wuhan[J]. Chin. J. Geophys. V48 (4), 751–758. Nava, B., Coisson, P., Radicella, S.M., 2008. A new version of the NeQuick ionosphere electron density model. J. Atmos. Sol.-Terr. Phys. 70 (15), 1856–1862. Norman R.J., 2003. An Inversion Technique for obtaining Quasi-Parabolic layer parameters from VI Ionograms. Radar Conference. Proceedings of the International. IEEE Press, Australia. Pezzopane, M., Scotto, C., 2005. The INGV software for the automatic scaling of foF2 and MUF(3000)F2 from ionograms: a comparison with ARTIST 4.01 from Rome data. J. Atmos. Sol. Terr. Phys. 67 (12), 1063–1073. Pezzopane, M., Scotto, C., 2008. A method for automatic scaling of F1 critical frequencies from ionograms. Radio Sci. 43, RS2S91, http://dx.doi.org/10.1029/ 2007RS003723. Pezzopane, M., Scotto, C., 2010. Highlighting the F2 trace on an ionogram to improve autoscala performance. Comput. Geosci. 36 (9), 1168–1177, http://dx. doi.org/10.1016/j.cageo.2010.01.010. Piggott, W.R., Rawer, K., 1972. U.R.S.I. Handbook of Ionogram Interpretation and Reduction. World Data Center A for Solar-Terrestrial Physics NOAA. Boulder, Colorado, pp. 17–18. Rao, N.N., 1974. Inversion of sweep-frequency sky-wave backscatter leading edge for quasi-parabolic ionospheric layer parameters. Radio Sci. 9 (10), 845–847, http://dx.doi.org/10.1029/RS009i010p00845. Reinisch, B.W., Xueqin, H., 1983. Automatic calculation of electron density profiles from digital ionograms: 3. Processing of bottomside ionograms. Radio Sci. 18 (3), 477–492, http://dx.doi.org/10.1029/RS018i003p00477. Reinisch, B.W., Huang, X., 2000. Redefining the IRI F1 layer profile. Adv. Space Res. 25 (1), 81–88. Scotto, C., 2009. Electron density profile calculation technique for autoscala ionogram analysis. Adv. Space Res. 44 (6), 756–766, http://dx.doi.org/10.1016/ j.asr.2009.04.037. Scotto, C., Pezzopane, M., 2002. A software for automatic scaling of foF2 and MUF (3000) F2 from ionograms. In: Proceedings of URSI 2002. Maastricht, August 2002, pp. 17–24. Scotto, C., Pezzopane, M., 2008. Removing multiple reflections from the F2 layer to improve Autoscala performance. J. Atmos. Sol. Terr. Phys. 70 (15), 1929–1934, http://dx.doi.org/10.1016/j.jastp.2008.05.012. Scotto, C., Pezzopane, M., Zolesi, B., 2012. Estimating the vertical electron density profile from an ionogram: On the passage from true to virtual heights via the target function method. Radio Sci. 47, RS1007, http://dx.doi.org/10.1029/ 2011RS004833. Song, J., Zhao, Z.Y., Zhou, C., Chen, G., 2011. Inversion of HF sweep-frequency backscatter ionograms. Chin. J. Geophys. V54 (8), 1953–1959, http://dx.doi.org/ 10.3969/j.issn.0001-5733.2011.08.002. (in Chinese).

C. Jiang et al. / Journal of Atmospheric and Solar-Terrestrial Physics 107 (2014) 20–29

Sun, Y., Zhou, G.H., Zhou, L.C., Shi, P.F., 2000. Fast template matching algorithm based on the projection. J. Shanghai Jiaotong Univ. V34 (5), 702–704. (cnki: ISSN:31-1466.0.2000-05-049). Titheridge, J.E., 1985. Ionogram analysis with the generalised program Polan, Report UAG-93. World Data Center A for Solar-Terrestrial Physics. Boulder, Colorado. Titheridge, J.E., 1988. The real height analysis of ionograms: A generalized formulation. Radio Sci. 23 (5), 831–849, http://dx.doi.org/10.1029/RS023i005p00831. Vesnin A.M., Ratovsky K.G., Klimenko M.V., Klimenko V.V., 2011. Vertical Sounding Ionogram Technique based on Modification of Initial Electron Density profile.13th

29

Ionospheric Effects Symp. IES2011. Alexandria VA 22308-1943. May 17–19, 2011: abstracts. 2011. P. A028. Wang, X YX.Y., Ji, S. YS.Y., Xi, D LD.L., 2004. EoFs analytic representation of the ionospheric electron density profile. Chin. J. Radio Sci. 19 (5), 560–564. (in Chinese). Zabotin, N.A., Wright, J.W., Zhbankov, G.A., 2006. NeXtYZ: Three-dimensional electron density inversion for dynasonde ionograms. Radio Sci. 41, RS6S32, http://dx.doi.org/10.1029/2005RS003352.