CN115168808A - Method for analyzing drought and waterlogging sudden turning event characteristic quantity - Google Patents
Method for analyzing drought and waterlogging sudden turning event characteristic quantity Download PDFInfo
- Publication number
- CN115168808A CN115168808A CN202210810680.8A CN202210810680A CN115168808A CN 115168808 A CN115168808 A CN 115168808A CN 202210810680 A CN202210810680 A CN 202210810680A CN 115168808 A CN115168808 A CN 115168808A
- Authority
- CN
- China
- Prior art keywords
- duration
- function
- intensity
- drought
- copula
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 74
- 230000006870 function Effects 0.000 claims abstract description 124
- 241000039077 Copula Species 0.000 claims abstract description 112
- 238000009826 distribution Methods 0.000 claims abstract description 108
- 238000005315 distribution function Methods 0.000 claims abstract description 84
- 238000013461 design Methods 0.000 claims abstract description 33
- 238000001556 precipitation Methods 0.000 claims abstract description 18
- 238000009825 accumulation Methods 0.000 claims abstract description 7
- 238000004364 calculation method Methods 0.000 claims description 30
- 230000008569 process Effects 0.000 claims description 22
- 238000011160 research Methods 0.000 claims description 14
- 230000001419 dependent effect Effects 0.000 claims description 11
- 230000008859 change Effects 0.000 claims description 8
- 230000002265 prevention Effects 0.000 claims description 7
- 230000001186 cumulative effect Effects 0.000 claims description 6
- 239000011159 matrix material Substances 0.000 claims description 6
- 238000012216 screening Methods 0.000 claims description 6
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 6
- 238000003657 Likelihood-ratio test Methods 0.000 claims description 5
- 238000007476 Maximum Likelihood Methods 0.000 claims description 5
- 238000011156 evaluation Methods 0.000 claims description 4
- 230000002123 temporal effect Effects 0.000 claims description 4
- 238000000342 Monte Carlo simulation Methods 0.000 claims description 3
- 239000000654 additive Substances 0.000 claims description 3
- 230000000996 additive effect Effects 0.000 claims description 3
- 238000013459 approach Methods 0.000 claims description 3
- 150000001875 compounds Chemical class 0.000 claims description 3
- 238000010276 construction Methods 0.000 claims description 3
- 230000005251 gamma ray Effects 0.000 claims description 3
- 238000010606 normalization Methods 0.000 claims description 3
- 230000004044 response Effects 0.000 claims description 3
- 238000005070 sampling Methods 0.000 claims description 3
- 238000000528 statistical test Methods 0.000 claims description 3
- 230000000007 visual effect Effects 0.000 claims description 3
- 238000012546 transfer Methods 0.000 claims 2
- 238000004458 analytical method Methods 0.000 description 7
- 238000012360 testing method Methods 0.000 description 7
- 238000004891 communication Methods 0.000 description 6
- 238000007689 inspection Methods 0.000 description 6
- 230000009286 beneficial effect Effects 0.000 description 3
- 238000004590 computer program Methods 0.000 description 3
- 238000000605 extraction Methods 0.000 description 3
- 230000000977 initiatory effect Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 238000010792 warming Methods 0.000 description 2
- 230000036461 convulsion Effects 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000007726 management method Methods 0.000 description 1
- 230000007704 transition Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
- G06Q10/00—Administration; Management
- G06Q10/06—Resources, workflows, human or project management; Enterprise or organisation planning; Enterprise or organisation modelling
- G06Q10/063—Operations research, analysis or management
- G06Q10/0635—Risk analysis of enterprise or organisation activities
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
- G06Q50/00—Information and communication technology [ICT] specially adapted for implementation of business processes of specific business sectors, e.g. utilities or tourism
- G06Q50/10—Services
- G06Q50/26—Government or public services
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A10/00—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE at coastal zones; at river basins
- Y02A10/40—Controlling or monitoring, e.g. of flood or hurricane; Forecasting, e.g. risk assessment or mapping
Landscapes
- Engineering & Computer Science (AREA)
- Business, Economics & Management (AREA)
- Physics & Mathematics (AREA)
- Human Resources & Organizations (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Tourism & Hospitality (AREA)
- Economics (AREA)
- Data Mining & Analysis (AREA)
- Strategic Management (AREA)
- Educational Administration (AREA)
- Pure & Applied Mathematics (AREA)
- Operations Research (AREA)
- Marketing (AREA)
- General Business, Economics & Management (AREA)
- Mathematical Physics (AREA)
- Entrepreneurship & Innovation (AREA)
- Mathematical Optimization (AREA)
- Development Economics (AREA)
- Mathematical Analysis (AREA)
- Computational Mathematics (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Evolutionary Biology (AREA)
- Bioinformatics & Computational Biology (AREA)
- Life Sciences & Earth Sciences (AREA)
- Probability & Statistics with Applications (AREA)
- Quality & Reliability (AREA)
- Algebra (AREA)
- Game Theory and Decision Science (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Health & Medical Sciences (AREA)
- General Health & Medical Sciences (AREA)
- Primary Health Care (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
The invention discloses a method for analyzing the characteristic quantity of drought and waterlogging emergency turning events. Firstly, collecting precipitation data, calculating a SWAP index, and extracting intensity and duration characteristics; then checking whether the time sequences of the two characteristics have consistency; sometimes, adopting probability distribution model fitting, and evaluating by using SBC and AIC criteria, otherwise adopting generalized accumulation model fitting, and finally obtaining an edge distribution function of intensity and duration; then, checking the consistency of the two edge distribution functions by adopting a Copula function likelihood ratio method; then, a Copula joint distribution function with constant parameters or time-varying parameters is adopted to simulate to obtain a joint distribution function with two characteristics; and finally, calculating the duration and intensity or the combined recurrence period of the ' AND ' and ' according to the combined distribution function, and calculating a combined design value of the duration and the intensity by adopting the most possible combination method, wherein the design value is the risk threshold. The method can evaluate the risk of the characteristic value of the drought and waterlogging emergency, and has very important significance for analyzing and preventing the drought and waterlogging emergency.
Description
Technical Field
The invention relates to the field of application of hydrology and water resources, in particular to a method for analyzing characteristic quantity of drought and waterlogging emergency turning events.
Background
The drought and flood sudden turning refers to the condition that a certain area or a river basin is continuously drought in the early stage and intensive rainfall occurs in the later stage, so that rivers swell and urban waterlogging and the like are changed from a drought state to a flood state in a short time, and the coexistence of drought and flood extreme events in the short time is reflected. The stability of the climate system is reduced under the background of global warming, increased extreme rainfall and enhanced water-gas circulation. Under the influence of climate change and human activities, the condition of uneven distribution of rainfall space and time is further aggravated, the frequency and the intensity of extreme disasters such as global drought, flood and the like are continuously increased, the frequency and the intensity of drought-flood turning events are also continuously increased in the world, and more students at home and abroad are attracted by wide attention in recent years.
Because the drought and waterlogging rush turning event and the characteristic change thereof have space-time discontinuity, the frequency analysis of the drought and waterlogging rush turning event characteristic is beneficial to formulating flood control and drought control measures and planning, provides a theoretical basis for preventing the drought and waterlogging rush turning disaster for the region with unbalanced water resource space-time distribution, and has important significance in regional water resource planning and management. Multiple studies show that the intensity, duration, frequency and other characteristics of drought and flood disasters have obvious variation trends driven by inconsistent climate backgrounds (such as global warming), so that when the frequency analysis is carried out on the drought and flood disasters, whether variation points exist in univariate time series or multivariate dependencies can be diagnosed. In the case of consistency, estimating a constant parameter of the variable; in the case of non-consistency, time-varying parameters of the variables are estimated. However, current research lacks consideration of inconsistent changes in drought and flood cutting events. Under the changing environment background, the inconsistency existing in the attributes of the detected drought and flood events needs to be discussed and judged, and the multivariate frequency analysis of the drought and flood rapid transition is carried out under the consideration of the potential inconsistency, so that the method has very important significance for the research and the prevention and the treatment of the meteorological disasters.
Disclosure of Invention
The method provided by the invention considers the time-varying statistical characteristics of the hydrological sequence, constructs a characteristic quantity analysis method for the drought and flood jerking events under the inconsistent background, and analyzes the characteristic quantity of the drought and flood jerking events so as to obtain more reliable risk prevention information, timely evaluates the risk of the characteristic value of the drought and flood jerking events, and has very important significance for the analysis and prevention of the drought and flood jerking disasters.
The technical scheme of the invention is as follows:
a method for analyzing the characteristic quantity of drought and flood sudden turn events comprises the following steps:
s1: determining a research area, and collecting daily scale precipitation data in the area;
s2: calculating a standardized weighted average rainfall SWAP index based on rainfall data of the drainage basin, identifying the drought and waterlogging turning process of the drainage basin, and extracting intensity and duration characteristics of drought and waterlogging turning events;
s3: adopting ADF (unit root inspection) to respectively inspect whether the time series of the intensity characteristic and the duration characteristic have consistency;
when the time sequence consistency is established, fitting intensity characteristics and duration characteristics by adopting a probability distribution model with constant parameters, and obtaining an intensity edge distribution function and a duration edge distribution function based on the evaluation of an SBC (Schwarz Bayesian Information Criterion) Criterion and an AIC (Akaike Information Criterion) Criterion;
when the time sequence consistency is not established, adopting a Generalized Accumulation Model (GAMLSS) to respectively fit the intensity characteristics and the duration characteristics to obtain an intensity edge distribution function and a duration edge distribution function;
s4: adopting a Copula-based likelihood-ratio test (CLR test) to test whether a consistency hypothesis between the intensity edge distribution function and the duration edge distribution function is established or not;
when the consistency is assumed to be set, simulating by adopting a Copula joint distribution function with constant parameters to obtain a joint distribution function of duration and intensity;
when the consistency assumption is not satisfied, simulating by adopting a Copula joint distribution function with time-varying parameters to obtain a joint distribution function with duration and intensity;
s5: and calculating a combined recurrence period of duration and intensity or and based on a combined distribution function of duration and intensity, and adopting a Most-likely combination method (Most-likely weight function) to obtain a combined design value of duration and intensity, and taking the design value as a risk threshold value for disaster control of drought and waterlogging jerk events.
Further, the calculation process of the normalized weighted average precipitation SWAP index in step S2 is as follows:
the SWAP index is an index on the premise that the drought and waterlogging state of the day is influenced by the early drought and waterlogging state and the rainfall of the day, the calculation method is to assume that a weighted average rainfall (WAP) sequence of the same day obeys Gamma distribution, fit the weighted average rainfall sequence through the Gamma distribution, and finally perform standard normalization processing on the cumulative probability distribution, wherein the calculation formula of the weighted average rainfall sequence is as follows:
w n ＝(1-a)a n-1
of the WAP i Representing a weighted mean precipitation sequence, P n Denotes the precipitation on day n, w n Represents P n A is a parameter of the weight decay with time, a n The weight of the nth day is an attenuation parameter along with time, N is the number of days affected by early precipitation, and a =0.9 and N =45 are selected according to the local soil-water system in the research area.
Further, in step S2, based on the calculated SWAP index, a running theory is used to identify and extract the drought-waterlogging and rush-turning process of the drainage basin, and the specific steps are as follows:
setting index critical value X for a certain drought and waterlogging index time sequence 0 、X 1 And a time length threshold value D 0 、D 1 And D 2 When the index is continuously lower than X 0 Time of not less than D 0 When drought occurs, it is considered to be drought; when the index is continuously higher than X 1 Time of not less than D 1 Then, the flood is considered to appear; when the time interval between drought and flood is not more than D 2 When the vehicle is in motion, the vehicle is considered to form a drought and waterlogging rush turning event; selecting a pair of upper bound of drought grade and lower bound of waterlogging grade according to drought and waterlogging index division standardCorresponding numerical value as X 0 And X 1 I.e. X 0 = -1 and X 1 =1; determining a recognition duration threshold D based on a daily index by referring to a run-length theory recognition process of various types of sudden turning event recognition research 0 ＝D 1 Time length threshold D of rush hour between drought state and flood state for =10 days 2 And (5) acquiring duration and intensity of each drought and waterlogging turning event according to the extracted drought and waterlogging turning events.
Further, in step S3, when the time sequence consistency is established, the concrete process of fitting to obtain the intensity edge distribution function and the duration edge distribution function is as follows;
selecting a probability distribution model which is commonly used in meteorological hydrological research as an alternative model, wherein the probability distribution model specifically comprises 7 probability distribution models of an index (Exponential, exp), a Gamma (Gamma), a Generalized Gamma (Gamma 3) with three parameters, a Normal (Normal), a Log-Normal (lornorm), a Logistic (logarithm, lorm), a Logistic (logigis, logis) and a Weibull (Weibull), fitting is carried out on the severity characteristic and the duration characteristic through the 7 alternative models respectively, and the fitting goodness of the 7 alternative models is evaluated through an SBC (SBC) criterion;
the calculation formula of SBC is as follows:
SBC＝GD+log(ξ)×df＝-2 log L(θ)+log(ξ)×df
in the formula, GD is global fitting deviation, df is number of degrees of freedom, theta and xi respectively represent estimated values and parameter numbers of probability distribution parameters, and L (-) represents a log-likelihood function value;
in order to avoid over-fitting of distribution, the goodness of fit of the 7 alternative models is evaluated by considering the AIC criterion;
the calculation formula for AIC is as follows:
AIC＝2ξ-2 log L(θ)
and finally, selecting a distribution model which enables the calculated values of the AIC and the SBC to be minimum simultaneously from all the tested candidate models as a relatively optimal model, and using the relatively optimal model as a strength edge distribution function or a duration edge distribution function.
Further, for the alternative normal and lognormal probability distribution models, the definition interval of the variable is [ - ∞, + ∞ ], and the intensity feature and duration feature to be tested are both greater than 0, so the definition interval of the 2 alternative models needs to be corrected to a positive interval, and the corrected cumulative probability is expressed by the following formula:
where f (x | θ) is a probability density function and θ is a probability distribution parameter.
Further, in step S3, when the time sequence consistency is not established, the concrete process of obtaining the intensity edge distribution function and the duration edge distribution function by fitting is as follows;
fitting the edge distribution using a Generalized Additive Models (GAMLSS) that can comprehensively consider the location, scale and shape parameters, in which the response variable Y is assumed
Individual independent observed value
Obeying probability distribution f (y) i |Θ i )；
Wherein, theta i ＝(μ i ,σ i ,v i ,τ i ) The vector is composed of probability distribution parameters, mu and sigma are position and scale parameters, and v and tau are shape parameters, and represent skewness and kurtosis of distribution respectively;
establishing a relation between the probability distribution parameters and the explanatory variables thereof through a connection function, wherein the expression is as follows:
namely:
in the formula, g k (. -) represents a join function; theta k Mu, sigma, v, tau are lengths of
A column vector of (a);
X k is a known matrix, corresponds
Interpretation variable values for the time periods; k represents the sequence number of the sequence to be analyzed;
is the parameter vector to be solved with dimension J k ；Z jk Is a known design matrix; gamma ray jk Is a column vector composed of random variables.
In the invention, the identified Drought and flood sharp turn occurrence time (Drought onset) is used as an explanatory variable, a functional relation between the duration of the sharp turn and the intensity probability distribution parameter and the Drought and flood sharp turn occurrence time is established, and the probability distribution parameter is expressed as a cubic spline function of the Drought and flood sharp turn initiation time by adopting a semi-parameter accumulation form of a GALSS model, thereby being beneficial to capturing the possible nonlinear relation. The method takes the maximum likelihood function as a target, calculates the time-varying parameters of the alternative probability distribution, and screens the optimal intensity marginal probability function and duration marginal probability function from 7 alternative distribution models according to the SBC criterion and the AIC criterion.
Further, in step S4, a Copula function likelihood ratio method is used to check whether the consistency assumption between the intensity edge distribution function and the duration edge distribution function is true or not as follows;
the basic assumption of the Copula function likelihood ratio method is that the type of the Copula function is fixed and is unchanged, and the parameter theta is detected c Analyzing the variable points of the multivariate hydrological sequence dependent structural strength;
known as y 1 ,y 2 ,...,y d Is a d-dimensional random variableA time series of observations, wherein y i ＝(y 1,i ,y 2,i ,...,y n,i ) I =1, 2.. D, the population where the observations lie obeys a probability distribution F (y) 1 |θ c,1 ),F(y 2 |θ c,2 ),...,F(y n |θ c,n )，θ c,i The Copula function parameters corresponding to different time periods, in the Copula function likelihood ratio test, the original assumption is that a multivariable dependent structure has no variable point (i.e. all observed values obey the same probability distribution), and the expression formula is as follows:
H 0 :θ c,1 ＝θ c,2 ＝...＝θ c,n ＝η 0
alternative hypothesis H 1 Comprises the following steps: presence of k * And is
So that theta c,1 ＝θ c,2 ＝...＝θ c,k* ＝η 1 And θ c,k* ＝...＝θ c,n ＝η 2 While eta 1 ≠η 2 (ii) a Constructing likelihood ratio Lambda based on Copula function k When a is k If the value is less than the given threshold value, the original hypothesis H is rejected 0 And considering that the multivariate dependent structure has a variable point k * ；
In the formula, L n (. H) represents a likelihood function for the entire sample sequence; l is k (. A) and
likelihood functions before and after the check point k, respectively; u. of i Representing the edge distribution probability of the bivariate sequence; c (-) represents the Copula probability density function;
and
are Copula function parameters η, respectively 0 、η 1 And η 2 Maximum likelihood estimate of Λ, rewrite Λ k In logarithmic form, we obtain:
further, the statistical test quantity in the form of construction logarithms is as follows:
when Z is n Greater than a threshold corresponding to a given level of significance, rejecting the original hypothesis H 0 Likelihood ratio statistics
Obeying an asymptotic distribution:
where z represents a given threshold, and when z approaches infinity, h = l = [ ln (n)] 3/2 ；
h. l and O represent parameters of asymptotic distribution; p is the number of parameters of the Copula function where there is one change point; will be provided with
Substituting the formula, if the calculated P value is smaller than the threshold, rejecting the original hypothesis under the significance level of the threshold, and if the consistency hypothesis is not established, otherwise, establishing the consistency hypothesis.
Further, in step S4, the joint distribution function of duration and intensity established according to the Copula function is expressed as:
P(D,S)＝C(u D ,u S |θ c )
in the formula, D and S respectively represent intensity and duration; u. u D ,u S Representing intensity and durationThe edge distribution of (2);
similar to the GAMLSS model, in order to estimate the time-varying parameters of Copula, explanation variables are introduced for describing the functional relationship between the time-varying parameters of Copula and the explanation variables;
g(θ c )＝β 0 +Xβ 1
wherein g (θ) c ) When the copula parameter is theta c A join function of (a);
θ c ＝(θ c,1 ,θ c,2 ,...,θ c,n ) T is a column vector composed of Copula time-varying parameters;
X＝(X 1 ,X 2 ,...,X n ) T n observations representing an explanatory variable; beta is a 0 And β represents the parameter vector to be solved;
the following expression of the connection function is given according to the condition whether the consistency assumption of the relation between duration and intensity is true or not:
(1) When the two-variable relation keeps consistent, the Copula parameter is a constant value which does not change along with time, and the connection function is expressed as:
g(θ c )＝constant
(2) When the consistency assumption is broken, the Copula parameter is a function of time, thus introducing the time of occurrence of the drought event as an explanatory variable, the join function is rewritten as:
g(θ c )＝β 0 +time·β 1
furthermore, a marginal function inference method is adopted, and a Copula parameter theta is solved by maximizing a log-likelihood function c,t The formula is as follows:
in the formula u D,i Intensity edge distribution representing the ith time period; u. of S,i An edge distribution representing duration of the ith time period; d i Indicating the intensity of the ith time period; s i Represents the duration of the ith time period; theta.theta. D Parameter representing intensity edge distributionCounting; theta S A parameter representing a duration edge distribution; theta D,i A parameter representing intensity edge distribution for the ith time period; theta S,i A parameter representing a temporal edge distribution of the ith time period; c (-) is a probability density function of Copula; f. of D (. And f) S (. H) probability density functions representing duration, intensity edge distributions, respectively; in the marginal function inference method (IFM), since the maximum values of the second and third terms on the right side of the equal sign of the formula are obtained by fitting the edge distribution of duration and intensity, the optimal beta can be obtained by maximizing the first term on the right side of the equal sign 0 And beta 1 And will be beta 0 And beta 1 Substituting the estimated value of the Copula function into the estimated value of the Copula function;
in the formula, time i Indicating the occurrence time of the ith period.
According to the invention, common Copula functions can be divided into an ellipse (elliptic) Copula, an Archimedean (Archimedean) Copula, an Extreme Value (Extreme Value) Copula, a mixed Copula and the like according to the difference of dependent structures among variables under description limit conditions, and the ellipse Copula and the Archimedean Copula are widely applied in the hydrological field at present. The invention adopts an elliptic Copula (Gaussian Copula) and four Archimedes Copula functions (Clayton Copula, frank Copula, gumbel Copula and Joe Copula) as well as five Copula functions as alternative Copula functions of alternative single parameters. And after fitting, performing goodness-of-fit inspection on the selected Copula function. According to a Bayesian Information Criterion (BIC), a function with the minimum BIC value is selected as an optimal joint distribution Copula function, and a joint distribution function of duration and intensity is established.
Further, in step S5, the process of calculating the duration and intensity "or" and "combined recurrence period is as follows:
univariate characteristic P (x is more than or equal to x) of drought and waterlogging sharp turn * ) The recurrence period of (c) is calculated as follows:
wherein T is the recurrence period of drought and waterlogging sudden-turn event, F (x) * ) Is x * An event probability function; interval is the time interval of all drought and waterlogging turning events, and the unit is year, namely the time length from the beginning of a certain drought and waterlogging turning event to the beginning of the next drought and waterlogging turning event; if the duration is N year Annual sampling to obtain N alter The interval calculation formula of the jerky event is as follows:
the multivariate drought and flood rush turn frequency analysis has a plurality of combined joint recurrence periods, including 'OR' the joint recurrence period and 'the joint recurrence period', OR 'the joint recurrence period' means that one OR more random variables of the drought and flood rush turn events exceed a given threshold value to obtain the recurrence period, and the recurrence period is represented by an OR operator '. V'; the 'AND' joint recurrence period 'refers to the recurrence period calculated when the random variables of the drought-waterlogging AND sudden-turning events exceed a given threshold, so that the' AND 'recurrence period' is also called 'joint exceeding recurrence period' or 'co-occurrence recurrence period' AND is represented by an AND operator 'reverse-inverted-V';
the calculation formula of the reappearance period under different visual angles is as follows;
under univariate viewing angle:
in the formula, T S 、T D Duration, intensity, respectively, of the reconstruction period under univariate viewing angle, F S (s) and F D (d) Respectively representing duration and intensityThe probability of birth;
under a bivariate viewing angle:
in the formula (I), the compound is shown in the specification,
representing the "or" joint recurrence period under both variable perspectives of duration and intensity,
denotes the "and" joint recurrence period, F, at both time and intensity variables S,D (s, d) is the joint occurrence probability of duration and intensity;
based on the various calculations described above, the duration and intensity "or" and "joint recurrence period is obtained.
Further, in step S5, the most probable combination method is adopted to calculate the joint design value of duration and intensity as follows:
the joint design value of drought and flood racing events refers to a multivariable combination (u, v) when the combined racing density f (u, v) composed of the racing characteristics reaches a maximum value at a certain design recurrence period level ml ,v ml ) I.e., the combination is most likely to occur, and therefore solved by constructing the following equation:
[u ml ,v ml ]＝argmaxf(u,v)
f(u,v)＝c(u,v)·f(u)·f(v)
wherein f (u, v) is a joint probability density function of double variables of sharp turn, c (u, v) is a probability density function corresponding to joint Copula distribution, f (u) and f (v) are edge density functions of variables, [ u, v ] ml ,v ml ]Combining the joint design values of the multiple variables;
the calculation steps for obtaining the joint design value by adopting the most probable combination method are as follows:
(1) Simulating a plurality of variable combinations (u) by adopting a Monte Carlo method based on a Copula joint distribution function ml ,v ml ) The number of the groups is 100000;
(2) According to the calculation formula of the ' OR ' and ' recurrence period, the joint recurrence period corresponding to the bivariate combination of the last step is calculated;
(3) According to the determined recurrence period level, screening to meet the precision requirement (e = 10) -4 ) Wherein the variable combination with the maximum joint probability density value is the joint design value.
Compared with the prior art, the invention has the following beneficial effects:
the invention provides a method for analyzing drought and flood jerking event characteristic quantities based on a non-uniformity theory, which is characterized in that a time-varying parameter Copula function is used for establishing combined distribution between duration and intensity with strong dependency of drought and flood jerking characteristics under a non-uniformity background, and a combined design value of the drought and flood jerking characteristics under each recurrence period is calculated and is used as a risk threshold value of the drought and flood jerking event, so that the risk of the drought and flood jerking event characteristic values can be further evaluated, and the method has very important significance for meteorological disaster research and prevention and control.
Drawings
FIG. 1 is a schematic flow chart of the present invention.
Detailed Description
The drawings are for illustrative purposes only and are not to be construed as limiting the patent; for the purpose of better illustrating the present embodiments, certain elements of the drawings may be omitted, enlarged or reduced, and do not represent the size of an actual product; it will be understood by those skilled in the art that certain well-known structures in the drawings and descriptions thereof may be omitted. The positional relationships depicted in the drawings are for illustrative purposes only and are not to be construed as limiting the present patent.
Example 1:
as shown in fig. 1, a method for analyzing the characteristic quantity of drought, waterlogging and sudden turning events comprises the following steps:
s1: determining a research area, and collecting daily scale precipitation data in the area;
s2: calculating a standardized weighted average rainfall SWAP index based on rainfall data of the drainage basin, identifying the drought and waterlogging turning process of the drainage basin, and extracting intensity and duration characteristics of drought and waterlogging turning events;
s3: adopting ADF (unit root inspection) to respectively inspect whether the time series of the intensity characteristics and the duration characteristics have consistency;
when the time sequence consistency is established, fitting intensity characteristics and duration characteristics by adopting a probability distribution model with constant parameters, and obtaining an intensity edge distribution function and a duration edge distribution function based on the evaluation of an SBC (Schwarz Bayesian Information Criterion) Criterion and an AIC (Akaike Information Criterion) Criterion;
when the time sequence consistency is not established, adopting a Generalized Accumulation Model (GAMLSS) to respectively fit the intensity characteristics and the duration characteristics to obtain an intensity edge distribution function and a duration edge distribution function;
s4: adopting a Copula-based likelihood-ratio test (CLR test) to test whether the consistency hypothesis between the intensity edge distribution function and the duration edge distribution function is established;
when the consistency is assumed to be set, simulating by adopting a Copula joint distribution function with constant parameters to obtain a joint distribution function of duration and intensity;
when the consistency assumption is not satisfied, simulating by adopting a Copula joint distribution function with time-varying parameters to obtain a joint distribution function with duration and intensity;
s5: calculating the duration and intensity or the combined recurrence period of 'and' based on the combined distribution function of the duration and the intensity, and adopting a Most-probable combination method (Most-likelyweight function) to obtain a combined design value of the duration and the intensity, wherein the design value is used as a risk threshold value for disaster prevention and control of drought and flood sudden-turn events.
The meteorological data used in the embodiment can be derived from a national weather science data center product, namely a Chinese ground climate data daily value data set (V3.0), the selected river basin can be a Zhujiang river basin, and the daily data of precipitation of 75 meteorological sites in the Zhujiang river basin in 1961-2020 is collected.
In this embodiment, the calculation process of the normalized weighted average precipitation SWAP index in step S2 is as follows:
the SWAP index is an index on the premise that the drought and waterlogging state of the day is influenced by the early drought and waterlogging state and the rainfall of the day, the calculation method is to assume that a weighted average rainfall (WAP) sequence of the same day obeys Gamma distribution, fit the weighted average rainfall sequence through the Gamma distribution, and finally perform standard normalization processing on the cumulative probability distribution, wherein the calculation formula of the weighted average rainfall sequence is as follows:
w n ＝(1-a)a n-1
in the form of WAP i Representing a weighted mean precipitation sequence, P n Denotes the precipitation on day n, w n Is represented by P n A is a parameter of the weight decay with time, a n The weight of the nth day is an attenuation parameter along with time, N is the number of days affected by early precipitation, and a =0.9 and N =45 are selected according to the local soil-water system in the research area.
In the implementation, in step S2, based on the calculated SWAP index, the drought-waterlogging and rush-turning process of the drainage basin is identified and extracted by using a run-length theory, and the specific steps are as follows:
setting an index critical value X for a certain drought-waterlogging index time sequence 0 、X 1 And a time length threshold value D 0 、D 1 And D 2 When the index is continuously lower than X 0 Time of not less than D 0 When drought occurs, it is considered to occur; when the index is continuously higher than X 1 Time of not less than D 1 When the flood occurs, the flood is considered to occur; when the time interval between drought and flood is not more than D 2 When the vehicle is in motion, the vehicle is considered to form a drought and waterlogging rush turning event; according to drought-waterlogging index division standard, selecting numerical value corresponding to drought grade upper bound and waterlogging grade lower bound as X 0 And X 1 I.e. X 0 = -1 and X 1 =1; determining a recognition duration threshold D based on a daily index by referring to a run-length theory recognition process of various types of sudden turning event recognition research 0 ＝D 1 =10 days, rush hour duration threshold D between drought and flood conditions 2 And (5) acquiring duration and intensity of each drought and waterlogging turning event according to the extracted drought and waterlogging turning events.
In the present embodiment, in step S3, when the time sequence consistency is established, the concrete process of obtaining the intensity edge distribution function and the duration edge distribution function by fitting is as follows;
selecting a probability distribution model which is commonly used in meteorological hydrological research as an alternative model, wherein the probability distribution model specifically comprises 7 probability distribution models of an index (Exponential, exp), a Gamma (Gamma), a Generalized Gamma (Gamma 3) with three parameters, a Normal (Normal), a Log-Normal (lornorm), a Logistic (logarithm, lorm), a Logistic (logigis, logis) and a Weibull (Weibull), fitting is carried out on the severity characteristic and the duration characteristic through the 7 alternative models respectively, and the fitting goodness of the 7 alternative models is evaluated through an SBC (SBC) criterion;
the calculation formula of SBC is as follows:
SBC＝GD+log(ξ)×df＝-2 log L(θ)+log(ξ)×df
in the formula, GD is global fitting deviation, df is the number of degrees of freedom, theta and xi respectively represent the estimated value and the number of parameters of probability distribution parameters, and L (-) represents a log-likelihood function value;
in order to avoid over-fitting of distribution, the fitting goodness of the 7 alternative models is evaluated by considering the AIC criterion;
the calculation formula for AIC is as follows:
AIC＝2ξ-2logL(θ)
and finally, selecting a distribution model which enables the calculated values of the AIC and the SBC to be minimum simultaneously from all the candidate models passing the test as a relatively optimal model, and using the relatively optimal model as an intensity edge distribution function or a duration edge distribution function.
In this embodiment, for the alternative normal and lognormal probability distribution models, the defined interval of the variables is [ - ∞, + ∞ ] and the intensity and duration characteristics to be tested are both greater than 0, so that the defined interval of the 2 alternative models is modified to a positive interval, and the modified cumulative probability is expressed by the following formula:
in the formula, f (x | θ) is a probability density function, and θ is a probability distribution parameter.
In the present embodiment, in step S3, when the time sequence consistency is not established, the concrete process of obtaining the intensity edge distribution function and the duration edge distribution function by fitting is as follows;
fitting the edge distribution using a Generalized Additive Models (GAMLSS) that can comprehensively consider the location, scale and shape parameters, in which the response variable Y is assumed
Individual independent observed value
Obeying probability distribution f (y) i |Θ i )；
Wherein, theta i ＝(μ i ,σ i ,v i ,τ i ) The vector is composed of probability distribution parameters, mu and sigma are position and scale parameters, and v and tau are shape parameters, and represent skewness and kurtosis of distribution respectively;
establishing a relation between the probability distribution parameters and the explanatory variables thereof through a connection function, wherein the expression is as follows:
namely:
in the formula, g k (. Cndot.) represents a join function; theta k Mu, sigma, v, tau are lengths of
A column vector of (a);
X k is a known matrix, corresponds
Interpretation variable values for the time periods; k represents the sequence number of the sequence to be analyzed;
is the parameter vector to be solved with dimension J k ；Z jk Is a known design matrix; gamma ray jk Is a column vector composed of random variables.
In the implementation, the identified Drought and flood rush turn occurrence time (Drought onset) is used as an explanatory variable, a functional relation between the rush turn duration and intensity probability distribution parameters and the Drought and flood rush turn occurrence time is established, and the probability distribution parameters are expressed as a cubic spline function of the Drought and flood rush turn initiation time by adopting a semi-parameter accumulation form of a GALSS model, so that the possible nonlinear relation can be captured. The invention takes a maximum likelihood function as a target, calculates time-varying parameters of alternative probability distribution, and screens optimal intensity marginal probability functions and duration marginal probability functions from 7 alternative distribution models according to an SBC (statistical distribution based) criterion and an AIC (empirical mode classification) criterion.
In the present embodiment, in step S4, a Copula function likelihood ratio method is used to check whether the consistency assumption between the intensity edge distribution function and the duration edge distribution function is true or not as follows;
the basic assumption of the Copula function likelihood ratio method is that the type of the Copula function is fixed and is unchanged, and the parameter theta is detected c Analyzing the variable points of the multivariate hydrological sequence dependent structural strength;
known as y 1 ,y 2 ,...,y d Is a time series of d-dimensional random variable observations, where y i ＝(y 1,i ,y 2,i ,...,y n,i ) I =1, 2.. D, the population where the observations lie obeys a probability distribution F (y) 1 |θ c,1 ),F(y 2 |θ c,2 ),...,F(y n |θ c,n )，θ c,i The Copula function parameters corresponding to different time periods are adopted, in the Copula function likelihood ratio test, the original assumption is that a multivariate dependent structure has no variable point (that is, all observed values obey the same probability distribution), and the expression formula is as follows:
H 0 :θ c,1 ＝θ c,2 ＝...＝θ c,n ＝η 0
alternative hypothesis H 1 Comprises the following steps: presence of k * And is
So that theta c,1 ＝θ c,2 ＝...＝θ c,k* ＝η 1 And theta c,k* ＝...＝θ c,n ＝η 2 While eta 1 ≠η 2 (ii) a Constructing likelihood ratio Lambda based on Copula function k When a is k Below a given threshold, rejecting the original hypothesis H 0 And considering that the multivariate dependent structure has a variable point k * ；
In the formula, L n (. O) a likelihood function representing the entire sequence of samples; l is k (. A) and
likelihood functions before and after the check point k, respectively; u. of i Representing the edge distribution probability of the bivariate sequence; c (-) represents the Copula probability density function;
and
are Copula functions, respectivelyParameter eta 0 、η 1 And η 2 Maximum likelihood estimate of, rewrite of Λ k In logarithmic form, we obtain:
further, the statistical test quantity in the form of construction logarithms is as follows:
when Z is n Above a threshold corresponding to a given level of significance, reject the original hypothesis H 0 Likelihood ratio statistics
Obey an asymptotic distribution:
where z represents a given threshold, and as z approaches infinity, h = l = [ ln (n)] 3/2 ；
h. l and O represent parameters of asymptotic distribution; p is the number of parameters of the Copula function where there is one change point; will be provided with
Substituting the formula, if the calculated P value is smaller than the threshold, rejecting the original hypothesis under the significance level of the threshold, and if the consistency hypothesis is not established, otherwise, establishing the consistency hypothesis.
In this embodiment, in step S4, the joint distribution function of duration and intensity established according to the Copula function is expressed as:
P(D,S)＝C(u D ,u S |θ c )
in the formula, D and S respectively represent intensity and duration; uD, uS represents the edge distribution of intensity and duration;
similar to the GALSS model, in order to estimate the time-varying parameters of the Copula, explanation variables are introduced for describing the functional relationship between the time-varying parameters of the Copula and the explanation variables;
g(θ c )＝β 0 +Xβ 1
wherein g (θ) c ) Indicates when the copula parameter is θ c A connection function of (a);
θ c ＝(θ c,1 ,θ c,2 ,...,θ c,n ) T is a column vector composed of Copula time-varying parameters;
X＝(X 1 ,X 2 ,...,X n ) T n observations representing an explanatory variable; beta is a 0 And β represents the parameter vector to be solved;
the following expression of the connection function is given according to the condition whether the consistency assumption of the relation between duration and intensity is true or not:
(1) When the two-variable relationship is consistent, the Copula parameter is a constant value which does not change with time, and the connection function is expressed as:
g(θ c )＝constant
(2) When the consistency assumption is broken, the Copula parameter is a function of time, thus introducing the time of occurrence of the drought event as an explanatory variable, the join function is rewritten as:
g(θ c )＝β 0 +time·β 1
further, a marginal function inference method is adopted, and a Copula parameter theta is solved by maximizing a log-likelihood function c,t The formula is as follows:
in the formula u D,i Representing an intensity edge distribution for the ith period; u. of S,i An edge distribution representing duration of the ith time period; d i Indicating the intensity of the ith time period; s i Represents the duration of the ith time period; theta D A parameter representing intensity edge distribution; theta S A parameter representing a temporal edge distribution; theta D,i A parameter representing intensity edge distribution for the ith time period; theta S,i A parameter representing a temporal edge distribution of the ith time period; c (-) is a probability density function of Copula; f. of D (. And f) S (. H) probability density functions representing duration, intensity edge distributions, respectively; in the marginal function inference method (IFM), since the maximum values of the second and third terms on the right side of the equal sign of the formula are obtained by fitting the edge distribution of duration and intensity, the optimal beta can be obtained by maximizing the first term on the right side of the equal sign 0 And beta 1 And will be beta 0 And beta 1 Substituting the estimated value of the Copula function into the estimated value of the Copula function;
in the formula, time i Indicating the occurrence time of the ith period.
According to the invention, common Copula functions can be divided into an ellipse (elliptic) Copula, an Archimedean (Archimedean) Copula, an Extreme Value (Extreme Value) Copula, a mixed type Copula and the like according to the difference of dependent structures among variables under the description limit condition, and the ellipse Copula and the Archimedean Copula are widely applied in the hydrology field at present. The invention adopts an elliptic Copula (Gaussian Copula) and four Archimedes Copula functions (Clayton Copula, frank Copula, gumbel Copula and Joe Copula) as well as five Copula functions as alternative Copula functions of alternative single parameters. And after fitting, the selected Copula function needs to be subjected to goodness-of-fit inspection. According to Bayesian Information Criterion (BIC), a function with the minimum BIC value is selected as an optimal joint distribution Copula function, and a joint distribution function of duration and intensity is established.
In this embodiment, in step S5, the process of calculating the duration and intensity "or" and "joint recurrence period is as follows:
univariate characteristic P (x is more than or equal to x) of drought and flood sudden turning * ) The recurrence period of (c) is calculated as follows:
wherein T is the recurrence period of drought and waterlogging sudden-turn event, F (x) * ) Is x * An event probability function; the interval is the average time interval of all drought and flood turning events, and the unit is year, namely the time length from the beginning of a certain drought and flood turning event to the beginning of the next drought and flood turning event; if the duration is N year Annual sampling to obtain N alter The interval calculation formula of the scene sharp turn event is as follows:
the multivariate drought and water-logging jerky frequency analysis has a plurality of combined joint recurrence periods, including 'OR' the joint recurrence period and 'the joint recurrence period' OR 'the joint recurrence period' refers to the recurrence period calculated by one OR more random variables of drought and water-logging jerky events exceeding a given threshold value, and is represented by an OR operator V-V; the 'AND' joint recurrence period 'refers to the recurrence period calculated when the random variables of the drought-waterlogging AND sudden-turning events exceed a given threshold, so that the' AND 'recurrence period' is also called 'joint exceeding recurrence period' or 'co-occurrence recurrence period' AND is represented by an AND operator 'reverse-inverted-V';
the calculation formula of the reappearance period under different visual angles is as follows;
under univariate viewing angle:
in the formula, T S 、T D Duration, intensity, respectively, for the reconstruction period under univariate viewing angle, F S (s) and F D (d) Respectively representing the probability of occurrence of duration and intensity;
under a bivariate viewing angle:
in the formula (I), the compound is shown in the specification,
representing the "or" joint recurrence period under both variable perspectives of duration and intensity,
denotes the "and" joint recurrence period, F, at both time and intensity variables S,D (s, d) is the joint occurrence probability of duration and intensity;
the duration and severity "or" and "combined recurrence period is derived based on the various calculations described above.
In this embodiment, in step S5, the process of obtaining the joint design value of duration and intensity by using the most probable combination method is as follows:
the joint design value of drought and flood racing events refers to a multivariable combination (u, v) when the combined racing density f (u, v) composed of the racing characteristics reaches a maximum value at a certain design recurrence period level ml ,v ml ) I.e., the combination is most likely to occur, and therefore solved by constructing the following equation:
[u ml ,v ml ]＝argmax(u,v)
f(u,v)＝c(u,v)·f(u)·f(v)
wherein f (u, v) is a joint probability density function of sharp turn bivariate, c (u, v) is a probability density function corresponding to joint Copula distribution, f (u) and f (v) are edge density functions of variables, [ u, v ] ml ,v ml ]Combining the joint design values of the multiple variables;
the calculation steps for obtaining the joint design value by adopting the most probable combination method are as follows:
(1) Simulating a plurality of variable combinations (u) by adopting a Monte Carlo method based on a Copula joint distribution function ml ,v ml ) The number of the groups is 100000;
(2) According to the calculation formula of the ' OR ' and ' recurrence period, the joint recurrence period corresponding to the bivariate combination of the last step is calculated;
(3) And screening to meet the precision requirement according to the determined recurrence period level (e = 10) -4 ) Wherein the variable combination with the maximum joint probability density value is the joint design value.
The invention provides a method for analyzing drought and flood jerking event characteristic quantities based on a non-uniformity theory, which is characterized in that a time-varying parameter Copula function is used for establishing combined distribution between duration and intensity with strong dependency of drought and flood jerking characteristics under a non-uniformity background, and a combined design value of the drought and flood jerking characteristics under each recurrence period is calculated and is used as a risk threshold value of the drought and flood jerking event, so that the risk of the drought and flood jerking event characteristic values can be further evaluated, and the method has very important significance for meteorological disaster research and prevention and control.
Example 2:
the embodiment provides a device for analyzing the characteristic quantity of drought and flood sudden turning events, which comprises an acquisition module, a characteristic extraction module, a screening module, a detection module and an output module which are sequentially in communication connection;
wherein, the collection module: the system is used for collecting daily scale precipitation data in a region;
a feature extraction module: calculating a standardized weighted average rainfall SWAP index based on the data of the acquisition module, identifying the drought and waterlogging rush-turn process of the drainage basin, and extracting the intensity and duration characteristics of the drought and waterlogging rush-turn event;
a screening module: the ADF (unit root inspection) is adopted to respectively inspect whether the time series of the intensity characteristics and the duration characteristics have consistency or not according to the intensity characteristics and the duration characteristics obtained by the characteristic extraction module;
when the time sequence consistency is established, fitting the intensity characteristic and the duration characteristic by adopting a probability distribution model with constant parameters, and obtaining an intensity edge distribution function and a duration edge distribution function based on the evaluation of the SBC (Schwarz Bayesian Information Criterion) Criterion and the AIC (Akaike Information Criterion) Criterion;
when the time sequence consistency is not established, adopting a Generalized Accumulation Model (GAMLSS) to respectively fit the intensity characteristics and the duration characteristics to obtain an intensity edge distribution function and a duration edge distribution function;
a checking module: the edge distribution function module is used for adopting a Copula-based likelihood ratio method (CLR test) to test whether a consistency hypothesis between the intensity edge distribution function and the duration edge distribution function is established or not according to the intensity edge distribution function and the duration edge distribution function obtained by the screening module;
when the consistency is assumed to be set, simulating by adopting a Copula joint distribution function with constant parameters to obtain a joint distribution function of duration and intensity;
when the consistency assumption is not satisfied, simulating by adopting a Copula joint distribution function with time-varying parameters to obtain a joint distribution function with duration and intensity;
an output module: and calculating the combined recurrence period of the duration and the intensity or the 'and' based on the combined distribution function of the duration and the intensity obtained by the inspection module, and obtaining a combined design value of the duration and the intensity by adopting a Most-probable combination method (Most-likely weight function), and finally outputting the design value as a risk threshold value for preventing and controlling the drought and waterlogging sudden-turn events.
The modules can be integrated in a system for implementation, and the result is output to a display screen through an output port for display.
Example 3:
the embodiment provides an electronic device, which comprises a processor, a communication interface, a memory and a communication bus, wherein the processor and the communication interface are used for completing mutual communication by the memory through the communication bus; a memory for storing a computer program; and a processor for executing the program stored in the memory to implement the method for analyzing the characteristic quantity of the drought/flood emergency in embodiment 1.
The present embodiment also provides a computer-readable storage medium, in which a computer program is stored, and when being executed by a processor, the computer program implements the method for analyzing the characteristic quantity of the drought/flood rush turning event.
It should be understood that the above-described embodiments of the present invention are merely examples for clearly illustrating the present invention and are not intended to limit the embodiments of the present invention. Other variations and modifications will be apparent to persons skilled in the art in light of the above description. And are neither required nor exhaustive of all embodiments. Any modification, equivalent replacement, and improvement made within the spirit and principle of the present invention should be included in the protection scope of the claims of the present invention.
Claims (10)
1. A method for analyzing drought and waterlogging turning incident characteristic quantity is characterized by comprising the following steps:
s1: determining a research area, and collecting daily scale precipitation data in the area;
s2: calculating a standardized weighted average rainfall SWAP index based on rainfall data of the drainage basin, identifying a drought and waterlogging rush-turn process of the drainage basin, and extracting intensity and duration characteristics of drought and waterlogging rush-turn events;
s3: adopting ADF to respectively check whether the time series of the intensity characteristic and the duration characteristic have consistency; when the time sequence consistency is established, fitting the intensity characteristic and the duration characteristic by adopting a probability distribution model with constant parameters, and obtaining an intensity edge distribution function and a duration edge distribution function based on the evaluation of the SBC criterion and the AIC criterion;
when the time sequence consistency is not established, fitting the intensity characteristics and the duration characteristics by adopting a generalized accumulation model respectively to obtain an intensity edge distribution function and a duration edge distribution function;
s4: adopting a Copula function likelihood ratio method to check whether a consistency assumption between the intensity edge distribution function and the duration edge distribution function is established or not;
when the consistency is assumed to be set, simulating by adopting a Copula joint distribution function with constant parameters to obtain a joint distribution function of duration and intensity;
when the consistency assumption is not satisfied, simulating by adopting a Copula joint distribution function with time-varying parameters to obtain a joint distribution function with duration and intensity;
s5: and calculating the combined recurrence period of duration and intensity or 'and' based on the combined distribution function of duration and intensity, and calculating a combined design value of duration and intensity by adopting a most possible combination method, and taking the design value as a risk threshold value for disaster prevention and control of drought and waterlogging emergency.
2. The method for analyzing the characteristic quantity of the drought/flood sudden-turn event according to claim 1, wherein the calculation process of the normalized weighted average precipitation SWAP index in the step S2 is as follows:
the SWAP index is an index on the premise that the drought and waterlogging state of the day is influenced by the early drought and waterlogging state and the rainfall of the day, the calculation method is to assume that the weighted average rainfall sequence of the same day obeys Gamma distribution, fit the weighted average rainfall sequence through the Gamma distribution, and finally perform standard normalization treatment on the cumulative probability distribution, wherein the calculation formula of the weighted average rainfall sequence is as follows:
w n ＝(1-a)a n-1
of the WAP i Representing a weighted mean precipitation sequence, P n Denotes the precipitation on day n, w n Is represented by P n A is a parameter of the weight decay with time, a n The weight of the nth day is a decay parameter along with time, N is the number of days of influence of early precipitation, and a and N are selected according to the local soil-water system in the research area, wherein a =0.9 and N =45 are selected.
3. The method for analyzing the drought and flood rush-turn incident characteristic quantity according to the claim 2, wherein in the step S2, the drought and flood rush-turn process of the drainage basin is identified and extracted by adopting a run-length theory based on the calculated SWAP index, and the specific steps are as follows:
setting index critical value X for drought and waterlogging index time sequence 0 、X 1 And a time length threshold value D 0 、D 1 And D 2 When the index is continuously lower than X 0 Time of not less than D 0 When drought occurs, it is considered to occur; when the index is continuously higher than X 1 Time of not less than D 1 When the flood occurs, the flood is considered to occur; when the time interval between drought and flood is not more than D 2 In time, the drought and waterlogging sudden turning event is considered to be formed; according to the drought-waterlogging index division standard, selecting a numerical value corresponding to the drought level upper bound and the waterlogging level lower bound as X 0 And X 1 I.e. X 0 = -1 and X 1 =1; determining a recognition duration threshold D based on a day index 0 ＝D 1 Time length threshold D of rush hour between drought state and flood state for =10 days 2 And (5) acquiring duration and intensity of each drought and waterlogging rush transfer event according to the extracted drought and waterlogging rush transfer events in 8 days.
4. The method for analyzing the drought and flood racing event characteristic quantity according to claim 1, wherein in step S3, when the time series consistency is established, the concrete process of fitting the intensity edge distribution function and the duration edge distribution function is as follows;
selecting 7 probability distribution models of exponential, gamma and three-parameter generalized gamma, normal, lognormal, logistic and Weibull as alternatives, respectively fitting intensity characteristics and duration characteristics through the 7 alternative models, and evaluating the goodness of fit of the 7 alternative models through an SBC (SBC) criterion;
the calculation formula of SBC is as follows:
SBC＝GD+log(ξ)×df＝-2log L(θ)+log(ξ)×df
in the formula, GD is global fitting deviation, df is the number of degrees of freedom, theta and xi respectively represent the estimated value and the number of parameters of probability distribution parameters, and L (-) represents a log-likelihood function value;
meanwhile, the goodness of fit of the 7 alternative models is evaluated through an AIC criterion;
the calculation formula for AIC is as follows:
AIC＝2ξ-2log L(θ)
and finally, selecting a distribution model which enables the calculated values of the AIC and the SBC to be minimum simultaneously from all the tested candidate models as a relatively optimal model, and using the relatively optimal model as a strength edge distribution function or a duration edge distribution function.
5. The method as claimed in claim 4, wherein the defined interval of the candidate normal and lognormal probability distribution models for the variables is [ - ∞, + ∞ ], and the intensity and duration features to be tested are both greater than 0, so that the defined interval of the 2 candidate models is modified to a positive interval, and the modified cumulative probability is expressed by the following formula:
in the formula, f (x | θ) is a probability density function, and θ is a probability distribution parameter.
6. The method for analyzing the drought and flood racing event characteristic quantity according to claim 4, wherein in step S3, when the time series consistency is not established, the specific process of fitting the intensity edge distribution function and the duration edge distribution function is as follows;
fitting the edge distribution by using a generalized additive model in which the response variable Y is assumed
Independent observed value
Obeying probability distribution f (y) i |Θ i )；
Wherein, theta i ＝(μ i ,σ i ,v i ,τ i ) The vector is composed of probability distribution parameters, mu and sigma are position and scale parameters, and v and tau are shape parameters, and represent skewness and kurtosis of distribution respectively;
establishing the relation between the probability distribution parameters and the interpretation variables thereof through a connection function, wherein the expression is as follows:
namely:
in the formula, g k (. -) represents a join function; theta k Mu, sigma, v, tau are lengths of
A column vector of (a);
X k is a known matrix, corresponds
Interpretation variable values for the time periods; k represents the sequence number of the sequence to be analyzed;
7. The method for analyzing the drought and flood racing event characteristic quantity according to claim 6, wherein in step S4, a Copula function likelihood ratio method is adopted to check whether the consistency assumption between the intensity edge distribution function and the duration edge distribution function is true or not as follows;
the basic assumption of the Copula function likelihood ratio method is that the type of the Copula function is fixed and is unchanged, and the parameter theta is detected c Analyzing the variable points of the multivariate hydrological sequence dependent structural strength;
knowing y 1 ，y 2 ，...，y d Is a time series of d-dimensional random variable observations, where y i ＝(y 1，i ，y 2，i ，…，y n，i ) I =1, 2.. D, the population where the observations lie obeys a probability distribution F (y) 1 |θ c，1 )，F(y 2 |θ c，2 )，...，F(y n |θ c，n )，θ c，i The Copula function parameters corresponding to different time periods, in the Copula function likelihood ratio test, the original assumption is that the multivariate dependent structure has no variable point, and the expression formula is as follows:
H 0 ：θ c，1 ＝θ c，2 ＝…＝θ c，n ＝η 0
alternative hypothesis H 1 Comprises the following steps: presence of k * And is provided with
So that theta c，1 ＝θ c，2 ＝…＝θ c.k* ＝η 1 And theta c，k* ＝…＝θ c，n ＝η 2 While eta 1 ≠η 2 (ii) a Constructing likelihood ratio Lambda based on Copula function k When Λ k Below a given threshold, rejecting the original hypothesis H 0 And considering that the multivariate dependent structure has a variable point k * ；
In the formula, L n (. H) represents a likelihood function for the entire sample sequence; l is k (. A) and
respectively before and after the check point kBut the function; u. u i Representing the edge distribution probability of the bivariate sequence; c (-) represents the Copula probability density function;
and
are Copula function parameters η, respectively 0 、η 1 And η 2 Maximum likelihood estimate of Λ, rewrite Λ k In logarithmic form, we obtain:
further, the statistical test quantity in the form of a construction log is as follows:
when Z is n Above a threshold corresponding to a given level of significance, reject the original hypothesis H 0 Likelihood ratio statistics
Obey an asymptotic distribution:
where z represents a given threshold, and as z approaches infinity, h = l = [ ln (n)] 3/2 ；
h. l and O represent parameters of asymptotic distribution; p is the number of parameters of the Copula function where there is one change point; will be provided with
Substituting the formula into the formula, if the calculated P value is less than the threshold value, the P value is significant at the threshold valueAnd rejecting the original hypothesis under the consistency level, wherein the consistency hypothesis does not hold, and otherwise, the consistency hypothesis holds.
8. The method for analyzing the characteristic quantity of the drought, flood and sudden turn events as claimed in claim 7, wherein in step S4, the joint distribution function of duration and intensity established according to Copula function is expressed as:
P(D，S)＝C(u D ，u S |θ c )
in the formula, D and S respectively represent intensity and duration; u. of D ，u S Edge distributions representing intensity and duration;
in order to estimate a time-varying parameter of Copula, introducing an explanation variable for describing a functional relation between the Copula time-varying parameter and the explanation variable;
g(θ c )＝β 0 +Xβ 1
wherein g (θ) c ) When the copula parameter is theta c A connection function of (a);
θ c ＝(θ c，1 ，θ c，2 ，...，θ c，n ) T is a column vector composed of Copula time-varying parameters;
X＝(X 1 ，X 2 ，...，X n ) T n observations representing an explanatory variable; beta is a 0 And β represents the parameter vector to be solved;
according to the consistency assumption of the relation between duration and intensity, giving out the expression of the connection function in different cases:
(1) When the two-variable relationship is consistent, the Copula parameter is a constant value which does not change with time, and the connection function is expressed as:
g(θ c )＝constant
(2) When the consistency assumption is broken, the Copula parameter is a function of time, thus introducing the time of occurrence of the drought event as an explanatory variable, the join function is rewritten as:
g(θ c )＝β 0 +time·β 1
further, a marginal function inference method is adopted, and log likelihood is maximizedFunction, solving Copula parameter θ c，t The formula is as follows:
in the formula u D，i Representing an intensity edge distribution for the ith period; u. of S，i An edge distribution representing duration of the ith time period; d i Representing the intensity of the ith time period; s i Represents the duration of the ith time period; theta D A parameter representing intensity edge distribution; theta.theta. S A parameter representing a temporal edge distribution; theta D，i A parameter representing intensity edge distribution for the ith time period; theta.theta. S，i A parameter representing a duration edge distribution of the ith time period; c (-) is a probability density function of Copula; f. of D (. And f) S () probability density functions representing duration, intensity edge distributions, respectively; in the marginal function inference method, because the maximum values of the second term and the third term on the right side of the equal sign of the formula are obtained by fitting the edge distribution of duration and intensity, the optimal beta can be obtained by maximizing the first term on the right side of the equal sign 0 And beta 1 And is beta 0 And beta 1 Substituting the estimated value of the Copula function into the estimated value of the Copula function;
in the formula, time i Indicating the occurrence time of the ith period.
9. The method for analyzing the characteristic quantity of the drought or flood sudden turn event according to the claim 8, wherein the step S5, the process of calculating the duration and intensity or the combined recurrence period of "and" is as follows: univariate characteristic P (x is more than or equal to x) of drought and waterlogging sharp turn * ) The recurrence period of (c) is calculated as follows:
wherein T is the recurrence period of drought and waterlogging sudden turn event, F (x) * ) Is x * An event probability function; the interval is the average time interval of all drought and flood turning events, and the unit is year, namely the time length from the beginning of a certain drought and flood turning event to the beginning of the next drought and flood turning event; if the duration is N year Annual sampling to obtain N alter The interval calculation formula of the scene sharp turn event is as follows:
the calculation formula of the reappearance period under different visual angles is as follows;
under univariate viewing angle:
in the formula, T S 、T D Duration, intensity, respectively, of the reconstruction period under univariate viewing angle, F S (s) and F D (d) Respectively representing the probability of occurrence of duration and intensity;
under a bivariate viewing angle:
in the formula (I), the compound is shown in the specification,
representing the "or" joint recurrence period under both variable perspectives of duration and intensity,
denotes the "and" joint recurrence period, F, at both time and intensity variables S，D (s, d) is the joint occurrence probability of duration and intensity;
the duration and severity "or" and "combined recurrence period is derived based on the various calculations described above.
10. The method for analyzing the characteristic quantity of the drought/flood sudden turn event according to the claim 9, wherein in the step S5, the most probable combination method is adopted to calculate the joint design value of duration and intensity as follows:
the combined design value of the drought and flood sudden-turn events refers to a multivariable combination (u, v) when the sudden-turn combined density f (u, v) composed of sudden-turn characteristics reaches a maximum value at a certain design recurrence period level ml ，v ml ) I.e., the combination is most likely to occur, and therefore, solved by constructing the following equation:
[u ml ，v ml ]＝argmax f(u，v)
f(u，v)＝c(u，v)·f(u)·f(v)
wherein f (u, v) is a joint probability density function of double variables of sharp turn, c (u, v) is a probability density function corresponding to joint Copula distribution, f (u) and f (v) are edge density functions of variables, [ u, v ] ml ，v ml ]Combining the joint design values of the multiple variables;
the calculation steps for obtaining the joint design value by adopting the most probable combination method are as follows:
(1) Simulating a plurality of variable combinations (u) by adopting a Monte Carlo method based on a Copula joint distribution function ml ，v ml ) The number of the groups is 100000;
(2) According to the calculation formula of the ' OR ' and ' recurrence period, the joint recurrence period corresponding to the bivariate combination of the last step is calculated;
(3) Root of Chinese scholar treeAccording to the determined level of the recurrence period, the screening meets the precision requirement (e = 10) -4 ) Wherein the variable combination with the maximum joint probability density value is the joint design value.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210810680.8A CN115168808A (en) | 2022-07-11 | 2022-07-11 | Method for analyzing drought and waterlogging sudden turning event characteristic quantity |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210810680.8A CN115168808A (en) | 2022-07-11 | 2022-07-11 | Method for analyzing drought and waterlogging sudden turning event characteristic quantity |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN115168808A true CN115168808A (en) | 2022-10-11 |
Family
ID=83492762
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210810680.8A Pending CN115168808A (en) | 2022-07-11 | 2022-07-11 | Method for analyzing drought and waterlogging sudden turning event characteristic quantity |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN115168808A (en) |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115329610A (en) * | 2022-10-17 | 2022-11-11 | 中山大学 | Method, device and equipment for identifying drought and waterlogging emergency turn based on soil moisture |
| CN116502891A (en) * | 2023-04-28 | 2023-07-28 | 西安理工大学 | Determination method of snow-drought dynamic risk |
| CN116663392A (en) * | 2023-04-23 | 2023-08-29 | 中国水利水电科学研究院 | Method, device, equipment and medium for simulating two-dimensional composite event of hydrological disaster |
-
2022
- 2022-07-11 CN CN202210810680.8A patent/CN115168808A/en active Pending
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115329610A (en) * | 2022-10-17 | 2022-11-11 | 中山大学 | Method, device and equipment for identifying drought and waterlogging emergency turn based on soil moisture |
| CN116663392A (en) * | 2023-04-23 | 2023-08-29 | 中国水利水电科学研究院 | Method, device, equipment and medium for simulating two-dimensional composite event of hydrological disaster |
| CN116663392B (en) * | 2023-04-23 | 2024-02-20 | 中国水利水电科学研究院 | Method, device, equipment and medium for simulating two-dimensional composite event of hydrological disaster |
| CN116502891A (en) * | 2023-04-28 | 2023-07-28 | 西安理工大学 | Determination method of snow-drought dynamic risk |
| CN116502891B (en) * | 2023-04-28 | 2024-03-29 | 西安理工大学 | Determination method of snow-drought dynamic risk |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN115168808A (en) | Method for analyzing drought and waterlogging sudden turning event characteristic quantity | |
| Gustafson | Quantifying landscape spatial pattern: what is the state of the art? | |
| Chiew et al. | Assessing the adequacy of catchment streamflow yield estimates | |
| CN106250905A (en) | A kind of real time energy consumption method for detecting abnormality of combination colleges and universities building structure feature | |
| CN111832607B (en) | Bridge disease real-time detection method based on model pruning | |
| Cappelaere et al. | Hydrologic process simulation of a semiarid, endoreic catchment in Sahelian West Niger. 2. Model calibration and uncertainty characterization | |
| CN109325888A (en) | A kind of students ' behavior prediction technique based on artificial neural network | |
| CN116597350A (en) | Flotation process fault early warning method based on BiLSTM predictive deviation | |
| CN115100819A (en) | Landslide hazard early warning method and device based on big data analysis and electronic equipment | |
| CN110738208A (en) | efficient scale-normalized target detection training method | |
| CN116760017A (en) | Prediction method for photovoltaic power generation | |
| CN115879594A (en) | Urban settlement population distribution trend prediction method based on geographic detector | |
| CN113487069B (en) | Regional flood disaster risk assessment method based on GRACE daily degradation scale and novel DWSDI index | |
| CN109544410A (en) | Cell source of houses value parameter estimation method and device | |
| CN111625525B (en) | Environment data repairing/filling method and system | |
| CN108876210A (en) | A kind of recognition methods, system and the device of land use and change causal structure | |
| CN114595425A (en) | Method for diagnosing and analyzing inconsistent mutation points of runoff relation of drainage basin rainfall | |
| CN114595764A (en) | Method and system for acquiring influence degree of urban factors on inland inundation disaster loss | |
| CN113807587A (en) | Integral early warning method and system based on multi-ladder-core deep neural network model | |
| CN110674471A (en) | Debris flow easiness prediction method based on GIS (geographic information System) and Logistic regression model | |
| Fan et al. | The impact of physiographic factors upon the probability of slides occurrence: a case study from the Kaoping River Basin, Taiwan | |
| CN112331342A (en) | Disease risk grade evaluation method based on gridding covariate factors | |
| CN110083930A (en) | The construction method and device of shale weathering index | |
| CN110751398A (en) | Regional ecological quality evaluation method and device | |
| CN117421985B (en) | Plant diffusion capacity assessment and early warning method |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination |