1 INTRODUCTION

Neutral wind in the middle and upper atmosphere, animportant area for satellite launch, plays a dominant role in the space environment. Not only is it an important parameter of studying the dynamic processes and energy transfer in the middle and upper atmosphere, but also has a great effect on the orbits of aerospace satellites. Therefore, neutral wind measurements make great contributions to the studying of dynamics, energetics and photochemical coupling of atmosphere for constructing atmospheric global scale models, and provide knowledge to ensure satellite safety and long-life operation.

Currently, Doppler-shift measurements of atmospheric emission lines, such as airglow and aurora, have been used to measure atmospheric wind velocity and direction using ground-based and satellite instruments. Doppler wind measurements from satellites providing global scale coverage have been made using either Fabry-Perot interferometers or stepped Michelson interferometers. The former can measure windfield information based on the Fourier Transform of equal inclination interference formed by Fabry-Perotetalon. The latter, also known as stepped Fourier transform spectrometers (FTS), use a phase shift which can be determined by measuring four interferogram points to retrieve wind velocity. Compared to Stepped Michelson interferometer, Fabry-Perot interferometer has the advantages of higher sensitivity, higher spectral resolving power and simpler structure. Further advantage is that Fabry-Perot interferometers have no moving parts, which results in a higher stability, reliability and sensitivity^[1].

To date, wind observations using Fabry-Perot interferometers from the ground have been made for several decades by many scientists^{[2, 3, 4, 5]} while wind observations from the space dates from 1982, in which Fabry-Perot interferometeron board Dynamics Explorer 2 (DE-2)^{[6, 7, 8, 9]} was launched. It obtained global measurements of the mesospheric and thermospheric wind in the altitude range 100~140 km using the 557.7 nm emission with an accuracy of ~10 m/s. Subsequently, the High Resolution Doppler Imager (HRDI)^{[10, 11, 12, 13]} on board the Upper Atmospheric Research Satellite (UARS) was launched in 1991. It observed the emission and absorption lines of molecular oxygen airglow using low-, medium-, and high-resolution etalons for wind velocity, temperature, and volume emission rate retrieval for investigating the dynamics of the stratosphere (10~40 km), mesosphere and lower thermosphere (50~120 km). The goal of HRDI is to measure the vector winds with an accuracy of 5 m/s. The TIMED Doppler Interferometer (TIDI)^{[14, 15, 16]} has been used for investigating the dynamics and energetics of the Earth’s mesosphere and lower-thermosphere (MLT) since 2001. It supplies vector wind field, temperature and minimum minor constituent of the atmosphere by observing airglow emission lines of 557.7 nm, 630.0 nm, and O₂(0-0) with a vertical resolution of 2 km and accuracies that approach 3 m/s.

In China, wind measurements and studies of the middle and upper atmosphere both from the ground and space have been carried out in recent years. Yuan et al. made continuous measurements of wind velocity and temperature using the Fabry-Perot interferometer^[17] and meteor radar^{[18, 19]} at National Astronomical Observatories (Xinglong site). The retrieval results are in good agreement with Horizontal Wind Model of 1993 (HWM1993). Zhao et al.^[20], Li et al.^[21], and E et al.^[22] studied the data process of wind retrieval using the method of determining complete interference fringe center proposed firstly by Kubota. Wang et al.^[23] derived complete interference fringe center and radius to retrieve wind velocity using the circle fitting with least squares method. Han et al.^{[1, 24, 25]} just studied the method of thinning incomplete interference fringes, but did not describe the method of determining the incomplete fringe center and radius.

To reduce costs of wind measurements, incomplete interference fringes, instead of complete circle, are becoming more common in application (e.g. TIDI). In this work, the center and radius of incomplete interference fringes were determined using the iterative method of Gaussian fitting and Taubincircle fitting. The methodology of wind retrieval is presented in Section 2. In Section 3, we describe the determination of center and radius of interference fringes. In Section 4, the retrieval results and analysis of wind retrieval are given and the conclusions are shown in the last section (Section 5).

2 METHODOLOGY

The offsets of interference fringe peaksformed by FPI from the Doppler shift of the airglow emissions are determined to obtain the wind velocity. After the airglow lights pass through the front optics, they are converted into nearly parallel lights and then undergo multiple reflections in FPI etalon, and finally are focused onto CCD by a focusing device, which forms a group of concentric interference circles.

The transmission function of FPI etalon is

where T is the transmission of etalon, R is reflectivity of inner surface film of etalon, and δ = Φ+2χ is the phase offset with χ phase offset resulting from light reflection in etalon. In Eq.(1), TFPI reaches a maximum value when interferential lights are in the same phase (δ = 2mπ, m, interference orders), while it reaches a minimum value when the phases of interferential lights are inverse (δ = (2m+1)π). For the negligible effect of χ on δ ^[26], δ can be described as where t is the etalon spacing, μ is the refractive index, θ is the reflection angle, and λ is the wavelength of incident lights. In Eq.(2), θ ≈ tan θ = r/f while θ is close to zero. Therefore, Eq.(2) can be transformed by the Taylor series expansion as Where r is the radius of the interference fringe, and f is the focal length. Due to r $ \ll $ f, the transmission function of FPI can also be described as Airy function As described above, the transmitted lights reach maximum when phase offset δ = 2mπ Therefore, interference order m is As the Doppler shift △λ is produced by the incident light wavelength, λ is substituted by λ+△λ and the radius r in Eq.(6) becomes r_λ+△λ Combining Eq.(6) and Eq.(7), the wind velocity can be derived from the offsets of fringe radius based on the shifts of incident light wavelength ($\frac{1}{{\lambda + \Delta \lambda }}{\text{ = }}\frac{1}{\lambda }\left( {1 + \frac{v}{c}} \right)$, here c is light velocity)

3 DATA PROCESSING

Two forms of interferograms imaged on FPI CCD are complete interference fringes (Fig.1) and incomplete interference fringes (Fig.2, fan-shaped in this paper). The center of the former form can be determined by an iterative method based on symmetry coordinates of fringe peaks with the center as symmetry point in horizontal and vertical directions alternately^{[5, 27]}. Then each fringe peak corresponding to radius in the sliced image data is obtained after annular summing by Gaussian function fitting with a nonlinear least-square method^{[28, 29]}. The center of the latter form can be determined by the perpendicular bisectors of uniformly-spaced sampling points atinterference fringe peaks, which, however, is with high randomness of point choice and low precision of point position. Therefore, the method is only applied to the thinning fringes.

In this paper, a method of center and radius determination of the incomplete interference fringes for wind retrieval is proposed. The procedures for data processing are as follows:

(1) Data pre-processing. The noise of the data must be reduced for the center and radius determination using mean filter and wiener weighting adaptive filter.

(2) The center and radius determination of interference fringes. Coordinates of fringe peaksare derived from Gaussian fitting based on the firstly guessed fringe center. Then, new center and radius of fringes are calculated using circle fitting proposed by Taubin^[30].

(3) The final center and radius are determined by using iterative method using Gaussian fitting and circle fitting alternately.

3.1 Data Pre-Processing

Interferograms of airglow formed on CCD are easy to be influenced and contaminated by noise for its low intensity. Besides cosmic rays, the noise, performing as a singular value or sudden change, comes from stray lights, including light leak of instruments and disturbance of environments. It directly affects the precision of center and radius determination. Therefore, data pre-processing must be carried out before wind retrieval including mean filtering^[31] and wiener weighting adaptive filtering^[32]. The intensity distributions of FPI interference fringes without pre-processing and after pre-processing are shown in Fig.3 and Fig.4, respectively. The former distribution is contaminated seriously by the noise, while the latter exhibits the interference intensity distribution well.

3.2 Center and Radius Determination of Interference Fringes

Procedures for determination of interference fringe center and radius are as follows:

(1) The coordinate (C_x, C_y) of the center is estimated firstly by manual operation.

(2) Interference fringes are sliced by several lines with equal angle-interval (Fig.5).

(3) Fringe peaks in the sliced image data are fitted by Gaussian function with a nonlinear least-square method^{[4, 33, 34, 35]}. Gaussian function is given

where r represents the radius of fringe. Constants A, B, and P are determined by the fitting. Then, the coordinates of peaks can be obtained with center coordinate and fringe radius. The lower intensity, the stronger noise effect, the first inner fringe was used for the calculation of first inner fringe is applied to the determination for its strongest intensity and most fitting points.

(4) New center and radius of fringes are calculated using Taubincircle fitting based on the peak coordinates derived in step 3.

(5) The new center coordinates are compared with the original or previous one. If one of the x-axis and y-axis offsets is greater than 0.5 pixels, the procedure goes back to step 2. If both of them are less than 0.5 pixels, the iterative procedure is stopped and the final center and radius of fringes are obtained.

4 ANALOG SIMULATION OF FPI AND RESULTS ANALYSIS OF WIND RETRIEVAL

The FPI system was simulated by ZEMAX using non-sequence mode. The simulation parameters are given in Table1.

In the simulation, the size of CCD is 13.3 mm×13.3 mm with 1024×1024, and pixel size is 13 μm×13 μm, which is used commonly in space instruments. Because of wide existence and stronger intensity in F layer of ionosphere, 630.0 nm airglow is in wide spread use to retrieve wind of the upper atmosphere. Therefore, 630.0nm was chosen as a central wavelength (λ₀) corresponding to 0 m/s. the intensity distribution of 630.0 nm airglow simulation after data pre-processing is shown in Fig.6. Since the common wind of the middle and upper atmosphere has a velocity from several scores of meters per second to hundreds of meters per second, the central wavelength of 629.99979 corresponding to ~100 m/s wind was chosen for simulation (Fig.7). For the purpose of validation of retrieving low wind velocity, the central wavelength of 629.999979 corresponding to ~10 m/s wind was also simulated (Fig.8).

The simulation data of 630.0 nm was processed using the method described in Section 3. The center coordinate of 630.0 nm airglow ZEMAX simulation (λ₀=630.0 nm) is (413.3283, 408.5913), and the radius (r) is 1.6579 mm filled with 127.5308 pixels shown in Fig.6 (red curve). The center coordinate of 629.99979 nm (λ₁=629.99979 nm) is (413.3571, 409.3311), and the radius (r_λ+△λ¹) is 1.7267 mm filled with 132.8231 pixels shown in Fig.7. The center coordinate of 629.999979 nm (λ₁=629.999979 nm) is (411.9925, 412.0979), and the radius (r_λ+△λ²) is 1.6633 mm filled with 127.9429 pixels shown in Fig.8. Then, two pairs of wavelengths (λ₀ and λ₁, λ₀ and λ₂) are substituted in Eq. 1 $\frac{1}{{\lambda + \Delta \lambda }}{\text{ = }}\frac{1}{\lambda }\left( {1 + \frac{v}{c}} \right)$ to derive theoretical wind velocity, while another two pairs of radius (r and r_λ+△λ¹, r and r_λ+△λ²) are substituted in Eq.(8) to obtain retrieval wind velocity. The results are compared in Table2.

The absolute errors (AbsE) of retrieval results are both less than 3 m/s, which meet he requirement of 5 m/s wind retrieval error of the middle and upper atmosphere. The relative errors of retrieval results are -2.98% (row 2 in Table2) and -24.67% (row 3 in Table2), which demonstrate a higher retrieval precision for high wind velocity, but a relative lower retrieval precision for low wind velocity. That may be due to CCD pixel size. According to Eq.(8) and Table1, an error of 1 μm radius determination can lead to 2.34 m/s deviation of wind velocity. Therefore, the pixel size 13 μm is not small enough for low wind retrieval. Therefore, sub-pixel technology will be taken into account in further studies.

5 CONCLUSIONS

In this paper, a new method of determination of center and radius of interference fringes is proposed. The method combines the intensity distribution (Gaussian fitting) and the spatial distribution (Circle fitting) of fringes, which can improve the precision of wind retrieval. After denoising processing of the fringe intensity including mean filtering and adaptive filtering, cross sections are taken for the interference fringes and fringe peak coordinates are acquired according to the Gauss distribution of the intensity. Then, new center coordinates and radius are estimated by using circle fitting with the least square method. Finally, center coordinate and radius can be determined through the iterative operation. The new method was applied to Fabry-Perot interferometer data processing for wind velocity calculation of the middle and upper atmosphere. The errors of the calculation were less than 3 m/s, which demonstrated preliminarily that the method was suitable for the wind retrieval. However, for the purpose of higher retrieval precision of low wind velocity, sub-pixel of CCD image will be taken into account in the future work.

ACKNOWLEDGMENTS

This work was supported by the National Natural Science Foundation of China (41005013).