Extraction of Quasi-Monochromatic Gravity Waves from an Airglow Imager Network

An algorithm has been developed to isolate the gravity waves (GWs) of different scales from airglow images. Based on the discrete wavelet transform, the images are decomposed and then reconstructed in a series of mutually orthogonal spaces, each of which takes a Daubechies (db) wavelet of a certain scale as a basis vector. The GWs in the original airglow image are stripped to the peeled image reconstructed in each space, and the scale of wave patterns in a peeled image corresponds to the scale of the db wavelet as a basis vector. In each reconstructed image, the extracted GW is quasi-monochromatic. An adaptive band-pass filter is applied to enhance the GW structures. From an ensembled airglow image with a coverage of 2100 km × 1200 km using an all-sky airglow imager (ASAI) network, the quasi-monochromatic wave patterns are extracted using this algorithm. GWs range from ripples with short wavelength of 20 km to medium-scale GWs with a wavelength of 590 km. The images are denoised, and the propagating characteristics of GWs with different wavelengths are derived separately.

Since Peterson and Adams [12] utilized the all-sky airglow imagers (ASAIs) to observe GWs from induced perturbations of airglow, ASAIs have been widely used to observe GWs. Based on the temporal and spatial information provided by ASAIs, quasi-monochromatic GWs can be studied to reveal the relationship between the wave scale and the propagation characteristics. Due to the limited field of view of ASAIs, most previous studies have focused on the GWs of wavelengths shorter than 100 km [13]. On the other hand, some studies have focused on small-scale ripples from broken GWs, revealing the instability of the local atmosphere [14]. Airglow sensors onboard low-earth orbit satellites

GW Events
Two GW events (shown in Figure 1) can be spotted easily at the same time on the night of 21 December 2014 (UTC): one propagated from northeast to southwest, lasting from 15:20 p.m. to 18:46 p.m. (UTC), and another wave propagated from north to south, lasting from 17:36 p.m. to 18:37 p.m. (UTC). An assembled image of 2100 pixels × 1200 pixels (one pixel denotes one kilometer) is shown in Figure 1, exhibiting an area of 2100 km × 1200 km. The blank area at the right lower corner is due to the breakdown of an ASAI at Rongcheng Station. No convective weather can be found around the center derived from the waveform. As thunderstorms are very rare in winter in Northern China, this GW event was possibly triggered by a frontal system or jet [6]. The airglow image matrix is 2100 km × 1200 km. The all-sky airglow imager (ASAI) stations are marked with yellow stars. The two GW events are marked by orange and green, respectively. The blank area at the right lower corner is due to the breakdown of the ASAI at Rongcheng Station. The left lower corner is out of the effective radius of Xinxiang and Shuozhou station, so that the area is blank.

Discrete Wavelet Transform (DWT)
DWT is an ideal algorithm to separate waves of different scales from the raw GWs with wide spectrum by peeling the wave patterns from a large to small scale. DWT divides the function space of object function F(x) into a series of mutually orthogonal subspaces. Information of different scales in F(x) are reconstructed and expressed separately in these subspaces with corresponding basis vectors. The function space L 2 (R) of F(x) is a square integrable [27] and can be decomposed as follows: The airglow image matrix is 2100 km × 1200 km. The all-sky airglow imager (ASAI) stations are marked with yellow stars. The two GW events are marked by orange and green, respectively. The blank area at the right lower corner is due to the breakdown of the ASAI at Rongcheng Station. The left lower corner is out of the effective radius of Xinxiang and Shuozhou station, so that the area is blank.

Discrete Wavelet Transform (DWT)
DWT is an ideal algorithm to separate waves of different scales from the raw GWs with wide spectrum by peeling the wave patterns from a large to small scale. DWT divides the function space of object function F(x) into a series of mutually orthogonal subspaces. Information of different scales in F(x) are reconstructed and expressed separately in these subspaces with corresponding basis vectors. The function space L 2 (R) of F(x) is a square integrable [27] and can be decomposed as follows: where the direct sum symbol ⊕ means the union of spaces. V indicates the scale space, used to construct an approximation to F(x). W is the wavelet space, used to construct waves representing the difference between the approximation and actual function. j is the scale factor. According to Equation (1), F(x) is divided into an approximation with lower frequency and wave of higher frequency. Subsequently, in the same way, the approximation with a lower frequency in scale space V j can be further decomposed into a new approximation of a lower frequency in spaces V j−1 and a new wave of higher frequency in W j−1 [27]: Equation (1) can then be expanded as follows: where all of the subspaces of V 0 and W j are orthogonal to each other. W j0 is the wavelet subspace with the smallest scale s 0 and the finest resolution. The scale of a basis vector of W j is 2 j 0 −j s 0 . In order of scale from small to large, information is peeled from F(x) and reconstructed in the wavelet subspace W j of the corresponding scale. The component f j (x) of F(x) in the scale space V j can be denoted by a linear combination: where a is the scale coefficient. The scale function φ(x) is the basis vector of V j : where k is the translation factor. W j is the wavelet subspace, with wavelet function ψ j,k (x) as its basis vector. Component g j (x) of F(x) in the wavelet space W j can be also decomposed into a linear combination: where b is the wavelet coefficient and In our work, Daubechies 45 (db 45) wavelet [28] was adopted in DWT; 45 indicates the maximum vanishing moment allowed in the MATLAB function library. A large vanishing moment leads to a smooth fitting, but at the cost of poor compact support in the frequency domain and more computing cost. In the db 45 wavelet, neither φ(x) in Equation (5) nor ψ(x) in Equation (7) had an analytic expression (shown in Figure 2).  Through DWT, F(x) is expanded as the sum of the linear combinations in the scale subspace and a number of wavelet subspaces: In short, Equation (8) is rewritten as follows: where f0 is the approximation to In the image processing, the two-dimensional image matrix is decomposed by DWT along three directions-vertical, horizontal, and diagonal. Then, the two-dimensional DWT expansion of different scales is reconstructed based on the DWT expansions along three directions [27]. The reconstructed two-dimensional DWT expansions are a series of images peeled from the original image, with both vertical and horizontal basis vector scales at 0 0 2 j j s − pixels. In the following image processing, s0 is 2 pixels. We used n (n = j0 -j + 1) to define a peeled layer. The wave patterns reconstructed in the peeled image at the nth layer with a basis vector at the scale of 2 n pixels is a quasi-monochromatic wave with a wavelength between the 2 n and 2 n+1 pixels. However, considering that the wavelet function (shown in Figure 2b) with a large vanishing moment is not a pulse, the wavelength in the nth layer may slightly exceed the range of 2 n to 2 n+1 , especially when n is small.

Enhancing the Wave Pattern
The wave patterns in a peeled image are limited in a narrow spectrum, so that they can be enhanced by a scan kernel with band-pass filtering. The kernel is a square with a side length of triple Through DWT, F(x) is expanded as the sum of the linear combinations in the scale subspace and a number of wavelet subspaces: where In short, Equation (8) is rewritten as follows: where f 0 is the approximation to F(x) by scale function with basis vector φ(x) at scale 2 j 0 s 0 , and g j is the DWT expansion with basis vector ψ(x) at scale 2 j 0 − j s 0 .
In the image processing, the two-dimensional image matrix is decomposed by DWT along three directions-vertical, horizontal, and diagonal. Then, the two-dimensional DWT expansion of different scales is reconstructed based on the DWT expansions along three directions [27]. The reconstructed two-dimensional DWT expansions are a series of images peeled from the original image, with both vertical and horizontal basis vector scales at 2 j 0 − j s 0 pixels. In the following image processing, s 0 is 2 pixels. We used n (n = j 0 -j + 1) to define a peeled layer. The wave patterns reconstructed in the peeled image at the nth layer with a basis vector at the scale of 2 n pixels is a quasi-monochromatic wave with a wavelength between the 2 n and 2 n+1 pixels. However, considering that the wavelet function (shown in Figure 2b) with a large vanishing moment is not a pulse, the wavelength in the nth layer may slightly exceed the range of 2 n to 2 n+1 , especially when n is small.

Enhancing the Wave Pattern
The wave patterns in a peeled image are limited in a narrow spectrum, so that they can be enhanced by a scan kernel with band-pass filtering. The kernel is a square with a side length of triple the basis vector of the layered image, going through the image with a stride of the basis vector. For example, a window of 48 pixels × 48 pixels slides on the image of the fourth layer (n = 4) with Atmosphere 2020, 11, 615 6 of 14 the stride of 16 pixels. In each step, the region covered by the filter (Figure 3a) was transferred to the frequency domain by two-dimensional fast Fourier transform, so that two symmetric frequency peaks could be revealed. Then, the spectrum of the covered region ( Figure 3c) was filtered by multiplying the following filtering function (shown in Figure 4): where u and v are the zonal and meridional frequency, respectively. u max and v max are the frequencies of the maximum intensity. Due to the symmetry of spectrum, two peak points can be spotted at (u max , v max ) and (−u max , −v max ), respectively. After frequency filtering, the noise was removed in the frequency domain. Then, the product (Figure 3d) of the original spectrum and H was transferred back to the image domain by two-dimensional inverse fast Fourier transform. Figure 3b shows the wave patterns became more visible after the wave enhancement.
Atmosphere 2020, 11, x FOR PEER REVIEW 6 of 14 the basis vector of the layered image, going through the image with a stride of the basis vector. For example, a window of 48 pixels × 48 pixels slides on the image of the fourth layer (n = 4) with the stride of 16 pixels. In each step, the region covered by the filter (Figure 3a) was transferred to the frequency domain by two-dimensional fast Fourier transform, so that two symmetric frequency peaks could be revealed. Then, the spectrum of the covered region ( Figure 3c) was filtered by multiplying the following filtering function (shown in Figure 4): where u and v are the zonal and meridional frequency, respectively. umax and vmax are the frequencies of the maximum intensity. Due to the symmetry of spectrum, two peak points can be spotted at (umax, vmax) and (−umax, −vmax), respectively. After frequency filtering, the noise was removed in the frequency domain. Then, the product ( Figure  3d) of the original spectrum and H was transferred back to the image domain by two-dimensional inverse fast Fourier transform. Figure 3b shows the wave patterns became more visible after the wave enhancement.

Multiple Scale Analysis of GW Events
The OH airglow image ( Figure 1) captured by the ASAI network contained GWs with a broad spectrum, which could be isolated separately through DWT. Considering that the assembled airglow image is 2100 pixels × 1200 pixels, the max n was set to 9. In other words, the largest scale of the basis vector was 2 9 . Figure 5 shows a series of layered images reconstructed with a basis vector at the scale of 2 n (n = 1, 2, …, 9) pixels, with one pixel representing one kilometer.
In each peeled image, the dominant wave pattern was much stronger than the noise caused by the background, clouds, or GWs of other scales, so that denoising could be performed by enhancing the strongest spectrum. For better visibility, all of the peeled images in Figure 5 were processed by the filtering kernel mentioned in Section 3.3 in order to reduce the noise and thus enhanced the wave patterns. Considering the spectrum peak of the various GWs in different regions of the peeled image, a local spectrum peak centered band-pass filter is required. An adaptive filtering kernel was applied for denoising, scanning through the peeled image region by region and making band-pass filtering with the peak of each local spectrum.
Waves with horizontal wavelengths of less than 12 km ere too weak to be observed by ASAI [29], so the first three layers (n = 1, 2, 3) were superposed together (Figure 5a), containing noises and a few vague ripples. Ripples with wavelengths of 25 ± 2 km, possibly from breaking GWs, were marked in Figure 5b. These ripples were perpendicular to the local GWs shown in the following sub-figures (Figure 5c-f), indicating dynamical instability [14].
The peeled images of the fifth to ninth layers are represented in Figure 5c-g, respectively. The GWs at different scales are illustrated in the corresponding images and their wavelengths are listed in Table 2. The medium-scale GW of the 591 ± 75 km wavelength (Figure 5g) was rarely found on ground observation [18], because its wavelength was beyond the field of view of a single ASAI and can be only captured by the ASAI network. By extracting the medium-scale GW from various GW spectrums, DWT offers a new way to measure the parameters of medium-scale GW, especially the temporal information absent in low-earth orbit satellite observations.

Multiple Scale Analysis of GW Events
The OH airglow image (Figure 1) captured by the ASAI network contained GWs with a broad spectrum, which could be isolated separately through DWT. Considering that the assembled airglow image is 2100 pixels × 1200 pixels, the max n was set to 9. In other words, the largest scale of the basis vector was 2 9 . Figure 5 shows a series of layered images reconstructed with a basis vector at the scale of 2 n (n = 1, 2, . . . , 9) pixels, with one pixel representing one kilometer.
In each peeled image, the dominant wave pattern was much stronger than the noise caused by the background, clouds, or GWs of other scales, so that denoising could be performed by enhancing the strongest spectrum. For better visibility, all of the peeled images in Figure 5 were processed by the filtering kernel mentioned in Section 3.3 in order to reduce the noise and thus enhanced the wave patterns. Considering the spectrum peak of the various GWs in different regions of the peeled image, a local spectrum peak centered band-pass filter is required. An adaptive filtering kernel was applied for denoising, scanning through the peeled image region by region and making band-pass filtering with the peak of each local spectrum.
Waves with horizontal wavelengths of less than 12 km ere too weak to be observed by ASAI [29], so the first three layers (n = 1, 2, 3) were superposed together (Figure 5a), containing noises and a few vague ripples. Ripples with wavelengths of 25 ± 2 km, possibly from breaking GWs, were marked in Figure 5b. These ripples were perpendicular to the local GWs shown in the following sub-figures (Figure 5c-f), indicating dynamical instability [14].

Multiple Scale Analysis of GW Events
The OH airglow image (Figure 1) captured by the ASAI network contained GWs with a broad spectrum, which could be isolated separately through DWT. Considering that the assembled airglow image is 2100 pixels × 1200 pixels, the max n was set to 9. In other words, the largest scale of the basis vector was 2 9 . Figure 5 shows a series of layered images reconstructed with a basis vector at the scale of 2 n (n = 1, 2, …, 9) pixels, with one pixel representing one kilometer.
In each peeled image, the dominant wave pattern was much stronger than the noise caused by the background, clouds, or GWs of other scales, so that denoising could be performed by enhancing the strongest spectrum. For better visibility, all of the peeled images in Figure 5 were processed by the filtering kernel mentioned in Section 3.3 in order to reduce the noise and thus enhanced the wave patterns. Considering the spectrum peak of the various GWs in different regions of the peeled image, a local spectrum peak centered band-pass filter is required. An adaptive filtering kernel was applied for denoising, scanning through the peeled image region by region and making band-pass filtering with the peak of each local spectrum.
Waves with horizontal wavelengths of less than 12 km ere too weak to be observed by ASAI [29], so the first three layers (n = 1, 2, 3) were superposed together (Figure 5a), containing noises and a few vague ripples. Ripples with wavelengths of 25 ± 2 km, possibly from breaking GWs, were marked in Figure 5b. These ripples were perpendicular to the local GWs shown in the following sub-figures (Figure 5c-f), indicating dynamical instability [14].
The peeled images of the fifth to ninth layers are represented in Figure 5c-g, respectively. The GWs at different scales are illustrated in the corresponding images and their wavelengths are listed in Table 2. The medium-scale GW of the 591 ± 75 km wavelength (Figure 5g) was rarely found on ground observation [18], because its wavelength was beyond the field of view of a single ASAI and can be only captured by the ASAI network. By extracting the medium-scale GW from various GW spectrums, DWT offers a new way to measure the parameters of medium-scale GW, especially the temporal information absent in low-earth orbit satellite observations. A new assembled image ( Figure 6) with clearer wave patterns than the original one ( Figure 1) could be obtained by superposing all of the peeled and denoised images together. To derive the wave parameters, we selected the data along the green arrow in Figure 6, along which the waveforms were very clear in each peeled image. Except in layer 9 (Figure 5g), all of the measurements were made along this line in the corresponding images (shown in Table 2). The relative fluctuation Ir is defined as follows: max avg 0avg r I I I I − = (12) where Imax denotes the grayscale of the pixels at the wave peak. Iavg and I0avg are the mean value of the corresponding peeled image (shown in Figure 5) and the original image (Figure 1), respectively. As discussed in Section 3.2, the data extracted from the original image were the wave pattern only and the background was still left in the original image, so that I0avg was used as the denominator to compare the relative fluctuation between the different scales. Then, the horizontal propagation speeds were measured by tracking the movement of the wave peaks (Figure 7). The periods ere represented by the time interval between the peaks of relative fluctuation, shown in Figure 8. The peeled images of the fifth to ninth layers are represented in Figure 5c-g, respectively. The GWs at different scales are illustrated in the corresponding images and their wavelengths are listed in Table 2. The medium-scale GW of the 591 ± 75 km wavelength (Figure 5g) was rarely found on ground observation [18], because its wavelength was beyond the field of view of a single ASAI and can be only captured by the ASAI network. By extracting the medium-scale GW from various GW spectrums, DWT offers a new way to measure the parameters of medium-scale GW, especially the temporal information absent in low-earth orbit satellite observations. A new assembled image ( Figure 6) with clearer wave patterns than the original one ( Figure 1) could be obtained by superposing all of the peeled and denoised images together. To derive the wave parameters, we selected the data along the green arrow in Figure 6, along which the waveforms were very clear in each peeled image. Except in layer 9 (Figure 5g), all of the measurements were made along this line in the corresponding images (shown in Table 2). The relative fluctuation I r is defined as follows: where I max denotes the grayscale of the pixels at the wave peak. I avg and I 0avg are the mean value of the corresponding peeled image (shown in Figure 5) and the original image (Figure 1), respectively. As discussed in Section 3.2, the data extracted from the original image were the wave pattern only and the background was still left in the original image, so that I 0avg was used as the denominator to compare the relative fluctuation between the different scales. Then, the horizontal propagation speeds were measured by tracking the movement of the wave peaks (Figure 7). The periods ere represented by the time interval between the peaks of relative fluctuation, shown in Figure 8.
Atmosphere 2020, 11, x FOR PEER REVIEW 9 of 14 Figure 6. Assembled airglow image after denoising. The axes above and to the right represent the position. The lower and left axes indicate the longitude and latitude, respectively. The ASAI stations are marked with yellow points. This image is obtained by superposing all of the images in Figure 5, showing wave patterns with more visibility than Figure 1. The green arrow indicates the propagating direction of the local GWs, which is perpendicular to most wave surfaces it crosses.  Table 2, its relative fluctuation is not measured for comparison. Figure 6. Assembled airglow image after denoising. The axes above and to the right represent the position. The lower and left axes indicate the longitude and latitude, respectively. The ASAI stations are marked with yellow points. This image is obtained by superposing all of the images in Figure 5, showing wave patterns with more visibility than Figure 1. The green arrow indicates the propagating direction of the local GWs, which is perpendicular to most wave surfaces it crosses.
Atmosphere 2020, 11, x FOR PEER REVIEW 9 of 14 Figure 6. Assembled airglow image after denoising. The axes above and to the right represent the position. The lower and left axes indicate the longitude and latitude, respectively. The ASAI stations are marked with yellow points. This image is obtained by superposing all of the images in Figure 5, showing wave patterns with more visibility than Figure 1. The green arrow indicates the propagating direction of the local GWs, which is perpendicular to most wave surfaces it crosses. N/A * * The characters of the GW of 9th layer are measured along its own propagating direction in the upper left part of the image. Since the medium-scale GW of 9th layer is from a different source from other GWs listed in Table 2, its relative fluctuation is not measured for comparison.

Discussion
Based on DWT, we extracted the wave patterns of different scales form the airglow image, ranging from ripples with wavelengths of 25 ± 2 km to medium-scale GWs with wavelengths of 591 ± 75 km. To our knowledge, the medium-scale GW in Figure 5g is the first report of continuous observation on the GWs with wavelengths of hundreds of kilometers. With the temporal information from ground-based observations, we could fill the gap of satellite observations of medium-scale GWs.
The waves are displayed in the corresponding peeled images separately, so that the propagation characteristics of the GW at different scales can be measured accurately without mutual interference. Table 2 shows that both the period and horizontal speed had positive correlations with the horizontal wavelength, except those for the horizontal speed of the 7th layer, which had a decreasing value. We need more statistical analysis in order to determine whether the obvious large amplitude of the GW in the 7th layer causes this exception.
The extraction based on DWT also provides a new method of noise reduction. The quasi-monochromatic GW was dominant in the frequency domain, so that the noise and other interferences could be effectively removed from the scan of the adaptive band-pass filtering kernel. An airglow image with more visible GWs could be obtained by superposing the peeled images together after denoising. A comparison of spectra and waveforms along the measuring line before (in Figure 1) and after (in Figure 6) denoising are shown in Figure 9. All of the grayscales of the pixels measured in the original image ( Figure 1) were subtracted by their average to simply remove the mean background noise, leading to the missing peak at the zero point in the spectrum before denoising (Figure 9a). In the spectrum after denoising (Figure 9b), the noise of the high frequency and the frequency between the marked peaks were removed, besides the mean background noise. The marked peaks in Figure 9b corresponded to the frequency of the quasi-monochromatic GWs in

Discussion
Based on DWT, we extracted the wave patterns of different scales form the airglow image, ranging from ripples with wavelengths of 25 ± 2 km to medium-scale GWs with wavelengths of 591 ± 75 km. To our knowledge, the medium-scale GW in Figure 5g is the first report of continuous observation on the GWs with wavelengths of hundreds of kilometers. With the temporal information from ground-based observations, we could fill the gap of satellite observations of medium-scale GWs.
The waves are displayed in the corresponding peeled images separately, so that the propagation characteristics of the GW at different scales can be measured accurately without mutual interference. Table 2 shows that both the period and horizontal speed had positive correlations with the horizontal wavelength, except those for the horizontal speed of the 7th layer, which had a decreasing value. We need more statistical analysis in order to determine whether the obvious large amplitude of the GW in the 7th layer causes this exception.
The extraction based on DWT also provides a new method of noise reduction. The quasi-monochromatic GW was dominant in the frequency domain, so that the noise and other interferences could be effectively removed from the scan of the adaptive band-pass filtering kernel. An airglow image with more visible GWs could be obtained by superposing the peeled images together after denoising. A comparison of spectra and waveforms along the measuring line before (in Figure 1) and after (in Figure 6) denoising are shown in Figure 9. All of the grayscales of the pixels measured in the original image ( Figure 1) were subtracted by their average to simply remove the mean background noise, leading to the missing peak at the zero point in the spectrum before denoising (Figure 9a). In the spectrum after denoising (Figure 9b), the noise of the high frequency and the frequency between the marked peaks were removed, besides the mean background noise. The marked peaks in Figure 9b corresponded to the frequency of the quasi-monochromatic GWs in layers 5, 6, 7, and 8, respectively. The peaks around 20 Hz indicate the mixed waves in layer 5 (Figure 5c). Due to the 1/f noise, whose power spectral density was inversely proportional to frequency, the peak values of the low frequency were higher than that of the high frequency. Figure 9c,d represents the waveform before and after denoise, from which we could see the waveform becomes smoother after denoise. layers 5, 6, 7, and 8, respectively. The peaks around 20 Hz indicate the mixed waves in layer 5 (Figure 5c). Due to the 1/f noise, whose power spectral density was inversely proportional to frequency, the peak values of the low frequency were higher than that of the high frequency. Figure  9c,d represents the waveform before and after denoise, from which we could see the waveform becomes smoother after denoise. Figure 9. Comparison before and after denoising. (a,b) The spectra before and after denoising, respectively. The horizontal axis is the frequency, and the vertical axis is the intensity. The intensity indicates the unnormalized proportion of the wave in the corresponding frequency band. (c,d) The waveforms before and after denoising, respectively. The horizontal axis is the distance from the starting point of the measuring line and the vertical axis denotes grayscale.

Conclusions
The algorithm based on DWT worked effectively for isolating the GWs of different scales in airglow images. The original airglow image was peeled into a series of image layers, taking the db wavelets of the different scales as their basis vectors, respectively. According to the order of the basis vector from small to large, the GWs were reconstructed in the image layer with the corresponding scale. In each peeled image, the reconstructed GW was a quasi-monochromatic wave.
The algorithm isolates the quasi-monochromatic GWs of different scales according to their wavelengths, expanding the research object from the strongest quasi-monochromatic GW to GWs of multiple scales in the airglow images. In addition, an airglow image with clearer wave patterns could be obtained via reconstructing the peeled images denoised by an adaptive band-pass filter kernel. Based on the observation of the ASAI network, the extraction algorithm can represent the temporal and spatial information of GWs of various scales, ranging from small-scale GWs of tens of kilometers, to medium-scale GWs of hundreds of kilometers. The denoise algorithm makes the wave patterns more visible in the image, which will be helpful in the auto detection of GW by machine learning [30]. We will also extend this method to satellite images (visible infrared imaging radiometer suit day night band). Compared with the wave patterns, the ground artificial light, which causes a serious disturbance in the satellite image, has the characteristics of a small scale and high spatial frequency. The denoise algorithm based on DWT provides a potential way to remove the disturbances in the satellite image, and makes the GWs in the image more visible for auto detection by machine learning.
Supplementary Materials: The following are available online at www.mdpi.com/xxx/s1, Video S1: wave_4.DivX.avi. One-hour evolution of the GW isolated in layer 4 (shown in Figure 5b). Video S2: wave_5.DivX.avi. One-hour evolution of the GW isolated in layer 5 (shown in Figure 5c). Video S3: Figure 9. Comparison before and after denoising. (a,b) The spectra before and after denoising, respectively. The horizontal axis is the frequency, and the vertical axis is the intensity. The intensity indicates the unnormalized proportion of the wave in the corresponding frequency band. (c,d) The waveforms before and after denoising, respectively. The horizontal axis is the distance from the starting point of the measuring line and the vertical axis denotes grayscale.

Conclusions
The algorithm based on DWT worked effectively for isolating the GWs of different scales in airglow images. The original airglow image was peeled into a series of image layers, taking the db wavelets of the different scales as their basis vectors, respectively. According to the order of the basis vector from small to large, the GWs were reconstructed in the image layer with the corresponding scale. In each peeled image, the reconstructed GW was a quasi-monochromatic wave.
The algorithm isolates the quasi-monochromatic GWs of different scales according to their wavelengths, expanding the research object from the strongest quasi-monochromatic GW to GWs of multiple scales in the airglow images. In addition, an airglow image with clearer wave patterns could be obtained via reconstructing the peeled images denoised by an adaptive band-pass filter kernel. Based on the observation of the ASAI network, the extraction algorithm can represent the temporal and spatial information of GWs of various scales, ranging from small-scale GWs of tens of kilometers, to medium-scale GWs of hundreds of kilometers. The denoise algorithm makes the wave patterns more visible in the image, which will be helpful in the auto detection of GW by machine learning [30]. We will also extend this method to satellite images (visible infrared imaging radiometer suit day night band). Compared with the wave patterns, the ground artificial light, which causes a serious disturbance in the satellite image, has the characteristics of a small scale and high spatial frequency. The denoise algorithm based on DWT provides a potential way to remove the disturbances in the satellite image, and makes the GWs in the image more visible for auto detection by machine learning.