Adaptive Inversion of the Curie Isothermal Surface in the Ordos Basin and Surrounding Regions Based on Multiple Constraints and Bayesian Optimization

1. Introduction

The thermal structure of the lithosphere provides a fundamental basis for understanding the composition, rheological properties, tectonic deformation, and deep geodynamic processes of the continental lithosphere, and is of considerable significance for geothermal resource assessment, tectonic stability evaluation, and studies of geodynamic evolution [1,2]. Surface heat-flow and geothermal-gradient observations directly reflect the thermal state of the shallow crust and therefore constitute important constraints on lithospheric thermal structure [3]. However, in tectonically complex regions, heat-flow observations are spatially uneven, and their ability to characterize large-scale deep thermal structure continuously is limited by both sparse sampling and local disturbances [4,5]. Crustal temperature generally increases with depth. When magnetic minerals reach their Curie temperature, ferromagnetism gradually transitions to paramagnetism, resulting in a marked decrease in crustal magnetization. The lower boundary of the effective magnetic layer inferred from magnetic anomalies is commonly referred to as the Curie Isothermal Surface (CIS). The CIS is jointly influenced by crustal temperature, magnetic mineral composition, and the structure of the magnetic basement, and thus represents an important geophysical interface for investigating lithospheric thermal structure [5,6,7,8].

Magnetic-anomaly spectral analysis is a commonly used method for estimating the CIS. Its basic principle is to use the relationship between the slopes of the power spectrum over different wavenumber ranges and the burial depths of magnetic sources to estimate the top and centroid depths of the magnetic layer, from which the depth to the bottom of the magnetic layer is further derived [9,10,11]. However, magnetic anomalies result from the superposition of magnetic bodies with different spatial scales, burial depths, and magnetization intensities, leading to non-uniqueness in their spectral characteristics [12,13]. In practical applications, the conventional centroid method and its modified variants generally require manual parameter selection, causing the inversion results to be influenced by empirical judgment and thereby reducing their objectivity [14].

In recent years, automated optimization methods have demonstrated strong global search capabilities in geophysical inversion and spatial prediction, particularly for optimization problems involving computationally expensive objective functions, nonlinear and multimodal parameter spaces, and discrete variables [15,16,17]. To reduce the influence of subjective parameter selection in conventional spectral inversion, this study proposes an adaptive CIS inversion method based on multiple constraints and Bayesian optimization. Within the de-fractal spectral framework, a Tree-structured Parzen Estimator (TPE) is introduced to jointly optimize the fractal index, window size, and wavenumber fitting intervals for the top and centroid depths of the magnetic layer, thereby enabling adaptive searching of the spectral parameters.

2. Data

2.1. Magnetic Anomaly Data

This study uses the global gridded magnetic anomaly dataset EMAG2v3 released by the National Centers for Environmental Information (NCEI). The dataset integrates shipborne, airborne, satellite, and geomagnetic observatory measurements. For the shipborne and airborne magnetic survey data, the effects of the main geomagnetic field and external magnetic fields were removed using Comprehensive Model 4, while the long-wavelength magnetic anomaly components were supplemented using the satellite lithospheric magnetic field model MF7. EMAG2v3 has a spatial resolution of 2′. In this study, the product referenced to a uniform altitude of 4 km was used and resampled to a spatial resolution of 0.1°. The spatial distribution of magnetic anomalies in the study area is shown in Figure 1.

The magnetic anomalies in the study area exhibit distinct spatial zonation. The main part of the Ordos Basin is characterized predominantly by broad positive and negative anomaly belts, whereas the amplitudes and spatial patterns of magnetic anomalies are more complex in the surrounding tectonic regions. Previous studies have shown that the sedimentary cover of the Ordos Basin is weakly magnetic and that regional magnetic anomalies are mainly controlled by the magnetic structure of the crystalline basement, with different magnetic anomaly belts showing good spatial correspondence with the regional basement structure [18,19]. Therefore, magnetic anomalies provide prior information on spatial variations in the magnetic structure of the study area and are used in this study as one of the constraints for CIS inversion parameter optimization.

2.2. Measured Heat-Flow Data

The measured heat-flow data were obtained from the China Terrestrial Heat Flow Database [20,21,22,23,24], as shown in Figure 2. In this study, the heat-flow observations were not interpolated to generate a continuous heat-flow field; instead, they were used only as constraint points at their corresponding spatial locations for objective-function calculation and result evaluation. To make full use of the limited measured heat-flow data, a five-fold cross-validation scheme was adopted. In each fold, 80% of the heat-flow observations were used as the optimization set for objective-function calculation and Bayesian parameter search, while the remaining 20% were used as the validation set to evaluate the inversion results. After completion of the five-fold cross-validation, Bayesian optimization was rerun using all measured heat-flow observations, and the optimal parameters obtained from the full dataset were adopted.

2.3. Other Geophysical Data

To investigate the relationships between the CIS, regional crustal structure, and present-day tectonic activity, Bouguer gravity anomaly, Moho depth, topography, and earthquake data were incorporated as auxiliary datasets. The Bouguer gravity anomaly data were obtained from the global gravity field model WGM2012 [25,26]. Moho depth, upper-crustal thickness, middle-crustal thickness, and topographic data were derived from the CRUST1.0 global crustal model [27]. The earthquake dataset mainly consists of magnitude records for events occurring within the study area since 1900 [28,29,30].

3. Methods

3.1. De-Fractal Spectral Method

Crustal magnetization exhibits fractal characteristics in the horizontal direction, and its power spectrum decays with wavenumber according to a power law [31,32]. The observed power spectrum E_F(k) can be expressed as the product of the equivalent power spectrum E_R(k) under a random magnetization model and the fractal term [33]:

$$E_F\left(k\right)=E_R\left(k\right)·k^{-\alpha }$$

(1)

Where k is the radial wavenumber and α is the fractal index. Let Z_t denote the depth to the top of the magnetic layer, Z₀ the centroid depth, and Z_b the CIS depth. Z_t and Z₀ are estimated separately using the approximately linear relationships over different wavenumber ranges, for which the logarithmic spectra vary approximately linearly with k:

$$\begin{gathered} {lnE}_R\left(k\right)\approx C_t-2k{Z}_t, k\in [k_{t1},k_{t2}] \\ ln\left({\displaystyle \frac{E_R\left(k\right)}{k}}\right)\approx C_0-2k{Z}_0, k\in [k_{01},k_{02}] \end{gathered}$$

(2)

Least-squares linear fitting is performed for lnE_R(k) over the wavenumber interval [k_t1,k_t2] and for ln(E_R(k)/k) over [k₀₁,k₀₂]:

$$Z_t=-{\displaystyle \frac{s_t}{2}}, Z_0=-{\displaystyle \frac{s_0}{2}}$$

(3)

Where s_t and s₀ are the corresponding fitted slopes. The CIS depth Z_b is then calculated as:

$$Z_b=2Z_0-Z_t$$

(4)

3.2. TPE Optimizer and Sampling Strategy

In CIS inversion, different combinations of the fractal index, window size, and wavenumber fitting intervals can all affect the inversion results. This optimization problem is characterized by high dimensionality, nonlinearity, and discrete variables, making Bayesian optimization well suited for parameter searching [30]. In this study, the Tree-structured Parzen Estimator (TPE) implemented in the Optuna framework was adopted as the sampler [31]. Let x denote the parameter combination to be optimized and y the corresponding objective-function value. TPE estimates p(x|y) and thereby constructs a probabilistic model of the objective function p(y|x). Based on the results of previous trials, the samples are divided into a better-performing set and a non-optimal set according to an objective-function threshold y*, and their probability density functions are modeled separately:

$$p\left(x|y\right)=\left\{\begin{array}{c}l\left(x\right),y\le y\ast \\ g\left(x\right),y>y\ast \end{array}\right.$$

(5)

Where l(x) represents the parameter distribution corresponding to the better-performing samples, whereas g(x) represents the parameter distribution of the remaining samples. At each iteration, TPE selects the candidate parameter set that maximizes l(x)/g(x), i.e., the parameter combination that is more likely to produce a lower loss value:

$$x_{n+1}=argmax{\displaystyle \frac{l\left(x\right)}{g\left(x\right)}}$$

(6)

3.3. Construction of the Objective Function

Bayesian optimization was used to search for the fractal parameter, moving-window size, and spectral fitting wavenumber ranges. The parameter vector to be optimized is defined as:

$$x=\left(\alpha ,W,k_{t1},k_{t2},k_{01},k_{02}\right)$$

(7)

Where α is the fractal index, W is the window size, and [k_t1,k_t2] and [k₀₁,k₀₂] are the spectral fitting wavenumber ranges for the top and centroid depths of the magnetic layer, respectively. The composite objective function is defined as:

$$L={\lambda }_g{L}_g+{\lambda }_p{L}_p+{\lambda }_s{L}_s+{\lambda }_u{L}_u+{\lambda }_m{L}_m$$

(8)

Where L_g is the magnetic-anomaly prior constraint term, L_p is the fitting term for the CIS converted from heat flow, L_s is the spatial-distribution constraint term, L_u is the uncertainty constraint term, and L_m is the valid-point coverage penalty term. λ denotes the weight assigned to the corresponding constraint term. During construction of the objective function, constraints are imposed only at locations with measured heat-flow data. After the parameters for all blocks have been determined, the CIS is calculated over the entire study area.

3.3.1. Magnetic-Anomaly Prior Constraint Term L_g

To preserve the information on deep magnetic structure reflected by the magnetic anomalies, a three-class prior representing deep, intermediate, and shallow CIS depths was established on the basis of the magnetic anomalies within each block. The 33rd and 66th percentiles of magnetic anomalies across all blocks were used as classification thresholds to distinguish the deep, intermediate, and shallow classes, yielding the prior CIS depth $Z_{b,i}^{piror}$ for each block:

$$Z_b^{prior}\left\{\begin{array}{c}{Z}_b^{deep}\\ Z_b^{middle}\\ Z_b^{upper}\end{array}\right.$$

(9)

The corresponding parameters were determined from the 33rd and 66th percentiles of the CIS values converted from measured surface heat flow. For the i th constraint point, the measured surface heat flow was converted into the corresponding reference CIS depth $Z_{b,i}^{ref}$ using a conductive radiogenic heat-production model [36]. The relationship between surface heat flow and CIS depth is given by [37]:

$$q_s={\displaystyle \frac{K\left(T_c-T_s\right)}{Z_b-Z_0}}+{\displaystyle \frac{H_0{h}_r^2{e}^{-Z_b{/h}_r}}{Z_b}}+H_0{h}_r{e}^{-Z_b{/h}_r}$$

(10)

Where K is the volumetric vertical thermal conductivity of the magnetic layer and is set to 2.62 W/(m℃). T_c and T_s are the Curie temperature and surface temperature, respectively, and are set to 580℃ and 15℃. H₀ is the surface heat-production rate, set to 2.85μW/m³, and h_r is the characteristic depth of heat production, set to 13.8 km [7]. According to the difference between $Z_{b,i}^{ref}$ and the prior CIS depth $Z_{b,i}^{piror}$ of the corresponding block, the reliability weight of each observation point is defined as:

$${\mathsf{\omega }}_{\mathrm{i}}=\mathrm{m}\mathrm{a}\mathrm{x}\left[0.15,{\displaystyle \frac{1}{1+{\left(\frac{\left|Z_{b,i}^{ref}-Z_{b,i}^{piror}\right|}{5}\right)}^2}}\right]$$

(11)

For a given block, let the median CIS depth of all valid inversion points be $\widetilde{Z_b^{inv}}$, The magnetic-anomaly prior constraint term is then defined as:

$${\mathrm{L}}_{\mathrm{g}}=\left\{\begin{array}{c}max\left(0,Z_{low}^{deep}-\widetilde{Z_b^{inv}}\right)+{\eta }_{deep}\left|\widetilde{Z_b^{inv}}-Z_b^{deep}\right|, Lower CIS\hfill \\ \left|\widetilde{Z_b^{inv}}-Z_b^{middle}\right|, Middle CIS\\ max\left(0,\widetilde{Z_b^{inv}}{-Z}_{up}^{upper}\right)+{\eta }_{upper}\left|\widetilde{Z_b^{inv}}-Z_b^{upper}\right|, Upper CIS\hfill \end{array}\right.$$

(12)

Where $Z_{low}^{deep}$ and $Z_{up}^{upper}$ are the lower and upper bounds of the prior CIS depth, respectively, defined as the 5th and 95th percentiles of the heat-flow-derived reference CIS values. ${\eta }_{deep}$ and ${\eta }_{upper}$ are soft-constraint coefficients. Because there is no strict one-to-one relationship between magnetic-anomaly amplitude and CIS depth, this term is implemented as a soft penalty to reduce the weight of parameter combinations that strongly deviate from the regional magnetic-anomaly background.

3.3.2. Heat-Flow-Derived Curie Isothermal Surface Fitting Term L_p

The CIS fitting term for the i th constraint point is represented by the weighted root mean square error:

$$L_P={\left[{\displaystyle \frac{{\sum }_{i=1}^{N_v}{\omega }_i{\left(Z_{b,i}^{ref}-Z_{b,i}^{inv}\right)}^2}{{\sum }_{i=1}^{N_v}{\omega }_i}}\right]}^{1/2}$$

(13)

Where N_v is the number of valid inversion points.

3.3.3. Spatial-Distribution Constraint Term L_s

Experimental tests showed that when only pointwise RMSE is used as the objective function, the overall error may remain small while local depth relationships and spatial variations become unreasonable. Therefore, the spatial-distribution constraint was constructed from three aspects: pairwise depth differences, spatial correlation, and relative depth ordering:

$$L_s={\gamma }_1{L}_{pair}+{\gamma }_2{L}_{corr}+{\gamma }_3{L}_{sign}$$

(14)

where L_pair is the pairwise depth-difference error term, Lcorr is the spatial-correlation constraint term, L_sign is the relative depth-ordering constraint term, and γ denotes the corresponding weighting coefficients. These terms are expressed as:

$$\begin{gathered} L_{pair}={\left[{\displaystyle \frac{\sum _{\left(i,j\right)\in P}{\omega }_{ij}{\left[\left(Z_{b,i}^{cal}-Z_{b,j}^{cal}\right)-\left(Z_{b,i}^{ref}-Z_{b,j}^{ref}\right)\right]}^2}{\sum _{\left(i,j\right)\in P}{\omega }_{ij}}}\right]}^{1/2} \\ {L}_{corr}={\sigma }_{ref}max\left(\mathrm{0,1}-\rho \right) \\ {L}_{sign}=r_{mis}max\left({\sigma }_{ref},1\right) \\ {r}_{mis}={\displaystyle \frac{1}{N_r}}\sum _{\left(i,j\right)\in P}I\left[sgn\left(Z_i^{inv}-Z_j^{inv}\right)\ne sgn\left(Z_i^{ref}-Z_j^{ref}\right)\right] \end{gathered}$$

(15)

where P denotes the set of valid point pairs, ω_ij=1/(d_ij+0.05) is the weighting factor, and d_ij is the distance between two points. ρ is the Pearson correlation coefficient between the inverted and constrained values, and σ_ref is the standard deviation of the CIS constraint values and inversion values. r_mis represents the proportion of point pairs for which the relative shallow–deep or high–low relationship is reversed, and I is the indicator function.

3.3.4. Uncertainty Constraint Term L_u

The CIS uncertainty ΔZ_b at each valid point was obtained by propagating the standard errors of the spectral fitting slopes, and the uncertainty constraint term L_u was constructed as:

$$\begin{gathered} L_u={\lambda }_{u,1}\overline{△Z_b}+{\lambda }_{u,2}\overline{\mathrm{m}\mathrm{a}\mathrm{x}\left(△Z_b-\mathrm{5,0}\right)}+{\lambda }_{u,3}\overline{\mathrm{m}\mathrm{a}\mathrm{x}\left(△Z_b-\mathrm{10,0}\right)} \\ △Z_b=\sqrt{{\left(2△Z_0\right)}^2+{\left(△Z_t\right)}^2} \end{gathered}$$

(16)

Where Z_b is the mean uncertainty of all valid points and λ_u,i denotes the corresponding weighting coefficient. Equation (16) imposes an additional penalty on results with uncertainties exceeding 5 km and a stronger penalty on those exceeding 10 km.

3.3.5. Valid-Point Coverage Penalty Term L_m

Let N be the total number of constraint points. The valid-point coverage penalty term is expressed as:

$$L_m=1-{\displaystyle \frac{N_v}{N}}$$

(17)

When certain parameter combinations result in insufficient data within the moving window, spectral fitting failure, or inversion results outside the reasonable range, the number of valid points N_v decreases and L_m increases accordingly. This encourages TPE to preferentially select parameter combinations capable of covering a larger number of observation points.

4. Experimental Testing and Algorithm Validation

4.1. Determination of Objective-Function Weights and Sensitivity Analysis

A combination of nested spatial cross-validation and TPE was adopted to jointly optimize the weights of the individual constraint terms and their internal subterms in the objective function. In the outer TPE optimization, 24 candidate weight combinations were generated simultaneously within the predefined ranges and initially screened using two-fold spatial cross-validation. Based on the spatial coordinates of the reference points, the study area was divided into two spatial subregions. In each fold, one subregion was used as the training region, within which block-wise Bayesian inversion was rerun, while the reference points in the other subregion were used to evaluate predictive performance. The inner Bayesian optimization was used only to search for the fractal parameter, moving-window size, and spectral fitting wavenumber ranges, with their search ranges kept unchanged, thereby maintaining independence between weight selection and physical-parameter inversion.

After the initial screening, the five candidate weight combinations with the best overall evaluation performance were selected and, together with the initial weight combination, further evaluated using three-fold spatial cross-validation. The final objective-function weights were determined according to the combined performance of prediction error, mean uncertainty, and valid-point coverage in the test regions. Through the above joint optimization and spatial cross-validation procedure, the optimal weight combination applicable to the study area was obtained within the predefined parameter ranges, as listed in Table 1.

4.2. Objective-Function Ablation Experiment Design

A small-area synthetic spectral experiment was designed to evaluate the ability of the proposed method to recover the spatial morphology, local anomalies, and overall depth trends of the CIS. A 4°×4° synthetic region was constructed, within which a true Curie isotherm depth field $Z_b^{ture}$ containing both regional variation trends and local anomalous perturbations was generated. The corresponding magnetic anomalies were then produced from the prescribed true CIS field. In the experiment, the fractal index was set to α=1.2. Random noise perturbations were incorporated during spectral construction by adding Gaussian random noise with an amplitude equivalent to 14% of the signal standard deviation, in order to simulate disturbances caused by observational errors, local anomalies, and instability in spectral estimation during practical spectral analysis. Twenty percent of the reference points were randomly sampled from the synthetic true CIS field and used as reference depths corresponding to simulated heat-flow constraints for calculation of the Bayesian-optimization objective function.

To evaluate the effects of different levels of parameter optimization and constraint information on the inversion results, a conventional fixed-parameter method and three progressively optimized schemes were designed. In the conventional method, the fractal index, window size, and fitting wavenumber intervals were fixed at α=1, windows_size=3°, k_t1=0.03 km^-1, k_t2=0.08 km^-1, k₀₁=0.005 km^-1, k₀₂=0.055 km^-1. This parameter combination represents commonly used empirical fitting-band ranges in the conventional centroid method or modified centroid method [38]. In M1, the fractal index and window size were fixed, and only the fitting wavenumber intervals for Z_t and Z₀ were treated as parameters to be optimized. In M2, the fractal index, window size, and the two fitting wavenumber intervals were jointly optimized, with the fitting error relative to the heat-flow-derived reference CIS used as the primary optimization criterion. In M3, magnetic-anomaly priors, spatial-distribution consistency, inversion uncertainty, and valid-point coverage constraints were further introduced on the basis of M2, constituting the complete adaptive inversion method proposed in this study. The number of Bayesian-optimization trials was set to 150 for each of the three optimization schemes.

4.3. Synthetic Experiment Results

The spatial distributions of the experimental results (Figure 3) show that the conventional fixed-parameter method can recover the regional variation trend of the Curie isotherm depth, but tends to underestimate amplitudes or produce spatial offsets in areas of local shallow and deep anomalies. In comparison, the proposed method, through joint optimization of the fractal index, window size, and fitting wavenumber ranges, more effectively recovers the regional depth variations and local anomalous features of the synthetic CIS. Areas with large errors are markedly reduced, indicating that the adaptive fitting-band selection strategy improves inversion stability at different spatial locations. The overall error of the proposed method is reduced, and errors near the boundaries of local anomalies are also decreased, indicating a better balance between spectral stability and preservation of local structural features.

Table 2 presents a quantitative comparison among the conventional fixed-parameter method, the method in which only the wavenumber intervals were optimized, the full-parameter Bayesian optimization method, and the final method proposed in this study. For the conventional fixed-parameter method, the RMSE and MAE are 1.24 km and 1.03 km, respectively, while the mean relative error and median relative error are 3.12% and 2.88%, respectively. The method that optimizes only the wavenumber intervals does not yield improved results, with the RMSE increasing to 1.37 km. This indicates that, when the fractal index and window size are fixed, optimizing only the fitting wavenumber ranges for Z_t and Z₀ cannot stably improve the inversion performance.

When the fractal index, window size, and wavenumber fitting ranges are jointly included in Bayesian optimization, the RMSE of M2 decreases to 0.98 km and the MAE decreases to 0.79 km, indicating that joint multi-parameter optimization can effectively reduce errors associated with empirical parameter selection. After further introducing magnetic-anomaly priors, spatial-distribution consistency, inversion uncertainty, and valid-point coverage constraints, the RMSE and MAE of the complete method decrease to 0.92 km and 0.74 km, respectively, representing reductions of approximately 25.8% and 28.2% relative to the conventional method. Meanwhile, the mean relative error decreases from 3.12% to 2.23%, and the median relative error decreases from 2.88% to 1.88%, indicating that the proposed method outperforms the conventional fixed-parameter method in terms of both overall accuracy and error level.

5. Curie Isothermal Surface in the Ordos Basin

5.1. Block Division of the Ordos Region

The study area spans the Ordos Basin and several surrounding tectonic units, which exhibit pronounced differences in magnetic basement structure, magnetic anomaly characteristics, and deep crustal structure. Applying a uniform set of spectral parameters to the entire study area would make it difficult to accommodate the distinct magnetic-anomaly spectral characteristics of different tectonic regions. Therefore, prior to parameter optimization, the study area was divided into multiple blocks, allowing spectral parameters to be searched independently within each block. This reduces the influence of regional structural differences on the selection of uniform parameters and provides a spatial basis for the subsequent adaptive inversion.

The block division follows different criteria for the basin interior and its surrounding regions, as shown in Figure 4. In tectonically active areas outside and around the Ordos Basin, block boundaries were primarily determined according to the regional tectonic framework, major tectonic units, and the spatial distribution of large-scale faults. The interior of the Ordos Basin is relatively tectonically stable, although differences in the crystalline basement remain among different regions. Because the sedimentary cover of the basin is weakly magnetic and the regional magnetic anomalies are mainly controlled by the magnetic structure of the crystalline basement, the basin interior was further subdivided with reference to previously reported basement structures and magnetic-anomaly zoning characteristics [19]. On this basis, subsequent TPE parameter optimization was conducted independently for each block to obtain parameter combinations adapted to the magnetic-anomaly spectral characteristics of different regions.

5.2. Inversion Results and Parameter Characteristics

To evaluate the five-fold cross-validation results, the inverted CIS at each validation point was substituted into Eq. (10), and heat flow was calculated using the same thermophysical parameters and then compared with the measured heat flow. The RMSE was 35.9462±1.2174 mW/m², and the MAE was 32.4984±0.6750 mW/m². After completion of the cross-validation, Bayesian optimization was rerun using all measured heat-flow points, and the parameters obtained from full-data optimization were used to perform the final regional inversion. The resulting CIS and associated parameter distributions are shown in Figure 5.

Figure 5 shows the distributions of the CIS, optimal window size, fractal index, and inversion uncertainty in the Ordos Basin and surrounding regions, together with representative spectral fitting results for the top and centroid depths of the magnetic layer. The CIS is generally deeper within the main body of the Ordos Basin, with a local deep zone developed near the Hetao region in the north. In parts of the central and northern basin, the CIS reaches 36–40 km. In contrast, the Qilian Fold System, northeastern Tibetan Plateau, Qinling Orogen, Weihe Graben, Taihang tectonic belt, and North China Plain are characterized by relatively shallow CIS values, locally ranging from approximately 22 to 30 km. This spatial pattern indicates pronounced differences in the CIS and deep thermal state between the stable interior of the Ordos Block and the surrounding tectonic belts.

Compared with inversion using uniform parameters throughout the entire study area, block-wise adaptive inversion better preserves depth variations among different tectonic units and reduces the smoothing effect of uniform parameters across tectonic boundaries and local anomalies. Figure 5(b) shows that the window size varies among tectonic units, indicating that a fixed window cannot simultaneously account for regional-scale deep information and local spatial resolution, whereas adaptive window selection allows the computational scale to be adjusted according to the magnetic-anomaly spectral characteristics of different blocks. Figure 5(c) shows that the fractal index varies among different basement units within the Ordos Basin, and pronounced differences are also observed between the basin interior and the surrounding orogenic belts, rift zones, and the eastern North China Craton. The fractal index reflects the decay characteristics of the magnetic-anomaly power spectrum, and its spatial variation indicates significant regional differences in the distribution of magnetization intensity, the complexity of magnetic sources, and spatial correlation.

Figure 5(d) shows that the CIS uncertainty is less than 5 km over most of the study area. Uncertainties are generally lower within the main part of the Ordos Basin, indicating relatively stable magnetic-anomaly spectra and good linear characteristics within the selected wavenumber intervals. Higher uncertainties are mainly concentrated near the southwestern and southeastern tectonic boundaries of the study area. The southwestern margin contains relatively few constraint points, whereas magnetic anomalies vary strongly in the southeastern tectonic boundary region. The superposition of spectral contributions from magnetic sources at different depths and spatial scales is therefore more complex, increasing the difficulty of identifying linear fitting segments. In the representative spectra shown in Figure 5(e) and 5(f), both the low-wavenumber and medium- to high-wavenumber ranges exhibit relatively clear near-linear relationships, indicating that Bayesian optimization can identify suitable fitting intervals for estimating Z₀ and Z_t from the candidate spectral ranges.

5.3. Comparison with Previous Results

Different inversion methods differ in their treatment of magnetic-anomaly data, assumptions regarding magnetic-layer models, and parameter-selection strategies. Consequently, the resulting CIS models are not completely consistent in terms of spatial resolution, degree of smoothing, and anomaly amplitude. Figure 6 presents the CIS results obtained using four methods, and Table 3 lists the corresponding evaluation metrics. To ensure comparability among the different CIS models, all four CIS models were sampled at the same heat-flow observation locations, and surface heat flow was calculated using the same thermophysical parameters and heat-conduction model.

The RMSE and MAE obtained using the method proposed in this study are close to those of the regional inversion by Gao et al. (2015) and lower than those of Li et al. (2017) and Zhou et al. (2024). The correlation coefficient is higher than those of Gao et al. (2015) and Zhou et al. (2024) and is generally comparable to that of Li et al. (2017). These results indicate that the proposed method improves the consistency between the inversion results and the spatial variation in heat flow while maintaining a regional-scale error level comparable to previous studies.

The global CIS reference model established by Li et al. (2017) using EMAG2 data was obtained by integrating and averaging results from multiple fixed windows. The resulting model exhibits good regional continuity and can stably characterize thermal-structure features at global and continental scales. In the Ordos region, this model shows the general pattern of a relatively deep CIS within the basin and a shallower CIS in the surrounding regions, which is broadly consistent with the results of this study. The present results exhibit more pronounced lateral variations along the southwestern, southern, and eastern margins of the Ordos Basin, particularly within several tectonic transition zones. This difference is related to the block-wise parameter optimization adopted for regional tectonic units and the incorporation of measured heat-flow constraints in this study. The similarity between the two results indicates that the present inversion preserves the stability of the regional background field, whereas the local differences reflect the different spatial-scale emphases of global-scale models and detailed regional inversions.

Gao et al. (2015) obtained the Curie point depth distribution in the Ordos region using a continual-model inversion method. Their results show that the interior of the Ordos Block is dominated by relatively deep values exceeding 32 km, in good agreement with the results of this study. Both studies indicate a deeper bottom boundary of the magnetic layer beneath the basin interior and a shallower boundary in the surrounding regions. The continual-model results exhibit a relatively complete high-value zone within the main body of the basin and further reveal local differences among different basement units within the basin. These differences reflect the distinct responses of continuous-interface recovery and block-wise spectral inversion to magnetic-anomaly information at different spatial scales.

The adaptive-window method proposed by Zhou et al. (2024) selects the optimal window according to the stability of the CIS with changing window size, effectively alleviating the difficulty of simultaneously capturing deep information and local spatial resolution with a fixed window. Their results exhibit relatively abundant local variations in some tectonic boundary zones and block-junction areas, reflecting the ability of adaptive windows to respond to magnetic structures at different scales. The overall variation trends along the western and southern margins of the basin and in the eastern transition zone are generally consistent with those obtained in this study, although the local gradients of variation are larger. This difference arises because, in addition to the window size, the present method simultaneously adjusts the fractal index and the fitting wavenumber intervals for the top and centroid depths of the magnetic layer, while heat-flow and spatial-structure constraints are used to comprehensively screen different parameter combinations. The resulting CIS therefore preserves variations in tectonic transition zones while maintaining relatively good continuity between adjacent regions. Profiles of the four results along A–A′ are shown in Figure 7.

5.4. Regional Differences Between the Curie Isothermal Surface and Crustal Structure

Comparison along the A–A′ profile (Figure 7.) shows that the CIS derived in this study exhibits pronounced spatial differences relative to Moho depth and Bouguer gravity anomalies. To further examine the regional distribution of these differences across the entire study area, the relationships between the CIS and regional crustal structure were analyzed in combination with Bouguer gravity anomalies, Moho depth, and upper- and middle-crustal thicknesses, as shown in Figure 8. The Bouguer gravity anomaly generally increases from west to east. The northeastern Tibetan Plateau and the Qilian Fold System are dominated by strong negative anomalies, whereas the main part of the Ordos Basin is characterized mainly by moderate negative anomalies, which gradually transition eastward to weak negative and positive anomalies toward the Taihang Mountains and North China Plain. Correspondingly, crustal thickness also exhibits a pronounced east–west variation. Moho depth is greatest beneath the northeastern Tibetan Plateau and the Qilian Fold System, where the total crustal thickness reaches approximately 48–54 km. The crust beneath the main part of the Ordos Basin is approximately 40–45 km thick, whereas farther east beneath the Taihang Mountains and North China Plain, crustal thickness gradually decreases to approximately 30–39 km. The thicknesses of the upper and middle crust also generally decrease from west to east, locally exceeding 20 km in the western tectonic region but becoming distinctly thinner in the eastern region. These characteristics indicate that Bouguer gravity anomalies and crustal thickness show good regional-scale correspondence across the study area, jointly reflecting the structural transition from the thickened crust of the northeastern Tibetan Plateau to the thinner crust beneath the North China Plain.

Within the main part of the Ordos Basin, the CIS generally shows good spatial correspondence with crustal structure. The central and northern parts of the basin are characterized by relatively thick crust and a generally deep CIS with good spatial continuity. Toward the eastern basin and the Taihang Mountains–Shanxi Rift, crustal thickness gradually decreases and the CIS also becomes shallower overall. In contrast, a more pronounced vertical structural difference occurs near the boundary between the northeastern Tibetan Plateau and the western Ordos Basin. In this region, Moho depth generally exceeds 48 km, and the upper and middle crust are also significantly thickened, whereas the CIS is mainly located at approximately 25–32 km, markedly shallower than the Moho and close to the base of the upper-crustal depth range. This feature indicates that, although the western tectonic region has undergone substantial crustal thickening, the effective magnetic layer has not increased synchronously with total crustal thickness and remains mainly confined to the upper crust and the upper part of the middle crust. This pronounced separation between the CIS and crustal thickness provides an important basis for further analysis of the thermal state and tectonic processes within the thickened crust of this region.

6. Discussion

The results of this study show pronounced differences in CIS distribution between the main part of the Ordos Basin and its eastern and western margins, reflecting the influences of different tectonic processes, including stable cratonic evolution, plateau–craton interaction, and reworking of the eastern North China Craton, on the regional crustal thermal state and magnetic structure, as illustrated in Figure 9.

The CIS in the central part of the Ordos Basin is generally deep, locally reaching 36–40 km, and corresponds spatially to a relatively deep Moho and smooth Bouguer gravity anomalies. This region developed on Archean–Paleoproterozoic crystalline basement, and large active faults and Phanerozoic magmatic activity are relatively weak within the basin, with present-day tectonic deformation concentrated mainly along its margins. Magnetotelluric and seismic imaging results indicate that the interior of the Ordos Block is generally characterized by high resistivity and high seismic velocities, suggesting relatively stable crustal and lithospheric structures and the preservation of overall stable crustal and lithospheric characteristics [39,40]. Such conditions are favorable for maintaining a relatively low regional geothermal gradient and therefore a deep CIS. In addition, the relatively intact crystalline basement and lower-crustal structure preserved within the basin provide a material basis for a thick effective magnetic layer. Therefore, the deep CIS in the central basin mainly reflects the relatively low thermal state and thick effective magnetic layer formed during the long-term evolution of a stable craton.

The CIS beneath the eastern Ordos region and the Taihang Mountains–Shanxi Rift is distinctly shallower than that beneath the central basin. Since the Late Mesozoic, westward subduction of the Paleo-Pacific Plate and related geodynamic processes have been widely regarded as an important geodynamic background for lithospheric reworking and destruction of the eastern North China Craton [41]. Slab rollback and associated asthenospheric upwelling promoted thermochemical modification of the lithospheric mantle, resulting in thinning and replacement of the ancient lithosphere [42]. Cenozoic extension of the Shanxi Rift further facilitated regional lithospheric adjustment and changes in the deep thermal state [43], increasing the regional geothermal gradient and causing relative shallowing of the bottom boundary of the effective magnetic layer. Meanwhile, lithospheric thinning reduced crustal thickness, resulting in overall Moho uplift and a corresponding increase in Bouguer gravity anomalies. Therefore, the shallow CIS in the eastern study area reflects the combined effects of lithospheric thinning of the eastern North China Craton, enhanced deep thermal conditions, and Cenozoic extension under the geodynamic background of Paleo-Pacific subduction.

The Qinling Orogen records multiple stages of convergence and collision between the North China Craton and the Yangtze Craton. Lithospheric roots, crustal detachment structures, and pronounced lateral velocity variations revealed by deep seismic reflection and seismic tomography indicate that the ancient collision-related suture zone continues to control crustal deformation and deep structure along the southern margin [44,45]. Seismic imaging has identified an upper-mantle low-velocity anomaly extending eastward from the northeastern Tibetan Plateau between the Ordos Block and the Yangtze Craton. It has been suggested that the rigid lithospheric roots beneath the Ordos Block and Yangtze Craton jointly restrict the eastward migration of asthenospheric material, forcing it to undergo channelized flow along the Qinling tectonic belt [46,47]. The relatively high-velocity structure in the middle and lower crust of the western Qinling may impede further eastward migration of low-velocity crustal material, whereas low-velocity anomalies in the lithospheric mantle beneath the Weihe Graben and eastern Qinling reflect thermochemical modification of the upper mantle associated with asthenospheric activity [48]. Three-dimensional gravity and magnetic inversion of the Weihe Basin indicates that regional heat flow does not show an obvious correspondence with Moho uplift, and that Yanshan granitic bodies buried beneath the basin sediments and their radiogenic heat production may constitute an important crustal heat source [49]. Differences in basement composition and magnetic mineral content on the two sides of the Weihe Fault also influence the position of the bottom boundary of the effective magnetic layer. Therefore, the relatively shallow CIS along the southern margin of the Ordos Basin may be jointly influenced by inherited paleocollisional structures, channelized asthenospheric flow and upper-mantle modification, tectonic activity of the Weihe Graben, and radiogenic heat production from granitic crust.

The northern Ordos Basin and the southern margin of the Hetao region form a local deep-CIS zone, with values reaching 36–40 km in some areas. High-resolution seismic imaging indicates that the northern Ordos region is generally characterized by relatively high seismic velocities and weak lateral variations, suggesting that relatively stable cratonic crustal and lithospheric structures are still preserved in this area [50,51]. Recent ambient-noise imaging further reveals a NE-trending high-velocity zone in the middle crust of the central–northern Ordos region that corresponds spatially to the magnetic anomalies and represents mafic–ultramafic basement structures formed during the amalgamation of ancient microcontinental blocks [52]. The relatively low thermal state, together with the well-preserved and strongly magnetic crystalline basement, may favor the development of a thick effective magnetic layer, consistent with the deep CIS obtained in the northern part of the study area. Receiver-function and seismic imaging results show that the upper crust beneath the Hetao–Yinchuan rift system has undergone extension, whereas the middle and lower crust exhibit local thickening, and residual crustal roots formed during Late Mesozoic intracontinental orogenesis are still preserved beneath parts of the rift [53,54]. Evidence for localized lithospheric thinning and mantle thermal modification beneath the Hetao Rift–Yinshan tectonic belt indicates that the northern margin has undergone varying degrees of reactivation [55]. Therefore, the deep CIS beneath the northern basin and southern Hetao margin may be related to the stable magnetic basement and residual thickened crust, whereas local variations near the rift axis and tectonic boundaries may be associated with mantle thermal perturbations and lithospheric modification.

The CIS beneath the western Ordos region, northeastern Tibetan Plateau, and Qilian Fold System is markedly shallower than the Moho. In this region, the Moho locally reaches depths greater than 48–50 km and corresponds to pronounced negative Bouguer gravity anomalies, whereas the CIS is mainly 25–32 km deep, indicating a substantial difference between total crustal thickness and effective magnetic-layer thickness. The northeastern Tibetan Plateau continues to expand toward the east and northeast. Under the obstruction of the rigid Ordos Block, the left-lateral strike-slip motion along the Haiyuan Fault gradually weakens eastward and is transformed into thrust shortening and tectonic uplift in the Liupanshan region. Crustal shortening, thrust stacking, and middle- to lower-crustal thickening jointly contribute to Moho deepening and lower Bouguer gravity anomalies in the western region [51,56]. Seismic imaging has identified low shear-wave velocities and anisotropic structures beneath the northeastern Tibetan Plateau and adjacent Liupanshan region, suggesting that ductile flow in the middle and lower crust contributes to outward plateau expansion and crustal thickening [40,57]. Other studies, however, indicate pronounced lateral heterogeneity in crustal deformation in this region. Thickened high-velocity lower-crustal structures are present beneath both Liupanshan and the southwestern margin of the Ordos Block, suggesting that crustal thickening along the western margin may be accommodated jointly by upper-crustal thrust shortening, middle- to lower-crustal stacking, and localized ductile deformation rather than by a single continuous crustal-flow mechanism [51,58]. Deformation of the middle and lower crust, ductile shearing, crust–mantle interaction, and localized thermal perturbations within the boundary zone may jointly enhance the thermal state of the middle and lower crust, preventing the bottom of the effective magnetic layer from deepening synchronously with Moho deepening caused by tectonic crustal thickening.

Figure 10 shows the spatial distributions of surface heat flow converted from the inverted CIS using the heat-conduction model and earthquakes with M≥3 in the study area since 1900. Seismicity is generally concentrated around the margins of the Ordos Block and along major active tectonic belts, whereas relatively few earthquakes occur within the main body of the basin. This corresponds to the deeper CIS, lower converted heat flow, and weaker seismic activity in the central and northern parts of the basin. In contrast, numerous earthquakes are distributed along the western and southwestern margins of the basin, particularly along the northeastern Tibetan Plateau and the Qilian tectonic belt, showing good spatial correspondence with the shallow CIS and relatively high converted heat flow. In particular, along the western margin of the Ordos Basin, Moho depth can exceed 48–50 km, whereas the CIS is mainly 25–32 km, showing a pronounced separation between crustal thickness and the bottom boundary of the effective magnetic layer. These characteristics indicate that, during the interaction between the continued expansion of the northeastern Tibetan Plateau and the rigid Ordos Block, crustal thickening did not result in synchronous deepening of the effective magnetic layer. The northeastern Tibetan Plateau and the Qilian Fold System are simultaneously characterized by relatively high CIS-derived heat flow and concentrated seismic activity on the centennial timescale, in marked contrast to the lower heat flow and weaker seismicity within the main body of the Ordos Basin.

7. Conclusions

This study addresses the dependence of the fractal index, window size, and spectral fitting intervals on manual selection in conventional magnetic-anomaly spectral methods by developing a Bayesian adaptive inversion method under multiple constraints and applying it to CIS inversion in the Ordos Basin and surrounding regions. The main conclusions are as follows:

(1) Based on de-fractal spectral analysis, the fractal index, moving-window size, and fitting wavenumber intervals for the top and centroid depths of the magnetic layer were incorporated into a unified optimization parameter space and jointly searched using TPE, thereby replacing empirical manual parameter selection with adaptive joint determination. Surface heat flow, magnetic-anomaly priors, spatial distribution, uncertainty, and valid-point coverage were further introduced as constraints during parameter optimization, transforming parameter selection from single spectral fitting into comprehensive optimization constrained by multiple sources of information and allowing adaptation to the differences in magnetic-anomaly spectra among different tectonic units. The weight-sensitivity analysis and spatial cross-validation results show that the different constraint terms provide a reasonable balance among inversion accuracy, spatial stability, and valid-data coverage.

(2) Synthetic experiments show that optimizing only a subset of the spectral parameters does not consistently improve the inversion results, whereas joint optimization of the fractal index, window size, and fitting wavenumber intervals can significantly reduce errors associated with parameter selection. With the further inclusion of multiple constraints, the recovery of regional CIS variations and local anomalies is further improved. The RMSE and MAE of the proposed method are 0.92 km and 0.74 km, respectively, representing reductions of 25.8% and 28.2% relative to the conventional fixed-parameter method. The mean and median relative errors also decrease simultaneously, indicating that joint multi-parameter optimization and multi-constraint screening improve both the accuracy and stability of spectral inversion.

(3) The CIS in the Ordos Basin and surrounding regions is generally approximately 22–40 km deep. Relatively deep values occur in the central and northern parts of the basin, reaching 36–40 km in some areas, whereas the northeastern Tibetan Plateau, Qilian Fold System, Qinling–Weihe tectonic belt, and Taihang Mountains–Shanxi Rift are generally characterized by shallower values of 22–30 km. The relatively deep CIS beneath the main body of the basin corresponds to the well-preserved crystalline basement, relatively stable crustal structure, and lower thermal state. The shallow CIS in the eastern study area shows good spatial correspondence with lithospheric thinning of the eastern North China Craton and the extensional setting of the rift system. In the western study area, the northeastern Tibetan Plateau and Qilian Fold System exhibit a more pronounced vertical structural contrast: the Moho reaches depths of 48–50 km, whereas the CIS is mainly 25–32 km, indicating that tectonically thickened crust has not developed an effective magnetic layer that increases synchronously with total crustal thickness. This region is also characterized by relatively high heat flow and concentrated seismic activity, indicating a pronounced enhancement of the thermal state of the middle and lower crust during crustal thickening along the western margin of the Ordos Basin.