Method for hydrological elastic attribution of a watershed based on chain mechanism decomposition

Abstract

The application discloses a kind of based on chain mechanism decomposition's watershed hydrology elastic attribution method, comprising: obtaining the meteorological hydrology sequence of target watershed, underlying surface snow sequence and climate state variable;Based on meteorological hydrology sequence and underlying surface snow sequence, extract the intermediate mechanism index of characterizing physical transmission path;Determine the empirical sensitivity of each intermediate mechanism index to snow change;Intermediate mechanism index and climate state variable are input into the watershed parameter model constructed in advance, and obtain partial sensitivity;Two kinds of sensitivity and basic environmental scaling factor are combined to calculate explicit mechanism elastic component;Based on explicit mechanism elastic component, calculate residual elastic component;Aggregation explicit mechanism elastic component and residual elastic component, output hydrology elastic attribution result.The application realizes the quantitative decoupling and complete attribution of multiple physical conduction channels, improves the objectivity and precision of complex hydrology response analysis.

Description

Technical Field

[0001] This invention relates to the field of hydrological, water resources and climate change response analysis technology, and in particular to a watershed hydrological elastic attribution method based on chain mechanism decomposition. Background Technology

[0002] Climate change and the evolution of underlying surface characteristics have profoundly impacted the hydrological cycle and the spatiotemporal allocation of water resources in watersheds. Quantitatively assessing the physical driving effects of changes in underlying surface elements on watershed runoff is of significant technical value for revealing the mechanisms of hydrological variability, predicting the evolution trends of extreme hydrological events, and supporting the assessment of watershed water resource carrying capacity. Accurately identifying the independent contributions of different hydrological driving factors is a crucial prerequisite for constructing high-fidelity hydrophysical models and realizing the scientific allocation of water resources in complex watersheds.

[0003] Current hydrological response analyses of underlying surface changes are mostly based on macroscopic statistical models using long-term meteorological and hydrological observation data. Existing analytical methods typically use single macroscopic regression coefficients or overall coefficients to summarize the overall impact of changes in individual underlying surface elements such as vegetation, snow cover, or land use on runoff. When processing watershed data under complex climatic backgrounds, the main approach is to establish statistical mappings directly between macroscopic underlying surface indicators and water storage capacity or runoff generation efficiency to extract overall hydrological elasticity values. However, when the watershed contains multiple concurrent physical processes such as radiative energy fluctuations and solid-liquid phase transitions, existing methods based on single statistical mappings struggle to effectively distinguish the independent roles of different physical driving mechanisms in the overall response. Furthermore, when faced with limited sample sizes or the fusion of multi-source heterogeneous meteorological data, the calculation results often face high uncertainty risks.

[0004] Existing methods struggle to balance the fine decoupling of physical processes with the rigor of quantitative attribution when analyzing the impact of complex underlying surface evolution on the hydrological cycle. Faced with complex hydrological environments involving the coupling and superposition of multiple physical mechanisms, breaking free from the black-box nature of macroscopic statistical parameters, effectively isolating and quantifying the independent contributions of each internal physical driving channel, while simultaneously ensuring the objectivity and robustness of complex attribution calculations, is a common technical challenge urgently needing to be addressed in the field of hydrological evolution attribution. It is necessary to investigate a method that can improve the accuracy and adaptability of variability response analysis under complex meteorological and hydrological scenarios. Summary of the Invention

[0005] The purpose of this invention is to provide a watershed hydrological elastic attribution method based on chain mechanism decomposition, in order to solve at least one of the above-mentioned problems in the prior art.

[0006] Technical solution: A watershed hydrological elastic attribution method based on chain mechanism decomposition, comprising:

[0007] Obtain meteorological and hydrological sequences, underlying snow cover sequences, and climate state variables for the target watershed;

[0008] Based on meteorological and hydrological sequences and underlying snow cover sequences, mediating mechanism indicators characterizing physical transport pathways are extracted.

[0009] Determine the empirical sensitivity of each mediation mechanism indicator to changes in underlying snow cover sequence;

[0010] By inputting the mediating mechanism indicators and climate state variables into a pre-built watershed parameter model, the partial sensitivity of watershed water storage parameters to each mediating mechanism indicator can be obtained.

[0011] By combining empirical sensitivity, partial sensitivity, and basic environmental scaling factor determined based on meteorological and hydrological sequences, the explicit mechanism elastic components corresponding to each physical transport path are calculated.

[0012] The residual elastic components of the target watershed are calculated based on the explicit mechanism elastic components.

[0013] By aggregating the explicit mechanism elastic components and residual elastic components, the hydrological elastic attribution results of the target watershed for snow cover changes are output.

[0014] Beneficial effects: This invention achieves quantitative decoupling and complete attribution of multiple physical conduction channels, improving the objectivity and accuracy of complex hydrological response analysis. Attached Figure Description

[0015] Figure 1 This is a schematic diagram of the overall process of a watershed hydrological elastic attribution method based on chain mechanism decomposition provided in an embodiment of this application.

[0016] Figure 2 This is a schematic diagram of the process for obtaining the water storage parameter sequence of the corresponding target watershed provided in the embodiments of this application.

[0017] Figure 3 This is a schematic diagram of the empirical sensitivity process for determining the sensitivity of various mediating mechanism indicators to changes in the underlying snow cover sequence, provided in an embodiment of this application.

[0018] Figure 4 This is a schematic diagram of the process for calculating the residual elastic components of the target watershed based on the explicit mechanism elastic components provided in the embodiments of this application.

[0019] Figure 5 This is a schematic diagram of the uncertainty quantification process for each explicit mechanism elastic component provided in the embodiments of this application. Detailed Implementation

[0020] Example 1: A general flowchart of a watershed hydrological elastic attribution method based on chain mechanism decomposition is provided, as follows: Figure 1 As shown, the main steps include the following:

[0021] Step 101: Obtain the meteorological and hydrological sequence, underlying snow cover sequence, and climate state variables of the target watershed.

[0022] In this embodiment, the target watershed typically refers to a high-altitude or mid-to-high latitude watershed with glacial and hydrological processes. The meteorological and hydrological series encompasses long-term, continuous daily or monthly observation data, specifically including but not limited to precipitation series, runoff series, potential evapotranspiration series, and net radiation series. The underlying snow cover series is a continuous temporal record reflecting the proportion of the watershed surface covered by snow, typically based on remote sensing inversion products, such as MODIS snow cover products. Climate state variables, for example, are represented by multi-year average drought indices, used to characterize the long-term average dry and wet climate background of the watershed. The aforementioned multi-source data provides comprehensive data support for subsequent analysis of the complex physical mechanisms by which underlying surface changes affect runoff.

[0023] In this embodiment, the climate state variable is used to characterize the long-term average climate or underlying surface baseline state of the watershed. In one implementation, the climate state variable can be implemented using a multi-year average drought index; in another implementation, the climate state variable can also be implemented using a long-term average snow cover percentage characterizing the underlying snow cover baseline level. It is understood that different physical transport pathways can select the state regulation quantity that best matches their physical mechanism as the corresponding climate state variable.

[0024] In one optional implementation, the state regulation variables corresponding to different physical transfer paths can be different. For energy transfer paths and hydrothermal time-series correlation paths, the multi-year average drought index Φ can be used as the state regulation variable to characterize the modulation effect of climate dryness / wetness background on their sensitivity. For water input intensity paths, the long-term average snow cover ratio fs of the watershed can be used as the state regulation variable to characterize the modulation characteristics of snow cover degree on the liquid water input intensity effect. Those skilled in the art can select appropriate state regulation variables based on the actual watershed characteristics.

[0025] Step 102: Based on meteorological and hydrological sequences and underlying snow cover sequences, extract mediating mechanism indicators that characterize physical transfer pathways.

[0026] Specifically, existing technologies often treat the impact of snow cover changes on runoff as a single, indistinguishable black-box path. To overcome this limitation, this scheme introduces a mechanism mediating layer between underlying surface characteristics and watershed water storage capacity. Based on acquired meteorological and hydrological observation data, a set of index variables is constructed to quantify different physical driving mechanisms. The mediating mechanism indicators correspond to the energy redistribution path caused by snow cover changes, the path of changes in liquid water input intensity, and the synchronous change path of water and energy in the intra-annual temporal distribution. By extracting the mediating mechanism indicators, a single macroscopic impact can be decomposed into multiple independent transmission channels with clear physical meaning.

[0027] Step 103: Determine the empirical sensitivity of each mediation mechanism indicator to changes in the underlying snow cover sequence.

[0028] Furthermore, after obtaining continuous sequences of mediation mechanism indicators and underlying snow cover sequences, statistical regression methods are used to establish the time-series response relationship between the two. Empirical sensitivity refers to the statistical rate of change of a specific mediation mechanism indicator for each unit change in snow cover percentage within a given historical climate context of the watershed. For example, a robust sensitivity slope can be obtained by extracting the mean sequences of adjacent time windows and calculating the differences between the mean sequences through linear regression, filtering out random fluctuations in interannual climate.

[0029] Step 104: Input the mediating mechanism indicators and climate state variables into the pre-constructed watershed parameter model to obtain the partial sensitivity of watershed water storage parameters to each mediating mechanism indicator.

[0030] Based on this, the pre-built watershed parameter model is a nonlinear estimation model trained offline. This model not only includes the main effect terms of each mediating mechanism indicator, but also introduces the interaction term of the product between the mediating mechanism indicators and the climate state variables.

[0031] After inputting the current watershed data into the model, the calculated partial sensitivity reflects the partial derivative response of the watershed storage parameters caused by a unit change in a single mediating mechanism indicator, while keeping other conditions constant. Introducing climate state variables for interactive regulation ensures that the partial sensitivity exhibits continuous spatial heterogeneity across watersheds with varying degrees of aridity and wetness, consistent with the regional differentiation patterns in hydrology.

[0032] Step 105: Combine empirical sensitivity, partial sensitivity, and basic environmental scaling factor determined based on meteorological and hydrological sequences and underlying snow cover sequences to calculate the explicit mechanism elastic components of each physical transport path.

[0033] A chain-like reconstruction logic is used to aggregate the aforementioned separated sensitivity features. The basic environmental scaling factor is used to characterize the comprehensive scaling coefficient of runoff elasticity when the watershed water storage parameters change by a unit, calculated through the evapotranspiration change path. It is a constant determined by the watershed's long-term hydrothermal balance ground state and the snow cover runoff mean state.

[0034] According to the chain rule of differentiation, for any explicit physical transmission path, multiplying its corresponding empirical sensitivity and partial sensitivity by the scaling factor of the underlying environment yields the independent contribution of that single physical path to the total elasticity, i.e., the explicit mechanism elastic component. This decoupled calculation of the above product structure enables the quantitative separation of complex physical mechanisms.

[0035] Step 106: Calculate the residual elastic components of the target watershed based on the explicit mechanism elastic components.

[0036] The extracted mediating mechanism indicators may not fully cover all underlying surface factors affected by snow cover changes, such as permafrost degradation and vegetation succession. It is necessary to perform fallback quantification on implicit path features not explicitly captured. Based on the obtained explicit mechanism elastic components, by comparing with the overall hydrological response trend of the watershed, the remaining contribution unexplained by the model, i.e., the residual elastic component, is separated. The calculation of this component ensures the mathematical completeness of the entire mechanism decomposition framework.

[0037] Step 107: Aggregate the explicit mechanism elastic components and residual elastic components, and output the hydrological elastic attribution results of the target watershed for snow cover changes.

[0038] All calculated explicit and residual elastic components are aggregated. This aggregation not only includes the overall elasticity values ​​of snow cover changes on watershed runoff but also outputs a quantitative attribution map, clearly indicating the percentage contribution of each independent transmission path, such as energy, water intensity, and temporal sequence, to the total impact. The resulting hydrological elasticity attribution results provide a refined physical mechanism basis for assessing water resource carrying capacity and scheduling water conservancy projects in high-altitude cold watersheds.

[0039] Example 2 describes the underlying mathematical derivation and numerical solution process for obtaining the watershed water storage parameter sequence and determining the basic environmental scaling factor.

[0040] In one possible implementation, when obtaining the water storage parameter sequence for the corresponding target watershed, such as Figure 2 As shown, it specifically includes:

[0041] Step 201: Extract historical precipitation, runoff, and potential evapotranspiration data from the meteorological and hydrological sequence.

[0042] In this embodiment, the meteorological and hydrological sequence serves as the input basis for calculating the water and energy balance of the entire watershed. Specifically, historical precipitation observations are typically obtained from actual measurements by rain gauges or meteorological stations distributed within the target watershed, or extracted based on gridded meteorological reanalysis data. Historical runoff observations represent the actual water yield at the watershed outlet section. Historical potential evapotranspiration observations are usually calculated using standard meteorological data combined with the Penman-Montes formula, used to characterize the maximum possible evaporation under conditions of sufficient atmospheric moisture supply. Obtaining the above-mentioned basic data sequence provides the necessary hydrological boundary conditions for subsequent inversion of underlying surface characteristic parameters.

[0043] Step 202: Using historical precipitation and runoff data, calculate the corresponding historical actual evapotranspiration.

[0044] Furthermore, based on the principle of long-term closed-loop water balance at the watershed scale, water consumption within the watershed can be inferred. Specifically, without considering inter-watershed water transfer from deep groundwater and relatively small interannual variations in water storage, the difference between the input historical precipitation observations and the output historical runoff observations can be used as the actual historical evapotranspiration. This calculation process eliminates the complex intermediate runoff generation and collection links, directly obtaining the total actual water evaporation and loss at the watershed scale.

[0045] In some alternative implementations, if there are known large-scale water diversion projects or significant groundwater extraction activities in the target watershed, additional correction terms can be introduced during the estimation process to compensate for the above-mentioned differences, so as to improve the estimation accuracy of the actual evapotranspiration over the years.

[0046] Step 203: Take the historical precipitation observations, historical potential evapotranspiration observations, and the estimated historical actual evapotranspiration as known constants, and substitute them into the target water-heat balance transcendental equation to construct the iterative objective function.

[0047] In this embodiment, the target hydrothermal balance transcendental equation preferably adopts the parameterized Boudico equation. This equation establishes a nonlinear constraint relationship between actual evapotranspiration, precipitation, and potential evapotranspiration. The equation includes an unknown that reflects the comprehensive characteristics of the watershed soil, vegetation, and topography, namely the watershed water storage parameter. This watershed water storage parameter is usually located in a power-law position, making it difficult to obtain an analytical solution from the equation through algebraic transformations. Therefore, iterative objective functions need to be constructed using this watershed water storage parameter as the independent variable.

[0048] g(n) = 1 + Φ - (1 + Φ n ) (1 / n) -E / P;

[0049] Where g(n) is the iterative objective function, n is the water storage parameter of the basin, Φ is the drought index, which is specifically calculated as the ratio of the observed potential evapotranspiration to the observed precipitation over the years, E is the actual evapotranspiration over the years, and P is the observed precipitation over the years. In the following text, when P is not marked with a year, it refers to the long-term average precipitation by default.

[0050] By constructing the above iterative objective function, the original complex inversion problem is transformed into the problem of finding the zeros of the function.

[0051] Step 204: The Newton-Raphson numerical iteration algorithm is used to differentiate and successively approximate the objective function until the preset convergence criterion is met, and the actual water storage parameters of the corresponding years are retrieved to form a water storage parameter sequence.

[0052] Specifically, the Newton-Raphson numerical iterative algorithm uses the first derivative of the function at the current estimation point to guide the next search direction. An initial value for the watershed storage parameter is set, for example, 2.0. The value of the objective function and its first derivative at the estimation point are calculated, and the estimation point is updated using the iterative formula. The specific iterative formula is expressed as follows:

[0053] n k+1 =n k -g(n k ) / g'(n k );

[0054] Where, n k+1 Let n be the estimated watershed water storage parameters for the (k+1)th iteration. k Let g(n) be the estimated water storage parameter for the k-th iteration. k Let g'(n) be the value of the iterative objective function at the k-th estimation point. k ) is the first derivative of the iterative objective function at the k-th estimation point.

[0055] Specifically, g'(n k Let be the first derivative of the iterative objective function g(n) at the current estimation point. Taking the derivative of the Boudico transcendental equation with respect to n, we obtain:

[0056] g'(n)=(1 / P)*Ω func (n,Φ,P);

[0057] Among them, Ω func The partial derivative function obtained in step (2) of this method is given by P, where P represents the historical precipitation observation value. The relationship between the two indicates that the calculation of the watershed water storage parameter inversion process and the basic environmental scaling factor share the same analytical derivative result, and Φ represents the drought index of the current year. The first derivative value can be calculated based on the above analytical derivative formula or by using the numerical difference method.

[0058] Continue the iterative update process described above, monitoring the absolute difference between two consecutive estimates. Stop the iteration when this absolute difference is less than a preset convergence criterion. The convergence criterion can be set to 10. -6 At this point, the latest estimate is taken as the actual watershed water storage parameter for that year. The above operation is repeated for each year in the observation sequence, eliminating outlier years that do not meet the physical boundary of water-heat balance, thus compiling a continuous water storage parameter sequence. This numerical iterative algorithm possesses quadratic convergence characteristics and can effectively extract underlying surface feature parameters hidden in macroscopic hydrological observation data.

[0059] In one possible implementation, determining the underlying environment scaling factor includes:

[0060] Step (1) extract the long-term average precipitation and potential evapotranspiration from the meteorological and hydrological series, and calculate the drought characteristic constant of the target watershed using the water-heat balance equation. The drought characteristic constant is the drought index.

[0061] Based on this, in addition to retrieving water storage parameters year by year, it is also necessary to determine the climate baseline state of the entire target watershed. Specifically, multi-year averages of precipitation and potential evapotranspiration data in the meteorological and hydrological series are calculated to extract stable long-term average precipitation and long-term average potential evapotranspiration. The ratio of long-term average potential evapotranspiration to long-term average precipitation is calculated to obtain the drought characteristic constant. This constant characterizes the overall dry and wet environment of the target watershed on a multi-year scale, determines the baseline position of the watershed's hydrological response on the Budico curve, and is a dimensionless static climate variable.

[0062] Step (2): Solve the partial derivative function of the target actual evapotranspiration with respect to the watershed water storage parameter for the basic equation of watershed water heat balance that characterizes the balance relationship between precipitation, actual evapotranspiration and potential evapotranspiration.

[0063] Furthermore, to establish an analytical relationship between changes in underlying surface area and changes in evapotranspiration, calculus calculations were performed on the fundamental equations of the watershed's water-heat balance. The target actual evapotranspiration was expressed as a composite function of long-term average precipitation, drought characteristic constant, and watershed water storage parameters.

[0064] Using the chain rule and the logarithmic rule, the first-order partial derivatives of the target actual evapotranspiration with respect to the watershed storage parameters are obtained, yielding the partial derivative function. The specific calculation formula is as follows:

[0065] Ω func (n,Φ,P)=P*(1+Φ n ) (1 / n) *[ln(1+Φ n ) / n 2 -(Φ n*ln(Φ)) / (n*(1+Φ n ))];

[0066] Among them, Ω func (n,Φ,P) is the partial derivative function obtained by solving, P is the long-term average precipitation, Φ is the drought characteristic constant, n is the watershed water storage parameter, and ln() is the natural logarithm function.

[0067] The analytical differentiation process described above makes the sensitivity relationship implicit in the nonlinear equation explicit, providing a computational tool for quantifying the water redistribution caused by changes in the underlying surface.

[0068] Step (3) Substitute the long-term average precipitation and drought characteristic constant into the partial derivative function to calculate the absolute change in the actual evapotranspiration of the target when the water storage parameters of the watershed change by a unit; combine the long-term average snow cover ratio and long-term average annual runoff of the target watershed to perform runoff elastic scale conversion on the absolute change to obtain the basic environmental scaling factor.

[0069] According to one aspect of this application, the step may also involve substituting the long-term average precipitation, drought characteristic constant, and multi-year average watershed water storage parameters determined based on meteorological and hydrological sequences into the partial derivative function to calculate the absolute change in the target actual evapotranspiration when the watershed water storage parameters change by a unit; and performing runoff elastic scale conversion by combining the long-term average snow cover ratio and long-term average annual runoff of the target watershed to obtain the basic environmental scaling factor.

[0070] The long-term average precipitation, drought characteristic constant, and multi-year average watershed water storage parameters extracted in the previous steps are used as input variables and substituted into the partial derivative function obtained from the analytical solution above for numerical calculation. The calculated output is the partial derivative value Ω of the target actual evapotranspiration with respect to the watershed water storage parameters. i .

[0071] Specifically, the absolute change Ω of the actual evapotranspiration is calculated when the watershed storage parameters change by a unit. i Ω i The basic environmental scaling factor is obtained by combining the long-term average snow cover percentage and long-term average annual runoff of the target watershed according to the following formula:

[0072] A i =-(fs i / Q i )×Ω i ;

[0073] Among them, A i Based on the environment scaling factor, fs i Q represents the long-term average snow cover percentage of the target watershed. i Ω represents the long-term average annual runoff of the target watershed.i Let fs be the partial derivative of the target actual evapotranspiration with respect to the watershed storage parameters, i.e., the absolute change in actual evapotranspiration when the watershed storage parameters change by a unit. i Q i Ω i All are positive values, A i The value is negative, and its absolute value reflects the runoff elasticity amplitude transformed by the unit change in watershed water storage parameters through its influence on evapotranspiration.

[0074] The partial derivative Ωi reflects the specific depth at which, under the current predetermined watershed climate conditions, the actual evapotranspiration will increase or decrease by one unit for every increase or decrease in the watershed's water storage parameter. For a real watershed, increased water storage capacity leads to greater water evaporation.

[0075] Based on this, the basic environment scaling factor A i The partial derivative value Ω i Coupling with the long-term homogeneous characteristics of watershed snow cover and runoff is equivalent to determining the core transformation center in the entire chain-like elastic decomposition architecture, so that the sensitivity of intermediate mechanisms in different dimensions can be uniformly converted into the final runoff elastic contribution.

[0076] Example 3 describes the specific process of extracting mediating mechanism indicators to characterize multiple independent physical transport paths and obtaining the phase characteristics of bottom precipitation.

[0077] In one possible implementation, the physical transfer path includes the energy transfer path, the water input intensity path, and the hydrothermal temporal correlation path; correspondingly, the mediating mechanism indicators include the energy mechanism indicator reflecting the characteristics of net radiation during the snow season, the liquid water input intensity indicator reflecting the characteristics of the average daily liquid water input during the year, and the hydrothermal asynchrony indicator reflecting the characteristics of the time misalignment between energy supply and water supply.

[0078] Specifically, when the proportion of snow cover changes, the hydrological response within the target watershed exhibits a multi-dimensional physical transmission mechanism.

[0079] For energy transfer pathways, an energy mechanism index is extracted to quantify the change in net radiation flux caused by changes in surface albedo. This index is obtained based on the calculation of the arithmetic mean of daily net radiation over a continuous period when the snow water equivalent is greater than 0.

[0080] I E =1 / N snow *∑(R n (t));

[0081] Among them, I E As an indicator of energy mechanisms, N snowR represents the number of days in the snow season, i.e., the total number of consecutive days in which the snow water equivalent is greater than 0. ∑() is a function that performs a summation operation on the data for each day within the snow season. n (t) represents the net surface radiation on day t.

[0082] Based on the water input intensity path, liquid water input intensity index is extracted to quantify the average daily input of liquid water composed of rainfall and snowmelt.

[0083] I LWI =V LWI / D LWI ;

[0084] Among them, I LWI V is the input intensity index for liquid water. LWI D represents the total annual input of liquid water. LWI This refers to the number of days with liquid water input per year.

[0085] Annual number of days of liquid water input D LWI The calculation is based on the criterion that the corresponding daily liquid water input LWI(t) is greater than a preset minimum effective daily input threshold. The minimum effective daily input threshold can be determined according to the climate zone of the target watershed and the accuracy of the observation data. In an optional implementation, the minimum effective daily input threshold can be set to 0.1 mm. The specific value of this threshold can be adaptively adjusted according to the precipitation characteristics of the actual watershed.

[0086] For the hydrothermal time-series correlation path, the hydrothermal asynchrony index is extracted to measure the degree of temporal deviation between the peak time of water supply and the peak time of energy supply within the year.

[0087] The above three indicators constitute the mediating mechanism layer connecting the underlying surface characteristic parameters and watershed runoff changes.

[0088] In one possible implementation, acquiring the liquid water input features corresponding to the path characterizing the moisture input intensity specifically includes:

[0089] Step 301: Obtain the daily precipitation sequence and daily average temperature sequence from the meteorological and hydrological sequence.

[0090] In this embodiment, daily precipitation and daily average temperature sequences are the fundamental input variables for watershed precipitation phase analysis. Specific methods for obtaining these sequences include extracting daily measured records from a ground-based meteorological observation network, or processing meteorological reanalysis datasets using spatial interpolation algorithms to generate gridded time series covering the target watershed. The temporal resolution of these sequences is uniformly configured to a daily scale to capture the transient characteristics of precipitation events.

[0091] Step 302: Configure the precipitation and snowfall thresholds for precipitation temperature, and perform proportional mapping based on the daily average temperature sequence and the precipitation and snowfall thresholds to divide the daily precipitation sequence into daily precipitation sequence and daily snowfall sequence.

[0092] Specifically, a dual-threshold temperature method is used to separate the solid and liquid components of mixed precipitation. The preset snowfall threshold can be set to 0 degrees Celsius, and the preset rainfall threshold can be set to 2 degrees Celsius. The corresponding rainfall proportion coefficient is calculated based on the specific position of the daily average temperature between these two thresholds. When the daily average temperature is less than or equal to the snowfall threshold, the corresponding rainfall proportion coefficient is assigned a value of 0. When the daily average temperature is greater than or equal to the rainfall threshold, the corresponding rainfall proportion coefficient is assigned a value of 1. When the daily average temperature is between the snowfall and rainfall thresholds, a linear interpolation model is used to calculate the rainfall proportion coefficient.

[0093] f rain (t)=(T(t)-T snow ) / (T rain -T snow );

[0094] Among them, f rain (t) is the rainfall proportion coefficient on day t, T(t) is the average daily temperature on day t, and T snow T is the snowfall threshold. rain This is the rainfall threshold.

[0095] Based on the extracted rainfall proportion coefficient, the daily precipitation series is decomposed into a time series of two independent variables.

[0096] P rain (t)=f rain (t)*P(t);

[0097] Among them, P rain P(t) represents the daily rainfall on day t, and P(t) represents the daily precipitation on day t.

[0098] P snow (t)=[1-f rain [(t)]*P(t);

[0099] Among them, P snow (t) represents the daily snowfall on day t.

[0100] Step 303: Obtain the snow water equivalent observation sequence of the target watershed, perform daily mass balance calculation based on the daily snowfall sequence and the snow water equivalent observation sequence, and extract the daily snow melt sequence that represents the amount of snow melting.

[0101] Furthermore, solid snowfall that falls to the ground forms snow storage within the watershed, and its melting process relies on continuous mass conservation calculations. The snow water equivalent observation sequence is extracted from multi-source fused hydrological data and combined with the previously generated daily snowfall sequence to calculate the daily meltwater output. Considering that sublimation loss or observation equipment errors may cause negative values ​​in physical quantities, a maximum value function is used to non-negatively truncate the output results. The specific mass balance formula is:

[0102] M(t) = max(0, SWE(t-1) + P snow (t)-SWE(t));

[0103] Where M(t) is the daily snowmelt amount on day t, max() is the maximum value function, SWE(t-1) is the snow water equivalent of the previous day, and P snow SWE(t) represents the daily snowfall on day t, and SWE(t) represents the snow water equivalent on that day.

[0104] Through the above mass balance calculations, the accurate time points and scale characteristics of the release of liquid water from melting snow were extracted.

[0105] Step 304: Merge the daily rainfall sequence and the daily snowmelt sequence to obtain the daily average liquid water input sequence used to calculate the characteristic index corresponding to the path of water input intensity.

[0106] According to one aspect of this application, the step may also involve merging the daily rainfall sequence and the daily snowmelt sequence to obtain a daily liquid water input sequence; based on the daily liquid water input sequence, calculating the total annual liquid water input and the number of days with annual liquid water input, and calculating the liquid water input intensity index corresponding to the water input intensity path.

[0107] After obtaining the two liquid water sources, direct rainfall and snowmelt, the two are added and aggregated on the same time dimension.

[0108] LWI(t)=P rain (t)+M(t);

[0109] Where LWI(t) is the amount of liquid water input on day t.

[0110] This liquid water input sequence forms the basis for calculating the total annual liquid water input and determining the number of days with annual liquid water input. Based on this sequence, the aforementioned liquid water input intensity index is calculated, providing a computational framework for quantifying the inhibitory characteristics of water input concentration on runoff efficiency.

[0111] Example 4 provides a specific quantitative extraction method for hydrothermal asynchrony index characterizing the time misalignment of energy supply and water supply, and provides a basic implementation method based on the difference of a single percentile for hydrothermal correlation paths and an optimized implementation method based on the difference of the full distribution morphology.

[0112] One possible implementation involves determining the hydrothermal asynchrony index, including the following steps:

[0113] Step 401: Obtain the liquid water input sequence and net surface radiation sequence of the target watershed based on meteorological and hydrological sequences. Specifically, obtain the liquid water input sequence and net surface radiation sequence of the target watershed based on meteorological and hydrological sequences and underlying snow cover sequences.

[0114] In this embodiment, the liquid water input sequence can be obtained by adding and merging the rainfall sequence and snowmelt sequence within a predetermined time scale. The net surface radiation sequence is directly extracted based on radiation flux monitoring data from meteorological and hydrological sequences. The above two sets of sequences provide the basic time-series variables required to quantify the temporal matching degree of water input and energy input.

[0115] In some specific embodiments, a basic implementation method is provided for calculating the hydrothermal asynchrony index based on the time difference of a single cumulative midpoint, and the following steps are performed accordingly.

[0116] Step 402: Determine the first date when the cumulative amount of liquid water input sequence reaches half of the total annual input, and the second date when the cumulative amount of surface net radiation sequence reaches half of the total annual radiation.

[0117] Specifically, the daily average liquid water input sequence is calculated by accumulating the data daily in ascending order within the year. The day number at which the accumulated value first exceeds or equals 50% of the total liquid water input for the corresponding year is recorded, and this day number is assigned as the first date. Similarly, the same daily ascending accumulation operation is performed on the daily net surface radiation sequence. The day number at which the accumulated value first exceeds or equals 50% of the total radiation for the corresponding year is recorded, and this day number is assigned as the second date. The values ​​of the first and second dates are integers ranging from 1 to 365.

[0118] Step 403: Calculate the absolute time difference between the first date and the second date as an indicator of water-thermal asynchrony.

[0119] After obtaining the median dates of the two dates representing 50% cumulative progress, perform the difference operation and extract the absolute value of the result.

[0120] I τ =|t LWI_50 -t Rn_50 |;

[0121] Among them, Iτ As an index of hydrothermal asynchrony, t LWI_50 For the first date, t Rn_50 For the second date, the symbol || represents the absolute value operation.

[0122] The physical basis of the aforementioned implementation method is established under predetermined climatic constraints. For most mid-to-high latitude snow-covered watersheds in the Northern Hemisphere, the concentrated water input period driven by snowmelt usually precedes the summer period when net radiation energy peaks. Using absolute time difference calculations can intuitively quantify the time-series misalignment between peak water and energy supply. However, this calculation method, which relies on single percentile differences, lacks sensitivity to the overall distribution characteristics of the indicator sequence and may lead to misjudgments when dealing with complex watersheds with multi-peak distributions.

[0123] Based on the above embodiments, an alternative implementation method is provided as a preferred method for calculating the hydrothermal asynchrony index, to overcome the limitations of a single percentile metric, and the following steps are performed accordingly:

[0124] Step 402a: Construct the first normalized cumulative distribution function of the liquid water input sequence on the time axis and the second normalized cumulative distribution function of the net surface radiation sequence on the time axis.

[0125] For each target year within the target watershed observation sequence, the proportion of the cumulative liquid water input from the first day of each year to the current time point is calculated daily to the total liquid water input for that year. This set of proportions is used to form the first normalized cumulative distribution function.

[0126] F LWI (t)=∑ d=1 t (LWI(d)) / ∑ d=1 365 (LWI(d));

[0127] Among them, F LWI (t) represents the value of the first normalized cumulative distribution function on day t, LWI(d) represents the liquid water input on day d, and ∑ d=1 t () represents the summation operation from the starting point of time to day t, ∑ d=1 365 () indicates a summation operation covering the entire year.

[0128] Furthermore, when processing radiation energy parameters, a condition determination operation is applied to the daily net surface radiation data, retaining and accumulating only the positive values ​​to eliminate radiation deficit interference during non-evaporation-driven periods.

[0129] F Rn (t)=∑d=1 t (max(R n (d),0)) / ∑ d=1 365 (max(R n (d),0));

[0130] Among them, F Rn (t) represents the value of the second normalized cumulative distribution function on day t, R n (d) represents the net surface radiation on day d, and max() is the conditional operation for taking the maximum value, used to constrain the lower limit of the input to 0. The range of the two cumulative distribution functions mentioned above is strictly limited to the closed interval between 0 and 1.

[0131] Step 403a: Extract the deviation distribution values ​​of the first normalized cumulative distribution function and the second normalized cumulative distribution function at each corresponding time node, calculate the sum of the deviation distribution values ​​for the entire time period as the distribution integral distance, and use the distribution integral distance as the hydrothermal asynchrony index.

[0132] After obtaining two normalized cumulative distribution functions, for each discrete day node on the annual timeline, the absolute difference between the corresponding values ​​of the two functions is extracted, and this absolute difference is defined as the deviation distribution value. The deviation distribution values ​​extracted from all time nodes throughout the year are summed to obtain the parameter value equivalent to the discretized one-dimensional Wasserstein distance, i.e., the distribution integral distance.

[0133] I τ_W =∑ t=1 365 (|F LWI (t)-F Rn (t)|);

[0134] Among them, I τ_W F is a hydrothermal asynchrony index in the form of distributed integral distance. LWI (t) represents the value of the first normalized cumulative distribution function on day t, F Rn (t) represents the value of the second normalized cumulative distribution function on day t, ∑ t=1 365 () indicates that the values ​​at each time point throughout the year are summed.

[0135] This preferred implementation is equivalent to the geometric area integral enclosed by the water accumulation curve and the energy accumulation curve on a two-dimensional time plane. When the liquid water input process and the net surface radiation process evolve completely synchronously within the year, the calculated result of the distribution integral distance is 0.

[0136] The integral distance, which incorporates all distribution morphology information, improves the accuracy of the index in recognizing multimodal input pulses and complex sequence widths. It can more reliably assess the degree of reduction in evapotranspiration efficiency due to hydrothermal misalignment in mid-latitude mixed winter-rainfall-spring-melting watersheds.

[0137] Example 5 mainly describes the specific implementation process of statistically extracting the empirical sensitivity of each indicator to snow cover changes at each watershed scale after obtaining multidimensional mechanism indicators.

[0138] One possible implementation involves determining the empirical sensitivity of each mediation mechanism indicator to changes in the underlying snow cover sequence, such as... Figure 3 As shown, it includes the following steps:

[0139] Step 501: According to the preset window length and sliding step size, the underlying snow accumulation sequence and each mediation mechanism indicator are divided into time series, and multiple overlapping time windows are extracted.

[0140] In this embodiment, a sliding window technique is applied to the time series data of each target watershed to construct an analysis sample. Specifically, the effective data year span for the target watershed is set to the first and last year, the preset window length is set to the predetermined number of years for smoothing interannual climate fluctuations, and the preset sliding step size is set to the number of years the window slides.

[0141] The preset window length should meet the following constraints: it should be long enough to smooth interannual climate random fluctuations, while ensuring that the number of differential samples meets the minimum significance requirements for statistical regression. The specific value can be determined through conventional statistical significance tests based on the length of the effective observation sequence and the climate variability characteristics of the target watershed. In this embodiment, the preset window length is set to W years, and the preset sliding step size is set to 1 year.

[0142] This sequence partitioning method discretizes continuous long-term observation sequences into a set of multiple overlapping subsequences containing temporal evolution information.

[0143] In one embodiment, the total number of overlapping time windows extracted based on the above parameters can be calculated by dividing the difference between the total number of years and the window length by the sliding step size, rounding down, and then adding 1.

[0144] Step 502: Calculate the mean data within each overlapping time window. By performing difference processing on the mean data of adjacent overlapping time windows, obtain the index sliding difference sequence of each mediation mechanism indicator and the snow accumulation sliding difference sequence of the underlying snow accumulation sequence.

[0145] Furthermore, within each extracted overlapping time window, an arithmetic mean is calculated for the data covering all years within that window. Taking any mediation mechanism indicator as an example, the average level of the multi-year observations of that indicator within the current time window is calculated.

[0146] I k_ij =1 / W*∑ y (I k_iy );

[0147] Among them, I k_ij Let W be the mean of the mediation mechanism index for the i-th target watershed within the j-th time window, and W be the preset window length. y () indicates that the summation operation is performed on the data of each year contained in the j-th time window, I k_iy This represents the actual value of the intermediary mechanism indicator for the corresponding year.

[0148] After calculating the internal mean of all sliding windows, the mean of two adjacent overlapping time windows on the time axis is subtracted.

[0149] Δ I_k_ij =I k_i (j+1)-I k_ij ;

[0150] Where, Δ I_k_ij For a single difference sample in the corresponding index moving difference sequence, I k_i (j+1) represents the mean of the data in the (j+1)th time window, I k_ij Let be the mean of the data in the j-th time window. Using the same mean calculation and difference processing logic for adjacent windows, the underlying snow cover ratio data are processed synchronously to extract the snow cover slip difference sequence. The difference processing constrains the analysis object to the dimension of the relative rate of change of variables.

[0151] Step 503: Regress the snow cover sliding difference sequence using the index sliding difference sequence of each mediation mechanism indicator to determine the empirical sensitivity of each mediation mechanism indicator.

[0152] After obtaining multiple sets of corresponding difference samples, for each target watershed, a linear regression equation is constructed using the snow cover sliding difference sequence as the independent variable and the index sliding difference sequences of each of the mediation mechanism indicators as the dependent variable.

[0153] Δ I_k_ij =c k_i *Δ fs_ij +d k_i +η k_ij ;

[0154] Where, Δ I_k_ijThe difference values ​​of the moving difference sequence of the index for the target watershed, Δ fs_ij Let c be the difference value of the snow slide difference sequence for the target watershed. k_i d is the empirical sensitivity coefficient of the independent variable to the dependent variable. k_i η is the intercept term of the linear regression equation. k_ij This is the random error term.

[0155] The unknown parameters in the above linear regression equation are estimated and solved using the least squares algorithm. The empirical sensitivity coefficient obtained from the solution is the empirical sensitivity of the corresponding target watershed.

[0156] Specifically, for the energy mechanism index, the estimated empirical sensitivity characterizes the specific change in net radiation during the snow season for each unit change in the underlying snow cover sequence, in units of W / m². 2 For the liquid water input intensity index, the estimated empirical sensitivity characterizes the specific change in liquid water input intensity per unit change in the underlying snow cover sequence, in mm / day. For the hydrothermal asynchrony index, the estimated empirical sensitivity characterizes the specific change in the number of asynchronous hydrothermal days per unit change in the underlying snow cover sequence, in days. This regression fitting process condenses the spatiotemporal response relationship into a deterministic statistical constant characterizing historical climate driving features.

[0157] Example 6 describes the construction of nonlinear watershed parameter models and the method for extrapolating partial sensitivity.

[0158] In one possible implementation, the watershed parameter model is pre-built through the following offline steps:

[0159] Step 601: Obtain historical sample sets for multiple watersheds. The historical sample sets include pre-extracted watershed water storage parameter change samples, as well as corresponding mediation mechanism indicator change samples and climate state variable samples.

[0160] In this embodiment, the offline training process requires the collection of a large amount of spatially heterogeneous watershed data. By applying a sliding window operation to long-term observation data of multiple known watersheds, the difference sequence between adjacent windows is extracted as a variation sample.

[0161] Building upon this, to ensure the adequacy of the model during training and to prevent unmodeled confounding factors from interfering with the regression coefficients, the historical sample set, in addition to including samples of changes in the aforementioned mediating mechanism indicators, also incorporates samples of changes in leaf area index and precipitation intensity as control variables. A complete nonlinear watershed parameter model equation is then constructed using these multivariate samples.

[0162] Δ n =β _0 +β E *ΔI_E +β LWI *Δ I_LWI +β τ *Δ I_τ +β LAI *Δ LAI +β Pi *Δ Pi +β E_Φ *Δ I_E *Φ+β LWI_fs *Δ I_LWI *fs+β τ_Φ *Δ I_τ *Φ+ε;

[0163] Where, Δ n For the sample of watershed water storage parameter changes, β _0 For the intercept term, β E For the univariate main effect characteristics of the energy mechanism, Δ I_E For the sample of changes in energy mechanism indicators, β LWI For the univariate main effect characteristics of the intensity mechanism, Δ I_LWI For the sample of changes in the intensity mechanism index, β τ For the univariate main effect characteristics of the time-series mechanism, Δ I_τ For the time-series mechanism indicator change sample, β LAI Δ is the sensitivity coefficient of leaf area index to watershed water storage parameters. LAI For the sample of leaf area index changes, β Pi Δ is the sensitivity coefficient of precipitation intensity to watershed water storage parameters. Pi For the sample of precipitation intensity variation, β E_Φ Φ is the interaction coefficient between energy and arid climate state, Φ is a climate state variable in the form of an aridity index, and β is... LWI_fs β is the interaction coefficient between intensity and snow cover percentage, fs is the long-term average snow cover percentage sequence sample of the watershed, and β τ_Φ ε is the interaction coefficient between the time series and the arid climate state, and ε is the model residual.

[0164] The above formula introduces multi-dimensional parameter interaction features through the product term, enabling the model to handle spatial heterogeneity.

[0165] Step 602: Based on the regularized regression algorithm, feature screening is performed on the historical sample set to determine a subset of significant predictive variables.

[0166] Specifically, multicollinearity may exist in a candidate variable pool containing multiple interaction terms and control variables. A lasso algorithm is used for the first step of variable selection. This regularized regression algorithm forces the coefficients of redundant features that contribute little to the prediction of the target dependent variable to be compressed to 0 by introducing a penalty term for the absolute value of the regression coefficients into the least squares loss function. After selection, all feature variables with non-zero regression coefficients are extracted, forming a subset of significant predictive variables. This feature selection process effectively reduces the dimensionality of the model, establishing a stable low-dimensional search space for subsequent constrained optimization solutions.

[0167] Step 603: Under the preset physical constraints, namely the preset directional inequality physical constraints, least squares regression fitting is performed using a subset of significant predictive variables to determine the model coefficients of the watershed parameter model.

[0168] After obtaining a significant subset of predictor variables, a restricted regression optimization is performed using a standard quadratic programming solver. To prevent the statistical fitting process from generating counterintuitive parameters that violate hydrophysical common sense, prior knowledge about the direction of the mechanism is converted into equivalent mathematical inequalities and applied to the solution boundary of the objective function.

[0169] The preset physical constraints include directional inequality constraints corresponding to each physical transfer path; the directional inequality constraints are configured as follows:

[0170] Specifically, the partial sensitivity corresponding to the forced energy transfer pathway maintains a positive effect under all climatic conditions to characterize the enhancement effect of energy increase on the watershed's water storage capacity. This constraint requires that the energy partial sensitivity be greater than 0 when the drought index reaches both the minimum and maximum values ​​of the sample.

[0171] β E +β E_Φ *Φ min >0; β E +β E_Φ *Φ max >0;

[0172] Where, β E As a univariate main effect characteristic of the energy mechanism, β E_Φ Φ is the interaction coefficient between energy and arid climate state. min Φ represents the minimum drought index in the historical sample set. max This represents the maximum drought index value in the historical sample set.

[0173] The partial sensitivity corresponding to the forced water input intensity path remains negative at all underlying snow cover levels, characterizing the inhibitory effect of increased input intensity on the watershed's water storage capacity.

[0174] βLWI +β LWI_fs *fs min <0; β LWI +β LWI_fs *fs max <0;

[0175] Where, β LWI For the univariate main effect characteristics of the intensity mechanism, β LWI_fs fs is the interaction coefficient between intensity and snow cover percentage. min fs is the minimum percentage of snow cover in the historical sample set. max This represents the maximum percentage of snow cover in the historical sample set.

[0176] The partial sensitivity corresponding to the forced hydrothermal time-series correlation path remains negative under all climatic conditions, which characterizes the inhibitory effect of hydrothermal time-series misalignment on the watershed's water storage capacity.

[0177] β τ +β τ_Φ *Φ min <0;

[0178] β τ +β τ_Φ *Φ max <0;

[0179] Where, β τ For the univariate main effect characteristics of the time-series mechanism, β τ_Φ This is the interaction coefficient between the time series and the arid climate state.

[0180] By solving for the minimum value of the sum of squared residuals within the strict boundaries of the above set of inequalities, the globally fixed model coefficient values ​​are obtained, and the training and construction of the nonlinear watershed parameter model is completed.

[0181] In one possible embodiment, the watershed parameter model includes univariate main effect features corresponding to each mediating mechanism indicator, as well as parameter interaction features corresponding to the combination of each mediating mechanism indicator and climate state variables.

[0182] In other words, the watershed parameter model includes the univariate main effect characteristics of each mediating mechanism indicator, as well as the parameter interaction characteristics of the combination of each mediating mechanism indicator with the state regulating variable characterizing the watershed climate or underlying surface state.

[0183] One possible implementation, when obtaining the partial sensitivity of watershed storage parameters to various mediating mechanism indicators, specifically includes:

[0184] By using a pre-built watershed parameter model, we can extract the univariate main effect characteristics of each mediating mechanism indicator, as well as the sensitive adjustment magnitude of the interaction characteristics of each parameter under the current climate state variable.

[0185] After the model is trained offline, for the predetermined online simulation watershed, the corresponding climate state variables and underlying surface average characteristic parameters are extracted. The corresponding univariate main effect features are retrieved as base values ​​from the pre-built watershed parameter model. The actual values ​​of the climate state variables of the target watershed are multiplied by the corresponding interaction coefficients to calculate the sensitive moderating magnitude of the parameter interaction characteristics.

[0186] Furthermore, the sensitivity moderating magnitudes of the main effects of each univariate and the corresponding interaction features of each parameter are linearly superimposed to obtain the partial sensitivity of each mediation mechanism indicator.

[0187] The extracted baseline values ​​are summed with the sensitivity adjustment range. Taking the energy mechanism as an example, its univariate main effect characteristics are added to the sensitivity adjustment range under the predetermined drought index to obtain the predetermined effective partial sensitivity of the energy mechanism in the target watershed.

[0188] S E =β E +β E_Φ *Φ;

[0189] Among them, S E For the partial sensitivity corresponding to the energy transfer path, β E This is a univariate main effect characteristic, β E_Φ Φ is the interaction coefficient between energy and arid climate state, and Φ is the climate state variable of the current watershed.

[0190] During the online simulation phase, the mediation mechanism index is used to specify the mechanism channels corresponding to the univariate main effect features and parameter interaction features retrieved in the watershed parameter model, and the climate state variable is used to calculate the sensitive adjustment magnitude of the parameter interaction features. In the implementation method of calculating the residual elastic components using the model fitting residual method, the changes in the mediation mechanism index and the climate state variable in each time window are also used as inputs to the watershed parameter model to generate the corresponding model fitting prediction values.

[0191] The above overlay operation realizes the adaptive smooth mapping of the output parameters of a single global model with regional climate characteristics, giving the derived partial sensitivity spatial resolution and physical rationality.

[0192] Example 7 describes the process of performing chain reconstruction to calculate the explicit mechanism elastic components after obtaining multidimensional sensitivity features, and provides a basic implementation method based on model fitting residual extraction, as well as a preferred implementation method based on closed interpolation method to extract complete residuals.

[0193] In one possible implementation, the calculation of the explicit mechanism elasticity components corresponding to each physical transmission path includes:

[0194] Step 701: For any target physical transfer path among all physical transfer paths, extract the empirical sensitivity and the corresponding partial sensitivity of the mediation mechanism index corresponding to the target physical transfer path in the target watershed.

[0195] In this embodiment, the impact of underlying surface changes on watershed runoff is decomposed into multiple parallel and independent physical transfer paths. For any one of these target physical transfer paths, two key sensitivity dimension features need to be extracted from the preprocessing module.

[0196] Specifically, an empirical sensitivity generated by watershed-by-watershed sliding window regression is extracted. This empirical sensitivity reflects the statistical response slope of a predetermined mediating mechanism indicator to changes in underlying snow cover sequences over historical time. Simultaneously, a partial sensitivity derived from a nonlinear watershed parameter model combined with current climate state variables is extracted. This partial sensitivity characterizes the partial derivative contribution of a single mediating mechanism indicator change to the overall watershed water storage parameter, assuming all other conditions remain frozen. By aggregating these two types of sensitivities from different dimensions, computational factors are provided for constructing a complete chain-like differential transmission logic.

[0197] Step 702: Multiply the empirical sensitivity and partial sensitivity corresponding to the target physical transfer path by the basic environment scaling factor to obtain the explicit mechanism elastic component of the corresponding target physical transfer path.

[0198] Furthermore, based on the chain rule of composite functions in calculus, the extracted local sensitivities are cascaded and multiplied along the physical conduction direction. To map the changes in water storage parameters to the final runoff elasticity dimension, a basic environmental scaling factor needs to be introduced as a common transformation coefficient.

[0199] The empirical sensitivity and partial sensitivity corresponding to the target physical transfer path are continuously multiplied with the scaling factor of the underlying environment. The resulting product is the independent contribution of the target physical transfer path to the total elasticity, i.e., the explicit mechanism elastic component.

[0200] ε fs_k_i =A i *S k_i *c k_i ;

[0201] Where, ε fs_k_i The explicit mechanism elastic component for prescribing the physical transport path for the target watershed, A i S is the basic environmental scaling factor determined by the long-term hydro-climatic ground state of the target watershed. k_i For the partial sensitivity corresponding to this physical transfer path, c k_i This represents the empirical sensitivity corresponding to the physical transmission path.

[0202] For the energy transfer path, the moisture input intensity path, and the hydrothermal temporal correlation path, the above multiplication operation is performed respectively to obtain the energy mechanism elastic component, the intensity mechanism elastic component, and the temporal mechanism elastic component. The above operation process quantitatively decomposes the originally indistinguishable overall elasticity into multiple clear physical conduction contributions.

[0203] In some specific embodiments, a basic implementation method is provided for extracting the additional underlying surface influence based on the fitting residual of the watershed parameter model, which involves performing the following steps.

[0204] When calculating the residual elastic components of the target watershed based on explicit mechanism elastic components, such as Figure 4 As shown, it specifically includes:

[0205] Step S1: Obtain the model fitting prediction value of the watershed parameter model that belongs to the same watershed parameter model as the generated explicit mechanism elastic component, and extract the actual observed value of the parameters of the target watershed.

[0206] In other words, the changes in the intermediate mechanism indicators and climate state variables corresponding to each time window are substituted into the watershed parameter model to obtain the model fitting prediction values ​​that belong to the same watershed parameter model as the generated explicit mechanism elastic components. Based on the meteorological and hydrological series, the actual observed values ​​of parameters of the target watershed are extracted through water balance inversion.

[0207] Specifically, the aforementioned three explicit mediation mechanism indicators may not fully capture all physical processes related to changes in underlying snow cover, such as underlying surface-related factors like permafrost degradation and changes in soil properties that may accompany snow reduction. To quantify the contribution of these omitted factors, the watershed parameter model used to generate the partial sensitivity was retrieved, and the model-fitted predictions output by this model within each historical time window were extracted. Simultaneously, the actual observed values ​​of the parameters obtained from inversion based on actual meteorological and hydrological data were acquired.

[0208] Step S2: Calculate the numerical difference between the actual observed values ​​of each parameter and the corresponding model-fitted predicted values, and generate a serialized fitting residual sequence.

[0209] On the same time dimension, the actual observed values ​​of the parameters for each time window are subtracted from the corresponding model-fitted predicted values. This subtraction operation separates all residual fluctuations not explained by the main characteristics of the watershed parameter model. The residual fluctuations calculated for each time window are then concatenated and recombined in chronological order to form a one-dimensional sequence of fitted residuals.

[0210] Step S3: Determine the residual sensitivity coefficient of the fitted residual sequence to the change in the underlying snow cover sequence.

[0211] After obtaining the fitted residual sequence, a univariate regression equation is constructed between the fitted residual sequence and the underlying snow cover sequence. The slope of the independent variable in this regression equation is estimated using the least squares method, and the extracted slope value is the residual sensitivity coefficient. This coefficient reflects the statistical law that the residual part not explained by the model changes synchronously with the change in the snow cover ratio.

[0212] Step S4: Multiply the residual sensitivity coefficient by the basic environment scaling factor to obtain the residual elastic component of the corresponding target watershed.

[0213] Using the same dimensional transformation logic as for the explicit mechanism, the obtained residual sensitivity coefficient is multiplied by the aforementioned basic environmental scaling factor. The calculated product is the residual elastic component. If the absolute value of this residual elastic component accounts for a small proportion of the total elasticity, it verifies the sufficiency of the aforementioned three explicit physical transfer path decompositions; if it accounts for a large proportion, it indicates the existence of important omitted physical coupling paths in the target watershed.

[0214] Based on the above embodiments, an alternative basic implementation method is provided as a preferred implementation method for calculating the residual elastic components and ensuring mathematical complete closure, so as to overcome the decomposition non-closure error caused by introducing other control variables in the watershed parameter model, and the following steps are performed accordingly.

[0215] Optionally, when calculating the residual elastic components of the target watershed, the following may be included:

[0216] Step S01 involves obtaining the water storage parameter sequence for the corresponding target watershed based on meteorological and hydrological sequences, and determining the overall sensitivity of the water storage parameter sequence to changes in the underlying snow cover sequence. Specifically, this step can also involve obtaining the water storage parameter sequence for the corresponding target watershed based on meteorological and hydrological sequences through annual water balance calculations and iterative solutions to the water-heat balance equations, and determining the overall sensitivity of the water storage parameter sequence to changes in the underlying snow cover sequence.

[0217] Specifically, when nonlinear watershed parameter models include control variables such as leaf area index or precipitation intensity, the contributions of these control variables as intermediaries transmitting the influence of underlying snow accumulation are absorbed by the model. This results in the sum of the four elastic components calculated directly using the fitting residuals not equaling the observed total elasticity, thus compromising the mathematical completeness of the decomposition. To address this issue, the closed-interference method is used to redefine the residuals. A linear regression model is directly established between the actual observed sequences of watershed water storage parameters and the underlying snow accumulation sequences, and the overall total sensitivity slope between the two is calculated. This overall total sensitivity encompasses the total effect of all known and unknown transmission paths.

[0218] Step S02: Extract the product of the empirical sensitivity and the partial sensitivity corresponding to each physical transmission path, and use it as the explicit mechanism attribution contribution of each physical transmission path.

[0219] Furthermore, for all explicit physical transfer pathways analyzed above, intermediate results were obtained by multiplying their empirical sensitivity and partial sensitivity. These intermediate results represent the changes in water storage parameters caused by a single pathway, and are defined as the explicit mechanism attribution contribution.

[0220] Step S03: Subtract the sum of the explicit mechanism attribution contributions of each physical transmission path from the overall total sensitivity, and extract the closure difference residual for the corresponding target watershed.

[0221] The explicit mechanism attribution contributions of all physical transport paths are summed, and this summation is subtracted from the overall sensitivity obtained from the aforementioned calculation. The resulting difference is defined as the closed-loop residual.

[0222] α closure_i =γ i -∑ k (S k_i *c k_i );

[0223] Where, α closure_i γ represents the closure difference residual of the target watershed. i For the overall total sensitivity, ∑ k () indicates that the attribution contribution of each explicit physical transport path is summed, S k_i For partial sensitivity, c k_i It is based on experience sensitivity.

[0224] The closed-loop interpolation residual forcibly absorbs the snow-related portion of the model fit residual, as well as the implicit path contributions passed by all control variables.

[0225] Step S04: Combine the closure difference residual with the basic environment scaling factor to obtain the residual elastic component of the corresponding target watershed.

[0226] The calculated closure difference residuals are multiplied by the basic environmental scaling factor, and a scaling transformation is performed. The result of the transformation is the final residual elastic component.

[0227] The closed-loop residual is defined based on the arithmetic difference. When the elastic components of each explicit mechanism are summed with the residual elastic component, the sum is necessarily strictly equal to the total elasticity of the target watershed. It does not require the assumption of absolute independence between variables and achieves 100% mathematical completeness of the decomposition across all climatic backgrounds. This preferred implementation eliminates systematic decomposition errors in the framework while retaining the residual diagnostic function.

[0228] Example 8 provides a chain error propagation and uncertainty quantification method based on the first-order Taylor expansion principle. That is, after completing the point estimation elasticity decomposition, the multi-stage independent estimation error is propagated to each explicit mechanism elasticity component to output the confidence interval in an advanced evaluation process.

[0229] According to one aspect of this application, such as Figure 5 As shown, the method further includes an uncertainty quantification step corresponding to each explicit mechanism elastic component:

[0230] Step 801: Obtain the first estimated variance of each mediation mechanism indicator during the stage of determining empirical sensitivity, and the second estimated variance during the stage of obtaining partial sensitivity.

[0231] Specifically, the first estimated variance of each mediation mechanism indicator is obtained in the regression fitting process for determining empirical sensitivity, and the second estimated variance is obtained in the regression fitting process for obtaining partial sensitivity.

[0232] In the computational framework of the elastic component, the sensitivity parameter originates from two independent statistical inference stages. The first estimated variance reflects the uncertainty of empirical sensitivity in watershed-by-watershed time-series regression, calculated by extracting the residual variance of the snow cover sequence regression from the mediation mechanism index and combining it with the sum of squared deviations of the sample independent variables. The second estimated variance is obtained by extracting the coefficient covariance matrix generated by constrained least squares fitting, reflecting the uncertainty of partial sensitivity in global spatial regression. Specifically, taking the partial sensitivity of the energy mechanism as an example, its second estimated variance is composed of the variance of the main effect coefficient, the variance of the climate interaction coefficient multiplied by the square of the climate state variable, and the covariance of both.

[0233] Step 802: Based on the first-order Taylor expansion principle, the first estimated variance and the second estimated variance are calculated independently along the corresponding empirical sensitivity and partial sensitivity in a chain error propagation to obtain the target approximate variance of each explicit mechanism elastic component.

[0234] Because the estimation models for empirical sensitivity and partial sensitivity differ significantly in structure and data organization, they statistically approximate the assumption of independent distribution. Utilizing this independence condition, a first-order Taylor expansion theory is applied to propagate statistical errors for the product formula structure of the explicit mechanism elastic components. First, a common scaling factor determined by the watershed climate mean is defined; the coefficient of variation of this factor is negligible compared to the sensitivity parameter. The variances of the two stages are multiplied by the squared terms of the corresponding sensitivity estimates for the other stage, summed, and then multiplied by the squared term of the common scaling factor to obtain the final variance propagation result.

[0235] Var Ε_k_i =A i2 *(c k_i 2 *Var S_k_i +S k_i 2 *Var c_k_i );

[0236] Among them, Var Ε_k_i To represent the target approximate variance corresponding to the explicit mechanism elastic component, Ai is the basic environment scaling factor, and c is... k_i For empirical sensitivity, Var S_k_i For the second estimate of variance, S k_i For biased sensitivity, Var c_k_i This is the first estimated variance.

[0237] The above chain-like error propagation formula shows that the uncertainty of the explicit mechanism elastic component is contributed by the variances of two stages: the variance Var of the first stage (empirical sensitivity estimation). c_k_i Squared amplification of the partial sensitivity estimate; Var of the variance in the second stage (partial sensitivity estimate) S_k_i The variance of the elastic component is amplified by the square of the empirical sensitivity estimate; the sum of the two contributions is then scaled by the square of the basic environmental scaling factor to obtain the final target approximate variance. Engineers can substitute the regression fitting results for specific watersheds into the above formula for numerical calculation.

[0238] Step 803: Construct corresponding confidence intervals based on the target approximate variance of each explicit mechanism elastic component to assess the statistical significance of the quantitative contribution rate of each physical transmission path.

[0239] Specifically, confidence intervals are constructed based on the target approximate variance of each explicit mechanism elastic component, and these confidence intervals are incorporated into the hydrological elasticity attribution results to assess the statistical significance of the quantitative contribution rate of each physical transport path.

[0240] After obtaining the approximate variance of the target, the standard deviation parameter is extracted by performing a square root operation. Combined with the set normal distribution quantiles, an interval range for hypothesis testing is constructed. Preferably, a significance level of 0.05 is set, and the corresponding quantile of 1.96 is extracted to construct a 95% confidence interval. The corresponding explicit mechanism elasticity component point estimates are subtracted from and added to by 1.96 times the standard deviation, respectively, to form the lower and upper limits of the interval.

[0241] Furthermore, it is determined whether the constructed confidence interval contains a value of 0. If the confidence interval does not contain a value of 0, then the contribution of the corresponding physical transfer path in the target watershed is determined to be statistically significant.

[0242] This uncertainty quantification scheme utilizes the chain-based architecture's unique multiplication expression structure to solve the problem of the lack of confidence measurement in traditional hydrological elasticity analysis, providing a reliable model diagnostic boundary for watersheds with low data signal-to-noise ratios.

[0243] According to one aspect of this application, internal consistency verification primarily involves a comparison and judgment of three physical laws. For the first type of target watershed where a decrease in underlying snow cover leads to a decrease in annual runoff, the absolute value of the elastic component of the energy mechanism is calculated and determined whether it is simultaneously greater than the absolute values ​​of the elastic components of the intensity mechanism and the time-series mechanism, to verify the dominant position of the energy mechanism in this type of watershed. For the second type of target watershed where a decrease in underlying snow cover leads to an increase in annual runoff, the sum of the absolute values ​​of the elastic components of the intensity mechanism and the time-series mechanism is calculated and determined whether it is greater than the absolute value of the elastic component of the energy mechanism. A threshold test is performed on the residual contribution rate to verify the structural sufficiency of the three-mechanism decomposition.

[0244] When examining the residual contribution rate, residual contribution rate data for all target watersheds are extracted. The established examination criterion is to count the proportion of watersheds with an absolute residual contribution rate less than a predetermined threshold to the total number of watersheds, and determine whether this proportion is greater than the preset target pass rate. Preferably, the predetermined threshold is set to 30%, and the target pass rate is set to 0.7.

[0245] P res =1 / N w *∑(1(|η res_i |<0.3));

[0246] Among them, P res To satisfy the residual threshold condition, the proportion of watersheds, N w Let ∑() be the total number of watersheds participating in the verification, ∑() be a function that performs statistical summation on watersheds that meet the predetermined conditions, 1() be an indicator function that takes a value of 1 when the internal conditions are met, and a value of 0 otherwise, and η be the total number of watersheds participating in the verification. res_i Let || be the residual contribution rate of the i-th watershed, and || denotes the operation of extracting the absolute value.

[0247] When calculating the residual contribution rate, the absolute value of the residual elastic component of the target watershed is divided by the sum of the absolute values ​​of the elastic components of each explicit mechanism and the absolute values ​​of the residual elastic components to obtain the corresponding residual contribution rate:

[0248] η res_i =|ε residual_i | / (∑ k (|ε fs_k_i |)+|ε residual_i |);

[0249] Where, η res_i ε represents the residual contribution rate of the i-th watershed.residual_i For the residual elastic component, ε fs_k_i For each explicit mechanism elastic component.

[0250] For external prediction validation, the hold-out method is used to verify the generalization ability of the nonlinear watershed parameter model. All acquired target watersheds are randomly divided into a training set containing a predetermined proportion of samples and a test set containing the remaining samples. Preferably, the training set contains 70% of the samples, and the test set contains 30% of the samples. Offline regression fitting of the aforementioned nonlinear watershed parameter model is performed on the training set. On the test set, the climate state variables and sensitivity estimation results of the test set watersheds are extracted, and the mechanism elasticity components and total elasticity prediction values ​​of each test watershed are predicted using the model coefficient values ​​fixed in the training set.

[0251] To quantify prediction accuracy, the sign agreement rate and root mean square error are extracted as quantitative evaluation indicators. The sign agreement rate is calculated based on the predicted elasticity and the reference elasticity on the test set.

[0252] SAR=1 / N test *∑(1(sgn(ε fs_i )==sgn(ε fs_obs_i )));

[0253] Wherein, SAR is the sign consistency rate, which characterizes the accuracy of elastic direction prediction, and its value is limited to 0~1, N test To test the total number of watersheds, ∑() is the summation function, 1() is the indicator function, and sgn() is the function to extract the sign of the numerical value, returning a positive sign when the input value is greater than 0 and a negative sign when the input value is less than 0. ε fs_i ε represents the elasticity of the model's predicted output. fs_obs_i The reference elasticity is calculated independently by the verification mechanism.

[0254] Furthermore, the root mean square error between the predicted elasticity on the test set and the reference elasticity is calculated to evaluate the accuracy of the model's prediction of the elasticity magnitude.

[0255] RMSE = sqrt(1 / N) test *∑((ε fs_i -ε fs_obs_i ) 2 ));

[0256] Where RMSE is the root mean square error, sqrt() is the square root operation, and N test Let ε be the total number of watersheds in the test set, ∑() be the summation function executed on the test set, and ε be the total number of watersheds in the test set. fs_i ε represents the elasticity of the model's predicted output. fs_obs_i For reference flexibility.

[0257] To obtain benchmark data independent of the three-mechanism decoupling framework when calculating the reference elasticity, the total sensitivity slope of the target watershed over the full time window is acquired. The reference elasticity is then obtained by multiplying the total sensitivity slope by the baseline environment scaling factor.

[0258] ε fs_obs_i =A i *γ s_i ;

[0259] Where, ε fs_obs_i To verify the reference elasticity of the required input, A i Based on the environment scaling factor, γ s_i The slope represents the overall sensitivity.

[0260] The calculation and judgment of the above model verification and evaluation mechanism completed the closed-loop verification of the reliability of the hydrological elastic decomposition scheme throughout the entire process.

[0261] Taking a mid-latitude seasonal snow-covered watershed as an example, the multi-year average precipitation P is approximately 600 mm, the multi-year average potential evapotranspiration EP0 is approximately 720 mm, the drought index Φ is approximately 1.2, the multi-year average snow cover percentage fs is approximately 0.25, and the multi-year average annual runoff Q is approximately 180 mm. The meteorological and hydrological sequence of this watershed spans 30 years.

[0262] Along the energy transfer pathway, for every 0.1 unit decrease in snow cover percentage, net radiation (IE) during the snow season increases. The empirical sensitivity cE is negative, representing the inverse response of net radiation decrease to increased snow cover. This inhibits runoff by increasing evapotranspiration and consuming more water, with the elastic component of the corresponding energy mechanism pointing in the same direction as the runoff change. Along the water input intensity pathway, reduced snow cover shortens the concentrated snowmelt period, increasing the daily average liquid water intensity (ILWI). The empirical sensitivity cLWI is negative, representing the inverse response of decreased liquid water input intensity to increased snow cover. This leads to a decrease in watershed water storage efficiency and a partial sensitivity S. LWI A negative value indicates a promoting effect on runoff. Considering the elastic components of the three pathways, the energy mechanism contributes the most to the total elasticity of the watershed, consistent with its characteristic of being a Class I watershed where reduced snow cover leads to reduced runoff.

[0263] In the example watershed described above, the total elasticity value obtained through chain mechanism decomposition is consistent with the reference elasticity value independently calculated using the direct regression method of watershed water storage parameter sequence on snow cover sequence, verifying the mathematical completeness of the decomposition framework. The sum of the elasticity components of the three explicit physical paths accounts for the main part of the total elasticity, and the absolute contribution rate of the residual elasticity component is less than 30% of the total elasticity, indicating that the three extracted mechanism paths have high structural sufficiency.

[0264] The results of the above embodiments demonstrate that the method of the present invention can achieve multi-channel quantitative decoupling of complex hydrological responses. It is understood that the specific contribution ratios of each pathway may vary depending on the geographical location, climate type, and snow cover characteristics of the watershed.

[0265] In summary, the watershed hydrological elastic attribution method based on chain mechanism decomposition of the present invention includes: acquiring meteorological and hydrological sequences, underlying snow cover sequences, and climate state variables of the target watershed; extracting multi-dimensional mediation mechanism indicators characterizing multiple independent physical transfer paths based on the meteorological and hydrological sequences and underlying snow cover sequences, wherein the physical transfer paths include at least energy transfer paths, water input intensity paths, and hydrothermal temporal correlation paths; constructing time-series observation samples based on the multi-dimensional mediation mechanism indicators and underlying snow cover sequences, and estimating the empirical sensitivity of each multi-dimensional mediation mechanism indicator to snow cover changes; and inputting the multi-dimensional mediation mechanism indicators and climate state variables into... The data are fed into a pre-constructed nonlinear watershed parameter model to obtain the partial sensitivity of watershed water storage parameters to various multi-dimensional mediating mechanism indicators. Chain decomposition and reconstruction are performed by combining empirical sensitivity, partial sensitivity, and a basic environmental scaling factor determined based on meteorological and hydrological sequences to calculate the explicit mechanism elastic components corresponding to each of the aforementioned physical transmission paths. The implicit underlying surface influence path features not captured by the explicit mechanism elastic components are extracted, and the residual elastic components of the corresponding target watershed are calculated. The explicit mechanism elastic components and residual elastic components are aggregated to output the total elasticity of the target watershed runoff to snow cover change and the quantitative contribution rate of each physical transmission path.

[0266] This invention introduces a multidimensional mediation mechanism index and combines it with the chain rule of calculus to elastically decouple the macroscopic black box of runoff into independent physical conduction contributions of energy, moisture intensity, and hydrothermal timing. This achieves a leap from superficial statistics to quantitative decoupling of physical mechanisms. It solves the problem that existing macroscopic statistical models struggle to distinguish the impact of underlying surface changes on the specific physical conduction of runoff.

[0267] By employing integral deviation distance based on the normalized cumulative distribution function for quantification, the ability to capture temporal misalignment is broadened from a single extreme point to a full-time distribution pattern, reducing the asynchronous assessment error in multimodal hydrological environments. This solves the problem that traditional single-point time difference methods are insufficient to accurately characterize the hydrothermal misalignment features under complex multi-peak precipitation patterns.

[0268] A nonlinear model incorporating climate interaction characteristics was constructed and directional inequality constraints were applied, enabling sensitivity inference to adapt to the arid and wet climate background of the watershed and ensuring the objectivity and rationality of sensitivity reasoning.

[0269] The residuals are extracted using the closed-loop difference method, which eliminates the systematic residual stripping omissions caused by traditional statistical fitting and elevates the summation of biased components to a mathematically complete closed state.

[0270] Based on the first-order Taylor expansion principle, multi-stage chain error propagation calculation is performed, which gives the attribution results a dynamic confidence interval boundary, transforming the single-value estimation lacking statistical support into a robust evaluation system with strict significance testing.

[0271] The preferred embodiments of the present invention have been described in detail above. However, the present invention is not limited to the specific details in the above embodiments. Within the scope of the technical concept of the present invention, various equivalent transformations can be made to the technical solutions of the present invention, and these equivalent transformations all fall within the protection scope of the present invention.

Claims

1. A watershed hydrological elastic attribution method based on chain mechanism decomposition, characterized in that, include: Obtain meteorological and hydrological sequences, underlying snow cover sequences, and climate state variables for the target watershed; Based on meteorological and hydrological sequences and underlying snow cover sequences, mediating mechanism indicators characterizing physical transport pathways are extracted. Determine the empirical sensitivity of each mediation mechanism indicator to changes in underlying snow cover sequence; By inputting the mediating mechanism indicators and climate state variables into a pre-built watershed parameter model, the partial sensitivity of watershed water storage parameters to each mediating mechanism indicator can be obtained. By combining empirical sensitivity, partial sensitivity, and basic environmental scaling factor determined based on meteorological and hydrological sequences, the explicit mechanism elastic components corresponding to each physical transport path are calculated. The residual elastic components of the target watershed are calculated based on the explicit mechanism elastic components. By aggregating the explicit mechanism elastic components and residual elastic components, the hydrological elastic attribution results of the target watershed for snow cover changes are output.

2. The method according to claim 1, characterized in that, The physical transfer pathways include energy transfer pathways, water input intensity pathways, and hydrothermal temporal correlation pathways; correspondingly, the mediating mechanism indicators include energy mechanism indicators reflecting the characteristics of net radiation during the snow season, liquid water input intensity indicators reflecting the characteristics of the average daily input of liquid water in the year, and hydrothermal asynchrony indicators reflecting the characteristics of the time misalignment between energy supply and water supply.

3. The method according to claim 2, characterized in that, Determine the hydrothermal asynchrony index, including: The liquid water input sequence and the net surface radiation sequence of the target watershed were obtained based on meteorological and hydrological sequences. The first normalized cumulative distribution function of the liquid water input sequence on the time axis and the second normalized cumulative distribution function of the net surface radiation sequence on the time axis were constructed respectively. The deviation distribution values ​​of the first normalized cumulative distribution function and the second normalized cumulative distribution function at each corresponding time node are extracted. The sum of the deviation distribution values ​​over the entire time period is calculated as the distribution integral distance, and the distribution integral distance is used as the hydrothermal asynchrony index.

4. The method according to claim 1, characterized in that, Determine the empirical sensitivity of each mediation mechanism indicator to changes in underlying snow cover sequence, including: According to the preset window length and sliding step size, the underlying snow accumulation sequence and various mediation mechanism indicators are divided into time series, and multiple overlapping time windows are extracted. The mean data within each overlapping time window is calculated separately. By differentiating the mean data of adjacent overlapping time windows, the sliding difference sequence of each mediation mechanism indicator and the snow accumulation sliding difference sequence of the underlying snow accumulation sequence are obtained. The sliding difference sequences of each mediation mechanism indicator were used to perform regression fitting on the snow cover sliding difference sequence to determine the empirical sensitivity of each mediation mechanism indicator.

5. The method according to claim 1, characterized in that, The watershed parameter model includes the univariate main effect characteristics of each mediating mechanism indicator, as well as the parameter interaction characteristics of the combination of each mediating mechanism indicator and climate state variables. The partial sensitivity of watershed water storage parameters to various mediating mechanism indicators was obtained, including: The univariate main effect characteristics of each mediating mechanism indicator and the sensitive adjustment magnitude of each parameter interaction characteristic under the current climate state variable are extracted using a pre-constructed watershed parameter model. By linearly superimposing the sensitivity moderating magnitudes of the main effects of each univariate with the corresponding interaction features of each parameter, the partial sensitivity of each mediation mechanism indicator is obtained.

6. The method according to claim 1, characterized in that, Calculate the explicit mechanism elastic components for each physical transport path, including: For any target physical transfer path among all physical transfer paths, the empirical sensitivity and the corresponding partial sensitivity of the mediation mechanism index corresponding to the target physical transfer path in the target watershed are extracted respectively. Multiplying the empirical sensitivity and partial sensitivity corresponding to the target physical transfer path by the scaling factor of the basic environment yields the explicit mechanism elastic component of the target physical transfer path.

7. The method according to claim 1, characterized in that, The residual elastic components of the target watershed are calculated based on the explicit mechanism elastic components, including: Obtain the model fitting prediction values ​​of the same watershed parameter model as the explicit mechanism elastic component, and extract the actual observed values ​​of the parameters of the target watershed. Calculate the numerical difference between the actual observed values ​​of each parameter and the corresponding model-fitted predicted values, and generate a serialized sequence of fitting residuals. Determine the residual sensitivity coefficient of the fitted residual sequence to changes in the underlying snow cover sequence; Multiplying the residual sensitivity coefficient by the basic environment scaling factor yields the residual elastic component of the corresponding target watershed.

8. The method according to claim 1, characterized in that, When calculating the residual elastic components of the target watershed, the following are included: Based on meteorological and hydrological sequences, the water storage parameter sequences of the corresponding target watershed are obtained, and the overall sensitivity of the water storage parameter sequences to changes in the underlying snow cover sequence is determined. Extract the product of the empirical sensitivity and partial sensitivity corresponding to each physical transmission path as the explicit mechanism attribution contribution of each physical transmission path. The sum of the explicit mechanism attribution contributions of each physical transmission path is subtracted from the overall total sensitivity, and the closure difference residual for the corresponding target watershed is extracted. By combining the closure difference residual with the basic environmental scaling factor, the residual elastic component of the corresponding target watershed is obtained.

9. The method according to claim 1, characterized in that, The watershed parameter model is pre-built through the following offline steps: Historical sample sets from multiple watersheds were obtained, including pre-extracted samples of watershed water storage parameter changes, as well as corresponding samples of mediating mechanism indicators and climate state variables. The historical sample set is used to screen features based on a regularized regression algorithm to determine a subset of significant predictive variables; Under the pre-defined physical constraints, least squares regression fitting is performed using a subset of significant predictor variables to determine the model coefficients of the watershed parameter model.

10. The method according to claim 1, characterized in that, Determine the base environment scaling factor, including: Long-term average precipitation and potential evapotranspiration are extracted from meteorological and hydrological sequences, and drought characteristic constants of the target watershed are calculated by combining them with the water-heat balance equation. For the fundamental equation of watershed water-heat balance that characterizes the balance between precipitation, actual evapotranspiration and potential evapotranspiration, the partial derivative function of the target actual evapotranspiration with respect to the watershed water storage parameter is solved. Substituting the long-term average precipitation and drought characteristic constant into the partial derivative function, the absolute change in the actual evapotranspiration of the target area is calculated when the watershed water storage parameters change by a unit. Combining the long-term average snow cover ratio and long-term average annual runoff of the target watershed, the absolute change is converted to a runoff elastic scale to obtain the basic environmental scaling factor.

Images