Satellite wind measurements represent an invaluable contribution to the description of the flow field over the oceans. Conventional cloud-tracking techniques suffer from the inability to simultaneously determine wind speed and height. Currently, the uncertainty in the independently calculated heights is the major factor limiting the accuracy of cloud motion winds. Near-simultaneous multiangle imagery from the multiangle imaging spectroradiometer (MISR) forms the basis of a unique method able to simultaneously retrieve cloud motion and height. The coupled motion and height parallaxes can be unscrambled from three properly selected multiangle views through a purely geometric, stereoscopic approach. Results based on simulated data indicate that for a mesoscale domain the average along-track and cross-track horizontal wind components may be obtained with an accuracy as good as 3–4 m s−1, and 1–2 m s−1, respectively, and with a corresponding height error of 300–400 m. The technique also possesses a limited capability to distinguish between low and high features moving at different velocities in a multilayer cloud field.
Detailed knowledge of the flow field (winds) is of paramount importance for specifying the current state and predicting the future state of the atmosphere. While the relatively dense surface observation network provides an adequate description of the large-scale horizontal wind fields over the continents, conventional measurements are extremely sparse over the vast oceanic regions. Remote sensing techniques seem to offer the only hope to fill this data gap. Space-based platforms are particularly promising with the possibility of achieving global coverage.
Winds may be derived from satellite measurements by a variety of active and passive techniques [for a comprehensive review see Isaacs et al. (1986) and Kidder and Vonder Haar (1995)]. Surface winds can be determined from the pattern of sunglint over the ocean (e.g., Wylie et al. 1981), from the observed microwave emissivities of the ocean surface (e.g., Wilheit and Chang 1980), and from microwave radar backscatter from the ocean (e.g., Cardone et al. 1983). Winds at higher atmospheric levels can be deduced from measuring the Doppler shift in emission and absorption lines in the visible and near-infrared portion of the spectrum (Burrage et al. 1996; Gault et al. 1996) or by tracking the motion of features in satellite imagery.
Undoubtedly, the tracking technique is best suited to operational purposes and enjoys the most popularity in the forecasting community. This method estimates the horizontal wind by determining the vector difference of the location of a feature in successive images. The technique may be applied in any spectral region in which distinctive features may be identified, such as the 15-μm CO2 band (Menzel et al. 1983) or the 6.7-μm water vapor band (Laurent 1993; Velden et al. 1997). In the majority of cases, however, the tracers are clouds in the visible and infrared channels, and hence the estimated winds are called cloud motion winds (CMWs). For operational purposes geostationary satellite imagery is used exclusively, although the feasibility of obtaining CMWs from overlapping polar-orbiter data has also been demonstrated for Nimbus-II High Resolution Infrared Radiometer by Shenk and Kreins (1970) and for TIROS-N Advanced Very High Resolution Radiometer by Warren and Turner (1988).
A shortcoming inherent in conventional cloud-tracking algorithms is the inability to determine the characteristic level of the winds simultaneously with the wind speed and direction. Cloud heights must be assigned by independent means, among which the most widely used (Nieman et al. 1993) are the infrared window (brightness temperature) technique, CO2 slicing, and the water vapor intercept method. These techniques heavily rely on radiative transfer calculations and require ancillary data such as cloud emissivities and an atmospheric temperature profile, the uncertainties of which may lead to height errors as large as 1–3 km. The improper height assignment appears to severely limit not only the CMWs’ impact on the quality of numerical weather forecasts but also their utility for mesoscale vorticity and divergence calculations by introducing fictitious horizontal gradients in the cloud motion data. Apart from the polar-orbiter tandem satellite concept of Lorenz (1983) and the yet experimental asynchronous stereo height and motion analysis developed by Campbell (1996) there are no techniques available that would remedy this problem by inherently combining the wind and height calculations.
This paper investigates the feasibility of a unique approach based on multiangle imagery to be obtained by the MISR instrument that enables the simultaneous retrieval of cloud height and motion. The retrieval results and error analysis presented are based on simulated data and will provide benchmarks for the evaluation of the wind product derived from actual measurements once MISR is fully operational. The technique is purely geometric in nature and uses a triplet of multiangle views to determine not only the horizontal wind vector but also the characteristic cloud height over a mesoscale domain. It is shown that certain view angle combinations result in higher accuracy cloud motion retrievals than others. Due to the near-simultaneous characteristic of the images (multiangle views are obtained within a few minutes) and to an efficient automatic image matching algorithm, cloud tracking is expected to be more reliable and much faster than that achievable with the traditional cross-correlation algorithm (Leese et al. 1971) in half-hourly geostationary satellite imagery.
The paper is structured as follows: first a brief overview of the MISR instrument and the wind retrieval problem are given in section 2. Section 3 investigates the separability of cloud height and motion with the help of a simplified geometrical model, then discusses the actual retrieval algorithm in detail including the stereo matchers used. An error analysis follows in section 4. Retrieval test results based on simulated imagery are presented in section 5. Finally, section 6 summarizes and concludes the paper.
2. MISR instrument and wind retrieval overview
The multiangle imaging spectroradiometer (MISR) was launched on 18 December 1999 aboard Terra, the first of the Earth Observing System spacecraft. The instrument is designed to measure the sunlight reflected by earth with pushbroom sensors oriented at nine different angles along track (see Fig. 1). One camera (designated An) points toward the nadir, one bank of four cameras (designated Af, Bf, Cf, and Df) points in the forward direction, and the other bank of four cameras (designated Aa, Ba, Ca, and Da) points in the aft direction. The nominal view angles, relative to the earth’s surface, are 0°, 26.1°, 45.6°, 60.0°, and 70.5° for An, Af/Aa, Bf/Ba, Cf/Ca, and Df/Da, respectively. In order to compensate for earth’s rotation, MISR cameras also contain small amounts of nominal side-looking angles (0°, ±1.0°, ±1.7°, and ±2.7° for An, A, B, C, and D cameras, respectively). Each of the nine cameras measures radiances in four spectral bands centered at 446, 558, 672, and 866 nm (blue, green, red, and near-infrared). The Terra satellite flies in a near-polar, sun-synchronous, 705-km descending orbit with a 1045 equator crossing time and has a 16-day global coverage repeat cycle. The combination of instrument geometry and orbital characteristics produces an overlap swath width of 360 km and allows that a point on the earth is observed at all nine angles within a 7-min interval. In the cross-track direction, the ground-projected instantaneous field of view (GIFOV) and sample spacing is 275 m for the off-nadir cameras, and 250 m for the nadir camera. The along-track GIFOV depends on view angle, ranging from 214 m in the nadir to 707 m at the most oblique angle (Df/Da). Nevertheless, the along-track sample spacing is still 275 m in all cameras.
The MISR wind retrieval algorithm operates on bidirectional reflectance (BRF) fields projected to the WGS84 reference ellipsoid (NIMA 1997), which approximates the earth’s surface. The method is sketched in Fig. 2. Viewing a cloud from two different directions at two different times results in an image disparity (AB2) that is a combination of a disparity due to height only (AB1) and a disparity due to motion only (B1B2). Stereophotogrammetrical height retrievals minimize the motion parallax by synchronizing the image pairs within a few seconds, while for cloud motion calculations in geostationary satellite imagery the height parallax becomes negligible because of the large distance between the satellite and the cloud. The technique presented in this paper, however, aims at unscrambling both cloud velocity and height from the measured total disparity. As shown in the next section, doing so requires a minimum of three near-simultaneous multiangle views.
In general, clouds at different altitudes may move with different horizontal and vertical velocities, and the corresponding 3D velocity ideally should be determined for each cloud feature. Such complexity would be very difficult to handle in an operational environment, and some simplifications are therefore introduced. First, the vertical cloud motion is neglected during the 7-min interval within which a cloud feature is imaged by all nine cameras so that only the horizontal velocities are determined. Second, even though the MISR wind retrieval calculates the velocity vectors of individual clouds, it is the average motion of an area of clouds that is finally reported. The individual velocity vectors representing a mesoscale region of 70.4 km × 70.4 km (256 × 256 pixels) are sorted into a 2D histogram and the two most common average velocities, corresponding to high and low features, are determined. This allows the separation of cloud from ground in the case of a broken cloud field, or of lower from higher clouds in a two-layer situation. Identifying additional layers would be difficult because the corresponding peaks are basically at the noise level due to the relative sparseness of the retrieved wind histogram.
3. Algorithm description
a. Separability of cloud motion and height
The optimal choice of wind retrieval cameras does not need the complexity of the complete wind retrieval algorithm. We therefore first describe how this choice is made using a simpler algorithm that ignores the cross-track components of the look vectors. This takes advantage of the much stronger stereo effect and correlation between cloud motion and cloud height that takes place in the along-track direction due to the large changes in viewing zenith angle from one camera to the next. In this simpler model, we can also assume a spherical earth, a circular orbit, and ignore the earth’s rotation.
Consider a cloud moving at a constant speed at a constant altitude above the surface (see Fig. 3). The cloud is viewed by two different cameras at times t1 and t2, respectively. During the interval t2 − t1 the cloud traveled the distance Sc described by
where υc is the velocity of the cloud in the along-track direction; R is earth’s radius; h is the cloud height; δ1 and δ2 are the angles between the initial radial line at time t0 and the radial lines to the surface locations x1 and x2, respectively; and γ1 and γ2 are the angles between the radial lines to the cloud and the corresponding surface locations x1 and x2. Since h ≪ R (h is typically less than 20 km) (1) can approximately be rewritten as
In (2) the camera view angles θ1 and θ2 are known, and the imaging times and surface locations are determined from stereo matching, so there are only two unknowns, υc and h. Thus two equations are needed to solve for cloud velocity and height. This means that at least three images (or two stereo image pairs) with different view zenith angles are required to separate cloud motion and height. Equation (2) can be generalized for a triplet of cameras:
This nonhomogeneous linear system can be solved only if its determinant
is nonzero. If detA equals zero, the above two linear equations are dependent on each other and cloud motion and height are inseparable. This singularity occurs for the triplets Df-An-Da, Cf-An-Ca, Bf-An-Ba, and Af-An-Aa that are symmetric to the nadir view, or for airborne measurements (AirMISR) where the flight line is straight, but does not occur, in general, for spaceborne observations on a curved orbit.
b. Camera selection for wind retrieval
Separation of cloud motion and height requires a triplet of images leading to [bu1002][xe3](93) = 84 possible camera combinations, out of which only the 80 asymmetric (nonsingular) triplets are potential candidates. Due to the varying magnitude of the corresponding determinant given in (4) some camera combinations are better suited to the retrieval than others. Combinations with the largest absolute determinant values are less prone to encounter numerical singularity and yield smaller errors in the retrieved cloud motion and height.
The determinant is plotted in Fig. 4 for different camera triplets. There are only 40 unique absolute values due to symmetry about the fore and aft directions. Considering only the value of the determinant suggests that the Df-Bf-Ca (or Da-Ba-Cf) combination, which would have the largest variety in stereo parallax and longest time difference between cameras, is best suited to wind retrieval. However, there are additional factors, such as the performance of stereo matching, that should also be taken into account. At lower sun elevation angles, either the forward cameras look at the sunlit side and the aft cameras view the shadowed side of objects or vice versa, and so matching the image of a forward camera to that of an aft camera may prove difficult. Therefore, preference is given to triplets consisting of either forward or aftward cameras only (shown as hatched bars in Fig. 4). The first two triplets meeting this criterion are Df-Cf-An (or Da-Ca-An) and Df-Bf-An (or Da-Ba-An) with respective determinant values of approximately 53 and 50. Both combinations have determinants of approximately the same magnitude. The expectation, however, is that stereo matching will perform better on triplets whose intermediate camera is farthest from both the other two cameras, that is, on the Df-Bf-An (or Da-Ba-An) triplet. Stereo is likely to work well on the D and C images, but can yield larger errors when matching the less correlated C and An images. Thus, the Df-Bf-An triplet and its aft pair the Da-Ba-An triplet are the first choices for wind retrieval. Of course, the operational camera selection algorithm has some flexibility built into it, so when the ideal set of cameras is not available the next best set will be chosen.
c. 3D ray intersection
Previously an elementary model was introduced with several simplifying assumptions, such as a nonrotating spherical earth, circular orbit, and look vectors and cloud motion having no cross-track components. Obviously, look vectors, except in the middle of the swath, have cross-track components, and the rotation of the earth, the nonsphericity of the reference surface, and the proper orbital characteristics all have to be taken into account when the look rays are calculated. Clouds also have velocity components in the cross-track direction. Here a general ray intersection algorithm is introduced that treats the 3D nature of the problem properly. The only limitation of the method discussed below is that, as before, vertical cloud motion is neglected and horizontal cloud motion is assumed to be constant over a mesoscale domain of 70.4 km × 70.4 km.
Consider a cloud moving at a constant horizontal speed, vc, at a constant altitude h above the earth (see Fig. 5). The cloud feature is imaged at times t1, t2, and t3, when its locations are Q1, Q2, and Q3, respectively. The cloud is projected to surface points P1, P2, and P3. The unit vectors of the projecting look rays are a, b, and c, respectively, and λ1, λ2, and λ3 are the corresponding scale factors of the look rays for them to intersect with the cloud. Then, the cloud motion vector, the conjugate look rays, and the surface disparities form a closed loop described by
It is preferable to represent these vector relationships in a coordinate system in which the condition of zero vertical wind can easily be imposed. Therefore, (5) and (6) are expressed in a local Cartesian coordinate system whose z axis is aligned with the zenith direction at image point P3, and having its origin at P3. The x axis and the y axis point toward the north and east, respectively. This particular coordinate system is called the local north (LN) coordinate system. The angle between the groundtrack of the satellite and the LN varies along the orbit. At lower and midlatitudes the x axis and the y axis of the LN system more or less correspond to the along-track and cross-track directions, respectively, while at higher latitudes the roles of the axes are reversed. The terms in (5) and (6) can be expressed in the LN system as P1 = (x1, y1, z1), P2 = (x2, y2, z2), P3 = (x3, y3, z3), a = (ax, ay, az), b = (bx, by, bz), c = (cx, cy, cz), and vc = (υcx, υcy, υcz). Substituting these into (5) and (6) with the condition of zero vertical wind (υcz = 0) yields the following set of equations:
where (7) and (8) are the generalizations of (3) to three dimensions. Because there are three unknowns—λ3, υcx, υcy—and four equations the problem is overdetermined, and therefore it is solved by least squares (Press et al. 1992). This arises because the small stereo effect in the cross-track direction makes the cross-track wind component almost independent of cloud height, so that the third image has little effect on this component. The height of the cloud is the z component of the look ray connecting the cloud location Q3 and the surface point P3; that is,
The above calculations provide the cloud height and the horizontal wind in the LN coordinate system assigned to P3, so υcx and υcy refer to the local north–south and east–west wind components, respectively.
d. Calculation of average mesoscale winds
Ray intersection returns the velocity vector of an individual cloud. The average cloud motion wind for a mesoscale domain of 70.4 km × 70.4 km is determined in the following manner.
First, stereo matching is performed on the 275-m resolution red band imagery of the wind retrieval camera triplet. This provides a set of matching point triplets, each corresponding to a given image feature.
For each member of a matched triplet the corresponding time tag, geolocation, and look vector are obtained from previously computed datasets, which constitute the input parameters for the 3D ray intersection algorithm.
- Step 2 is repeated for every matching point triplet. The LN system in which the wind and the cloud height are calculated will be different for every point triplet. Therefore, all wind velocities are transformed into a common LN system whose origin is in the center of the mesoscale region of interest. Then the cloud velocity vectors are sorted into a 2D histogram, where each bin of the histogram corresponds to north–south (N–S) and east–west (E–W) wind speeds. The default bin widths are 6 m s−1 in both directions. Next, the two most populated bins are determined, and the average velocities, 〈vc〉1, 〈vc〉2, and median heights, hmedian1, hmedian2, are calculated in each of the two bins. Then the average of the median heights is obtained:
e. Stereo matchers
A key component of any feature-tracking method is the stereo matcher that automatically identifies a given feature in successive images. In this study two different techniques have been tested for this purpose: the nested max (NM) and the multipoint matcher (M2) algorithms [see Diner et al. (1999) for a detailed description].
For operational MISR wind retrievals, NM is the designated stereo matcher. It is a very efficient feature-based matcher that quickly matches relatively few features with a high degree of confidence. In the course of testing NM on simulated images and also on AirMISR data, matches were typically found for 1% of all pixels within the scene with a fairly even spatial distribution. For both a target and a search window the algorithm first finds the local maxima in the signal and then the local maxima within this string of numbers and so on, and thus builds up a hierarchy of the brightest spots within the scene. Finally, each “bright spot” within the target window is matched with candidate points within the search window and the best match is determined. Here NM is favored over the cross-correlation technique because the latter would be prohibitively slow considering the very large MISR search windows that have to allow for both the wind and height disparities.
The M2 area-based matcher was used solely for research purposes. It calculates a matching metric value between a target patch and a search patch for every possible location of the search patch within a window. This metric is computed by taking all the BRF values in each patch, subtracting the mean BRF within the patch from each pixel, and normalizing by the difference in the maximum and minimum BRFs. Then, the absolute difference between these values in the target patch and the corresponding values in the comparison patch, averaged over the area of the patches and normalized by an uncertainty estimate, is tested against a threshold. The location with the smallest metric value above the threshold is then established as the final position of the search patch. The M2 is considerably slower than NM but has larger coverage (theoretically almost every pixel can be matched) and is expected to yield slightly more accurate results.
4. Error analysis
Most of the analysis below is based on the simplified model introduced in section 3a and assumes perfect navigation and stereo matching. In that simple model the look vectors have no cross-track component and thus the stereo effect is in the along-track direction only. Navigational and stereo-matching errors are discussed at the end of the section. Errors of the test runs in section 5b using the 3D ray intersection algorithm are in good agreement with these error estimates, justifying the use of a simpler approach.
a. Effect of pixel quantization
1) Along-track wind and cloud height
The effect of finite pixel size generates some error in the conjugate image location. Whenever a cloud feature is projected to the surface, quantization registers it to the center of the relevant pixel, giving rise to a maximum error of (plus or minus) half the pixel dimension.
Let us consider a single cloud first. The effect of pixel quantization on the accuracy of the velocity and height retrievals and the correlation between the velocity and height errors are shown in Fig. 6 for the default Df-Bf-An camera triplet and for the MISR pixel size of 275 m. The cloud height is 2 km and the cloud speed varies from 0 to 50 m s−1. The retrievals fluctuate around the truth depending on how close the surface projection of the cloud corresponds to a pixel center. In the worst case, the velocity and height retrieval errors are on the order of 10 m s−1 and 700–800 m, respectively. The velocity and height errors are positively correlated, that is, an overestimation of the true wind speed causes an overestimation of the true height and vice versa. Also, the height and velocity errors are related roughly linearly, a 1 m s−1 error in the along-track velocity produces an 70–80-m height error.
For a mesoscale cloud field, cloud-top heights vary within a certain range depending on the cloud type (stratiform, cumuliform, etc.). Therefore, the effect of pixel quantization will in general be different for different individual clouds within the scene, and the retrieval algorithm will produce a distribution of velocities corresponding to a distribution of pixel quantization error associated with variable cloud-top height. A best estimate of mean cloud speed is therefore determined by averaging the velocities in the most heavily populated bin of the wind speed histogram.
The results are summarized in Fig. 7a for a field of 500 clouds moving with the same speed but with cloud tops varying randomly between 1 and 4 km. The errors are generally smaller than for a single cloud (cf. Fig. 6a) and the retrievals follow a step function with the computed velocities clustering about specific values. This is due to the fairly even distribution of the winds in the most populated bin resulting in average velocities close to the bin center values. A more detailed explanation of the observed grouping of the retrieved velocities is given in section 5b.
While we assume that the cloud field advects with a certain average speed, the individual clouds may show some fluctuation about the mean velocity. Combined with cloud-top height variations this disseminates the winds within the most populated bin much more evenly, which, in turn, yields smoother average values. This is illustrated by choosing a normal distribution with a 3 m s−1 standard deviation about the mean speed, the results from which are plotted in Fig. 7b. Compared to Fig. 7a, the retrieved velocities show less scatter, the “plateaus” of the step function are more distinct.
In summary, for a single cloud the errors in the calculated cloud speed and height due to the 275-m MISR pixel size are about 10 m s−1 and 750 m, respectively. When dealing with a field of clouds with cloud-top height and/or velocity distributions, the retrieved velocities become quantized into steps of 6 m s−1, yielding a retrieval error of ±3 m s−1. A speed error of that magnitude translates into a height error of around 200–300 m (see Fig. 6c).
2) Cross-track wind
Since the cross-track stereo effect is much smaller than the along-track effect, the simplified model can also be used to estimate the expected error in this component. Assuming perfect registration, the maximum error in the measured cross-track distance traveled by a single cloud is one pixel, corresponding to a velocity error of about 2.5 m s−1. This is much smaller than the along-track error of 10 m s−1 for a single cloud. When a cloud field is considered, averaging over fluctuations will further reduce this error. The overall error in the resultant wind is thus dominated by the along-track component.
b. Errors caused by the vertical motion of clouds
One of the fundamental assumptions of the wind retrieval algorithm is of no vertical cloud motion. When this assumption breaks down the retrieved winds are biased. This wind bias is much more significant in the along-track direction, where cloud motion and height are highly correlated, than in the cross-track direction. The along-track speed bias due to vertical cloud motion can easily be estimated within the framework of the simplified model introduced in section 3a.
For a single cloud imaged by the forward Df-Bf-An camera triplet, for example, we obtain a bias in retrieved along-track cloud motion of ∓6 m s−1 (positive velocity is in the direction of the satellite orbit) for each 1 m s−1 of, respectively, upward or downward vertical cloud velocity.
The above results suggest that the along-track speed error can be considerable for rising cumulus towers. Observations of isolated cumulus bubbles show that the rate of rise of cloud tops can be as high as 8–10 m s−1 for a period of a few hundred seconds (e.g., Malkus and Scorer 1955; Warner 1977). Vertical cloud motion of such magnitude can lead to errors of a few tens of meters per second in the along-track wind. However, the average horizontal velocity of a cloud field as determined by MISR’s wind retrieval algorithm is expected to be much less affected by this type of error. Observations show that, at each instant, only about 3% of the Tropics is covered by active cumulus clouds. Furthermore, they also show that only about 1% of the area covered by active clouds is covered by undiluted convective drafts (Cotton and Anthes 1989; Black et al. 1994). Thus, it is likely that only a small portion of the calculated winds will be affected by large along-track speed errors due to vertical cloud development. Therefore, the peak of the 2D histogram of retrieved winds, though slightly smeared, will still correspond to the cloud field’s average advection velocity.
c. Stereo-matching errors
A limited stereo matcher testing was completed for both NM and M2. Test data mostly comprised of clear-sky, surface datasets but some cloudy scenes generated with a fractal cloud generator (Barker and Davies 1992) and a Monte Carlo radiative transfer code (Várnai 2000) were included as well. The results indicate that the stereo matchers used are very accurate, both the cross-track and along-track disparity error distributions were centered at zero with relatively few outliers (C. M. Moroney 2000, personal communication). The M2 proved to be slightly more accurate than NM, with the disparity errors having a sharper peak at zero. A point of concern was the increase of MISR’s footprint with view angle, which renders the oblique D and B images slightly blurrier than the An image and inevitably introduces small errors in stereo matching. These errors, however, are random and thus expected to have a negligible effect when averaged over a number of matched triplets in a domain.
In line with the above, stereo-matching errors did not seem to be a limiting factor in our wind retrieval tests. It is noted though, that both the stereo matcher and the wind retrieval tests were limited in terms of solar zenith angle and range of spatial contrast. It is understood that for extremely oblique sun angles and featureless cloud fields the stereo-matching schemes may fail or produce large errors. A comprehensive evaluation of stereo matcher performance, however, will be possible only once real measurements are available.
d. Georectification errors
The MISR product requirements call for a nadir camera geolocation accuracy of ±250 m in both the cross-track and along-track directions, and a spatial coregistration of all nine cameras with an uncertainty of ±250 m cross track and ±500 m along track. The wind retrieval algorithm is relatively insensitive to absolute image geolocation errors, due to the fact that the camera look vectors and time tags change very slowly within a mesoscale domain. Accurate coregistration of the different viewing directions, however, is crucial to obtain unbiased disparities and winds. Preliminary testing of the ellipsoid projection using simulated nominal navigation data with added measurement errors showed that the coregistration errors were generally less than 0.6 pixel in the along-track direction and 0.3 pixel in the cross-track direction (Jia Zong 2000, personal communication). This accuracy, which is significantly better than the stated product requirement and which guarantees high quality wind retrievals, is expected to carry over to actual MISR measurements.
a. Test data generation
In the absence of actual measurements, the MISR wind retrieval algorithm was tested on two simulated datasets. One of them is a clear-sky dataset provided by the Jet Propulsion Laboratory (JPL). This dataset covering a portion of central Mexico was generated from real satellite measurements by the MISRSIM simulation software (Lewicki et al. 1994). Since the data include topographical features only, the wind retrieval algorithm should return, within a certain margin of error, zero wind and the height of topography.
In order to test the retrieval of cloud-track winds, cloudy scenes were generated in the following manner. First, a cloud height field of 256 × 256 pixels (70.4 km × 70.4 km) was created with a fractal cloud generator (Barker and Davies 1992). The cloud field was assumed to consist of 256 × 256 prisms. Each one of these prisms had a base area of 275 m × 275 m and a height drawn from the previously generated fractal cloud height field. An arbitrary brightness field was then assigned to the assumed Lambertian cloud tops, with higher tops being brighter. This cloud field with all the prisms being at the same cloud-base height was assumed to move with a constant horizontal velocity. Finally, a brightness value was assigned to each pixel of a 256 × 256 pixel region on the surface ellipsoid by means of a ray-tracing algorithm (ellipsoid projection). The input to the ray tracing consisted of the geodetic latitude and longitude, and the corresponding camera look vectors and time tags of each pixel. In the course of the ray tracing, the brightness value was linearly interpolated from neighboring values, whenever the side of a cloud column was hit instead of its top.
The main limitations of the cloudy test data stem from assuming isotropic reflection. While cloud tops might be considered as Lambertian reflectors at small solar zenith angles, reflection of sunlight is definitely anisotropic at small sun elevations. Thus, forward-looking views can considerably differ from aftward-looking views in oblique sun situations. This and several other effects are not accounted for in our data simulation method. Even the above simplified cloudy test scenes, however, are suitable for studying the limitations of the wind retrieval algorithm arising from the geometry of the problem.
b. Wind retrieval tests
1) Clear-sky test run: Retrieval of zero winds
The MISR wind retrieval algorithm was first tested on a simulated clear-sky mosaic produced from several Landsat Thematic Mapper tiles by JPL. When the algorithm is applied to this dataset, it is expected to retrieve zero winds and the height of topography. The dataset consists of three blocks and covers a portion of central/southern Mexico. Each block is made up of 512 lines (140.8 km) with 2048 samples (563.2 km) in each line. The groundtrack of the satellite does not run exactly in the north–south direction but is slightly tilted eastward. For any pixel, the solar zenith angle falls within the range of (47 ± 4)°, while the solar azimuth angle is in the range of (327 ± 4)°, which corresponds approximately to the NW direction. The actual radiance data extend longitudinally and latitudinally from 99° to 103°W, and from 18° to 21°N, respectively. Areas with no valid data are depicted by white color. An expanded view of the radiance field is shown in Fig. 8. In this image the blocks are apparently unaligned and shifted with respect to each other; this, however, has no effect on the retrievals. Since wind retrieval is performed over a 70.4 km × 70.4 km (256 × 256 pixels) mesoscale domain, 30 such regions, represented by black squares in Fig. 8, were cut out of the dataset. Individual cuts are numbered 1–30 from top to bottom and from west to east.
The results of the wind retrieval test are summarized in Figs. 9 and 10. As shown in Fig. 8, cuts 1–5 and cuts 26–30 contain variable amount of invalid, spurious data. Therefore, the performance of the wind retrieval algorithm might be degraded over the top five and the bottom five regions. In fact, NM could not find a matching point triplet in cuts 1, 2, and 30, thus the wind retrieval returned no results. Here M2 yielded false retrievals for cuts 1–4. This provides 27 and 26 usable, realistic retrievals for NM and M2, respectively.
The results in Fig. 9 show that in most cases the error of the retrieval was less than 3 m s−1 in both directions and, as expected, the E–W wind component (approximately the cross-track component) is determined more accurately than the N–S component (approximately the along-track component).
The obvious grouping of the wind speeds is an artifact of the binning algorithm. Instead of having a single central bin centered at 0 m s−1 the histogram has four central bins centered at ±3 m s−1 in both the along-track and cross-track directions. The most populated bin is thus forced to be one of the four central bins yielding nonzero average velocities. Due to the relatively small number of domains only the central bins centered at respective along-track and cross-track velocities of (−3 m s−1, +3 m s−1) and (+3 m s−1, +3 m s−1) were retrieved as most populated bins. For a larger number of retrievals the distribution of most populated bins among the four central bins would have been fairly even.
It is noted that using M2 as stereo matcher yields slightly more accurate results. This is mainly because M2 gives more matching triplets than NM. The increase in accuracy is insignificant, however, and since NM is much faster than M2 it is the first choice for operational purposes.
As mentioned in section 3d, the primary reason for calculating the parameter hwind is to separate the high feature bin and the low feature bin describing features at substantially different altitudes (multilayer clouds, or cloud above ground). In this case, however, both most populated bins refer to the surface, and thus hwind corresponds to the surface as well. Experience shows that the matching point triplets found by NM or M2 are distributed fairly evenly over the domain of interest. Therefore, hwind is expected to approximate the true median height of the domain. This comparison is shown in Fig. 10 for NM and M2. The algorithm was able to capture the general outline of the height variation throughout the three-block test scene in both cases. For most cuts the difference between hwind and hmed was on the order of 200–300 m, but in a few cases it was as large as 600–700 m. For this particular set of retrievals the positive bias in the heights is caused by the mainly positive along-track velocity errors (see Fig. 9), which are positively correlated with the height errors (see Fig. 6c). For a larger number of retrievals this bias is expected to mostly disappear.
2) Cloud motion retrieval: Single cloud layer
In this section the retrieval algorithm is tested on simulated cloud fields for a wide range of cloud velocities. The cloudy test data were generated as described before with the geodetic latitude and longitude, the camera look vectors and the time tags being taken from the three-block Mexican dataset.
First a single cloud layer with a median height of 2.4 km was considered. The cloud cover was 100%, that is, the surface was not visible. The clouds were moving in the NE direction, thus the N–S and E–W velocity components were identical and varied from 0 to 50 m s−1 with 1 m s−1 increments. The retrieved cloud motion components are plotted in Fig. 11 for NM. In this and subsequent figures the solid lines represent the true cloud speeds, while the plus signs mark the retrieved values. The N–S and E–W directions approximately correspond to the along-track and cross-track directions, respectively.
The calculated along-track wind components are strongly quantized. The retrieved velocities are concentrated around 3, 9, 15, 21, 27, 33, 39, 45, and 51 m s−1. This is an artifact of the current wind-binning algorithm. As described in section 3d, the cloud velocity vectors are sorted into a 2D histogram with bin widths of 6 m s−1 in both directions. The bins have speeds falling into the range of 0–6, 6–12, 12–18, 18–24 m s−1, etc. The characteristic cloud motion velocity is determined by averaging the speeds in the most populated bin.
Obviously, the wind retrieval results plotted in Fig. 11 correspond to the bin centers of the 2D histogram. For all true wind speeds that fall into the same histogram bin, the retrieved wind speed is the center value of that bin. This yields an accuracy of ±3 m s−1 in the along-track direction. It was shown previously that for a single cloud the error in the along-track component due to pixel quantization can be as high as 10 m s−1 (see Fig. 6). This means that in the case of a distribution of cloud-top heights, the measured along-track components may scatter within a wide range around the true value, which in turn yields that the speeds in the most populated bin are distributed fairly evenly and vary from the lower to the upper limit of that bin. Therefore, the average of these speeds is approximately the center value of the bin. This effect is much less pronounced in the cross-track direction, since the wind error for a single cloud (1–2 m s−1) and hence the spread of the measured winds around the truth is much smaller. Not surprisingly, the retrieved cross-track components follow the true values very closely with a typical error of 1–2 m s−1.
The root-mean-square errors (rmse’s) corresponding to the retrievals using NM and M2 are rmseNM = 2.2 m s−1 and rmseM2 = 1.8 m s−1, respectively. Compared to NM, results obtained with M2 show less fluctuation and the step function character of the retrievals is more distinct. The general accuracy of the calculated winds, however, has not increased dramatically by using M2 instead of NM.
The height retrievals corresponding to the winds plotted in Fig. 11 are shown in Fig. 16a. The computed median heights are scattered within ±300 m about the true median height of the cloud field (solid line). Here M2 yielded cloud heights of similar accuracy.
3) Cloud motion retrieval: Broken cloud field above ground
In the previous section, the retrieval algorithm was tested on a single cloud layer with individual clouds moving at exactly the same speed. There are, however, several realistic situations when different features move at different velocities, for example, a multilayer cloud field with each layer having its own advection velocity. Identifying all the occurring cloud velocities operationally would be too ambitious a task, therefore the MISR wind retrieval was designed to calculate only the two most common velocities. In order to test this feature of the retrieval, the algorithm was run on simulated scenes containing partially cloudy land. Simulated broken cloud fields were combined with surface reflectance data from the three-block Mexican dataset. As before, clouds were moving at the same speed in the NE direction. As described in section 3d, the algorithm determines the average wind speeds of the two most populated bins of the wind vector histogram and then labels them as low or high depending on their corresponding median heights. The expectation was that one of the two retrieved bins would refer to surface features (zero wind), while the other would correspond to the advecting clouds.
The computed winds for two different cloud fraction values are shown in Figs. 12–15, for both NM and M2. In these figures, the x axis is the true cloud velocity. The solid lines represent the truth, with the horizontal line at 0 m s−1 retrieved wind referring to the surface and the diagonal line referring to the true cloud velocities. The plus signs (Vx,low, Vy,low) and the circles (Vx,high, Vy,high), respectively, denote the low feature bin and the high feature bin. The corresponding median heights are plotted in Fig. 16 for NM (M2 yielded similar heights). Here the retrieved low and high values are marked by boxes and filled boxes, respectively, while the solid lines refer to the true median cloud and surface heights.
Figures 12, 13, and Fig. 16c plot the retrieved speeds and heights, respectively, for a cloud cover of 20%. The median cloud-top and surface height was 2.9 and 1.1 km, respectively. With NM as stereo matcher (Fig. 12), up to a true wind speed of approximately 30 m s−1, the low feature bin corresponds to the surface in most cases, while the high feature bin refers to the clouds. At speeds close to 0 m s−1 the peaks corresponding to the ground and clouds are basically merged together causing a large scatter in the calculated heights at those speeds. At speeds >6 m s−1 the peak referring to the clouds separates from the steady peak referring to surface features, which is also reflected in the retrieved height values. At cloud speeds >30 m s−1 the surface signal was not picked up and so even the low feature bin corresponds to the clouds usually with large errors. The explanation for this stems from the fact that at larger speeds the clouds move considerably from one image to another. Hence the common surface area seen by all three cameras becomes very small. For a low-coverage feature matcher such as NM, the surface point triplets thus become increasingly difficult to identify as the wind velocity increases.
Applying M2, as opposed to NM, as stereo matcher results in sharper histogram peaks (Fig. 13). This probably is due to the larger number of measured wind vectors, which makes the separation of the modes easier. There is much less fluctuation in the retrieved wind speeds and the steps in the along-track component are clearer. Notice that M2, being a high-coverage area matcher, was able to pick up the surface signal even at large cloud speeds. This was also evident in the retrieved heights, which are not shown here for M2. In both of the above cases the error in the along-track and cross-track components was, as before, 3–4 m s−1 and 1–2 m s−1, respectively, while the median heights were retrieved with an uncertainty of about ±300 m.
When the cloud fraction is sufficiently large, no parts of the surface are seen by all three cameras. This is depicted in Figs. 14, 15, and Fig. 16b, where the cloud cover was 60% and the median cloud-top height was 2.7 km. Both peaks correspond to the clouds, with the high feature bin usually overestimating and the low feature bin underestimating the true velocities. Due to the correlation between the speed and height errors (see Fig. 6c), this yields a respective positive/negative bias in the high/low height values.
The above demonstrated that it is possible to determine the two largest peaks in the wind vector histogram that refer to features at different altitudes moving at different velocities. The successful separation of the peaks requires them to be sufficiently apart from each other. The M2 always provided slightly better results than NM, underlying the importance of having as many matching point triplets to work with as possible. In the case of an extensive high cloud field, the lower altitude features, such as the ground terrain or a low cloud layer, might be obscured (especially in the most oblique Df/Da images) and thus only the high altitude features are likely to be identified.
6. Summary and conclusions
The foregoing analysis and use of simulated data demonstrates the feasibility of retrieving cloud-tracked winds with the multiangle imaging spectroradiometer. The unique multiangle views provided by the instrument allowed for the development of a purely geometric enhanced stereo technique that, unlike traditional methods, computes not only the horizontal cloud motion vectors but their characteristic altitudes as well. Cloud motion and height are shown to be separable from the total parallaxes on an orbital scale with a curved orbit. The simultaneous retrieval of cloud motion and height requires at least three images. Not all the possible camera triplets, however, are equally suitable for retrieval. Certain combinations are less sensitive to the inevitable errors in the measured disparities than others. Taking into account the performance of stereo matching as well, the Df-Bf-An and Da-Ba-An camera triplets appear to be best suited to cloud motion retrieval.
A general 3D ray intersection algorithm was introduced to calculate cloud motion and height simultaneously. Over a mesoscale domain of 70.4 km × 70.4 km, the algorithm determines the average velocities of low- and high-altitude features and assigns them to characteristic heights. The main limitations of the method stem from neglecting vertical cloud motion and assuming a constant horizontal cloud advection over the domain. Under some atmospheric conditions involving intense convection, frontal wind shear, etc., these assumptions may fail and the retrieved winds for these cases will likely be unreliable.
Error analysis reveals that the uncertainty in the retrieved velocity due to pixel quantization is much larger in the along-track direction, where retrieved cloud motion and height are highly correlated, than in the cross-track direction. For a single cloud, the errors in the calculated along-track and cross-track speeds are on the order of 10 m s−1, and 3 m s−1, respectively. This translates into a height error of 700–800 m. For a mesoscale cloud field, due to fluctuations in cloud height and speed the errors in the average along-track and cross-track velocities and in cloud height should be reduced to 3–4 m s−1, 1–2 m s−1, and 300–400 m, respectively.
The wind retrieval algorithm was tested on simulated clear-sky and cloudy datasets. For stereo-matching purposes, two different algorithms were assessed: (i) an extremely fast one with sparse coverage (NM), and (ii) a more time-consuming one with relatively dense coverage (M2). For the clear-sky dataset provided by the Jet Propulsion Laboratory, the wind algorithm retrieved the zero wind within the expected rmse of 3–4 m s−1. Cloudy scenes were generated with a ray-tracing algorithm assuming Lambertian reflection. A wide range of cloud velocities was successfully retrieved with the expected errors; the rmse was 2–3 m s−1. Median heights were retrieved with an accuracy of ±300 m both for clear-sky and for cloudy scenes. The algorithm’s capability of determining the two most common velocities was demonstrated by running the code on scenes containing broken cloud fields as well as surface features. For a modest cloud amount the peaks corresponding to the terrain and the clouds could clearly be distinguished. For a larger cloud cover, however, only the high-level features could be tracked, since low-altitude features were hidden from view. In general, the magnitudes of retrieval errors were as expected. Here NM and M2 produced results of comparable accuracy in most cases, except from the “multilayer” test, when M2 yielded less noisy, much more unambiguous retrievals than NM, due to a larger number of tracked features.
It is emphasized here that the above results should be interpreted with care, since they are based on a limited set of tests. All the test experiments assumed perfect image navigation. These assumptions eliminated certain sources of error that will inevitably occur in real data. Therefore, the obtained retrieval accuracy should be considered as a theoretical upper limit.
Future work should aim at improving the quality of the calculated winds. The observed quantization of the retrieved winds, especially in the case of the along-track component, is mainly an artifact of the wind-binning algorithm. The accuracy of the retrieved along-track components could be enhanced by simply reducing the width of the histogram bins. A smaller bin width, however, means a greater number of less populated bins, which may render the determination of the histogram peak more difficult. So there is really a trade-off here between a small bin size having small theoretical error and the ability to find the mode of the histogram. Once real data are available, different binning strategies should be tested that may lead to less quantization and hence smaller overall errors.
This work, which is based on the M.S. thesis of Á. Horváth, was supported by the Jet Propulsion Laboratory (JPL) of the California Institute of Technology under Contract 960489. Our thanks to Robert Ando of JPL for providing some of the test data, and to C. M. Moroney and T. Várnai for many helpful discussions. It is also a pleasure to acknowledge the significant contributions of D. J. Diner, J. Zong, and J. P. Muller to the Algorithm Theoretical Basis Documentation, which prompted this study.
Corresponding author address: Ákos Horváth, Department of Atmospheric Sciences, The University of Arizona, Tucson, AZ 85721-0081.