CN116127247B - Probability risk analysis and calculation method for coupling multiple seismic source models - Google Patents

Probability risk analysis and calculation method for coupling multiple seismic source models Download PDF

Info

Publication number
CN116127247B
CN116127247B CN202310111806.7A CN202310111806A CN116127247B CN 116127247 B CN116127247 B CN 116127247B CN 202310111806 A CN202310111806 A CN 202310111806A CN 116127247 B CN116127247 B CN 116127247B
Authority
CN
China
Prior art keywords
seismic
formula
earthquake
source
fault
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.)
Active
Application number
CN202310111806.7A
Other languages
Chinese (zh)
Other versions
CN116127247A (en
Inventor
潘华
李昌珑
程江
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
INSTITUTE OF GEOPHYSICS CHINA EARTHQUAKE ADMINISTRATION
Original Assignee
INSTITUTE OF GEOPHYSICS CHINA EARTHQUAKE ADMINISTRATION
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by INSTITUTE OF GEOPHYSICS CHINA EARTHQUAKE ADMINISTRATION filed Critical INSTITUTE OF GEOPHYSICS CHINA EARTHQUAKE ADMINISTRATION
Priority to CN202310111806.7A priority Critical patent/CN116127247B/en
Publication of CN116127247A publication Critical patent/CN116127247A/en
Application granted granted Critical
Publication of CN116127247B publication Critical patent/CN116127247B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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
    • 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)
  • Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • General Physics & Mathematics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Algebra (AREA)
  • Computational Mathematics (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Mathematical Analysis (AREA)
  • Testing Of Devices, Machine Parts, Or Other Structures Thereof (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

The invention relates to the technical field of earthquake engineering, and discloses a probability risk analysis and calculation method for coupling multiple earthquake source models, which comprises a fault source risk calculation method and a surface source risk calculation method, wherein the surface source risk calculation method specifically differentiates a potential earthquake source area into sub-sources and distributes earthquake activity parameters; the fault source risk calculation method is to expand the fault trace into a fault plane, gridding the fault plane and then obtain the earthquake fracture plane according to the magnitude-area empirical formula. And finally, coupling seismic risk results of different source models by adopting a full probability calculation method. The invention is coupled with the probability risk analysis and calculation method of various seismic source models, and can describe the seismic vibration influence of different potential seismic source areas around the field point more finely.

Description

Probability risk analysis and calculation method for coupling multiple seismic source models
Technical Field
The invention relates to the technical field of seismic engineering, in particular to a probability risk analysis and calculation method for coupling multiple seismic source models.
Background
The latest fifth generation earthquake motion parameter demarcation diagram is taken as a potential earthquake focus area model, an earthquake activity parameter model and matched earthquake motion attenuation relation parameters as basic input in China, and the probability earthquake risk analysis is carried out by taking the earthquake activity as a poisson model, so that the method is a main means for carrying out earthquake risk analysis, engineering field area safety evaluation and earthquake insurance model development. The fifth generation diagram adopts a three-level potential seismic source region model formed by a seismic statistics region, a background potential seismic source region and a construction potential seismic source region, and builds a corresponding seismic activity model. The potential source area is a collection in the macroscopic earthquake of the earthquake, the earthquake is regarded as a point instead of an earthquake fracture surface, the method is completely feasible for small and medium earthquakes or far-field large earthquakes, but the fracture surface scale of the earthquake with larger magnitude is larger, and the influence of the large earthquake fracture surface scale on the earthquake risk can be reflected by more reasonable consideration of the earthquake as the fracture surface when the earthquake risk is calculated. At present, in the domestic work of earthquake division, earthquake cell division, earthquake safety evaluation, earthquake insurance and the like, a single surface source model is mainly adopted, and the corresponding earthquake risk calculation is carried out by considering the earthquake as a point.
The method mainly uses a program-CPSHA 90 matched with a five-generation diagram in China, can map a potential seismic source region model and a seismic activity model, and is combined with a seismic motion prediction model to calculate the seismic risk, and the seismic risk result is that comprehensive products such as historical seismic data, geophysical data, seismic data and the like are fully utilized, so that the method has scientificity and advancement. At present, the domestic earthquake risk analysis program calculates the differential of a potential earthquake focus area into a point source, which is applicable to small and medium earthquakes, but is not reasonable for large earthquakes, especially near-field large earthquakes. In software such as openquate, the earthquake event is considered as a fracture surface instead of a point in the earthquake risk, and the assumption that the earthquake fracture surface has the position information of the spatial parameters such as trend, tendency, depth and the like instead of a single point is more accurate in expressing the earthquake risk. The domestic seismic risk research and application also develop from a single surface source model to a multi-model source direction, so that foreign seismic risk analysis software cannot be directly applied, and the foreign seismic risk analysis software or programs have certain limitations, such as not matching with a five-generation graph model and being inconvenient to derive the intermediate process quantity of calculation. Therefore, developing a set of earthquake risk algorithm which meets the requirements of domestic medium and small earthquakes on a regional map surface source model and simultaneously adopts multi-model coupling of a fault source model for large earthquakes is particularly important, and the point source and the fracture surface model can be flexibly selected to carry out risk calculation, so that the method has great help for more accurately describing the risk. In the next generation of earthquake zone diagrams, a potential earthquake source zone dividing technical route adopts a multi-type potential source mode, so that a dangerous algorithm for coupling point sources and a fracture surface model has certain reference significance for development of the zone diagrams.
The software for performing earthquake risk calculation at home is a program matched with a fifth generation earthquake motion parameter zoning chart, and the program can only calculate a model with a potential earthquake source zone as a surface source and consider the earthquake as a point. With the accumulation and development of fault data at present, a single surface source potential source region model cannot meet the requirements of the facets of earthquake activity and danger and is large in the cracking area of a large earthquake, and the large earthquake is regarded as a point to be calculated in the calculation process, so that a probability danger analysis calculation method for coupling multiple source models is needed.
Disclosure of Invention
The invention aims to provide a probability risk analysis and calculation method for coupling multiple seismic source models, which can describe the seismic influence of different potential seismic source areas around a field point more precisely.
The invention is realized in the following way:
the invention provides a probability risk analysis and calculation method for coupling multiple seismic source models, which comprises a fault source risk calculation method and a surface source risk calculation method, wherein the surface source risk calculation method is specifically executed according to the following steps:
S 1 differentiating the potential seismic source area into sub-sources, distributing seismic activity parameters, and calculating probability seismic risk by using a full probability formula;
S 2 setting the seismic event in the surface source of the seismic statistic region as a random event, wherein the seismic magnitude m and the distance r from a field point are random variables, and are related, the joint probability density function is f (m, r), nz seismic statistic regions are divided in the region altogether, and the probability that the field point causes A to be more than or equal to a is shown as formula (1) -formula (11) when the secondary seismic magnitude m and the distance r which randomly occur in the kth seismic statistic region are seismic;
wherein P (A is greater than or equal to a|m, r) is random uncertainty caused by the variance of the decay relation fit;
wherein f (m) k As for the magnitude distribution function in the kth seismic statistics area, f (r|m) k is the seismic space distribution function, and Nks potential seismic source areas are assumed to be shared in the kth seismic statistics area and are discretized, and the magnitude is divided into Nm grades; magnitude m j Representative ofThen formula (2) becomes as formula (3);
wherein s is ki Representing the ith potential source zone within the kth seismic statistics zone: as As ki Representing an ith potential source zone area within the kth seismic statistics zone; (s) ki |m j ) Representing the ith potential source zone m within the kth seismic statistics zone j The occurrence probability of the earthquake; let gn=a be the magnitude frequency relationship in the kth seismic statistics area k -b k m and n are the number of earthquakes greater than or equal to the magnitude m, a k ，b k Coefficients for the G-R relationship; as shown in formula (4);
wherein m is 0 For the lower limit of the magnitude, m uzk Is the upper limit of the magnitude of the kth seismic statistics area, beta k ＝b k ln10, formula (5);
p(s) ki ∣m j ) Abbreviated as f ki,mj Then formula (6) is obtainable by formulas (3) - (5);
let P (A is greater than or equal to a) k ＝p k And assuming that the kth seismic statistic region is m in unit time period 0 The number of times of occurrence of earthquakes above the level meets the average value of v mk0 According to the random selection characteristic of the Poisson distribution, the occurrence times of the earthquake motion A is larger than or equal to the a event at the siteAlso consistent with an average incidence of p k v m0k Poisson distribution of (a) as in formula (7);
the earthquake motion overrun probability of all the earthquakes to the field points in the unit time period of the kth earthquake statistical region can be deduced according to the probability theory, and the earthquake motion overrun probability is shown as a formula (8);
the influence of the earthquake in all the earthquake statistical areas around the field point on the field point is synthesized according to the full probability theorem, as shown in the formula (9);
substituting formula (8) into formula (9) to obtain formula (10);
substituting the formula (6) into the formula (10) to obtain the overrun probability of the site earthquake motion parameter A exceeding the given value a as shown in the formula (11), wherein the result calculated by the formula (11) is the earthquake motion influence of all potential earthquake source areas around the site on the site;
S 3 the non-point source risk calculation method is specifically implemented according to the following steps of;
S 3.1 importing a potential seismic source region model and seismic activity parameters and coordinate data of a required site;
S 3.2 reading longitude and latitude coordinates of field points, counting the number of the field points, presetting a PGA threshold value in an overrun probability curve, and presetting the overrun probability curveSetting the PGA threshold value range between 1 and 2000gal, wherein the attenuation relation is divided into different vibration levels and attenuation directions, and the attenuation relation model is shown as a formula (12);
lgY＝A+BW+Clg(R+ EM De) (12)
Wherein Y represents a seismic parameter, a seismic acceleration or velocity, M represents a surface wave magnitude, R represents a midrange, A, B, C, D, and E are regression coefficients;
S 3.3 gridding potential seismic source areas and distributing occurrence rates;
S 3.4 calculating the mean earthquake motion and annual overrun probability.
Further, in step S 3.4 Firstly, calculating the distance between grid points, namely the ground vibration attenuation distance; converting the coordinate system to a local distance coordinate system which takes a point source as an origin, takes the attenuation long axis direction as an x 'axis and takes the attenuation short axis direction as a y' axis;
respectively calculating the parameters of the attenuation long axis and the attenuation short axis along the connection line of the field point and the grid point to obtain the possible maximum value and the possible minimum value of the earthquake motion; calculating the average value of the influence of grid points on the site vibration according to the elliptical attenuation relation by adopting a dichotomy;
taking the mean earthquake motion parameter as a normal distribution mean value, constructing normal distribution by taking the standard deviation in the attenuation relation parameter as the normal distribution standard deviation, and respectively calculating the probability of a group of thresholds of the preset earthquake motion parameter, wherein the probability is shown as a formula (13);
wherein P is t ' A is equal to or greater than a, and is the t-year overrun probability after peak acceleration correction; f (f) 1 (z) a probability density function of normal distribution with a mean value of 0 and a standard deviation of sigma as shown in formula (14);
the relation between the annual overrun probability P1 (Y > Y) and the overrun probability Pt (Y > Y) of t years exists in the formulas (12) and (13), the annual overrun probability can be converted into an arbitrary annual overrun probability which is generally 50 years or 100 years, and then the earthquake motion parameters under the specified probability, such as 2%, 10% and 53%, are obtained according to the interpolation method. As in formula (15) -formula (16);
P 1 (Y＞y)＝1-[1-P t (Y＞y)] 1/t (15)
P t (Y＞y)＝1-[1-P 1 (Y＞y)] t (16)
further, the fault source risk calculation method generates a seismic fracture surface according to the calculated magnitude for each discretized sub-source, and the seismic fracture surface calculation method specifically comprises the following steps:
S 4.1 : firstly, reading longitude and latitude coordinates of field points, counting the number of the field points, and exceeding a pga threshold value in a probability curve and a seismic attenuation relation parameter of a fault source; the attenuation relation mode is shown as a formula (17);
lgY(M，R)＝A+BM+Clg(R+De EM ) (17)
S 4.2 : computed region tomographic data comprising the following parameters: longitude and latitude coordinates of each turning point of the fault trace, upper and lower depths of a pregnant earthquake region, fault sliding angle, starting earthquake level, upper earthquake level limit, inclination angle, sliding speed, b value and earthquake level gear step length;
S 4.3 calculating fault activity parameters;
S 4.4 constructing a gridding for the fault plane;
S 4.5 calculating the breaking area and the grid points;
S 4.6 calculating fault fracture face distance R rup 。
Further, step S 4.3 The method is as follows;
b value adopts the value b of the seismic zone where the fault is located in the relation of magnitude-frequency of the fault, and G-R relation a value on the fault source is calculated by using the relation of the fault sliding rate and the seismic incidence;
and then, according to the relation of the magnitude and the frequency, the annual incidence rate is distributed on each magnitude grade, and the median value of the magnitude grade is used as the calculated magnitude grade of the magnitude grade.
Further, it is characterized by S 4.5 Firstly, converting a surface wave magnitude Ms into a moment magnitude Mw, then converting the magnitude to be calculated into a fracture surface size according to a magnitude-fracture area empirical formula, selecting different empirical formulas for different fault types when the fracture surface area is calculated, and then obtaining the length L (M) and the width W (M) of the fracture surface according to the aspect ratio of the seismic fracture surface; as in formula (18) -formula (19);
MW=0.86.Ms+0.5gMs < 7.0 formula (18)
Mw=1.28.times.MS-2.42 Ms. Gtoreq.7.0 formula (19)
Then respectively calculating the grid number required by the earthquake fracture surface in the trend direction and the trend direction, wherein the spreading direction of the earthquake fracture surface is consistent with the spreading direction of the faults, if the earthquake fracture surface cannot be rounded, rounding the grid number, the rounding purpose is to facilitate the fracture surface to slide on the fault grid surface,and->Respectively representing the number of grid points required by the length and width directions of the fracture surface;and->Representing the number of times the fracture surface can slide along the fault strike and trend directions, respectively, wherein delta is the mesh size, as in formulas (20) - (23);
if the length or width direction of the earthquake fracture surface is larger than the trend length or width of the fracture surface, taking the length or width of the fracture surface as the length or width of the fracture surface, and then recalculating the length or width according to the redetermined length or width under the condition of ensuring that the fracture area is unchanged; when the length and width directions of the fracture surface are larger than those of the fracture surface, the whole fracture surface is directly selected as the fracture surface.
Further, in step S 4.6 In, calculate fault fracture face distance R rup Comprises a point-by-point calculation method and a parsing calculation method;
further, the present invention provides a computer-readable storage medium having stored thereon a computer program which, when executed by a main controller, implements a method as described in any of the above.
Compared with the prior art, the invention has the beneficial effects that:
1. the invention uses the potential source area as the grid point, uses the grid point as the differential unit and uses the source distance as the attenuation distance, and uses the elliptical attenuation mode to calculate the earthquake vibration parameter; for a fracture surface model, dispersing a potential source into equidistant grid surfaces, floating the obtained seismic fracture surface on the potential source according to the magnitude and area relation and the aspect ratio, and obtaining the seismic vibration parameter in a linear attenuation mode by taking the fracture surface distance as an attenuation distance; and then obtaining a point source and a fracture surface source according to a full probability formula, and integrating the point source and the fracture surface source to obtain the multi-source model danger to the field point. The multi-source model risk is more reasonable than the method of adopting a single point source mode to describe the risk. The influence of the large earthquake fracture scale on the danger can be reflected by adopting an earthquake fracture surface model algorithm. The Fenwei seismic zone is selected as the application, wherein the potential seismic source zone of a fifth generation seismic parameter demarcation graph and a seismic activity model are used as inputs for 4.0-6.5 seismic class, the potential source effect on the site seismic vibration is calculated by using a point source model, and the fault source mode is adopted for the above-6.5 seismic class to calculate the fault effect on the site seismic vibration. The point source part inherits the characteristics of seismic regions in China in the potential source division and the seismic activity parameters, and the fracture surface model part applies more detailed fault data. The algorithm for coupling the point source and the fracture surface model can describe the earthquake motion influence of different potential earthquake focus areas around the field point more finely.
Drawings
In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings that are needed in the embodiments will be briefly described below, it being understood that the following drawings only illustrate some examples of the present invention and therefore should not be considered as limiting the scope, and other related drawings may be obtained according to these drawings without inventive effort for a person skilled in the art.
FIG. 1 is a flow chart of a fault source hazard method of the present invention;
FIG. 2 is a schematic diagram of a potential source area grid of the present invention;
FIG. 3 is a schematic view of the distance parameters of the present invention;
FIG. 4 is a coordinate system conversion schematic diagram of the present invention;
FIG. 5 is a schematic representation of a fracture surface of the present invention;
FIG. 6 is a schematic view of the burst interval of the present invention;
fig. 7 is a flowchart of a method of the present invention for a non-point source risk.
Detailed Description
For the purpose of making the objects, technical solutions and advantages of the embodiments of the present invention more apparent, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention, and it is apparent that the described embodiments are some embodiments of the present invention, but not all embodiments. All other embodiments, based on the embodiments of the invention, which are apparent to those of ordinary skill in the art without inventive faculty, are intended to be within the scope of the invention. Thus, the following detailed description of the embodiments of the invention, as presented in the figures, is not intended to limit the scope of the invention, as claimed, but is merely representative of selected embodiments of the invention. All other embodiments, based on the embodiments of the invention, which are apparent to those of ordinary skill in the art without inventive faculty, are intended to be within the scope of the invention.
Referring to fig. 1-7, a method for calculating probability risk analysis of coupling multiple seismic source models includes a fault source risk calculation method and a surface source risk calculation method, wherein the surface source risk calculation method is specifically implemented according to the following steps:
S 1 differentiating the potential seismic source area into sub-sources, distributing seismic activity parameters, and calculating probability seismic risk by using a full probability formula;
S 2 setting the seismic event in the surface source of the seismic statistic region as a random event, wherein the seismic magnitude m and the distance r from a field point are random variables, and are related, the joint probability density function is f (m, r), nz seismic statistic regions are divided in the region altogether, and the probability that the field point causes A to be more than or equal to a is shown as formula (1) -formula (11) when the secondary seismic magnitude m and the distance r which randomly occur in the kth seismic statistic region are seismic;
wherein P (A is greater than or equal to a|m, r) is random uncertainty caused by the variance of the decay relation fit;
wherein f (m) k As for the magnitude distribution function in the kth seismic statistics area, f (r|m) k is the seismic space distribution function, and Nks potential seismic source areas are assumed to be shared in the kth seismic statistics area and are discretized, and the magnitude is divided into Nm grades; magnitude m j Representative ofThen formula (2) becomes as formula (3);
wherein s is ki Representing the ith potential source zone within the kth seismic statistics zone: as As ki Representing an ith potential source zone area within the kth seismic statistics zone; (s) ki |m j ) Representing the ith potential source zone m within the kth seismic statistics zone j The occurrence probability of the earthquake; let gn=a be the magnitude frequency relationship in the kth seismic statistics area k -b k m and n are the number of earthquakes greater than or equal to the magnitude m, a k ，b k Coefficients for the G-R relationship; as shown in formula (4);
wherein m is 0 For the lower limit of the magnitude, m uzk Is the upper limit of the magnitude of the kth seismic statistics area, beta k ＝b k ln10, formula (5);
p(s) ki ∣m j ) Abbreviated as f ki,mj Then formula (6) is obtainable by formulas (3) - (5);
let P (A is greater than or equal to a) k ＝p k And assuming that the kth seismic statistic region is m in unit time period 0 The number of times of occurrence of earthquakes above the level meets the average value of v mk0 Is a poisson distribution, and the earthquake occurs at the site according to the random selection characteristic of the poisson distributionThe number of events of the movement A being more than or equal to a also accords with the average incidence rate p k v m0k Poisson distribution of (a) as in formula (7);
the earthquake motion overrun probability of all the earthquakes to the field points in the unit time period of the kth earthquake statistical region can be deduced according to the probability theory, and the earthquake motion overrun probability is shown as a formula (8);
the influence of the earthquake in all the earthquake statistical areas around the field point on the field point is synthesized according to the full probability theorem, as shown in the formula (9);
substituting formula (8) into formula (9) to obtain formula (10);
substituting the formula (6) into the formula (10) to obtain the overrun probability of the site earthquake motion parameter A exceeding the given value a as shown in the formula (11), wherein the result calculated by the formula (11) is the earthquake motion influence of all potential earthquake source areas around the site on the site;
S 3 the non-point source risk calculation method is specifically implemented according to the following steps of;
S 3.1 importing a potential seismic source region model and seismic activity parameters and coordinate data of a required site;
S 3.2 reading longitude and latitude coordinates of field points, counting the number of the field points, exceeding a preset PGA threshold value in a probability curve, and exceedingThe range of the preset PGA threshold value in the probability curve is between 1 and 2000gal, the attenuation relation is divided into different vibration levels and attenuation directions, and the attenuation relation model is shown as a formula (12);
lgY＝A+BW+Clg(R+De EM ) (12)
Wherein Y represents a seismic parameter, a seismic acceleration or velocity, M represents a surface wave magnitude, R represents a midrange, A, B, C, D, and E are regression coefficients; specifically as shown in Table 1;
TABLE 1 Point Source model attenuation relation parameter Format
S 3.3 Gridding potential source regions and assigning occurrences as in FIG. 2; meshing a potential seismic source area within a field point 200km range according to the size of 0.2 degrees by 0.2 degrees, and meshing the potential seismic source area by taking each grid point as a point source; distributing the occurrence rate of each earthquake level on the potential earthquake source area to each grid point according to the annual average occurrence rate, the earthquake level-frequency distribution relation and the spatial distribution function;
S 3.4 as shown in figure 3, calculating the average earthquake motion and annual overrun probability, and calculating the distance between a grid point and a field point, wherein the distance is the earthquake center distance, namely the earthquake motion attenuation distance, and the attenuation distance becomes the earthquake source distance if the earthquake source depth is considered; as shown in fig. 4, the coordinate system is converted into a local distance coordinate system with the origin of the point source, the x 'axis of the attenuation major axis direction and the y' axis of the attenuation minor axis direction.
In the present embodiment, in step S 3.4 Firstly, calculating the distance between grid points, namely the ground vibration attenuation distance; converting the coordinate system to a local distance coordinate system which takes a point source as an origin, takes the attenuation long axis direction as an x 'axis and takes the attenuation short axis direction as a y' axis;
and then respectively calculating attenuation major axis and attenuation minor axis parameters along the connection line of the field point and the grid point to obtain the possible maximum value and minimum value of the earthquake motion
Calculating the average value of the influence of grid points on the ground vibration of the field point according to the elliptical attenuation relation by adopting a dichotomy, and firstly calculating the length axis value of the ellipse according to the average value of the maximum value and the minimum value in the steps; and substituting the field point coordinates into an elliptic relational expression, if the result is larger than 1 and the field point falls outside the ellipse, if the result is smaller than 1 and the field point falls inside the ellipse, selecting the median value of the latest twice earthquake motion values according to the position of the ellipse, and continuing to calculate until the difference value of the latest twice earthquake motion values is smaller than a threshold value, and determining that the calculation effect is achieved.
The condition threshold value set in the dichotomy is not suitable to be too small, if the too small calculation is too large, even the program is possibly circulated, the threshold value is set to be 0.5gal or 1gal through repeated checking, so that the accuracy of the calculation result can be ensured, and the calculation speed can be improved. If the threshold is too large, the accuracy of the result of the overrun probability curve is lost, and if the threshold is too small, the calculation step is increased, and the calculation efficiency is reduced
The mean seismic vibration parameter is taken as a normal distribution mean value, the standard deviation in the attenuation relation parameter is the normal distribution standard deviation to construct normal distribution, the probability of a group of threshold values of the preset seismic vibration parameter is calculated respectively, the calculation of the part is also called seismic vibration parameter uncertainty correction, a front correction algorithm is adopted for standard deviation decomposition when seismic hazard result decomposition is carried out conveniently, and the general correction range is 3 times of standard deviation interval. Thus obtaining the annual overrun probability curve result of the non-point source to the field point, as shown in the formula (13);
wherein P is t ' A is equal to or greater than a, and is the t-year overrun probability after peak acceleration correction; f (f) 1 (z) a probability density function of normal distribution with a mean value of 0 and a standard deviation of sigma as shown in formula (14);
the relation between the annual overrun probability P1 (Y > Y) and the overrun probability Pt (Y > Y) of t years exists in the formulas (12) and (13), the annual overrun probability can be converted into an arbitrary annual overrun probability which is generally 50 years or 100 years, and then the earthquake motion parameters under the specified probability, such as 2%, 10% and 53%, are obtained according to the interpolation method. As in formula (15) -formula (16);
P 1 (Y＞y)＝1-[1-P t (Y＞y)] 1/t (15)
P t (Y＞y)＝1-[1-P 1 (Y＞y)] t (16)
in this embodiment, the fault source risk calculation method generates a seismic fracture surface for each sub-source after discretization according to the calculated magnitude, and the seismic fracture surface calculation method specifically includes the following steps:
S 4.1 : firstly, reading longitude and latitude coordinates of field points, and counting the number of the field points, a pga threshold value in an overrun probability curve and a seismic attenuation relation parameter of a fault source, wherein the table 2 is shown; the attenuation relation mode is shown as a formula (17);
lgY(M，R)＝A+BM+Clg(R+De EM ) (17)
TABLE 2 damping relationship parameter format for fracture surface model
A B C D E Standard deviation sigma
Ms≤6.5 1.812 0.555 1.935 0.956 0.462 0.377
Ms＞6.5 2.389 0.474 1.935 0.956 0.462 0.377
S 4.2 : computed region tomographic data comprising the following parameters: longitude and latitude coordinates of each turning point of the fault trace, upper and lower depths of a pregnant earthquake region, fault sliding angle, starting earthquake level, upper earthquake level limit, inclination angle, sliding speed, b value and earthquake level gear step length;
S 4.3 : calculating fault activity parameters;
S 4.4 constructing a gridding for the fault plane; firstly, the read fault data is a projection trace of a fault on the earth surface, which is defined as an intersection point between the fault surface and the earth surface, the trace of the fault needs to be expanded into a fault plane according to trend, dip angle and the like, and then subsequent calculation is carried out.
(1) And converting the longitude and latitude coordinate system into a distance coordinate system, selecting the minimum value of longitude and latitude according to the longitude and latitude coordinate of each fault trace, forming a reference point O (Lonmin, latmin), and then moving the O point downwards for a distance of between pregnancy areas/sin (dip). This ensures that the fault trace and the whole fault plane interval after expansion into planes fall entirely in the first quadrant for calculation
(2) And under a distance coordinate system, adopting a grid with the size of 2km x 2km to grid-divide the fault plane. Firstly, calculating the slope between connecting lines at two ends of a fault as a fault trace trend, translating the fault trace from top to bottom in a seismic inoculation depth along a direction perpendicular to the fault trace trend, and constructing a fault surface angle with an inclination angle equal to the inclination angle. It should be noted that for segmented fault traces, the segmentation cannot be performed separately to ungrid the traces of different segments are not on the same plane. And calculating the slope and the length of the sectional trace, supplementing the next section of the sectional trace with insufficient grid size, and finally selecting to round and round when the trend length and the trend length on the whole fault plane are insufficient in grid size, wherein the grid number of the obtained grid plane is an integer in the trend and trend aspects. For characteristic earthquake faults, one fault is an earthquake pregnancy range without grid; for complex faults, the fault is divided into different areas according to boundary parameters, and then the fault is gridded according to a simple fault method.
S 4.5 Calculating the breaking area and the grid points;
S 4.6 calculating fault fracture face distance R rup 。
In the present embodiment, step S 4.3 The method is as follows; b value adopts the value b of the seismic zone where the fault is located in the relation of magnitude-frequency of the fault, and G-R relation a value on the fault source is calculated by using the relation of the fault sliding rate and the seismic incidence;
and then, according to the relation of the magnitude and the frequency, the annual incidence rate is distributed on each magnitude grade, and the median value of the magnitude grade is used as the calculated magnitude grade of the magnitude grade.
In the present embodiment, in step S 4.5 Firstly converting a surface wave magnitude Ms into a moment magnitude Mw, then converting the magnitude to be calculated into a fracture surface size according to a magnitude-fracture area empirical formula, selecting different empirical formulas for different fault types when the fracture surface area is calculated, such as a data relation of a fault type shown in Table 3, and then obtaining the length L (M) and the width W (M) of the fracture surface according to the aspect ratio of the earthquake fracture surface; as in formula (18) -formula (19);
Mw=0.86×Ms+0.59Ms < 7.0 formula (18)
Mw=1.28.Ms-2.42 Ms.gtoreq.7.0 formula (19)
Table 3 data relationship of fault types
Then respectively calculating the grid number required by the earthquake fracture surface in the trend direction and the trend direction, wherein the spreading direction of the earthquake fracture surface is consistent with the spreading direction of the faults, if the earthquake fracture surface cannot be rounded, rounding the grid number, the rounding purpose is to facilitate the fracture surface to slide on the fault grid surface,and->Respectively representing the number of grid points required by the length and width directions of the fracture surface;and->Representing the number of times the fracture surface can slide along the fault strike and trend directions, respectively, wherein delta is the mesh size, as in formulas (20) - (23);
if the length or width direction of the seismic fracture surface is greater than the strike length or width of the fracture surface, the length or width of the fracture surface is taken as the length or width of the fracture surface, and the width or length is recalculated according to the redetermined length or width under the condition of keeping the fracture area unchanged; when the length and width directions of the fracture surface are larger than those of the fracture surface, the whole fracture surface is directly selected as the fracture surface. Wherein, the liquid crystal display device comprises a liquid crystal display device,
a) Representing fault trace expansion into fault plane grid-connected meshing;
b) Representing a single earthquake at a fault level;
c) Indicating that the earthquake is floating on the fault plane;
in the sliding seismic fracture surface, as shown in the figure b), generally selecting the corners of the fracture surface as initial positions, then sliding the seismic fracture surface in the trend direction and the trend direction of the fracture surface in a stepping mode respectively, and determining the total sliding number;
/>
in the present embodiment, in step S 4.6 In, calculate fault fracture face distance R rup Comprises a point-by-point calculation method and a parsing calculation method; the point-by-point calculation method comprises the steps of directly calculating the distance from each grid vertex coordinate to a field point, and then selecting the minimum value of the distance to be the fracture face distance. The calculation is convenient and quick.
The analysis calculation method includes two cases, namely, the shortest distance from a point to a limited plane on a space on a horizontal plane is divided into two cases, if the point is located in a section where the limited plane is projected on the horizontal plane along a normal vector, the shortest distance from the point to the limited plane is the distance between the point and the point projected on the plane along the normal vector of the limited plane, for example, site1 in fig. 5; secondly, if the field point is not in the projection area, the distance from the field point to the limited plane is the minimum value of the shortest distances from the field point to four side lines of the limited plane, namely, the shortest distances from the field point to four line segments are respectively calculated, and then the minimum value is taken, for example, site2 in fig. 6.
It should be noted that if the shortest distance from the field point to the seismic fracture surface is calculated on a segmented fault according to the method, if the seismic fracture surface is not a plane across multiple fault segments, i.e. the fracture surface is not a plane but a fracture surface, the fracture surface needs to be divided into different surfaces for recalculation.
In this embodiment, the present invention provides a computer-readable storage medium having stored thereon a computer program which, when executed by a host controller, implements a method as described in any of the above.
The above description is only of the preferred embodiments of the present invention and is not intended to limit the present invention, and various modifications and variations may be made to the present invention by those skilled in the art. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention should be included in the protection scope of the present invention.

Claims (9)

1. A probability risk analysis and calculation method for coupling multiple seismic source models is characterized in that: the method comprises a fault source risk calculation method and a non-point source risk calculation method, wherein the non-point source risk calculation method is specifically executed according to the following steps:
S 1 : differentiating the potential seismic source area into sub-sources, distributing seismic activity parameters, and calculating probability seismic risk by using a full probability formula;
S 2 : setting the seismic event in the surface source of the seismic statistics area as a random event, wherein the seismic magnitude m and the distance r from a field point are random variables, the seismic magnitude m and the distance r from the field point are related, the joint probability density function is f (m, r), nz seismic statistics areas are divided in the area altogether, and the probability that the A is more than or equal to a is shown as the formula (1) -formula (11) when the secondary seismic magnitude m and the distance r which randomly occur in the kth seismic statistics area are seismic;
P(A≥a) k ＝∫ m ∫ r f(m，r) k ·P(A≥a|m，r)d r d m (1)
Wherein P (A is larger than or equal to a|m, r) is random uncertainty caused by fitting variance of the attenuation relation;
P(A≥a) k ＝∫∫ m f(m) k ·f(r|m) k ·P(A≥a|m，r)d r d m (2)
Wherein f (m) k As for the magnitude distribution function in the kth seismic statistics area, f (r|m) k is the seismic space distribution function, and Nks potential seismic source areas are assumed to be shared in the kth seismic statistics area and are discretized, and the magnitude is divided into Nm grades; magnitude m j Representative ofThen formula (2) becomes as formula (3);
wherein s is ki Representing the ith potential source zone within the kth seismic statistics zone: as As ki Representing an ith potential source zone area within the kth seismic statistics zone; (s) ki |m j ) Representing the ith potential source zone m within the kth seismic statistics zone j The occurrence probability of the earthquake; let gn=a be the magnitude frequency relationship in the kth seismic statistics area k -b k m and n are the number of earthquakes greater than or equal to the magnitude m, a k ，b k Coefficients for the G-R relationship; as shown in formula (4);
wherein m is 0 For the lower limit of the magnitude, m uzk Is the upper limit of the magnitude of the kth seismic statistics area, beta k ＝b k ln10，
As shown in formula (5);
p(s) ki ∣m j ) Abbreviated as f ki,mj Then formula (6) is obtainable by formulas (3) - (5);
let P (A is greater than or equal to a) k ＝p k And assuming that the kth seismic statistic region is m in unit time period 0 The number of times of occurrence of earthquakes above the level meets the average value of v mk0 According to the random selection characteristic of the Poisson distribution, the number of times of occurrence of the earthquake motion A is larger than or equal to a event at the site also accords with the average occurrence rate p k v m0k Poisson distribution of (a) as in formula (7);
the earthquake motion overrun probability of all the earthquakes to the field points in the unit time period of the kth earthquake statistical region can be deduced according to the probability theory, and the earthquake motion overrun probability is shown as a formula (8);
the influence of the earthquake in all the earthquake statistical areas around the field point on the field point is synthesized according to the full probability theorem, as shown in the formula (9);
substituting formula (8) into formula (9) to obtain formula (10);
substituting the formula (6) into the formula (10) to obtain the overrun probability of the site earthquake motion parameter A exceeding the given value a as shown in the formula (11), wherein the result calculated by the formula (11) is the earthquake motion influence of all potential earthquake source areas around the site on the site;
S 3 : the non-point source risk calculation method is specifically executed according to the following steps of;
S 3.1 : importing a potential seismic source region model and seismic activity parameters and coordinate data of a required field point;
S 3.2 : reading longitude and latitude coordinates of field points, counting the number of the field points, exceeding a preset PGA threshold value in a probability curve, and dividing an attenuation relation into different vibration levels and attenuation directions, wherein the attenuation relation model is shown as a formula (12);
lgY＝A+BM+Clg(R+ EM De) (12)
Wherein Y represents a seismic parameter, a seismic acceleration or velocity, M represents a surface wave magnitude, R represents a midrange, A, B, C, D, and E are regression coefficients;
S 3.3 gridding potential seismic source areas and distributing occurrence rates;
S 3.4 calculating the mean earthquake motion and annual overrun probability.
2. The method of claim 1, wherein in step S 3.4 Firstly, calculating the distance between grid points, namely the ground vibration attenuation distance; converting the coordinate system to a local distance coordinate system which takes a point source as an origin, takes the attenuation long axis direction as an x 'axis and takes the attenuation short axis direction as a y' axis;
respectively calculating the parameters of the attenuation long axis and the attenuation short axis along the connection line of the field point and the grid point to obtain the possible maximum value and the possible minimum value of the earthquake motion; calculating the average value of the influence of grid points on the site vibration according to the elliptical attenuation relation by adopting a dichotomy;
taking the mean earthquake motion parameter as a normal distribution mean value, constructing normal distribution by taking the standard deviation in the attenuation relation parameter as the normal distribution standard deviation, and respectively calculating the probability of a group of thresholds of the preset earthquake motion parameter, wherein the probability is shown as a formula (13);
wherein P is t ' A is equal to or greater than a, and is the t-year overrun probability after peak acceleration correction; f (f) 1 (z) a probability density function of normal distribution with a mean value of 0 and a standard deviation of sigma as shown in formula (14);
3. the method for calculating the probability risk analysis of coupling multiple seismic source models according to claim 1, wherein the fault source risk calculation method generates a seismic fracture surface for each sub-source after discretization according to the calculated magnitude, and the method for calculating the seismic fracture surface is specifically executed by the following steps:
S 4.1 firstly, reading longitude and latitude coordinates of field points, and counting the number of the field points, a pga threshold value in an overrun probability curve and earthquake motion attenuation relation parameters of a fault source; the attenuation relation mode is shown as a formula (17);
lgY(M，R)＝A+BM+Clg(R+ EM De) (17) S 4.2 Calculating regional fault data comprising the following parameters: longitude and latitude coordinates of each turning point of the fault trace, upper and lower depths of a pregnant earthquake region, fault sliding angle, starting earthquake level, upper earthquake level limit, inclination angle, sliding speed, b value and earthquake level gear step length;
S 4.3 calculating fault activity parameters;
S 4.4 constructing a gridding for the fault plane;
S 4.5 calculating the breaking area and the grid points;
S 4.6 calculating fault fracture face distance R rup 。
4. The method of claim 3, wherein the step of S 4.3 The method is as follows;
b value adopts the value b of the seismic zone where the fault is located in the relation of magnitude-frequency of the fault, and G-R relation a value on the fault source is calculated by using the relation of the fault sliding rate and the seismic incidence;
and then, according to the relation of the magnitude and the frequency, the annual incidence rate is distributed on each magnitude grade, and the median value of the magnitude grade is used as the calculated magnitude grade of the magnitude grade.
5. The method of claim 1, wherein in step S 3.2 The preset PGA threshold value in the overrun probability curve is in the range of 1-2000 gal.
6. The method of claim 3, wherein S 4.5 Firstly, converting a surface wave magnitude Ms into a moment magnitude Mw, then converting the magnitude to be calculated into a fracture surface size according to a magnitude-fracture area empirical formula, selecting different empirical formulas for different fault types when the fracture surface area is calculated, and then obtaining the length L (M) and the width W (M) of the fracture surface according to the aspect ratio of the seismic fracture surface; as in formula (18) -formula (19);
Mw=0.86×Ms+0.59Ms < 7.0 formula (18)
Mw=1.28×Ms-2.42 Ms.gtoreq.7.0 formula (19).
7. The method of claim 6, wherein S 4.5 In the middle, the grid number required by the earthquake fracture surface in the trend direction and the trend direction is calculated respectively, the spreading direction of the earthquake fracture surface is consistent with the spreading direction of the fault,and->Respectively representing the number of grid points required by the length and width directions of the fracture surface;and->Representing the number of times the fracture surface can slide along the fault strike and trend directions, respectively, wherein delta is the mesh size, as in formulas (20) - (23);
8. the method of claim 6, wherein in step S 4.6 In, calculate fault fracture face distance R rup Including point-by-point calculations and analytical calculations.
9. A computer readable storage medium, on which a computer program is stored, characterized in that the program, when executed by a main controller, implements the method according to any of claims 1-8.
CN202310111806.7A 2023-02-14 2023-02-14 Probability risk analysis and calculation method for coupling multiple seismic source models Active CN116127247B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310111806.7A CN116127247B (en) 2023-02-14 2023-02-14 Probability risk analysis and calculation method for coupling multiple seismic source models

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310111806.7A CN116127247B (en) 2023-02-14 2023-02-14 Probability risk analysis and calculation method for coupling multiple seismic source models

Publications (2)

Publication Number Publication Date
CN116127247A CN116127247A (en) 2023-05-16
CN116127247B true CN116127247B (en) 2023-08-18

Family

ID=86306177

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310111806.7A Active CN116127247B (en) 2023-02-14 2023-02-14 Probability risk analysis and calculation method for coupling multiple seismic source models

Country Status (1)

Country Link
CN (1) CN116127247B (en)

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103955620A (en) * 2014-05-13 2014-07-30 中国地质大学（北京） Engineering site earthquake hazard analysis method considering effect of potential earthquake source orientations
CN109001805A (en) * 2018-05-24 2018-12-14 深圳大学 A kind of earthquake source inverting Uncertainty Analysis Method, storage medium and server
CN109375253A (en) * 2018-12-13 2019-02-22 中国地震局地球物理研究所 Ground motion parameter evaluation method based on whole seismic structure maximum credible earthquakes
CN110334482A (en) * 2019-05-25 2019-10-15 成都理工大学 Mud-rock flow Dynamic flexibility evaluation method after shake based on material resource activity intensity
CN113187556A (en) * 2021-05-06 2021-07-30 西南交通大学 Tunnel earthquake risk analysis method under fault dislocation based on total probability
CN113268852A (en) * 2021-04-14 2021-08-17 西南交通大学 Monte Carlo simulation-based earthquake landslide probability risk analysis method
CN114675325A (en) * 2022-03-13 2022-06-28 中国地震台网中心 Mindlin solution-based method for estimating permanent displacement of seismic surface fracture
CN115362391A (en) * 2020-03-06 2022-11-18 吉奥奎斯特系统公司 Marine seismic imaging

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103955620A (en) * 2014-05-13 2014-07-30 中国地质大学（北京） Engineering site earthquake hazard analysis method considering effect of potential earthquake source orientations
CN109001805A (en) * 2018-05-24 2018-12-14 深圳大学 A kind of earthquake source inverting Uncertainty Analysis Method, storage medium and server
CN109375253A (en) * 2018-12-13 2019-02-22 中国地震局地球物理研究所 Ground motion parameter evaluation method based on whole seismic structure maximum credible earthquakes
CN110334482A (en) * 2019-05-25 2019-10-15 成都理工大学 Mud-rock flow Dynamic flexibility evaluation method after shake based on material resource activity intensity
CN115362391A (en) * 2020-03-06 2022-11-18 吉奥奎斯特系统公司 Marine seismic imaging
CN113268852A (en) * 2021-04-14 2021-08-17 西南交通大学 Monte Carlo simulation-based earthquake landslide probability risk analysis method
CN113187556A (en) * 2021-05-06 2021-07-30 西南交通大学 Tunnel earthquake risk analysis method under fault dislocation based on total probability
CN114675325A (en) * 2022-03-13 2022-06-28 中国地震台网中心 Mindlin solution-based method for estimating permanent displacement of seismic surface fracture

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
考虑震群活动的川南工业区地震危险性分析;李昌珑 等;《地震地磁观测与研究》;377-379 *

Also Published As

Publication number Publication date
CN116127247A (en) 2023-05-16

Similar Documents

Publication Publication Date Title
EP2869096B1 (en) Systems and methods of multi-scale meshing for geologic time modeling
Shah et al. Modelling the effects of spatial variability in rainfall on catchment response. 1. Formulation and calibration of a stochastic rainfall field model
CN102944174A (en) Point cloud data processing method and system
CN102434210B (en) Method and system for monitoring underground engineering portrait information and monitoring information safely
CN106529755A (en) Mine geological resource reserve management method
CN110276732B (en) Mountain area point cloud cavity repairing method considering topographic characteristic line elements
CN107610021A (en) The comprehensive analysis method of environmental variance spatial and temporal distributions
CN104021267A (en) Geological disaster liability judgment method and device
CN104317772A (en) Method of quick geometric processing for Monte-Carlo particle transport on basis of spatial grid partitioning
CN107119657A (en) A kind of view-based access control model measures foundation ditch monitoring method
EP4036836A1 (en) Method, system for calculating operation acres of agricultural machinery
CN116359981A (en) Quick determination method for travel time of earthquake waves
CN103793939A (en) Local increasing type curved-surface reconstruction method of large-scale point cloud data
CN116127247B (en) Probability risk analysis and calculation method for coupling multiple seismic source models
Baehmann et al. Adaptive multiple-level h-refinement in automated finite element analyses
CN104794332A (en) Uncertainty analysis method for high-rise building wind response analysis models
CN117557681A (en) High-precision topographic map generation method and device based on multi-source mapping data
CN107507179B (en) Rock-soil mass quantitative analysis method based on GOCAD
CN112070129A (en) Ground settlement risk identification method, device and system
CN115455706A (en) Regional rock mass quality evaluation method considering unloading fracture effect and related assembly
CN110070234B (en) Earthquake landslide personnel death number prediction method and application thereof
US20150193707A1 (en) Systems and Methods for Estimating Opportunity in a Reservoir System
CN113970782A (en) Big data processing high and steep rock mass mining slope stability sound wave detection evaluation system
Bakula et al. The role of structural lines extraction from high-resolution digital terrain models in the process of height data reduction
CN113607920B (en) Reconstruction analysis method, experimental device and medium for sedimentary basin by magma bottom wall

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
GR01 Patent grant
GR01 Patent grant