Study of aircraft icing warning algorithm based on millimeter wave radar

In order to avoid accidents due to aircraft icing, an algorithm for identifying supercooled water was studied. Specifically, a threshold method based on millimeter wave radar, lidar, and radiosonde was used to retrieve the coverage area of supercooled water and a fuzzy logic algorithm was used to classify the observed meteorological targets. The macrophysical characteristics of supercooled water could be accurately identified by combing the threshold method with the fuzzy logic algorithm. In order to acquire microphysical characteristics of supercooled water in a mixed phase, the unimodal spectral distribution of supercooled water was extracted from a bimodal or trimodal spectral distribution of a mixed phase cloud, which was then used to retrieve the effective radius and liquid water content of supercooled water by using an empirical formula. These retrieved macro- and micro-physical characteristics of supercooled water can be used to guide aircrafts during takeoff, flight, and landing to avoid dangerous areas.


Introduction
Clouds consist of water droplets or ice crystals suspending in the air (Liu et al., 2014(Liu et al., , 2015Wang et al., 2016a, b, c). Non-freezing water droplets can be called supercooled water when the ambient temperature is below the freezing point. Such water is very unstable and can easily freeze into ice due to slight vibration. If an aircraft with surface temperature slightly below 0°C collides with supercooled water, ice can form on the aircraft surface, which is a phenomenon known as aircraft icing. The ambient temperature, liquid water content, and diameter of supercooled water droplets are important factors in aircraft icing (Ahrens, 1982;Eagleman, 1985).
Many aircraft accidents have occurred due to supercooled water. During the mid-1980s, Federal Aviation Administration reported 12 serious accidents due to aircraft icing in which turboprop engines stopped suddenly under the influence of supercooled water. In four of the accidents, the aircraft went out of control and fell to the ground. On November 21, 2004, a CRJ200 airplane of China Eastern Airlines crashed near Baotou airport when there was less than a minute after the takeoff; experts speculated that the cause of the crash was aircraft icing. On June 3, 2006, a military transport plane of Chinese Air Force crashed somewhere in Anhui; after investiga-tion, experts established the direct cause of the accident as the frequent trips of the aircraft across frozen regions. In order to ensure flight safety of aircraft in icing meteorological conditions, protective equipment should be used. The existing aircraft icing protection systems can be divided into two categories (Wu, 2012): anti-icing systems and de-icing systems. Anti-icing systems mainly include liquid anti-icing technology, electric anti-icing technology, hot gas anti-icing technology, and evaporation anti-icing technology. De-icing systems mainly include pneumatic de-icing technology and electrical pulse de-icing technology. The aircraft icing phenomenon has been extensively studied in developed countries such as United States, France, Germany, and Canada (Addy Jr. et al., 1997;Gent et al., 2000;Bragg et al., 2005).
If icing intensity is weak and the coverage area of supercooled water is small, aircrafts can fly safely using anti-icing and deicing equipment; if icing intensity is strong and the coverage area of supercooled water is large, aircrafts should stay away from the icing area by changing flight altitude or flight path. Therefore, the analysis and forecasting of supercooled water (Serke et al., 2014;Yang et al., 2014;Alexandrov et al., 2016), which causes aircraft icing, is a high-value task in aviation meteorology.
In this paper, we propose an algorithm for determining the macro-and micro-physical parameters of supercooled water. The retrieved results can be used to judge whether the aircraft can fly through a supercooled water area. These warning parameters can be retrieved by using a combination of a threshold method, a fuzzy logic algorithm, and power spectral density based on the data of millimeter wave radar, lidar, and radiosonde. The retrieved results mainly include the coverage area, effective droplet radius, and liquid water content of supercooled water. The results can be used to avoid the accident risk of aircrafts due to icing.

Coverage area of supercooled water identified by the threshold method
A cloud with a lidar backscatter coefficient greater than 5 × 10 -5 m sr -1 can be regarded as a supercooled water layer (Luke et al., 2010). Khain et al. (2001) thought that the radar reflectivity factor of a cloud that includes only droplets should be less than -20 dBZ. Based on extensive experiments with 7 yr of multi-sensor observations, Shupe (2007) found that the spectral width of millimeter wave radar exceeds 0.4 m s -1 for a cloud consisting of su-percooled water and ice particles, while this value is below 0.4 m s -1 for a cloud consisting of only supercooled water or ice particles. Note that although lidar can identify supercooled water layers of a cloud, the lidar signal is attenuated strongly; thus, lidar can only identify the bottom of a supercooled water layer. In order to better identify the coverage area of supercooled water, millimeter wave radar and radiosonde data should be used. An algorithm for using this data is shown below: 1) The temperature of supercooled water ranging from -40 to 0°C can be determined by using sounding data; 2) In the temperature region ranging from -40 to 0°C, when the lidar backscatter coefficient of a cloud is greater than 5 × 10 -5 m sr -1 , the cloud can be regarded as a supercooled water layer; 3) In the temperature region ranging from -40 to 0°C, the region of a supercooled water layer that cannot be detected by lidar due to attenuation could be identified by the spectral width and linear depolarization ratio (LDR) of millimeter wave radar. Specifically, the spectral width should exceed 0.4 m s -1 and the value of LDR should be below -15 dBZ; 4) The entire area identified by lidar and millimeter wave radar can be labelled as supercooled water area.

Cloud phase type classified by a fuzzy logic algorithm
Compared with neural network approach (Inoue, 1987;Lee et al., 1990;Shi and Qu, 2002), a fuzzy logic algorithm mainly includes four processes (Liu and Chandrasekar, 2000): 1) fuzzification, 2) inference of rule, 3) aggregation, and 4) defuzzification. A fuzzy logic algorithm for identifying the hydrometeor type is based on converting the matrix of radar data to the matrix of particle phase type using conversion rules described by the membership function (Al-Sakka et al., 2013). In this paper, four measurements, namely, radar reflectivity factor, Doppler velocity, spectral width, and temperature, were used as input variables of the fuzzy logic algorithm. Trapezoidal membership functions were used to fuzzify the four measurements. The different membership functions were constructed by using the threshold of Shupe (2007). The hydrometeor types and output variables of the fuzzy classifier system are listed in Table 1. Table 1 shows that the output of the fuzzy classifier system is one of the hydrometeor types: (1) snow, (2) ice, (3) mixed, (4) liquid, (5) drizzle, and (6) rain. Figure 1 shows the block scheme of the fuzzy classifier. The input parameters of the fuzzy classifier are radar reflectivity factor, Doppler velocity, spectral width, and temperature, and the output parameter is the hydrometeor type.

DECEMBER 2017
Next, we describe the four steps in the classification of hydrometeors: 1) Fuzzification The process of fuzzification involves conversion of the four input parameters into fuzzy sets. This conversion process mainly relies on membership function. For each input parameter (radar reflectivity factor, Doppler velocity, spectral width, and temperature), there are six membership functions corresponding to the six hydrometeor types. The membership functions are defined as MBF ij (fuzzy sets); subscripts i and j stand for the input parameter and the particle type, respectively. There are many forms of membership functions, and an asymmetric trapezoidal function was chosen as the basic form of membership functions used in this paper. The trapezoidal shape is determined by coefficients X 1 , X 2 , X 3 , and X 4 , and the membership function is expressed as follows: where T stands for membership, and x is the considered parameter (radar reflectivity factor, Doppler velocity, spectral width, or temperature). Determination of the coefficients of the trapezoidal function (X 1 , X 2 , X 3 , and X 4 ) is the key to the classification result of the fuzzy logic algorithm. These coefficients would not greatly influence the retrieved results due to ambiguity of fuzzy logic algorithm, that is, the tiny change of these parameters is not enough to affect the retrieved results. The threshold values for the different particle types, in terms of the radar reflectivity factor (Z h ), Doppler velocity (v D ), spectral width (w D ), and temperature (T), are given in Table 2 (Peng et al., 2011).
2) Inference of rule By constructing the membership functions successfully, the fuzzy parameters can be obtained by using the  fuzzy logic algorithm; rules are used to describe the complex relationships between input and output fuzzy variables. The process of deducing the "strength" of these fuzzy parameters from the strength of the antecedents is called rule inference. The most commonly used inference methods are correlation product and correlation stack, which are shown below: where P ij stands for the membership, subscripts i and j are the parameter and the particle type, and A i is a weight coefficient.
The results from the two inference methods have no significant differences. In this study, we selected the method of correlation stack and set A i = 1, which is shown in Eq. (3).
3) Aggregation This step involves aggregating S j computed by the inference method and setting the largest S j as Q, that is, Q = S jmax .

4) Defuzzification
The final results obtained by the fuzzy logic algorithm are the index values shown in Table 1, from which the particle classification can be identified.

Microphysical parameters of supercooled water retrieved by Doppler power spectral density
In order to acquire the microphysical parameters of different hydrometeor particles forming a hybrid cloud, the noise level in Doppler power spectra density must be computed first. The actual power spectral densities of different particles can be obtained after eliminating the noise power spectral density, using the following basic principle: Assuming that particles inside the radar effective irradiation volume have the same attributes, the power spectral density of these particles in still air can be regarded as a single pulse function. However, actual cloud particles detected by radar have different attributes and the air has a vertical velocity; thus, Doppler power spectral densities detected by radar are results of convolution of signal power spectral densities and the air turbulence power spectral density (Gossard, 1994;Gossard et al., 1997). Because air turbulence can be regarded as a Gaussian signal, the result of convolution is still a Gaussian function. The power spectral density of individual particles can be expressed by a Gaussian distribution and those of a variety of particles can be expressed by the addition of these Gaussian distributions. For a 35-GHz millimeter wave radar, there is velocity ambiguity and spectral overlap in the case of large raindrops, and the power spectral density is not Gaussian. In order to avoid such complications and ignore the influence of updraft and downdraft on particle movement, we analyze only stratiform clouds. In general, the bimodal spectrum of a stratiform cloud can be regarded as a spectrum with one mode due to liquid droplets and the other due to ice particles (Kollias et al., 2011).
The Gaussian distribution form of individual particles can be written as follows: where S 0 is maximum signal power, v is the Doppler velocity, v 0 is the corresponding velocity, and σ is the spectral width. Because the output power spectral density contains noise, we should remove the noise first to retrieve the radar reflectivity factor and other physical parameters using power spectral density.
Assuming that the noise of the radar system is N(v), the measured output of power spectral density can be expressed by the addition of meteorological signal power and noise power: There are several methods to remove the noise from the measured power spectral density. Battan determined a fixed noise threshold for separating signal from noise (Battan, 1964). Donaldson (1967) and Sekhon and Srivastava (1971) utilized a noise threshold below the peak to separate the signal of meteorological targets from radar system noise. However, these methods of removing noise are problematic because noise is not fixed and its value changes randomly. Hildebrand and Sekhon (1974) treated noise as Gaussian white noise and proposed an algorithm to separate noise from the detected signal. Currently, this method has been widely used and has the following steps: 1) The Doppler power spectral density is smoothed in three points.
2) The smoothed Doppler power spectral density values are arranged in the order of smallest to largest; the new sequence of Doppler power spectral density values can be defined as A i (i = 1, 2, 3, ..., n).
3) The average P n and variance Q n of the new sequence A i is computed by using following formulas: 4) A new variable R n is defined: When R n > 1, the corresponding P n value can be regarded as noise power.

Case study
The data of millimeter wave radar (35-GHz Copernicus), lidar (CT75K), and radiosonde (http://weather. uwyo.edu/upperair/sounding.html) taken during 1600-2400 UTC 11 January 2011 at Chilbolton Observatory (United Kingdom) were used to retrieve the coverage area of supercooled water. Table 3 shows the specifications of 35-GHz Copernicus and lidar, which can be found as follows.
Figures 2a-d show the radar reflectivity factor, Doppler velocity, spectral width, and linear depolarization ratio of the 35-GHz millimeter wave radar, respectively. Figure 3 shows the backscatter coefficient of lidar, and Fig. 4 shows the sounding temperature data at 2400 UTC.
According to the algorithm for identification of supercooled water, which was introduced in Section 2.1, the retrieved coverage area of supercooled water based on data of Figs. 2-4 is shown in Fig. 5. Figure 5 shows that supercooled water is present at the height of 1 km during 1600-1700 UTC, the area of the supercooled water layer is gradually increasing during 2000-2400 UTC, and the maximum coverage of supercooled water occurs during 2200-2400 UTC, with a coverage area from 1-to 5-km heights.
The fuzzy logic method introduced in Section 2.2 was Backscattering sensitivity~10 -7 m -1 sr -1 at low altitude  Fig. 6. Figure 6 shows that the bottom of the cloud is almost entirely liquid water, with the height of liquid water about 1.3 km during 1600-1700 UTC. After 2200 UTC, the coverage area of liquid water ranges in height from 1 to 5 km, and other areas of the cloud consist almost entirely of ice particles. This conclusion is in agreement with the retrieved result shown in Fig. 5. Moreover, the cloud after 2000 UTC is a mixed-phase cloud containing ice particles and liquid water droplets.
In order to obtain microphysical parameters of supercooled water in a mixed-phase cloud, the Doppler power spectral density of supercooled water should be extracted from that of the mixed-phase cloud. Figure 7 shows the power spectral densities detected by 35-GHz millimeter wave radar.
The power spectral densities at 2000, 2100, and 2200 UTC have obvious bimodal spectral distribution in the case of a coded mode, which implies that there are two different hydrometeors (ice particle and liquid water). These power spectral densities at different heights are shown in Fig. 8.
Note that the Doppler power spectral densities in Fig. 8 contain noise power. We eliminated the Gaussian white noise as described in Section 2.3 to obtain the final spectral power containing the information of only meteorological particles. Figure 9 shows the original and denoised power spectral densities at different heights and times. Figure 9 shows that the bimodal nature of the power spectral density becomes clearer after denoising; the left part of the bimodal power spectral distribution stands for ice particles and the right part stands for supercooled water. The velocity of ice particles is negative, indicating falling velocity (toward the radar), and the velocity of supercooled water is positive, indicating rising velocity (away from the radar). There are two modes of pulse width in the 35-GHz Copernicus millimeter wave radar, uncoded mode and coded mode. The measured velocity of the uncoded mode ranges from -5.3633 to 5.3214 m s -1 and that of the coded mode ranges from -2.6816 to 2.6397 m s -1 .
Assuming that the particle size distribution is normal, it can be expressed as    where N is particle number concentration, r is the particle radius, r 0 is the median radius, x = lnr, x 0 = lnr 0 , and σ is the spectral width of log spectra.
When radar reflectivity factors of a mixed phase, ice phase, and liquid phase are computed based on their Doppler power spectral densities, the effective radius (r e ) and liquid water content (LWC) are given by the following expressions (Shupe et al., 2004): The effective radius and LWC computed by using results shown in Figs. 9b, 9d, and 9f are listed in Tables 4  and 5, respectively. Table 4 shows that the retrieved effective radii of the mixed phase are smaller than those of the ice phase, with a minimum relative error of 165.34%, and larger than those of the liquid phase, with a minimum relative error of 57.69%. Table 5 shows that the retrieved LWC values of the mixed phase are larger than those of the ice phase or the liquid phase. The minimum relative error is 94.67% for the ice phase and 91.77% for the liquid phase. Based on above analysis, the microphysical parameters retrieved for the mixed phase are completely different from those of the ice phase or the liquid phase.
Based on above analysis, the macro-and micro-physical parameters of supercooled water, such as coverage area, effective radius, and LWC, can be retrieved by millimeter wave radar, lidar, and radiosonde measurements. The classification results can provide reference information about whether an aircraft can fly through a supercooled water layer and the flight paths to avoid supercooled water.

Summary
The macro-and micro-physical parameters of a supercooled water layer were extracted from millimeter wave    radar, lidar, and radiosonde data. A threshold method was used to retrieve the coverage area of supercooled water and a fuzzy logic algorithm was used to classify the observed hydrometeor particles. The macrophysical characteristics of supercooled water can be confirmed by comparing the results of the threshold method and the  fuzzy logic algorithm. To obtain the microphysical characteristics of supercooled water, a method of extracting the unimodal power spectral density of supercooled water from the power spectral density of a mixed phase after denoising was proposed. This method can be used to retrieve the effective radius and LWC of supercooled water. A case study verified that the coverage of supercooled water at different heights and times can be identified by the proposed method and that the microphysical parameters of the mixed water-ice phase are completely different from those of the constituent phases.