Environ. Sci. Technol. 2005, 39, 2247-2254
Novel Algorithm for Tomographic Reconstruction of Atmospheric Chemicals with Sparse Sampling W I M V E R K R U Y S S E † A N D L O R I A . T O D D * ,‡ Beckman Laser Institute, University of CaliforniasIrvine, Irvine, California, and Department of Environmental Sciences and Engineering, University of North Carolina at Chapel Hill, Chapel Hill, North Carolina 27599
Numerical studies were performed to evaluate a new air monitoring method for reconstructing chemical exposures and source emissions, based upon optical remote sensing (ORS) and computed tomography (CT). With an ORSCT system, two-dimensional maps of chemical concentrations can be created that have good spatial and temporal resolution. The mathematical algorithm used to compute the distribution is critical for accurate and useable reconstructions of the concentrations. In this research, a novel reconstruction method was tested that uses maximum likelihood expectation maximization (MLEM) combined with two techniques called grid-translation and multi-grid (GT-MG). To evaluate this method, computer simulations were performed using 120 test maps of varying complexity and a simulated ORS system with four instruments and a total of 40 path-integrated measurements. The results were quantitatively compared with two previously used reconstruction methods (single-grid and gridtranslation). Results using the GT-MG method were dramatically improved over previously used methods. Quantitatively, peak exposure errors were reduced by up to 85% and artifacts were dramatically minimized.
Introduction Techniques for monitoring chemicals in air are typically limited to point sampling devices that obtain time weighted averaged or real-time concentration measurements. For source monitoring or human exposure assessment, these devices are usually placed at a limited number of locations throughout an area or on a subset of people from a larger exposed population. These point sampling measurements provide information only about chemical concentrations at a specific location or a person during the sampling time interval. Environmental scientists then use these spatially and/or temporally limited sample results to estimate concentrations in the unsampled areas between the point samples or for the unsampled people in the population. These fixed measurements cannot be used to discern chemical generation or transport patterns unless they are combined with air modeling techniques. A chemical monitoring method has been developed that is a major departure from the traditional methods and generates pollutant maps that allow two-dimensional visualization of chemicals in air. This chemical imaging system fills an important gap in currently used sampling tech* Corresponding author phone: (919) 219-5603; fax: (919) 9664711; e-mail:
[email protected]. † University of CaliforniasIrvine. ‡ University of North Carolina at Chapel Hill. 10.1021/es035231v CCC: $30.25 Published on Web 02/15/2005
2005 American Chemical Society
niques: the maps allow visualization of spatial and temporal changes in concentrations and chemical species across an entire area with high temporal resolution. The maps are generated using an environmental computer-assisted tomography (CAT) scanning system that couples the chemical detection technology of optical remote sensing (ORS) with the mapping capabilities of computerassisted tomography. In theory, a variety of different ORS methods could be used including tunable diode laser and UV-DOAS. The real-world application of the technology used in this research used open-path Fourier transform infrared (OP-FTIR) spectroscopy to measure chemical concentrations across an area. In this system, several OP-FTIR spectrometers transmit beams of infrared light along open beam paths across a space (up to ∼1 km) to simultaneously measure and identify a wide range of multiple contaminants directly in the atmosphere to below part per million (ppm) or part per trillion (ppt) levels. A wide variety of chemicals can be measured as long as they are relatively strong absorbers in the infrared region of the electromagnetic spectrum. Multiple rotating OP-FTIR spectrometers can be set up around the periphery of an area to obtain a noninvasive network of intersecting path-integrated concentrations. Tomographic inversion algorithms are used to transform the network of concentrations into a two-dimensional spatial distribution of concentrations. Every time a new open-path measurement is obtained, an updated two-dimensional concentration map can be generated. When the maps are linked together they can be used to visualize spatial and temporal changes in concentrations and chemical species across an entire space and over time. Outdoors, this method is particularly important for quantifying chemical air emissions and emission rates. An environmental CAT scanning system could be applied to volatile organic compounds, pesticides, and gases from fugitive sources such as chemical and agricultural waste lagoons, landfills, and spray fields. These are large open applications that are difficult to assess with currently used point sampling techniques. Unless many point samples are used, environmental scientists must assume homogeneity over the large areas. Indoors, in a workplace, this system could be used to evaluate ventilation patterns and map exposures to chemicals in the air. CAT has been has been widely applied in medicine and the biological sciences (1, 2). The first discussions of applying CAT to air pollution were in 1979 when it was suggested that a theoretical laser-based system be used on an urban scale (3, 4). It was not until 1990 that CAT was proposed to measure pollutants in indoor air (5). This technology has now been evaluated in theoretical and experimental indoor (6- 11) and outdoor field studies (12-14). For practical reasons, the typical number of pathintegrated measurements used in CAT applied to environmental imaging is much smaller (10,000). This relatively sparse sampling is due to several factors including the relatively low speed and high cost of remote sensing technology coupled with spatially and temporally fluctuating chemical concentrations in air. In addition to sparse sampling, the variety of indoor and outdoor field conditions places challenges on the spatial arrangement of optical rays (configuration) and usually results in nonsymmetrical configurations. Therefore, algorithms developed for medicine are not usually optimal or directly applicable to reconstructing chemical pollutants in air. VOL. 39, NO. 7, 2005 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
2247
FIGURE 1. Configuration of optical rays and OP-FTIR spectrometers used in the experiments. The squares in the corners (four) represent the spectrometers and the solid lines represent the path of the optical rays (10 for each source). The dashed lines indicate the width of one of the rays. Several algorithms have been investigated for reconstructing concentrations from a limited number of measurements. In particular, iterative algebraic reconstruction techniques have been found to be successful for generating chemical concentration maps: algebraic reconstruction technique (ART), multiplicative algebraic reconstruction technique (MART), and maximum likelihood expectation maximization (MLEM) (15-17). Several nonalgebraic reconstruction methods have also been evaluated with good success (12, 18-20). The algebraic reconstruction techniques are quick and allow real-time reconstructions of concentrations; however, as the concentration distribution becomes complex relative to the configuration, artifact production increases. This paper describes an evaluation of a novel iterative reconstruction method to reconstruct chemical concentrations in air that provides a considerable improvement in both the accuracy of reconstructing complex concentration distributions and the minimization of artifact production. This method does not use any a priori information on chemical distributions and uses MLEM combined with two techniques called grid-translation and multi-grid (GT-MG). To evaluate this method, computer simulations were performed using 120 test maps of varying complexity and a simulated optical remote sensing configuration with four instruments and a total of 40 path-integrated measurements. The reconstructed maps were evaluated based upon quantification of concentrations over the entire map and at specific peak locations, and based upon visual production of smooth distributions and artifacts.
Methods Test Concentration Data. To evaluate the reconstruction method (described below), 120 test maps were generated by randomly distributing up to six concentration peaks on a 40 × 40 grid. These maps were similar to those that have been used in previous studies (15). Of the 120 maps, 20 maps had 1 peak, 38 maps had 2 peaks, 28 maps had 3 peaks, 18 maps had 4 peaks, 14 maps had 5 peaks, and 2 maps had 6 peaks. A 40 × 40 grid could represent an area that is 40 m on one side, which could be used for reconstruction of chemicals outdoors. The concentration peaks were defined approximating a bi-variate Gaussian source (eq 1)
C(x, y) ) P2248
9
(x-x0)2+(y-y0)2 σ2
(1)
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 39, NO. 7, 2005
where C(x, y) is the concentration in a grid cell, P is the peak height (0 to 40 ppm), x0 and y0 are the location coordinates (0 to 40), and σ is the width of the peak (2.8, 4.2, 5.7, or 7.1). Optical Remote Sensing Configuration. In all of the experiments, a configuration (or arrangement) of open-path rays was used that simulated four OP-FTIR spectrometers (one in each corner of the area) with 10 optical rays per OP-FTIR spectrometer (40 total rays), see Figure 1. The simulated ray width is depicted by dashed lines for one ray. The rays were distributed such that they were separated by angles of 15°. The inputs used by a reconstruction algorithm, from an actual optical remote sensing system (ORS), are the individual chemical concentrations (ray sums) measured by each openpath ray. When using an OP-FTIR spectrometer, the concentrations are usually expressed as parts per millionmeters (ppm-m), which is the path-integrated concentration along the length of a one-meter ray. This is converted to ppm by dividing this number by the actual length of the ray. With a simulated system, the algorithm uses ray sums (reprojection values) that are calculated from the test maps per eq 2 N×M
∑AC
ij i
pj )
i)1
(2)
W
where pj is the ray sum for ray j, Aij is the fractional area of the i-th grid cell intercepted by the j-th measured ray sum, Ci is the grid cell concentration (ppm), N × M is the number of grid cells (40 × 40), and W is the width of the ray. Traditional MLEM Reconstruction Algorithm. The maximum likelihood with expectation maximization (MLEM) tomographic reconstruction algorithm was used as the basic reconstruction algorithm in this study to reconstruct the simulated ray sums that would have been achieved using an optical remote sensing device (17, 21-22). To reconstruct a concentration map, an idealized area is broken into an N × M square or rectangular grid of cells. This reconstruction grid has a fixed, single basis resolution (N × M). The concentration in each cell is assumed to be homogeneous and nonnegative, and is updated with each iteration. In this paper we will refer to this method as single grid (SG). The initial concentration values used in the reconstructions are called “smart guesses” and are determined based upon eq 3
∑pj
j
j
Ci )
(3)
ki
where p j j is the average concentration for pj (equaling the ray sum value divided by the length of the ray) and ki is the number of times cell Ci is hit by a ray. For each MLEM iteration, the concentrations Cqi are updated using eq 4
Cq+1 i
)
Cqi
∑
Aij
Aijp ˆj
∑p j
q
(4)
j
j
where p ˆ j is the set of measured (or simulated) ray sums, and Aij is the intersection areas of ray j with reconstruction grid cell Ci. The reprojection values after q iterations, pqj , are calculated with eq 2, now using the concentration values Cqi and intersection areas corresponding to the reconstruction grid. As in previous studies, a fixed number of iterations (50)
FIGURE 2. Example of generation of eight additional reconstruction grids from a basis grid (a) of 2 × 2 cells. Translation in the x-direction of the basis grid by 1/3and 2/3 of a basis grid cell dimension results in the grids shown in (b) and (c), respectively. Parts (d), (e), and (f) result from translation in the y-direction by 1/3 the basis grid cell dimension, combined with an x-direction translation by 0, 1/3, and 2/ , respectively. Similarly, the grids in (g), (h), and (i) result from 3 a y-direction translation of 2/3, combined with the x-direction translations. were used for all of the reconstructions. This number was found to provide reasonable reconstructions for most maps. In almost all cases, major improvements were achieved within 10 iterations; improvement was minimal at 40-50 iterations. Grid Translation Method. One limitation in the traditional method of choosing a fixed single reconstruction grid resolution is that the location of the edges of each of the grid cells in relation to a particular peak position and shape and optical ray position could adversely impact reconstruction accuracy. A grid translation (GT) method was developed to be used with MLEM to improve upon these limitations (23). With the GT method, reconstructions are performed on a number of different grids rather than on a fixed single grid. The basis grid is the same as used in the single grid method and the additional grids are formed by translating, in the two orthogonal directions described by the edges of the grid cells, over distances equal to a fraction of a single basis grid cell size. In the example of Figure 2, a basis grid of 2 × 2 grid cells is used to generate eight additional grids by translating over distances equal to multiples of 1/3 the size of a basis grid cell. After the determined number of MLEM iterations on each of the translated grids, the resulting “interim” concentration maps of each of the grids are combined into one final map of higher resolution. For the results of the GT method presented here, 50 iterations of MLEM were used for each grid. The combined maps match the original test maps significantly better than maps reconstructed using only a single grid. Moreover, peak concentration value estimates are more accurate with the GT method when compared to the SG method. For the grid translation-multi grid (GT-MG) method (described below), the GT portion was altered. Rather than waiting until all of the individual maps reached 50 iterations to create the combined higher resolution map, a higher resolution map was created after each individual iteration. This higher resolution map was then used as input for the starting values for the second (next) iteration for each translated map. For the results of the GT method in this paper, we used a basis grid of 8 × 8 cells and translation distances equal to multiples of 1/7 the width of a basis grid cell; this resulted in
FIGURE 3. Interim reconstruction maps resulting from the GT-MG method corresponding to basis grids with resolutions of 2 × 2 cells up to 10 × 10 cells (M2 to M10, respectively). a total of 49 translated grids. The smaller the translation distance, the higher the number of combined maps. As the number of combined maps increases, resolution and reconstruction accuracy increases. However, at some point, the improvement achieved by using more translated grids (smaller translation distances) becomes less dramatic while the computational cost continues to increase. Grid Translation Combined with Multi Grid (GT-MG). The new method presented in this paper combines the grid translation (GT) method with a technique referred to as the multi-grid (MG) method (24-26). The grid translation-multi grid (GT-MG) method starts with a very coarse resolution 2 × 2 basis grid. Grid translation is then used to generate additional grids starting with this 2 × 2 grid. After the determined number of translations are performed using MLEM, a final combined higher resolution map (M2) is obtained for the initial 2 × 2 basis grid (Figure 3). VOL. 39, NO. 7, 2005 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
2249
TABLE 1. Grid Translation Distances and the Resulting Total Numbers of Grids Used for Each Basis Grid Resolution basis grid cell resolution
basic translation distance (units of basis grid cell size)
total number of grids in set
2×2 3×3 4×4 5×5 6×6 7×7 8×8 9×9 10 × 10
1/8 1/8 1/8 1/8 1/7 1/6 1/5 1/4 1/4
64 64 64 64 49 36 25 16 16
Next, the basis grid resolution is increased to 3 × 3 cells. The initial concentration values used in the 3 × 3 grid are calculated from the final combined map generated from the 2 × 2 basis grid (M2). Concentrations are selected from the cells in M2 that are closest in location to the cells in the new basis grid. Iterations are then performed on the new set of grids, resulting in a new final combined map (M3). Figure 3 shows the improvement in reconstruction using this method up to a basis resolution of 10 × 10 (M2 to M10). This research used an initial basis grid resolution of 2 × 2 cells that was incrementally increased up to 10 × 10. Translation distances and corresponding number of grids used for each basis grid resolution are listed in Table 1. For the GT-MG method, the number of iterations was not fixed at 50; it could vary up to a maximum of 50 iterations. To determine the number of iterations, a statistic was calculated called projection distance (PD) that measured how closely the reconstructed ray sums matched the original ray sums (see eq 5) (27). A projection data distance of zero implies a perfect match. The program stopped iterating when the PD was smaller than 0.5%, or a maximum of 50 iterations was reached. R
q
PD )
∑(pˆ - p j
j)1
q 2 j
)
(5)
R
where R is the total number of rays, pˆj is the original ray sum for ray j, and pqj is the reconstructed ray sum (reprojection value) after q iterations. Evaluating Reconstruction Quality. The GT-MG method was evaluated quantitatively in several ways. A conventional image quality measure, called nearness, was used to describe the discrepancy between the original test maps and the reconstructed test maps (9, 28). Nearness evaluated errors over all the grid cells in the map (eq 6).
Nearness )
x
40×40
∑ (C
/ i
- Ci)2
i
(6)
40×40
∑
(C/i
-C
/
2
avg)
i
where C/i and Ciare the true value and the estimated value, respectively, for the i-th cell in the map, and C/avg is the mean concentration of all the grid cells in the “original” map. The lower the nearness value, the better the match. A second evaluation parameter used was exposure error percent (eq 7) which evaluated how well average concentrations in an area around a peak are reconstructed. In industrial hygiene, this could reflect the accuracy of performing a worstcase exposure assessment for a worker or of tracking the movement of volatile organic chemicals in air in a workspace. 2250
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 39, NO. 7, 2005
FIGURE 4. Mean nearness values shown in quartiles for SG, GT, and GT-MG methods, categorized by map complexity. The box shows the errors in the middle two quartiles and the whiskers show the top and bottom of the quartiles. In the environmental field, this could reflect the accuracy of measuring chemical air emissions and emission rates from fugitive sources such as chemical and agricultural waste lagoons, landfills, and spray fields. AOI
AOI
∑C - ∑C
Exposure error % )
|
/ i
i
i
AOI
∑
C/i
|
i
× 100%
(7)
i
where C/i and Ci are the true value and the estimated value, respectively, for the i-th cell in the map, and AOI is the area of interest (cells) around the peak being evaluated. For this study, the area of interest was all the cells within a peak width σ (see eq 1). It can be anticipated, on the basis of sampling theory, that well-sampled peaks (hit by relatively many rays) will be reconstructed more accurately than peaks that are sampled by relatively few rays. To investigate the impact of sampling density on the accuracy of peak concentration reconstruction, using the current configuration, peaks were categorized by the number of rays that crossed the peaks. Three categories were used: peaks sampled by more than 9 rays (16 peaks), peaks sampled by 5-9 rays (187 peaks), and peaks sampled by fewer than 5 rays (131 peaks).
Results and Discussion Overall Image Quality. For the 120 test maps, regardless of the number of peaks in the test map, the average nearness using the GT-MG method was statistically smaller (better) than the average nearness using the conventional single grid or grid translation reconstruction method (P < 0.05) (Figure 4 and Table 2). Nearness deteriorated as the number of peaks in the maps increased using any of the three methods (p < 0.05). On the basis of the test maps used in this research, the GT-MG method reduced the nearness values by more than 50% of the values previously obtained using the single grid method for all categories of map complexity. This improvement is important because nearness is a measure of the accuracy of reconstruction over the entire reconstructed map; therefore, it represents reconstruction of peak heights, shapes, and the production of artifacts. Figure 5 shows reconstructions obtained using the SG, GT, and GT-MG methods for four different test maps with two, three, five, and six peaks (5a, 5b, 5c, and 5d, respectively). The maps reconstructed using the GT-MG method (right column) are much closer to the original test maps in peak
FIGURE 5. Original test maps (1st column) and corresponding reconstructed test maps using the SG (2nd column), GT (3rd column), and GT-MG (4th column) methods.
TABLE 2. Mean, Median, and Standard Deviation of Nearness nearness
mean
median
SD
SG 1 GT 1 GT-MG 1 SG 2 GT 2 GT-MG 2 SG 3 GT 3 GT-MG 3 SG 4 GT 4 GT-MG 4 SG 5/6 GT 5/6 GT-MG 5/6
0.63 0.32 0.09 0.57 0.38 0.16 0.63 0.42 0.19 0.68 0.47 0.25 0.72 0.57 0.37
0.62 0.27 0.09 0.54 0.38 0.14 0.63 0.41 0.19 0.72 0.45 0.24 0.72 0.56 0.36
0.15 0.11 0.05 0.11 0.09 0.07 0.14 0.09 0.06 0.13 0.09 0.08 0.09 0.05 0.09
height and shape compared with the SG (2nd column) and the GT method (3rd column). Using the GT-MG method, artifacts are dramatically reduced or eliminated; this results in concentration distribution shapes that are smooth and peak locations that can be easily identified. The dramatic improvement in both the nearness and in the visual appearance of the concentration maps is very important for real-time assessment of chemical concentrations in environmental and industrial applications such as locating chemical leaks and evaluating human exposures.
FIGURE 6. Absolute exposure error percentages shown in quartiles for SG, GT, and GT-MG methods, categorized by map complexity. The box shows the errors in the middle two quartiles and the whiskers show the top and bottom of the quartiles. Exposure Error Percent. Regardless of the number of peaks in the original test map, the mean absolute exposure error percent using the GT-MG method was statistically smaller (better) than the error percent using the conventional single grid reconstruction method as well as the GT method (P < 0.05) for all categories of number of peaks per map (Figure 6 and Table 3). Using the GT-MG method, the mean exposure error was reduced by up to 85% from using the SG method, when only a few peaks were present in the maps. VOL. 39, NO. 7, 2005 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
2251
TABLE 3. Mean, Median, and Standard Deviation of Exposure Errors exposure error
mean
median
SD
SG 1 GT 1 GT-MG 1 SG 2 GT 2 GT-MG 2 SG 3 GT 3 GT-MG3 SG 4 GT 4 GT-MG 4 SG 5/6 GT 5/6 GT-MG 5/6
22.1 6.8 2.9 20.6 14.1 7.9 29.3 21.8 14.6 30.9 22.1 17.2 34.1 28.5 23.3
12.4 5.7 2.5 17.6 10.2 6.0 26.7 20.9 13.5 27.8 23.1 16.3 35.6 30.4 26.4
20.3 7.3 1.8 13.5 10.0 6.8 17.8 11.2 9.7 13.7 9.7 8.9 16.1 14.3 12.3
FIGURE 8. Histograms of exposure error percentages for all peaks (n ) 334), stratified by the number of rays by which each peak area is sampled: (a) SG method, (b) GT method, and (c) GT-MG. The dashed line is the distribution of exposure errors for areas hit by more than 9 rays (n ) 131), the thin solid line is for areas hit by 5-9 rays (n ) 187), and the thick solid line is for areas hit by 10 or more rays (n ) 16).
FIGURE 7. Histograms of exposure error percentages for all peaks (n ) 334) in all 120 test maps for the SG method (a), GT method (b), and GT-MG (c). The histograms were made using exposure error ranges of 20%. Positive errors correspond to overestimation of the average concentration in the peak area, and negative values correspond to underestimation. 2252
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 39, NO. 7, 2005
For one to two peaks, the mean reconstructed concentrations were within 8% of the actual concentrations. To investigate if exposure was under- or overestimated, the relative exposure errors for all the individual peaks (n ) 334) in all 120 test maps were plotted as histograms (Figure 7). For all three methods, the exposure errors were generally negative; therefore, the average concentrations around the peaks were underestimated. The distribution of errors for the SG method (7a) was the broadest. The tighter the distribution of errors, the greater the confidence and predictability of the peak concentrations. The distribution of errors was then categorized by the number of rays that sampled each peak (Figure 8). As the number of sampling rays increased, the distribution of errors narrowed and decreased (were closer to 0%) for all three methods. With the GT-MG method, the exposure errors for the peaks sampled
FIGURE 9. Original test map (1st column) and corresponding reconstructed test maps using the SG, GT, and GT-MG methods for a complex distribution with narrow and overlapping peaks.
FIGURE 10. Original test map (1st column) and corresponding reconstructed test maps using the GT-MG method for a complex distribution with narrow peaks with double the number of optical rays.
FIGURE 11. Original test map (1st column) and corresponding reconstructed test maps using the SG, GT, and GT-MG methods for a “realistic” distribution. by more than 9 rays decreased to (10%. As expected, reconstruction accuracy increased as the number of rays sampling each peak increased. An under-sampled peak is flattened out in the reconstruction process and consequently the exposure in the peak area is underestimated. The distributions of exposure errors in Figure 7 are only for the current set of test maps and optical remote sensing configuration. Additional simulations were performed using a denser configuration with 60 rays and the shapes of the distributions in Figure 8 were similar to the distributions obtained using 40 rays. In terms of reconstruction accuracy, the GT-MG method represents a significant improvement over the SG and GT method. While the GT-MG method is an important improvement over the previous methods, the effectiveness of an algorithm is limited by the density of the optical remote sensing configuration, the complexity of the concentration profile, and the speed of data acquisition compared with the rate of plume movement. With the test maps used in this study, the GT-MG and the other reconstruction methods were less effective in reconstructing individual peaks in complex profiles with over four peaks. However, the GT-MG method was always an improvement over the previously used reconstruction methods. With the configuration used in this research study, reconstruction accuracy decreased for maps with multiple very narrow peaks or multiple overlapping peaks. For multiple peaks, the very narrow peaks may fall between the rays and a denser configuration may be required to differentiate between overlapping peaks. This is not a drawback of the reconstruction method; it is due to a low sampling frequency relative to the spatial frequencies of the
concentration maps. Figure 9 shows reconstructions of a map with five overlapping peaks using the SG, GT, and GTMG methods. Visually, the GT-MG map gave the best reconstruction; however, two of the five peaks were not well differentiated. The SG, GT, and GT-MG methods had nearness values of 0.796, 0.636, and 0.518, respectively. When the number of rays in the configuration was doubled, nearness improved by over 30%. Figure 10 shows the improved reconstruction and peak differentiation for the GT-MG method using a higher density configuration. A two-dimensional dynamic flow model was used to create maps of single pollutants in a room.29 This model was used to represent “realistic” chemical distributions. Some of the peaks represent potentially difficult scenarios to reconstruct, using the configuration used in this research, because the peaks were very narrow. For example, in Figure 11, the location of the peak was accurately reconstructed; however, the nearness was high at 0.547, 0.424, and 0.307, for the SG, GT, and GT-MG methods, respectively. The exposure errors were also high at 61%, 58%, and 46%, for the SG, GT, and GT-MG methods, respectively. When the configuration was doubled in density the reconstruction improved visually and nearness improved by 20%. The results presented in this paper show that high-quality reconstruction maps can be obtained with minimal artifacts at moderate computational cost when an algebraic iterative method is combined with grid-translation and multi-grid techniques. Peaks in maps with multiple concentration profiles could easily be identified and the associated exposure errors were low. Considerable reductions in computation time can be achieved by reducing the number of increases VOL. 39, NO. 7, 2005 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
2253
in basis grid resolution with only moderate deterioration of image quality or by using fewer translated grids at each basis grid resolution. Optimization of the code may further reduce computation time. These improvements in tomographically reconstructing concentration profiles using relatively sparsely sampled remote sensing systems provide industrial hygienists and environmental scientists with a powerful tool for evaluating human exposures, air flow patterns, emission measurements, and validating fluid flow models. In addition, these improvements may enable simpler configurations to be used in areas where lower resolution is acceptable; this would significantly reduce the hardware cost and setup time. As in the medical field, a natural extension of this work in the future is expanding it from two dimensions to three dimensions. Presently, in two dimensions, a single plane of interest must be selected for monitoring. Outdoors, this system is usually placed adjacent to the emission source and indoors it is placed near the breathing zone of a worker. As ORS technology becomes faster, cheaper, and smaller, a system can be set up to monitor chemicals through a volume (three dimensions) of air allowing us to completely visualize chemical flows and emissions. For example, sensors could be placed around the periphery of a large indoor work area to map concentration flow around workers who are seated or standing. While ORS technology is the limiting factor for this expansion, algorithms would need to be modified as well. This work is underway.
(11)
(12)
(13)
(14)
(15)
(16)
(17)
(18)
(19)
Acknowledgments This material is based upon work supported by the National Science Foundation under grant 00011385. We would like to thank Dr. Kathleen Mottus for her technical assistance with the manuscript.
Literature Cited (1) Housfield, G. N. Computerized Transverse Axial Scanning (Tomography). 1. Description of the system. Br. J. Radiol. 1973, 46, 1016. (2) DeRosier, D. J.; Klug, A. Reconstruction of three-dimensional structures from electron micrographs. Nature 1968, 217, 130. (3) Byer, R. L.; Shepp, L. A. Two-dimensional remote air-pollution monitoring via tomography. Appl. Optics Lett. 1979, 4, 75. (4) Wolfe, D. C.; Byer, R. L. Model studies of laser absorption computed tomography for remote air pollution measurement. Appl. Optics 1979, 21, 1165. (5) Todd, L. A.; Leith, D. Remote sensing and computed-tomography in industrial hygiene. Am. Ind. Hyg. Assoc. J. 1990, 51, 224. (6) Yost, M. G.; Gadgil, A. J.; Drescher, A. C.; Zhou, Y.; Simonds, M. A.; Levine, S. P.; Nazaroff, W. W.; Saison, P. A. Imaging indoor tracer-gas concentrations with computed-tomography: experimental results with a remote-sensing FTIR system. Am. Ind. Hyg. Assoc. J. 1994, 55, 395. (7) Drescher, A. C.; Park, D. Y.; Yost, M. G.; Gadgil, A. J.; Levine, S. P.; Nazaroff, W. W. Stationary and time-dependent indoor tracergas concentration profiles measured by OP-FTIR remote sensing and SBFM-computed tomography. Atmos. Environ. 1997, 31, 727. (8) Bhattacharya, R.; Todd, L. A. Spatial and temporal visualization of gases and vapours in air using computed tomography. Numerical studies. Ann. Occup. Hyg. 1997, 41, 105. (9) Samanta, A.; Todd, L. A. Mapping air contaminants indoors using a prototype computed tomography system. Ann. Occup. Hyg. 1996, 40, 675. (10) Todd, L. A.; Farhat, S. K.; Mottus, K. M.; Mihlan, G. J. Experimental evaluation of an environmental CAT scanning system for
2254
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 39, NO. 7, 2005
(20)
(21)
(22)
(23)
(24) (25) (26)
(27) (28)
(29)
mapping chemical in air in real-time. Appl. Occup. Environ. Hyg. 2001, 16, 45. Fischer, M. L.; Price, P. N.; Thatcher, T. L.; Schwalbe, C. A.; Craig, M. J.; Wood, E. E.; Sextro, R. G.; Gadgil, A. J. Rapid measurements and mapping of tracer gas concentrations in a large indoor space. Atmos. Envir. 2001, 35, 2837. Cuccoli, F.; Facheris, L.; Tanelli, S.; Giuli, D. Infrared tomographic system for monitoring the two-dimensional distribution of atmospheric pollution over limited areas. IEEE Trans. Geosci. Remote Sens. 2000, 38, 1922. Piper, A.; Todd, L.; Mottus, K. A field study using open-path FTIR spectroscopy to measure and map air emissions from volume sources. Field Anal. Chem. Technol. 1999, 3, 69. Todd, L. A.; Ramanathan, M.; Mottus, K.; Katz, R.; Dodson, A.; Mihlan, G. Measuring chemical emissions using open-path Fourier transform infrared (OP-FTIR) spectroscopy and computer assisted tomography. Atmos. Environ. 2001, 35, 1937. Todd, L. A.; Ramachandran, G. Evaluation of optical sourcedetector configurations for tomographic reconstruction of concentration in indoor air. Am. Ind. Hyg. Assoc. J. 1994, 55, 1133. Todd, L. A.; Bhattacharyya, R. Tomographic reconstruction of air pollutants: evaluation of measurement geometries. Appl. Optics 1997, 30, 7678. Samanta, A.; Todd, L. A. Mapping chemicals in air using an environmental CAT scanning system: evaluation of algorithms. Atmos. Environ. 2000, 34, 699. Drescher, A. C.; Gadgil, A. J.; Price, P. N.; Nazaroff, W. W. Novel approach for tomographic reconstruction of gas concentration distribution in air: Use of smooth basis functions and simulated annealing. Atmos. Environ. 1996, 30, 929. Price, P. N.; Fischer, M. L.; Gadgil, A. J.; Sextro, R. G. An algorithm for real-time tomography of gas concentrations using prior information about spatial derivatives. Atmos. Environ. 2001, 35, 2827. Park, D. Y.; Fessler, J. A.; Yost, M. G.; Levine, S. P. Tomographic reconstruction of tracer gas concentration profiles in a room with the use of a single OP-FTIR and two iterative algorithms ART and PWLS. J. Air Waste Manage. Assoc. 2000, 50, 357. Tsui, B. M. W.; Zhao, X. D.; Frey, E. C.; Gulberg, G. T. Comparison between ML-EM and WLS-CG algorithms for spect imagereconstruction. IEEE Trans. Nucl. Sci. 1991, 38, 1766. Shepp, L. A.; Vardi, Y. Maximum likelihood reconstruction in positron emission tomography. IEEE Trans. Med. Imag. 1982, 1, 113. Verkruysse, W.; Todd, L. A. Improved Method (“Grid Translation”) For Mapping Environmental Pollutants Using A TwoDimensional Cat Scanning System. Atmos. Environ. 2004, 38, 1801. Fedorenko, R. P. On the speed of convergence of an iterations process. USSR Comput. Math Phys. 1964, 4, 227. Brandt, A. Multi-level adaptive solutions to boundary-value problems. Math. Comput. 1977, 31, 333. Pan, T. S.; Yagle, A. E. Numerical study of multigrid implementations of some iterative image reconstruction algorithms. IEEE Trans. Med. Imag. 1991, 10, 572. Herman, G. T.; Lent, A. Iterative reconstruction algorithms. Comput. Biol. Med. 1976, 6, 273. Herman, G. T.; Lent, A.; Rowland, S. ART: Mathematics and applications (a report on the mathematical foundations and on the applicability to real data of the algebraic reconstruction techniques). J. Theor. Biol. 1973, 42, 1. Yamamoto, T.; Sparks, L. E. USEPA RVMPC Version 1.0; EPA: Washington, DC, 1994.
Received for review November 4, 2003. Revised manuscript received December 22, 2004. Accepted December 22, 2004. ES035231V
