Derivation of Cloud-Free-Region Atmospheric Motion Vectors from FY-2E Thermal Infrared Imagery

Since the launch of the first meteorological experimental satellite, TIROS-1, in 1960, satellites have been widely used in global meteorological observation. Data and inversion products, such as atmospheric motion vectors (AMVs), derived by tracing clouds or water vapor, have long been studied and verified. Cloud motion vectors (CMVs) were successfully derived as early as in the late 1960s, by detecting the motion of cloud in visible images using traditional film animation technology (Izawa and Fujita, 1968). Later, (Endlich et al., 1971) proposed a method to calculate CMVs using a pattern recognition technique, and (Leese et al., 1971) derived CMVs using the cross-correlation coefficient of brightness temperature in infrared images, which enabled the automatic computation of CMVs. Water vapor motion vectors are derived from geostationary satellite water vapor imagery with a similar procedure (Stewart et al., 1985; Velden et al., 1997). Countries in Europe, along with America and China, have set up regulations in the operational calculation and quality control of AMVs, and these quality-controlled AMV data can improve the initial conditions of numerical weather prediction (NWP) (Schmetz et al., 1993; Nieman et al., 1997; Xu and Zhang, 2006). Velden et al. (1992, 1998) and Le Marshall et al. (1994, 2002) indicate that a positive effect can be obtained if quality-checked CMVs are assimilated into NWP for numerical hurricane track forecasting. (Hou and Xu, 2006) revealed the relationship between upper-tropospheric circulation and rain bands in China through the use of observed precipitation data and GMS-5 CMVs from 1998 to 2002.

AMVs derived from geostationary satellite imagery have evolved into an important data source of meteorological information (Menzel, 2001). However, AMVs for the lower troposphere in cloud-lacking areas (cAMVs) are in need of further development. CMVs inversed from visible and IR imagery can only be obtained within cloud regions, and water vapor motion vectors inversed from water vapor imagery are mainly available for the upper troposphere, because the weighting functions of water vapor channels around 6.7 μm are typically within the 250-600 hPa (Velden et al., 1997).

Progress has been made in this field in the last few years, with a focus on water vapor and aerosols in cloud-free areas (Guo et al., 2010; Zhang et al., 2013). For instance, the "Split Window Difference (WD)" method has been adopted to suppress the influence of surface temperature fluctuation and enhance the weak, spatial texture signals of water vapor and aerosols (Zhan et al., 2012; Wang et al., 2014a). The result for the clear region surrounding a typhoon, reported by (Zhang et al., 2015), showed the RMSE of WD-derived cAMVs with respect to the NCEP 800-hPa wind field to be between 5.4 and 7.0 m s^-1. Meanwhile, (Yang et al., 2014) proposed the "Temporal Difference (TD)" method, which calculates the "temporal difference" of brightness temperature from two adjoining IR1 (or IR2) images for the same purpose as WD. (Wu et al., 2014) obtained cAMVs with the TD method and assimilated the quality-checked cAMVs into NWP, revealing that cAMVs can improve the initial field of NWP and then decrease the error in the prediction of the cyclone track. Most recently, through a combination of the WD and TD methods, we (Wang et al., 2014b; Zhao et al., 2015) proposed the "Second Order Difference (SD)" algorithm, with which the surface temperature fluctuation interferences can be suppressed by twice as much, and the weak signals of water vapor and aerosols can be strengthened.

This paper, as an extension to the preliminary study of (Wang et al., 2014b), introduces the difference principles and algorithm schemes (i.e., SD, WD and TD) for extracting AMVs in cloud-free areas from FY-2E IR imagery (section 2), and then compares their performances in a set of controlled experiments (section 3) and case studies (section 4). Conclusions are given in section 5.

The spectral radiance in the thermal infrared channel is commonly represented by \begin{equation} \label{eq1} I_\lambda=\varepsilon_\lambda B_\lambda (T_{\rm s})\tau_\lambda(0)+\int_0^\infty{B_\lambda(T(z))}\left[\frac{\partial\tau_\lambda(z)}{\partial z}\right]dz ,(1) \end{equation} where the subscript Λ is used to specify a certain channel, such as "1" for IR channel around 10.3-11.5 μm and "2" for IR channel 11.6-12.8 μm, ε_Λ is the surface emissivity, T _s is surface temperature (units: K), B_{Λ(T s)} is the surface blackbody radiance, B_Λ(T(z)) represents the atmospheric blackbody radiance at the height z, and τ_Λ(0) represents the transmittance of the whole atmospheric layer——often abbreviated as τ_Λ. The expression ∂τ_{Λ(z)/∂ z} is defined as a weighting function, i.e., \begin{equation} \label{eq2} W_\lambda(z)=\frac{\partial\tau_\lambda(z)}{\partial z} ,(2) \end{equation} where τ_Λ(z) is the transmittance from the height z to the satellite, and is represented by \begin{equation} \label{eq3} \tau_\lambda(z)=\exp\left\{-\int_z^\infty{k_\lambda(z')\rho_\alpha(z')dz}'\right\} , (3)\end{equation} where k_Λ is the mass coefficient of absorption at wavelength Λ, and ρ_α is the mass of absorbing water vapor and aerosols. Aerosols are considered together with water vapor in this study because aerosols increase the opacity in the split window channels effectively, and can be taken as a tracer in addition to water vapor for extracting AMVs in cloud-free areas from FY-2E IR imagery.

Since the thermal radiance exhibits good linearity with the temperature in the vicinity of the thermal infrared split window (Prabhakara et al., 1974), Eq. (1)——in terms of brightness temperature (TB)——becomes \begin{equation} \label{eq4} {\rm TB}_\lambda=\varepsilon_\lambda \tau_\lambda T_{\rm s}+\int_0^\infty{T(z)}W_\lambda(z)dz .(4) \end{equation}

Based on the above theory, (Zhang et al., 2013) simulated the brightness temperatures of FY-2E split window channels with MODTRAN 4 (Berk et al., 1999) for the three reference atmospheric profiles for the tropics (15^°N), midlatitude summer (45^°N, July) and midlatitude winter (45^°N, January). In the computation, surface temperature (T _s) ranged from 298 to 302 K for summer, and 271 to 275 K for midlatitude winter; water vapor content (WV) was between 0.5 and 5.5 g cm^-2; and aerosol optical depth (AOD) was between 0.1 and 2.0, the aerosol type was set to "Rural (VIS = 5 km)", and the surface emissivity was assumed to be 0.98. To convert AOD to VIS for computations with MODTRAN, the relationship between AOD and VIS (horizontal meteorological visibility) was used as follows (He et al., 2003): \begin{equation} \label{eq5} {\rm VIS}=\frac{(1-b){\rm AOD}}{a{\rm AOD}} .(5) \end{equation} Here, the coefficients a and b are equal to 0.1202 and 0.2974 for spring and summer, and 0.1419 and 0.1377 for autumn and winter, respectively.

One of the approximations to Eq. (5) for split window channels is \begin{equation} \label{eq6} {\rm TB}_\lambda=a_{0,\lambda}+a_{T_{\rm s},\lambda} T_{\rm s}+a_{{\rm WV},\lambda}{\rm WV}+a_{{\rm AOD},\lambda}{\rm AOD} ,(6) \end{equation} where a_0,Λ, a_{T s,Λ}, a _WV,Λ and a _AOD,Λ are regression coefficient and can be obtained with regression analysis on the simulated brightness temperature database, as given in Table 1. Since a_{T s,Λ}, a _WV,Λ and a _AOD,Λ express the sensitivity of TB_Λ to T _s, WV and AOD, the thresholds for ST, WV and AOD increment observations with FY-2E thermal infrared channels, with a sensitivity threshold (NEdT) of 0.2 K, would be as those given by the last three columns in Table 1. It is quite common, except in dry winters, for each of the thresholds to occur within 30 minutes in a cloud-free area (Zhan et al., 2012; Yang et al., 2014), and therefore FY-2E thermal infrared channel imagery can be used for WV and AOD texture tracking, as long as the influence from the surface feature variation can be alleviated.

Surface interference, which is a main error for AMV retrieval, is caused by surface temperature variation and emissivity homogeneity. The optimum surface in this sense is ocean surface, as compared with land. To alleviate the influence from the surface feature variation, especially for land surface, we adopt the difference principle; that is, ∆TB (temporal difference of TB between time t₂ and t₁), TBD (split window difference of TB between channels 1 and 2), and ∆TBD (the second-order difference of TB). They are defined as follows: \begin{eqnarray} \label{eq7} \Delta {\rm TB}_\lambda&=&{\rm TB}_{\lambda,t_2}-{\rm TB}_{\lambda,t_1} ,(7)\\[1mm] \label{eq8} {\rm TBD}_t&=&{\rm TB}_{1,t}-{\rm TB}_{2,t} ,(8)\\[1mm] \label{eq9} \Delta {\rm TBD}&=&{\rm TBD}_{t_2}-{\rm TBD}_{t_{\rm 1}}=\Delta{\rm TB}_{\rm 1}-\Delta{\rm TB}_{2} . (9)\end{eqnarray} Assuming that surface emissivity, ε, changes little in a short period of time (e.g., 30 min) in a cloud-free area, and the difference of ε between channels 1 and 2 is negligible, after substituting Eq. (5) into the above three equations, one has \begin{eqnarray} \label{eq10} \Delta {\rm TB}_\lambda&=&\varepsilon(\tau_\lambda T_{\rm s})|_{t_1}^{t_2}+\left.\left[\int_0^\infty {T(z)} W_\lambda(z)dz\right]\right|_{t_1}^{t_2} ,(10)\\[1mm] \label{eq11} {\rm TBD}&=&\varepsilon(\tau_1-\tau_2)T_{\rm s}+\int_0^\infty{T(z)}[W_1(z)-W_2(z)]dz ,(11)\\[1mm] \label{eq12} \Delta {\rm TBD}&=&\varepsilon[(\tau_1\!-\!\tau_2)T_{\rm s}]|_{t_1}^{t_2}\!+\!\left.\left\{\int_0^\infty{T(z)}[W_1(z)-W_2(z)]dz\right\}\right|_{t_1}^{t_2} ,(12)\nonumber\\ \end{eqnarray} where ∆τ=τ₁-τ₂ represents the transmittance difference between the two split window channels. It can be seen that TB, ∆TB, TBD and ∆TBD each consist of two parts: the first part depends on the surface characteristics (surface emissivity, surface temperature) multiplied by the transmittance of the whole atmosphere; while the second part is only determined by the status of the atmosphere (temperature and absorption coefficient of the absorbing molecules and aerosols). In order to extract the signal of atmospheric motion in a cloud-free area, we hope that the spatial texture in imagery as AMV tracers depends only on the status of the atmosphere, rather than on the surface characteristics. Results from calculation show that it is common in nature that either τ_Λ|t1t2, ∆τ or ∆τ|_t1^t₂ is less than τ_Λ. This indicates that, in the sense of alleviating the surface influence, the differential imagery, such as ∆TB, TBD and ∆TBD, should be better than the original images for tracing the atmospheric textures in a cloud-free area for atmospheric motion calculation, and the ∆TBD imagery is expected to be the best. The methods associated with ∆TB, TBD and ∆TBD have been respectively named as the TD method (Yang et al., 2014), the WD method (Zhan et al., 2012), and the SD method (Wang et al., 2014b).

Figure 1.  Algorithm schemes for the SD, WD and TD methods: (a) scheme for WD; (b) scheme for WD with quality control; (c) scheme for TD with IR₂; (d) scheme for TD with quality control; (e) scheme 1 for SD; (f) scheme 2 for SD. In the figure, IR_Λ,t indicates an image taken by an IR channel at wavelength Λ at time t.

The algorithm schemes for the above methods are shown in Fig. 1, in which a circle with "-" means computing the difference between two images, a circle with "R" means computing the maximum correlation between two images, and a circle with "QC" means quality control in terms of spatial and temporal consistency between two AMV fields. One can see that the WD method works with four images observed with split window channels at two different times, and two more images observed at the third time are required for QC. The TD method can use either of the two IR channels, but IR₂ is better than IR₁, based on a comprehensive consideration according to Table 1, especially for tracking water vapor textures. Therefore, the TD method works with only three images observed sequentially at three different times. However, the method does need all six images as well if QC is required. Two schemes can be designed for SD according to the definition of the second-order differentiation. They lead to the same result, but Scheme 1——as shown in Fig. 1e——is more time-saving than Scheme 2 (Fig. 1f). Either of the two schemes needs at least six images to work, and so we can see that two more images observed at the fourth time are required for QC.

In order to be conducive to the height assignment of AMV retrievals and better understanding the surface interference situation, a sensitivity analysis of FY-2E weighting height using the specific spectral response function has been performed according to the theory, as in Eqs. (10)-(12). As shown in our previous study on height assignment and the sensitivity of spectral radiance to the distribution of water vapor for different types of situation (Yang et al., 2014), the peak level of the water vapor weighting function based on Eq. (10) is at approximately 850 hPa under the tropical and midlatitude summer atmosphere, and at 800 hPa under the U.S. standard atmosphere. The results from Eqs. (11)-(12) have similar features. This implies that the surface interference has been alleviated and the height assigned to the derived AMVs in cloud-free regions is approximately 850 hPa under the tropical and midlatitude summer atmosphere, and 800 hPa under the U.S. standard atmosphere; that is, low-level winds.

A simplified flow chart for the simulation experiment is shown in Fig. 2. The supposed AMV field follows the 850 hPa wind field and is taken as representative of the true AMVs. The surface temperature fluctuation in 30 min is sequentially set to 0, 0.5, 1.0 and 1.5 K (Hong et al., 2014) to simulate the random variation of surface temperature.

Figure 2.  Distribution of meteorological factors for the simulation experiment over tropical ocean: (a) surface temperature (units: K); (b) water vapor content (units: g cm^-2); (c) AOD; (d) supposed 850-hPa AMVs (average speed was set at 10 m s^-1). The area shown by Fig. 3d is the same as the interested region shown by the black square inside Fig. 3a, b and c.

The simulated area was a cloud-free area of 400 km × 400 km over the tropical ocean in a northwesterly wind controlled region. The ST, atmospheric WV, AOD, and supposed 850 hPa AMVs for the simulation experiment, are shown in Fig. 3. The SD, WD and TD imagery simulated with the standard tropical atmospheric profile with MODTRAN 4 is shown in Fig. 4. Though the simulated SD, WD and TD imagery were influenced by many elements, such as the vertical distribution pattern of water vapor and aerosols, and the algorithm for processing the convergence and divergence caused by the wind field, the difference caused by the uncertainty in these elements was negligible compared with the major texture features in SD, WD and TD imagery. By tracking the motion of textures with the maximum correlation method (Leese et al., 1971), the AMVs derived with the SD, WD and TD methods are shown in Fig. 5, for four values of ST fluctuation (0, 0.5, 1.0 and 1.5 K).

Figure 3.  Flow chart for the simulation experiment.

Error statistics, such as bias (BIAS), mean vector difference (MVD), standard deviation (STD) and RMSE, were calculated based on the following definitions (Nieman et al., 1997): \begin{eqnarray} \label{eq13} {\rm VD}_i&=&[(U_i-U_{ri})^2+(V_i-V_{ri})^2]^{1/2} ,(13)\\[1mm] \label{eq14} {\rm BIAS}&=&\frac{1}{N}\sum_{i=1}^N[(U_i^2+V_i^2)^{1/2}-(U_{ri}^2+V_{ri}^2)^{1/2}] ,(14)\\[1mm] \label{eq15} {\rm MVD}&=&\frac{1}{N}\sum_{i=1}^N{\rm VD}_i ,(15)\\[1mm] \label{eq16} {\rm STD}&=&\left[\frac{1}{N}\sum_{i=1}^N({\rm VD}_i-{\rm MVD})^2\right]^{1/2} ,(16)\\[1mm] \label{eq17} {\rm RMSE}&=&({\rm MVD}^2+{\rm STD}^2)^{1/2} , (17)\end{eqnarray} where (U_i,V_i) represents the truth of AMV u- and v-components, as shown in Fig. 3d; and (U_ri,V_ri) represents their retrievals, as in each of the 12 panels in Fig. 5. Table 2 presents the results. According to Table 2, as well as the difference between Fig. 3d and each of the 12 panels in Fig. 5, it can be seen that, although further quality control (Xu et al., 2008) can be applied to the retrievals, all three methods worked well. SD was generally the best of the three methods, and the quality of the retrieved AMVs decreased as ST fluctuation increased.

Figure 4.  The SD, WD and TD imagery simulated with the standard tropical atmospheric profile with MODTRAN4.

Figure 5.  AMVs derived with the SD (left-hand panels), WD (middle panels) and TD (right-hand panels) methods, with ST fluctuation values of (a-c) 0, (d-f) 0.5, (g-i) 1.0 and (j-l) 1.5 K.

For deriving AMVs in a cloud-free area with satellite observations in practice, there are generally four steps that are followed (Zhang et al., 2015): cloud masking (Ma et al., 2007); TD, WD or SD image production; a tracing procedure with maximum cross-correlation (Wang and Zeng, 1996); and quality control (Yang et al., 2014; Xu et al., 2008). The height assignment is simply completed with 800 hPa, as discussed in section 2.

Figure 6.  Case 0000 UTC 4 Aug 2012: (a) FY-2E IR1 channel image with a dash-lined square inside to show the almost cloud-free conditions over China's eastern ocean region; (b) NCEP wind field (800 hPa); (c) NSMC CMVs; (d) TD-derived AMVs overlain on the TD image; (e) WD-derived AMVs overlain on the WD image; (f) SD-derived AMVs overlain on the SD image. The dash-lined square in each of panels indicates the same geographic area for comparison. The cloud-free areas are white in panel (a) and grey in panels (d), (e) and (f) as seen in the dash-lined squares.

Figure 7.  As in Fig. 6 but for the 0600 UTC 29 July 2011 case.

Figure 8.  As in Fig. 6 but for the 0600 UTC 17 March 2013 case.

Results from studying many cases were similar, three of which we present here as examples:

The first case was 0000 UTC 4 August 2012 (Fig. 6). Figure 6a shows the FY-2E IR1 imagery (pixel resolution: 5 km × 5 km) and the dash-lined square inside shows almost cloud-free conditions over China's eastern ocean region. Figure 6b is the NCEP 800-hPa wind field, Fig. 6c is the CMVs produced by National Satellite Meteorological Center of China (NSMC), and Figs. 6d-f show the AMVs obtained with TD, WD and SD, respectively. It can be seen that the cloud-free area is between typhoon Haikui and Damray, lying to the south of the subtropical high, and the wind field inside is in general closely consistent with the NCEP data (Fig. 6b). Moreover, SD performed better than any of the other methods in representing the true situation.

Figure 7 shows the results for a case on 29 July 2011. On the west of a cold front associated with a low-pressure system is a cloud-lacking, dry-tongue area (highlighted by the dash-lined square), where the air at low levels flowed counterclockwise (Fig. 7b). The airflow patterns shown in Figs. 7d-f are generally consistent with that shown in Fig. 7b, and the SD-derived field looks better than any of the others in terms of continuity and consistency.

Lastly, Fig. 8 is an example of tracing sand aerosol for AMVs in a cloud-free area in North China. The IR1 image observed by FY-2E at 0600 UTC 17 March 2013 shows clearly (Fig. 8a) that the area in the dash-lined square was almost free of cloud for cloud tracking, but full of northwesterly wind (Fig. 8b). Therefore, cloud motion vectors were rare (Fig. 8c), but the TD, WD and SD methods were able to retrieve a number of northwesterly wind vectors (Figs. 8d-f) that were consistent with the NCEP data.

The above three cases indicate that the different methods for AMV inversion can provide additional low-level airflow information in cloud-free areas, where the operational CMV method fails. Table 3 presents the error statistics of AMV retrievals with respect to NCEP FNL data. The accuracy is acceptable, as compared with the accuracy of the cloud motion vectors. In general, the SD method is superior to the TD and WD methods in terms of consistency and continuity.

According to the results of the present study, the weak signals representing the spatial texture of water vapor and aerosols in cloud-free areas shown in FY-2E split window images can be strengthened and traced for AMVs in the lower troposphere around 800 hPa.

Surface temperature interference when identifying the spatial texture of water vapor and aerosol must be suppressed, and so algorithms for retrieving AMVs [e.g., TD, WD and SD, based on ∆TB (temporal difference of brightness temperature imagery between time t₂ and t₁), TBD (split window difference of brightness temperature imagery between channels 1 and 2), and ∆TBD (the second-order difference of brightness temperature imagery with respect to both time and wavelength)] have been proposed. The simulation analysis in the present study demonstrated that all of the algorithms, though not perfect in alleviating the influence of surface temperature fluctuation in some cases, generally work well in retrieving AMVs, and the second-order difference algorithm seems to be the best (as compared to either the split window difference or temporal difference methods).

The cases studies reported in this paper illustrate that AMVs can be obtained in certain areas where the traditional cloud-tracking method fails. The accuracy with respect to NCEP 800 hPa reanalysis data was found to be acceptable, as compared with the accuracy of the cloud motion vectors.

The algorithms proposed can take either aerosol or water vapor, or both, as a whole tracer for AMV retrieval, but would fail in cases where both aerosol and water vapor are not available and suitable for tracing. Further study on the performance of the difference algorithms in an operational system, as well as verifications of the AMVs produced operationally for cloud-free regions and the evaluation of their impact on NWP models, is also needed.