Terrestrial heat ecosystem service coupling value evaluation method

Abstract

The invention discloses a geothermal ecosystem service coupling value evaluation method, and relates to the technical field of big data analysis and processing, and the method comprises the steps: obtaining the five indexes of geothermal output, underground water level, vegetation coverage, tourism income and carbon emission reduction at a plurality of moments of a target region, determining the edge distribution model of each monitoring time sequence, and obtaining the value of each monitoring time sequence; forming an edge distribution model set; each record in the basic monitoring data set is converted into a five-dimensional sample containing five types of index standardized values through a cumulative distribution function of the corresponding edge distribution model, a posterior sample set is obtained, and a Copula-Bayesian mixture model is constructed; and utilizing Monte Carlo simulation to generate five-dimensional uniform random numbers, calculating coupling value expectation and risk indexes based on the comprehensive value sample set, and obtaining a geothermal ecosystem service coupling value evaluation result. According to the method, the accuracy, the robustness and the interpretability of geothermal-ecological coupling relation evaluation are remarkably improved.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of big data analysis and processing, and particularly relates to a geothermal ecosystem service coupling value evaluation method. BACKGROUND

[0002] With the continuous expansion of geothermal development, its influence on groundwater systems, vegetation coverage, soil temperature field and regional ecological service capacity gradually attracts the attention of researchers and managers. In the existing technical system, the geological characteristics of the reservoir, the temperature distribution of the thermal reservoir or the geothermal production capacity are usually studied by geothermal resource evaluation methods, and the values of ecological elements such as carbon sink, water conservation, tourism and cultural services are studied by ecosystem service value evaluation methods. However, most studies often deal with geothermal development and ecological service changes separately, and lack systematic measurement means for the complex coupling relationship between the two. Due to the significant volatility of indicators such as geothermal production, groundwater level, vegetation coverage, tourism income and carbon emission reduction in time, and the significant nonlinear correlation characteristics in change trend, extreme event response, etc., the existing technology is difficult to accurately depict the linkage mechanism between the geothermal system and the ecological system. Many existing methods rely on linear regression, Pearson correlation coefficient, and multi-index analytic hierarchy process to comprehensively evaluate the relationship between geothermal development and ecosystem services. In these methods, the interaction between indicators is often simplified as a linear relationship, ignoring the fact that when a certain indicator is at an extreme high value or an extreme low value, the response of other indicators may show mutation or strong coupling characteristics. For example, when the groundwater level drops rapidly, the vegetation coverage may decrease significantly faster, and the tourism income may decrease at the same time due to the degradation of the landscape, but the linear correlation coefficient is difficult to capture this highly synchronized tail linkage. In addition, traditional multi-index comprehensive evaluation methods usually rely on human-set weights, and it is difficult to maintain the stability of the results when the data shows strong fluctuations. SUMMARY

[0003] The purpose of the present application is to provide a geothermal ecosystem service coupling value evaluation method, which realizes the dynamic depiction of the coupling relationship between geothermal and ecological systems, and can accurately handle the nonlinear dependence and tail linkage characteristics between indicators. Through the construction of edge distribution model set and standardized sample set, the monitoring data of different dimensions are uniformly processed in the same probability space. Through the Markov chain Monte Carlo method, the dependence parameters are continuously updated, so that the model can remain adaptive with new monitoring data. Through Monte Carlo simulation, a set of comprehensive value samples is generated, and the coupling value expectation and risk indicators are calculated, so that the evaluation results reflect both the long-term average income level and the extreme scenario risk. The present application significantly improves the accuracy, robustness and interpretability of the evaluation of geothermal-ecological coupling relationship, and provides a reliable quantitative basis for geothermal development strategy formulation and ecosystem service management.

[0004] To address the aforementioned technical problems, this invention provides a method for assessing the coupled value of geothermal ecosystem services, comprising the following steps: acquiring five types of indicators—geothermal output, groundwater level, vegetation cover, tourism revenue, and carbon emission reduction—at multiple times in a target area; combining these five types of indicators at the same time point to form a basic monitoring dataset; constructing five types of monitoring time series based on the basic monitoring dataset; performing parameter estimation and goodness-of-fit tests on each monitoring time series to determine the marginal distribution model for each monitoring time series, forming a set of marginal distribution models; and converting each record in the basic monitoring dataset into a five-dimensional sample containing standardized values ​​of the five types of indicators using the cumulative distribution function of the corresponding marginal distribution model, thus obtaining... A standardized sample set is used to select the target Copula structure model from the pre-defined Copula structure candidate set. Based on the standardized sample set, the dependency parameters are updated Bayesianly using the Markov chain Monte Carlo method to obtain the posterior sample set and construct a Copula-Bayesian mixture model. Five-dimensional uniform random numbers are generated using Monte Carlo simulation and input into the Copula-Bayesian mixture model to obtain standardized joint samples. Five types of index samples are recovered through the inverse function of the marginal distribution model set. The comprehensive value sample is calculated according to the pre-defined benefit-cost rule. Based on the comprehensive value sample set, the expected coupling value and risk index are calculated to obtain the geothermal ecosystem service coupling value assessment results.

[0005] Furthermore, the acquisition of the basic monitoring dataset includes: setting up at least 5 geothermal production monitoring points, at least 5 groundwater level monitoring points, and at least 5 vegetation coverage monitoring sample areas within the target area; setting up at least 1 tourism revenue metering point at the main tourist entrance; setting up at least 1 carbon emission reduction calculation point at the entrance and exit of the geothermal energy supply facility; continuously collecting geothermal production monitoring values, groundwater level monitoring values, vegetation coverage monitoring values, tourism revenue monitoring values, and carbon emission reduction monitoring values ​​for at least 365 days at a time interval of 1 day; and combining the above 5 monitoring values ​​collected on the same natural day into 1 record to form a basic monitoring dataset containing at least 365 records.

[0006] Furthermore, in constructing the set of marginal distribution models for geothermal ecosystem services coupling, candidate marginal distribution models include normal distribution models, log-normal distribution models, and gamma distribution models. For each monitoring time series, maximum likelihood estimation is performed on the normal distribution model, log-normal distribution model, and gamma distribution model respectively to obtain the corresponding maximum likelihood objective function value. The Kolmogorov-Smirnov test is used to obtain the goodness-of-fit statistic. Among the normal distribution model, log-normal distribution model, and gamma distribution model, the model with the largest maximum likelihood objective function value and the smallest Kolmogorov-Smirnov test statistic is selected as the marginal distribution model for that monitoring time series.

[0007] Furthermore, in the process of constructing the standardized sample set, the geothermal production monitoring value of each record in the geothermal ecosystem service coupled basic monitoring dataset is input into the cumulative distribution function of the geothermal production marginal distribution model to obtain the standardized geothermal production value; the groundwater level monitoring value is input into the cumulative distribution function of the groundwater level marginal distribution model to obtain the standardized groundwater level value; the vegetation cover monitoring value is input into the cumulative distribution function of the vegetation cover marginal distribution model to obtain the standardized vegetation cover value; the tourism revenue monitoring value is input into the cumulative distribution function of the tourism revenue marginal distribution model to obtain the standardized tourism revenue value; and the carbon emission reduction monitoring value is input into the cumulative distribution function of the carbon emission reduction marginal distribution model to obtain the standardized carbon emission reduction value. The standardized sample vectors are constructed in the order of geothermal production, groundwater level, vegetation cover, tourism revenue, and carbon emission reduction, and all five-dimensional standardized sample vectors are stored in the order of records to form a five-dimensional standardized sample set.

[0008] Furthermore, in the process of selecting the target Copula structure model, the preset candidate Copula structure set includes at least the Clayton Copula structure model, the Gumbel Copula structure model, and the Gaussian Copula structure model. For each Copula structure model, the range of dependent parameter values ​​is set to a closed interval from 0.1 to 10.0. Within this closed interval, 10 equally spaced initial parameter values ​​are selected. For each initial parameter value, the log-likelihood value of all samples is calculated using the corresponding Copula structure model and the five-dimensional standardized sample set, and the sum is used to obtain the log-likelihood objective function value. The maximum value of the log-likelihood objective function value of all initial parameter values ​​corresponding to the same Copula structure model is taken as the initial evaluation value of the Copula structure model. Among the Clayton Copula structure model, the Gumbel Copula structure model, and the Gaussian Copula structure model, the Copula structure model with the largest initial evaluation value is selected as the target Copula structure model.

[0009] Furthermore, the method for constructing the prior distribution of the dependent parameters of the target Copula structural model includes: setting the initial value point of the parameter that maximizes the initial evaluation value of the target Copula structural model as the center value of the prior distribution of the dependent parameters; setting 10% of the length of the range of dependent parameter values ​​as the standard deviation of the prior distribution of the dependent parameters; constructing a normal prior distribution with the center value and the standard deviation as parameters; and maintaining the dependent parameter values ​​within a closed interval of 0.1 to 10.0 by a limiting method within the interval where the dependent parameter values ​​are less than 0.1 or greater than 10.0.

[0010] Furthermore, the Markov chain Monte Carlo method performs a multi-step Bayesian update of the dependency parameters, which includes: setting the total number of iterations of the Markov chain to 10,000 steps, defining the first 2,000 steps as the burning stage, using the initial parameter value as the current value of the dependency parameter in the first step, and starting from the second step, in each iteration, using a normally distributed random number with a mean of 0.0 and a standard deviation of 0.1 as a perturbation value, adding it to the current value of the dependency parameter from the previous step to obtain a candidate value of the dependency parameter, setting the candidate value of the dependency parameter less than 0.1 to 0.1, and setting the candidate value of the dependency parameter greater than 10.0 to 10.0, and using the candidate value of the dependency parameter and the five-dimensional standardized sample set... The log-likelihood values ​​of all samples are calculated and summed using the target Copula structural model. This sum of log-likelihoods is then added to the log density values ​​of the candidate dependent parameters under the prior distribution of the dependent parameters to obtain the log-posterior objective function value corresponding to the candidate dependent parameter value. The log-posterior objective function value obtained in the same way from the previous step of the dependent parameter is used for comparison. The acceptance probability is calculated based on the difference between the two values. Uniformly distributed random numbers with values ​​between 0 and 1 are generated. The current value of the dependent parameter in the current iteration step is determined based on the acceptance probability and the uniformly distributed random numbers. After 10,000 iterations, the current values ​​of the dependent parameters other than those in the burning stage are used to form the dependent parameter posterior sample set.

[0011] Furthermore, the geothermal ecosystem service coupling value assessment process also includes: collecting geothermal production monitoring values, groundwater level monitoring values, vegetation coverage monitoring values, tourism revenue monitoring values, and carbon emission reduction monitoring values ​​in the second and subsequent batches under the same conditions as the monitoring point layout and time intervals; converting each batch of monitoring values ​​into a new five-dimensional standardized sample set; using the arithmetic mean of the dependent parameter posterior sample set as the new dependent parameter prior distribution center value; using the sample standard deviation of the dependent parameter posterior sample set as the new dependent parameter prior distribution standard deviation; constructing a new dependent parameter prior distribution; and using the Markov chain Monte Carlo method to perform dependent parameter Bayesian update on the new five-dimensional standardized sample set to obtain the updated dependent parameter posterior sample set, which is then combined with the target Copula structural model to form a target Copula-Bayes hybrid model with dynamic dependent parameters.

[0012] Furthermore, during the Monte Carlo simulation, the number of simulation rounds was set to 10,000. In each round of simulation, a random number generation algorithm was used to generate five independent uniformly distributed random numbers with values ​​ranging from 0 to 1. These five independent uniformly distributed random numbers were input into the Copula-Bayes mixture model. Under the joint distribution constraints of the target Copula structure model, a set of five-dimensional standardized joint samples was obtained. The five components of the five-dimensional standardized joint samples were then input into the inverse functions of the geothermal production edge distribution model, the groundwater level edge distribution model, the vegetation cover edge distribution model, the tourism revenue edge distribution model, and the carbon emission reduction edge distribution model, respectively, to obtain geothermal production samples, groundwater level samples, vegetation cover samples, tourism revenue samples, and carbon emission reduction samples.

[0013] Furthermore, the calculation method for the comprehensive value sample of geothermal ecosystem services includes: multiplying the geothermal production sample by the unit price of geothermal energy revenue to obtain geothermal energy revenue; multiplying the tourism revenue sample by the tourism industry surcharge to obtain comprehensive tourism revenue; multiplying the carbon emission reduction sample by the unit price of carbon emission reduction revenue to obtain emission reduction revenue; comparing the groundwater level sample with the groundwater level safety threshold; when the groundwater level sample is lower than the groundwater level safety threshold, calculating the ground subsidence control cost by multiplying the difference between the groundwater level safety threshold and the groundwater level sample by the unit cost of ground subsidence control; when the groundwater level sample is higher than or equal to the groundwater level safety threshold, setting the ground subsidence control cost to 0; multiplying the vegetation cover sample by the ecosystem service gain coefficient to obtain ecosystem service gain; and adding the geothermal energy revenue, comprehensive tourism revenue, emission reduction revenue, and ecosystem service gain, then subtracting the ground subsidence control cost to obtain the comprehensive value sample. The unit price of geothermal energy revenue is 30.0 monetary units per megawatt-hour, the tourism industry surcharge is 1.2, the unit price of carbon emission reduction revenue is 60.0 monetary units per ton, the groundwater level safety threshold is the arithmetic mean of groundwater level monitoring values ​​during the reference period minus 0.5 meters, and the ecosystem service gain coefficient is the ecosystem service value corresponding to 10.0 monetary units per unit of vegetation cover. In a sample set containing 10,000 comprehensive value samples, the comprehensive value samples are sorted in ascending order. The sample value with position number 500 after sorting is used as the first risk value indicator, the arithmetic mean of comprehensive value samples with position numbers 1 to 500 after sorting is used as the second risk value indicator, and the arithmetic mean of the sample set is used as the expected comprehensive value indicator. The first risk value indicator, the second risk value indicator, and the expected comprehensive value indicator are combined to form the result of the geothermal ecosystem service coupling value assessment.

[0014] This invention provides a method for assessing the coupled value of geothermal ecosystem services, offering the following advantages: First, it comprehensively reflects the nonlinear dependency structure and tail-synchronous fluctuation characteristics of multiple indicators over time, maintaining stable identification of the coupled state under complex environmental changes. Second, it uses a standardized sample set to unify the dimensions of five types of indicators, enabling linked analysis of different ecological and energy indicators within the same probability space, avoiding error accumulation caused by dimensional differences in traditional methods. Third, through a Bayesian update strategy based on Markov chain Monte Carlo methods, it dynamically adjusts dependency parameters when new monitoring data is input, allowing the joint distribution model to continuously adapt to various natural and engineering disturbances such as seasonal changes, load changes, and groundwater recharge changes, thereby improving the model's reliability in long-term assessments. Fourth, it utilizes Monte Carlo simulation to generate a multi-scenario comprehensive value sample set, enabling the calculation of coupled value expectation and risk indicators in the same process, ensuring that the assessment results simultaneously consider average levels and extreme scenarios, providing quantitative basis for the selection of geothermal development strategies and the formulation of ecological protection measures. Unlike traditional methods that typically provide only a single evaluation value, this invention not only reveals the comprehensive value of the coupled system under normal conditions but also accurately measures the degree of risk exposure under adverse scenarios, thus possessing greater practical value in the integrated management of complex energy-ecosystems. Furthermore, the process of this invention is scalable, allowing for flexible adjustment of the edge distribution model type and risk quantile settings based on monitoring conditions in different regions. This enables the method to adapt to diverse geothermal utilization patterns and ecosystem structures, thereby maintaining robustness and applicability across various application scenarios. Attached Figure Description

[0015] Figure 1 The present invention provides a pairwise correlation heatmap of the raw data of five types of monitoring indicators provided in this embodiment. Figure 2 A schematic diagram of the joint probability density surface of geothermal output and tourism revenue in physical space, provided for an embodiment of the present invention; Figure 3 This is a schematic diagram of a standardized joint distribution structure constructed using the GumbelCopula structural model, provided for an embodiment of the present invention. Detailed Implementation

[0016] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0017] A method for assessing the coupled value of geothermal ecosystem services includes the following steps: acquiring five types of indicators—geothermal output, groundwater level, vegetation cover, tourism revenue, and carbon emission reduction—at multiple times in a target area; combining these five types of indicators at the same time to form a basic monitoring dataset; constructing five types of monitoring time series based on the basic monitoring dataset; performing parameter estimation and goodness-of-fit tests on each monitoring time series to determine the marginal distribution model for each monitoring time series, forming a set of marginal distribution models; and converting each record in the basic monitoring dataset into a five-dimensional sample containing standardized values ​​of the five types of indicators using the cumulative distribution function of the corresponding marginal distribution model, obtaining a standardized sample set. The target Copula structure model is selected from the pre-set Copula structure candidate set, and its dependency parameters are updated Bayesianly using the Markov chain Monte Carlo method based on the standardized sample set to obtain the posterior sample set and construct the Copula-Bayes mixture model. Five-dimensional uniform random numbers are generated using Monte Carlo simulation and input into the Copula-Bayes mixture model to obtain standardized joint samples. The five types of index samples are recovered through the inverse function of the marginal distribution model set. The comprehensive value sample is calculated according to the pre-set benefit-cost rule. Based on the comprehensive value sample set, the expected coupling value and risk index are calculated to obtain the evaluation result of the coupling value of geothermal ecosystem services.

[0018] In one implementation, the geothermal ecosystem service coupling value assessment method was applied to a geothermal heating demonstration park. The park covers approximately 10 square kilometers and includes geothermal well sites and heat exchange stations, as well as residential areas, landscaped green spaces, and hot spring tourism facilities. To ensure the stability of subsequent statistical fitting and the integrity of the basic monitoring dataset, this implementation uses a continuous 365-day monitoring period with a 1-day time resolution to simultaneously monitor five indicators: geothermal production, groundwater level, vegetation cover, tourism revenue, and carbon emission reduction.

[0019] In this implementation, monitoring points are first deployed within the target area, and the monitoring method is defined. Geothermal production monitoring points are located at the junction of all production wells, equipped with digital flow meters and temperature sensors. The system automatically accumulates the daily hot water flow rate at 15-minute sampling intervals and converts it into equivalent heat energy. The accumulated equivalent heat energy from 00:00 to 24:00 each day is used as the geothermal production monitoring value for that day, measured in megawatt-hours (MWh). Five groundwater level monitoring points are deployed along the perimeter of the geothermal well site and downstream. A pressure level gauge is installed in each monitoring well. The water level reading at 12:00 local time is selected as the groundwater level monitoring value for that well each day, measured in meters. The arithmetic mean of the water level monitoring values ​​from the five monitoring wells is then calculated and used as the groundwater level monitoring value for that day. Five vegetation cover monitoring sample areas were set up in the target area according to land use type. Aerial imaging was conducted daily at 10:00 AM using drones, with a resolution of 0.1 meters. Image processing software was used to calculate the ratio of vegetation pixels to total pixels in each sample area. The arithmetic mean of the ratios from the five sample areas was calculated to obtain the daily vegetation cover monitoring value, which ranges from 0 to 1. The tourism revenue metering point was set at the ticket system at the entrance of the hot spring resort. The ticket system recorded daily ticket revenue and secondary consumption settlement data within the park. The sum of these two figures was used as the daily tourism revenue monitoring value, expressed in monetary units. The carbon emission reduction calculation point was set up inside the geothermal heat exchange station. Based on the daily standard coal consumption replaced by the geothermal heating system, the carbon emission reduction was calculated using an energy balance program according to the carbon dioxide emission coefficient corresponding to each ton of standard coal, yielding the daily carbon emission reduction monitoring value, expressed in tons.

[0020] refer to Figure 1This chart comprehensively reveals the linear strength and direction of pairwise interactions among five key variables: geothermal production, groundwater level, vegetation cover, tourism revenue, and carbon emission reduction. The chart uses a 5x5 matrix, with each row and column corresponding to one of the five indicators. Each cell in the matrix is ​​filled with different shades of color and labeled with a specific Pearson correlation coefficient value, ranging from -1 to +1. Darker colors indicate a stronger absolute correlation, meaning a closer connection between the variables; lighter colors indicate a correlation closer to zero, meaning the variables tend to be linearly independent. Observing the main diagonal of the matrix reveals a series of dark squares with positive values ​​of +1, representing perfect correlation between each indicator and itself—a statistically inevitable result and forming the baseline reference for the heatmap. The area outside the main diagonal shows the interaction relationships between different indicators. For example, in the cell where geothermal production and carbon emission reduction intersect, a darker, higher-valued square is observed, indicating a significant positive correlation between the two. From a technical perspective, this is because increased geothermal production directly replaces more fossil fuel consumption, thus linearly increasing carbon emission reductions. This strong correlation provides data support for setting high-intensity dependency parameters in the model. Similarly, a positive correlation coefficient is shown at the intersection of geothermal production and tourism revenue, reflecting the driving effect of geothermal resource development on related industries such as hot spring tourism. Although the correlation strength may be slightly lower than that of carbon emission reductions, it is sufficient to prove that the two cannot be treated as independent variables. Furthermore, heat maps may also contain negative or weak correlations. For example, groundwater levels may show a negative correlation with other indicators, meaning that excessive increases in geothermal production may put downward pressure on groundwater levels. Although reinjection technology can mitigate this trend, the negative correlation in the statistical data can keenly capture this potential ecological constraint. Vegetation cover and other industrial indicators may appear in lighter colors, suggesting that their changes are more influenced by seasons and climate, and their direct linear relationship with geothermal production is weaker. However, this does not rule out the existence of a non-linear dependency structure. Figure 1 The heatmap analysis in this embodiment strongly demonstrates the necessity of using a multidimensional joint distribution evaluation method: since there are a large number of non-zero correlation coefficients with varying strengths in the matrix, if these five types of indicators are simply treated as independent variables and superimposed for calculation, the organic connection within the system will inevitably be severed, leading to a systematic bias in the evaluation results.

[0021] Under the above monitoring configuration, one geothermal production monitoring value, one groundwater level monitoring value, one vegetation coverage monitoring value, one tourism revenue monitoring value, and one carbon emission reduction monitoring value can be obtained for each day of the entire monitoring period. To ensure that each record in the subsequent statistical analysis corresponds to the same physical time point, this implementation combines the five monitoring values ​​obtained within the same natural day into one record in a fixed order: date, geothermal production monitoring value, groundwater level monitoring value, vegetation coverage monitoring value, tourism revenue monitoring value, and carbon emission reduction monitoring value. For example, on the 10th day of monitoring, if the geothermal production monitoring value is 520.0, the groundwater level monitoring value is 23.4, the vegetation coverage monitoring value is 0.62, the tourism revenue monitoring value is 350,000.0, and the carbon emission reduction monitoring value is 145.0, then the corresponding record for that day would be "Day 10, 520.0, 23.4, 0.62, 350,000.0, 145.0". By combining 365 days of data in this way, a basic monitoring dataset containing 365 records can be formed. By combining five types of indicators at the same time into one record, it can be ensured that the correlation between the five types of indicators is based on the true time-series correspondence when constructing the joint distribution. This maintains the coupling structure as much as possible in a statistical sense, and can more accurately reflect the synchronous change characteristics of geothermal systems and ecosystem services compared to processing individual samples of each indicator separately.

[0022] After the basic monitoring dataset was constructed, five types of monitoring time series were built based on it. Specifically, all geothermal production monitoring values ​​from the 365 records were arranged in chronological order to obtain a geothermal production monitoring time series; all groundwater level monitoring values ​​were arranged in the same chronological order to obtain a groundwater level monitoring time series; and vegetation cover monitoring time series, tourism revenue monitoring time series, and carbon emission reduction monitoring time series were obtained in the same way. Each monitoring time series contains 365 monitoring values ​​sorted by time. This chronological arrangement allows the long-term statistical characteristics, seasonal variations, and frequency of extreme events for each indicator to be fully captured within the same series, avoiding the loss of temporal structure information that might be caused by random shuffling.

[0023] When constructing marginal distribution models for each monitoring time series, this implementation method selects the normal distribution model, log-normal distribution model, and gamma distribution model as candidate marginal distribution models for each monitoring time series. The normal distribution model is chosen because, under stable operating conditions, indicators such as geothermal output and groundwater level tend to fluctuate around a certain central value, and the probability of deviating from the central value is approximately symmetrical in the positive and negative directions. The normal distribution model can effectively describe this phenomenon of "high probability near the mean and low probability far from the mean". The log-normal distribution model is chosen because tourism revenue and carbon emission reduction are generally non-negative values, and there may be a few maxima in actual monitoring, such as a significant increase in heating demand during peak tourist seasons or abnormally cold weather. The log-normal distribution, by performing a logarithmic transformation on the data before applying the normal distribution, can naturally handle the "right skewness" and "long tail" characteristics. The gamma distribution model was chosen to improve the ability to characterize indicators that are non-negative and have a strong tendency to concentrate near low values, while having a softer long tail in high value regions. For example, when geothermal production has not yet reached full capacity or is limited by supply and demand scheduling, there will be a large number of low and medium production values ​​and a small number of high production peaks in the monitoring time series. The gamma distribution has good fitting flexibility in this scenario.

[0024] For each monitoring time series, the parameters of three candidate marginal distribution models are estimated using the maximum likelihood estimation method. Taking the geothermal production monitoring time series as an example, this series contains 365 monitoring values, with the values ​​for the first 5 days being 480.0, 510.0, 495.0, 530.0, and 505.0, respectively. For the normal distribution model, the arithmetic mean of the 365 monitoring values ​​is first calculated. For example, assuming the arithmetic sum of the 365 monitoring values ​​is 182500.0, the arithmetic mean is 182500.0 divided by 365, yielding approximately 500.0. Then, the square of the difference between each monitored value and the arithmetic mean is calculated, and the 365 squared values ​​are summed. Assuming the sum of squares is 925000.0, 925000.0 is divided by 365 to get approximately 2534.25. The square root of 2534.25 is then taken to get approximately 50.34. 500.0 and 50.34 are used as the central location parameter and dispersion parameter of the normal distribution model, respectively. For the log-normal distribution model, the natural logarithm is first taken for each of the 365 monitoring values ​​in the geothermal production monitoring time series. For example, the natural logarithm of 480.0 is approximately 6.1738, and the natural logarithm of 510.0 is approximately 6.2344, and so on. After taking the logarithm for all monitoring values, the arithmetic mean and standard deviation are calculated for the 365 logarithmic values. Assuming the arithmetic mean is approximately 6.215 and the standard deviation is approximately 0.095, 6.215 and 0.095 can be regarded as the location parameter and discrete parameter in the logarithmic space. These two parameters have a one-to-one correspondence with the scale and shape in the log-normal distribution model. For the gamma distribution model, the shape parameter and scale parameter can be found through numerical iteration, such that, under this parameter combination, the theoretical mean and variance of the gamma distribution are as close as possible to the sample arithmetic mean of 500.0 and the sample variance of 2534.25 of the geothermal production monitoring time series. For example, we can enumerate the shape parameters in the range of 5.0 to 20.0 with a step size of 0.5, calculate the corresponding candidate values ​​of the scale parameter for each shape parameter, so that the product of the shape parameter and the scale parameter equals 500.0, and then calculate the theoretical variance of the gamma distribution for each set of shape parameters and scale parameters. The set of parameters with the smallest absolute value of the difference between the theoretical variance and the sample variance is regarded as the parameter estimate of the gamma distribution model.

[0025] After obtaining the parameter estimation results for each candidate marginal distribution model, the maximum likelihood objective function value is calculated for each candidate marginal distribution model. Taking the normal distribution model as an example, for each monitoring value in the geothermal production monitoring time series, the probability density value of the monitoring value is calculated using the normal distribution probability density function with a center location parameter of 500.0 and a dispersion parameter of 50.34. Then, the natural logarithm of each of the 365 probability density values ​​is taken, and the 365 natural logarithm values ​​are summed to obtain the log-likelihood value of the normal distribution model. For example, if the calculated log-likelihood value is equal to −1230.5, then this value is taken as the maximum likelihood objective function value of the normal distribution model. For the log-normal distribution model, the probability density value is calculated for each geothermal production monitoring value in the logarithmic space with the arithmetic mean of the logarithmic values ​​(6.215) and the standard deviation (0.095) as parameters, and the log-likelihood value is calculated in the aforementioned manner. Assume that the obtained log-likelihood value is −1185.2. For the gamma distribution model, using the shape and scale parameters obtained from the aforementioned enumeration as input, the gamma distribution probability density value is calculated for each monitoring value, and the natural logarithm is summed. The assumed log-likelihood value is -1178.9. By comparing -1230.5, -1185.2, and -1178.9, it can be seen that the gamma distribution model has the largest log-likelihood value. Therefore, at the level of maximum likelihood objective function value, the gamma distribution model has the best fitting effect on the geothermal production monitoring time series.

[0026] To further evaluate the goodness of fit of the candidate marginal distribution models, this implementation method performs the Kolmogorov–Smirnov test on each candidate marginal distribution model. The specific steps are as follows: First, the 365 monitoring values ​​in the geothermal production monitoring time series are sorted in ascending order, denoted as the 1st to the 365th sorted values. For each sorted value, the empirical cumulative distribution value is calculated by dividing the sorted value's index in the sorting sequence by 365. For example, the empirical cumulative distribution value of the 73rd sorted value is 73 divided by 365, approximately 0.2. Then, based on the candidate marginal distribution models and their parameter estimation results, the theoretical cumulative distribution value of the corresponding sorted value under the theoretical cumulative distribution is calculated. For example, under the gamma distribution model, the theoretical cumulative distribution value corresponding to the 73rd sorted value might be 0.19. For each sorted value, the absolute difference between the empirical cumulative distribution value and the theoretical cumulative distribution value is calculated, and the maximum value among the 365 absolute differences is used as the Kolmogorov–Smirnov test statistic. Assuming the maximum difference is 0.12 under the normal distribution model, 0.09 under the log-normal distribution model, and 0.06 under the gamma distribution model, then the gamma distribution model can be considered to have the best goodness of fit in the Kolmogorov–Smirnov test. By simultaneously comparing the maximum likelihood objective function value and the Kolmogorov–Smirnov test statistic, we can comprehensively measure the model's performance in terms of fitting accuracy and distribution shape consistency, thus avoiding selection bias that may result from relying on only one indicator.

[0027] For groundwater level monitoring time series, vegetation cover monitoring time series, tourism revenue monitoring time series, and carbon emission reduction monitoring time series, this implementation method adopts the same parameter estimation and goodness-of-fit test process as the geothermal production monitoring time series. That is, maximum likelihood estimation is performed on the three candidate marginal distribution models for each monitoring time series to obtain their respective maximum likelihood objective function values. Then, the corresponding goodness-of-fit statistics are calculated through the Kolmogorov-Smirnov test. The candidate marginal distribution model with the largest maximum likelihood objective function value and the smallest Kolmogorov-Smirnov test statistic for the monitoring time series is selected as the marginal distribution model for the monitoring time series. For example, in practical calculations, groundwater level monitoring time series may better conform to a normal distribution model, while vegetation cover monitoring time series may better conform to a beta distribution model. In optional implementations where improved detail is needed at both the low and high value ends, a beta distribution model can be added to the original three candidate marginal distribution models. Tourism revenue monitoring time series and carbon emission reduction monitoring time series typically exhibit right-skewed characteristics. The log-normal distribution model has a larger maximum likelihood objective function value and a smaller Kolmogorov–Smirnov test statistic, making it more likely to be selected as the marginal distribution model for these monitoring time series. After completing the model selection for the five monitoring time series, the marginal distribution models for geothermal production, groundwater level, vegetation cover, tourism revenue, and carbon emission reduction are combined to form a set of marginal distribution models, laying the foundation for the subsequent construction of a standardized sample set and a Copula–Bayesian mixture model.

[0028] In an optional implementation, to take into account the variation characteristics at different time scales, the monitoring period can be extended to 730 days, or the time resolution can be shortened from 1 day to 12 hours. With the shortened time resolution, each record can represent the average geothermal production, average groundwater level, average vegetation cover, tourism revenue, and carbon emission reduction over a continuous 12-hour period. This method can capture intraday fluctuation characteristics, but the construction method of the basic monitoring dataset and monitoring time series still maintains the overall structure of "acquiring five types of indicators at multiple times in the target area and combining the five types of indicators at the same time into a record".

[0029] In one implementation, the basic monitoring dataset and the set of marginal distribution models have been constructed according to the preceding steps. The basic monitoring dataset contains 365 records, each corresponding to one calendar day, including monitoring values ​​for geothermal production, groundwater level, vegetation cover, tourism revenue, and carbon emission reduction. For the five types of monitoring time series, corresponding marginal distribution models have been selected through maximum likelihood estimation and goodness-of-fit tests, and their respective parameter estimates have been obtained, forming the set of marginal distribution models.

[0030] In this implementation, each record in the basic monitoring dataset is first converted into a five-dimensional standardized sample vector containing standardized values ​​of five categories of indicators using the cumulative distribution function of the corresponding marginal distribution model, thus obtaining a standardized sample set. Specifically, taking a single day's record as an example, assume the geothermal production monitoring value is 520.0, the groundwater level monitoring value is 23.4, the vegetation coverage monitoring value is 0.62, the tourism revenue monitoring value is 350,000.0, and the carbon emission reduction monitoring value is 145.0. The aforementioned steps have determined that the marginal distribution model of the geothermal production monitoring time series is a normal distribution model, with an arithmetic mean of 500.0 and a standard deviation of 50.0. For the geothermal production monitoring value of 520.0, first subtract 500.0 from 520.0 to obtain 20.0, then divide 20.0 by 50.0 to obtain 0.4. Using 0.4 as the independent variable of the standard normal distribution, the corresponding cumulative probability found in the cumulative distribution function table of the standard normal distribution is approximately 0.6554, which is the standardized value of the geothermal production. The meaning of this treatment is that, under the marginal distribution model of geothermal production, the probability of randomly observing a geothermal production of no more than 520.0 is about 0.6554. This maps the original value in megawatt-hours to a dimensionless probability position in the interval between 0 and 1, so that the subsequent modeling of the dependent structure is no longer affected by the dimension.

[0031] For a groundwater level monitoring value of 23.4, assuming the marginal distribution model of the groundwater level monitoring time series is also a normal distribution model, with an arithmetic mean of 23.0 and a standard deviation of 0.8, we first subtract 23.0 from 23.4 to get 0.4, then divide 0.4 by 0.8 to get 0.5. Looking up the cumulative probability in the standard normal distribution cumulative distribution function table, the corresponding cumulative probability is approximately 0.6915. This value is used as the standardized groundwater level value. This can be interpreted as 23.4 being slightly higher than the center position in the groundwater level marginal distribution model; therefore, its standardized value is close to 0.7, indicating that the groundwater level on that day was relatively high.

[0032] For a vegetation cover monitoring value of 0.62, assuming the marginal distribution model of the vegetation cover monitoring time series is a beta distribution model in the interval between 0 and 1, and the shape parameters obtained through the aforementioned fitting steps are approximately 3.5 and 2.5. Since the cumulative distribution function of the beta distribution is relatively complex, this embodiment uses the cumulative distribution function calculation routine provided in the numerical calculation software to numerically integrate the independent variable 0.62, obtaining a cumulative probability of approximately 0.68, which is then used as the standardized value of vegetation cover. Because the beta distribution is a type of distribution that flexibly describes the skewed shape in the interval between 0 and 1, through the cumulative distribution function mapping, "vegetation cover of 0.62" can be transformed into "approximately the 68th percentile in the current long-term distribution context," which is beneficial for comparing relative levels between different indicators.

[0033] For the tourism revenue monitoring value of 350,000.0, assuming the marginal distribution model of the tourism revenue monitoring time series is a log-normal distribution model, the previous steps have already taken the natural logarithm of each day's tourism revenue monitoring value, resulting in an arithmetic mean of approximately 12.5 and a standard deviation of approximately 0.6. First, calculate the natural logarithm of 350,000.0, which is approximately 12.765. Then, subtract 12.5 from 12.765 to get 0.265, and divide 0.265 by 0.6 to get approximately 0.4417. Looking up the cumulative probability in the standard normal distribution cumulative distribution function table, the corresponding cumulative probability is approximately 0.6703. This value is used as the standardized value of tourism revenue. The effect of this is that the originally large-scale tourism revenue data is compressed into the 0-1 interval through logarithmic transformation and mapping with the normal cumulative distribution function, smoothing out the influence of maxima and helping to improve the numerical stability of the model.

[0034] For the carbon emission reduction monitoring value of 145.0, assuming the marginal distribution model of the carbon emission reduction monitoring time series is a gamma distribution model, and by back-calculating the sample mean and sample variance, the shape parameter is approximately 98.65 and the scale parameter is approximately 5.07. Since the cumulative distribution function of the gamma distribution involves incomplete integration, this implementation method calls the cumulative distribution function interface of the gamma distribution in the scientific computing library, using 145.0 as the independent variable and 98.65 and 5.07 as parameters, to perform numerical calculations, obtaining a cumulative probability of approximately 0.52. This value is used as the standardized value of the carbon emission reduction. After this processing, different absolute values ​​of carbon emission reduction are mapped to probability position values ​​in the interval between 0 and 1, which is beneficial for subsequent unified processing of the correlation between different indicators.

[0035] Using the above method, the five original monitoring values ​​for that day can be converted into five standardized values ​​between 0 and 1, namely 0.6554, 0.6915, approximately 0.68, 0.6703, and approximately 0.52. These five values ​​are then combined into a five-dimensional standardized sample vector in the order of standardized geothermal production, groundwater level, vegetation cover, tourism revenue, and carbon emission reduction. This process is repeated for each of the 365 records in the basic monitoring dataset, resulting in 365 five-dimensional standardized sample vectors, which are stored chronologically to form a five-dimensional standardized sample set. This method transforms the original multidimensional, dimensional monitoring data into a dimensionless standardized sample set with uniformly distributed margins, thus concentrating all the complexity on dependency modeling. This facilitates the use of Copula structures to specifically characterize the coupling patterns between geothermal production, groundwater level, vegetation cover, tourism revenue, and carbon emission reduction.

[0036] After the standardized sample set is constructed, a target Copula structure model is selected from a pre-set candidate Copula structure set. In this embodiment, the pre-set candidate Copula structure set includes the Clayton Copula structure model, the Gumbel Copula structure model, and the Gaussian Copula structure model. The Clayton Copula structure model emphasizes the lower tail linkage, that is, when the multidimensional standardized samples are simultaneously at small values, it can amplify the correlation between variables, making it suitable for describing the common risk of multiple indicators falling to low levels simultaneously (e.g., reduced geothermal production, declining groundwater levels, and reduced vegetation cover). The Gumbel Copula structure model emphasizes the upper tail linkage, that is, when multiple indicators are simultaneously close to high levels, it can strengthen the probability of joint occurrence, making it suitable for describing situations where "high geothermal production, high tourism revenue, and high carbon emission reduction" occur simultaneously. The Gaussian Copula structure model reflects a structure with overall correlation and relatively mild tail linkage, making it suitable for scenarios with approximately linear correlations. This embodiment selects the most suitable structure from the above-mentioned structures for the characteristics of the target region through data-driven selection.

[0037] The specific selection process is as follows: For each Copula structural model, the range of dependent parameter values ​​is set to a closed interval from 0.1 to 10.0. Within this closed interval, 10 initial parameter values ​​are selected at equal intervals: 0.1, 1.2, 2.3, 3.4, 4.5, 5.6, 6.7, 7.8, 8.9, and 10.0. For a fixed initial parameter value, using this initial parameter value and the standardized sample set, the joint density value of each five-dimensional standardized sample vector is calculated using the corresponding Copula structural model. For numerical stability considerations, this implementation first calculates the natural logarithm of the joint density value of each five-dimensional standardized sample vector, and then sums the natural logarithmic joint density values ​​of all five-dimensional standardized sample vectors to obtain the log-likelihood objective function value corresponding to the initial parameter value. For example, under the Clayton Copula structural model, with an initial dependency parameter value of 2.3, the sum of the logarithmic joint density values ​​of the 365 five-dimensional standardized sample vectors is assumed to be −410.7; under the Gumbel Copula structural model, with an initial dependency parameter value of 3.4, the sum of the logarithmic joint density values ​​is −395.2; and under the Gaussian Copula structural model, with an initial dependency parameter value of 2.3, the sum of the logarithmic joint density values ​​is −402.5. This implementation calculates the corresponding log-likelihood objective function value for each of the 10 initial parameter values ​​for each Copula structural model, and selects the largest value as the initial evaluation value for that Copula structural model. For example, after traversal, the initial evaluation value of the Clayton Copula structural model is -405.3, with an initial dependency parameter value of 1.8; the initial evaluation value of the Gumbel Copula structural model is -390.1, with an initial dependency parameter value of 3.6; and the initial evaluation value of the Gaussian Copula structural model is -398.7, with an initial dependency parameter value of 2.5. Since -390.1 is higher than -405.3 and -398.7, the Gumbel Copula structural model is chosen as the target Copula structural model, and 3.6 is used as the initial value of the dependency parameter. This "coarse search" approach quickly finds a starting point with a good overall fit within a general range of different structures and strengths, avoiding the long convergence time caused by the Markov chain Monte Carlo method starting with completely unreasonable parameters.

[0038] After determining the target Copula structure model, a prior distribution needs to be constructed for the dependency parameters. In this embodiment, the length of the dependency parameter's value interval is 10.0 minus 0.1, which equals 9.9. Ten percent of 9.9 is set as the standard deviation of the prior distribution, resulting in 0.99. The initial value of the aforementioned dependency parameter, 3.6, is used as the center value of the prior distribution. A unimodal symmetric prior distribution is constructed based on this center value and this standard deviation. In actual implementation, grid points can be divided within the interval from 0.1 to 10.0 with a fixed step size. The distance between each grid point and the center value is substituted into the commonly used normal shape calculation process to obtain the relative weight of each grid point. Then, all weights are normalized so that their sum over the interval from 0.1 to 10.0 is 1, thus obtaining the discretized dependency parameter prior distribution. For values ​​less than 0.1 or greater than 10.0, this embodiment does not generate them in the grid, thus naturally keeping the dependency parameter values ​​within a predetermined closed interval. By constructing a prior distribution centered at 3.6, we can assign a preference to "parameters near the current optimum are more likely" in subsequent Bayesian updates, while still retaining the ability to explore other parameter values ​​far from 3.6.

[0039] After determining the prior distribution of the dependency parameters, this implementation uses the Markov chain Monte Carlo method based on a standardized sample set to perform Bayesian updates on the dependency parameters, obtaining a posterior sample set and constructing a Copula-Bayes mixture model. The Markov chain Monte Carlo method employs a random walk-based accept-reject sampling strategy. First, the total number of iterations in the Markov chain is set to 10,000 steps, with the first 2,000 steps defined as the burning phase, used to gradually bring the chain closer to the high posterior density region from its initial position. In the first iteration, the current value of the dependency parameters is set to 3.6, which is the initial value of the dependency parameters obtained based on the log-likelihood objective function mentioned above. Starting from the second iteration, in each iteration, a random number generation algorithm is first called to generate a perturbation value that conforms to a normal distribution with a mean of 0.0 and a standard deviation of 0.1. For example, in the second iteration, the generated perturbation value is 0.08, which is added to the current value of the dependency parameters of 3.6 from the previous step to obtain a candidate value of 3.68 for the dependency parameters. If the candidate value of the dependency parameter is less than 0.1, it is directly replaced with 0.1; if the candidate value of the dependency parameter is greater than 10.0, it is directly replaced with 10.0. In this embodiment, 3.68 is between 0.1 and 10.0, so no amplitude limiting is required.

[0040] After obtaining a candidate dependency parameter value of 3.68, the joint density value of each five-dimensional standardized sample vector is calculated using this candidate value and the standardized sample set through the objective GumbelCopula structure model. The natural logarithm of each joint density value is then taken, and the 365 natural logarithmic joint density values ​​are summed to obtain the log-likelihood sum corresponding to the dependency parameter candidate value. This log-likelihood sum is assumed to be −388.4. Subsequently, based on the previously constructed dependency parameter prior distribution, the prior probability density value corresponding to the grid point close to 3.68 is found in the grid. Assuming the natural logarithm of this probability density value is −1.2, −388.4 is added to −1.2 to obtain −389.6, which is regarded as the log-posterior objective function value corresponding to the dependency parameter candidate value. Similarly, for the current value of the dependency parameter 3.6 in the previous iteration, the above calculation process is repeated to obtain the sum of the log-likelihoods corresponding to the current value of the dependency parameter as -390.1, and the natural logarithm of the prior probability density value as -1.1. Adding these two together yields -391.2, which is used as the log-posterior objective function value corresponding to the current value of the dependency parameter. Then, the difference between the two is calculated, i.e., subtracting -391.2 from -389.6 to obtain 1.6. Inputting 1.6 into the exponential function yields approximately 4.95. Since 4.95 is greater than 1, the acceptance probability is set to 1. Next, a random number generation algorithm is called to generate a uniformly distributed random number between 0 and 1, for example, 0.37. Since 0.37 is less than 1, in this implementation, the candidate value of the dependency parameter 3.68 is accepted as the new current value of the dependency parameter in the second iteration.

[0041] In step 3 and subsequent iterations, the process of "adding perturbations to obtain candidate values ​​for dependency parameters, calculating the sum of log-likelihoods, combining the prior distribution to obtain the log-posterior objective function value, calculating the difference and obtaining the acceptance probability through an exponential function, generating uniformly distributed random numbers and deciding whether to accept or reject" is repeated. When the acceptance probability is less than 1, if the generated uniformly distributed random number is greater than the acceptance probability, the current value of the dependency parameter from the previous iteration is retained unchanged in the current step. This "rejection" operation allows the Markov chain to slowly wander near the high-density region of the posterior distribution, which is beneficial for detailed sampling in local regions. Through multiple acceptances and rejections, the value of the dependency parameter will dynamically fluctuate around the high-density region of the posterior distribution within the range supported by the prior distribution. After completing 10,000 iterations, the first 2,000 steps are discarded as the current values ​​of the dependency parameters corresponding to the burning stage, and the current values ​​of the dependency parameters from each of the remaining 8,000 steps are collected to form a posterior sample set of dependency parameters. For example, the values ​​in the posterior sample set of dependent parameters may be concentrated between 2.8 and 4.5, with an arithmetic mean of approximately 3.7 and a sample standard deviation of approximately 0.4, indicating that after considering the standardized sample set and prior information, the upper tail coupling strength among the five indices in the target region tends to be at a moderately strong level.

[0042] In this implementation, the target GumbelCopula structural model and the posterior sample set of dependency parameters are combined to construct a Copula-Bayesian mixture model. Specifically, when generating standardized joint samples, a random number generation algorithm is first used to randomly select an integer between 1 and 8000. This integer is used as an index to retrieve the corresponding dependency parameter value from the posterior sample set of dependency parameters. Then, this dependency parameter value is substituted into the target GumbelCopula structural model, and a set of five-dimensional standardized joint samples is generated using an existing multidimensional Copula random number generation algorithm. In this way, a large number of slightly different GumbelCopula structural models are superimposed with the posterior sample set of dependency parameters as weights to form a joint distribution "hybrid" under parameter uncertainty. Compared with a fixed Copula model that only uses a single dependency parameter point estimate, the Copula-Bayesian mixture model can more fully reflect the dependency structure uncertainty inherent in the data. When there are nonlinear and multi-tailed characteristics between geothermal production, groundwater level, vegetation cover, tourism revenue, and carbon emission reduction, it can provide a more robust characterization of the joint distribution.

[0043] In an optional implementation, the pre-defined candidate Copula structure set can include a Student distribution Copula structure model in addition to the Clayton Copula, Gumbel Copula, and Gaussian Copula structures to better describe tail coupling when extreme covariances are particularly significant. The dependency parameters can be extended to a parameter vector containing two or more components, controlling the overall correlation strength and tail correlation strength respectively. In the Markov chain Monte Carlo method, the standard deviation of the perturbation values ​​can be dynamically adjusted in the first few hundred iterations. For example, the standard deviation can be appropriately increased when there are too many consecutive acceptances and appropriately decreased when there are too many consecutive rejections, to maintain the acceptance rate between 0.2 and 0.5, thereby improving sampling efficiency.

[0044] In this embodiment, the basic monitoring dataset, the marginal distribution model set, and the Copula-Bayesian mixture model have been constructed according to the aforementioned implementation method. The goal is to use Monte Carlo simulation to generate five-dimensional uniform random numbers, input them into the Copula-Bayesian mixture model to obtain standardized joint samples, then recover the five types of indicator samples through the inverse function of the marginal distribution model set, calculate the comprehensive value sample according to the preset benefit-cost rule, and finally calculate the coupled value expectation and risk index based on the comprehensive value sample set to obtain the geothermal ecosystem service coupled value assessment result.

[0045] In this implementation, the number of Monte Carlo simulation rounds is first determined. To strike a balance between computational load and statistical accuracy, the number of simulation rounds is set to 10,000. The reason for choosing 10,000 rounds is that the larger the number of rounds, the closer the comprehensive value sample set is to the true theoretical distribution, and the more stable the coupled value expectation and risk indicators are, but the computation time also increases accordingly; 10,000 rounds can usually keep the estimation error within an acceptable range for engineering applications, and at the same time, the calculation can be completed on a conventional server within a reasonable time.

[0046] In each round of simulation, five independent uniformly distributed random numbers with values ​​between 0 and 1 are generated consecutively using a random number generation algorithm. For example, in the first round of simulation, the five generated independent uniformly distributed random numbers are 0.37, 0.52, 0.80, 0.62, and 0.41. These five values ​​are independent of each other, representing unbiased random points within a unit interval. While directly using these five values ​​can yield five types of index samples through the inverse function of the marginal distribution model set, it cannot reflect the actual linkage between geothermal production, groundwater level, vegetation cover, tourism revenue, and carbon emission reduction. Therefore, these five independent uniformly distributed random numbers need to be input into the Copula-Bayes mixture model, allowing the existing dependency structure to transform them, resulting in standardized joint samples that include dependencies.

[0047] In the first round of simulation, a dependency parameter sample value is randomly selected from the posterior sample set of dependency parameters. For example, the sample with index 2500 is randomly selected, corresponding to a dependency parameter value of 3.7. This dependency parameter value, along with five independent uniform random numbers, is input into the target Copula structure model. A multidimensional transformation is performed in the Copula-Bayes mixture model to obtain a set of five-dimensional standardized joint samples. For example, the transformed five-dimensional standardized joint samples are 0.58, 0.65, 0.73, 0.81, and 0.49. Compared with the original 0.37, 0.52, 0.80, 0.62, and 0.41, these five values ​​still fall between 0 and 1, but they are no longer independent of each other. Instead, they reflect the characteristics learned through the Copula-Bayes mixture model, such as "simultaneous high upper tail" and "interconnected at a medium level." This approach allows all dependency structures to be processed in the Copula layer without changing the marginal distribution shape of each indicator, making the construction of the joint distribution highly flexible.

[0048] After obtaining the five-dimensional standardized joint sample, the five index samples are recovered using the inverse function of the marginal distribution model set. Specifically, the first component of the five-dimensional standardized joint sample, 0.58, is input into the inverse function of the geothermal production marginal distribution model. Assuming the marginal distribution model selected in the aforementioned steps for the geothermal production monitoring time series is a gamma distribution model with a shape parameter of approximately 98.65 and a scale parameter of approximately 5.07, the gamma distribution quantile calculation function from the scientific computing library is called, using 0.58 as the cumulative probability and 98.65 and 5.07 as parameters, to obtain the corresponding geothermal production sample. For example, a calculation result of 504.0 indicates that, under the long-term distribution background, the probability of geothermal production being less than or equal to 504.0 is approximately 58%. Therefore, 504.0 can be considered a simulation result "slightly higher than the long-term median level".

[0049] The second component of the five-dimensional standardized joint sample, 0.65, is input into the inverse function of the groundwater level marginal distribution model. Assuming the marginal distribution model of the groundwater level monitoring time series is a normal distribution model with an arithmetic mean of 23.0 and a standard deviation of 0.8, we first look up the standardized value corresponding to the cumulative probability of 0.65 in the standard normal distribution quantile table. This standardized value is approximately 0.385. Then, we multiply 0.385 by 0.8 to get 0.308, and add 0.308 to 23.0 to get 23.308. For practical applications, to retain three decimal places, this value can be rounded to 23.308 and used as the groundwater level sample. In this way, the groundwater level sample not only reflects the randomness of the current simulation cycle but also maintains consistency with the long-term monitoring distribution, making the subsequent risk characterization more closely resemble actual hydrological conditions.

[0050] The third component, 0.73, of the five-dimensional standardized joint sample is input into the inverse function of the vegetation cover marginal distribution model. Assuming the marginal distribution model of the vegetation cover monitoring time series is a beta distribution model with shape parameters of approximately 3.5 and 2.5, respectively, using the beta distribution quantile calculation function and taking 0.73 as the cumulative probability, the vegetation cover sample is approximately 0.68. This indicates that in the long-term vegetation cover distribution, 0.68 is approximately at the 73rd percentile, representing a relatively good vegetation condition.

[0051] The fourth component of the five-dimensional standardized joint sample, 0.81, is input into the inverse function of the tourism revenue marginal distribution model. Assuming the marginal distribution model of the tourism revenue monitoring time series is a log-normal distribution model with an arithmetic mean of 12.5 and a standard deviation of 0.6 in the log space, we first look up the standardized value corresponding to the cumulative probability of 0.81 in the standard normal distribution quantile table. This standardized value is approximately 0.877. Multiplying 0.877 by 0.6 gives 0.526. Adding 0.526 to 12.5 gives 13.026. Then, we perform an exponential operation on this logarithmic value to obtain the tourism revenue sample. The calculated tourism revenue sample is approximately 451,000.0 monetary units, indicating that tourism revenue is at a long-term high level in this simulation.

[0052] The fifth component, 0.49, of the five-dimensional standardized joint sample is input into the inverse function of the marginal distribution model of carbon emission reductions. Assuming the marginal distribution model of the carbon emission reduction monitoring time series is also a gamma distribution model with a shape parameter of approximately 100.0 and a scale parameter of approximately 1.6, the gamma distribution quantile calculation function is used, and 0.49 is used as the cumulative probability to calculate the carbon emission reduction sample, which is approximately 158.0 tons. This indicates that under the long-term carbon emission reduction distribution, the simulation result is close to the median level.

[0053] Through the above process, a set of restored samples of the five categories of indicators can be obtained in the first round of simulation: geothermal production sample 504.0, groundwater level sample 23.308, vegetation coverage sample 0.68, tourism revenue sample 451000.0, and carbon emission reduction sample 158.0. In subsequent simulation rounds, the above steps are repeated. For each simulation round, five new independent uniformly distributed random numbers are generated, transformed using a Copula-Bayesian mixture model to obtain new five-dimensional standardized joint samples, and then the five categories of indicator samples are restored using the inverse function of the marginal distribution model set. After completing 10,000 simulation rounds, 10,000 sets of geothermal production, groundwater level, vegetation coverage, tourism revenue, and carbon emission reduction samples are obtained, providing input for the calculation of the comprehensive value samples.

[0054] refer to Figure 2This graph visually reflects the dependency patterns between standardized variables after marginal distribution transformation, particularly highlighting the nonlinear coupling between standardized geothermal output and standardized tourism revenue. In the graph, the horizontal axis represents standardized geothermal output, and the vertical axis represents standardized tourism revenue, both strictly confined to a closed interval between zero and one. This means that regardless of whether the original data is measured in megawatt-hours or monetary units, it is transformed into its cumulative probability position within its respective marginal distribution under this coordinate system, thus eliminating the interference of dimensional differences on correlation assessment. The numerous black dots scattered throughout the graph represent simulated sample points generated after updating parameters using the Markov chain Monte Carlo method. Observing the distribution characteristics of the scattered points reveals that these sample points do not uniformly fill the entire square area, but rather exhibit a clear asymmetric clustering pattern. In the lower left corner of the graph, the sample points are relatively sparse and dispersed, indicating that when both geothermal output and tourism revenue are at low levels, the correlation between the two is relatively weak; that is, low geothermal output does not necessarily correspond to extremely low tourism revenue, and the two exhibit a certain degree of independence in the low-value region. However, moving the viewpoint towards the upper right corner of the graph reveals a clear convergence of sample points towards the diagonal region, with a significant increase in density. This high-density, peak-like distribution in the upper right corner is a typical characteristic of the Gumbel-Copula structural model, known as upper-tail correlation. Physically, it reveals a key ecosystem service coupling law: when local heat production remains at extremely high levels, it is often accompanied by extremely high levels of tourism revenue; the probability of both occurring simultaneously at extreme values ​​is far higher than the product of their individual probabilities. In addition to the scatter plot, the graph also includes several closed contour lines. These contour lines represent the level set of standardized joint probability densities, with values ​​gradually increasing from the outside in. The shape of the contour lines further confirms the upper-tail correlation characteristic, exhibiting a sharp, outward-protruding shape in the upper right corner and a smoother, more rounded shape in the lower left corner. This contour line distribution clearly reveals the range of variable combinations at different joint probability levels. For example, the region within the innermost contour lines represents the most frequently occurring state combinations in this coupled system. For the evaluator, Figure 2 This study not only validated that the selected Copula model effectively captures the asymmetric dependencies between variables, but more importantly, it provides an intuitive basis for subsequent risk analysis. If a conventional normal distribution or Gaussian Copula were used, the graph would present a centrally symmetric ellipse, failing to describe the synergistic enhancement effect in this high-value range, leading to either an overestimation or underestimation of the synergistic benefits of geothermal energy and tourism. Therefore, Figure 2 This invention establishes the technical advantages of handling nonlinear and asymmetric ecosystem service value coupling assessment and proves the necessity of introducing the Copula Bayesian mixture model.

[0055] In each simulation round, based on the above five types of indicator samples, a comprehensive value sample is calculated according to a preset revenue-cost rule. In this embodiment, the unit price of geothermal energy revenue is 30.0 monetary units per megawatt-hour, the tourism industry surcharge coefficient is 1.2, the unit price of carbon emission reduction revenue is 60.0 monetary units per ton, the groundwater level safety threshold is the arithmetic mean of groundwater level monitoring values ​​during the reference period minus 0.5 meters, and the ecosystem service gain coefficient is the ecosystem service value corresponding to 10.0 monetary units per unit of vegetation cover. Assuming the arithmetic mean of groundwater level monitoring values ​​during the reference period is 23.0 meters, then the groundwater level safety threshold is 23.0 minus 0.5, i.e., 22.5 meters. Assuming the unit cost of ground subsidence remediation is 100,000.0 monetary units per meter, this means that for every 1.0 meter drop in groundwater level, 100,000.0 monetary units need to be invested in ground subsidence remediation.

[0056] Taking the first round of simulation as an example, the geothermal production sample is 504.0, so the geothermal energy revenue equals 504.0 multiplied by 30.0, resulting in 15120.0 monetary units. The tourism revenue sample is 451000.0, so the comprehensive tourism revenue equals 451000.0 multiplied by 1.2, resulting in 541200.0 monetary units. The carbon emission reduction sample is 158.0, so the emission reduction revenue equals 158.0 multiplied by 60.0, resulting in 9480.0 monetary units. The vegetation cover sample is 0.68, so the ecosystem service gain equals 0.68 multiplied by 10.0, resulting in 6.8 monetary units. The groundwater level sample is 23.308, which is higher than the groundwater level safety threshold of 22.5, therefore the land subsidence control cost in this round of simulation is 0.

[0057] Summarizing the aforementioned benefits and costs, the comprehensive value sample equals geothermal energy benefits plus tourism benefits plus emission reduction benefits plus ecosystem service gains, and then subtracts the cost of land subsidence remediation. Substituting the values, the comprehensive value sample equals 15120.0 plus 541200.0 plus 9480.0 plus 6.8, minus 0, resulting in 15120.0 plus 541200.0 equaling 556320.0. Adding 9480.0 gives 565800.0, and adding 6.8 gives 565806.8 monetary units. This value is the comprehensive value sample corresponding to the first round of simulation. A larger value indicates that in this round of simulation, the comprehensive benefits from geothermal energy supply, tourism development, carbon emission reduction, and ecosystem services are higher, while the remediation costs from negative impacts such as land subsidence are lower.

[0058] In the second and subsequent simulations, the same multiplication and addition operations and threshold comparisons are performed on the corresponding five categories of indicator samples to generate comprehensive value samples round by round. After completing 10,000 simulations, a sample set containing 10,000 comprehensive value samples is obtained. This sample set reflects various scenarios in which the coupled comprehensive value of geothermal ecosystem services may emerge in the future under the current geothermal development strategies and ecological environment conditions. Statistical analysis of this sample set can yield the expected coupled value and risk indicators.

[0059] Specifically, the sample set containing 10,000 composite value samples is first sorted in ascending order. After sorting, each composite value sample is assigned a position number, ranging from position number 1 to position number 10,000. The composite value sample value at position number 500 after sorting is used as the first value-at-risk (VAT) indicator. For example, if the composite value sample value at position number 500 after sorting is 360,000.0 currency units, then the first VAT indicator is 360,000.0 currency units. This value indicates that, under the current strategy, the composite value is lower than 360,000.0 currency units in approximately 5% of scenarios, and higher than 360,000.0 currency units in the remaining 95% of scenarios. Therefore, this indicator can be used to characterize the lower limit of the composite value under more unfavorable scenarios.

[0060] Next, the arithmetic mean of the values ​​of the 500 composite value samples numbered from 1 to 500 after sorting is calculated to obtain the second value-at-risk (VAT) index. For example, if the arithmetic sum of these 500 composite value samples is 170,000,000.0 units, then dividing 170,000,000.0 by 500 yields 340,000.0 units, which is used as the second VAT index. Unlike the first VAT index, which only focuses on a single point at the "5th percentile," the second VAT index is an average representation of the composite value level of the 500 worst-case scenarios. It reflects the typical level of composite value under persistently adverse conditions and is therefore more suitable for describing the comprehensiveness of extreme risks.

[0061] Finally, the arithmetic mean of the values ​​from 10,000 comprehensive value samples is calculated to obtain the expected comprehensive value index. For example, if the arithmetic sum of the 10,000 comprehensive value samples is 520,000,000.0 monetary units, then dividing 520,000,000.0 by 10,000 yields 52,000.0 monetary units, which is used as the expected comprehensive value index. The expected comprehensive value index corresponds to the average level of comprehensive value under all possible scenarios and is mainly used to measure the overall profitability of geothermal ecosystem service coupling value in long-term operation.

[0062] In this implementation, the first risk value index, the second risk value index, and the comprehensive value expectation index are collectively used as the output of the geothermal ecosystem service coupled value assessment. Decision-makers can use the comprehensive value expectation index to determine the average long-term return of a geothermal development and ecological protection strategy, the first risk value index to determine the minimum return of the strategy under worst-case scenarios, and the second risk value index to determine the average loss under several unfavorable scenarios, thereby weighing multiple alternative strategies.

[0063] In an optional implementation, the number of simulation rounds can be adjusted from 10,000 to 20,000 to further improve the estimation accuracy of the coupled value expectation and risk indicators; the location number of the first risk value indicator can be adjusted from 500 to 100, corresponding to the worst 1% of the comprehensive value samples after ranking, for a more conservative assessment of extreme adverse scenarios; the calculation of the second risk value indicator can be done by averaging the comprehensive value samples from location number 1 to location number 1000 to balance the impact of individual extreme samples. When these optional implementations are adopted, the overall process of generating five-dimensional uniform random numbers using Monte Carlo simulation, obtaining standardized joint samples by inputting them into a Copula-Bayes mixture model, recovering the five types of indicator samples through the inverse function of the marginal distribution model set, and calculating the comprehensive value samples according to the preset benefit-cost rule remains unchanged. Only the sample size and risk quantile settings are adjusted, thereby allowing the same geothermal ecosystem service coupled value assessment framework to be used flexibly under different risk preferences and accuracy requirements.

[0064] refer to Figure 3The graph contains three dimensions: the horizontal axis represents geothermal output in megawatt-hours; the vertical axis represents tourism revenue in monetary units; and the vertical axis represents the joint probability density value, used to measure the likelihood of a specific combination of geothermal output and tourism revenue. The entire surface presents a mountain-like structure rising from the bottom, with the total volume below the surface mathematically strictly equal to one, encompassing all possible scenarios. A closer look at the surface's shape reveals that it is not a standard cone or dome, but rather a ridge extending along a specific direction. The highest point of the surface, the apex of the mountain, corresponds to the most likely numerical combination of geothermal output and tourism revenue, typically representing the system's rated operating state or long-term average level. Extending outwards from the peak, the surface height gradually decreases, but the rate of decrease varies significantly in different directions. Along the diagonal direction from low geothermal output and low revenue to high geothermal output and high revenue, the slope of the surface is relatively gentle, forming a distinct ridge. The ridge's orientation profoundly reveals the positive correlation between two variables: as geothermal production increases, the high-value area of ​​the joint probability density shifts towards the direction of increased tourism revenue. This implies that in actual observations, high-yield geothermal extraction activities often occur simultaneously with high tourism economic returns, exhibiting a macroscopic trend of synergistic growth. Conversely, perpendicular to the ridge, the surface height rapidly drops to near zero, indicating an extremely low probability of divergent scenarios such as "extremely high geothermal production paired with extremely low tourism revenue" or "extremely low geothermal production paired with extremely high tourism revenue." Furthermore, the surface's edge contour and cross-sectional shape are actually determined by their respective edge distribution models. Integrating along the horizontal axis yields a projection curve that conforms to the characteristics of a gamma distribution, exhibiting a right-skewed long-tailed shape; integrating along the vertical axis yields a log-normal distribution. Figure 3 The curved surface does not exhibit perfect symmetry; its tail extends further in the upper right corner, which is consistent with... Figure 2 Corresponding to the upper-tail correlation structure in the text, in physical space, this manifests as a more significant fluctuation range in the system under extremely favorable conditions (such as the peak heating season coinciding with the peak tourist season) compared to under extremely unfavorable conditions. By constructing and displaying this joint probability density surface, this embodiment allows decision-makers to intuitively see the system's "comfort zone" (near the peak) and potential "risk zones" or "high-yield zones" (edge ​​areas). This three-dimensional visualization transforms abstract statistical coupling parameters into concrete probabilistic topographic maps, enabling resource allocation and capacity planning for geothermal ecosystems to be based on a precise quantification of joint uncertainties, avoiding the potential for biased decisions caused by relying on a single indicator.

[0065] The present invention has been described in detail above. Specific examples have been used to illustrate the principles and implementation methods of the invention. The descriptions of the embodiments above are merely for the purpose of helping to understand the method and core ideas of the present invention. It should be noted that those skilled in the art can make various improvements and modifications to the present invention without departing from its principles, and these improvements and modifications also fall within the protection scope of the claims of the present invention.

Claims

1. A method for assessing the coupled value of geothermal ecosystem services, characterized in that, Includes the following steps: Five indicators—geothermal output, groundwater level, vegetation cover, tourism revenue, and carbon emission reduction—were acquired at multiple times in the target area. These five indicators at the same time were combined into a basic monitoring dataset. Based on this dataset, five monitoring time series were constructed, and parameter estimation and goodness-of-fit tests were performed on each time series to determine their marginal distribution models, forming a set of marginal distribution models. Each record in the basic monitoring dataset was transformed into a five-dimensional sample containing standardized values ​​of the five indicators using the cumulative distribution function of the corresponding marginal distribution model, resulting in a standardized sample set. This set was then used within a pre-defined Copula structure. A target Copula structural model is selected from the selection set, and its dependency parameters are updated Bayesianly using the Markov chain Monte Carlo method based on the standardized sample set to obtain the posterior sample set and construct a Copula-Bayes mixture model. Five-dimensional uniform random numbers are generated using Monte Carlo simulation and input into the Copula-Bayes mixture model to obtain standardized joint samples. The five types of index samples are recovered through the inverse function of the marginal distribution model set. The comprehensive value sample is calculated according to the preset benefit-cost rule. Based on the comprehensive value sample set, the expected value of coupling and risk indicators are calculated to obtain the evaluation results of the coupling value of geothermal ecosystem services.

2. The method according to claim 1, characterized in that, The acquisition of the basic monitoring dataset includes: setting up at least 5 geothermal production monitoring points, at least 5 groundwater level monitoring points, and at least 5 vegetation coverage monitoring sample areas within the target area; setting up at least 1 tourism revenue metering point at the main tourist entrance; and setting up at least 1 carbon emission reduction calculation point at the entrance and exit of the geothermal energy supply facility. Geothermal production monitoring values, groundwater level monitoring values, vegetation coverage monitoring values, tourism revenue monitoring values, and carbon emission reduction monitoring values ​​are continuously collected for at least 365 days at a time interval of 1 day. The above 5 monitoring values ​​collected on the same natural day are combined into 1 record to form a basic monitoring dataset containing at least 365 records.

3. The method according to claim 1, characterized in that, In constructing the set of marginal distribution models for geothermal ecosystem services, candidate marginal distribution models include normal distribution models, log-normal distribution models, and gamma distribution models. For each monitoring time series, maximum likelihood estimation is performed on the normal distribution model, log-normal distribution model, and gamma distribution model respectively to obtain the corresponding maximum likelihood objective function value. The Kolmogorov-Smirnov test is used to obtain the goodness-of-fit statistic. Among the normal distribution model, log-normal distribution model, and gamma distribution model, the model with the largest maximum likelihood objective function value and the smallest Kolmogorov-Smirnov test statistic is selected as the marginal distribution model for that monitoring time series.

4. The method according to claim 1, characterized in that, In the process of constructing the standardized sample set, the geothermal production monitoring value of each record in the geothermal ecosystem service coupling basic monitoring dataset is input into the cumulative distribution function of the geothermal production marginal distribution model to obtain the standardized geothermal production value. The groundwater level monitoring value is input into the cumulative distribution function of the groundwater level marginal distribution model to obtain the standardized groundwater level value. The vegetation cover monitoring value is input into the cumulative distribution function of the vegetation cover marginal distribution model to obtain the standardized vegetation cover value. The tourism revenue monitoring value is input into the cumulative distribution function of the tourism revenue marginal distribution model to obtain the standardized tourism revenue value. The carbon emission reduction monitoring value is input into the cumulative distribution function of the carbon emission reduction marginal distribution model to obtain the standardized carbon emission reduction value. The standardized sample vectors are constructed in the order of geothermal production, groundwater level, vegetation cover, tourism revenue, and carbon emission reduction. All five-dimensional standardized sample vectors are stored in the order of records to form a five-dimensional standardized sample set.

5. The method according to claim 1, characterized in that, In selecting the target Copula structure model, the pre-set candidate Copula structure set includes at least the Clayton Copula, Gumbel Copula, and Gaussian Copula structures. For each Copula structure model, the range of dependent parameter values ​​is set to a closed interval from 0.1 to 10.

0. Within this closed interval, 10 equally spaced initial parameter values ​​are selected. For each initial parameter value, the log-likelihood value of all samples is calculated using the corresponding Copula structure model and the five-dimensional standardized sample set, and the sum is used to obtain the log-likelihood objective function value. The maximum value of the log-likelihood objective function value of all initial parameter values ​​corresponding to the same Copula structure model is taken as the initial evaluation value of the Copula structure model. Among the Clayton Copula, Gumbel Copula, and Gaussian Copula structures, the Copula structure model with the largest initial evaluation value is selected as the target Copula structure model.

6. The method according to claim 1, characterized in that, The method for constructing the prior distribution of dependent parameters of the target Copula structural model includes: setting the initial value of the parameter that maximizes the initial evaluation value of the target Copula structural model as the center value of the prior distribution of dependent parameters; setting 10% of the length of the range of dependent parameter values ​​as the standard deviation of the prior distribution of dependent parameters; constructing a normal prior distribution with the center value and the standard deviation as parameters; and maintaining the dependent parameter values ​​within a closed interval of 0.1 to 10.0 by means of amplitude limiting within the interval where the dependent parameter values ​​are less than 0.1 or greater than 10.

0.

7. The method according to claim 1, characterized in that, The Markov chain Monte Carlo method performs a multi-step Bayesian update of the dependency parameters, which includes: setting the total number of iterations of the Markov chain to 10,000 steps, defining the first 2,000 steps as the burning stage, and using the initial parameter values ​​as the current values ​​of the dependency parameters in the first step. Starting from the second step, in each iteration, a normally distributed random number with a mean of 0.0 and a standard deviation of 0.1 is used as a perturbation value and added to the current value of the dependency parameter in the previous step to obtain a candidate value of the dependency parameter. Candidate values ​​of dependency parameters less than 0.1 are set to 0.1, and candidate values ​​of dependency parameters greater than 10.0 are set to 10.

0. The candidate values ​​of the dependency parameters and the five-dimensional standardized sample set are then used to... The target Copula structural model calculates and sums the log-likelihood values ​​of all samples. This sum of log-likelihoods is then added to the log density value of the candidate dependency parameter under the prior distribution of the dependency parameter to obtain the log-posterior objective function value corresponding to the candidate dependency parameter. The log-posterior objective function value obtained in the same way from the previous step of the current dependency parameter is used for comparison. The acceptance probability is calculated based on the difference between the two values. Uniformly distributed random numbers with values ​​between 0 and 1 are generated. The current dependency parameter value of the current iteration step is determined based on the acceptance probability and the uniformly distributed random numbers. After 10,000 iterations, the current dependency parameter values ​​except for the burning stage are used to form the dependency parameter posterior sample set.

8. The method according to claim 1, characterized in that, The process of assessing the coupled value of geothermal ecosystem services also includes: collecting monitoring values ​​of geothermal output, groundwater level, vegetation coverage, tourism revenue, and carbon emission reduction in the second and subsequent batches under the same conditions as the monitoring point layout and time intervals; converting each batch of monitoring values ​​into a new five-dimensional standardized sample set; using the arithmetic mean of the posterior sample set of dependent parameters as the new central value of the prior distribution of dependent parameters; using the sample standard deviation of the posterior sample set of dependent parameters as the new standard deviation of the prior distribution of dependent parameters; constructing a new prior distribution of dependent parameters; and performing Bayesian updates of dependent parameters on the new five-dimensional standardized sample set using the Markov chain Monte Carlo method to obtain the updated posterior sample set of dependent parameters. This updated posterior sample set of dependent parameters is then combined with the target Copula structural model to form a target Copula-Bayes hybrid model with dynamic dependent parameters.

9. The method according to claim 1, characterized in that, In the Monte Carlo simulation, the number of simulation rounds was set to 10,000. In each round, a random number generation algorithm was used to generate five independent uniformly distributed random numbers with values ​​between 0 and 1. These five independent uniformly distributed random numbers were input into the Copula-Bayes mixture model. Under the joint distribution constraints of the target Copula structure model, a set of five-dimensional standardized joint samples was obtained. The five components of the five-dimensional standardized joint samples were input into the inverse functions of the geothermal production edge distribution model, the groundwater level edge distribution model, the vegetation cover edge distribution model, the tourism revenue edge distribution model, and the carbon emission reduction edge distribution model, respectively, to obtain geothermal production samples, groundwater level samples, vegetation cover samples, tourism revenue samples, and carbon emission reduction samples.

10. The method according to claim 1, characterized in that, The calculation method for the comprehensive value sample of geothermal ecosystem services includes: multiplying the geothermal production sample by the unit price of geothermal energy revenue to obtain geothermal energy revenue; multiplying the tourism revenue sample by the tourism industry surcharge to obtain comprehensive tourism revenue; multiplying the carbon emission reduction sample by the unit price of carbon emission reduction revenue to obtain emission reduction revenue; comparing the groundwater level sample with the groundwater level safety threshold; when the groundwater level sample is lower than the groundwater level safety threshold, the ground subsidence control cost is calculated by multiplying the difference between the groundwater level safety threshold and the groundwater level sample by the unit cost of ground subsidence control; when the groundwater level sample is higher than or equal to the groundwater level safety threshold, the ground subsidence control cost is set to 0; multiplying the vegetation cover sample by the ecosystem service gain coefficient to obtain ecosystem service gain; and adding the geothermal energy revenue, comprehensive tourism revenue, emission reduction revenue, and ecosystem service gain, and then subtracting the ground subsidence control cost to obtain the comprehensive value sample. The unit price of geothermal energy revenue is 30.0 monetary units per megawatt-hour, the tourism industry surcharge is 1.2, the unit price of carbon emission reduction revenue is 60.0 monetary units per ton, the groundwater level safety threshold is the arithmetic mean of groundwater level monitoring values ​​during the reference period minus 0.5 meters, and the ecosystem service gain coefficient is the ecosystem service value corresponding to 10.0 monetary units per unit of vegetation cover. In a sample set containing 10,000 comprehensive value samples, the comprehensive value samples are sorted in ascending order. The sample value with position number 500 after sorting is used as the first risk value indicator, the arithmetic mean of comprehensive value samples with position numbers 1 to 500 after sorting is used as the second risk value indicator, and the arithmetic mean of the sample set is used as the comprehensive value expectation indicator. The first risk value indicator, the second risk value indicator, and the comprehensive value expectation indicator are used together as the result of the geothermal ecosystem service coupling value assessment.

Images