- Email: [email protected]
Impact of sensor characteristics on source characterization for dispersion modeling
Measurement 44 (2011) 802–814
Contents lists available at ScienceDirect
Measurement journal homepage: www.elsevier.com/locate/measurement
Impact of sensor characteristics on source characterization for dispersion modeling Luna M. Rodriguez a,b,c, Sue Ellen Haupt a,b,c,⇑, George S. Young b a
P.O. Box 30, Applied Research Laboratory, The Pennsylvania State University, State College, PA 16804-0030, United States 503 Walker Building, Meteorology Department, The Pennsylvania State University, University Park, PA 16802, United States c P.O. Box 3000, The National Center for Atmospheric Research/Research Applications Laboratory, Boulder, CO 80307-3000, United States b
a r t i c l e
i n f o
Article history: Received 9 September 2010 Received in revised form 7 January 2011 Accepted 20 January 2011 Available online 2 February 2011 Keywords: Sensor thresholds Sensor saturation Dispersion monitoring Dispersion modeling Source characterization Source term estimation Genetic algorithm
a b s t r a c t An accidental or intentional release of hazardous chemical, biological, radiological, or nuclear material into the atmosphere obligates responsible agencies to model its transport and dispersion in order to mitigate the effects. This modeling requires input parameters that may not be known and must therefore be estimated from sensor measurements of the resulting concentration field. The genetic algorithm (GA) method used here has been successful at back-calculating not only these source characteristics but also the meteorological parameters necessary to predict the contaminants subsequent transport and dispersion. This study assesses the impact of sensor thresholds, i.e. the sensor minimum detection limit and saturation level, on the ability of the algorithm to back-calculate modeling variables. The sensitivity of the back-calculation to these sensor constraints is analyzed in the context of an identical twin approach, where the data is simulated using the same Gaussian Puff model that is used in the back-calculation algorithm in order to analyze sensitivity in a controlled environment. The solution is optimized by the GA and further tuned with the Nelder–Mead downhill simplex algorithm. For this back-calculation to be successful, it is important that the sensor capture the maximum concentrations. Ó 2011 Elsevier Ltd. All rights reserved.
1. Introduction When a toxic agent is released into the atmosphere, government agencies are tasked with tracking and forecasting its transport and dispersion [1–4]. These scenarios include both accidental and intentional release of any chemical, biological, radiological, or nuclear (CBRN) contaminant that can be a threat to human and ecological health. Although tracking the contaminant can be accomplished with arrays of sensors, forecasting the transport and dispersion requires input variables describing the source of the release and the meteorological conditions. Sometimes, however, there is inadequate a priori source or meteorological information to perform those predictions and it becomes necessary to characterize the source ⇑ Corresponding author. Tel.: +1 303 497 2763; fax: +1 303 497 8386. E-mail addresses: [email protected] (L.M. Rodriguez), haupts2@ asme.org (S. Ellen Haupt), [email protected] (G.S. Young). 0263-2241/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.measurement.2011.01.014
of an airborne contaminant from the sensor field measurements of the resulting concentration field. Such source term estimation involves back-calculating the source location and emission rate history. There has been extensive work on methods for back-calculating source characteristics. For example, Thompson et al. [5] use a random search algorithm with simulated annealing and Rao [6] reviews methods that include adjoint and tangent linear models, Kalman filters, and variational data assimilation, among others. Some previous work uses genetic algorithms (GAs) to optimize source characteristics [7–13]. In addition to characterizing the source, some of these studies also back-calculate the meteorological variables required for transport and dispersion modeling, including wind speed, wind direction, and atmospheric stability. Several of the aforementioned papers examine the effect of adding noise to the data to simulate errors in the sensor data as well as the uncertainty inherent in atmospheric turbulence. Most studies, however, do not take into account the sensor
803
L.M. Rodriguez et al. / Measurement 44 (2011) 802–814
constraints, such as the minimum detection threshold and saturation levels. The sensor characteristics that we seek to evaluate are the detection and saturation levels, not only because they exist in real sensors but also to assess how they influence Sensor Data Fusion (SDF) algorithms, specifically our method discussed in Section 3. The understanding of sensor limitations is vital for the military and homeland defense agencies given that sensors used in the field vary in their saturation and detection capabilities. These capabilities can vary from simple binary (on/off) sensors to more complex sensors whose detection and saturation capabilities depend on the gain settings used. For example, the range specifications for the Digital Fast Response Photo-Ionization Detector (digiPID) sensors used in the FUsing Sensor Information from Observing Networks (FUSION) Field Trial (FFT07) are displayed in Table 1. At the lowest gain setting of 0, the sensor saturation is 1000 ppm. At this setting, the response is fast (0.25 ms) but noisy [14]. For higher gain settings, saturation occurs at lower concentration values, but noise is reduced. The threshold level for the propylene released in that experiment is 50 ppb. Due to expense and difficulty of calibration, some networks may include sensors that only supply binary data. For example, the Joint Warning and Reporting Network (JWARN) is a US DoD system that employs binary sensors for early warning and reporting functions [15]. It is precisely such systems that may also be used to back-calculate source variables. If sensors with a high threshold are not able to detect a contaminant it is obvious that many lives may be in danger; however, if the sensors are only able to detect over a certain threshold this again may affect how SDF algorithms interpret the spread of the puff and consequently the source characteristics. Fig. 1 is an example of how a saturation level of one order of magnitude smaller than the maximum concentration detected and a detection level three orders of magnitude smaller than the maximum concentration level detected truncates the concentration field. Much of the information is lost. In real applications SDF algorithms will use the truncated concentration field detected from these sensors. It is important that we understand these highly truncated scenarios even when using sensors with a large dynamic range. In a scenario where concentrations are high, if these sensors become saturated immediately, the concentration field will appear binary to the SDF algorithm. Also, if sensors become saturated the concentration gradient is altered; as a result we lose information about the path and spread of the contaminant puff, which can also affect the retrieval of source characteristics.
Table 1 Sensor specifications for Digital Fast Response Photo-Ionization Detector (digiPID) used in FFT07. Frequency response: 50 Hz. Detection limit: 50 ppb (propylene). Gain setting
Gain factor
Full scale conc. (PPM)
Integration time (MSEC)
Number of integrations/sample
0 1 2 3
1 4 16 32
1000 250 62 31
0.25 1.00 4.00 8.00
64 16 4 2
The goal of the present study is to determine the sensor performance characteristics required to properly back-calculate the source parameters and the wind direction. This is done by first creating ‘‘true data’’ using a Gaussian Puff model and then using a genetic algorithm (GA) coupled with the same Gaussian Puff model to optimize the source and meteorological parameters for agreement with these simulated observations: this is known as the identical twin approach. By applying thresholds to the concentration data, sensor detection and saturation constraints are modeled. In this way the sensitivity of the source characterization to the sensor performance can be evaluated in a controlled and systematic manner. For this study, the source characteristics (x, y) location of the release and source strength) plus the wind direction are back-calculated by a GA, and then fine-tuned to the global optimum using the Nelder–Mead downhill simplex algorithm (NMDS) [16]. The GA coupled with the NMDS will be referred to as a Hybrid GA. The process is repeated for varying numbers of sensors to estimate how many sensors are required to perform the back-calculation for each thresholding scenario. This sensitivity analysis is needed because sensor thresholds limit the amount of information contained in the data. Moreover, because the algorithm works by comparing model predictions to the observations it will perform best if it is aware of these thresholds so that it can model the observed data instead of an ideal distribution. To a larger extent, this study assesses the way that sensor design impacts utility. Specifically, does a low saturation level or high detection level prevent the data from being used for this source characterization? The sensitivity analysis of sensor constraints including data descriptions and procedures appear in Section 2. Section 3 presents the results. Section 4 is a case study using a more refined model to create the test data. Section 5 summaries the conclusions and discusses the implications of the study. 2. Model construction We construct an identical twin experiment with Gaussian Puff model output being used to back-calculate source strength, location of the release, and wind direction. We first describe how the concentration data is generated, then discuss the application of thresholds to the data to simulate the sensor constraints, and finally describe the Hybrid GA used to do the back-calculation. 2.1. The Gaussian Puff model The Gaussian Puff model is used to produce the synthetic concentration observations:
Cr ¼
Q Dt
exp
ðxr UtÞ2 2r2x
!
ð2pÞ1:5 rx ry rz !" ! y2r ðzr He Þ2 exp exp 2r2y 2r2z !# ðzr þ He Þ2 þ exp 2r2x
ð1Þ
804
L.M. Rodriguez et al. / Measurement 44 (2011) 802–814
Fig. 1. Concentration pattern generated with Gaussian Puff Model over three time steps on a 16 16 receptor grid. The panel on the left shows the concentration with no thresholds and the panel on the right shows the same receptor grid with a saturation and detection level of 1 1012, 1 1016 and detection level.
where Cr is the concentration at receptor r, xr is the coordinate aligned downwind from the emission source, yr is the crosswind coordinate, zr is the vertical coordinate, Q is the emission rate of the source, Dt is duration of the release, t is the time since the release, U is the wind speed, He is the effective height of the puff center, and (rx, ry, rz) are the dispersion coefficients. These latter values are computed according to the fit by McMullen [17] with coefficients tabulated in Beychok [18]. In this study we assume horizontal isotropy so that rx and ry are equivalent. Neutral stability is assumed. We use Dt = 1 s. Note that QDt, the total mass released, is a conserved quantity. We produce the synthetic data over five time steps, each separated by 6 min. We consider seven different grid sizes on a 16 km2 domain (see Table 2) in order to study how large the sensor network Table 2 Characteristics of grids evaluated on a constant 16 km2 domain. Grid size
Number of grid points
Spacing between grid points (km)
22 44 66 88 10 10 12 12 16 16
4 16 36 64 100 144 256
16.00 5.33 3.20 2.29 1.78 1.45 1.07
must be to accurately determine the source characteristics for each set of sensor constraints. In order to model a realistic situation we do not assume that the wind direction is known a priori. Therefore the sensor grid is constructed to surround the source and is never located on the centerline of the puff. As a result, many of the sensors will not be in the path of the contaminant puff. 2.2. Data generation Twin data puffs are created by means of Eq. (1). This is done by first defining the parameters that we wish to backcalculate (wind direction, emission rate, and the x and y values of the source location). Then the ‘‘correct’’ values of these parameters are inserted into Eq. (1) to generate concentration observations at the grid points. The values of those observations are then clipped at specified thresholds to simulate detection and saturation levels. Fig. 2 shows how the concentration strength varies over three time steps for a 16 16 sensor grid before the data is clipped to simulate the saturation and detection levels of the sensors for two different wind directions examined (180° and 225°, i.e. along and diagonal to the grid axes). The minimum detection level is specified with respect to the maximum concentration value and any data under that level is set to zero. Thus, for a 1 1016 cutoff,
L.M. Rodriguez et al. / Measurement 44 (2011) 802–814
805
Fig. 2. Concentration pattern over three time steps on a 16 16 receptor grid. The panel on the left shows the concentration with a 180° wind direction and the panel on the right for a 225° wind direction.
anything smaller than 1 1016 of the maximum concentration is changed to 0.0, and similarly for the 1 1012, 1 108, 1 104, 1 103, 1 102, 1 101, 5 101, 1 100 detection threshold levels. Saturation level is defined similarly for this study. Any value over a specified percentage of the maximum concentration is clipped to that level. For example, for a saturation level of 1 1000 of the maximum concentration, 1.0 kg m3 is used as the cutoff value but for a saturation of 5 101 of the maximum concentration, 0.5 kg m3 is the cutoff. Examples of the detection and saturation cutoff levels are provided in Fig. 3 where the maximum concentration in this illustration is 1 kg m3. After the data is clipped it serves as our ‘‘true observations’’. 2.3. Solution procedure After our ‘‘true observations’’ are generated, we use the Hybrid GA to back-calculate from them the source parameters and the wind direction. The GA begins with a random population of ‘‘first guess’’ values for these variables taken from the ranges shown in Table 3. These trial solutions, known as chromosomes, are used to predict the downwind contaminant concentration via Eq. (1). The results are then compared to our ‘‘true observations’’ by means of the cost function:
ffi P5 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi PTR 2 t¼1 r¼1 ðlog10 ðC r þ eÞ log10 ðRr þ eÞÞ cost function ¼ PTR P5 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 t¼1 r¼1 ðlog10 ðRr þ eÞÞ ð2Þ where Cr is the concentration as predicted by the dispersion model given by (1), Rr is the observation data value at receptor r, TR is the total number of receptors, and e is a constant added to avoid logarithms of zero. Thus, the cost function is formulated to measure the difference (as an L2 norm) of the log of the concentrations predicted with the trial variables from those that are measured by the virtual sensors. The predictions are clipped in the same manner as the ‘‘true observations’’ as described above. GAs work by evaluating an initial population of trial solutions via the cost function then selecting the highest ranking individuals to reproduce, forming a new generation of trial solutions by the GA operators of mating and mutation [19]. These ‘‘offspring’’ solutions are then evaluated in turn and the process iterates (Fig. 4). The real-valued chromosomes used for mating are determined by a rank weighted roulette wheel selection [19]. The mating procedure blends each of the parameter values of the original two chromosomes to create new chromosomes [9], and then a fixed percentage of the population is subject to mutation. Mutation is needed to explore the complete
806
L.M. Rodriguez et al. / Measurement 44 (2011) 802–814
Fig. 3. Data fit to a Gaussian distribution. The maximum concentration normalized to 1 kg m3. The lines indicate the censored levels of 1 1000, 5 1001, 1 1001, 1 1002, 1 1003, 1 1004, 1 1008, 1 1012, and 1 1016.
Table 3 Variable thresholds used to populate the GA. Parameter
Minimum value
Maximum value
Location (x, y) (m) Source strength (kg s1) Wind direction (°)
8000 0 0
8000 20 360
solution space. In mutation, each selected variable is replaced by a random number chosen from a uniform distribution between the values of each parameter from Table 2. The real-valued GA used here is further described in Haupt and Haupt [19]. We use a population of 40 chromosomes, 640 iterations, and a mutation rate of 0.32; these parameters were shown to lead to convergence by Kuroki et al. [20]. Since a GA is initialized with random values, each application of the algorithm will produce a different path toward the global minimum. Therefore, we evaluate the results of an ensemble of 200 optimizations. From this ensemble of 200 optimizations we apply a bootstrapping technique from which we take the median of the ten optimizations 1000 times. Our final result is the mean of the 1000 medians and for error analysis we also evaluate the standard deviation (SDev) of the 1000 medians. This ap-
proach serves to eliminate outliers resulting from the occasional failure of the optimization, and thus, represents a ‘‘typical’’ application. We select 10 as the sample size because in operational use it achieves protection against outlier contamination without the incurring excessive computational cost, as seen in Kuroki et al. [20]. The mean absolute error (MAE) of the ensemble medians is used to evaluate the success of the algorithm in estimating each individual variable. These 200 runs were executed for each of seven grid sizes for each of the 16 detection/saturation combinations and the bootstrapping statistics computed as described above. In this manner, we are able to thoroughly investigate the method’s behavior for a wide range of potential configurations. 3. Results To assess the results an error score is computed for the model for each configuration of saturation and detection levels. Error scores combine the results for the four variables that characterize our source (wind direction, strength, and location x, y) into a single metric that allows a broader representation of the results. They are computed following the method described in Long et al. [13] and summarized in Appendix A. Similar to the MAE used for
Fig. 4. Schematic of genetic algorithm with parameter configuration.
L.M. Rodriguez et al. / Measurement 44 (2011) 802–814
the individual variables, these error scores are constructed so that smaller error scores indicate better results. We consider solutions with error scores less than or equal to 0.1 to be very successful [13]. Figs. 5–7 show the results in terms of the mean and their standard deviation of the error scores for the 1000 ten-run bootstrap samples discussed in Section 2.3. These results indicate that the larger grid sizes generally have a better retrieval quality, but when thresholding the data too severely even those fail. Specifically, in Fig. 5, the vertical contours indicate that the algorithm’s success has little dependence on the minimum detection level. In that case, the number of sensors in the field is the primary constraint as long as there is a two order of magnitude range between the saturation and detection level. On the other hand, if the saturation level is not within an order of magnitude of the maximum concentration value, then there is insufficient information to successfully perform the back-calculation. Thus, the dependence on saturation level is substantially greater than that on the minimum threshold level given that the Hybrid GA is only able to retrieve correct source characteristics when there are at least three orders of magnitude between the maximum concentration and the detection level as shown with Fig. 6. Fig. 7 also shows that even using a
807
range of four orders of magnitude in the lower portion of the clipped Gaussian is not sufficient information to obtain a good solution. This result implies that sensors for such applications should be chosen so that the maximum concentrations are accurately captured if the data are to be used for source characterization. Similar results (not shown) were obtained for other cases with differing wind directions. The mean error score plots in Figs. 5 and 6 contradict the expected monotonic behavior with increasing number of sensors; that is, that the overall skill should improve with grid resolution. For the saturation level at 0.1 of the maximum concentration (Fig. 6), the sole variable that behaved better using a 6 6 grid than an 8 8 grid was the alongwind direction variable, y. The monotonic decrease does occur in Fig. 5 (saturation level of 1.0). However, the SDev is higher in the region of the distorted error score contours, showing more uncertainty in those results. In Fig. 6 we observe vertical error score contours, but there is a local minimum around the 6 6 grid. The SDev is quite low in that region indicating low uncertainty. Upon inspection we determined that the source of this inconsistency was not the clipping of the Gaussian concentrations. In Fig. 6, the mean error score for the 6 6 and 8 8 grids behaves
Fig. 5. Error score and SDev for a constant 1 1000 saturation level.
Fig. 6. Error score and SDev for a constant 1 1001 saturation level.
808
L.M. Rodriguez et al. / Measurement 44 (2011) 802–814
Fig. 7. Error score and SDev for a constant 1 1012 saturation level.
Fig. 8. Error score and SDev for a constant 10 10 grid with the saturation level (ceiling) varying on the abscissa and the detection level (floor) varying on the ordinate.
Fig. 9. Error score and SDev for a constant 4 4 grid with the saturation level (ceiling) varying on the abscissa and the detection level (floor) varying on the ordinate.
inconsistently, yet the SDev of the mean error score indicates certainty in the result, implying that the nonmonotonic behavior is a real phenomenon. To assess why
this is occurring, we considered the puff progression across the domain. We note that the sampling times (every 5 min) with a 5 m/s wind speed results in observations every
L.M. Rodriguez et al. / Measurement 44 (2011) 802–814
1800 m, which is roughly half the grid spacing (3200 m) for the 6 6 grid; therefore, the sampling times for the 6 6 grid resulted in observing the puff every other sampling time, but it observed the same portion of the puff. In contrast, the sampling in the 8 8 grid was at about 80% of the spacing (2286 m) resulting in observing different portions of the puff each sampling time, which caused an incorrect estimate of the along-wind direction variable, y. Given
809
enough observation times, a slower wind speed, a different source location, or any combinations of these would likely change this behavior and the 8 8 grid would eventually correctly identify the along-wind parameter, y. Error scores were also plotted for each grid size with the saturation level varying on the abscissa and the detection level varying on the ordinate (Figs. 8 and 9). Fig. 8 represents the typical scenario where the back-calculation was
Fig. 10. Results for all grid sizes (colored lines) and detection levels (x-axis) of the saturation level at 1 100 of the maximum concentration value. Panel a shows the MAE of h (wind direction, 180°), panel b indicates the MAE of the source strength (1 kg s1), and panels c and d show the MAE of the location (0, 0) (in m) for x and y respectfully. All of the MAEs are computed from the mean of 1000 medians bootstrapped from an ensemble of 200 optimizations. We also evaluate the standard deviation (SDev) of the 1000 medians. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
810
L.M. Rodriguez et al. / Measurement 44 (2011) 802–814
Fig. 11. Results for all grid sizes (colored lines) and detection levels (x-axis) of the saturation level at 1 1012 of the maximum concentration value. Panel a shows the MAE of h (wind direction, 180°), panel b indicates the MAE of the source strength (1 kg s1), and panels c and d show the MAE of the location (0, 0) (in m) for x and y respectfully. All of the MAEs are computed from the mean of 1000 medians bootstrapped from an ensemble of 200 optimizations. We also evaluate the standard deviation (SDev) of the 1000 medians. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
successful for a large grid (10 10) while Fig. 9 illustrates a smaller (4 4) grid where the data quantity becomes insufficient for the algorithm to fully optimize the unknown input variables. These plots show the dependency on grid size,1 but more importantly, they reveal the limita1
Note that an observation is taken at each grid point.
tions of the algorithm when the data is thresholded by the sensor limitations. Fig. 8 documents the skill of every combination of detection level (floor) combined with every possible saturation level (ceiling) that exceeds it for the 10 10 grid. Note that results are quite similar for the grids with even more sensors (not shown). In that figure the Hybrid GA is able to obtain good results when using
L.M. Rodriguez et al. / Measurement 44 (2011) 802–814
a saturation level of 1 1000 and 5 1001 in conjunction with any detection level lower than 1 1002. As the detection level increases, the scenario loses information with the limiting case being that of a binary sensor represented by the diagonal. As the scenario approaches the binary case, we see a breakdown in skill. Fig. 9 plots results for the 4 4 grid, the first resolution at which results begin to degrade. The breakdown in skill for all ceiling and floor combinations is evident, but even in this case, some skill exists for the largest dynamic ranges. Figs. 10 and 11 plot the aforementioned MAE values for each of the solution variables so that we can further analyze the results. Each figure reflects a specific saturation level (full maximum of 1 1001 in Fig. 10 and 1 1012 for Fig. 11) while varying detection levels appear across the
811
abscissa. Varying grid sizes are indicated by line color and marker. In these figures the results for each parameter (wind direction, source strength, and x, y release location) are plotted separately. For the wind direction we are seeking a value of 180°, for the source strength a value of 1 kg/s, and for the source release location of (x, y) = (0, 0) m. The results for the 225° wind direction case are similar, and therefore, not shown here. Fig. 10 indicates that the Hybrid GA retrieves the correct source characteristics (that is, has low MAE) for all the detection levels using the 10 10, 12 12, and 16 16 observation grids at the 1 1000 saturation level but the coarser grids were not as good and often exhibited an inconsistency. The results for saturation levels 5 1001 and 1 1001 (not shown) displayed similar results. When the
Fig. 12. SCIPUFF setup at first time step with source location shown as blue square. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Fig. 13. Error score and SDev after using Hybrid GA for a 10 10 grid with the saturation level (ceiling) varying on the abscissa and the detection level (floor) varying on the ordinate.
812
L.M. Rodriguez et al. / Measurement 44 (2011) 802–814
Table 4 GA and Hybrid GA MAE for a 10 10 grid and a 1.0E00 saturation level. Note that the error values here are not bound. Mean absolute error for a 10 10 grid and a 1.0E00 saturation level Floor
Wind direction
GA solutions 1.0E00 5.0E01 1.0E01 1.0E02 1.0E03 1.0E04 1.0E+08
1.68 0.95 2.03 2.74 13.37 53.57 80.91
Hybrid GA solutions 1.0E00 0.57 5.0E01 3.79 1.0E01 7.91 1.0E02 9.78 1.0E03 180.82 1.0E04 45.10 1.0E+08 142.64
Q
X
Y
Error score
0.89 0.88 0.82 0.76 0.80 0.74 0.73
12.64 5.41 9.37 17.33 35.25 7.02 7.15
1.90 2.02 17.79 1.13 3.27 191.25 209.13
0.08 0.05 0.03 0.03 0.07 0.52 0.57
0.93 1.02 1.01 1.00 1.01 0.71 0.72
3.11 24.14 37.63 59.31 47.23 96.86 152.40
0.62 2.65 28.48 1.82 0.43 84.57 102.67
0.12 0.40 0.39 0.37 0.96 0.64 0.68
saturation level is further reduced to 1 1002, 1 1003, 1 1004, 1 1008, 1 1012 (see Fig. 11), 1 1016 no grid size produces reliable results. These findings confirm that thresholding the data too severely eliminates so much information that retrieval quality declines substantially: thus, greater dynamic range in sensors leads to more accurate inversion for the source characteristics and meteorological variables. This dynamic range produces the greatest accuracy if it extends to the maximum concentration as is illustrated in Fig. 7. 4. Case study To demonstrate the impact of more realistic data limitations and to move away from the identical twin experimental setup, a 5 min scenario was created using the Second-order Closure Integrated Puff model [21]. The dimensions of the grid, quantity of sensors, source location, and emission rate for the scenario created were modeled after the FFT07 field experiment introduced in Section 1.
The grid consists of 10 10 sensors arranged evenly in a 0.5 km 0.5 km domain with the source located at the bottom (southern edge) in the horizontal center of the domain at (0.25, 0). The emission rate used was 0.1 kg/s with a wind speed of 2 m/s, the boundary layer type was a Simple Diurnal with the release at 3PM which produces an unstable stability class. Fig. 12 displays the resulting concentration pattern at the first time step, 1 min after the release. This scenario was used to create our concentration data. It is important to note that the data generated with SCIPuff includes only six orders of magnitude between the maximum and minimum concentrations, unlike the Gaussian generated data that included 300 orders of magnitude between the maximum and minimum, i.e. the internal floating point precision of the programming language used (MATLAB). The tests are limited to sensor floors six orders of magnitude below the peak concentrations. Our Hybrid GA back-calculation model was run on this SCIPuff-created data. The dispersion model for the backcalculation remained the Gaussian Puff model described in Section 2. The GA solution search space was decreased to match the needs of the problem: for the emission rate the search space was reduced to 0–5 kg/s and for the domain to 0–500 m. The range for the wind direction remained the same. A similar ensemble configuration of the Hybrid GA with the Gaussian Puff model was run – bootstrapping and taking the mean of 1000 medians of 10 realizations out of 100 optimizations. Fig. 13 plots the resulting error scores for all possible combinations of detection thresholds (floors) and saturation limits (ceilings) with the maximum concentration normalized to 1. In this more realistic SCIPuff scenario the horizontal error score contours indicate that the saturation level does not influence the back-calculation as long as we can capture the lowest detection level. This finding contrasts to the results shown in the comparable Fig. 8, where for the smallest detection levels the contours were more vertical, indicating sensitivity to saturation level. However, as detection level (floor) increases beyond 1 1002, the detection level becomes more important in agreement with this SCIPuff case (Fig. 13). Therefore, when
Fig. 14. Error score and SDev after GA for a 10 10 grid with the saturation level (ceiling) varying on the abscissa and the detection level (floor) varying on the ordinate.
L.M. Rodriguez et al. / Measurement 44 (2011) 802–814
using more realistic data with a lower dynamic range, the most important feature to capture is the spread of the puff, which is revealed by the full range of concentrations. This finding supports the conclusion above that having access to the largest available dynamic range enhances the ability of our algorithm to estimate the source term. We note, however, that the error score values in Fig. 13 are substantially higher than for the identical twin case of Fig. 8. In fact, upon further examination, we established that the NMDS algorithm produced a worse solution than the GA alone. Table 4 lists the MAEs for each individual variable as well as the error scores for each detection level maintaining a constant saturation level of 1 1000 for both the GA alone and the Hybrid GA that includes the NMDS. For this case study, the GA alone produces much better results (lower MAEs and lower error scores) than when the NMDS is included. The NMDS algorithm does not allow specific bounds on the variables; therefore, we often obtained solutions that were out of our domain. Fig. 14 depicts the error scores when only the GA was used. It is obvious that the detection level becomes less important in agreement with our identical twin results. This analysis indicates that it may be more appropriate to apply only the GA for mixed model applications.
5. Conclusions This study has sought to determine to what extent sensor thresholds limit the feasibility of back-calculating source characteristics and wind direction from sensor field measurements. It has augmented the literature on source term estimation for CBRN data by considering the impact of more realistic sensor characteristics. The Hybrid GA method used here (with NMDS) is successful in back-calculating source characteristics and wind direction with identical twin data that has been thresholded to include detection limits and saturation levels. These thresholds, if too restrictive, eliminate so much information that retrieval quality degrades severely. The identical twin experiments indicate that inversion is most successful if each sensor can detect the maximum concentrations (ie. has a high saturation level), which means that the most effective sensors have this characteristic. The case study using SCIPuff to generate the sensor data demonstrates that one must be careful in applying algorithms to different data types. In fact, the GA alone produced better results than the Hybrid GA that includes the NMDS. For such cases of differing data types, it is critical to use a very robust optimization method that does not rely on a good first guess solution. For sensor data inversion using the Gaussian model with GA optimization both sensor detection and saturation limits are equally important. In contrast, when the SCIPuff model is used with GA optimization the detection limit is more important. Thus, for the Gaussian model sensor dynamic range is the key to success while for SCIPuff both dynamic range and a detection threshold low enough to sample the puff periphery are essential. In conclusion, it is important to select an accurate optimization technique as well as model saturation and detection levels in sensors,
813
specifically when using a GA coupled with an AT&D model to estimate source parameters. This work demonstrates that it is important that the sensors have sufficient dynamic range in order to use the data effectively for source term estimation. From the case study results using SCIPuff, we see that the necessary range varies according to the problem. In addition, a sufficient number of sensors must be available to define the centerline and spread of the plume. Once again, having the full range of concentrations allows the algorithm to better quantify the unknown source variables. However, a robust algorithm (such as a GA) still shows skill when the dynamic range is reduced. This work modeled generic sensor characteristics rather than emulating any particular sensor in order to retain generality. It points out the importance of considering sensor characteristics when developing new numerical techniques that use their data. It also demonstrates that information needs of the intended retrieval mission must be considered when specifying sensors for a monitoring network. Studies such as this one that use known modeled solutions as input data are an important step since this approach allows separation of the sensor characteristics from other sources of error, including sensor calibration, atmospheric turbulence, and inadequacies of the atmospheric transport and dispersion models. Such considerations should be included in future work. Of particular interest is whether higher frequency sampling would improve our ability to compensate for lack of dynamic range. Note that typical sensors sample at a sufficient frequency to detect the puff or plume. To sample more frequently would only further resolve the noise and would therefore not be helpful. The results demonstrated here should be confirmed with concentration data measured during a field experiment and we are currently assessing these issues using the FFT07 data. We note that this current study uses a single method for source characterization to assess the sensor constraints. Other methods should also be considered when evaluating the limits of sensor usage. Note that one can also assimilate the available sensor data directly into atmospheric transport and dispersion models or, equivalently, apply sensor data fusion to those models to improve contaminant concentration predictions without recovering source information [22,23]. Such work is complementary to the current study. In fact, one can view the methods presented here as an example of using sensor data fusion techniques to retrieve the source parameters. In summary, the current study indicates the importance of designing sensor networks for chem-bio applications with knowledge of the sensor constraints and of how those constraints will limit the future use of their data. Additionally, algorithms that use such data must be designed with knowledge of those constraints to prevent failure in a realistic situation. For FFT07 we will use the lessons learned to implement sensor noise thresholds. Acknowledgements This work was supported by the Defense Threat Reduction Agency under Grants W911NF-06-C-0162, 01-03-D-0010-0012, and by the Bunton–Waller Fellowship.
814
L.M. Rodriguez et al. / Measurement 44 (2011) 802–814
Support for professional development was supplied by the Significant Opportunities in Atmospheric Research and Science Program. Dr. Randy L. Haupt provided the framework for the genetic algorithm and Christopher Allen and Kerrie Long contributed to the source characterization code. We would like to thank Kerrie Long, Dave Swanson, Arthur Small, Anke Beyer-Lout, Andrew Annunzio, and Yuki Kuroki for helpful discussions. Finally, we thank two anonymous reviewers for some very helpful comments that helped us refine our analysis and strengthen our results. Appendix A. Definition of error scores The error score is used to assess the accuracy of the solutions found by the GA and the Nelder–Mead search. The ‘‘error’’ is based on the following equations.
Eloc ¼ 1:25 104 dist
ðA:1Þ
Edir ¼ minððjhfou hact jÞ; ð360 jhfou hact jÞÞ=90
ðA:2Þ
Eact Efou Efou Eact Estr ¼ max ; Eact Emin Emax Eact
ðA:3Þ
where dist is the distance between the GA computed source location and the actual source location, (0, 0), in meters, hfou is the wind direction found by the GA, hact is the actual wind direction, 180°, Eact is the actual source strength of 1.0 kg/s, Efou is the strength found by the GA, Emin and Emax are the minimum and maximum values allotted for the GA described in Table 2. Each parameter is normalized on a scale of [0, 1], where 1 is the absolute worst possible solution and 0 is the best. The final error score is the mean of the Eloc, Edir,and Estr. The error scores are evaluated with the results from only the GA and then with results from the Hybrid GA. It is important to note that the Hybrid GA is not bound and the error score may give us a value large than 1, however in this work we bound the error between 0 and 1. References [1] Homeland Security Office, National Strategy for Homeland Security, 2007.. [2] National Research Council, Tracking and Predicting the Atmospheric Dispersion of Hazardous Material Releases. Implications for Homeland Security. The National Academies Press, 2003, 114 pp.
[3] Defence Research and Development Canada, CRTI-IRTC Annual Report 2005–2006, ISSN 1912-0974.. [4] European Union Committee, Preventing Proliferation of Weapons of Mass Destruction: The EU Contribution. Report with Evidence, House of Lords, the Stationery Office Limited, London, 2005. . [5] L.C. Thomson, B. Hirst, G. Gibson, S. Gillespie, P. Jonathan, K.D. Skeldon, M.J. Padgett, An improved algorithm for locating a gas source using inverse methods, Atmospheric Environment 41 (2007) 1128–1134. [6] S. Rao, Source estimation methods for atmospheric dispersion, Atmospheric Environment 41 (2007) 6963–6973. [7] S.E. Haupt, A demonstration of coupled receptor/dispersion modeling with a genetic algorithm, Atmospheric Environment 39 (2005) 7181–7189. [8] S.E. Haupt, G.S. Young, C.T. Allen, Validation of a receptor/dispersion model with a genetic algorithm using synthetic data, Journal of Applied Meteorology 45 (2006) 476–490. [9] C.T. Allen, G.S. Young, S.E. Haupt, Improving pollutant source characterization by optimizing meteorological data with a genetic algorithm, Atmospheric Environment 41 (2006) 2283–2289. [10] C.T. Allen, S.E. Haupt, G.S. Young, Source characterization with a receptor/dispersion model coupled with a genetic algorithm, Journal of Applied Meteorology and Climatology 46 (2007) 273– 287. [11] S.E. Haupt, G.S. Young, C.T. Allen, A genetic algorithm method to assimilate sensor data for a toxic contaminant release, Journal of Computers 2 (2007) 85–93. [12] S.E. Haupt, R.L. Haupt, G.S. Young, A mixed integer genetic algorithm used in chem-bio defense applications, Journal of Soft Computing 15 (2011) 51–59. [13] K.J. Long, S.E. Haupt, G.S. Young, Assessing sensitivity of source term estimation, Atmospheric Environment 44 (2010) 1558–1567. [14] FUsing Sensor Information from Observing Networks (FUSION) Field Trial 2007 (FFT07), 2010.
Contents lists available at ScienceDirect
Measurement journal homepage: www.elsevier.com/locate/measurement
Impact of sensor characteristics on source characterization for dispersion modeling Luna M. Rodriguez a,b,c, Sue Ellen Haupt a,b,c,⇑, George S. Young b a
P.O. Box 30, Applied Research Laboratory, The Pennsylvania State University, State College, PA 16804-0030, United States 503 Walker Building, Meteorology Department, The Pennsylvania State University, University Park, PA 16802, United States c P.O. Box 3000, The National Center for Atmospheric Research/Research Applications Laboratory, Boulder, CO 80307-3000, United States b
a r t i c l e
i n f o
Article history: Received 9 September 2010 Received in revised form 7 January 2011 Accepted 20 January 2011 Available online 2 February 2011 Keywords: Sensor thresholds Sensor saturation Dispersion monitoring Dispersion modeling Source characterization Source term estimation Genetic algorithm
a b s t r a c t An accidental or intentional release of hazardous chemical, biological, radiological, or nuclear material into the atmosphere obligates responsible agencies to model its transport and dispersion in order to mitigate the effects. This modeling requires input parameters that may not be known and must therefore be estimated from sensor measurements of the resulting concentration field. The genetic algorithm (GA) method used here has been successful at back-calculating not only these source characteristics but also the meteorological parameters necessary to predict the contaminants subsequent transport and dispersion. This study assesses the impact of sensor thresholds, i.e. the sensor minimum detection limit and saturation level, on the ability of the algorithm to back-calculate modeling variables. The sensitivity of the back-calculation to these sensor constraints is analyzed in the context of an identical twin approach, where the data is simulated using the same Gaussian Puff model that is used in the back-calculation algorithm in order to analyze sensitivity in a controlled environment. The solution is optimized by the GA and further tuned with the Nelder–Mead downhill simplex algorithm. For this back-calculation to be successful, it is important that the sensor capture the maximum concentrations. Ó 2011 Elsevier Ltd. All rights reserved.
1. Introduction When a toxic agent is released into the atmosphere, government agencies are tasked with tracking and forecasting its transport and dispersion [1–4]. These scenarios include both accidental and intentional release of any chemical, biological, radiological, or nuclear (CBRN) contaminant that can be a threat to human and ecological health. Although tracking the contaminant can be accomplished with arrays of sensors, forecasting the transport and dispersion requires input variables describing the source of the release and the meteorological conditions. Sometimes, however, there is inadequate a priori source or meteorological information to perform those predictions and it becomes necessary to characterize the source ⇑ Corresponding author. Tel.: +1 303 497 2763; fax: +1 303 497 8386. E-mail addresses: [email protected] (L.M. Rodriguez), haupts2@ asme.org (S. Ellen Haupt), [email protected] (G.S. Young). 0263-2241/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.measurement.2011.01.014
of an airborne contaminant from the sensor field measurements of the resulting concentration field. Such source term estimation involves back-calculating the source location and emission rate history. There has been extensive work on methods for back-calculating source characteristics. For example, Thompson et al. [5] use a random search algorithm with simulated annealing and Rao [6] reviews methods that include adjoint and tangent linear models, Kalman filters, and variational data assimilation, among others. Some previous work uses genetic algorithms (GAs) to optimize source characteristics [7–13]. In addition to characterizing the source, some of these studies also back-calculate the meteorological variables required for transport and dispersion modeling, including wind speed, wind direction, and atmospheric stability. Several of the aforementioned papers examine the effect of adding noise to the data to simulate errors in the sensor data as well as the uncertainty inherent in atmospheric turbulence. Most studies, however, do not take into account the sensor
803
L.M. Rodriguez et al. / Measurement 44 (2011) 802–814
constraints, such as the minimum detection threshold and saturation levels. The sensor characteristics that we seek to evaluate are the detection and saturation levels, not only because they exist in real sensors but also to assess how they influence Sensor Data Fusion (SDF) algorithms, specifically our method discussed in Section 3. The understanding of sensor limitations is vital for the military and homeland defense agencies given that sensors used in the field vary in their saturation and detection capabilities. These capabilities can vary from simple binary (on/off) sensors to more complex sensors whose detection and saturation capabilities depend on the gain settings used. For example, the range specifications for the Digital Fast Response Photo-Ionization Detector (digiPID) sensors used in the FUsing Sensor Information from Observing Networks (FUSION) Field Trial (FFT07) are displayed in Table 1. At the lowest gain setting of 0, the sensor saturation is 1000 ppm. At this setting, the response is fast (0.25 ms) but noisy [14]. For higher gain settings, saturation occurs at lower concentration values, but noise is reduced. The threshold level for the propylene released in that experiment is 50 ppb. Due to expense and difficulty of calibration, some networks may include sensors that only supply binary data. For example, the Joint Warning and Reporting Network (JWARN) is a US DoD system that employs binary sensors for early warning and reporting functions [15]. It is precisely such systems that may also be used to back-calculate source variables. If sensors with a high threshold are not able to detect a contaminant it is obvious that many lives may be in danger; however, if the sensors are only able to detect over a certain threshold this again may affect how SDF algorithms interpret the spread of the puff and consequently the source characteristics. Fig. 1 is an example of how a saturation level of one order of magnitude smaller than the maximum concentration detected and a detection level three orders of magnitude smaller than the maximum concentration level detected truncates the concentration field. Much of the information is lost. In real applications SDF algorithms will use the truncated concentration field detected from these sensors. It is important that we understand these highly truncated scenarios even when using sensors with a large dynamic range. In a scenario where concentrations are high, if these sensors become saturated immediately, the concentration field will appear binary to the SDF algorithm. Also, if sensors become saturated the concentration gradient is altered; as a result we lose information about the path and spread of the contaminant puff, which can also affect the retrieval of source characteristics.
Table 1 Sensor specifications for Digital Fast Response Photo-Ionization Detector (digiPID) used in FFT07. Frequency response: 50 Hz. Detection limit: 50 ppb (propylene). Gain setting
Gain factor
Full scale conc. (PPM)
Integration time (MSEC)
Number of integrations/sample
0 1 2 3
1 4 16 32
1000 250 62 31
0.25 1.00 4.00 8.00
64 16 4 2
The goal of the present study is to determine the sensor performance characteristics required to properly back-calculate the source parameters and the wind direction. This is done by first creating ‘‘true data’’ using a Gaussian Puff model and then using a genetic algorithm (GA) coupled with the same Gaussian Puff model to optimize the source and meteorological parameters for agreement with these simulated observations: this is known as the identical twin approach. By applying thresholds to the concentration data, sensor detection and saturation constraints are modeled. In this way the sensitivity of the source characterization to the sensor performance can be evaluated in a controlled and systematic manner. For this study, the source characteristics (x, y) location of the release and source strength) plus the wind direction are back-calculated by a GA, and then fine-tuned to the global optimum using the Nelder–Mead downhill simplex algorithm (NMDS) [16]. The GA coupled with the NMDS will be referred to as a Hybrid GA. The process is repeated for varying numbers of sensors to estimate how many sensors are required to perform the back-calculation for each thresholding scenario. This sensitivity analysis is needed because sensor thresholds limit the amount of information contained in the data. Moreover, because the algorithm works by comparing model predictions to the observations it will perform best if it is aware of these thresholds so that it can model the observed data instead of an ideal distribution. To a larger extent, this study assesses the way that sensor design impacts utility. Specifically, does a low saturation level or high detection level prevent the data from being used for this source characterization? The sensitivity analysis of sensor constraints including data descriptions and procedures appear in Section 2. Section 3 presents the results. Section 4 is a case study using a more refined model to create the test data. Section 5 summaries the conclusions and discusses the implications of the study. 2. Model construction We construct an identical twin experiment with Gaussian Puff model output being used to back-calculate source strength, location of the release, and wind direction. We first describe how the concentration data is generated, then discuss the application of thresholds to the data to simulate the sensor constraints, and finally describe the Hybrid GA used to do the back-calculation. 2.1. The Gaussian Puff model The Gaussian Puff model is used to produce the synthetic concentration observations:
Cr ¼
Q Dt
exp
ðxr UtÞ2 2r2x
!
ð2pÞ1:5 rx ry rz !" ! y2r ðzr He Þ2 exp exp 2r2y 2r2z !# ðzr þ He Þ2 þ exp 2r2x
ð1Þ
804
L.M. Rodriguez et al. / Measurement 44 (2011) 802–814
Fig. 1. Concentration pattern generated with Gaussian Puff Model over three time steps on a 16 16 receptor grid. The panel on the left shows the concentration with no thresholds and the panel on the right shows the same receptor grid with a saturation and detection level of 1 1012, 1 1016 and detection level.
where Cr is the concentration at receptor r, xr is the coordinate aligned downwind from the emission source, yr is the crosswind coordinate, zr is the vertical coordinate, Q is the emission rate of the source, Dt is duration of the release, t is the time since the release, U is the wind speed, He is the effective height of the puff center, and (rx, ry, rz) are the dispersion coefficients. These latter values are computed according to the fit by McMullen [17] with coefficients tabulated in Beychok [18]. In this study we assume horizontal isotropy so that rx and ry are equivalent. Neutral stability is assumed. We use Dt = 1 s. Note that QDt, the total mass released, is a conserved quantity. We produce the synthetic data over five time steps, each separated by 6 min. We consider seven different grid sizes on a 16 km2 domain (see Table 2) in order to study how large the sensor network Table 2 Characteristics of grids evaluated on a constant 16 km2 domain. Grid size
Number of grid points
Spacing between grid points (km)
22 44 66 88 10 10 12 12 16 16
4 16 36 64 100 144 256
16.00 5.33 3.20 2.29 1.78 1.45 1.07
must be to accurately determine the source characteristics for each set of sensor constraints. In order to model a realistic situation we do not assume that the wind direction is known a priori. Therefore the sensor grid is constructed to surround the source and is never located on the centerline of the puff. As a result, many of the sensors will not be in the path of the contaminant puff. 2.2. Data generation Twin data puffs are created by means of Eq. (1). This is done by first defining the parameters that we wish to backcalculate (wind direction, emission rate, and the x and y values of the source location). Then the ‘‘correct’’ values of these parameters are inserted into Eq. (1) to generate concentration observations at the grid points. The values of those observations are then clipped at specified thresholds to simulate detection and saturation levels. Fig. 2 shows how the concentration strength varies over three time steps for a 16 16 sensor grid before the data is clipped to simulate the saturation and detection levels of the sensors for two different wind directions examined (180° and 225°, i.e. along and diagonal to the grid axes). The minimum detection level is specified with respect to the maximum concentration value and any data under that level is set to zero. Thus, for a 1 1016 cutoff,
L.M. Rodriguez et al. / Measurement 44 (2011) 802–814
805
Fig. 2. Concentration pattern over three time steps on a 16 16 receptor grid. The panel on the left shows the concentration with a 180° wind direction and the panel on the right for a 225° wind direction.
anything smaller than 1 1016 of the maximum concentration is changed to 0.0, and similarly for the 1 1012, 1 108, 1 104, 1 103, 1 102, 1 101, 5 101, 1 100 detection threshold levels. Saturation level is defined similarly for this study. Any value over a specified percentage of the maximum concentration is clipped to that level. For example, for a saturation level of 1 1000 of the maximum concentration, 1.0 kg m3 is used as the cutoff value but for a saturation of 5 101 of the maximum concentration, 0.5 kg m3 is the cutoff. Examples of the detection and saturation cutoff levels are provided in Fig. 3 where the maximum concentration in this illustration is 1 kg m3. After the data is clipped it serves as our ‘‘true observations’’. 2.3. Solution procedure After our ‘‘true observations’’ are generated, we use the Hybrid GA to back-calculate from them the source parameters and the wind direction. The GA begins with a random population of ‘‘first guess’’ values for these variables taken from the ranges shown in Table 3. These trial solutions, known as chromosomes, are used to predict the downwind contaminant concentration via Eq. (1). The results are then compared to our ‘‘true observations’’ by means of the cost function:
ffi P5 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi PTR 2 t¼1 r¼1 ðlog10 ðC r þ eÞ log10 ðRr þ eÞÞ cost function ¼ PTR P5 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 t¼1 r¼1 ðlog10 ðRr þ eÞÞ ð2Þ where Cr is the concentration as predicted by the dispersion model given by (1), Rr is the observation data value at receptor r, TR is the total number of receptors, and e is a constant added to avoid logarithms of zero. Thus, the cost function is formulated to measure the difference (as an L2 norm) of the log of the concentrations predicted with the trial variables from those that are measured by the virtual sensors. The predictions are clipped in the same manner as the ‘‘true observations’’ as described above. GAs work by evaluating an initial population of trial solutions via the cost function then selecting the highest ranking individuals to reproduce, forming a new generation of trial solutions by the GA operators of mating and mutation [19]. These ‘‘offspring’’ solutions are then evaluated in turn and the process iterates (Fig. 4). The real-valued chromosomes used for mating are determined by a rank weighted roulette wheel selection [19]. The mating procedure blends each of the parameter values of the original two chromosomes to create new chromosomes [9], and then a fixed percentage of the population is subject to mutation. Mutation is needed to explore the complete
806
L.M. Rodriguez et al. / Measurement 44 (2011) 802–814
Fig. 3. Data fit to a Gaussian distribution. The maximum concentration normalized to 1 kg m3. The lines indicate the censored levels of 1 1000, 5 1001, 1 1001, 1 1002, 1 1003, 1 1004, 1 1008, 1 1012, and 1 1016.
Table 3 Variable thresholds used to populate the GA. Parameter
Minimum value
Maximum value
Location (x, y) (m) Source strength (kg s1) Wind direction (°)
8000 0 0
8000 20 360
solution space. In mutation, each selected variable is replaced by a random number chosen from a uniform distribution between the values of each parameter from Table 2. The real-valued GA used here is further described in Haupt and Haupt [19]. We use a population of 40 chromosomes, 640 iterations, and a mutation rate of 0.32; these parameters were shown to lead to convergence by Kuroki et al. [20]. Since a GA is initialized with random values, each application of the algorithm will produce a different path toward the global minimum. Therefore, we evaluate the results of an ensemble of 200 optimizations. From this ensemble of 200 optimizations we apply a bootstrapping technique from which we take the median of the ten optimizations 1000 times. Our final result is the mean of the 1000 medians and for error analysis we also evaluate the standard deviation (SDev) of the 1000 medians. This ap-
proach serves to eliminate outliers resulting from the occasional failure of the optimization, and thus, represents a ‘‘typical’’ application. We select 10 as the sample size because in operational use it achieves protection against outlier contamination without the incurring excessive computational cost, as seen in Kuroki et al. [20]. The mean absolute error (MAE) of the ensemble medians is used to evaluate the success of the algorithm in estimating each individual variable. These 200 runs were executed for each of seven grid sizes for each of the 16 detection/saturation combinations and the bootstrapping statistics computed as described above. In this manner, we are able to thoroughly investigate the method’s behavior for a wide range of potential configurations. 3. Results To assess the results an error score is computed for the model for each configuration of saturation and detection levels. Error scores combine the results for the four variables that characterize our source (wind direction, strength, and location x, y) into a single metric that allows a broader representation of the results. They are computed following the method described in Long et al. [13] and summarized in Appendix A. Similar to the MAE used for
Fig. 4. Schematic of genetic algorithm with parameter configuration.
L.M. Rodriguez et al. / Measurement 44 (2011) 802–814
the individual variables, these error scores are constructed so that smaller error scores indicate better results. We consider solutions with error scores less than or equal to 0.1 to be very successful [13]. Figs. 5–7 show the results in terms of the mean and their standard deviation of the error scores for the 1000 ten-run bootstrap samples discussed in Section 2.3. These results indicate that the larger grid sizes generally have a better retrieval quality, but when thresholding the data too severely even those fail. Specifically, in Fig. 5, the vertical contours indicate that the algorithm’s success has little dependence on the minimum detection level. In that case, the number of sensors in the field is the primary constraint as long as there is a two order of magnitude range between the saturation and detection level. On the other hand, if the saturation level is not within an order of magnitude of the maximum concentration value, then there is insufficient information to successfully perform the back-calculation. Thus, the dependence on saturation level is substantially greater than that on the minimum threshold level given that the Hybrid GA is only able to retrieve correct source characteristics when there are at least three orders of magnitude between the maximum concentration and the detection level as shown with Fig. 6. Fig. 7 also shows that even using a
807
range of four orders of magnitude in the lower portion of the clipped Gaussian is not sufficient information to obtain a good solution. This result implies that sensors for such applications should be chosen so that the maximum concentrations are accurately captured if the data are to be used for source characterization. Similar results (not shown) were obtained for other cases with differing wind directions. The mean error score plots in Figs. 5 and 6 contradict the expected monotonic behavior with increasing number of sensors; that is, that the overall skill should improve with grid resolution. For the saturation level at 0.1 of the maximum concentration (Fig. 6), the sole variable that behaved better using a 6 6 grid than an 8 8 grid was the alongwind direction variable, y. The monotonic decrease does occur in Fig. 5 (saturation level of 1.0). However, the SDev is higher in the region of the distorted error score contours, showing more uncertainty in those results. In Fig. 6 we observe vertical error score contours, but there is a local minimum around the 6 6 grid. The SDev is quite low in that region indicating low uncertainty. Upon inspection we determined that the source of this inconsistency was not the clipping of the Gaussian concentrations. In Fig. 6, the mean error score for the 6 6 and 8 8 grids behaves
Fig. 5. Error score and SDev for a constant 1 1000 saturation level.
Fig. 6. Error score and SDev for a constant 1 1001 saturation level.
808
L.M. Rodriguez et al. / Measurement 44 (2011) 802–814
Fig. 7. Error score and SDev for a constant 1 1012 saturation level.
Fig. 8. Error score and SDev for a constant 10 10 grid with the saturation level (ceiling) varying on the abscissa and the detection level (floor) varying on the ordinate.
Fig. 9. Error score and SDev for a constant 4 4 grid with the saturation level (ceiling) varying on the abscissa and the detection level (floor) varying on the ordinate.
inconsistently, yet the SDev of the mean error score indicates certainty in the result, implying that the nonmonotonic behavior is a real phenomenon. To assess why
this is occurring, we considered the puff progression across the domain. We note that the sampling times (every 5 min) with a 5 m/s wind speed results in observations every
L.M. Rodriguez et al. / Measurement 44 (2011) 802–814
1800 m, which is roughly half the grid spacing (3200 m) for the 6 6 grid; therefore, the sampling times for the 6 6 grid resulted in observing the puff every other sampling time, but it observed the same portion of the puff. In contrast, the sampling in the 8 8 grid was at about 80% of the spacing (2286 m) resulting in observing different portions of the puff each sampling time, which caused an incorrect estimate of the along-wind direction variable, y. Given
809
enough observation times, a slower wind speed, a different source location, or any combinations of these would likely change this behavior and the 8 8 grid would eventually correctly identify the along-wind parameter, y. Error scores were also plotted for each grid size with the saturation level varying on the abscissa and the detection level varying on the ordinate (Figs. 8 and 9). Fig. 8 represents the typical scenario where the back-calculation was
Fig. 10. Results for all grid sizes (colored lines) and detection levels (x-axis) of the saturation level at 1 100 of the maximum concentration value. Panel a shows the MAE of h (wind direction, 180°), panel b indicates the MAE of the source strength (1 kg s1), and panels c and d show the MAE of the location (0, 0) (in m) for x and y respectfully. All of the MAEs are computed from the mean of 1000 medians bootstrapped from an ensemble of 200 optimizations. We also evaluate the standard deviation (SDev) of the 1000 medians. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
810
L.M. Rodriguez et al. / Measurement 44 (2011) 802–814
Fig. 11. Results for all grid sizes (colored lines) and detection levels (x-axis) of the saturation level at 1 1012 of the maximum concentration value. Panel a shows the MAE of h (wind direction, 180°), panel b indicates the MAE of the source strength (1 kg s1), and panels c and d show the MAE of the location (0, 0) (in m) for x and y respectfully. All of the MAEs are computed from the mean of 1000 medians bootstrapped from an ensemble of 200 optimizations. We also evaluate the standard deviation (SDev) of the 1000 medians. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
successful for a large grid (10 10) while Fig. 9 illustrates a smaller (4 4) grid where the data quantity becomes insufficient for the algorithm to fully optimize the unknown input variables. These plots show the dependency on grid size,1 but more importantly, they reveal the limita1
Note that an observation is taken at each grid point.
tions of the algorithm when the data is thresholded by the sensor limitations. Fig. 8 documents the skill of every combination of detection level (floor) combined with every possible saturation level (ceiling) that exceeds it for the 10 10 grid. Note that results are quite similar for the grids with even more sensors (not shown). In that figure the Hybrid GA is able to obtain good results when using
L.M. Rodriguez et al. / Measurement 44 (2011) 802–814
a saturation level of 1 1000 and 5 1001 in conjunction with any detection level lower than 1 1002. As the detection level increases, the scenario loses information with the limiting case being that of a binary sensor represented by the diagonal. As the scenario approaches the binary case, we see a breakdown in skill. Fig. 9 plots results for the 4 4 grid, the first resolution at which results begin to degrade. The breakdown in skill for all ceiling and floor combinations is evident, but even in this case, some skill exists for the largest dynamic ranges. Figs. 10 and 11 plot the aforementioned MAE values for each of the solution variables so that we can further analyze the results. Each figure reflects a specific saturation level (full maximum of 1 1001 in Fig. 10 and 1 1012 for Fig. 11) while varying detection levels appear across the
811
abscissa. Varying grid sizes are indicated by line color and marker. In these figures the results for each parameter (wind direction, source strength, and x, y release location) are plotted separately. For the wind direction we are seeking a value of 180°, for the source strength a value of 1 kg/s, and for the source release location of (x, y) = (0, 0) m. The results for the 225° wind direction case are similar, and therefore, not shown here. Fig. 10 indicates that the Hybrid GA retrieves the correct source characteristics (that is, has low MAE) for all the detection levels using the 10 10, 12 12, and 16 16 observation grids at the 1 1000 saturation level but the coarser grids were not as good and often exhibited an inconsistency. The results for saturation levels 5 1001 and 1 1001 (not shown) displayed similar results. When the
Fig. 12. SCIPUFF setup at first time step with source location shown as blue square. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Fig. 13. Error score and SDev after using Hybrid GA for a 10 10 grid with the saturation level (ceiling) varying on the abscissa and the detection level (floor) varying on the ordinate.
812
L.M. Rodriguez et al. / Measurement 44 (2011) 802–814
Table 4 GA and Hybrid GA MAE for a 10 10 grid and a 1.0E00 saturation level. Note that the error values here are not bound. Mean absolute error for a 10 10 grid and a 1.0E00 saturation level Floor
Wind direction
GA solutions 1.0E00 5.0E01 1.0E01 1.0E02 1.0E03 1.0E04 1.0E+08
1.68 0.95 2.03 2.74 13.37 53.57 80.91
Hybrid GA solutions 1.0E00 0.57 5.0E01 3.79 1.0E01 7.91 1.0E02 9.78 1.0E03 180.82 1.0E04 45.10 1.0E+08 142.64
Q
X
Y
Error score
0.89 0.88 0.82 0.76 0.80 0.74 0.73
12.64 5.41 9.37 17.33 35.25 7.02 7.15
1.90 2.02 17.79 1.13 3.27 191.25 209.13
0.08 0.05 0.03 0.03 0.07 0.52 0.57
0.93 1.02 1.01 1.00 1.01 0.71 0.72
3.11 24.14 37.63 59.31 47.23 96.86 152.40
0.62 2.65 28.48 1.82 0.43 84.57 102.67
0.12 0.40 0.39 0.37 0.96 0.64 0.68
saturation level is further reduced to 1 1002, 1 1003, 1 1004, 1 1008, 1 1012 (see Fig. 11), 1 1016 no grid size produces reliable results. These findings confirm that thresholding the data too severely eliminates so much information that retrieval quality declines substantially: thus, greater dynamic range in sensors leads to more accurate inversion for the source characteristics and meteorological variables. This dynamic range produces the greatest accuracy if it extends to the maximum concentration as is illustrated in Fig. 7. 4. Case study To demonstrate the impact of more realistic data limitations and to move away from the identical twin experimental setup, a 5 min scenario was created using the Second-order Closure Integrated Puff model [21]. The dimensions of the grid, quantity of sensors, source location, and emission rate for the scenario created were modeled after the FFT07 field experiment introduced in Section 1.
The grid consists of 10 10 sensors arranged evenly in a 0.5 km 0.5 km domain with the source located at the bottom (southern edge) in the horizontal center of the domain at (0.25, 0). The emission rate used was 0.1 kg/s with a wind speed of 2 m/s, the boundary layer type was a Simple Diurnal with the release at 3PM which produces an unstable stability class. Fig. 12 displays the resulting concentration pattern at the first time step, 1 min after the release. This scenario was used to create our concentration data. It is important to note that the data generated with SCIPuff includes only six orders of magnitude between the maximum and minimum concentrations, unlike the Gaussian generated data that included 300 orders of magnitude between the maximum and minimum, i.e. the internal floating point precision of the programming language used (MATLAB). The tests are limited to sensor floors six orders of magnitude below the peak concentrations. Our Hybrid GA back-calculation model was run on this SCIPuff-created data. The dispersion model for the backcalculation remained the Gaussian Puff model described in Section 2. The GA solution search space was decreased to match the needs of the problem: for the emission rate the search space was reduced to 0–5 kg/s and for the domain to 0–500 m. The range for the wind direction remained the same. A similar ensemble configuration of the Hybrid GA with the Gaussian Puff model was run – bootstrapping and taking the mean of 1000 medians of 10 realizations out of 100 optimizations. Fig. 13 plots the resulting error scores for all possible combinations of detection thresholds (floors) and saturation limits (ceilings) with the maximum concentration normalized to 1. In this more realistic SCIPuff scenario the horizontal error score contours indicate that the saturation level does not influence the back-calculation as long as we can capture the lowest detection level. This finding contrasts to the results shown in the comparable Fig. 8, where for the smallest detection levels the contours were more vertical, indicating sensitivity to saturation level. However, as detection level (floor) increases beyond 1 1002, the detection level becomes more important in agreement with this SCIPuff case (Fig. 13). Therefore, when
Fig. 14. Error score and SDev after GA for a 10 10 grid with the saturation level (ceiling) varying on the abscissa and the detection level (floor) varying on the ordinate.
L.M. Rodriguez et al. / Measurement 44 (2011) 802–814
using more realistic data with a lower dynamic range, the most important feature to capture is the spread of the puff, which is revealed by the full range of concentrations. This finding supports the conclusion above that having access to the largest available dynamic range enhances the ability of our algorithm to estimate the source term. We note, however, that the error score values in Fig. 13 are substantially higher than for the identical twin case of Fig. 8. In fact, upon further examination, we established that the NMDS algorithm produced a worse solution than the GA alone. Table 4 lists the MAEs for each individual variable as well as the error scores for each detection level maintaining a constant saturation level of 1 1000 for both the GA alone and the Hybrid GA that includes the NMDS. For this case study, the GA alone produces much better results (lower MAEs and lower error scores) than when the NMDS is included. The NMDS algorithm does not allow specific bounds on the variables; therefore, we often obtained solutions that were out of our domain. Fig. 14 depicts the error scores when only the GA was used. It is obvious that the detection level becomes less important in agreement with our identical twin results. This analysis indicates that it may be more appropriate to apply only the GA for mixed model applications.
5. Conclusions This study has sought to determine to what extent sensor thresholds limit the feasibility of back-calculating source characteristics and wind direction from sensor field measurements. It has augmented the literature on source term estimation for CBRN data by considering the impact of more realistic sensor characteristics. The Hybrid GA method used here (with NMDS) is successful in back-calculating source characteristics and wind direction with identical twin data that has been thresholded to include detection limits and saturation levels. These thresholds, if too restrictive, eliminate so much information that retrieval quality degrades severely. The identical twin experiments indicate that inversion is most successful if each sensor can detect the maximum concentrations (ie. has a high saturation level), which means that the most effective sensors have this characteristic. The case study using SCIPuff to generate the sensor data demonstrates that one must be careful in applying algorithms to different data types. In fact, the GA alone produced better results than the Hybrid GA that includes the NMDS. For such cases of differing data types, it is critical to use a very robust optimization method that does not rely on a good first guess solution. For sensor data inversion using the Gaussian model with GA optimization both sensor detection and saturation limits are equally important. In contrast, when the SCIPuff model is used with GA optimization the detection limit is more important. Thus, for the Gaussian model sensor dynamic range is the key to success while for SCIPuff both dynamic range and a detection threshold low enough to sample the puff periphery are essential. In conclusion, it is important to select an accurate optimization technique as well as model saturation and detection levels in sensors,
813
specifically when using a GA coupled with an AT&D model to estimate source parameters. This work demonstrates that it is important that the sensors have sufficient dynamic range in order to use the data effectively for source term estimation. From the case study results using SCIPuff, we see that the necessary range varies according to the problem. In addition, a sufficient number of sensors must be available to define the centerline and spread of the plume. Once again, having the full range of concentrations allows the algorithm to better quantify the unknown source variables. However, a robust algorithm (such as a GA) still shows skill when the dynamic range is reduced. This work modeled generic sensor characteristics rather than emulating any particular sensor in order to retain generality. It points out the importance of considering sensor characteristics when developing new numerical techniques that use their data. It also demonstrates that information needs of the intended retrieval mission must be considered when specifying sensors for a monitoring network. Studies such as this one that use known modeled solutions as input data are an important step since this approach allows separation of the sensor characteristics from other sources of error, including sensor calibration, atmospheric turbulence, and inadequacies of the atmospheric transport and dispersion models. Such considerations should be included in future work. Of particular interest is whether higher frequency sampling would improve our ability to compensate for lack of dynamic range. Note that typical sensors sample at a sufficient frequency to detect the puff or plume. To sample more frequently would only further resolve the noise and would therefore not be helpful. The results demonstrated here should be confirmed with concentration data measured during a field experiment and we are currently assessing these issues using the FFT07 data. We note that this current study uses a single method for source characterization to assess the sensor constraints. Other methods should also be considered when evaluating the limits of sensor usage. Note that one can also assimilate the available sensor data directly into atmospheric transport and dispersion models or, equivalently, apply sensor data fusion to those models to improve contaminant concentration predictions without recovering source information [22,23]. Such work is complementary to the current study. In fact, one can view the methods presented here as an example of using sensor data fusion techniques to retrieve the source parameters. In summary, the current study indicates the importance of designing sensor networks for chem-bio applications with knowledge of the sensor constraints and of how those constraints will limit the future use of their data. Additionally, algorithms that use such data must be designed with knowledge of those constraints to prevent failure in a realistic situation. For FFT07 we will use the lessons learned to implement sensor noise thresholds. Acknowledgements This work was supported by the Defense Threat Reduction Agency under Grants W911NF-06-C-0162, 01-03-D-0010-0012, and by the Bunton–Waller Fellowship.
814
L.M. Rodriguez et al. / Measurement 44 (2011) 802–814
Support for professional development was supplied by the Significant Opportunities in Atmospheric Research and Science Program. Dr. Randy L. Haupt provided the framework for the genetic algorithm and Christopher Allen and Kerrie Long contributed to the source characterization code. We would like to thank Kerrie Long, Dave Swanson, Arthur Small, Anke Beyer-Lout, Andrew Annunzio, and Yuki Kuroki for helpful discussions. Finally, we thank two anonymous reviewers for some very helpful comments that helped us refine our analysis and strengthen our results. Appendix A. Definition of error scores The error score is used to assess the accuracy of the solutions found by the GA and the Nelder–Mead search. The ‘‘error’’ is based on the following equations.
Eloc ¼ 1:25 104 dist
ðA:1Þ
Edir ¼ minððjhfou hact jÞ; ð360 jhfou hact jÞÞ=90
ðA:2Þ
Eact Efou Efou Eact Estr ¼ max ; Eact Emin Emax Eact
ðA:3Þ
where dist is the distance between the GA computed source location and the actual source location, (0, 0), in meters, hfou is the wind direction found by the GA, hact is the actual wind direction, 180°, Eact is the actual source strength of 1.0 kg/s, Efou is the strength found by the GA, Emin and Emax are the minimum and maximum values allotted for the GA described in Table 2. Each parameter is normalized on a scale of [0, 1], where 1 is the absolute worst possible solution and 0 is the best. The final error score is the mean of the Eloc, Edir,and Estr. The error scores are evaluated with the results from only the GA and then with results from the Hybrid GA. It is important to note that the Hybrid GA is not bound and the error score may give us a value large than 1, however in this work we bound the error between 0 and 1. References [1] Homeland Security Office, National Strategy for Homeland Security, 2007.
[3] Defence Research and Development Canada, CRTI-IRTC Annual Report 2005–2006, ISSN 1912-0974.