Reconstructive Mapping from Sparsely-Sampled Groundwater Data Using Compressive Sensing

Journal of Geographic Information System, 2021, 13, 287-301 https://www.scirp.org/journal/jgis ISSN Online: 2151-1969 ISSN Print: 2151-1950 DOI: 10.4236/jgis.2021.133016 May 10, 2021 287 Journal of Geographic Information System Reconstructive Mapping from Sparsely-Sampled Groundwater Data Using Compressive Sensing T.-W. Lee1*, J. Y. Lee1, J. E. Park1, H. Bellerova2, M. Raudensky2 1Mechanical and Aerospace Engineering, Arizona State University, Tempe, AZ, USA 2Heat Laboratory, Brno University of Technology, Brno, Czech Republic Abstract Compressive sensing is a powerful method for reconstruction of sparsely-sampled data, based on statistical optimization. It can be applied t o a range of flow measurement and visualization data, and in this work we show the usage in groundwater mapping. Due to scarcity of water in many regions of the world, including southwestern United States, monitoring and management of groundwater is of utm ost importance. A complete mapping of groundwater is difficult since the monitored sites are far from one another, and thus the data sets are considered extremely “sparse”. To overcome this difficulty in complete mapping of groundwater, compressive sensing is an ideal tool, as it bypasses the classical Nyquist criterion. We show that compressive sensing can effectively be used for reconstructions of gro undwater level maps, by validating against data. This approach can have an impact on geographical sensing and information, as effective monitoring and management are enabled without constructing numerous or expensive measu rement sites for groundwater. Keywords Visualization Data, Compressive Sensing, Reconstruction, Mapping 1. Introduction Fresh water is critical to sustainable living, economic growth, social stability and public health [1]. For management and dispensing of the groundwater, accurate monitoring is of utmost importance. However, existing methods such as drilling are expensive and intrusive, and alternate methods such as gravimetry and seisHow to cite this paper: Lee, T.-W., Lee, J.Y ., Park, J.E., Bellerova, H. and Raudensky , M. (2021) Reconstructive Mapping from SparselySampled Groundwater Data Using Compressive Sensing . Journal of Ge ographic Information System , 13 , 287-301. https://doi.org/10.4236/jgis.2021.133016 Received: November 22, 2020 Accepted: May 7, 2021 Published: May 10, 2021 Copyright © 20 21 by author(s) and Scientific Research Publishing Inc. This work is licensed under the Creative Commons Attribution International License (CC BY 4.0). http://creativecommons.org/licenses/by/4.0/ Open Access T.-W. Lee et al. DOI: 10.4236/jgis.2021.133016 288 Journal of Geographic Information System mic sounding are also costly to perform, subject to errors and limited in range. Thus, a method for mapping complete distributions of groundwater based on sparse data sets would be a tremendous asset in groundwater management. Here, we introduce the use of compressed sensing (CS) for generating accurate contour maps of groundwater levels based on sparsely-sampled or incomplete data. Compressive sensing has been developed for signal/image processing technology to enable reconstruction of the source signal based on sparsely sampled data, far below the Nyquist criterion. We have applied this method for recovering turbulent flow structures in so-called “sub-grid” scales [2]. Here, we demonstrate that the compressive sensing can also be used for generating accurate groundwater level maps from sparsely sampled data, by using direct measurement data in Arizona and in other parts of the southwestern United States (California, Nevada, Utah and Colorado). The method is compared with conventional interpolation method, such as kriging, which demonstrates that compressive sensing outperforms by a large margin. Compressive sensing is a powerful method that allows for the actual signals to be recovered from far fewer samples than what have been possible according to the Nyquist-Shannon sampling theorem [3]. The Nyquist-Shannon sampling theorem states that if a function f ( t ) has no frequencies higher than υ Hz, then tracing its ordinates at a series of points (spaced at least by 1/2 υ seconds apart) determines f ( t ) perfectly. This requirement had limited the full mapping of data, unless very tightly-spaced data were sampled, either in space or time. Compressive sensing, in this regard, is a transcending method, initially developed in image processing, where a very small amount of data can be effectively used to reconstruct the full mapping [3]. Compressive sensing works in a way similar (but opposite) to image compression methods such as the JPEG protocol. A signal f ( x ), in general, can be expressed in a series of orthonormal functions, Ψ i (e.g. Fourier series of sinusoidal functions). () ( ) 1 n ii i fx c x ψ = = ∑ (1) For most signals, only a small number of the coefficients ci are of significant magnitude, and most of the others can be discarded with negligible loss of information. Image compression works by having the full data, f ( x ), and evaluating the coefficients, ci , a priori . Since the full image data are used during image compression, the entirety of f ( x ) is available for de-composition (Equation (1)) to find ci , through Fourier transform, for example. From the set of coefficients, only the most significant ci terms are retained, for later reconstruction, which is how image compression/reconstruction works. Compressive sensing takes advantage of this concept by attempting to find this coefficient set, ci , in the absence of full data. Therefore, an algorithm is required to “recover” ci , from under-sampled signal, f ( x ) in Equation (1), with the requirement that the coefficient matrix must have certain properties such as a small number of dominant terms. Due to under-sampled f ( x ), the equation or the system of equations such T.-W. Lee et al. DOI: 10.4236/jgis.2021.133016 289 Journal of Geographic Information System as Equation (1) is highly indeterminate with a constraint (a condition that the matrix must follow). For this type of problems, highly effective numerical methods such as l 1-optimization, gradient pursuit or other convex optimization algorithms have been developed [4]-[12]. Although Equation (1) is written in one-dimensional form for illustration purposes, the method is easily translatable to two-dimensions, albeit at larger matrix size and computational time. Here, we demonstrate that compressive sensing can be applied for constructing groundwater mapping with far fewer sampled data than classically thought possible (Nyquist criterion). Compressive sensing has now become quite established in many applications of image analyses and other areas [2] [4]-[12]. Compressive sensing is of course not without limits, and does need the following criteria to work, in applications such as groundwater mapping: 1) The desired groundwater network database is “compressible”, meaning that the discrete data (at sparsely-spaced probe locations) can be approximated as a series of orthogonal functions (e.g. sinusoidal function in a Fourier series). Based on working with most groundwater network data base, we find that this condition is satisfied since most continuous signals can be expressed reasonably accurately using truncated Fourier series. 2) The second condition relates to the coefficient matrix, ci in Equation (1) as an example. This matrix needs to have a relatively small number of dominant terms, so that the indeterminate matrix can be numerically found under the “compressibility” constraint. Again, tests with groundwater database suggest that this condition is also satisfied. Technical details can be found in the Appendix and in the references cited above. If the above criteria are met, then numerical algorithms such as the l 1-optimization or total variational optimization routines can be used to find highly under-determined data matrix based on sparse input data. The step-by-step computational procedures used in this work follow Boyd and Vandenberghe [13] and Candes et al. [3]. 2. Compressive Sensing for Reconstructive Groundwater Mapping For a demonstration of the compressive sensing technique, a sample MODFLOW image (256 by 256 pixel resolution or a total pixel count of 65,536) is used, for groundwater level in Fort Cobb Reservoir experimental watershed (780 km2) [14], as shown in Figure 1(a). The location is in southwestern Oklahoma, USA, and the region shown in Figure 1(a) shows a mapping for a portion of the Rush Spring aquifer [14]. The feed streams supply the groundwater, from which diffusion patterns are observable. First, we start with sampling of only 2025 pixel values (out of 65,536 total) V-shaped (pyramidal) sampling lines (top box in Figure 1(b), which is an extremely sparse pixel sampling of 3%. Even at such low sampling density, remarkably some baseline features of the groundwater distribution in the MODFLOW data are recovered, as shown in Figure 1(c). T.-W. Lee et al. DOI: 10.4236/jgis.2021.133016 290 Journal of Geographic Information System Figure 1. (a) Original MODFLOW contour map (256 by 256 pixels) of groundwater [14]; (b) Pyramidal sampling lines at various sampling rates (0.0308, 0.0761, 0.1492, and 0.2857, top to bottom); (c) Reconstructed images. The sampling density is progressively increased from top to bottom, in Figure 1(b), which results in rapid increase in the fidelity of the reconstructed data (Figure 1(c)). Visually, the original MODFLOW data are nearly fully recovered at sampling density of only 14%, nearly an order of magnitude less than the full data. Thus, Figure 1 demonstrates the utility of compressive sensing where sparsely sampled data can be used to map out the two-dimensional data in spite of the complex feed and distribution patterns for groundwater. Of course, practical monitoring of groundwater does not follow any specific geometrical patterns; nonetheless the compressive sensing operations work in non-geometrical sampling as long as above-mentioned criteria are met [3]. We can track the rate at which the original data can be recovered at different sampling rates, and also test different data sampling geometries (pyramidal, radial and Cartesian), for the same MODFLOW data. The results are shown in Figure 2, where the “error norm” is plotted as a function of the sampling rate. The error norm is the normalized pixel-by-pixel comparison of the reconstructed to original image, and the sampling rate (or density) is the ratio of the read pixels to the total number of pixels. T.-W. Lee et al. DOI: 10.4236/jgis.2021.133016 291 Journal of Geographic Information System Figure 2. The error norm (Equation (2)) as a function of the sampling rate for three sampling geometries: pyramidal, Cartesian and radial. The MODFLOW data in Figure 1 is used for this analysis. ( ) 2 Error Norm TRUE ij ij II−= ∑ (2) Iij is the pixel value, relative to the actual data, TRUE ij I . We can easily see that the pyramidal sampling geometry far outperforms the others, and using this sampling method error norm of 10 or less is possible at only 10% sampling rate, meaning that 1/10th of the data is necessary for a nearly full recovery of the data. Approximate reconstructions are possible at even less sampling rates. Thus, the ability to generate a full picture for the groundwater distribution is significantly enhanced using compressive sensing. This demonstration of the compressive sensing also points to the “compressibility” of groundwater data. As shown in Figure 1, the groundwater distribution is continuous (no step changes or spikes exist) and related to the adjacent data simply due to the physics of groundwater transport, regardless of strong effects of mountain ranges or other disruptive conditions. This means that the groundwater data can be decomposed using standard functions in Fourier series expansions, using a relatively small number of terms. 3. Results and Discussion For the groundwater level measurements, we use the USGS database, augmented by data provided by the state water resources management in the southwestern states (California, Nevada, Arizona, and Utah). We first focus on the groundwater resources in the southwestern United States, where the recent historical drought conditions make the accurate assessments and management important. T.-W. Lee et al. DOI: 10.4236/jgis.2021.133016 292 Journal of Geographic Information System Compressive sensing can be applied in a similar manner to any other regions where minimal amount of measurement data are available. The groundwater level data are typically generated using a pressure transducer or a floating device in registered wells. In each registered well, depth-to-well (the water head) is measured [15]. Data sets are obtained from USGS and respective water resource departments in Arizona, California, Nevada and Utah (ADWR, CDWR, NDWR and UDWR, respectively). The two pieces of information at each registered well, the depth-to-well and the site elevation (altitude), are converted to obtain the “groundwater level”, referenced to the sea level. For the site elevation, we use the NGVD29 (National Geodetic Vertical Datum of 1929), NAVD88 (National Geodetic Vertical Datum of 1988), or NED (National Elevation Dataset). The data are averaged over 0.25˚ latitude × 0.25˚ longitude areas, for data processing purposes ( i.e. to arrange the data in manageable array sizes, and also to average the fluctuations), as shown in Figure 3. It can be seen in Figure 3 that even with this segmentation and averaging, the data sampling is quite sparse, and there are large regions (white) where there aren’t any data available (e.g. some swaths in northeastern and southwestern Arizona). Thus, it would of high value to be able to track the groundwater in the entire area, through a minimally installed probe measurement. For validation of the compressive sensing algorithm, we start with the data set as shown in Figure 3(a), which shows the absolute groundwater level relative to the sea level for the year 2006, averaged to 0.25˚ latitude × 0.25˚ longitude cells. The spatial averaging is performed, to reduce the matrix size which exponentially reduces the computational time. Compressive sensing involves array operation, and the size of the array is one of the determining factors in reliable reconstruction. Again, white cells represent regions where no data are available. From this data set, we intentionally remove the data for 26 sites (cells) for a later comparison. The criterion for selecting these 26 sites is that there should be some nearby data locations, otherwise no validation is possible. Other than that, the selection was designed to be more or less evenly distributed. These sites are marked with inverted triangles in Figure 3(a), and also identified in Figure 3(b). The purpose is to test whether compressive sensing can recover the data at these sites, after “blinding” itself by removing these data from the analysis. Figure 4 shows the result of this test, where we have not only removed the data at the 26 sites shown in Figure 3(b), but also applied different sampling rates. The latter was done in order to find critical limits in the data sampling rate. Out of several sampling rates that we have applied, we exhibit three: 32 (unacceptable), 59 (acceptable) and 100% (full) sampling rates. The “city index” (the horizontal axis) represents the 26 sites shown in Figure 3(b), while the vertical axis is the “absolute” groundwater level relative to the sea level. At or below 32% sampling rate of the available data, the validation plot shows poor outcome where the compressive sensing data mostly tend to overestimate the groundwater level. However, above 59% sampling rate, the reconstruction becomes almost exact. Note that this 59% sampling rate is relative to the available data set, as T.-W. Lee et al. DOI: 10.4236/jgis.2021.133016 293 Journal of Geographic Information System (a) (b) Figure 3. (a) Groundwater absolute level relative to the sea level, averaged to 0.25 × 0.25 (longitude x latitude) resolution, inverted triangles are the 26 sites removed from the data for internal validation; and (b) Location and names of these sites. CA = California, NV = Nevada, UT = Utah, and AZ = Arizona. shown in Figure 3(a), which is already under-sampled. This means that even if we had only 59% of the data shown in Figure 3(a), minus the 26 data cells removed, we would still obtain accurate readings of the groundwater levels. Above this sampling rate, the data remain much the same, and overlap with the exact levels. Thus, Figure 4 shows that compressive sensing can convert sparse data set into a more complete spatial map of the groundwater levels. We can test the compressive sensing to perform this reconstruction over time, by applying it on an annual basis from 1992 to 2012, as shown in Figure 5. Here, we also compare T.-W. Lee et al. DOI: 10.4236/jgis.2021.133016 294 Journal of Geographic Information System (a) (b) (c) Figure 4. Validation of the compressive sensing reconstruction. (a) 32% sampling rate; (b) 59%; and (c) 100%. Circles denote measured data, and asterisks reconstructed data. Due to local altitude, the groundwater level in absolute ft can be high. the results with traditional methods of recovering missing data, e.g. Delaunay polygon, kriging, radial basis function, and inverse distance weighting methods. In hydrological monitoring, sparsity in the data is evidently a persistent problem, and several methods have been developed to “patch” this deficiency [15] [16]. After testing all of the above “interpolation” methods, kriging worked as well or slightly better (although still far from the actual data as shown in Figure 3) than the others, and therefore we only show the results of kriging method as open squares in Figure 5. Figure 5 shows that both the “under-sampled” (32%) compressive sensing and kriging deviate in a random manner from the actual (filled diamonds), and also from amply sampled (>59%) compressive sensing T.-W. Lee et al. DOI: 10.4236/jgis.2021.133016 295 Journal of Geographic Information System (a) (b) Figure 5. Point validation in (a) Cedar city, Utah and (b) Phoenix, Arizona from 1992 to 2012. data. When the sampling rate is above 59%, as in Figure 4, the data overlap exactly with the actual at all 26 locations and in time (only 2 sites are shown in Figure 5 as examples). Open triangles (>59%), dark squares (100%), and actual (filled diamonds) are all merged in Figure 5. Delaunay polygon, kriging, radial basis function, and inverse distance weighting methods have been traditionally used to compensate for lack of data [15] [16]. The main assumption is that water head can be interpolated in a linear or some other interpolation function. The Delaunay interpolation (or sometimes referred to as Thiessen polygon) is for defining regions based on a set of data points in the plane such that regions must be enclosed by a line midway between the measured point under consideration and surrounding points [17]. In geostatistics, kriging is an interpolation method for which the interplated values are modeled by prior covariance-governed Gaussian process, and is known as the best linear unbiased prediction method. Radial basis function uses a real-valued function whose value is dependent only on the distance from the origin. Using