## 1 Introduction  
Spatial interpolation is a computational technique used to predict attribute values at any location in a region. It involves filling the gaps between sampled points and converting discrete sample data into a continuous surface [1, 2, 3]. This method is based on [PERSON]'s [4] first law of geography and the principle of distance decay. Currently, spatial interpolation is widely used in various fields such as spatial analysis, 3D geological modeling, management of geographic data, smart city geological management, and spatial assessment [5, 6, 7, 8, 9, 10]. Recently, as China's coal mining industry moves towards intelligent mining practices [11; 12; 13], constructing high-precision 3D geological models has become a critical technical necessity. The quality and reliability of these models directly hinge upon the selection of a scientifically sound and appropriate spatial interpolation method [14; 15]. Spatial interpolation methods are broadly classified into three categories: deterministic interpolation, geostatistical interpolation, and machine learning interpolation. Deterministic interpolation relies on the similarity of sample points or surface smoothness to develop a predictive model. Common deterministic methods include inverse distance weighting (IDW) and radial basis function (RBF) [16]. [PERSON] et al. [17] employed IDW for 3D geological modeling of groundwater levels. [PERSON] et al. [18; 19] conducted a systematic analysis of IDW parameters and introduced the concept of \"Relative distance\" for 3D geological modeling. Addressing the issue of RBF stability, [PERSON] [20] proposed a method for 3D geological modeling utilizing scattered data. [PERSON] et al. [21] explored a method for reconstructing implicit curves and surfaces through RBF interpolation. Deterministic interpolation algorithms offer advantages in terms of simplicity, flexibility, and speed. However, they require a high density of spatially distributed sample points. Therefore, interpolation precision reduces when sample points are sparsely distributed. Geostatistical interpolation analyzes natural characteristics by referencing spatial or spatiotemporal data. Central to this field are kriging [22] and its variations, including ordinary kriging, simple kriging, probabilistic kriging, and directed kriging, along with their optimized algorithms [23; 24]. These methods utilize the semivariogram to characterize spatial autocorrelation and achieve optimal, linear, and unbiased estimations. [PERSON] et al. [25] developed a spatially interpolated precipitation prediction employing collaborative kriging. [PERSON] et al. [26] generated a multi-attribute 3D geological model of a high-geostress tunnel utilizing indicator kriging. [PERSON] et al. [27] evaluated the effectiveness of discrete smooth interpolation (DSI), kriging, and function interpolation for transparent surfaces by analyzing the resultant errors. [PERSON] et al. [28] built a 3D fault model utilizing weighted kriging. However, kriging methods necessitate stable, consistently distributed data, and when significant spatial variability exists between sampling points, the resulting interpolation is excessively smooth, failing to capture localized fluctuations. With advances in geostatistics, Bayesian theory and its derived methods are widely applied in fields such as engineering geology, spatiotemporal interpolation, and public health[PERSON] [PERSON] et al. [10; 29; 30; 31] precisely identified the soil-rock [PERSON] using Bayesian evidence theory. They inverted issues such as rainfall-induced slope instability based on Bayesian theory. In the field of spatiotemporal interpolation, Bayesian Maximum Entropy (BME) interpolation [32; 33] has become prominent in the field. In contrast to kriging, BME is knowledge-driven, capable of integrating diverse data sources with varying levels of precision under a rigorous and comprehensive theoretical framework. This enhances predictive precision and offers adaptability [34; 35][PERSON] [PERSON] et al.