A Hybrid Framework for Spatial Interpolation:
Merging Data-driven with Domain Knowledge

Abstract

Estimating spatially distributed information through the interpolation of scattered observation datasets often overlooks the critical role of domain knowledge in understanding spatial dependencies. Additionally, the features of these datasets are typically limited to the spatial coordinates of the scattered observation locations. In this paper, we propose a hybrid framework that integrates data-driven spatial dependency feature extraction with rule-assisted spatial dependency function mapping to augment domain knowledge. We demonstrate the superior performance of our framework in two comparative application scenarios, highlighting its ability to capture more localized spatial features in the reconstructed distribution fields. Furthermore, we underscore its potential to enhance nonlinear estimation capabilities through the application of transformed fuzzy rules and to quantify the inherent uncertainties associated with the observation datasets. Our framework introduces an innovative approach to spatial information estimation by synergistically combining observational data with rule-assisted domain knowledge.

Keywords: domain knowledge integration, spatial interpolation, spatial dependency correlation, neuro-fuzzy system

1 Introduction

Spatially dependent properties are widespread across various fields including subsurface resource exploitation [1, 2, 3, 4, 5], water resources management [6, 7], traffic engineering [8, 9], and environmental studies [10, 11]. These spatially varying attributes are typically recorded at a limited number of observed points or monitoring stations, often insufficient to represent the entire and multiscale spatial distribution accurately [12, 13, 5]. Due to the impracticality of monitoring every part of a domain to fully capture its attribute field, spatial interpolation is commonly employed to estimate values at unobserved locations, utilizing observed data and underlying spatial dependency correlations [14, 15].

Numerous spatial interpolation methods exist, as reviewed by Li and Heap [16], who suggest that these methods often rely on similar principles, such as inverse distance weighting [17] and Kriging [18]. The inverse distance weighting approach, which is built on a straightforward distance-based method, is generally unsuitable for spatially aggregated observations. Kriging, on the other hand, functions as an unbiased linear estimator under the assumption of a stationary normal distribution. Since it utilizes the conditional distribution of a Gaussian process, Kriging is essentially a specific case of Gaussian process regression [19, 20]. However, all these methods are founded on simplified assumptions, including stationarity, linearity, and normality, which limit their applicability in real-world scenarios [20, 21, 22]. In particular, their reliance on linear interpolation makes them inadequate for handling complex nonlinear spatial dependency relationships, often resulting in suboptimal estimations for continuous attribute fields. Therefore, there is a clear need for alternative approaches that can construct non-stationary, nonlinear, and non-Gaussian spatial dependency functions for spatial interpolation.

Given their versatile nonlinear approximation capabilities, the emerging data-driven approaches are promising to provide an alternative framework for spatial information modeling [11, 23, 24, 25, 26]. Various machine learning-based methods, such as random forests [27, 28, 29], artificial neural networks [30, 31, 32], radial basis function networks [21, 33], long short-term memory networks [34, 35], convolutional neural networks [36, 37], conditional generative adversarial networks [22, 38, 39], and ensemble learning [40, 41], have been successfully adapted and applied to spatial interpolation. However, without significant reengineering, these original machine learning methods are not inherently designed for spatial interpolation and often fail to account for the available spatial configuration information. In other words, directly applying these machine learning models may pose challenges and may not be suitable for most spatial interpolation scenarios [42, 43].

While methods involving convolutional neural networks and generative adversarial networks are well-suited for imagery or regularly gridded datasets, such as terrain elevation images, observational data are often irregularly sparse and scattered [44, 45]. Traditional deep learning methods may not be ideal for extracting spatial dependency features from such sparse observations, as they typically only rely on the spatial coordinates of the observation locations as input features [46]. To address this limitation, modified approaches have been developed. For example, Wu et al. proposed an inductive graph neural network Kriging model that better leverages distance information [47], while Ma et al. introduced a geo-layer into long short-term memory networks to integrate spatial correlations from monitoring stations [48]. In essence, effectively extracting spatial dependency features from irregularly sparse observed data necessitates considering both individual observations and their neighboring spatial configurations.

In addition to the data-driven approaches, the family of rule-based fuzzy inference systems can effectively model nonlinear functions by leveraging inherent reasoning paradigms [49, 50, 51, 52, 53]. Unlike data-driven methods, fuzzy inference systems can directly and broadly incorporate expert knowledge, making them valuable for enhancing spatial interpolation [52]. For example, the first law of geography, which states that “near things are more related than distant things” [54], can be articulated as a fuzzy IF-THEN rule—IF two spatial observations are closer, THEN their spatial dependency correlations are stronger. This rule underpins many existing interpolation methods that rely on concepts of distance and neighborhood. Similarly, there is a wealth of domain-specific knowledge related to both general spatial dependencies and particular interpolation scenarios. For instance, Yesilkanat et al. [55] demonstrated the successful application of domain expertise for spatial interpolation of environmental radioactivity using fuzzy IF-THEN rules, either collected or self-defined. However, their fuzzy rulesets were tailored to specific domains rather than forming a general rule base for spatial dependency extraction. Therefore, there is a need to develop a more general rule-based spatial interpolation framework that can effectively exploit and utilize latent domain knowledge. It is important to note that domain knowledge related to spatial dependency is often difficult to collect and represent in the form of fuzzy IF-THEN rules. Moreover, the lack of standardized approaches to transforming domain knowledge into rule bases limits the effective use of such knowledge. Consequently, it is both necessary and advantageous to systematically and automatically generate fuzzy rulesets to incorporate domain knowledge. Another advantage of rule-based fuzzy inference systems is their ability to tolerate inaccurate information. Observation data and monitoring records are often subject to errors from various sources [56], which can negatively impact interpolation accuracy. A robust spatial interpolation approach should account for these uncertainties associated with the observed data. Notably, fuzzy logic is well-suited for handling such inherent uncertainties, thereby enhancing the performance of spatial interpolation [57, 58].

In this paper, we present a hybrid framework merging data-driven and domain knowledge for interpolating spatially dependent properties. The framework extracts latent spatial dependency basis by leveraging observation data and neighboring information. By embedding a fuzzy inference system within our adaptive network, we can automatically transform domain knowledge into fuzzy rulesets. This framework capitalizes on the advantages of fuzzy reasoning, particularly its ability to tolerate inaccurate information and manage nonlinear spatially dependent functions. We validate our framework in two scenarios: subsurface formation parameter estimation and air quality mapping. Additionally, we conduct comparative studies to assess the estimation performance quantitatively and qualitatively against conventional interpolation techniques, including ordinary Kriging, inverse distance weighting, and Gaussian process regression.

2 Hybrid data-driven and rule-assisted learning framework

To address the limitations mentioned above, our framework combines the extraction of spatial dependency features from observation data with the transformation of latent domain knowledge into fuzzy rulesets. For constructing a fuzzy inference system, it is ideal to automatically convert domain knowledge into fuzzy IF-THEN rules. The Adaptive-Network-based Fuzzy Inference System (ANFIS) developed by J.S. Jang [49] provides a solid foundation for this task. However, the number of rules generated by ANFIS can become intractable as the number of input features increases. Specifically, if the number of membership functions assigned to each input dimension is fixed and the dimensionality increases linearly, the number of generated rules grows exponentially. This exponential growth can significantly limit the applicability of ANFIS when incorporating both spatial coordinates and relevant neighboring information as input features. Inspired by ANFIS, we propose an enhanced architecture to overcome this issue by input feature decomposition. As shown in Figure 1, our architecture integrates a data-driven approach with rule-based assistance. The network consists of an intrinsic input layer, a ruleset layer, a T-Norm operation layer, a normalization layer, a consequent layer, and a summation layer. Additionally, it includes a newly designed spatial dependency layer to exploit sparse observation data and their neighboring information. Essentially, our architecture retains the ANFIS learning mechanism for constructing fuzzy IF-THEN rules and approximating the estimation function.

Figure 1: The network architecture of the hybrid framework based on ANFIS

The detailed architecture of our hybrid framework is illustrated in Figure 2. It has two main components: 1) data-driven Spatial Dependency Basis (SDB) extraction, and 2) rule-assisted spatial dependency function approximation. The SDB serves as a critical link between these two components, significantly reducing the number of generated rules. In the following sections, we outline the main steps and their corresponding implementation.

Figure 2: The detailed architecture of the hybrid framework

2.1 Spatial Dependency Basis (SDB) extraction

We apply the reduced-rank approach for obtaining SDB [59]. The SDB is extracted from spatial observations to represent fixed spatial features. Rather than using the direct input features from these spatial observations, we use SDB as the input for the rule-based ANFIS. This approach aims to minimize the number of automatically generated fuzzy IF-THEN rules while fully leveraging the available spatial observations. As illustrated in Figure 2, the extraction of the SDB begins with constructing nearest neighboring spatial covariates.

2.1.1 Nearest neighboring spatial covariates

We use nearest neighboring spatial covariates [27, 42], which is the combination of the environmental covariates with the spatial covariates, to comprehensively describe the spatial dependency relationship from spatial observations. Compared to approaches that consider only spatial coordinates as input features, these combined covariates better characterize spatial correlations by fully accounting for spatial coordinates, observation data, and neighboring configurations. The neighboring configurations are established by constructing a neighboring graph for each observation location. We use nearest neighboring algorithm to select the m nearest neighbors for each observation location based on Euclidean distance, see Equation (1).

$$d_{ij}=\|\mathbf{Obs}_{i}-\mathbf{Obs}_{j}\|_{2}$$ (1)

where $\mathbf{Obs}_{i}$ and $\mathbf{Obs}_{j}$ are the $i$-th and $j$-th observed points, respectively. In 2D, $\mathbf{Obs}_{i}=(x_{i},y_{i})$, and in 3D, if necessary, $\mathbf{Obs}_{i}=(x_{i},y_{i},z_{i})$. The value $d_{ij}$ represents the Euclidean distance between $\mathbf{Obs}_{i}$ and $\mathbf{Obs}_{j}$.

As illustrated in Figure 3, the observation locations and their corresponding neighbors form a neighboring graph. This graph is then used to construct the associated features, allowing the features at a specific observation location to encompass both its own data and the information from its neighbors.

Figure 3: An illustration of Neighboring graph

We can express the neighboring graph using nearest neighboring spatial covariates to generate more valuable features for each observation. For a given observation $\text{Obs}_{i}=(x_{i},y_{i})$ with observation value of $\Phi_{i}$ and the $m$ nearest neighbors, the corresponding neighboring spatial covariates are written as:

$$\text{[}{Obs}_{i}]=[x_{i},y_{i},x_{i}^{1},y_{i}^{1},\Phi_{i}^{1},x_{i}^{2},y_{i}^{2},\Phi_{i}^{2},\dots,x_{i}^{j},y_{i}^{j},\Phi_{i}^{j},\dots,x_{i}^{m},y_{i}^{m},\Phi_{i}^{m}]_{1\times(3m+2)}$$ (2)

where $x_{i}$ and $y_{i}$ are the spatial coordinates of the $i$-th observed point, and $i=1,2,\dots,N$ if there are $N$ observations in total. For the case with $N$ observed points, $\text{Obs}_{1},\text{Obs}_{2},\dots,\text{Obs}_{i},\dots,\text{Obs}_{N}$, we can construct a matrix of nearest neighboring spatial covariates as shown in Equation (3). The superscript $j$ represents the $j$-th nearest neighbor of the $i$-th observed point, and $j=1,2,\dots,m$.

For each observation location, the dimension of its nearest neighbor spatial covariates is $3m+2$. Assigning two membership functions to each dimension results in the adaptive neural networks automatically generating $2^{3m+2}$ fuzzy IF-THEN rules. This can lead to a significant number of generated rules, causing computational challenges in tuning and updating the hyperparameters during training to accurately map the nonlinear spatial dependency function. To address this, we further extract primary features that can characterize the spatial correlations using SDB.

2.1.2 Reduced-rank basis extraction

We use Uniform Manifold Approximation and Projection (UMAP) [60] to decompose the nearest neighboring spatial covariates into several fixed reduced-rank bases specifically targeting spatial dependency basis. According to Equation (3), the dimension of the nearest neighboring Spatial Covariates (SC) is $3m+2$. If we want to reduce the dimensionality to $d$, we can write such reduction as in Equations (4) and (5) to extract the Spatial Dependency Basis (SDB).

$$SDB_{i}=\begin{bmatrix}SDB_{i}^{1}&SDB_{i}^{2}&\dots&SDB_{i}^{d}\end{bmatrix}\in\mathbb{R}^{d}$$ (5)

Where $i$ denotes the $i$-th observation location.

Then we construct a local relationship graph [61] based on the SC by converting the differences between neighbors into weights or probabilities using Equation (6).

$$W_{ij}=\exp\left(-\frac{d(SC_{i},SC_{j})-p_{i}}{\sigma_{i}}\right)$$ (6)

Where $d(SC_{i},SC_{j})$ is the distance between the $i$-th and $j$-th observation.

Apparently, a shorter distance implies a stronger spatial dependency relationship. Hence, $d(SC_{i},SC_{j})$ is considered as the spatial dependency between these two observations. $p_{i}$ is the local spatial dependency between the $i$-th observation and its adjacent neighbor. $\sigma_{i}$ is the local spatial dependency between the $i$-th observation and its $m$ nearest neighbors, which is regularized as $\sum_{j}W_{ij}=\log_{2}(m)$.

We then build an embedding spatial dependency basis which can preserve the structure of the local relationship graph between the nearest neighboring spatial covariates. This basis is optimized to minimize the difference between SC and SDB with respect to fuzzy set cross entropy in Equations (7) and (8).

$$C(W_{ij},\mu_{ij})=\sum_{ij}{}\left[W_{ij}\log\left(\frac{W_{ij}}{\mu_{ij}}\right)+(1-W_{ij})\log\left(\frac{1-W_{ij}}{1-\mu_{ij}}\right)\right]$$ (7)

$$\mu_{ij}=\left[1+a\cdot d(SDB_{i},SDB_{j})^{2b}\right]^{-1}$$ (8)

The local relationship graph between extracted SDB is characterized by the weights $\mu_{ij}$. $d(SDB_{i},SDB_{j})$ denotes the distance between the $i$-th observation location and the $j$-th observation location with respect to the dataset of SDB. $a$ and $b$ are adjusted by non-linear least squares fitting.

Further implementation of UMAP, such as the fuzzy topological representation and parameter calculation, can be found in [60, 61].

2.2 Spatial dependency estimation

Estimating a parameter field involves approximating the spatial dependency function using both the Spatial Dependency Basis (SDB) and domain knowledge, particularly in the form of automatically generated fuzzy IF-THEN rules. As illustrated in Figure 2, our rule-assisted adaptive network transforms domain knowledge related to spatial dependency into a rule base.

ANFIS combines the nonlinear function approximation capabilities of neural networks with the knowledge utilization strengths of fuzzy inference systems. In this context, we adapt ANFIS to model the nonlinear spatial dependency function by generating latent knowledge from the dataset of spatial dependency bases. For example, in a 4-dimensional case, the transformed knowledge in the form of IF-THEN rules can be expressed as:

IF $$SDB_{i1}\text{ is }\mathit{Rule}_{11}\text{ and }SDB_{i2}\text{ is }\mathit{Rule}_{21}\text{ and }SDB_{i3}\text{ is }\mathit{Rule}_{31}\text{ and }SDB_{i4}\text{ is }\mathit{Rule}_{41}$$ (9)

THEN $$\Phi_{i}=\overrightarrow{C_{i}}\cdot\overrightarrow{SDB_{i}}$$

where $\overrightarrow{SDB_{i}}=\left[SDB_{i1},SDB_{i2},SDB_{i3},SDB_{i4},1\right]$ includes the 4-dimensional spatial dependency basis and one as a constant. $C_{i}$ is the consequent parameter matrix to be discussed later, and $\Phi_{i}$ is the estimated value at an unknown location. $\mathit{Rule}_{11},\mathit{Rule}_{21},\mathit{Rule}_{31},\text{and }\mathit{Rule}_{41}$ are the corresponding defined linguistic labels of each spatial dependency basis, respectively.

A complete IF-THEN rule consists of both the antecedent clause in the IF part and the consequent clause in the THEN part. Typically, linguistic labels are often defined such as “Close,” “Medium,” and “Fair.” In this study, we refer to these linguistic labels as $\mathit{Rule}_{ij}$ for convenience. $\mathit{Rule}_{ij}$ represents the $j$-th linguistic label of the $i$-th spatial dependency basis. Additionally, it is characterized by a bell membership function as shown in Equation (10):

$$\mu_{\mathit{Rule}}=\frac{1}{1+\left[\left(\frac{s-g}{e}\right)^{2}\right]^{f}}$$ (10)

where $s$ is the input data into the rules-set layer (shown in Figure 2), $\mathit{Rule}$ is the linguistic label, and $e,f,g$ are premise parameters to be adjusted during the training process. The bell membership function is important as it can tolerate imprecise information by assigning a certain degree (between 0 and 1) of membership to each input $s$ based on fuzzy logic. This quantifies the uncertainties associated with the observation data.

Within our ANFIS architecture, the rules-set layer is predefined to set the number of membership functions (or the number of linguistic labels) directly, which determines the number of the finally generated IF-THEN rules. Each node in this layer is given by Equation (11):

$$O_{\mathit{l,k}}^{Rule-set}=\mu_{\mathit{Rule}_{l,k}}(x)$$ (11)

where $O_{\mathit{l,k}}^{Rule-set}$ is the $k$-th node or membership function of the $l$-th dimensional input spatial dependency basis that makes the $SDB_{l}$ fuzzy between 0 and 1, where $l=1,2,\dots,d$ and $d$ is the dimensionality of the extracted spatial dependency basis. In this study, we assign two membership functions for each input, so $k=1,2$. As a result, we can write the premise parameters in $\vec{P}=[e_{11},f_{11},g_{11},e_{12},f_{12},g_{12},\dots,e_{d1},f_{d1},g_{d1},e_{d2},f_{d2},g_{d2}]_{1\times 6d}$.

The T-Norm layer performs a multiplication operation on the signals from the previous layer. It outputs the weight $\omega_{t}$ or the firing strength of a rule, according to Equation (12):

$$O_{\mathit{t}^{TNorm}}=\omega_{t}=\mu_{\mathit{Rule}_{l_{1}k_{1}}}(s_{1})\cdot\mu_{\mathit{Rule}_{l_{2}k_{2}}}(s_{2})\cdot\mu_{\mathit{Rule}_{l_{3}k_{3}}}(s_{3})\cdot\mu_{\mathit{Rule}_{l_{4}k_{4}}}(s_{4})$$ (12)

where $O_{\mathit{t}}^{TNorm}$ is the $t$-th weight and $t=1,2,3,\dots,2^{d}$, $l_{1}\neq l_{2}\neq l_{3}\neq l_{4}$, and $k_{1},k_{2},k_{3},k_{4}=1,2$.

The normalized layer calculates the ratio of the $t$-th firing strength to the sum of all rules’ firing strengths using Equation (13):

$$O_{\mathit{t}}^{Normalized}=\overline{\omega_{t}}=\frac{\omega_{t}}{\sum\omega_{t}},\quad t=1,2,3,\dots,2^{d}$$ (13)

The consequent layer computes the result of the consequent clause of the IF-THEN rules following Equation (14). The parameters in this layer, known as consequent parameters, are adjusted during the training process using least squares estimation:

$$O_{\mathit{t}}^{Consequent}=\overline{\omega_{t}}f_{t}=\overline{\omega_{t}}(\overrightarrow{C}\cdot\overrightarrow{SDB})=\overline{\omega_{t}}(p_{1}SDB_{1}+p_{2}SDB_{2}+\dots+p_{d}SDB_{d}+p_{d+1})$$ (14)

where $\overrightarrow{SDB}=[SDB_{1},SDB_{2},\dots,SDB_{d},1]$ and $\overrightarrow{C}=[p_{1},p_{2},\dots,p_{d},p_{d+1}]$.

The summation layer calculates the overall output by summing all the incoming signals from the previous layer, as shown in Equation (15):

$$O_{\mathit{summation}}=\sum O_{t}^{Consequent}=\sum\overline{\omega_{t}}f_{t}$$ (15)

2.3 Implementation of the proposed hybrid learning framework

We outline the key implementation steps of our hybrid framework. First, we obtain the nearest neighboring spatial covariates by constructing neighboring graph for each observation point. We then use these covariates to extract the latent spatial dependency basis as the main input features for our rule-assisted adaptive networks. Next, we model and train the spatial dependency function using the rule-based ANFIS. Finally, we approximate the nonlinear spatial function, which allows us to perform spatial interpolation of a specific attribute field. Table 1 lists the corresponding pseudo code.

Table 1: Pseudo code of hybrid data-driven and rule-assisted learning procedure

Algorithm: Hybrid data-driven and rule-assisted learning procedure
Input: Observation record $\mathit{Obs}=\left[\mathit{Obs}_{1},\mathit{Obs}_{2},\dots,\mathit{Obs}_{i},\dots,\mathit{Obs}_{N}\right]$
Output: Premise parameters $\vec{P}$ and consequent parameters $\vec{C}$ of fuzzy IF-THEN rules for establishing the approximated spatial dependency function
Step 1: Obtain nearest neighboring Spatial Covariates (SC) for each observation point
Step 1.1: Construct neighboring graph for each observation point
Function: NeighboringGraph($Obs,m$)
For each observation $\mathit{Obs}_{i}\in Obs$:
# Return $\mathit{m}$ nearest neighbors of $\mathit{Obs}_{i}$ in all observations
Neighbors $\leftarrow$ NN($\mathit{Obs}_{i},Obs,m$)  # nearest neighbors algorithm
NeighborsMatrix.append(Neighbors)
Return NeighborsMatrix
Step 1.2: Build SC from NeighborsMatrix
Step 2: Extract SDB from SC using UMAP
# the extraction is performed by minimizing cross-entropy loss
Step 3: Model and train the spatial dependency function using rule-based ANFIS
Function: ANFIS($SDB,\Phi$)
Step 3.1: Update premise parameters by the gradient descent method
$\vec{P}\leftarrow\vec{P}-\eta\frac{\partial(\text{error})}{\partial\vec{P}}$  # adjust premise parameters
Step 3.2: Update consequent parameters by least squares estimation
$\vec{C}=(SDB^{T}SDB)^{-1}SDB^{T}\Phi$  # using Kalman filtering algorithm to calculate $\vec{C}$
Return $\vec{P}$ and $\vec{C}$

3 Applications

To evaluate the general applicability and practical performance of the proposed framework, we apply it to two different cases involving the estimation of spatially dependent properties using interpolation of scattered observation data. The first case addresses a classic challenge in subsurface formation characterization [1], where the goal is to estimate the distribution of formation properties, such as rock porosity and/or permeability, based on very scattered observations, and produce a continuous property distribution map for subsurface resource exploitation. This is a particularly relevant use case as only limited observation data are available spatially in general [14, 15]. The second case focuses on estimating the distribution of air pollutants by interpolating data collected from scattered observation stations.

3.1 Spatial interpolation for oil reservoir porosity mapping

We use an oil reservoir example as a representing case. We aim to generate the formation porosity distribution using limited porosity measurements at scattered spatial observation locations. The benchmark model is provided in Figure 4. This model contains 50-by-50-by-1 grid blocks.

Figure 4: Oil reservoir model and its porosity distribution

We irregularly sample the porosity field, as shown in Figure 4 to mimic the scattered observation locations with porosity measurement data, see Figure 5. We then randomly choose 100 observation locations as the input for the proposed hybrid learning framework. By modeling a proper spatial dependency function, we obtain the porosity estimation for unsampled locations in this model with 2500 grid blocks.

Figure 5: Scattered porosity data at randomly chosen observation locations

As shown in Figure 6, the estimated porosity distribution map closely recaptures the original spatial pattern in the benchmark model. The high porosity channel and the low porosity channel along the lower-left diagonal are apparently reconstructed. Additionally, several dispersed regions with significant porosity variations are also qualitatively estimated, such as the high porosity region in the upper-right corner and the low porosity zone in the lower-left corner.

Figure 6: Porosity distribution constructed using our proposed hybrid learning framework

We further compare the estimated porosity values with the actual benchmark values using line plots as shown in Figure 7. The red dotted line represents the estimated porosity for each grid block, while the black solid line shows the actual porosity. Clearly, our estimation aligns well with the true porosity data; however, upon closer inspection of the curves, some minor discrepancies between the actual and predicted porosity values appear. To address these differences, we perform a quantitative analysis of the estimation performance using some evaluation metrics.

The evaluation metrics include Mean Square Error (MSE), Root Mean Square Error (RMSE), Mean Absolute Error (MAE), Mean Absolute Percentage Error (MAPE), and coefficient of determination ($R^{2}$), as given in the following Equations.

MSE $$=\frac{1}{n}\sum_{i=1}^{n}(y_{i}-\hat{y}_{i})^{2}$$ (16)

RMSE $$=\sqrt{\frac{1}{n}\sum_{i=1}^{n}(y_{i}-\hat{y}_{i})^{2}}$$ (17)

MAE $$=\frac{1}{n}\sum_{i=1}^{n}|y_{i}-\hat{y}_{i}|$$ (18)

MAPE $$=\frac{1}{n}\sum_{i=1}^{n}\frac{|y_{i}-\hat{y}_{i}|}{\max(y_{i})}$$ (19)

$$R^{2} =1-\frac{\sum_{i=1}^{n}(y_{i}-\hat{y}_{i})^{2}}{\sum_{i=1}^{n}(y_{i}-\bar{y})^{2}}$$ (20)

Where $\hat{y}_{i}$ and $y_{i}$ are the estimated and actual porosity of the $i$-th grid, respectively. $i=1,2,3,\dots,n$, and $n$ is the total number of grid blocks in the model. The symbol $|\cdot|$ represents the absolute value, and $\bar{y}$ is the mean of actual porosities.

Figure 7: comparative curve between actual value and estimated predictions for each grid

The calculated metrics are given in Table 2. The errors are very low in terms of MSE, RMSE, MAE, and MAPE. The coefficient of determination ($R^{2}$) is high at 0.9283. Furthermore, we demonstrate the fitted straight line between estimated and actual porosities for a better illustration of $R^{2}$ in Figure 8. The regression line $y=x$ indicates that the estimation closely matches the actual porosity, underscoring the reliability of our proposed spatial interpolation method.

Table 2: The evaluation metrics for our proposed method

MSE	RMSE	MAE	MAPE	$R^{2}$
0.00011236	0.0106	0.0086	0.0391	0.9283

Figure 8: Interpolation performance of each grid using scatter plot along with regression straight line

3.2 Spatial interpolation for ozone value mapping

The second case is to estimate ozone value distribution based on sparse measurement data at scattered air quality monitoring stations. The ozone data used in this study is publicly available from the U.S. Environmental Protection Agency (EPA) Air Quality System (https://aqs.epa.gov/aqsweb/airdata/). Those ozone values are based on the U.S. national ambient air quality standards, specifically in the form of the “annual fourth-highest daily maximum 8-hour concentration”. In our application, we select the observed ozone data recorded at the same time across various monitoring stations as our test dataset, as illustrated in Figure 9. We then conduct spatial interpolation across the entire continental United States to estimate ozone levels in areas without direct observations. This allows for an assessment of air quality, even in regions with limited monitoring infrastructures.

The map generated using our proposed method captures the higher ozone areas in the west and lower areas in the east, as clearly and intuitively visualized in Figure 10. The reconstructed ozone distribution demonstrates the capabilities in capturing spatial variations, effectively distinguishing between highly contaminated regions and less affected areas. This is of practical use for representing the whole air quality situation in sparsely monitored regions, where fewer monitoring stations are available. The ability to reflect these variations in regions with limited data highlights the robustness and practical utility of our spatial interpolation method for air quality monitoring and assessment.

Figure 9: Ozone value at scattered monitoring station in continental United States

Figure 10: Ozone value distribution interpolated using our hybrid framework

In addition to visualizing the spatial distribution patterns of ozone levels, we also provide quantitative results using the evaluation metrics mentioned earlier. Unlike the first case, where actual values at unknown locations are available for comparison, the actual ozone values at non-monitoring stations in this scenario are not accessible. Therefore, we cannot directly compare the estimated values with actual measurements at those locations. We employ the 10-fold cross-validation to conduct the comparison. We divide the observed data into a training dataset and a validation dataset to assess the model’s performance. We then calculate the average values for MSE, RMSE, MAE, MAPE, and $R^{2}$ across the 10 validation datasets. The results of these evaluations are listed in Table 3.

Table 3: The evaluation metrics of 10-fold validation
datasets using our proposed method

MSE	RMSE	MAE	MAPE	$R^{2}$
0.00002601	0.0051	0.0037	0.0792	0.7626

4 Comparison and Discussion

We compare our hybrid framework with several established interpolation methods, including Inverse Distance Weight (IDW) [17], ordinary kriging [18], and gaussian process [62]. These three comparative methods are widely used for spatial interpolation.

4.1 Comparison in terms of evaluation metrics

Table 4 and Table 5 provide the respective comparative results. Our proposed method outperforms in terms of MSE, RMSE, MAE, MAPE, and $R^{2}$ for both cases.

Table 4: Comparative evaluation metrics of different interpolation methods for reservoir porosity mapping

Metrics	Ordinary Kriging	IDW	Gaussian Process	Our Method
MSE	0.00046976	0.00113792	0.00030674	0.00011236
RMSE	0.02167388	0.03373304	0.03373304	0.0106
MAE	0.01617806	0.02701021	0.01324166	0.0086
MAPE	0.15395796	0.2807828	0.13766530	0.0391
$R^{2}$	0.71519912	0.31011110	0.81403229	0.9283

Table 5: Comparative evaluation metrics of different interpolation methods for ozone value mapping

Metrics	Ordinary Kriging	IDW	Gaussian Process	Our Method
MSE	0.00003324	0.00002602	0.00002665	0.00002601
RMSE	0.00576572	0.00510121	0.00516196	0.0051
MAE	0.00447331	0.00411962	0.00411051	0.0037
MAPE	0.13360340	0.12824334	0.12676602	0.0792
$R^{2}$	0.29874051	0.45106761	0.43791603	0.7626

4.2 Comparison in terms of spatial feature reconstruction

In addition to quantitative analysis in terms of error and accuracy, we also compare the spatial distribution patterns from both global and local perspectives. The comparative spatial distribution maps of porosity are shown in Figure 11. It is evident that all three traditional methods can reconstruct the global trends and features. Global features here mean the prominent spatial distribution patterns that are easily recognizable across the entire map. These global patterns include the high porosity channel or region in the lower-left diagonal and upper-right direction, as well as the low porosity channel along the diagonal.

However, when it comes to capturing spatially local features, all three traditional methods fall short. Local features refer to the finer details and variations within smaller regions that contribute to the overall spatial heterogeneity. The inability of these methods to accurately reconstruct these local features highlights a limitation in their capacity to fully represent complex spatial variability. This comparison underscores the importance of methods that can address both global and local spatial patterns to provide a more comprehensive and detailed understanding of the spatial distribution. As demonstrated in Figure 11(d), the circled regions highlight the local features that reflect variations within the porosity field. Compared to other commonly used techniques, our proposed method captures and represents these local features more effectively, providing a more detailed and nuanced depiction of the spatial variability. This enhanced ability to preserve local variations underscores the strength of our approach in handling complex spatial patterns.

Figure 11: Comparative results of porosity interpolation generated from different methods

Similarly, the comparative spatial distribution maps of ozone values are shown in Figure 12. All methods can predict global trends to some extent. For instance, the highly polluted region in the southwest and the less polluted region in the east are consistently identified across all methods. However, when it comes to reconstructing local features, our proposed method outperforms the other three interpolation techniques, as indicated by the dashed circles.

It should be noted that the three commonly used interpolation techniques tend to produce overly smooth interpolations. While they may capture general trends, it often comes at the cost of losing important local variations within the spatial field. Excessively smooth interpolation can negatively impact the accuracy of spatial predictions by overlooking these critical local features or variations in the spatial dependency field. Our proposed method, by contrast, strikes a better balance between capturing global trends and preserving local details, leading to more accurate and nuanced spatial interpolations.

Figure 12: Comparative results of ozone interpolation generated from different methods

4.3 Sensitive analysis with respect to the number of nearest neighbors

Although our proposed method offers advantages over traditional techniques in terms of accuracy and the preservation of local features in both cases, its performance in interpolating ozone is slightly less effective compared to its success in reconstructing the porosity field. The primary difference between these two cases is the density of observed data within the interpolation space. The denser data in the porosity field allows for more precise predictions, whereas the sparser data in the ozone case makes it more challenging to capture local variations accurately.

Compared to the first porosity field case, the observation density in the ozone case is less uniform than in our randomly selected porosity data. This nonuniform observation density tends to be more sensitive to the number of nearest neighbors m when constructing nearest neighboring spatial covariates in our proposed method. In other words, the choice of the number of nearest neighbors has certain impact on the spatial distribution patterns and the estimation accuracy of the reconstructed attribute fields. This sensitivity underscores the importance of carefully selecting the number of nearest neighbors to ensure high quality spatial interpolation, particularly in cases with uneven observation densities.

To further illustrate this point, we focus on the ozone case to discuss the sensitivity of the number of nearest neighbors m. We compare and analyze different hyperparameter settings, as demonstrated in Figure 13. This comparison highlights how varying the number of nearest neighbors affects the spatial interpolation results, allowing us to better understand the impact of this parameter on the accuracy and reliability of our estimated predictions.

Figure 13: Reconstructed spatial distribution patterns with different number of nearest neighbors

Clearly, different values of m have significant impact on the interpolated spatial distribution and its corresponding accuracy. We choose to use 10 nearest neighbors to achieve the desired accuracy. However, it is important to note that the optimal value of m is highly dependent on the specific observed sampling data. Considerable time and effort are required to carefully tune the hyperparameter m to achieve higher interpolation quality. This tuning process is crucial for adapting the model to different datasets and ensuring reliable interpolation results.

Additionally, a larger value of m tends to preserve global spatial patterns, while a smaller value of m tends to emphasize local features. For instance, when m is set to 30, the estimated high-ozone value region in the middle appears more expansive. Conversely, when m is reduced to 5, the results reveal some sparse low-ozone value areas within the high-ozone value region. This occurs because a larger m means that a broader neighborhood is considered when estimating unknown locations, leading to results that are more influenced by higher attribute values within that neighborhood. On the other hand, a smaller m often overlooks the general trend due to its limited neighborhood size, resulting in a focus on local variations. This, however, makes the interpolations more susceptible to outliers or extreme values.

Therefore, neither a large nor a small value of m is ideal for constructing nearest neighboring spatial covariates. Careful tuning of the hyperparameter m is essential when implementing our hybrid data-driven and rule-assisted learning framework for spatial interpolation. Finding the right balance in m allows for a better representation of both global and local patterns, ensuring reliable and robust estimations.

5 Conclusions

In this study, we develop a hybrid data-driven and rule-assisted learning framework that combines spatial feature extraction with domain knowledge integration via fuzzy rule sets to enhance the interpolation of spatially dependent properties. By decomposing nearest neighbor spatial covariates into multiple spatial dependency bases and utilizing fuzzy IF-THEN rules within an adaptive network framework, our method accommodates imprecise information, considers data uncertainties, and improves spatial interpolation accuracy. Validated through applications in subsurface formation characterization and air quality assessment, our approach surpasses traditional techniques like ordinary Kriging, inverse distance weighting, and Gaussian processes by achieving lower error metrics and better capturing local spatial features. However, the method’s performance is sensitive to the number of nearest neighbors used, which influences the balance between global trends and local variations. While our approach is parametric and requires careful tuning, future research could explore non-parametric methods that dynamically adjust the number of nearest neighbors, offering a promising avenue for further improving spatial interpolation techniques.