Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://num-meth.srcc.msu.ru/english/zhurnal/tom_2005/ps/v6r124.ps Äàòà èçìåíåíèÿ: Tue Oct 18 17:53:56 2005 Äàòà èíäåêñèðîâàíèÿ: Mon Oct 1 22:59:52 2012 Êîäèðîâêà:

Numerical Methods and Programming, 2005, Vol. 6 (http://num­meth.srcc.msu.su) 253
UDC 520.88
CORRECTION OF SOLAR DOPPLERGRAMS' LARGE­SCALE FLATFIELD
I. Patrickeyev 1 , G. Mashnich 2 , A. Khlystova 2 , H. Zhang 3
The authors have proposed and compared three different methods of large­scale flatfield correction
(time averaging, Fourier filtering, and wavelet­based filtering). As an example, the real dopplergram of
the velocity field in a solar active region obtained at the Huairou Solar Observing Station of National
Astronomical Observatories of China has been processed. It is shown that any of the applied methods
can improve the visual quality of the image. A quantitative comparison of the proposed methods was
performed by evaluating the heterogeneity of the images.
1. Introduction. Modern astronomical observations widely use detectors based on a CCD­matrix. There
are many papers which discuss not only the features of CCD­detectors, but also the means of overcoming their
instrumental distortions. The most essential intensity variations introduced by measuring systems are caused
by the following factors: the non­uniform sensitivity of matrix pixels, the non­uniform light field because of the
optical irising, the interference or the dust particles on the camera lenses, and the filter defects. Moreover, the
intensity varies from frame to frame due to the variations of the diffuse light and the visibility conditions. In
the observations of the elongated bright objects such as the Sun, the images obtained outside the focal plane
are often used as a flatfield. In this case, the flatfield and the observation data are obtained at the different
configurations of an optical system. A calibration method of the spatial non­uniformity of the matrix is proposed
in [1]. In this paper, the construction of a flatfield is carried out as follows: the authors first determine the pixel
gain coefficients based on the images with known shift and then calculate the flatfield for the whole image by
an iteration technique.
Our methods of 2D FeI –5324.19 š A and Hfi Dopplergram processing have been used at the Huairou Solar
Observing Station [2, 3]. It was suggested that a flatfield can be described by a surface of the second order. Better
results were obtained by subtracting the flatfield of the image calculated by extrapolating and interpolating the
flatfield for a limited number of points of the unprocessed image obtained after large­scale averaging. A filtering
of the radio maps with the discrete Daubechi wavelets was discussed in [4] and an advantage of the wavelet
application over the traditional filtering technique was shown. An estimation of the ff­effect based on the velocity
field observations and its reliability was discussed in [5].
The main aim of this paper is to show the possibility of the flatfield correction in the obtained time series
of the images and to propose other methods of the large­scale flatfield correction based on the truncated Fourier
series and the expansion in terms of the continuous wavelet basis.
2. Methods of large­scale flatfield correction. Unfortunately, at present it is impossible to use the
flatfield obtained directly in an observation (the so­called ``zero frame''), because this will require updating of
the software installed on board the observing station. All the methods proposed in this paper do not use a zero
frame and are based on different hypotheses of the temporal or the spatial structure of the flatfield. Moreover,
all methods use the assumption that the observed image is formed by combining a large­scale flatfield with the
ideal velocity field. Therefore, we take into account an additive component only. However, it should be noted
that the proposed methods can also be adopted for a multiplicative component correction (by converting the
intensity to the logarithmic scale).
2.1. Time Averaging. Let us assume that the flatfield is a stationary structure. The definition ``stationary''
has a special meaning in the theory of the stochastic processes. In this paper we use this term in a common sense,
i.e., ``stationary'' means time invariable. This assumption looks reasonable at least in the case of an instrumental
flatfield. If the velocity field varies with time, then the flatfield and the velocity fields can be separated on the
basis of a different behavior of the fields in time. The one way to identify stationary structures is to collect frames
of the same active region and then to perform averaging. Such a mean frame involves structures reproduced
1 Institute of Continuum Mechanics, Russian Academy of Sciences, Korolyova 1, 614013, Perm, Russia; e­mail:
pat@icmm.ru
2 Institute of Solar and Terrestrial Physics, Russian Academy of Sciences, Lermontova 126, 664033, Irkutsk,
Russia; e­mail: mashnich@iszf.irk.ru, ancha­k@nm.ru
3 National Astronomical Observatories, Chinese Academy of Sciences, 100012, Beijing, China; e­mail:
hzhang@bao.ac.cn

254 Numerical Methods and Programming, 2005, Vol. 6 (http://num­meth.srcc.msu.su)
in each frame as well as the time averaged non­stationary structures. Therefore, if the velocity at a point is
governed by a random time function with zero mean value, then the point in the mean frame will have the zero
velocity value.
The mean frame can be found by the following formula: f i;j = 1
N
N \Gamma1
X
n=0
s (n)
i;j , where N is the number of the
frames, i; j are the indices of a pixel in the x­ and y­axes, respectively, s i;j are the frames, and (n) is a frame
number.
The corrected velocity field is obtained by subtraction of the mean frame from the observed velocity frame.
A degree of separation of the velocity field and flatfield increases with the number of frames; however, the
observation time cannot be too long. In practice, acceptable results are achieved even at N = 15 -- 20.
2.2. Fourier Filtering. The application of the Fourier transform to the flatfield filtering [6] relies on the
assumption that the flatfield and the velocity field have different spectral features. Since the flatfield has a
low­frequency spectrum and the velocity field spectrum involves mainly the middle and high frequencies, it is
possible to separate them in the frequency domain. The visual analysis of the frames shows that the typical
flatfield appears as a structure consisting of a few large spots or of the two or three wide strips extending across
the whole frame area. To delete the low­frequent structure, it is necessary to apply the 2D Fourier transform
and to replace the Fourier coefficients at the lowest frequencies by zeros. The acceptable results are achieved
after deleting the 2 or 3 lowest frequencies.
The 2D Fourier coefficients are obtained from the formula
bs kx ;ky = 1
N x N y
Nx \Gamma1
X
i=0
Ny \Gamma1
X
j=0
s i;j e \Gamma2ú
p
\Gamma1 (kxi=Nx+ky racj=Ny ) ;
where N x and N y are the x­ and y­image dimensions in the physical space; k x and k y are the coordinates in the
frequency domain.
2.3. Wavelet Filtering. A continuous wavelet transform [7] allows us to obtain almost the same result as
the Fourier transform. Following the initial assumption, the flatfield can be considered as a large­scale structure,
whereas the velocity field --- as a combination of the middle­scale and the small­scale structures. The 2D isotropic
wavelet transform coefficients are obtained according to the formula
Fig. 1. The original velocity field
e
s a;i 0 ;j 0 = a \Gamma1
Nx \Gamma1
X
i=0
Ny \Gamma1
X
j=0
s i;j / \Lambda
i i \Gamma i 0
a
;
j \Gamma j 0
a
j
;
where a, i 0 , and j 0 are the parameters of the scaling and the translation
along the x­ and y­axes, respectively, / is a wavelet, and \Lambda indicates complex
conjugation.
A flatfield image is obtained by the synthesis (inverse wavelet transform)
on the largest scales. The subtraction of the flatfield from the original image
yields the filtered image. In this paper, the Mexican Hat wavelet is used:
/(x; y) = (2 \Gamma x 2 \Gamma y 2 ) exp
\Gamma
\Gamma(x 2 + y 2 )=2
\Delta
:
Here the wavelet coefficients are calculated as a linear convolution of the original image by the wavelet
series whose scale varies from 1 pixel to the image size. The synthesis is made using the Morlet formula [8]. In
contrast to the Fourier­based method with a sharp frequency cut­off, the wavelet­based method gives smooth
borders of the frequency bandwidth [9].
3. Numerical results. All the methods described above were applied to the real Dopplergrams obtained
at the Huairou Solar Observing Station.
Each Dopplergram is obtained by the formula
V dop ¸
256
X
n=1
I red \Gamma
256
X
n=1
I blue
256
X
n=1
I red +
256
X
n=1
I blue
;
where I red and I blue are the measured intensities in the red and blue wings of the spectral line. Each pair is
separated in time by 20 ms.The time period of collection is about 40 s.

Numerical Methods and Programming, 2005, Vol. 6 (http://num­meth.srcc.msu.su) 255
Figure 1 illustrates a typical velocity field obtained at the Hfi wavelength – equal to 486.1 nm. Some
results of the large­scale correction are presented in Figures 2a, 3a, and 4a. The flatfields obtained by different
methods are given in Figures 2b, 3b, and 4b. Each filtered image is calculated as a result of the subtraction of
the corresponding flatfield from the velocity field 1, whereas each filtered image looks more or less homogeneous,
without large­scale structures.
Fig. 2. Time averaging: a) filtered image; b) flatfield
Fig. 3. Filtering in the frequency domain: a) filtered image; b) flatfield obtained with FFT
The results obtained by the mean frame method are presented in Figure 2. A filtered image (Figure 2a)
was obtained by subtracting the mean frame from the original image. The corresponding flatfield (Figure 2) is
the result of the averaging of 23 frames over a period of 9 minutes. Our comparison with the original image
shows that the flatfield distortions become almost invisible.
Figure 3 illustrates the results of the image processing with the Fourier transform. The large­scale structures
were filtered almost completely out; however, the flatfield image obtained with the FFT (Figure 3b) differs from
the mean frame (Figure 2b).
Fig. 4. Wavelet filtering: a) filtered image; b) flatfield obtained with wavelets
Some results of the wavelet­based filtering is presented in Figure 4. Similarly to the previous cases, the
filtering yields a visible improvement of the image quality, which can be seen as a more uniform distribution
of the intensity over the frame field. The flatfield image obtained by a synthesis of the large­scale wavelets
(Figure 4b) seems to be the smoothest one compared to the other flatfield images.
A simple quantitative estimation of the image homogeneity can be done by comparing the mean values in
the upper and lower halves of the image. It is expected that the perfectly homogeneous image has the same
mean value in both parts. In order to make a quantitative comparison of the above methods, the calculated
mean values of each half­image were normalized to the mean value of the whole image. The table shows the
deviation of the calculated values from the mean value of the whole image for the velocity field (Figure 1) and

256 Numerical Methods and Programming, 2005, Vol. 6 (http://num­meth.srcc.msu.su)
for each of the filtered images. The upper (lower) half­image of the original velocity field has the mean values
of about 5% smaller (larger) than the mean value of the whole image. The halves of the filtered images have
the mean values weakly deviating in the range from 2.3 % to 0.03%. Thus, each of the proposed methods can
reduce the image heterogeneity by subtracting the corresponding flatfield.
Deviations (in %) of the mean value calculated in halves of the image from
the mean value of the whole image
Image Upper half­image Lower half­image
Velocity field (Fig. 1) ­5.5 5.1
Time averaging (Fig. 2) ­0.1 ­0.4
Fourier filtering (Fig. 3) ­0.03 ­0.6
Wavelet filtering (Fig. 4) ­2.3 1.8
Actually, it is possible to apply more sophisticated methods of the estimation of the image homogeneity
(the calculation of the correlation radius, etc.). Perhaps, we do not need to use such methods if one can say
that the images have no large­scale structures. Indeed, the aim of the proposed methods is making an image
free of the instrumental artifacts (not making them completely homogeneous). As a result, now we have three
similar images; each image looks adequately enough for the further analysis of the velocity distribution.
4. Conclusion. From the visual estimation of the filtered images it may be concluded that the acceptable
results can be obtained using any of the proposed methods. The time averaging method yields acceptable results
in studying the structure evolution over the time intervals comparable with the observation time of a single
image series. The advantage of the Fourier­based and th wavelet­based methods is the possibility of processing
a single image, whereas the mean frame method requires a series of the images. If a zero frame is known, an
optimized filter should be used; this gives a further improvement of the Fourier­based method. The wavelet­
based filtering can also be improved by choosing the other wavelets for the analysis and the synthesis. For
example, the directional structures can be better filtered with the anisotropic wavelets [10].
Acknowledgements.The authors would like to thank V. V. Grechnev, D.D. Sokoloff, P. G. Frick, K. M. Ku­
zanyan for the fruitful discussions. This work was partially supported by the collaborative grants of RFBR and
NSFC 02--02--39027, RFBR and DFG 03--02--04031, the grant RFBR 03--02--16384, the State Support of Leading
Scientific Schools of the Russia Federation, the grant SS--733.2003.2, and the Exchange Program for Scientists'
visits between the Russian Academy of Sciences and the Chinese Academy of Sciences. The Russian authors
would like to thank the Huairou Solar Observing Station of the National Astronomical Observatories of China
for hospitality.
References
1. Kuhn J., Lin H., Loranz D. Gain calibrating nonuniform image­array data using only the image data // Publications
of the Astronomical Society of the Pacific. 1991. 103. 1097--1108.
2. Ai G., Li W., Zhang H. Formation of the FeI 5324.19 line in the Sun and theoretical calibration of solar magnetic
telescope // Chinese Journal of Astronomy and Astrophysics. 1992. 6. 129--136.
3. Zhang H., Ai G. The formation of the H­Beta line in the solar chromosphere magnetic field // Acta Astronomica
Sinica. 1986. 27. 286.
4. Chen P., Zhang X., Xiang S., Feng L., Reich W. Protruding structure buried in radio map by wavelet // Chinese
Physics Letters. 2000. 17, N 5. 388--389.
5. Zhang H., Bao S., Kuzanyan K. Twist of magnetic fields in solar active regions // Astronomy Reports. 2002. 46,
N 5. 424--434.
6. Bracewell R. Fourier transform and its application. New York: McGrow­Hill, 1965.
7. Farge M., Hunt J.C.R., Vassilicos J.C. Wavelets, fractals and Fourier transforms. Oxford: Clarendon Press, 1993.
8. Daubechies I. Ten lectures on wavelets. Philadelphia: SIAM, 1992.
9. Frick P., Beck R., Berkhuijsen E., Patrickeyev I. Scaling and correlation analysis of galactic images // Monthly
Notices of the Royal Astronomical Society. 2001. 327. 1145--1157.
10.Patrickeyev I., Fletcher A., Beck R., Berkhuijsen E., Frick P., Horrelou C. Anisotropic wavelet analysis of spiral arms
and magnetic fields in the galaxy M51 // The Magnetized Plasma in Galaxy Evolution (Proceedings of the Conference
held in Cracow, Poland, Sept. 27th ­ Oct. 1st, 2004. Edited by K. Chyzy, K. Otmianowska­Mazur, M. Soida, and
R.­J. Dettmar). Cracow, 2005. 130--135.
Received 09 September 2005