Hydrological model parameter correction method, device and system

Abstract

The application provides a hydrological model parameter correction method, device and system, the method comprises the following steps: obtaining the to-be-corrected parameters of a hydrological model, and constructing an initial background field parameter set; taking the initial background field parameter set as a current parameter set, and performing a parameter correction step: mapping the current parameter set to an observation space to obtain a simulated flow set; calculating a parameter perturbation matrix and a flow perturbation matrix; based on the measured flow data of the hydrological station in the assimilation window, the flow perturbation matrix, and taking a set weight vector as an independent variable, constructing an objective function; iteratively solving the objective function, and combining the parameter perturbation matrix to determine a corrected posterior parameter analysis field; based on the flow perturbation matrix and the posterior parameter analysis field, updating the current parameter set, and continuing to perform the parameter correction step to realize a parameter dynamic correction process. The application can use the measured flow data to accurately correct the model parameters in a cycle, which is beneficial to improving the accuracy of hydrological simulation and flood forecasting.

Description

Technical Field

[0001] This application relates to the field of hydrology and water resources management technology, and in particular to a method, apparatus and system for correcting hydrological model parameters. Background Technology

[0002] Existing Xin'anjiang models typically calibrate a fixed set of parameters based on historical data for watershed hydrological simulation and flood forecasting. However, actual watershed underlying surface conditions (such as soil moisture content distribution and vegetation interception capacity) exhibit significant spatiotemporal heterogeneity with seasonal changes and rainfall events. This simulation method based on fixed parameters not only struggles to adapt to the dynamic changes in watershed characteristics but also fails to effectively compensate for systematic errors in rainfall input and biases introduced by model structural generalization, resulting in a difficulty in achieving consistent simulation accuracy for different flood events.

[0003] To address the aforementioned issues, ensemble Kalman filtering (EnKF) and its variants have been introduced in this field. Although these methods can provide flow-dependent background error covariance estimates through ensemble forecasts, in practical applications, they are limited by computational resources, and the limited ensemble sample size easily introduces sampling errors. Furthermore, EnKF is essentially a sequential assimilation algorithm, which only uses observation information from a single moment for updates, resulting in temporal discontinuities in the analysis field within the assimilation window, making it difficult to guarantee the dynamic consistency of hydrological processes. Summary of the Invention

[0004] The purpose of this application is to provide a method, device and system for calibrating hydrological model parameters, which realizes a hydrological model parameter calibration method based on four-dimensional ensemble variation, and can use measured flow data to cyclically and accurately calibrate model parameters, which is beneficial to improving the accuracy of hydrological simulation and flood forecasting.

[0005] In a first aspect, this application provides a method for correcting hydrological model parameters. The method includes: obtaining the parameters to be corrected from the hydrological model and constructing an initial background field parameter set; using the initial background field parameter set as the current parameter set and performing parameter correction steps: mapping the current parameter set to the observation space to obtain a simulated flow set; calculating the parameter perturbation matrix of the current parameter set and the flow perturbation matrix of the simulated flow set; based on the measured flow data and flow perturbation matrix of hydrological stations within the assimilation window, constructing an objective function containing background prior constraints and flow observation fitting with the set weight vector as the independent variable; iteratively solving the objective function and determining the corrected posterior parameter analysis field by combining the parameter perturbation matrix; updating the current parameter set based on the flow perturbation matrix and the posterior parameter analysis field, and continuing to perform parameter correction steps to realize the dynamic parameter correction process.

[0006] Furthermore, the steps of obtaining the parameters to be calibrated from the hydrological model and constructing the initial background field parameter set include: extracting key physical parameters from the Xin'anjiang model as parameters to be calibrated; generating multiple sets of parameter vectors using the Monte Carlo sampling method based on the prior distribution information of the parameters to be calibrated; and mapping the multiple sets of parameter vectors to the computational space through a transformation function to obtain the initial background field parameter set that satisfies the Gaussian distribution assumption.

[0007] Furthermore, the above-mentioned step of mapping the current parameter set to the observation space to obtain the simulated flow set includes: inputting each set of parameter vectors in the current parameter set into the observation operator determined by the Xin'anjiang model, so as to map the parameter vectors from the computation space to the observation space and obtain the simulated flow set used to characterize the simulated flow process within the assimilation window.

[0008] Furthermore, the steps of calculating the parameter perturbation matrix of the current parameter set and the flow perturbation matrix of the simulated flow set include: calculating the mean background field parameters of the current parameter set and the mean flow of the simulated flow set; subtracting the mean background field parameters from the background field parameter vector corresponding to each set member in the current parameter set and multiplying by a normalization coefficient to obtain the parameter perturbation matrix; and subtracting the mean flow from the simulated flow vector corresponding to each set member in the simulated flow set and multiplying by a normalization coefficient to obtain the flow perturbation matrix.

[0009] Furthermore, the steps described above, based on the measured flow data and flow disturbance matrix of hydrological stations within the assimilation window, and using the set weight vector as the independent variable, to construct an objective function comprising two parts: background prior constraints and flow observation fitting, include: constructing an observation vector using the measured flow data of hydrological stations within the assimilation window; calculating the difference vector between the observation vector and the mean of the simulated flow sample within the assimilation window; constructing an observation error covariance matrix based on the error characteristics of the observation data within the assimilation window; and constructing an objective function comprising two parts: background prior constraints and flow observation fitting, based on the difference vector, the observation error covariance matrix, the flow disturbance matrix, and the regularization coefficient, using the set weight vector as the independent variable.

[0010] Furthermore, the steps described above, which construct an objective function comprising background prior constraints and flow observation fitting based on the difference vector, observation error covariance matrix, flow disturbance matrix, and regularization coefficients, with the set weight vector as the independent variable, include: The objective function J(w) is determined according to the following formula: J(w)=0.5γw T w+0.5(d Yb'w) T R 1 (d Yb'w); Where w represents the set weight vector; γ represents the regularization coefficient; d represents the difference vector; Yb' represents the flow disturbance matrix; R represents the observation error covariance matrix; and T represents the matrix transpose operation.

[0011] Furthermore, the steps described above for iteratively solving the objective function and determining the corrected posterior parameter analysis field by combining the parameter perturbation matrix include: iteratively solving the constructed objective function using a sequential least squares programming algorithm to obtain the optimal weight vector that minimizes the objective function; calculating the optimal correction amount of the parameters based on the parameter perturbation matrix and the optimal weight vector; and superimposing the optimal correction amount onto the mean of the background field parameters to obtain the corrected posterior parameter analysis field.

[0012] Furthermore, the steps described above for updating the current parameter set based on the flow perturbation matrix and the posterior parameter analysis field include: calculating the Hessian matrix of the objective function based on the flow perturbation matrix, the observation error covariance matrix, and the regularization coefficient; performing Cholesky decomposition on the Hessian matrix and constructing a transformation matrix for updating the perturbation weights based on the decomposition results; performing a linear transformation on the parameter perturbation matrix based on the transformation matrix to obtain the posterior perturbation matrix; superimposing the posterior perturbation matrix onto the posterior parameter analysis field to obtain the updated posterior parameter set; performing an inverse mapping on the analysis field members in the posterior parameter set using the inverse function of the transformation function to obtain a physically constrained posterior parameter set; and updating the current parameter set using the physically constrained posterior parameter set.

[0013] Secondly, this application also provides a hydrological model parameter calibration device, comprising: a set construction module for acquiring the parameters to be calibrated of the hydrological model and constructing an initial background field parameter set; a parameter calibration module for using the initial background field parameter set as the current parameter set and performing parameter calibration steps: mapping the current parameter set to the observation space to obtain a simulated flow set; calculating the parameter perturbation matrix of the current parameter set and the flow perturbation matrix of the simulated flow set; constructing an objective function containing background prior constraints and flow observation fitting, based on the measured flow data and flow perturbation matrix of hydrological stations within the assimilation window, with the set weight vector as the independent variable; iteratively solving the objective function and determining the posterior parameter analysis field after calibration, combined with the parameter perturbation matrix; and a set update module for updating the current parameter set based on the flow perturbation matrix and the posterior parameter analysis field, continuing to execute the parameter calibration steps, and realizing the dynamic parameter calibration process.

[0014] Thirdly, this application also provides a hydrological model parameter correction system, including a processor and a memory, wherein the memory stores computer-executable instructions that can be executed by the processor, and the processor executes the computer-executable instructions to implement the method of the first aspect above.

[0015] The hydrological model parameter calibration method, apparatus, and system provided in this application first obtain the parameters to be calibrated for the hydrological model and construct an initial background field parameter set. Then, using the initial background field parameter set as the current parameter set, the parameter calibration steps are performed: mapping the current parameter set to the observation space to obtain a simulated flow set; calculating the parameter perturbation matrix of the current parameter set and the flow perturbation matrix of the simulated flow set; based on the measured flow data and flow perturbation matrix of hydrological stations within the assimilation window, constructing an objective function containing background prior constraints and flow observation fitting with the set weight vector as the independent variable; iteratively solving the objective function and determining the posterior parameter analysis field after calibration by combining the parameter perturbation matrix; updating the current parameter set based on the flow perturbation matrix and the posterior parameter analysis field, and continuing to execute the parameter calibration steps to achieve a dynamic parameter calibration process. This application implements a hydrological model parameter calibration method based on four-dimensional ensemble variation, which can use measured flow data to cyclically and accurately calibrate model parameters, thus improving the accuracy of hydrological simulation and flood forecasting. Attached Figure Description

[0016] To more clearly illustrate the technical solutions in the specific embodiments of this application or the prior art, the drawings used in the description of the specific embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of this application. For those skilled in the art, other drawings can be obtained from these drawings without creative effort.

[0017] Figure 1 A flowchart of a hydrological model parameter correction method provided in this application embodiment; Figure 2 A schematic diagram showing the comparison between simulated flow rate and measured flow rate before and after calibration, provided for an embodiment of this application; Figure 3 A structural block diagram of a hydrological model parameter correction device provided in an embodiment of this application; Figure 4 This is a schematic diagram of the structure of an electronic device provided in an embodiment of this application. Detailed Implementation

[0018] The technical solutions of this application will be clearly and completely described below with reference to the embodiments. Obviously, the described embodiments are only some embodiments of this application, not all embodiments. Based on the embodiments of this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0019] While the four-dimensional variational (4DVar) method can theoretically solve the problem of time discontinuity through global constraints, its engineering implementation is extremely difficult. This method relies heavily on the development of adjoint models, and for hydrological models with complex logic and numerous nonlinear physical processes, the derivation and coding of adjoint models are extremely cumbersome, resulting in high development and maintenance costs, as well as enormous computational demands, making it difficult to meet the stringent timeliness requirements of real-time flood forecasting.

[0020] Based on this, embodiments of this application provide a hydrological model parameter calibration method, apparatus, and system, which realizes a hydrological model parameter calibration method based on four-dimensional ensemble variation. It can use measured flow data to cyclically and accurately calibrate model parameters, which is beneficial to improving the accuracy of hydrological simulation and flood forecasting.

[0021] To facilitate understanding of this embodiment, a hydrological model parameter correction method disclosed in this application embodiment will first be described in detail. Figure 1 A flowchart of a hydrological model parameter correction method provided in this application embodiment is shown. The method specifically includes the following steps: Step S1: Obtain the parameters to be corrected for the hydrological model and construct the initial background field parameter set; Physical parameters in the Xin'anjiang model that are sensitive to the runoff generation and confluence process were selected as parameters to be corrected.

[0022] Based on the prior information of each parameter, an initial set of background field parameters is generated using Monte Carlo sampling, and the sample mean of this set is defined as the background field of the parameter to be corrected.

[0023] Step S2: Using the initial background field parameter set as the current parameter set, perform the parameter correction step: Step S21: Map the current parameter set to the observation space to obtain the simulated flow set; Each member variable in the initial parameter set is input into the hydrological model to obtain the background simulated flow sequence within the assimilation window, thus obtaining the simulated flow set.

[0024] Step S22: Calculate the parameter perturbation matrix of the current parameter set and the flow perturbation matrix of the simulated flow set; The current parameter set and the simulated flow set are centered to determine the mean of the set samples. Then, the parameter perturbation matrix and the flow perturbation matrix are constructed using the set sample mean and the normalization coefficient.

[0025] Step S23: Based on the measured flow data and flow disturbance matrix of the hydrological stations within the assimilation window, construct an objective function that includes background prior constraints and flow observation fitting, using the set weight vector as the independent variable. When constructing the objective function, a regularization coefficient is introduced to adjust the weight ratio of the background term and the observation term, so as to approximate the background error covariance.

[0026] Step S24: Iteratively solve the objective function and determine the corrected posterior parameter analysis field by combining the parameter perturbation matrix; In this embodiment, the Sequential Least Squares Programming (SLSQP) algorithm can be used to solve the objective function, obtain the optimal weight vector, and then calculate the corrected posterior parameter analysis field.

[0027] Step S3: Based on the flow disturbance matrix and the posterior parameter analysis field, update the current parameter set and continue to execute the parameter correction step to realize the dynamic parameter correction process.

[0028] By decomposing the Hessian matrix of the objective function, a set of physically consistent posterior parameters is generated and used as the initial background field for the next time period to achieve rolling forecasting.

[0029] Step S1 above involves obtaining the parameters to be corrected for the hydrological model and constructing an initial set of background field parameters, specifically including: (1) Extract key physical parameters from the Xin'anjiang model as parameters to be corrected; Physical parameters sensitive to runoff generation and confluence processes in the Xin'anjiang model were selected as parameters to be corrected. These parameters include, but are not limited to: tensional water storage capacity (WM), the power of the soil water storage capacity curve (B), saturated hydraulic conductivity (Ks), free water storage capacity (SM), source delineation coefficients (KI, KG), and runoff coefficients (CS, CI, CG). Upper and lower limits and prior distributions conforming to hydrophysical laws were set for each parameter, and a preset duration before the assimilation window was set as the model warm-up period.

[0030] (2) Based on the prior distribution information of the parameters to be corrected, multiple sets of parameter vectors are generated using the Monte Carlo sampling method; Prior distribution information refers to an estimate of the parameter to be corrected based on existing knowledge and experience before parameter correction. For example, based on experience, it is believed that a certain parameter follows a normal distribution with a mean of miu and a variance of sigma; this normal distribution is the prior distribution. Based on the prior distribution information of the parameters, a Monte Carlo sampling method is used to generate N sets of parameter vectors.

[0031] (3) By transforming multiple sets of parameter vectors to the computation space, the initial background field parameter set that satisfies the Gaussian distribution assumption is obtained.

[0032] To address the physical boundary constraints of hydrological model parameters, a nonlinear bijective transformation function is constructed between the physical and computational spaces. This transforms the parameters to be corrected from the bounded parameter region. Unbounded open intervals of the entire real number field mapped to the computational space .

[0033] Forward transformation: ; In the formula, and These are the lower and upper limits of the physical threshold for the corresponding parameters, respectively.

[0034] The generated parameter vectors are mapped to the computational space through a transformation function, thereby forming an initial background field parameter set matrix that lies in the computational space and satisfies the Gaussian distribution assumption. .

[0035] Step S21 above maps the current parameter set to the observation space to obtain the simulated flow set, specifically including: Each parameter vector in the current parameter set is input into the observation operator determined by the Xin'anjiang model to map the parameter vector from the computation space to the observation space, thereby obtaining a set of simulated flow rates used to characterize the simulated flow process within the assimilation window.

[0036] In practice, the observation operator is used. The operator performs a mapping, which involves a mapping process from the computation space to the observation space, i.e. Here, M represents the runoff generation and concentration calculation process of the Xin'anjiang model. Each set of parameter vectors from the generated background field parameter set is input into the Xin'anjiang model. The model is run up to the current time step, and the simulated flow process within the assimilation window is output, thus obtaining the simulated flow set. ,Right now .

[0037] Further, step S22 above, calculating the parameter perturbation matrix of the current parameter set and the flow perturbation matrix of the simulated flow set, specifically includes: (1) Calculate the mean value of the background field parameters of the current parameter set and the mean value of the flow rate of the simulated flow rate set; Calculate the mean value of background field parameters in the background field parameter set. : ; Average flow rate of the simulated flow set : .

[0038] (2) Subtract the mean value of the background field parameters from the background field parameter vector corresponding to each set member in the current parameter set, and multiply by the normalization coefficient to obtain the parameter perturbation matrix; Construct the parameter perturbation matrix for the unbiased estimate: ; Background covariance matrix .

[0039] (3) Subtract the mean flow from the simulated flow vector corresponding to each member of the simulated flow set, and multiply by the normalization coefficient to obtain the flow disturbance matrix.

[0040] For a given background field state and tiny perturbations Using the first-order approximation principle of Taylor expansion: ; Ignore higher-order terms Rearranging the above Taylor expansion, we get: ; Each set member in the set Both can be seen as the mean A disturbance caused Then the above approximate relationship becomes: ; In ensemble methods, to construct the statistical background field error covariance, a statistical perturbation relative to the ensemble center must be defined. Theoretically, the Taylor expansion is a simulated value based on the background field mean. A linear approximation is performed. However, in the ensemble statistical calculation of this embodiment, the ensemble mean of the simulated flow is directly used as the centralization benchmark to ensure that the flow disturbance matrix satisfies the statistical property of zero mean. The simulated flow vector corresponding to each ensemble member is subtracted from the flow mean and multiplied by the normalization coefficient. ; The error introduced by this replacement is negligible when the model has low nonlinearity (i.e., small ensemble perturbation). In practice, the effectiveness of this linearization approximation is ensured by controlling the magnitude of the ensemble perturbation (e.g., using 15% relative error in this embodiment).

[0041] Furthermore, in step S23 above, based on the measured flow data and flow disturbance matrix of the hydrological stations within the assimilation window, an objective function is constructed using the set weight vector as the independent variable, comprising two parts: background prior constraints and flow observation fitting. Specifically, this includes: (1) Construct observation vectors using measured flow data from hydrological stations within the assimilation window; Collect measured flow data from hydrological stations within an assimilation window (e.g., the past 24 hours). Assume the assimilation window contains... At any moment (from) The measured flow values ​​at each time point are arranged in chronological order to obtain the measured flow data within the assimilation window, forming an observation vector. =[ ].

[0042] (2) Calculate the difference vector between the observed vector and the mean of the simulated flow sample within the assimilation window period; Calculate the measured flow vector within the assimilation window period. and the mean of simulated flow samples The difference vector between .

[0043] (3) Construct the observation error covariance matrix based on the error characteristics of the observation data within the assimilation window; Based on the number of observations within the assimilation window (i.e., the length of the observation vector) Construct a dimension of The square matrix is ​​used as the observation error covariance matrix. Based on the assumption that the observation errors at different times are independent, matrix R is set as a diagonal matrix, meaning that all off-diagonal elements are 0. The elements on the diagonal of the matrix correspond to the error variance of the flow measurement at each time point (e.g., set as a fixed empirical value or a dynamic value that varies with the flow rate). Then, the diagonal matrix is ​​directly inverted to obtain... This is used for the subsequent weighted calculation of the cost function.

[0044] (4) Based on the difference vector, observation error covariance matrix, flow disturbance matrix, and regularization coefficient, and with the set weight vector as the independent variable, construct an objective function that includes background prior constraints and flow observation fitting.

[0045] The objective function J(w) is determined according to the following formula: J(w)=0.5γw T w+0.5(d Yb'w) T R 1 (d Yb'w); Where w represents the set weight vector; γ represents the regularization coefficient; d represents the difference vector; Yb' represents the flow disturbance matrix; R represents the observation error covariance matrix; and T represents the matrix transpose operation.

[0046] Then the objective function with respect to the set weight vector The first derivative is: .

[0047] Further, step S24 above involves iteratively solving the objective function and, in conjunction with the parameter perturbation matrix, determining the corrected posterior parameter analysis field, specifically including: (1) The objective function is solved iteratively by using the sequential least squares programming algorithm to obtain the optimal weight vector that minimizes the objective function; Due to the set weight vector The optimizer operates on a set of parameters in the computational space. The optimization process does not require explicit setting of physical boundary constraints for the parameters (i.e., no boundary settings are needed). The optimizer can efficiently search across the entire real number domain, fundamentally avoiding the risk of parameter updates exceeding limits (such as negative soil moisture content). Iteration starts with the initial value of the zero vector, using gradient information to guide the search direction. Iteration stops when the change in the objective function value is less than a set threshold, thus obtaining the optimal weight vector that minimizes the objective function. Minimize the cost function J(w) and find an optimal weight vector, denoted as wopt, such that the adjusted parameters best fit the measured flow rate.

[0048] (2) Calculate the optimal correction amount of the parameters based on the parameter perturbation matrix and the optimal weight vector; Based on the normalized parameter perturbation matrix, the optimal weights are projected back into the original model parameter space, and the optimal parameter correction is calculated. .

[0049] (3) The optimal correction is superimposed on the mean of the background field parameters to obtain the corrected posterior parameter analysis field.

[0050] The parameter corrections are then added to the mean of the background parameters to obtain the optimized model parameter set. That is, the posterior parameter analysis field.

[0051] Furthermore, step S3 above updates the current parameter set based on the flow perturbation matrix and the posterior parameter analysis field, specifically including: (1) Calculate the Hessian matrix of the objective function based on the flow disturbance matrix, the observation error covariance matrix, and the regularization coefficient; Using the normalized perturbation matrix And linearized observation operators, to approximately calculate the Hessian matrix of the objective function ( The matrix characterizes the curvature of the objective function, reflecting the degree to which the observation information constrains the background error.

[0052] (2) Perform Cholesky decomposition on the Hessian matrix and construct a transformation matrix for updating the perturbation weights based on the decomposition results; right Perform Cholesky decomposition: ; In the formula, It is a lower triangular matrix.

[0053] Due to the Cholesky factor That is A matrix square root, therefore the transformation matrix It can be obtained directly by solving the system of trigonometric equations: ; If the matrix is ​​not strictly positive definite due to numerical issues, causing Cholesky decomposition to fail, then ridge regression regularization (applying a small positive increment to the diagonal) should be used first. ( Let represent the regularization coefficient, which is a preset, small positive real number. I is the identity matrix, with order . (To ensure consistency), it is forced to transform into a positive definite matrix before decomposition. The discretization of the prior set is "shrunken" and rotated to ensure that the generated posterior perturbation can correctly reflect the reduced parameter uncertainty (i.e., posterior error covariance) after fusing the observation data.

[0054] (3) Based on the transformation matrix, perform a linear transformation on the parameter perturbation matrix to obtain the posterior perturbation matrix; superimpose the posterior perturbation matrix onto the posterior parameter analysis field to obtain the updated posterior parameter set; The prior perturbation matrix is ​​linearly transformed using the transformation matrix T to obtain the posterior perturbation matrix, which is then superimposed onto the posterior parameter analysis field obtained in S43. The above generates an updated posterior parameter set containing N members. For the i-th member of the set, the calculation formula is: ; In the formula, Representative matrix The i-th column vector.

[0055] (4) Use the inverse function of the transformation function to perform inverse mapping on the analysis field members in the posterior parameter set to obtain a posterior parameter set with physical constraints; update the current parameter set through the posterior parameter set with physical constraints.

[0056] After completing the variational optimization within the computational space, in order for the assimilation results to drive the Xin'anjiang model for forecasting the next time period, the computational variables must be reduced to parameters with clear physical meaning. Based on the bijective symmetry of the aforementioned transformation function, its inverse function is used to analyze the posterior field members. Perform the reverse mapping. This step realizes the mapping and regression of data from the unbounded computational space to the bounded physical space, and its restoration formula is as follows: ; In the formula, For the i-th member of the posterior parameter set in the computational space, and These represent the lower and upper limits of the physical thresholds for the corresponding parameters. This is comprised of all members of the reconstructed posterior parameter set. This constitutes the final set of posterior parameters with physical constraints, the mean of which is... These are the final posterior analytical field parameters.

[0057] The method provided in this application is a hydrological model parameter calibration method based on four-dimensional ensemble variation, which realizes accurate correction of model parameters using measured flow data, and is conducive to improving the accuracy of hydrological simulation and flood forecasting.

[0058] Here is a specific example: This embodiment selects a typical small-to-medium-sized watershed in the humid region of southern my country as the research object. The watershed has a catchment area of ​​459.28 km2, the terrain is mainly mountainous and hilly, the climate is mild and humid, and the vegetation cover is good. It belongs to the typical Xin'anjiang model applicable area (i.e., the runoff generation area).

[0059] Given that this method has implemented a long-duration model preheating in the preceding steps, effectively eliminating random errors in the initial conditions, this embodiment adopts a 'state locking, parameter decoupling' strategy. Specifically, the hydrological state at the end of the preheating is regarded as a reliable deterministic boundary condition, and only the model physical parameters sensitive to the flow field are used as control variables for variational assimilation, thereby achieving stable convergence of the optimal solution while ensuring the consistency of the physical mechanism.

[0060] Step 1: Select 15 parameters from the Xin'anjiang model used for flood simulation as control variables to be corrected. Based on the hydrogeological characteristics of the watershed and prior experience, the physical value ranges for each parameter are preset. And prior means. To address the strict physical boundaries of hydrological parameters, a Logit space transformation is introduced: a bidirectional mapping mechanism is constructed from the physical space to the computational space. Before assimilation and update, the Logit function performs a forward transformation to map bounded parameters to unbounded variables. After the update, a function restores these bounded variables to ensure that the parameters always remain within the defined physical range.

[0061] Forward transformation:

[0062] In the formula, and These are the lower and upper limits of the physical threshold for the corresponding parameters, respectively.

[0063] Based on the watershed's hydrogeological characteristics and empirical values, the physical upper and lower bounds and prior mean values ​​of the parameters are set, as shown in Table 1 below: Table 1

[0064] The set of parameters is set to N=100 members, and prior values ​​are determined based on empirical values. The prior mean of the parameters is set based on empirical values ​​from the watershed, and a 15% relative error is set as the initial perturbation. Background field parameter sets are generated using random sampling with a Gaussian distribution. .

[0065] Step 2: Run the hydrological model to obtain simulated flow rates and construct disturbance matrices for parameters and flow rates. Input the above 100 sets of parameters into the Xin'anjiang model, using rainfall and evaporation data from June 7th to June 17th, 2015 (a typical flood event) as the driving data. Calculate the mean of the 100 flow process lines output by the model: ; Each simulated traffic with the mean Subtracting them, we obtain the flow perturbation matrix: ; Calculate the mean of the background field parameter set. Subtract the mean from each member parameter to construct the parameter perturbation matrix. ; Step 3: Obtain the measured flow data of the hydrological station during the same period Set the observation error covariance R, and based on the actual situation, set the observation error to 10% of the measured flow rate. Construct a system that includes background prior constraints and observation fitting terms.

[0066] Construct an objective function that includes background prior constraints and observation fitting terms: ; In the formula, flow vector and the average flow rate of the background field simulation The difference vector between ; This is the regularization coefficient, used to adjust the weight ratio between the background term and the observation term. In this embodiment, it is set to... The value is set to 0.1 (or between 0.01 and 1.0) to balance the parameter correction amount with the observation fitting accuracy.

[0067] Step 4: Calculate the gradient of the objective function with respect to w based on the perturbation matrix and observation error: ; The SLSQP (Sequential Least Squares) algorithm is used to search within the set space. This algorithm can efficiently utilize gradient information to find the optimal weight vector that minimizes the objective function. .

[0068] The obtained optimal weight vector Mapping back to the parameter space, we calculate the corrected optimal parameter mean, i.e., the posterior parameter analysis field: ; Step 5: Based on the set perturbation matrix and observation error, calculate the Hessian matrix of the objective function with respect to the weight vector.

[0069] Calculate the Hessian matrix of the objective function. Then, use the normalized perturbation matrix. And linearized observation operators, to approximately calculate the Hessian matrix of the objective function ( The matrix characterizes the curvature of the objective function, reflecting the degree to which the observation information constrains the background error.

[0070] right Perform Cholesky decomposition ( , (where the matrix is ​​a lower triangular matrix), we can obtain... The square root factor of the matrix. Based on the decomposition of the Hessian matrix, a transformation matrix for updating the perturbation weights is constructed. (i.e., satisfy) , which can be determined by This is efficiently obtained by solving a system of triangular equations. This transformation "shrinks" and rotates the discretization of the prior set, ensuring that the generated posterior perturbation can correctly reflect the reduced parameter uncertainty (i.e., posterior error covariance) after fusing the observation data.

[0071] Using the transformation matrix A linear transformation is performed on the prior perturbation matrix to obtain the posterior perturbation matrix, which is then superimposed on the mean of the posterior optimal parameters obtained in S43 to generate an updated posterior parameter set containing N members. For the i-th member of the set, the calculation formula is: ; In the formula, Representative matrix The i-th column vector represents the posterior perturbation of the i-th member.

[0072] When making forecasts for the next time period or outputting results, the calculated variables are restored to model parameters with clear physical meaning. Based on the bijective symmetry of the transformation function in step one, its inverse function is used to analyze the posterior field members. The inverse mapping is performed, and its restoration formula is as follows: ; In the formula, Let i be the i-th member of the posterior parameter set in the computational space. and These represent the lower and upper limits of the physical threshold for the corresponding parameters, respectively. The set consists of all the reconstructed posterior parameters. This constitutes the final set of posterior parameters with physical constraints, the mean of which is... These are the final posterior analytical field parameters.

[0073] Figure 2 This diagram illustrates the comparison between simulated and measured flow rates before and after correction using the method provided in this embodiment.

[0074] Based on the above method embodiments, this application also provides a hydrological model parameter correction device, see [link to relevant documentation]. Figure 3 As shown, the device includes: a set construction module 32, used to acquire the parameters to be corrected of the hydrological model and construct an initial background field parameter set; a parameter correction module 34, used to take the initial background field parameter set as the current parameter set and execute parameter correction steps: mapping the current parameter set to the observation space to obtain a simulated flow set; calculating the parameter perturbation matrix of the current parameter set and the flow perturbation matrix of the simulated flow set; based on the measured flow data and flow perturbation matrix of the hydrological station within the assimilation window, constructing an objective function containing two parts: background prior constraints and flow observation fitting, with the set weight vector as the independent variable; iteratively solving the objective function and determining the corrected posterior parameter analysis field by combining the parameter perturbation matrix; and a set update module 36, used to update the current parameter set based on the flow perturbation matrix and the posterior parameter analysis field, continue to execute the parameter correction steps, and realize the dynamic parameter correction process.

[0075] Furthermore, the aforementioned set construction module 32 is used to extract key physical parameters from the Xin'anjiang model as parameters to be corrected; based on the prior distribution information of the parameters to be corrected, a Monte Carlo sampling method is used to generate multiple sets of parameter vectors; and the multiple sets of parameter vectors are mapped to the computational space through a transformation function to obtain an initial set of background field parameters that satisfy the Gaussian distribution assumption.

[0076] Furthermore, the parameter correction module 34 is used to input each set of parameter vectors in the current parameter set into the observation operator determined by the Xin'anjiang model, so as to map the parameter vectors from the computation space to the observation space and obtain the simulated flow set used to characterize the simulated flow process within the assimilation window.

[0077] Furthermore, the parameter correction module 34 is used to calculate the mean value of the background field parameters of the current parameter set and the mean value of the flow of the simulated flow set; subtract the mean value of the background field parameters from the background field parameter vector corresponding to each set member in the current parameter set, and multiply by the normalization coefficient to obtain the parameter perturbation matrix; subtract the mean value of the flow from the simulated flow vector corresponding to each set member in the simulated flow set, and multiply by the normalization coefficient to obtain the flow perturbation matrix.

[0078] Furthermore, the parameter correction module 34 is used to construct an observation vector using the measured flow data of hydrological stations within the assimilation window; calculate the difference vector between the observation vector and the mean of the simulated flow sample within the assimilation window; construct an observation error covariance matrix based on the error characteristics of the observation data within the assimilation window; and construct an objective function containing background prior constraints and flow observation fitting, based on the difference vector, the observation error covariance matrix, the flow disturbance matrix, and the regularization coefficient, with the set weight vector as the independent variable.

[0079] Furthermore, the parameter correction module 34 described above is used to determine the objective function J(w) according to the following formula: J(w)=0.5γw T w+0.5(d Yb'w) T R 1 (d Yb'w); Where w represents the set weight vector; γ represents the regularization coefficient; d represents the difference vector; Yb' represents the flow disturbance matrix; R represents the observation error covariance matrix; and T represents the matrix transpose operation.

[0080] Furthermore, the parameter correction module 34 is used to iteratively solve the constructed objective function using a sequential least squares programming algorithm to obtain the optimal weight vector that minimizes the objective function; calculate the optimal correction amount of the parameters based on the parameter perturbation matrix and the optimal weight vector; and superimpose the optimal correction amount onto the mean of the background field parameters to obtain the corrected posterior parameter analysis field.

[0081] Furthermore, the aforementioned set update module 36 is used to calculate the Hessian matrix of the objective function based on the flow perturbation matrix, the observation error covariance matrix, and the regularization coefficient; perform Cholesky decomposition on the Hessian matrix, and construct a transformation matrix for updating the perturbation weights based on the decomposition result; perform a linear transformation on the parameter perturbation matrix based on the transformation matrix to obtain the posterior perturbation matrix; superimpose the posterior perturbation matrix onto the posterior parameter analysis field to obtain the updated posterior parameter set; perform inverse mapping on the analysis field members in the posterior parameter set using the inverse function of the transformation function to obtain a physically constrained posterior parameter set; and update the current parameter set using the physically constrained posterior parameter set.

[0082] The device provided in this application embodiment has the same implementation principle and technical effect as the aforementioned method embodiment. For the sake of brevity, any parts of the device embodiment not mentioned can be referred to the corresponding content in the aforementioned method embodiment.

[0083] This application also provides a hydrological model parameter correction system, such as... Figure 4 The diagram shows the structure of the system, which includes a processor 41 and a memory 40. The memory 40 stores computer-executable instructions that can be executed by the processor 41, and the processor 41 executes the computer-executable instructions to implement the above-described method.

[0084] exist Figure 4 In the illustrated embodiment, the system further includes a bus 42 and a communication interface 43, wherein the processor 41, the communication interface 43, and the memory 40 are connected via the bus 42.

[0085] The memory 40 may include high-speed random access memory (RAM) and may also include non-volatile memory, such as at least one disk storage device. Communication between this system network element and at least one other network element is achieved through at least one communication interface 43 (which can be wired or wireless), such as the Internet, wide area network, local area network, metropolitan area network, etc. The bus 42 may be an ISA (Industry Standard Architecture) bus, a PCI (Peripheral Component Interconnect) bus, or an EISA (Extended Industry Standard Architecture) bus, etc. The bus 42 can be divided into an address bus, a data bus, a control bus, etc. For ease of representation, Figure 4 The symbol is represented by a single double-headed arrow, but this does not mean that there is only one bus or one type of bus.

[0086] Processor 41 may be an integrated circuit chip with signal processing capabilities. In implementation, each step of the above method can be completed by the integrated logic circuitry in the hardware of processor 41 or by instructions in software form. Processor 41 can be a general-purpose processor, including a Central Processing Unit (CPU), a Network Processor (NP), etc.; it can also be a Digital Signal Processor (DSP), an Application Specific Integrated Circuit (ASIC), a Field-Programmable Gate Array (FPGA), or other programmable logic devices, discrete gate or transistor logic devices, or discrete hardware components. The general-purpose processor can be a microprocessor or any conventional processor. The steps of the method disclosed in the embodiments of this application can be directly implemented by a hardware decoding processor, or implemented by a combination of hardware and software modules in the decoding processor. The software modules can reside in random access memory, flash memory, read-only memory, programmable read-only memory, electrically erasable programmable memory, registers, or other mature storage media in the art. The storage medium is located in the memory, and the processor 41 reads the information in the memory and, in conjunction with its hardware, completes the steps of the method described in the foregoing embodiment.

[0087] This application also provides a computer-readable storage medium storing computer-executable instructions. When the computer-executable instructions are called and executed by a processor, the computer-executable instructions cause the processor to implement the above-described method. For specific implementation details, please refer to the foregoing method embodiments, which will not be repeated here.

[0088] The computer program products of the methods, apparatus, and electronic devices provided in the embodiments of this application include a computer-readable storage medium storing program code. The instructions included in the program code can be used to execute the methods described in the preceding method embodiments. For specific implementations, please refer to the method embodiments, which will not be repeated here.

[0089] Unless otherwise specifically stated, the relative steps, numerical expressions, and values ​​of the components and steps described in these embodiments do not limit the scope of this application.

[0090] If the aforementioned functions are implemented as software functional units and sold or used as independent products, they can be stored in a processor-executable, non-volatile, computer-readable storage medium. Based on this understanding, the technical solution of this application, in essence, or the part that contributes to the prior art, or a portion of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute all or part of the steps of the methods described in the various embodiments of this application. The aforementioned storage medium includes various media capable of storing program code, such as USB flash drives, portable hard drives, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical disks.

[0091] In the description of this application, it should be noted that the terms "center," "upper," "lower," "left," "right," "vertical," "horizontal," "inner," and "outer," etc., indicate the orientation or positional relationship based on the orientation or positional relationship shown in the accompanying drawings. They are used only for the convenience of describing this application and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation. Therefore, they should not be construed as limitations on this application. Furthermore, the terms "first," "second," and "third" are used for descriptive purposes only and should not be construed as indicating or implying relative importance.

[0092] Finally, it should be noted that the above-described embodiments are merely specific implementations of this application, used to illustrate the technical solutions of this application, and not to limit them. The protection scope of this application is not limited thereto. Although this application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that any person skilled in the art can still modify or easily conceive of changes to the technical solutions described in the foregoing embodiments, or make equivalent substitutions for some of the technical features, within the technical scope disclosed in this application. Such modifications, changes, or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of this application, and should all be covered within the protection scope of this application. Therefore, the protection scope of this application should be determined by the protection scope of the claims.

Claims

1. A method for correcting parameters of a hydrological model, characterized in that, The method includes: Obtain the parameters to be calibrated from the hydrological model and construct an initial set of background field parameters; Using the initial background field parameter set as the current parameter set, perform the parameter correction step: Mapping the current parameter set to the observation space yields a simulated flow set; Calculate the parameter perturbation matrix of the current parameter set and the flow perturbation matrix of the simulated flow set; Based on the measured flow data of hydrological stations within the assimilation window and the flow disturbance matrix, an objective function is constructed with the set weight vector as the independent variable, which includes background prior constraints and flow observation fitting. The objective function is solved iteratively, and the corrected posterior parameter analysis field is determined by combining the parameter perturbation matrix. Based on the flow disturbance matrix and the posterior parameter analysis field, the current parameter set is updated, and the parameter correction step is continued to realize the dynamic parameter correction process.

2. The method according to claim 1, characterized in that, The steps for obtaining the parameters to be calibrated from the hydrological model and constructing the initial set of background field parameters include: Key physical parameters were extracted from the Xin'anjiang model as parameters to be corrected. Based on the prior distribution information of the parameters to be corrected, multiple sets of parameter vectors are generated using the Monte Carlo sampling method; By transforming the multiple sets of parameter vectors to the computational space, an initial set of background field parameters that satisfies the Gaussian distribution assumption is obtained.

3. The method according to claim 1, characterized in that, The step of mapping the current parameter set to the observation space to obtain the simulated flow set includes: Each parameter vector in the current parameter set is input into the observation operator determined by the Xin'anjiang model to map the parameter vector from the computation space to the observation space, thereby obtaining a set of simulated flow rates used to characterize the simulated flow process within the assimilation window.

4. The method according to claim 1, characterized in that, The steps of calculating the parameter perturbation matrix of the current parameter set and the flow perturbation matrix of the simulated flow set include: Calculate the mean value of the background field parameters of the current parameter set, and the mean value of the flow rate of the simulated flow rate set; Subtract the mean value of the background field parameters from the background field parameter vector corresponding to each set member in the current parameter set, and multiply by the normalization coefficient to obtain the parameter perturbation matrix; Subtract the mean flow from the simulated flow vector corresponding to each member of the simulated flow set, and multiply by the normalization coefficient to obtain the flow perturbation matrix.

5. The method according to claim 4, characterized in that, Based on the measured flow data from hydrological stations within the assimilation window and the flow disturbance matrix, the steps for constructing an objective function comprising background prior constraints and flow observation fitting, using the set weight vector as the independent variable, include: Observation vectors are constructed using measured flow data from hydrological stations within the assimilation window; Calculate the difference vector between the observed vector and the mean of the simulated flow sample within the assimilation window period; Based on the error characteristics of the observed data within the assimilation window, construct the observation error covariance matrix; Based on the difference vector, the observation error covariance matrix, the flow disturbance matrix, and the regularization coefficient, an objective function is constructed with the set weight vector as the independent variable, which includes background prior constraints and flow observation fitting.

6. The method according to claim 5, characterized in that, Based on the difference vector, the observation error covariance matrix, the flow disturbance matrix, and the regularization coefficient, the step of constructing an objective function comprising background prior constraints and flow observation fitting, with the set weight vector as the independent variable, includes: determining the objective function according to the following formula. : ； in, Represents the set weight vector; Represents the regularization coefficient; Represents the difference vector; Represents the flow disturbance matrix; Represents the observation error covariance matrix; This represents the matrix transpose operation.

7. The method according to claim 5, characterized in that, The steps of iteratively solving the objective function and determining the corrected posterior parameter analysis field by combining the parameter perturbation matrix include: The constructed objective function is solved iteratively using a sequential least squares programming algorithm to obtain the optimal weight vector that minimizes the objective function; The optimal correction amount of the parameters is calculated based on the parameter perturbation matrix and the optimal weight vector. The optimal correction is superimposed on the mean of the background field parameters to obtain the corrected posterior parameter analysis field.

8. The method according to claim 7, characterized in that, The step of updating the current parameter set based on the flow perturbation matrix and the posterior parameter analysis field includes: Calculate the Hessian matrix of the objective function based on the flow disturbance matrix, the observation error covariance matrix, and the regularization coefficient; The Hessian matrix is ​​subjected to Cholesky decomposition, and a transformation matrix for updating the perturbation weights is constructed based on the decomposition results. Based on the transformation matrix, a linear transformation is performed on the parameter perturbation matrix to obtain the posterior perturbation matrix; The posterior perturbation matrix is ​​superimposed onto the posterior parameter analysis field to obtain the updated posterior parameter set. The analysis field members in the posterior parameter set are inversely mapped using the inverse function of the transformation function to obtain a posterior parameter set with physical constraints. The current parameter set is updated using the posterior parameter set with physical constraints.

9. A hydrological model parameter correction device, characterized in that, The device includes: The set construction module is used to obtain the parameters to be calibrated in the hydrological model and construct the initial background field parameter set; The parameter correction module is used to take the initial background field parameter set as the current parameter set and perform parameter correction steps: mapping the current parameter set to the observation space to obtain a simulated flow set; calculating the parameter perturbation matrix of the current parameter set and the flow perturbation matrix of the simulated flow set; based on the measured flow data of hydrological stations within the assimilation window and the flow perturbation matrix, constructing an objective function containing two parts: background prior constraints and flow observation fitting, with the set weight vector as the independent variable; iteratively solving the objective function and combining it with the parameter perturbation matrix to determine the corrected posterior parameter analysis field. The set update module is used to update the current parameter set based on the flow disturbance matrix and the posterior parameter analysis field, and continue to execute the parameter correction step to realize the dynamic parameter correction process.

10. A hydrological model parameter correction system, characterized in that, The system includes a processor and a memory, the memory storing computer-executable instructions executable by the processor, the processor executing the computer-executable instructions to implement the method of any one of claims 1 to 8.

Images