Earthquake prediction method and system based on multi-physics coupling and quantum optimization

Abstract

The invention discloses an earthquake prediction method and system based on multi-physics coupling and quantum optimization, and relates to the technical field of earthquake prediction, and the method comprises the following steps: obtaining earth surface deformation observation data, geomagnetic observation data, gravity gradient data, and infrasonic wave monitoring data, and carrying out the preprocessing; performing multi-physics field coupling calculation on the preprocessed earth surface deformation observation data, geomagnetic observation data, gravity gradient data and infrasonic wave monitoring data; optimizing the fault network by adopting a D-Wave quantum annealing technology to find a lowest energy state to obtain a quantum optimization result; and generating an earthquake probability distribution diagram, an earthquake magnitude-time window table and a fault stress field three-dimensional model, thereby greatly expanding the breadth and depth of data, and remarkably improving the comprehensiveness of earthquake prediction.

Description

Technical Field

[0001] The present invention relates to the technical field of earthquake prediction, and more specifically, it relates to an earthquake prediction method and system based on multi-physical field coupling and quantum optimization. Background Art

[0002] An earthquake is a highly destructive natural disaster that can have a huge impact on human life, property, and social development. Accurate earthquake prediction can provide crucial support for early disaster prevention and mitigation work and reducing disaster losses, so it has always been the focus and difficulty of research in the field of earth science. The background art of this earthquake prediction method based on multi-physical field coupling and quantum optimization is as follows:

[0003] Early earthquake prediction mainly relied on statistical analysis of historical earthquake data, such as studying the frequency distribution of the time, location, and magnitude of earthquakes, etc., to speculate on the possibility of future earthquakes. However, this method is only based on historical experience and does not deeply explore the physical mechanism of earthquake occurrence. Moreover, the occurrence of earthquakes is complex and uncertain, making the reliability of statistical results limited.

[0004] Subsequently, a prediction method based on geological structure was developed. By studying the distribution and activity characteristics of faults, etc., potential earthquake hazard areas are judged. However, this method can only determine the general area where an earthquake may occur, and it is difficult to accurately predict the specific time and magnitude of an earthquake, and it is also difficult to comprehensively and accurately grasp the details of the geological structure deep underground.

[0005] In addition, some methods based on single physical field observations, such as monitoring surface deformation, geomagnetic changes, gravity anomalies, or infrasound, etc., although they can capture some information related to earthquakes, since the Earth's interior is a complex multi-physical field coupling system, single physical field observations often cannot comprehensively reflect all the changes during the earthquake gestation process, easily leading to information loss and misjudgment. Summary of the Invention

[0006] Aiming at the deficiencies of the existing technology, the purpose of the present invention is to provide an earthquake prediction method and system based on multi-physical field coupling and quantum optimization.

[0007] To achieve the above purpose, the present invention provides the following technical solutions:

[0008] An earthquake prediction method based on multi-physical field coupling and quantum optimization, the method includes the following steps:

[0009] Obtain surface deformation observation data, geomagnetic observation data, gravity gradient data, and infrasound monitoring data and then perform preprocessing;

[0010] Perform multi-physical field coupling calculations on the preprocessed surface deformation observation data, geomagnetic observation data, gravity gradient data, and infrasound monitoring data;

[0011] Use D-Wave quantum annealing technology to optimize the fault network to find the lowest energy state and obtain the quantum optimization result;

[0012] Generate a seismic probability distribution map, a magnitude-time window table, and a three-dimensional model of the fault stress field.

[0013] Preferably, after obtaining the surface deformation observation data, geomagnetic observation data, gravity gradient data, and infrasound monitoring data, perform preprocessing, specifically including:

[0014] The preprocessing of the surface deformation observation data, geomagnetic observation data, gravity gradient data, and infrasound monitoring data is specifically to clean and convert the formats of the surface deformation observation data, geomagnetic observation data, gravity gradient data, and infrasound monitoring data.

[0015] Preferably, the multi-physical field coupling calculation specifically includes:

[0016] Obtain the crustal stress field inversion result by processing and calculating the surface deformation observation data. The crustal stress field inversion result is specifically the Coulomb stress change and the fault slip rate;

[0017] Obtain the electromagnetic disturbance index by processing and calculating the geomagnetic observation data;

[0018] Obtain the gravity gradient anomaly value by processing and calculating the gravity gradient data;

[0019] Obtain the energy flux by processing and calculating the infrasound monitoring data.

[0020] Preferably, generating a seismic probability distribution map, a magnitude-time window table, and a three-dimensional model of the fault stress field specifically includes the following steps:

[0021] Generate a seismic probability distribution map based on the multi-physical field coupling calculation results;

[0022] Generate a magnitude-time window table based on the quantum optimization result;

[0023] Generate a three-dimensional model of the fault stress field based on the crustal stress field inversion result.

[0024] Preferably, obtaining the Coulomb stress change and the fault slip rate by processing and calculating the surface deformation observation data specifically includes the following steps:

[0025] The surface deformation observation data includes the horizontal displacement rate vector v=(v x , v y) Geological fault strike θ, geological fault dip angle δ, geological fault slip direction φ, and average geological fault slip rate

[0026] Obtain the elastic modulus μ of the geological medium, the Poisson's ratio ν of the geological medium, the friction coefficient μ of the geological fault f and the geological pore pressure p p ;

[0027] Calculate the strain rate tensor through the first calculation formula

[0028] Calculate the strain rate tensor in the vertical direction through the second calculation formula

[0029] Calculate the stress rate tensor through the third calculation formula where δ ij is the Kronecker symbol;

[0030] Calculate the Coulomb stress change Δσ c =Δτ + μ f (Δσ n -Δp p ) through the fourth calculation formula, where Δτ is the shear stress change, Δσ c is the normal stress change, Δp n is the pore pressure change, p t, n are the unit vectors in the tangential and normal directions of the fault plane respectively; Calculate the fault slip rate

[0031] through the fifth calculation formula where A is the friction parameter; In the formula, A is the friction parameter;

[0032] Construct the least squares objective function α, β are weight coefficients, N is the number of observation points, is the m th fault segment slip rate, is the prior slip rate estimated based on geological or geophysical methods, M the number of fault segments; α is the weight coefficient controlling the prior slip rate constraint strength;

[0033] Discretize the crustal medium using the finite element method or the boundary element method, and solve for the optimal stress field in combination with the gradient descent method or Bayesian inversion. The specific equation is Gm = d, where G is the Green's function matrix (transfer function from deformation to stress), m is the vector of stress parameters to be solved, and d is the vector of observed data (GNSS strain + fault slip rate).

[0034] Preferably, the electromagnetic disturbance index is calculated by processing and calculating the geomagnetic observation data, specifically including the following steps:

[0035] The geomagnetic observation data includes the vertical component B z obs (t) of the geomagnetic field actually observed at time t, and the long-term average value of the vertical component of the geomagnetic field during periods of no significant geomagnetic activity

[0036] Calculate the geomagnetic disturbance anomaly fluctuation parameter ΔB of the vertical component through the sixth calculation formula Z (t);

[0037] Calculate the ionospheric disturbance parameter ΔTEC(t) through the seventh calculation formula ; where μ TEC and σ TEC are the mean and standard deviation of the sliding time window, and TEC(t) represents the total electron content;

[0038] Construct a cylindrical spatial domain with a preset radius centered on the epicenter, and interpolate the geomagnetic disturbance anomaly fluctuation parameter ΔB Z (t) and the ionospheric disturbance parameter ΔTEC(t) to the grid to ensure spatio-temporal alignment;

[0039] Use a sliding time window to calculate the covariance matrix of ΔB Z (t) and ΔTEC(t) Extract the principal component direction as the electromagnetic coupling weight (α, β); where σ 2 ΔB is the variance of ΔB Z (t), σ 2 ΔT is the variance of ΔTEC(t), and σ ΔBΔT is the covariance of ΔB Z (t) and ΔTEC(t);

[0040] Calculate the electromagnetic disturbance index EMI(t) through the eighth calculation formula EMI(t) = α·ΔB z (t) + β·ΔTEC(t), where the weight coefficients satisfy α 2 + β 2 = 1;

[0041] When EMI(t) > 3σ EMI , it is determined as a valid abnormal signal.

[0042] Preferably, the gravity gradient anomaly value is obtained by processing and calculating the gravity gradient data, which specifically includes the following steps:

[0043] The gravity gradient data includes the observed gravity gradient tensor T ij and the normal gravity gradient U ij ;

[0044] The gravity gradient anomaly value δΓ ij = T ij - U ij (i, j = x, y, z) is calculated to obtain the gravity gradient anomaly value δΓ ij ;

[0045] Among them, the normal gravity gradient U is calculated through the tenth calculation formula ij ; where G represents the universal gravitational constant, M represents the mass of the earth, R represents the equatorial radius, r represents the distance value from the calculation point to the center of the earth, n represents the order of the spherical harmonic function, m represents the degree of the spherical harmonic function, represents the cosine term of the mass distribution symmetric about the longitude direction, represents the sine term of the mass distribution antisymmetric about the longitude direction, λ represents the longitude of the calculation point, θ = 90° - geographical latitude, represents the value of the normalized associated Legendre function at the colatitude θ;

[0046] The observed gravity gradient tensor T is calculated through the eleventh calculation formula ij ; where Δa i represents the acceleration change in the i-axis direction, Δx j represents the displacement change in the j-axis direction, represents the spatial gradient of the acceleration in the j-axis direction, represents the spatial gradient of the acceleration in the i-axis direction.

[0047] Preferably, the energy flux is obtained by processing and calculating the infrasound monitoring data, which specifically includes the following steps:

[0048] The infrasound monitoring data includes the sound pressure fluctuation p(t), the air density ρ 0 , the sound speed c, and the infrasound characteristic frequency band bandwidth f;

[0049] The energy flux φ is calculated through the twelfth calculation formula ; where ρ 0 is the static density of the medium, c is the sound speed, T is the average time, and p(t) is the instantaneous sound pressure.

[0050] Preferably, the D-Wave quantum annealing technology is adopted to optimize the tomography network to find the lowest energy state and obtain the quantum optimization result, which specifically includes the following steps:

[0051] Calculate the Hamiltonian H through the thirteenth calculation formula wherein, σ Problem , represents the spin state of the th qubit, J i is the coupling strength between qubits i and j, and h ij is the local magnetic field strength of the ith qubit; i

[0052] Calculate the transverse magnetic field Hamiltonian H through the fourteenth calculation formula wherein, σ initial ; wherein, σ i x is the Pauli-x operator of the ith qubit, and Γ is the transverse magnetic field strength;

[0053] Calculate through the fifteenth calculation formula H(t) = A(t)H initial + B(t)H problem ; wherein, A(t) is the weight of the transverse magnetic field Hamiltonian, and B(t) is the weight of the Hamiltonian;

[0054] If the evolution is slow enough, it will always remain in the instantaneous ground state and finally converge to the ground state of H problem ;

[0055] During the annealing process, the quantum system crosses the energy barrier through the tunneling rate formula Γ tunnel = e -s / h to jump out of the local minimum; wherein, s is the action of the tunneling path, h is the reduced Planck constant;

[0056] When t = T, the transverse magnetic field is turned off. At this time, A(T) = 0, and the system Hamiltonian is completely dominated by H problem , and the spin state {σ i} of the qubit is measured at this time, that is, the solution with the lowest energy is obtained.

[0057] A seismic prediction system based on multi-physical field coupling and quantum optimization includes

[0058] A preprocessing module that preprocesses the surface deformation observation data, geomagnetic observation data, gravity gradient data, and infrasound monitoring data after obtaining them;

[0059] A calculation module that performs multi-physical field coupling calculations on the preprocessed surface deformation observation data, geomagnetic observation data, gravity gradient data, and infrasound monitoring data;

[0060] A processing module that uses D-Wave quantum annealing technology to optimize the fault network to find the lowest energy state and obtain a quantum optimization result;

[0061] A generation module that generates a seismic probability distribution map, a magnitude-time window table, and a three-dimensional model of the fault stress field.

[0062] Compared with the prior art, the present invention has the following beneficial effects:

[0063] In the present invention, surface deformation observation data, geomagnetic observation data, gravity gradient data, and infrasound monitoring data are obtained through this application. These data respectively reflect the state changes of the interior and surface of the earth in detail from multiple different physical field dimensions such as crustal movement, earth's magnetic field, gravity field, and sound waves. The surface deformation observation data can accurately present the displacement and deformation of the crust, the geomagnetic observation data can sensitively capture the subtle anomalies of the earth's magnetic field, the gravity gradient data helps to detect the changes in the distribution of underground substances, and the infrasound monitoring data can effectively capture the infrasound signals that may be related to earthquakes. By comprehensively processing and deeply analyzing these rich and diverse data, various physical phenomena involved in the earthquake gestation process can be understood comprehensively and multi-level. Compared with the earthquake prediction methods that only rely on a single data source, this solution expands the breadth and depth of the data and significantly improves the comprehensiveness of earthquake prediction. BRIEF DESCRIPTION OF THE DRAWINGS

[0064] Figure 1 It is a schematic flow chart of a seismic prediction method based on multi-physical field coupling and quantum optimization proposed by the present invention;

[0065] Figure 2 It is a schematic module diagram of a seismic prediction system based on multi-physical field coupling and quantum optimization proposed by the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0066] Refer to Figures 1 to 2 as shown.

[0067] The embodiments are used to further illustrate a seismic prediction method and system based on multi-physical field coupling and quantum optimization proposed by the present invention.

[0068] A seismic prediction method based on multi-physical field coupling and quantum optimization, the method includes the following steps:

[0069] After obtaining surface deformation observation data, geomagnetic observation data, gravity gradient data, and infrasound monitoring data, perform preprocessing;

[0070] Perform multi-physical field coupling calculations on the preprocessed surface deformation observation data, geomagnetic observation data, gravity gradient data, and infrasound monitoring data;

[0071] The D-Wave quantum annealing technology is used to optimize the fault network to find the lowest energy state, obtaining a quantum optimization result;

[0072] A seismic probability distribution map, a magnitude-time window table, and a three-dimensional model of the fault stress field are generated.

[0073] This application comprehensively obtains multi-source information such as surface deformation observation data, geomagnetic observation data, gravity gradient data, and infrasound monitoring data. Different types of data reflect the changes inside and on the surface of the Earth from multiple physical field perspectives. For example, surface deformation observation data reflects the movement state of the crust, geomagnetic observation data reflects the anomalies of the Earth's magnetic field, gravity gradient data can show the subtle changes in the gravity field, and infrasound monitoring data can capture infrasound signals that may be related to earthquakes. By comprehensively processing and analyzing these data, various physical phenomena during the earthquake gestation process can be comprehensively understood, greatly improving the comprehensiveness of prediction compared to earthquake prediction methods based on a single data source.

[0074] Multi-physical field coupling calculations are performed on the preprocessed multi-source data. Through different calculation methods, key information related to earthquakes is extracted from various types of data, such as the inversion result of the crustal stress field (Coulomb stress change and fault slip rate) calculated from surface deformation observation data, the electromagnetic disturbance index obtained from geomagnetic observation data, the gravity gradient anomaly value obtained from gravity gradient data, and the energy flux obtained from infrasound monitoring data. The information from these different physical fields is interrelated and mutually corroborative, which can more accurately reveal the internal mechanism of earthquake occurrence, thereby enhancing the accuracy of earthquake prediction and reducing false alarms and missed alarms.

[0075] The D-Wave quantum annealing technology is used to optimize the fault network to find the lowest energy state. The quantum annealing technology utilizes the special properties of quantum systems and can quickly find the lowest energy state in a complex fault network model. Compared with traditional optimization algorithms, the calculation efficiency is greatly improved. At the same time, by precisely optimizing the fault network model, a more accurate quantum optimization result is obtained, and then a more accurate magnitude-time window table is generated based on this, providing a more accurate time and magnitude estimate for earthquake prediction.

[0076] Finally, an earthquake probability distribution map, a magnitude-time window table, and a three-dimensional model of the fault stress field are generated. The earthquake probability distribution map shows the likelihood of earthquakes occurring in different regions in an intuitive graphical manner, which helps relevant departments plan and deploy disaster prevention and mitigation measures in advance; the magnitude-time window table clearly gives the magnitude range and time interval in which an earthquake may occur, providing key information for emergency preparedness; the three-dimensional model of the fault stress field can clearly present the spatial distribution of fault stress, helping researchers deeply study the process of earthquake gestation and occurrence. These visualized result outputs greatly facilitate decision-making and application in practical work such as earthquake research, monitoring and early warning, and disaster prevention and mitigation.

[0077] With the help of various professional observation devices, such as the Global Navigation Satellite System (GNSS) used to obtain the horizontal displacement rate vectors of surface observation points in surface deformation observation data, geomagnetic observatories to obtain geomagnetic observation data, gravity measurement instruments to obtain gravity gradient data, and infrasound monitoring arrays to obtain infrasound monitoring data. Before an earthquake occurs, due to the movement of materials inside the Earth and the change of physical fields, these data will show characteristics different from the normal state, which are important bases for earthquake prediction.

[0078] Based on various earthquake-related parameters obtained from multi-physics field coupling calculations, such as Coulomb stress change, electromagnetic disturbance index, gravity gradient anomaly value, energy flux, etc., through specific algorithms and models, these parameters are correlated and analyzed with the probability of earthquake occurrence to generate an earthquake probability distribution map. In this map, different colors or patterns represent the probability of earthquake occurrence in different regions, intuitively showing the spatial distribution of the likelihood of earthquake occurrence.

[0079] Using the results obtained by quantum optimization and combining the relevant information obtained from multi-physics field coupling calculations, through further data analysis and model prediction, the magnitude range and corresponding time interval in which an earthquake may occur are determined to generate a magnitude-time window table. This table provides important time and magnitude estimation information for earthquake monitoring and early warning, helping relevant departments make preparations in advance.

[0080] Based on the crustal stress field inversion results obtained from surface deformation data processing, that is, information such as Coulomb stress change and fault slip rate, using three-dimensional modeling technology, a three-dimensional model of the fault stress field is constructed. This model can intuitively display the stress distribution of the fault in three-dimensional space, including information such as the magnitude and direction of stress, providing a powerful tool for deeply studying the mechanism of earthquake gestation and occurrence.

[0081] After obtaining surface deformation observation data, geomagnetic observation data, gravity gradient data, and infrasound monitoring data, preprocessing is carried out, specifically including:

[0082] The preprocessing of surface deformation observation data, geomagnetic observation data, gravity gradient data, and infrasound monitoring data specifically refers to cleaning and format conversion of surface deformation observation data, geomagnetic observation data, gravity gradient data, and infrasound monitoring data.

[0083] Cleaning surface deformation observation data, geomagnetic observation data, gravity gradient data, and infrasound monitoring data can effectively remove noise and invalid data generated by instrument errors, environmental interference, etc. in various types of data, ensuring the accuracy and reliability of the data, and providing a high-quality data basis for subsequent precise earthquake analysis. Format conversion can unify the formats of data from different sources, eliminate processing obstacles caused by data format differences, enable different types of data to be smoothly integrated and interacted, and improve the efficiency and compatibility of data processing.

[0084] High-quality and uniformly formatted data can more accurately reflect the true changes in the internal physical fields of the Earth in subsequent multi-physical field coupling calculations and other analysis steps. This helps to more accurately reveal the physical mechanisms during earthquake gestation, and then generate more reliable earthquake probability distribution maps, magnitude-time window tables, and three-dimensional models of fault stress fields, ultimately enhancing the overall reliability of earthquake prediction.

[0085] Use a variety of professional observation equipment to collect different types of data. For example, obtain surface deformation observation data through the Global Positioning System (GPS) or Interferometric Synthetic Aperture Radar (InSAR) technology, which can reflect the displacement and deformation of the Earth's crust surface; use geomagnetic observatories to collect geomagnetic observation data for monitoring changes in the Earth's magnetic field; collect gravity gradient data with the help of gravity gradiometers to detect gravity field differences caused by changes in underground material distribution; obtain infrasound monitoring data through infrasound monitoring stations to capture infrasound signals that may be related to earthquakes.

[0086] Perform cleaning operations on the four types of data obtained respectively. For surface deformation observation data, there may be abnormal displacement values caused by satellite signal interference, GPS receiver errors, etc., and these noise points need to be removed through filtering algorithms, etc.; geomagnetic observation data may be affected by solar activities, nearby electromagnetic interference, etc. and show abnormal fluctuations, and data smoothing and other methods need to be used to eliminate the interference; gravity gradient data may have deviations due to instrument precision problems or terrain effects, and correction and outlier identification processing need to be carried out; infrasound monitoring data may be mixed with invalid signals such as environmental noise, and effective infrasound signals need to be screened out through spectrum analysis and other means to ensure the purity of the data.

[0087] Since the data formats generated by different data acquisition devices and systems vary, for example, surface deformation observation data may be in a specific Geographic Information System (GIS) format, geomagnetic observation data may be in text format with different recording methods, etc. Therefore, according to the requirements of subsequent data processing and analysis software, these data need to be converted into a unified format, such as a common database format or a general data file format (CSV, JSON, etc.), so that various types of data can be integrated and analyzed on the same platform.

[0088] The multi-physical field coupling calculation specifically includes:

[0089] By processing and calculating the surface deformation observation data, the inversion results of the crustal stress field are obtained. The inversion results of the crustal stress field are specifically the Coulomb stress change and the fault slip rate;

[0090] By processing and calculating the geomagnetic observation data, the electromagnetic disturbance index is obtained;

[0091] By processing and calculating the gravity gradient data, the gravity gradient anomaly value is obtained;

[0092] By processing and calculating the infrasound monitoring data, the energy flux is obtained.

[0093] The surface deformation observation data of this application includes the horizontal displacement rate vector of the surface observation point, the strike of the geological fault, the dip angle of the geological fault, the sliding direction of the geological fault, and the average sliding rate of the geological fault, etc. These data reflect the movement state of the earth's crust surface and the geometric and kinematic characteristics of the fault.

[0094] At the same time, parameters such as the elastic modulus of the geological medium, the Poisson's ratio of the geological medium, the friction coefficient of the geological fault, and the geological pore pressure need to be obtained. These parameters describe the physical properties of the earth's crust materials and are crucial for calculating the crustal stress field.

[0095] First, through a specific first calculation formula, combining surface deformation observation data and relevant parameters, the strain rate tensor is calculated, which reflects the distribution of the crustal deformation rate. Then, through the second calculation formula, the strain rate tensor in the vertical direction is further obtained to analyze the deformation in the vertical direction. Next, using the third calculation formula, based on information such as the strain rate tensor, the stress rate tensor is calculated. The stress rate tensor represents the rate of change of stress with time. On this basis, through the fourth calculation formula, combining shear stress change, normal stress change, pore pressure change, and the tangential and normal unit vectors of the fault plane, the Coulomb stress change is calculated. The Coulomb stress change is a key indicator for judging whether a fault is prone to slip. Then, through the fifth calculation formula, the fault slip rate is calculated by combining relevant parameters. Finally, a least-squares objective function is constructed, the crustal medium is discretized using the finite element method or the boundary element method, and the optimal stress field is solved by combining the gradient descent method or Bayesian inversion, so as to obtain the accurate Coulomb stress change and fault slip rate, and complete the inversion of the crustal stress field.

[0096] Through the sixth calculation formula, the geomagnetic observation data collected are processed to calculate the geomagnetic disturbance anomaly fluctuation parameter in the vertical component, which reflects the abnormal change in the vertical direction of the geomagnetism. Then, through the seventh calculation formula, combining information such as the mean and standard deviation of the sliding time window, the ionospheric disturbance parameter is calculated. The ionospheric disturbance is closely related to geomagnetic activity and may also be related to the earthquake gestation process.

[0097] With the epicenter as the center, a cylindrical spatial domain with a preset radius is constructed, and the geomagnetic disturbance anomaly fluctuation parameter and the ionospheric disturbance parameter are interpolated into the grid to ensure spatio-temporal alignment, so that data at different positions and times can be analyzed in a unified framework.

[0098] The covariance matrix of the geomagnetic disturbance anomaly fluctuation parameter and the ionospheric disturbance parameter is calculated using a sliding time window, and the principal component direction is extracted as the electromagnetic coupling weight. Then, through the eighth calculation formula, the electromagnetic disturbance index is calculated by combining information such as the weight coefficient. When this index meets certain conditions, it is determined as an effective abnormal signal for subsequent earthquake analysis.

[0099] Gravity gradient data includes the observed gravity gradient tensor and the normal gravity gradient. The observed gravity gradient tensor is obtained through actual measurement by a gravity gradiometer, which reflects the current change of the gravity field; the normal gravity gradient is a theoretical value calculated based on the basic physical parameters and geometric shape of the earth and is used as a reference standard.

[0100] The observed gravity gradient tensor is obtained by the eleventh calculation formula, and the normal gravity gradient is obtained by the tenth calculation formula. Then, the gravity gradient anomaly is obtained by subtracting the normal gravity gradient from the observed gravity gradient tensor by the ninth calculation formula. The change of the gravity gradient anomaly may be related to the change of underground material distribution, crustal deformation and stress adjustment, etc., and is one of the important reference indicators for earthquake prediction.

[0101] Infrasound monitoring data include information such as the magnitude of sound pressure fluctuations, air density, sound speed, and bandwidth of the characteristic frequency band of infrasound, etc. These data reflect the basic physical characteristics of infrasound and the properties of the propagation medium.

[0102] The energy flux is calculated using the twelfth calculation formula, combined with data such as the amplitude of sound pressure fluctuations, air density, sound speed, and bandwidth of the characteristic frequency band of infrasound waves. The change in energy flux carried by infrasound waves during propagation may be related to the release and propagation of energy inside the earth before an earthquake occurs. Therefore, energy flux can be used as an important parameter to determine the possibility of an earthquake.

[0103] The earthquake probability distribution map, magnitude-time window table and fault stress field three-dimensional model are formed, which specifically includes the following steps:

[0104] Generate earthquake probability distribution map based on multi-physics field coupling calculation results;

[0105] Generate a magnitude-time window table based on quantum optimization results;

[0106] A three-dimensional model of fault stress field is generated based on the inversion results of the crustal stress field.

[0107] The earthquake probability distribution map generated by this application presents the possibility of earthquakes in different regions in a visual way, which can enable scientific researchers, disaster management departments and the public to intuitively understand the spatial distribution of earthquake risks, facilitate advance planning and take targeted disaster prevention and mitigation measures, such as strengthening infrastructure earthquake resistance construction in high-probability areas and rationally planning densely populated areas.

[0108] The magnitude-time window table is generated based on the results of quantum optimization, which can accurately estimate the magnitude range and time interval of possible earthquakes. This provides key time and magnitude references for emergency response agencies to reserve supplies and formulate evacuation plans in advance, which helps to improve the efficiency and effectiveness of emergency response and reduce the losses caused by earthquakes.

[0109] The three-dimensional model of fault stress field is constructed based on the inversion results of crustal stress field, which can intuitively present the stress distribution state of fault in three-dimensional space. This provides an important tool for geologists and seismologists to conduct in-depth research on fault mechanical properties and earthquake breeding mechanisms, and helps promote the development of earthquake science theory, thereby improving the accuracy and reliability of earthquake prediction.

[0110] The three results cooperate with each other. From the possibility of earthquake occurrence, magnitude-time prediction to fault stress analysis, they provide comprehensive and integrated decision-making support for earthquake monitoring and early warning, disaster prevention and control, scientific research, etc., forming a relatively complete earthquake-related information system.

[0111] Results such as Coulomb stress changes, electromagnetic disturbance indices, gravity gradient anomalies, and energy fluxes obtained from multi-physics field coupling calculations reflect abnormal information related to earthquakes from different physical field perspectives. These information comprehensively reflect the physical state changes in various regions of the Earth's interior and are important bases for judging the possibility of earthquake occurrence.

[0112] Using specialized data analysis models and algorithms, combine the above multi-physics field coupling calculation results with historical earthquake data, geological structure information, etc. Through learning and analyzing a large amount of data, establish the correlation between each physical quantity and the earthquake occurrence probability. Then, according to these correlation relationships, calculate and evaluate the earthquake occurrence probability for each region.

[0113] Based on the calculated earthquake occurrence probabilities for each region, use Geographic Information System (GIS) technology to visually represent the probability values on the map in the form of different colors, patterns, or contour lines, etc., thus generating an earthquake probability distribution map.

[0114] The result of the lowest energy state obtained after optimizing the fault network by D-Wave quantum annealing technology reflects a stable state of the fault system and a potential energy release pattern. This result contains information such as the interaction relationship between faults and stress concentration areas, which is closely related to the magnitude and occurrence time of earthquakes.

[0115] Combined with historical earthquake data and seismological theories, conduct in-depth analysis of the quantum optimization results. By establishing a mathematical model between earthquake magnitude and occurrence time and the state parameters of the fault network, predict the magnitude range and corresponding time interval of possible earthquakes in the future for a period of time.

[0116] Sort and arrange the predicted magnitude range and time interval, and generate a magnitude-time window table in a certain format to clearly show the time range in which earthquakes of different magnitudes may occur.

[0117] Results such as Coulomb stress changes and fault slip rates obtained from crustal stress field inversion accurately depict the stress state and movement of faults at different positions. These data are the basis for constructing a three-dimensional model of the fault stress field and can reflect the mechanical properties and deformation characteristics of faults.

[0118] Using 3D modeling software and technology, integrate the stress data and fault geometric information in the crustal stress field inversion results. Model the fault in 3D space according to a certain scale and coordinate system. In the model, represent the magnitude and direction of stress through different colors, textures, or isosurfaces, etc., so as to generate a 3D model that can intuitively display the distribution of the fault stress field.

[0119] By processing and calculating the surface deformation observation data, obtain the Coulomb stress change and fault slip rate, which specifically include the following steps:

[0120] The surface deformation observation data includes the horizontal displacement rate vector v=(v x , v y ) of the surface observation points, the strike θ of the geological fault, the dip angle δ of the geological fault, the slip direction φ of the geological fault, and the average slip rate of the geological fault

[0121] Obtain the elastic modulus μ of the geological medium, the Poisson's ratio ν of the geological medium, the friction coefficient μ f of the geological fault, and the pore pressure p p ;

[0122] Calculate the strain rate tensor through the first calculation formula

[0123] Calculate the strain rate tensor in the vertical direction through the second calculation formula

[0124] Calculate the stress rate tensor through the third calculation formula where δ ij is the Kronecker symbol;

[0125] Calculate the Coulomb stress change Δσ c = Δτ + μ f (Δσ n - Δp p ) through the fourth calculation formula, and obtain the Coulomb stress change Δσ c ; where Δτ is the shear stress change, Δσ n is the normal stress change, Δp p is the pore pressure change, t, n are the tangential and normal unit vectors of the fault plane respectively;

[0126] Calculate the fault slip rate through the fifth calculation formula In the formula, A is the friction parameter;

[0127] Construct the least - squares objective function α and β are weight coefficients, N is the number of observation points, is the m slip rate of the ith fault segment, and M is the prior slip rate estimated based on geological or geophysical methods; α is the weight coefficient controlling the constraint intensity of the prior slip rate;

[0128] Discretize the crustal medium using the finite - element method or the boundary - element method, and solve for the optimal stress field by combining the gradient - descent method or Bayesian inversion. The specific equation is Gm = d, where G is the Green's function matrix (the transfer function from deformation to stress), m is the vector of stress parameters to be solved, and d is the vector of observed data (GNSS strain + fault slip rate).

[0129] Through the comprehensive processing of various surface deformation observation data and related geological medium parameters and a series of precise calculations in this application, the Coulomb stress change and fault slip rate can be accurately obtained, and then the crustal stress field state can be inverted. This provides key information for in - depth understanding of the stress distribution within the crust and fault activities, helps to accurately grasp the mechanical changes during earthquake gestation, and greatly improves the accuracy of stress field inversion compared with single - data or simple calculation methods.

[0130] Accurate Coulomb stress change and fault slip rate data can be used to evaluate the activity and stability of faults. The Coulomb stress change is an important indicator for judging whether a fault is likely to slip, and the fault slip rate directly reflects the movement state of the fault. By analyzing these data, fault areas with high seismic risk can be identified in advance, providing important basis for earthquake early warning and disaster prevention, and helping relevant departments to take targeted measures in advance to reduce the seismic disaster risk.

[0131] Accurate crustal stress field inversion results provide more reliable input parameters for earthquake prediction models. Combining with other physical - field data and seismological research results, key information such as the likelihood of earthquake occurrence and magnitude can be analyzed more comprehensively, thus improving the accuracy and reliability of earthquake prediction and providing strong support for the improvement of earthquake monitoring and early - warning systems.

[0132] The electromagnetic disturbance index is obtained by processing and calculating geomagnetic observation data, specifically including the following steps:

[0133] The geomagnetic observation data includes the vertical component B z obs (t) of the geomagnetic field actually observed at time t, and the long - term average value of the vertical component of the geomagnetic field during non - significant geomagnetic activity

[0134] Through the sixth calculation formula Calculate the geomagnetic disturbance anomaly fluctuation parameter ΔB of the vertical component Z (t);

[0135] Through the seventh calculation formula Calculate the ionospheric disturbance parameter ΔTEC(t); where μ TEC and σ TEC are the mean and standard deviation of the sliding time window, and TEC(t) represents the total electron content;

[0136] With the epicenter as the center, construct a cylindrical spatial domain with a preset radius, and interpolate the geomagnetic disturbance anomaly fluctuation parameter ΔB Z (t) and the ionospheric disturbance parameter ΔTEC(t) to the grid to ensure spatio-temporal alignment;

[0137] Use a sliding time window to calculate the covariance matrix of ΔB Z (t) and ΔTEC(t) Extract the principal component direction as the electromagnetic coupling weight (α, β); where, σ 2 ΔB is the variance of ΔB Z (t), σ 2 ΔT is the variance of ΔTEC(t), and σ ΔBΔT is the covariance of ΔB Z (t) and ΔTEC(t);

[0138] Calculate the electromagnetic disturbance index EMI(t) through the eighth calculation formula EMI(t) = α·ΔB z (t) + β·ΔTEC(t), and the weight coefficients satisfy α 2 + β 2 = 1;

[0139] When EMI(t) > 3σ EMI , it is determined as a valid anomaly signal.

[0140] This application comprehensively uses a variety of calculation methods and data processing means, and can accurately extract electromagnetic anomaly information related to earthquakes from geomagnetic observation data. By calculating the geomagnetic disturbance anomaly fluctuation parameter and the ionospheric disturbance parameter, and then combining covariance matrix analysis and specific formula calculation of the electromagnetic disturbance index, weak but important electromagnetic anomaly signals can be effectively identified, providing a more sensitive and accurate indicator for earthquake monitoring.

[0141] Quantifying the electromagnetic coupling relationship: By using a sliding time window to calculate the covariance matrix and extracting the principal component direction as the electromagnetic coupling weight, the complex relationship between geomagnetic and ionospheric disturbances is quantified. This enables a deeper and more precise understanding of the electromagnetic coupling effect, helps to study the interaction mechanism of the electromagnetic physical fields within the Earth during earthquake gestation, and provides a new perspective and data support for the study of earthquake genesis.

[0142] Set clear judgment criteria (when the number of EMI(t) > 3σ EMI it is judged as a valid abnormal signal), which can screen out the signals that truly have potential earthquake indication significance from a large amount of electromagnetic data, reduce misjudgment and noise interference, improve the reliability and credibility of earthquake prediction, provide more valuable information for earthquake early warning work, and enable relevant departments to make more scientific decisions and deploy response measures.

[0143] Processing and calculating the gravity gradient data to obtain the gravity gradient anomaly value, which specifically includes the following steps:

[0144] The gravity gradient data includes the observed gravity gradient tensor T ij and the normal gravity gradient U ij ;

[0145] Through the ninth calculation formula δΓ ij = T ij - U ij (i, j = x, y, z) to calculate the gravity gradient anomaly value δΓ ij ;

[0146] Among them, the normal gravity gradient U is calculated through the tenth calculation formula ij ; where G represents the gravitational constant, M represents the mass of the Earth, R represents the equatorial radius, r represents the distance value from the calculation point to the center of the Earth, n represents the order of the spherical harmonic function, m represents the degree of the spherical harmonic function, represents the cosine term of the mass distribution symmetric about the longitude direction, represents the sine term of the mass distribution antisymmetric about the longitude direction, λ represents the longitude of the calculation point, θ = 90° - geographical latitude, represents the value of the normalized associated Legendre function at the colatitude θ;

[0147] Through the eleventh calculation formula to calculate the observed gravity gradient tensor T ij ; where Δa i represents the acceleration change in the i-axis direction, Δx j represents the displacement change in the j-axis direction, represents the spatial gradient of the acceleration in the j-axis direction, Represents the spatial gradient of acceleration in the i-axis direction.

[0148] Through the processing and calculation of gravity gradient data, this application can accurately obtain gravity gradient anomaly values. These anomaly values can reflect the inhomogeneity of underground material distribution, changes in geological structures, and possible stress adjustments, etc. Compared with other detection methods, it can more sensitively capture subtle physical changes underground, providing high-precision data support for studying the internal structure of the Earth and the earthquake gestation environment.

[0149] Accurate gravity gradient anomaly values are one of the important reference indicators for earthquake prediction. When the underground material distribution or stress state changes, corresponding abnormal changes will occur in the gravity gradient, and these changes may be related to the gestation and occurrence of earthquakes. By analyzing gravity gradient anomaly values and combining other physical field data (such as surface deformation, geomagnetism, etc.), the possibility of earthquake occurrence can be comprehensively judged from multiple angles, improving the accuracy and reliability of earthquake prediction.

[0150] The detailed processing and calculation of gravity gradient data help to deepen the understanding of the Earth's gravity field and improve the Earth science theory. At the same time, these data also have important application values in the fields of mineral resource exploration, geological structure research, etc., and can provide strong technical support for the research and practice in related fields, promoting the development of Earth science and its related industries.

[0151] The energy flux is obtained by processing and calculating the infrasound monitoring data, specifically including the following steps:

[0152] The infrasound monitoring data includes the acoustic pressure fluctuation momentum p(t), air density ρ 0 , sound speed c, and the bandwidth f of the characteristic frequency band of infrasound;

[0153] The energy flux φ is calculated through the twelfth calculation formula ; where ρ 0 is the static density of the medium, c is the sound speed, T is the average time, and p(t) is the instantaneous acoustic pressure.

[0154] This application can accurately calculate the energy flux from infrasound monitoring data. The change in the energy flux of infrasound may be closely related to the energy release and propagation inside the Earth before an earthquake. By monitoring and analyzing the energy flux, the energy information related to earthquake gestation can be captured, providing a new perspective and data basis for earthquake prediction, and helping to detect signs of earthquake activities in advance.

[0155] As a quantitative indicator, the energy flux can be used to continuously monitor the energy changes carried by infrasound waves. Combining with other seismic monitoring data (such as surface deformation, geomagnetism, etc.), it can provide a more comprehensive understanding of the physical state of the seismic activity area, assist in judging the intensity, scope and development trend of seismic activities, improve the accuracy and timeliness of seismic monitoring, and provide strong support for earthquake early warning and disaster prevention.

[0156] The processing of infrasound monitoring data and the calculation of energy flux contribute to deepening the research on the characteristics of infrasound waves and their relationship with earthquakes. This not only has important applications in the field of seismology, but also can be extended to other research and application scenarios involving infrasound waves, such as atmospheric physics, geological disaster monitoring, etc., promoting the technological development and theoretical improvement of related fields.

[0157] Using D-Wave quantum annealing technology to optimize the fault network to find the lowest energy state to obtain the quantum optimization result, which specifically includes the following steps:

[0158] Calculate the Hamiltonian H through the thirteenth calculation formula where σ Problem represents the spin state of the i-th qubit, J i is the coupling strength between qubits i and j, and h ij is the local magnetic field strength of the i-th qubit; i

[0159] Calculate the transverse magnetic field Hamiltonian H through the fourteenth calculation formula where σ initial is the Pauli-x operator of the i-th qubit, and Γ is the transverse magnetic field strength; i x

[0160] Calculate through the fifteenth calculation formula H(t) = A(t)H initial + B(t)H problem where A(t) is the weight of the transverse magnetic field Hamiltonian and B(t) is the weight of the Hamiltonian;

[0161] If the evolution is slow enough, it will always remain in the instantaneous ground state and finally converge to the ground state of H problem ;

[0162] During the annealing process, the quantum system crosses the energy barrier through the tunneling rate formula Γ tunnel = e -s / h to jump out of the local minimum; where s is the action of the tunneling path, h is the reduced Planck constant;

[0163] When t = T, the transverse magnetic field is turned off. At this time, A(T) = 0, and the system Hamiltonian is completely determined by Hproblem Dominant, at this time, the spin state of the quantum bit {σ i}, that is, the solution with the lowest energy is obtained.

[0164] This application uses D-Wave quantum annealing technology to quickly find the lowest energy state in a complex fault network model, greatly improving the computational efficiency compared to traditional optimization algorithms. When dealing with large-scale, high-dimensional fault network optimization problems, traditional methods may require a lot of time and computing resources, while quantum annealing technology can complete the calculation in a shorter time, providing timely data support for earthquake research and prediction.

[0165] By constructing and evolving a specific Hamiltonian, this technology can accurately optimize the fault network. It can fully consider factors such as the coupling strength between quantum bits and the local magnetic field strength, so as to more accurately reflect the actual physical state of the fault system and obtain optimization results that are more in line with the actual situation, providing a powerful tool for in-depth research on fault activity and earthquake incubation mechanisms.

[0166] The quantum optimization results obtained can be used to generate important data such as magnitude-time window tables, providing more accurate information for earthquake prediction. More accurate fault network optimization results can help to more accurately estimate the magnitude and time of possible earthquakes, improve the accuracy and reliability of earthquake prediction, and better guide earthquake early warning and disaster prevention and mitigation work, reducing losses caused by earthquakes.

[0167] An earthquake prediction system based on multi-physics field coupling and quantum optimization, including

[0168] The preprocessing module obtains the surface deformation observation data, geomagnetic observation data, gravity gradient data, and infrasound monitoring data and then performs preprocessing;

[0169] The calculation module performs multi-physics field coupling calculation on the pre-processed surface deformation observation data, geomagnetic observation data, gravity gradient data, and infrasound monitoring data;

[0170] The processing module uses D-Wave quantum annealing technology to optimize the fault network and find the lowest energy state to obtain the quantum optimization result;

[0171] The generation module generates earthquake probability distribution maps, magnitude-time window tables and three-dimensional models of fault stress fields.

[0172] Through the description of the above embodiments, those skilled in the art can clearly understand that each embodiment can be implemented by means of software plus a necessary general hardware platform, and of course, it can also be implemented by hardware. Based on such an understanding, the essence of the above technical solution, or the part that contributes to the prior art, can be embodied in the form of a software product. This computer software product can be stored in a computer-readable storage medium, such as ROM / RAM, magnetic disk, optical disk, etc., and includes several instructions to enable a computer device (which can be a personal computer, a server, or a network device, etc.) to execute the methods described in each embodiment or some parts of the embodiments.

[0173] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and are not intended to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those of ordinary skill in the art should understand that they can still modify the technical solutions recorded in the foregoing embodiments, or perform equivalent replacements on some of the technical features; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. An earthquake prediction method based on multi-physics field coupling and quantum optimization, characterized in that: The method comprises the following steps: Obtain surface deformation observation data, geomagnetic observation data, gravity gradient data, and infrasound monitoring data and perform preprocessing; Perform multi-physics field coupling calculations on the pre-processed surface deformation observation data, geomagnetic observation data, gravity gradient data, and infrasound monitoring data; Use D-Wave quantum annealing technology to optimize the fault network and find the lowest energy state to obtain quantum optimization results; Generate earthquake probability distribution maps, magnitude-time window tables and three-dimensional models of fault stress fields.

2. The earthquake prediction method based on multi-physical field coupling and quantum optimization according to claim 1, characterized in that: After obtaining the surface deformation observation data, geomagnetic observation data, gravity gradient data, and infrasound monitoring data, preprocessing is performed, including: The preprocessing of the surface deformation observation data, geomagnetic observation data, gravity gradient data, and infrasound monitoring data specifically includes cleaning and format conversion of the surface deformation observation data, geomagnetic observation data, gravity gradient data, and infrasound monitoring data.

3. The earthquake prediction method based on multi-physical field coupling and quantum optimization according to claim 1, characterized in that: The multi-physics field coupling calculation specifically includes: The crustal stress field inversion result is obtained by processing and calculating the surface deformation observation data, and the crustal stress field inversion result is specifically the Coulomb stress change and the fault slip rate; The electromagnetic disturbance index is obtained by processing and calculating the geomagnetic observation data; The gravity gradient anomaly value is obtained by processing and calculating the gravity gradient data; The energy flux is obtained by processing and calculating the infrasound monitoring data.

4. The earthquake prediction method based on multi-physical field coupling and quantum optimization according to claim 3 is characterized in that: Generate earthquake probability distribution map, magnitude-time window table and fault stress field three-dimensional model, specifically including the following steps: Generate earthquake probability distribution map based on multi-physics field coupling calculation results; Generate a magnitude-time window table based on quantum optimization results; A three-dimensional model of fault stress field is generated based on the inversion results of the crustal stress field.

5. The earthquake prediction method based on multi-physical field coupling and quantum optimization according to claim 3, characterized in that: The Coulomb stress change and fault slip rate are obtained by processing and calculating the surface deformation observation data, which specifically includes the following steps: The surface deformation observation data includes the horizontal displacement rate vector v of the surface observation point = (v x , v y ), the strike of the geological fault θ, the dip angle of the geological fault δ, the slip direction of the geological fault φ and the average slip rate of the geological fault Obtain geological medium elastic modulus μ, geological medium Poisson's ratio ν, geological fault friction coefficient μ f and geological pore pressure p p ; By the first calculation formula The strain rate tensor is calculated By the second calculation formula The strain rate tensor in the vertical direction is calculated By the third calculation formula The stress rate tensor is calculated in, δ ij is the Kronecker symbol; The fourth calculation formula Δσ c =Δτ+μ f (Δσ n -Δp p ) to calculate the Coulomb stress change Δσ c ; where Δτ is the shear stress change, Δσ n is the normal stress change, Δp p is the pore pressure change, t, n are the unit vectors tangent and normal to the fault plane, respectively; By the fifth calculation formula Calculate the fault slip rate Where A is the friction parameter; Constructing the least squares objective function α, β are weight coefficients, N is the number of observation points, For the m The slip rate of the fault segment, is the a priori slip rate estimated based on geological or geophysical methods, M The number of fault segments; α is the weight coefficient that controls the strength of the prior slip rate constraint; The crustal medium is discretized by the finite element method or boundary element method, and the optimal stress field is solved by combining the gradient descent method or Bayesian inversion. The specific equation is Gm=d, where G is the Green's function matrix, m is the stress parameter vector to be determined, and d is the observed data vector.

6. The earthquake prediction method based on multi-physical field coupling and quantum optimization according to claim 3, characterized in that: The electromagnetic disturbance index is obtained by processing and calculating the geomagnetic observation data, which specifically includes the following steps: The geomagnetic observation data includes the vertical component B of the geomagnetic field actually observed at time t. z obs (t), long-term average value of the vertical component of the geomagnetic field when there is no significant geomagnetic activity By the sixth calculation formula The vertical component of the geomagnetic disturbance anomaly fluctuation parameter ΔB is calculated Z (t); By the seventh calculation formula The ionospheric disturbance parameter ΔTEC(t) is calculated; where μ TEC and σ TEC is the mean and standard deviation of the sliding time window, TEC(t) represents the total electron content; With the epicenter as the center, a cylindrical space domain with a preset radius is constructed, and the geomagnetic disturbance abnormal fluctuation parameter ΔB is calculated. Z (t) and ionospheric disturbance parameter ΔTEC(t) are interpolated to the grid to ensure spatial and temporal alignment; Calculate ΔB using a sliding time window Z The covariance matrix of (t) and ΔTEC(t) The principal component direction is extracted as the electromagnetic coupling weight (α, β); where σ 2 ΔB ΔB Z The variance of (t), σ 2 ΔT is the variance of ΔTEC(t), σ ΔBΔT ΔB Z (t) and the covariance of ΔTEC(t); By the eighth calculation formula EMI(t) = α·ΔB z (t)+β·ΔTEC(t) to calculate the electromagnetic disturbance index EMI(t), and the weight coefficient satisfies α 2 +β 2 =1; When EMI(t)>3σ EMI , it is judged as a valid abnormal signal.

7. The earthquake prediction method based on multi-physical field coupling and quantum optimization according to claim 3, characterized in that: The gravity gradient anomaly value is obtained by processing and calculating the gravity gradient data, which specifically includes the following steps: The gravity gradient data includes the observed gravity gradient tensor T ij and the normal gravity gradient U ij ; By the ninth calculation formula δΓ ij =T ij -U ij (i, j = x, y, z) to calculate the gravity gradient anomaly δΓ ij ; Among them, through the tenth calculation formula Calculate the normal gravity gradient U ij ; Where G represents the gravitational constant, M represents the mass of the earth, R represents the equatorial radius, r represents the distance from the calculation point to the center of the earth, n represents the order of the spherical harmonic function, and m represents the degree of the spherical harmonic function. represents the cosine term of mass distribution symmetrical to the longitude direction, represents the sine term of mass distribution in the longitude direction, λ represents the longitude of the calculation point, θ=90°-geographic latitude, represents the value of the normalized joint Legendre function at co-latitude θ; By the eleventh calculation formula The observed gravity gradient tensor T is calculated ij ; where Δa i represents the acceleration change in the i-axis direction, Δx j represents the displacement change in the j-axis direction, represents the spatial gradient of acceleration in the j-axis direction, Represents the spatial gradient of acceleration in the i-axis direction.

8. The earthquake prediction method based on multi-physical field coupling and quantum optimization according to claim 3, characterized in that: The energy flux is obtained by processing and calculating the infrasound monitoring data, which specifically includes the following steps: The infrasound monitoring data includes the sound pressure fluctuation p(t), air density ρ0, sound speed c and infrasound characteristic frequency band bandwidth f; By the twelfth calculation formula The energy flux φ is calculated; where ρ0 is the static density of the medium, c is the speed of sound, T is the averaging time, and p(t) is the instantaneous sound pressure.

9. The earthquake prediction method based on multi-physical field coupling and quantum optimization according to claim 8, characterized in that: The D-Wave quantum annealing technology is used to optimize the fault network to find the lowest energy state to obtain the quantum optimization result, which specifically includes the following steps: By the thirteenth calculation formula The Hamiltonian H is calculated Problem , where σ i represents the spin state of the qubit, J ij is the coupling strength between quantum bits i and j, h i is the local magnetic field strength of the i-th quantum bit; By the fourteenth calculation formula The transverse magnetic field Hamiltonian H is calculated initial ; Among them, σ i x The Pauli-x operator of the i-th quantum bit, Γ is the transverse magnetic field strength; By the fifteenth calculation formula H(t)=A(t)H initial +B(t)H problem ; Where A(t) is the transverse magnetic field Hamiltonian weight, and B(t) is the Hamiltonian weight; If the evolution is slow enough, it will always remain in the instantaneous ground state and eventually converge to H problem The ground state of During the annealing process, the quantum system passes through the tunneling rate formula Γ tunnel =e -s / h Cross the energy barrier and jump out of the local minimum; among them, s is the action of the tunneling path, h is the reduced Planck constant; When t = T, the transverse magnetic field is turned off, A(T) = 0, and the system Hamiltonian is completely determined by H problem Dominant, at this time, the spin state of the quantum bit {σ i }, that is, the solution with the lowest energy is obtained.

10. An earthquake prediction system based on multi-physical field coupling and quantum optimization, applied to an earthquake prediction method based on multi-physical field coupling and quantum optimization as claimed in any one of claims 1 to 9, characterized in that: include The preprocessing module obtains the surface deformation observation data, geomagnetic observation data, gravity gradient data, and infrasound monitoring data and then performs preprocessing; The calculation module performs multi-physics field coupling calculation on the pre-processed surface deformation observation data, geomagnetic observation data, gravity gradient data, and infrasound monitoring data; The processing module uses D-Wave quantum annealing technology to optimize the fault network and find the lowest energy state to obtain the quantum optimization result; The generation module generates earthquake probability distribution maps, magnitude-time window tables and three-dimensional models of fault stress fields.

Images