- Email: [email protected]
Systematic tracking of boundary layer aerosols with laser radar
Optics & Laser Technology 34 (2002) 177 – 185
www.elsevier.com/locate/optlastec
Systematic tracking of boundary layer aerosols with laser radar N.C. Parikha; ∗ , J.A. Parikhb a Department
of Physics and Earth Sciences, Central Connecticut State University, 1615 Stanley Street, New Britain, CT 06050, USA of Computer Science, Southern Connecticut State University, 501 Crescent Street, New Haven, CT 06515, USA
b Department
Received 4 April 2001; received in revised form 27 August 2001; accepted 25 October 2001
Abstract Atmospheric aerosols exhibit a high degree of variability in their properties and their spatial and temporal distribution. Mapping atmospheric aerosols is of great importance for monitoring the e3ects of anthropogenic activities and natural processes upon local air quality, as well as for producing accurate input for radiative transfer models. To help enhance our understanding of radiative, physical, chemical, and dynamic processes in the atmosphere, laser radar is used to provide systematic monitoring of the temporal evolution of the aerosol rich boundary layer. The Micro Pulse Lidar (Light Detection and Ranging) System is capable of long-term autonomous mapping of vertical aerosol structure. We present a new automated method for tracking boundary layer height from Micro Pulse Lidar data and show the results of this method for a continuous 16 h study conducted on the east coast of the United States. ? 2002 Published by Elsevier Science Ltd. Keywords: Lidar; Aerosols; Boundary layer
1. Introduction An understanding of the spatial and temporal distribution of aerosols is important for atmospheric studies including air quality assessment, pollution management, global climate change research, and precipitation and atmospheric dynamics. The in
0030-3992/02/$ - see front matter ? 2002 Published by Elsevier Science Ltd. PII: S 0 0 3 0 - 3 9 9 2 ( 0 1 ) 0 0 1 0 7 - 4
a Geiger mode avalanche photodiode photon counting detector [1]. The system is capable of continuous, autonomous data acquisition in both daytime and nighttime. Several MPL systems are currently in operation across the country. Each MPL system can generate large volumes of continuous data. For long-term atmospheric monitoring, it is often desirable to employ lidar data in conjunction with data from other instruments [2]. Therefore, to facilitate data analysis and intercomparison of results between instruments, an automated method of tracking boundary layer heights from MPL images is especially important. We have developed a methodology for tracking boundary layer heights with Micro Pulse Lidar data. Various deInitions of average boundary layer height have been used in the literature for derivation of the boundary layer height from lidar backscatter signals. Flamant et al. [3] deIned average boundary layer height as the base of the transition zone (the interface between the mixed layer and the free troposphere)—a deInition associated with a Irst derivative approach. Menut et al. [4] examined methods to calculate average boundary layer height deIned as the middle of the transition zone using second-order derivatives of backscatter signals. Boundary layer height over a continuous time span is calculated in this study as a continuous track based on the gradient of
178
N.C. Parikh, J.A. Parikh / Optics & Laser Technology 34 (2002) 177–185
Fig. 1. Boundary layer track determination
the backscatter signal subject to continuity and curvature constraints. The lidar measurements consist of photoelectrons per microsecond backscattered from a range r. Images formed as a function of time map the signal strength to color intensity at a given range and time. The vertical extent of the aerosol-rich boundary layer is often estimated visually from such images. These manual estimations can be time-intensive and images may be diMcult to interpret, especially in regions of low visual contrast. In addition, interpretation by di3erent individuals may result in subjective di3erences in estimated boundary layer heights. As illustrated in Fig. 1,we applied edge detection, followed by coarse boundary detection to these images. The result provided input for an active contour model of the boundary layer height. The active contour model minimizes an energy functional consisting of terms describing contour curvature, continuity and consistency with image edge points [5].
2. Experimental procedure As a single frequency lidar unit, the Micro Pulse Lidar transmits pulses of 523 nm light at 2500 Hz and measures backscattered photons as a function of time. The signal intensity is a3ected by both extinction and backscatter. The high-resolution instrument we used in this study yielded data with 30 m resolution. The data are corrected for the 1=r 2 range dependence and for instrument calibration factors including detector deadtime, detector afterpulse, and Ield of view overlap as described in Parikh et al. [6]. Each data point consists of a 60 s average of the signal returned from each range. An image of the temporal evolution of the atmosphere is generated from continuous data by plotting the corrected returned signal on a grid where the y-axis represents vertical range above the instrument and the x-axis represents the time of the observation. The signal strength can be represented by a false-color or gray-scale mapping.
Images obtained from a sequence of MPL measurements we conducted at a site o3 the east coast of the United States (near Atlantic City, NJ) over an unbroken 16 h period from 4 pm to 8 am the following morning are shown in Figs. 2a and b. This dataset will be used to illustrate our automated method for boundary layer height determination from MPL data. Sky conditions during the observation period were clear from 4 pm until approximately 10 pm, at which time scattered high-altitude clouds could be seen over the observation site. Shortly after midnight, a low cloud deck settled over the observation site and remained for the duration of the observation period. Visual observations of sky conditions during the period are consistent with sky condition reports from the ASOS Atlantic City, New Jersey Weather Station located approximately six miles southeast of the observation site. Surface wind speed and temperature over the observation period from the ASOS station (shown in Figs. 3a and b) were moderate to mild. The images shown in Figs. 2a and b depict byte-scaled corrected signals mapped into the integer interval between zero and 255. Strong signal returns from overlying dense clouds were mapped into a value of 255 and extremely weak signal returns were mapped into zero, thereby enhancing mid-range intensity values indicative of aerosol variations. 3. Edge detection Canny edge detection techniques [7] using MATLAB were applied to the byte-scaled, mid-range enhanced images described in the previous section to produce a binary edge image. In the Canny edge detection procedure we applied a horizontal Gaussian Ilter to the input images [8]. Non-maximal suppression then removed all locations which were not along a ridge of the gradient. The resulting image was then thresholded to produce a binary edge image in which all image points were either valued at zero (no edge was detected) or 1 (an edge was detected). Parameters that must be speciIed for the edge detection technique include the standard deviation of the Gaussian Ilter and the thresholds used to produce the binary edge image. For our experiments, we used a standard deviation of 1 and a threshold of 0.0125 for weak edges and a threshold of 0.0313 for strong edges. Weak edges were included in the Inal binary edge images if they were connected to strong edges. Figs. 4a and b show the binary edge images for each day overlaid on the original MPL dataset images. The edge detection operation picks out edges corresponding to locations of strong changes in intensity. 4. Coarse boundary detection Once edges have been identiIed it is necessary to discriminate between those edges which potentially correspond to locations representing areas of division between the aerosol
N.C. Parikh, J.A. Parikh / Optics & Laser Technology 34 (2002) 177–185
179
Fig. 2. MPL data (a) day 338, 4:00 –11:59 pm and (b) day 339, 12:00 –8:00 am.
rich boundary layer and the upper atmosphere, and all other edges. In this phase of the procedure, we constructed a coarse boundary consisting of a linked series of blocks that initialized and constrained the Inal boundary produced by the active contour model. The construction of the coarse block boundary required analysis and processing of the original byte-scaled images as well as of the binary edge images. The Irst step in the coarse boundary detection phase of the procedure was to identify blocks within the binary edge image with strong edge persistence over both the altitude
range and time frame represented by the length and width of the block respectively. For each center point of each block, a heuristic measure h was deIned to characterize block edge persistence. For our experiments, h was deIned as h = k(columnedge )
(1)
where k may be any positive constant and columnedge is the total number of columns containing one or more edge points. The next step, altitude block localization, involved locating, for each time frame, the highest altitude block
180
N.C. Parikh, J.A. Parikh / Optics & Laser Technology 34 (2002) 177–185
Fig. 3. Surface observations: (a) Wind speed and (b) Temperature.
with maximal edge persistence. The binary edge image was divided into a sequence of Ixed time frames corresponding to block width. For each time frame the values of the heuristic measure at the center point of each altitude-overlapping block within the time frame were compared to determine the uppermost block associated with the highest heuristic value for that frame. The co-located region in the original byte-scaled image was then checked for cloud amount. The block selected for the given time frame was rejected if more than half of the pixels in the co-located region were identiIed as clouds and the procedure continued until an appropriate block was found. Cloudy pixels were identiIed in the original image as pixels whose value exceeded a threshold value Tc . The rejection of a block in this step of the procedure did not necessarily exclude the block from inclusion in the Inal coarse boundary due to linking criteria applied in the subsequent step. The discontinuous block boundary resulting from altitude block localization served as an initialization for an iterative linking technique which yielded the Inal coarse boundary approximation. Adjustments to the current coarse boundary approximation at locations of discontinuity between boundary blocks were based on the values of the heuristic measure associated with the blocks in the given time frames as well as on continuity with neighboring blocks. The linking procedure terminated when all blocks in the coarse boundary approximation touched any neighboring blocks to the right or left at one or more points. The result of the analysis was a linked series of blocks which corresponded to a coarse estimate of upper altitude constraints on the track of the top of the aerosol rich boundary layer. Figs. 5a and b show the results of the coarse boundary detection technique as a block boundary layer track overlaid over the original byte-scaled data. These results were obtained using a 16 × 16 block size
corresponding to approximately 480 m in altitude and 16 min in time. 5. Active contour The coarse block boundary is reIned using an active contour model. The active contour model requires good initialization. While the midlines of each block in the block contour could be used for starting points, these points may not be edge points. Therefore we applied line integral minimization of an abbreviated energy functional, Einit , to the coarse block midline initialization to snap the contour to a set of edge points where possible, thereby providing an improved initialization. Here Einit is given by Einit = Eimg ds; where
Eimg (q) = 2(es − eq ) − 1
for q = s;
Eimg (s) = −2es − 1;
(2)
where s is the current contour point, q is a potential contour point at the same time and e is the value of the point in the binary edge image [9]. The result is a reIned initial contour estimate shown in Figs. 6a and b. The Inal step to obtain a smooth contour through edge and non-edge points to automatically track the boundary layer height is to minimize the full energy functional, E, constrained to altitude variations. E is given by E = (Econt + Ecurv + Eimg ) ds; where
Econt (q) = (yq − ysprevious )2 Ecurv (q) = (ysprevious − 2yq + ysnext )2 :
(3)
N.C. Parikh, J.A. Parikh / Optics & Laser Technology 34 (2002) 177–185
181
Fig. 4. Binary edge results overlaid on original image (a) day 338, 4:00 –11:59 pm and (b) day 339, 12:00 – 8:00 am.
Here, yq is the altitude of a potential contour point, q, for a given time, t. ysprevious is the altitude of the contour point at time t − 1. ysnext is the altitude of the contour point at time t + 1, and , , and are weights. Finite di3erences are used to approximate Irst and second derivative properties [5]. The energy functional minimizes the weighted sum of a continuity term (an estimate of the Irst derivative property), a curvature term (second derivative property), and an edge term. Note that regions containing cloud pixels may be included in the Inal boundary layer track, if doing so results in the overall minimization of the energy functional.
The Inal result, shown in Figs. 7a and b, is a contour which shows the derived boundary layer height tracked over the time span of the image.
6. Discussion of results The edge detection step, the Irst major component of the technique, was parameterized to yield an abundance of edges including several subtle edges that appeared in the Inal results. For this application the Canny edge detector
182
N.C. Parikh, J.A. Parikh / Optics & Laser Technology 34 (2002) 177–185
Fig. 5. Coarse block contour (a) day 338, 4:00–11:59 pm and (b) day 339, 12:00–8:00 am.
was found to be superior to other standard edge detection operators such as the Sobel operator. Strong edges that were selected by the Canny edge detector and rejected in subsequent steps included edges from cloud shadows, edges from clouds not associated with the top of the boundary layer, and edges demarcating internal structure within the aerosol layer. The result of the second major step of the technique was a connected sequence of Ixed-size regions that formed a coarse outline providing an upper constraint on the Inal boundary layer track. In this coarse, boundary
detection stage of the procedure physical constraints on the data, problem-speciIc expertise, and data-speciIc meteorological information can be incorporated. In addition, adaptation in this stage allows this technique to be extended to analysis of data from other lidar instruments. The Inal boundary layer track was produced using a customized active contour model. The iterative energy functional minimization procedure was terminated when no more than 10% of the contour points were changed in an iteration. Weights on edge, continuity, and curvature terms may be set depending on the degree of desired sensitivity
N.C. Parikh, J.A. Parikh / Optics & Laser Technology 34 (2002) 177–185
183
Fig. 6. Active contour initialization: (a) day 338, 4:00–11:59 pm and (b) day 339, 12:00–8:00 am.
of the contour to speciIc physical features such as aerosol plumes. 7. Conclusions An automated method for tracking the temporal variation of boundary layer height from Micro Pulse Lidar data has been developed. The method has been illustrated on a dataset that spans both day and nighttime conditions and clear and cloudy conditions. Complex internal structure within the
aerosol layer marked by strong edge intensities was also prevalent. The results match well with visual determinations of boundary layer height. Advantages of this automated method are objectivity, consistency, and
184
N.C. Parikh, J.A. Parikh / Optics & Laser Technology 34 (2002) 177–185
Fig. 7. Boundary layer track: (a) day 338, 4:00–11:59 pm and (b) day 339, 12:00–8:00 am.
well suited for adaptation to use with data from other lidar instruments.
The Irst author (NCP) would also like to thank her students, Jennifer Falciglia, Gabriel Ma3ei, and Stephen Price for their assistance in collection of the laser radar data used in this work.
Acknowledgments This work was supported by Connecticut State University research and release time grants. The authors wish to thank Frank Schmidlin at NASA Wallops Flight Facility for the use of the Micro Pulse Lidar system and Steven Palm at NASA Goddard Space Flight Center for helpful discussions.
References [1] Spinhirne JD. Micro pulse lidar. IEEE Trans Geosci Remote Sensing 1993;31(1):48–54. [2] Siebert P, Beyrich F, Gryning S, Jo3re S, Rasmussen A, Tercier P. Review and intercomparison of operational methods for the
N.C. Parikh, J.A. Parikh / Optics & Laser Technology 34 (2002) 177–185 determination of the mixing height. Atmos Environ 2000;34: 1001–27. [3] Flamant C, Pelon J, Flamant P, Durand P. Lidar determination of the entrainment zone thickness at the top of the unstable marine atmospheric boundary layer. Boundary-Layer Meteorol 1997;83: 247–84. [4] Menut L, Flamant C, Pelon J, Flamant P. Urban boundary-layer height determination from lidar measurements over the Paris area. Appl Opt 2000;38(6):945–54. [5] Williams DJ, Shah M. A fast algorithm for active contours and curvature estimation. Comput Vision Graphics Image Process 1992;55(1):14–26.
185
[6] Parikh NC, Schmidlin FJ, Campbell JR, Hlavka DL. Experimental characterization of laser radar instrumentation. J Indian Geophys Union 2000;4(1):91–5. [7] Canny J. A computational approach to edge detection. IEEE Pattern Anal Mach Intell 1986;8(6):679–98. [8] Jain R, Kasturi R, Schunck B. Machine vision. New York: McGraw-Hill, 1995. [9] Niedermeier A, Romaneessen E, Lehner S. Detection of coastlines in SAR images using wavelet methods. IEEE Trans Geosci Remote Sensing 2000;38(5):2270–81.
www.elsevier.com/locate/optlastec
Systematic tracking of boundary layer aerosols with laser radar N.C. Parikha; ∗ , J.A. Parikhb a Department
of Physics and Earth Sciences, Central Connecticut State University, 1615 Stanley Street, New Britain, CT 06050, USA of Computer Science, Southern Connecticut State University, 501 Crescent Street, New Haven, CT 06515, USA
b Department
Received 4 April 2001; received in revised form 27 August 2001; accepted 25 October 2001
Abstract Atmospheric aerosols exhibit a high degree of variability in their properties and their spatial and temporal distribution. Mapping atmospheric aerosols is of great importance for monitoring the e3ects of anthropogenic activities and natural processes upon local air quality, as well as for producing accurate input for radiative transfer models. To help enhance our understanding of radiative, physical, chemical, and dynamic processes in the atmosphere, laser radar is used to provide systematic monitoring of the temporal evolution of the aerosol rich boundary layer. The Micro Pulse Lidar (Light Detection and Ranging) System is capable of long-term autonomous mapping of vertical aerosol structure. We present a new automated method for tracking boundary layer height from Micro Pulse Lidar data and show the results of this method for a continuous 16 h study conducted on the east coast of the United States. ? 2002 Published by Elsevier Science Ltd. Keywords: Lidar; Aerosols; Boundary layer
1. Introduction An understanding of the spatial and temporal distribution of aerosols is important for atmospheric studies including air quality assessment, pollution management, global climate change research, and precipitation and atmospheric dynamics. The in
0030-3992/02/$ - see front matter ? 2002 Published by Elsevier Science Ltd. PII: S 0 0 3 0 - 3 9 9 2 ( 0 1 ) 0 0 1 0 7 - 4
a Geiger mode avalanche photodiode photon counting detector [1]. The system is capable of continuous, autonomous data acquisition in both daytime and nighttime. Several MPL systems are currently in operation across the country. Each MPL system can generate large volumes of continuous data. For long-term atmospheric monitoring, it is often desirable to employ lidar data in conjunction with data from other instruments [2]. Therefore, to facilitate data analysis and intercomparison of results between instruments, an automated method of tracking boundary layer heights from MPL images is especially important. We have developed a methodology for tracking boundary layer heights with Micro Pulse Lidar data. Various deInitions of average boundary layer height have been used in the literature for derivation of the boundary layer height from lidar backscatter signals. Flamant et al. [3] deIned average boundary layer height as the base of the transition zone (the interface between the mixed layer and the free troposphere)—a deInition associated with a Irst derivative approach. Menut et al. [4] examined methods to calculate average boundary layer height deIned as the middle of the transition zone using second-order derivatives of backscatter signals. Boundary layer height over a continuous time span is calculated in this study as a continuous track based on the gradient of
178
N.C. Parikh, J.A. Parikh / Optics & Laser Technology 34 (2002) 177–185
Fig. 1. Boundary layer track determination
the backscatter signal subject to continuity and curvature constraints. The lidar measurements consist of photoelectrons per microsecond backscattered from a range r. Images formed as a function of time map the signal strength to color intensity at a given range and time. The vertical extent of the aerosol-rich boundary layer is often estimated visually from such images. These manual estimations can be time-intensive and images may be diMcult to interpret, especially in regions of low visual contrast. In addition, interpretation by di3erent individuals may result in subjective di3erences in estimated boundary layer heights. As illustrated in Fig. 1,we applied edge detection, followed by coarse boundary detection to these images. The result provided input for an active contour model of the boundary layer height. The active contour model minimizes an energy functional consisting of terms describing contour curvature, continuity and consistency with image edge points [5].
2. Experimental procedure As a single frequency lidar unit, the Micro Pulse Lidar transmits pulses of 523 nm light at 2500 Hz and measures backscattered photons as a function of time. The signal intensity is a3ected by both extinction and backscatter. The high-resolution instrument we used in this study yielded data with 30 m resolution. The data are corrected for the 1=r 2 range dependence and for instrument calibration factors including detector deadtime, detector afterpulse, and Ield of view overlap as described in Parikh et al. [6]. Each data point consists of a 60 s average of the signal returned from each range. An image of the temporal evolution of the atmosphere is generated from continuous data by plotting the corrected returned signal on a grid where the y-axis represents vertical range above the instrument and the x-axis represents the time of the observation. The signal strength can be represented by a false-color or gray-scale mapping.
Images obtained from a sequence of MPL measurements we conducted at a site o3 the east coast of the United States (near Atlantic City, NJ) over an unbroken 16 h period from 4 pm to 8 am the following morning are shown in Figs. 2a and b. This dataset will be used to illustrate our automated method for boundary layer height determination from MPL data. Sky conditions during the observation period were clear from 4 pm until approximately 10 pm, at which time scattered high-altitude clouds could be seen over the observation site. Shortly after midnight, a low cloud deck settled over the observation site and remained for the duration of the observation period. Visual observations of sky conditions during the period are consistent with sky condition reports from the ASOS Atlantic City, New Jersey Weather Station located approximately six miles southeast of the observation site. Surface wind speed and temperature over the observation period from the ASOS station (shown in Figs. 3a and b) were moderate to mild. The images shown in Figs. 2a and b depict byte-scaled corrected signals mapped into the integer interval between zero and 255. Strong signal returns from overlying dense clouds were mapped into a value of 255 and extremely weak signal returns were mapped into zero, thereby enhancing mid-range intensity values indicative of aerosol variations. 3. Edge detection Canny edge detection techniques [7] using MATLAB were applied to the byte-scaled, mid-range enhanced images described in the previous section to produce a binary edge image. In the Canny edge detection procedure we applied a horizontal Gaussian Ilter to the input images [8]. Non-maximal suppression then removed all locations which were not along a ridge of the gradient. The resulting image was then thresholded to produce a binary edge image in which all image points were either valued at zero (no edge was detected) or 1 (an edge was detected). Parameters that must be speciIed for the edge detection technique include the standard deviation of the Gaussian Ilter and the thresholds used to produce the binary edge image. For our experiments, we used a standard deviation of 1 and a threshold of 0.0125 for weak edges and a threshold of 0.0313 for strong edges. Weak edges were included in the Inal binary edge images if they were connected to strong edges. Figs. 4a and b show the binary edge images for each day overlaid on the original MPL dataset images. The edge detection operation picks out edges corresponding to locations of strong changes in intensity. 4. Coarse boundary detection Once edges have been identiIed it is necessary to discriminate between those edges which potentially correspond to locations representing areas of division between the aerosol
N.C. Parikh, J.A. Parikh / Optics & Laser Technology 34 (2002) 177–185
179
Fig. 2. MPL data (a) day 338, 4:00 –11:59 pm and (b) day 339, 12:00 –8:00 am.
rich boundary layer and the upper atmosphere, and all other edges. In this phase of the procedure, we constructed a coarse boundary consisting of a linked series of blocks that initialized and constrained the Inal boundary produced by the active contour model. The construction of the coarse block boundary required analysis and processing of the original byte-scaled images as well as of the binary edge images. The Irst step in the coarse boundary detection phase of the procedure was to identify blocks within the binary edge image with strong edge persistence over both the altitude
range and time frame represented by the length and width of the block respectively. For each center point of each block, a heuristic measure h was deIned to characterize block edge persistence. For our experiments, h was deIned as h = k(columnedge )
(1)
where k may be any positive constant and columnedge is the total number of columns containing one or more edge points. The next step, altitude block localization, involved locating, for each time frame, the highest altitude block
180
N.C. Parikh, J.A. Parikh / Optics & Laser Technology 34 (2002) 177–185
Fig. 3. Surface observations: (a) Wind speed and (b) Temperature.
with maximal edge persistence. The binary edge image was divided into a sequence of Ixed time frames corresponding to block width. For each time frame the values of the heuristic measure at the center point of each altitude-overlapping block within the time frame were compared to determine the uppermost block associated with the highest heuristic value for that frame. The co-located region in the original byte-scaled image was then checked for cloud amount. The block selected for the given time frame was rejected if more than half of the pixels in the co-located region were identiIed as clouds and the procedure continued until an appropriate block was found. Cloudy pixels were identiIed in the original image as pixels whose value exceeded a threshold value Tc . The rejection of a block in this step of the procedure did not necessarily exclude the block from inclusion in the Inal coarse boundary due to linking criteria applied in the subsequent step. The discontinuous block boundary resulting from altitude block localization served as an initialization for an iterative linking technique which yielded the Inal coarse boundary approximation. Adjustments to the current coarse boundary approximation at locations of discontinuity between boundary blocks were based on the values of the heuristic measure associated with the blocks in the given time frames as well as on continuity with neighboring blocks. The linking procedure terminated when all blocks in the coarse boundary approximation touched any neighboring blocks to the right or left at one or more points. The result of the analysis was a linked series of blocks which corresponded to a coarse estimate of upper altitude constraints on the track of the top of the aerosol rich boundary layer. Figs. 5a and b show the results of the coarse boundary detection technique as a block boundary layer track overlaid over the original byte-scaled data. These results were obtained using a 16 × 16 block size
corresponding to approximately 480 m in altitude and 16 min in time. 5. Active contour The coarse block boundary is reIned using an active contour model. The active contour model requires good initialization. While the midlines of each block in the block contour could be used for starting points, these points may not be edge points. Therefore we applied line integral minimization of an abbreviated energy functional, Einit , to the coarse block midline initialization to snap the contour to a set of edge points where possible, thereby providing an improved initialization. Here Einit is given by Einit = Eimg ds; where
Eimg (q) = 2(es − eq ) − 1
for q = s;
Eimg (s) = −2es − 1;
(2)
where s is the current contour point, q is a potential contour point at the same time and e is the value of the point in the binary edge image [9]. The result is a reIned initial contour estimate shown in Figs. 6a and b. The Inal step to obtain a smooth contour through edge and non-edge points to automatically track the boundary layer height is to minimize the full energy functional, E, constrained to altitude variations. E is given by E = (Econt + Ecurv + Eimg ) ds; where
Econt (q) = (yq − ysprevious )2 Ecurv (q) = (ysprevious − 2yq + ysnext )2 :
(3)
N.C. Parikh, J.A. Parikh / Optics & Laser Technology 34 (2002) 177–185
181
Fig. 4. Binary edge results overlaid on original image (a) day 338, 4:00 –11:59 pm and (b) day 339, 12:00 – 8:00 am.
Here, yq is the altitude of a potential contour point, q, for a given time, t. ysprevious is the altitude of the contour point at time t − 1. ysnext is the altitude of the contour point at time t + 1, and , , and are weights. Finite di3erences are used to approximate Irst and second derivative properties [5]. The energy functional minimizes the weighted sum of a continuity term (an estimate of the Irst derivative property), a curvature term (second derivative property), and an edge term. Note that regions containing cloud pixels may be included in the Inal boundary layer track, if doing so results in the overall minimization of the energy functional.
The Inal result, shown in Figs. 7a and b, is a contour which shows the derived boundary layer height tracked over the time span of the image.
6. Discussion of results The edge detection step, the Irst major component of the technique, was parameterized to yield an abundance of edges including several subtle edges that appeared in the Inal results. For this application the Canny edge detector
182
N.C. Parikh, J.A. Parikh / Optics & Laser Technology 34 (2002) 177–185
Fig. 5. Coarse block contour (a) day 338, 4:00–11:59 pm and (b) day 339, 12:00–8:00 am.
was found to be superior to other standard edge detection operators such as the Sobel operator. Strong edges that were selected by the Canny edge detector and rejected in subsequent steps included edges from cloud shadows, edges from clouds not associated with the top of the boundary layer, and edges demarcating internal structure within the aerosol layer. The result of the second major step of the technique was a connected sequence of Ixed-size regions that formed a coarse outline providing an upper constraint on the Inal boundary layer track. In this coarse, boundary
detection stage of the procedure physical constraints on the data, problem-speciIc expertise, and data-speciIc meteorological information can be incorporated. In addition, adaptation in this stage allows this technique to be extended to analysis of data from other lidar instruments. The Inal boundary layer track was produced using a customized active contour model. The iterative energy functional minimization procedure was terminated when no more than 10% of the contour points were changed in an iteration. Weights on edge, continuity, and curvature terms may be set depending on the degree of desired sensitivity
N.C. Parikh, J.A. Parikh / Optics & Laser Technology 34 (2002) 177–185
183
Fig. 6. Active contour initialization: (a) day 338, 4:00–11:59 pm and (b) day 339, 12:00–8:00 am.
of the contour to speciIc physical features such as aerosol plumes. 7. Conclusions An automated method for tracking the temporal variation of boundary layer height from Micro Pulse Lidar data has been developed. The method has been illustrated on a dataset that spans both day and nighttime conditions and clear and cloudy conditions. Complex internal structure within the
aerosol layer marked by strong edge intensities was also prevalent. The results match well with visual determinations of boundary layer height. Advantages of this automated method are objectivity, consistency, and
184
N.C. Parikh, J.A. Parikh / Optics & Laser Technology 34 (2002) 177–185
Fig. 7. Boundary layer track: (a) day 338, 4:00–11:59 pm and (b) day 339, 12:00–8:00 am.
well suited for adaptation to use with data from other lidar instruments.
The Irst author (NCP) would also like to thank her students, Jennifer Falciglia, Gabriel Ma3ei, and Stephen Price for their assistance in collection of the laser radar data used in this work.
Acknowledgments This work was supported by Connecticut State University research and release time grants. The authors wish to thank Frank Schmidlin at NASA Wallops Flight Facility for the use of the Micro Pulse Lidar system and Steven Palm at NASA Goddard Space Flight Center for helpful discussions.
References [1] Spinhirne JD. Micro pulse lidar. IEEE Trans Geosci Remote Sensing 1993;31(1):48–54. [2] Siebert P, Beyrich F, Gryning S, Jo3re S, Rasmussen A, Tercier P. Review and intercomparison of operational methods for the
N.C. Parikh, J.A. Parikh / Optics & Laser Technology 34 (2002) 177–185 determination of the mixing height. Atmos Environ 2000;34: 1001–27. [3] Flamant C, Pelon J, Flamant P, Durand P. Lidar determination of the entrainment zone thickness at the top of the unstable marine atmospheric boundary layer. Boundary-Layer Meteorol 1997;83: 247–84. [4] Menut L, Flamant C, Pelon J, Flamant P. Urban boundary-layer height determination from lidar measurements over the Paris area. Appl Opt 2000;38(6):945–54. [5] Williams DJ, Shah M. A fast algorithm for active contours and curvature estimation. Comput Vision Graphics Image Process 1992;55(1):14–26.
185
[6] Parikh NC, Schmidlin FJ, Campbell JR, Hlavka DL. Experimental characterization of laser radar instrumentation. J Indian Geophys Union 2000;4(1):91–5. [7] Canny J. A computational approach to edge detection. IEEE Pattern Anal Mach Intell 1986;8(6):679–98. [8] Jain R, Kasturi R, Schunck B. Machine vision. New York: McGraw-Hill, 1995. [9] Niedermeier A, Romaneessen E, Lehner S. Detection of coastlines in SAR images using wavelet methods. IEEE Trans Geosci Remote Sensing 2000;38(5):2270–81.