Detection of gravity waves in Meteosat imagery by grating cell operators

ABSTRACT An effort was made to objectively identify regions of parallel stripes in images recorded by the Spinning Enhanced Visible and InfraRed Imager instrument on board the Meteosat Second Generation geostationary satellites. Such a grating pattern can signal the presence of gravity waves which eventually may lead to atmospheric turbulence, a hazard to be dealt with in the daily routine work of aviation meteorologists. A pattern detection algorithm, based on Gabor filters plus subsequently applied grating cell operators, is implemented, adapted and tested. The method is shown to be capable of identifying grating patterns despite the unfavourable relation between the scale of the sought waves and the spatial resolution of the imagery. Some phenomena producing similar patterns (hence similar responses of the algorithm) have been identified: marine stratocumulus and regular arrangements of mountain ranges and valleys. If present in the image, they have the potential to considerably impact the analyses. The results for the 7.3 µm water vapour channel are particularly encouraging in this context as filtering the unwanted alarms is straightforward in this case.


Introduction
Gravity waves may occur in the atmosphere for various reasons and be spotted by several different sensors; see e.g. Ehard et al. (2016) for comprehensive lists. That paper, as well as several of its precursors quoted therein, also substantiates the fact that the relevance of gravity waves is by no means constrained to the troposphere, where their origin is generally found, but their vertical propagation frequently has an eventual profound impact also on the state of the middle atmosphere.
The focus here is on the appearance of gravity waves in imagery from geostationary satellites. The easiest way to find nice examples is to monitor certain privileged areas in the lee of mountain ranges where those waves can be seen quite often in the form of parallel white stripes (Figure 1). For a recent paper on physical modeling of this particular type of gravity waves, see Sachsperger, Serafin, and Grubišić (2015). For the purpose of this paper, it is enough to know that what is reflected in Figure 1 are alternating zones of upward and downward air motion resulting from the gravity waves, and that this motion is apparently vigorous enough and the moisture supply in the affected layers is sufficient to allow the formation of clouds along the lines of upward motion.
The gravity waves may become unstable, breaking down into turbulence (e.g. Sharman & Lane, 2016) which eventually occurs in layers high above the tops of the mountain range (and the clouds) and may surprise the traveller aboard an aircraft as the notorious "clear-air turbulence" (CAT) in~10 km height. Yet, spotting lee waves may also alert one to adverse flight weather conditions at the ground, as Keller et al. (2015) demonstrated with their analysis of the meteorological conditions contributing to an aircraft accident during take-off.
In principle, it should be possible to master the last-mentioned situation with a trained person that subjectively evaluates all of the available meteorological material, among that being animations of satellite imagery of the region around the airport. A stripe pattern occurring in the vicinity of the airport would not be missed (to draw the right conclusion from this observation and the other material is quite another challenge). Yet for the CAT detection issue, not only does the local situation have to be monitored but also the weather along the flight routes followed by the different aircrafts. Depending on the airport / airline that one is responsible for, the area to be scanned subjectively may be considerable. Hence, the idea emerged that an aviation meteorologist may be assisted by an objective-analysis product which highlights areas with striped patterns indicative of gravity waves. Consequently, in the "Satellite Application Facility on Support to Nowcasting and Very Short-Range Forecasting" (referred to as "Nowcasting-SAF" hereafter), a product has been included into the portfolio to provide such an objective identification. The project has a focus on European operational meteorology and so the satellite input is that of the geostationary Meteosat Second Generation (MSG) series, a succession of four identical satellites positioned around 0°longitude, with the first one launched in 2002 and the last one expected to observe the Earth well into the 2020s. Images are provided by the main instrument, Spinning Enhanced Visible and InfraRed Imager (SEVIRI), every 15 min in 12 spectral channels from 0.4 to 13.4 µm, with 3 km spatial resolution in the sub-satellite point (except for the High-Resolution Visible channel (HRVIS), with 1 km).
Screening past work of stripe detection yielded a series of interesting papers by a pair of authors (Kruizinga & Petkov, 1995, 1999Petkov & Kruizinga, 1997), where they used Gabor filters to identify single stripes, and subsequently employed grating cell operators which looked for configurations of three stripes in equal distance. These papers feature nice demonstrations showing how important the second component is, whereas solely relying on Gabor filter responses, gradients or other measures of local brightness variability does not separate stripe features from other edges present in the image very well. Though it was immediately felt that the approach could be a major contribution to successful gravity wave detection in Meteosat imagery, it was likewise expected from the very beginning that some modifications will have to be made for the application pursued here, which is well different: parts of the Kruizinga-Petkov model were explicitly incorporated to mimic the response of visual neurons (of monkeys) in case of spotting parallel lines. Although the term "grating cell operator" refers to this biological origin rather than lee waves, it is retained here for simplicity. After all, it can be argued that the goal is still to simulate with the computer what grating cells conclude if confronted with striped patterns in satellite imagery.

Algorithm
The Kruizinga-Petkov model is based on Gabor filter functions g x; y ð Þ described as follows: (2) where the parameters have the following meaning: ; η ð Þ are the coordinates of the filter centre, x; y ð Þ are those of the pixels in the surrounding. The standard deviation σ of the Gaussian factor determines the size of the field considered around ; η ð Þ (cf. parameter Ω in the following Equation (4)). The eccentricity of the filter is steered by the spatial aspect ratio γ. The parameter λ, which is the wavelength of the cosine factor, corresponds to the spatial frequency of the periodic patterns detected best. The parameter Θ 2 0; π ½ Þ specifies the orientation of the normal to the parallel stripe zones (this normal is the axis x 0 of Equation (2)). Finally φ, the phase offset in the argument of the harmonic factor, determines the symmetry of the function. However, only symmetric filters with φ=0 and φ=π, respectively, will actually be used in the following. It helps in building an efficient numerical implementation to note that with a phase shift of π, the Gabor function values differ only in their signs, i.e. g ;η;λ;Θ;φ¼π x; y ð Þ ¼ Àg ;η;λ;Θ;φ¼0 x; y ð Þ. The filter response at any location ; η ð Þ is obtained by convolving g x; y ð Þ with the input image f x; y ð Þ in a vicinity Ω around ; η ð Þ: The values of g x; y ð Þ quickly approach 0 with increasing distance from the central point, so there is no use in convolving the image with a boundless filter. Here, an implementation was used which spans a square going 3σ in the four directions AE x 0 ; AEy 0 , with Ω being the smallest square in the original, unrotated image coordinate system that encloses that square with side-length 6σ.
Though interestingly there is no trace of it in the Kruizinga-Petkov papers, other literature on Gabor filters (e.g. Boukerroui, Noble, & Brady, 2006;Lourens, Barakova, Okuno, & Tsujino, 2005) leaves no doubt about the necessity of supplementary measures to account for the fact that the integral of g ;η;λ;Θ;φ x; y ð Þ even over infinite domains does not necessarily vanish so that the result of Equation (4) is not invariant against addition of a constant value to the input image. In other words, the results presented below would be different if the input to the algorithm were in degree Celsius instead of Kelvin. In order to avoid this, the approach discussed in detail in Lourens (1998) was adopted where all negative values resulting from Equation (1) are multiplied by a factor, which is obtained by summing the positive values of Equation (1) and dividing the sum by the negative sum of the negative values. In order to preserve invariance against the sign of the image input (this symmetry is indispensable as the notions of "bright stripes" and "dark stripes" in imagery outside the visible spectrum are just a matter of convention), this modification of negative filter coefficients was done only for φ=0. The relation g ;η;λ;Θ;φ¼π x; y ð Þ ¼ Àg ;η;λ;Θ;φ¼0 x; y ð Þ remained in force and provided the Gabor function values for φ ¼ π. Figure 2 shows displays of Gabor filter functions actually used in the present study. The one from panel (a) with φ=0 yields strong positive responses for stripes of high pixel values, with matching orientation, whereas the filter of panel (b) with φ ¼ π is sensitive to stripes constituted by low values compared to the surrounding.
It was found necessary to filter r ;η;λ;Θ;φ by withholding values below a threshold r min from further processing (in other words, demanding a certain minimum response). One interesting reason for that is the architecture of the used SEVIRI instrument which at each satellite revolution acquires three image lines at once using three sensors, which theoretically should be identical but in the real world of course have unavoidable subtle differences. Over very homogeneous surfaces (e.g. extensive cloudfree ocean scenes), otherwise inexplicable widespread signals were obtained. The fluctuations are too small to be noticed in everyday's meteorological practice but that an algorithm dedicated to detect regularities may react to such perfect periodicity vigorously comes as no surprise. The threshold r min was obtained from the convolution of the Gabor filter with an amplitude-adjusted copy of itself, with the idea being that the signal from a perfectly matching pattern with very small magnitude is assessed. This threshold magnitude was set to 2% reflectance for the visible channels used below; for the water vapour channel, 0.4 K brightness temperature was chosen. Another simple measure against capturing the 3-sensor signals pattern is to avoid the angle Θ ¼ π=2, which is the most susceptible to non-meteorological stripes in the scanning direction.
After having obtained r ;η;λ;Θ;φ for all ; η ð Þ, the subsequently applied grating cell operator inspects along a straight line, direction prescribed through the angle Θ, the values orthogonal to the orientation of the filter's bars. A quantity q ;η;Θ;λ , is determined at all ; η ð Þ: where n 2 Àn max ; . . . ; n max f g(n max thus defining how many adjacent stripes one wishes to find), ρ is "a threshold parameter with a value smaller than but near one" (Kruizinga & Petkov, 1995) and the auxiliary quantities m ;η;λ;Θ and M ;η;λ;Θ;n are computed as follows: with x b c representing the floor operator (nearest integer value smaller than or equal to x) and x d e designating the ceiling operator (nearest integer value larger than or equal to x). It is irrelevant whether the pixel at the centre of the grating pattern is dark or bright. Therefore: (1) in case r ;η;λ;Θ;φ¼0 is larger than the threshold r min , the test is run with φ n ¼ 0 for even n and φ n ¼ π for the uneven n; (2) if r ;η;λ;Θ;φ¼π is larger than r min , the test is run with φ n ¼ π for the even n and φ n ¼ 0 for the uneven values.
The test returns q ;η;Θ;λ =1, signaling a grating pattern, if alternating Gabor filter signals of the same preferred orientation Θ and spatial frequency 1=λ have similar magnitudes in intervals of length λ=2 along a line passing in direction Θ, having length n max λ and being centred on point ; η ð Þ. It is of course unrealistic to expect that all the stripes in the image look the same or thateven if such a situation were encountereda parameter like λ can be guessed with such precision that the convolution values r 0 ;η 0 ;λ;Θ;φ will be identical along the line. Therefore, a certain tolerance is introduced into the above formulas, via parameter ρ. A value of 0.9 was suggested (Kruizinga & Petkov, 1999;Petkov & Kruizinga, 1997). However, that this value worked well in the quoted papers is probably owing to the application to man-made artefacts (computer-generated stripe patterns, a photograph of a textile with stripe design).
Distinctly lower values of ρ are employed hereafter as one has to be less demanding on the similarity of the stripes in the real atmosphere (i.e. their delineation in the satellite imagery).

Visible channel
Experimentation with various sets of parameter combinations suggested the following one for application of the algorithm to MSG SEVIRI 0.6 µm imagery: γ=0.4, n max ¼ 5, ρ ¼ 0:47; σ ¼ 0:4 Ã λ; the tests were run with six values of λ=2,. . .,7 image pixels, and eight values of Θ, namely π/16, 3π/16, . . ., 15π/16. Figure 3 shows a representative example result of the pattern recognition with this setting, revisiting the case of Figure 1. The stripe patterns over Spain and Algeria are signalled at many points. Occasionally, the numerical criteria are satisfied also at other places though the density of signals is generally much lower. One notable and frequent source of false alarms (with a couple of signals being visible also in Figure 3(a)) are the cloudfree Alps. However, the notion of false alarm is applicable only from the viewpoint of the gravity wave search. The pattern recognition process is absolutely correct in signaling alternating black and white patterns there, simply reflecting the regular alignment of valleys and mountain ranges, with the latter covered by snow and showing high reflectivities in the visible spectrum, and the former being snowfree and dark. Because these are valid signals, there is no point in trying to make them vanish by adjusting algorithmic parameters (thereby inevitably reducing also the other correct signals one prefers to retain). It appears more adequate to deal with such regular patterns exhibited by land surface areas in cloudfree scenes (presence of snow is no criterionsuch patterns were captured during our experiments also elsewhere) by using a proper cloud mask (the Nowcasting-SAF has one for the SEVIRI instrument in its portfolio (Derrien & Le Gléau, 2005; Derrien & Le Gléau, 2010).
There is (at least) one more meteorological phenomenon producing disturbing features of certain regularity, though they actually do not have the right length/width ratio that a human would name them a grating or stripes: marine stratocumulus arranging in mesoscale cellular structures (e.g. Koren & Feingold, 2013) often yield quite a number of unwanted alarms (Figure 3(a) displays a few signals of this type around Sicily). These large areas of quite regular arrangements of white patches are reminiscent of the checkerboard patterns for which von der Heydt, Peterhans, and Dürsteler (1992) obtained strong response for monkeys' grating cells. As the Kruizinga-Petkov algorithm was created to model such cells, one cannot blame it for responding to those patterns. On the other hand, human beings clearly sense an optical difference between the stripe patterns produced by gravity waves and the stratocumulus patches. Therefore, there is certain hope that a supplementary concept to eliminate those signals can be identified in future.

High-resolution visible channel
The HRVIS channel aboard MSG offers threefold spatial resolution (1 km pixel in sub-satellite point). This is a highly relevant improvement here since the goal is to detect parallel cloud bands on the scale of a few kilometres, so very often at the limit of what can be reasonably resolved by the ordinary SEVIRI visible channels (in Shutts (1997), figures between 9 and 22 km are given for the wavelengths of actual lee clouds; the waves obtained in the modeling efforts of Sachsperger et al. (2015) are even shorter). As the algorithm demands that dark zones and bright zones mutually resemble each other, with signs reversed where applicable, the blurring of patterns through mixing cloudfree and cloudy areas in varying and unknown proportions clearly becomes an issue when affecting too many pixels (which makes a performance as documented in Figure 3(a) even the more respectable). Figure 3(b) shows the outcome of an experiment designed to isolate the impact of higher spatial resolution using the same algorithmic parameters as in Figure 3(a) for the 0.6 µm image except for the wavelengths which have been multiplied by 3 (in terms of pixels), i.e. λ=6, 9, 12, 15, 18, 21 HRVIS pixels. Although the set of sought-after waves thus has not changed, their sharper delineation clearly makes the test more often indicate a hit. This is true for the actual gravity waves butas had to be expectedalso for the signals from the Alps and the stratocumulus around Sicily. The new signals over mainland Italy correspond to actual gravity wave patternsthere, even the human eye requires the higher resolution to identify them with certainty and so does the algorithm.
The higher resolution allows to extend the test to wavelengths λ which were not meaningful for the coarser 0.6 µm channel. Figure 4 shows the increased number of signals obtained for a run with λ=2, . . ., 21 HRVIS pixels, with the other parameters unchanged. The further improvement in the visibility of the lee waves over Central Italy, due to the inclusion of the smallest λ in the test battery, deserves particular notice. In this figure, the whole lines for which the test was passed are drawn, not only their central pixels. In Figure 3, the lines were omitted for clarity but, in fact, a positive outcome of the test indicates the presence of a grating pattern of which all pixels along the line are part. They should consequently altogether be marked in the final graphical presentation. Figure 5 contrasts the analyses of the visible channels further by showing histograms of the wavelengths of the indicated gratings. The total number of hits for the 0.6 µm visible (VIS) channel and the HRVIS (across the wavelengths λ used in Figure 3(b)) were 393 and 3565, respectively. As anticipated, the ratio becomes most remarkable when considering the smallest wavelengths: only 6 hits were recorded for the analysis of the 0.6 µm channel with λ ¼ 2 VIS pixels = 6 HRVIS pixels, whereas the HRVIS analysis yielded 496 signals for that wavelength, complemented by additional 689 hits for wavelengths λ<2 VIS pixels.

Water vapour channel
For a 24-h product, other channels have to be considered, at least at night-time where information from the visible channels is not available. A number of authors used water vapour (WV) absorption channels to monitor gravity waves, and there is one strong argument in favour of using such a channel (Feltz, Bedka, Otkin, Greenwald, & Ackerman, 2009): the results of the vertical motions induced by gravity waves may be seen even if the atmospheric conditions are not sufficiently favourable for the development of clouds which is the prerequisite for identifying the gravity waves in visible and infrared bands. Feltz et al. (2009) stated that all three investigated WV bands (6.2, 6.7 and 7.3 µm) have their value, but from the two WV channels aboard MSG, 6.2 and 7.3 µm, the latter still seems to be the more promising choice, so our experiments focussed on that channel (hereafter referred to as WV 7.3). The contrast between the bright and dark zones of the gravity waves is less pronounced in WV 7.3 than in the visible channel; it was nevertheless found that most of the parameter settings could be left unchanged, with only the tolerance parameter being relaxed somewhat to ρ ¼ 0:39. Analysing the same case as in Figures 3 and 4, the WV 7.3 evaluation in Figure 6 provides a clear picture of three agglomerations with the lines of detection showing spatially coherent directions (Spain, Algeria/Tunisia and Sardinia, the waves over the latter actually being good examples for being seen best in the WV channel). There is moreover some contamination by spurious, comparatively erratic, signals from randomly matching minor water vapour fluctuations in thick cloud decks, but the countermeasure is fortunately straightforward: as the phenomenon searched for is normally not associated with the thick clouds contributing to the largest water vapour signals (and the channel is not illumination-dependent), a constant threshold should help. To be exact, gravity waves do occur over thick convective clouds as well and may in fact result in turbulence there (Sharman & Lane, 2016). However, as there are no prospects for seeing the waves of this kind in the WV 7.3 channel, the outlined reduction of the search area may be safely done. Figure 6 shows the effect of setting r ;η;λ;Θ;φ =0 for pixels with brightness temperature below −30°C, which makes all lines passing through such pixels ineligible for the grating finder. It was verified with subjective analyses that this threshold generally leaves the waves discernible in Meteosat imagery untouched. This filter at the same time partly cures the problems with marine stratocumulus, but even without it, already the different appearance of these clouds in WV7.3 results in notably fewer false alarms around Sicily in Figure 6 than in Figures 3(b) and 4. Disturbing regular-pattern signals from the Earth's surface are at this wavelength anyhow either heavily blurred or even invisible.

Conclusions
The combination of Gabor filter plus grating cell operator is capable of detecting the stripe patterns in satellite imagery reflecting the presence of gravity waves in the atmosphere. However, as regular patterns can have other origins as well (and the grating cell operator has every right to flag them), the pattern detection result needs to be combined with auxiliary meteorological information before it can be presented to the aviation meteorologist, and some ideas were already indicated (and even partly executed, Figure 6). Figures 3, 4 and 6 show a higher density of signals q ;η =1 in areas with gravity waves than elsewhere, and this can also be exploited to filter false alarms to some degree (note that subscripts Θ and λ were dropped here since both the mentioned figures and the final map assemble the hits obtained for any tested value of Θ and λ). To distinguish the occasional isolated signal from the agglomeration of responses in an actual gravity wave area, one could use a distanceweighted spatial integration, similar to what Kruizinga and Petkov (1999) suggested:  . Results of the grating pattern search for the Meteosat-10 7.3µm water vapour image recorded on 25 January 2014, 1200 UTC. Centres of lines along which the algorithm indicated grating patterns are marked with dots: they are cyan where discarded through the imposed threshold −30°C on "too cold" pixels, the red dots are the accepted ones. Other points along the lines marked in yellow (accepted) and pink (discarded).
where Ω and σ have analogous meanings as before but are otherwise unrelated to Equations (1) and (4). The meaning of q 0 is as follows: in the original Kruizinga-Petkov formulation, q=1 only for the central point of the line where a grating pattern is detected. Presenting Figure 4, it was argued that not only the central point should be considered but all points along the line that contributed to the test's hit. However, for a formula like Equation (8), setting q ¼ 1 for all points along the line will result in an artificial bias towards waves with larger λ. Hence, the newly introduced variable q 0 is set to 1 only for the central point and the two end points of the line where a hit was recorded. Thus, every hit equally contributes three pixels of value 1 regardless of λ, yet the extension of the phenomenon is still communicated to Equation (8). When the resulting density function is projected with a logistic function onto the (0,1) range, a sort of "probability of gravity wave" field is obtained which appears to constitute a suitable presentation for the aviation meteorologist (Figure 7). Also shown in that figure are the areas where gravity waves were suspected according to a subjective analysis done beforehand. As this analysis relied heavily on the HRVIS image, several small patches are indicated where waves cannot be detected by any means from WV7.3. By and large, however, the objective analysis of this single channel compares reasonably well with the result of a time-consuming subjective scrutiny of the available meteorological material, and there is still the option to include more channels, additional angles Θ and/or wavelengths λ in order to increase the detection rate further.
It is anticipated that the next release of the Nowcasting-SAF software package will offer such a gravity wave pattern detection, in its initial form probably on the basis of the WV7.3 channel being available 24 h a day, and the false alarms being apparently easy to filter. For a 1100 × 650-pixel region covering Europe, the computation time of such a WV7.3 analysis could be brought down to approximately 10 s on a single 3-GHz processor; as the effort essentially increases linearly with the number of pixels, the figure for a corresponding HRVIS scene is around 1.5 min. Since the Gabor filter convolutions and grating searches for the diverse Θ and λ are independent from each other, there should be ample potential for reducing those numbers further by employing parallel computing.
It was demonstrated how higher spatial resolution improves the visibility of the gravity wave to the algorithm so there is every reason to look forward to the future application of the presented algorithm to imagery of Meteosat's "Third Generation" (expected to become available at the beginning of the next decade) where the imager will provide all visible imagery with 1 km resolution (sub-satellite point), and infrared and water vapour imagery with 2 km.