A tomographic imaging method based on Kalman filtering

Abstract

The present invention discloses a tomography method based on Kalman filtering, characterized by comprising the following steps: obtaining seismic wave arrival information and its uncertainty; establishing an initial velocity model and its uncertainty; performing ray tracing based on the given initial velocity model to obtain the calculated travel time and ray path from each shot to each receiving point; forming a Kalman filter inversion equation based on the difference between the above-mentioned seismic wave arrival information and the calculated travel time, as well as the ray path, and solving the inversion equation to obtain an updated velocity model and the corresponding uncertainty; repeating the above steps until the obtained updated velocity model meets the accuracy requirements. By considering the seismic wave arrival information and the uncertainty of the initial velocity model in the inversion equation, the present invention can not only obtain an updated velocity model, but also obtain the uncertainty of the updated velocity model, which is helpful for subsequent geological interpretation and construction design.

Description

Technical Field

[0001] The present invention belongs to the technical field of geophysical inversion, and in particular relates to a method for performing travel-time tomography of underground media based on Kalman filtering. Background Art

[0002] Obtaining accurate velocity models and discontinuity distributions of subsurface media is a core issue in seismology, crucial for earthquake location, seismic attenuation imaging, microseismic monitoring, and subsequent earthquake prediction. Numerous theoretical and experimental studies have demonstrated that the physical properties of subsurface media change before and after earthquakes. Studying these changes in physical properties offers a new avenue for exploring earthquake prediction.

[0003] Seismic traveltime tomography is a technique for reconstructing the subsurface velocity structure using the observed travel time of seismic waves. After nearly 40 years of development, methods for velocity imaging based on seismic traveltime have made significant progress, including first-arrival tomography, combined first-arrival and reflection traveltime tomography, and multi-phase traveltime tomography. However, to date, almost all traveltime tomography methods can only provide a final inversion result. As is well known, geophysical inversion is a nonlinear inverse problem. Although the above methods have improved the inversion accuracy to varying degrees, none can fully accurately image the subsurface structure. Various factors, such as different shot arrangement and arrival time picking accuracy, can affect the inversion results. Therefore, the accuracy (or reliability) of the inversion results has become a bottleneck restricting the practical application of these inversion methods in subsequent address interpretation, construction design, and other practical applications.

[0004] Existing 4D time-lapse tomography methods combine seismic wave data received by monitoring stations at different time periods and, through continuous iterative optimization, bring model parameters close to the true solution. This method obtains model parameters for different time periods, thereby determining the variation of model parameters and providing the final results. However, it does not provide the uncertainty of the inversion results, making it difficult for personnel to conduct analysis and geophysical interpretation. To address these shortcomings, a solution is now provided.

[0005] 4D seismic traveltime tomography, as disclosed in the prior art, involves a process where, after given an initial model, the model parameters are continuously adjusted to ensure that the calculated seismic wave arrival times are sufficiently close to the observed arrival times. This process involves the following specific steps: 1. Reading the observed arrival times; 2. Given the initial model parameters; 3. Calculating the theoretical arrival times; 4. Calculating the traveltime residuals and determining whether these residuals meet the required accuracy. Uncertainty exists in each step of traveltime tomography, particularly in the input parameters of the initial velocity model and the observed arrival times, which has so far been difficult to effectively eliminate or avoid.

[0006] Therefore, those skilled in the art are committed to developing a tomography method based on Kalman filtering to address the above-mentioned deficiencies in the prior art. Summary of the Invention

[0007] In view of the above-mentioned defects of the prior art, the technical problem to be solved by the present invention is that the current prior art can only give the final inversion result, but cannot give the uncertainty of the inversion result, which is not convenient for staff analysis and geophysical interpretation. A method is proposed that can simultaneously give the underground velocity structure parameters and the uncertainty (or confidence level) of the velocity structure parameters.

[0008] To achieve the above object, the present invention provides a tomography method based on Kalman filtering, the specific steps of which are:

[0009] Step 1: Obtain the observation time information of the monitoring station and its uncertainty;

[0010] Step 2: Establish an initial velocity model and the uncertainty of the initial velocity model;

[0011] Step 3: Calculate the travel time of seismic waves and perform ray tracing based on the initial velocity model given in step 2 to obtain the arrival time and ray path of seismic waves from each shot point to each receiving point;

[0012] Step 4: The difference between the observed arrival time obtained in step 1 and the calculated arrival time in step 3, as well as the uncertainty of the observed arrival time and the uncertainty of the initial velocity model, constitutes a Kalman filter inversion equation. The inversion equation is solved to obtain an updated velocity model and the uncertainty of the updated velocity model.

[0013] Step 5: Repeat steps 3 and 4 until the updated velocity model meets the accuracy requirements;

[0014] Furthermore, in step 1, the observed time information includes the specific time when the seismic wave arrives at the monitoring station and the uncertainty of the observed time;

[0015] Furthermore, in step 2, the initial velocity model is: L kilometers wide and H kilometers deep, and is divided into dx meters * dz meters according to the velocity distribution according to each grid size, and the initial velocity model is described by a discrete grid of L / dx*H / dz;

[0016] Furthermore, in step 2, the establishment of the initial velocity model and the uncertainty of the initial velocity model are specifically as follows: when the initial velocity model is given, the uncertainty of the size of the initial velocity model on each discrete grid is given;

[0017] Furthermore, in step 3, the travel time calculation and ray tracing of seismic waves are performed based on the initial velocity model given in step 2. The travel time calculation and ray tracing method is a fast marching method, specifically:

[0018] When seismic waves propagate, they satisfy the Eikonal equation as shown in formula (1):

[0019]

[0020] in,

[0021] T(x, z) and V(x, z) are functions of the spatial coordinates (x, z), representing the travel time of the seismic wave from the source point 0 to the spatial point (x, z) and the velocity of the point, respectively;

[0022] Discretizing formula (1) yields:

[0023]

[0024] The travel time of point (i, j) is solved by equation (2). The specific process includes the x-direction solution process and the z-direction solution process:

[0025] Furthermore, the x-direction solution process includes:

[0026] If the travel time of points (i-1, j) and (i+1, j) is known, determine T(x i-1，j )<T(x i+1，j );

[0027] If T(x i-1，j )<T(x i+1，j ) condition is established, judge T(x i-2，j )<T(x i-1，j ) If the condition is met, then If the condition is not met, then

[0028] If T(x i-1，j )<T(x i+1，j ) condition is not met, judge T(x i-1，j )>T(x i+1，j ), if the condition is met, judge T(x i+2，j )<T(x i+1，j ),

[0029] If the condition is met, then

[0030] If the condition is not met, then

[0031] If the travel time of point (i-1, j) is known, judge T(x i-2，j )<T(x i-1，j), if the condition is met, then If the condition is not met, then

[0032]

[0033] If the travel time of point (i+1, j) is known, judge T(x i+2，j )<T(x i+1，j ); if the condition is met, then

[0034]

[0035] If the condition is not met, then

[0036] If none of the above conditions are met, then

[0037] The z-direction solution process includes:

[0038] If the travel time of points (i, j-1) and (i, j+1) is known, determine T(x i，j-1 )<T(x i，j+1 );

[0039] If T(x i，j-1 )<T(x i，j+1 ) condition is established, judge T(x i，j-2 )<T(x i，j-1 ) If the condition is met, then If the condition is not met, then

[0040] If T(x i，j-1 )<T(x i，j+1 ) condition is not met, judge T(x i，j-1 )>T(x i，j+1 ), if the condition is met, judge T(x i，j+2 )<T(x i，j+1 ),

[0041] If the condition is met, then

[0042] If the condition is not met, then

[0043] If the travel time of point (i, j-1) is known, judge T(x i，j-2 )<T(x i，j-1 ), if the condition is met, then If the condition is not met, then

[0044] If the travel time of point (i, j+1) is known, judge T(x i，j+2 )<T(x i，j+1); if the condition is met, then

[0045]

[0046] If the condition is not met, then

[0047] If none of the above conditions are met, then

[0048] in,

[0049]

[0050]

[0051] Furthermore, in step 3, the specific steps of calculating the seismic wave travel time are:

[0052] Step 3-1, assign the attribute value of all grid nodes to "far away point";

[0053] Step 3-2: Assign the travel time of the earthquake source point to 0, assign its attribute to "Completed Point", calculate the travel time of the four nodes on the grid where the earthquake source is located, and assign the attributes of the four surrounding nodes to "Points to be Updated";

[0054] Step 3-3: Find the point with the smallest travel time among the "points to be updated" attributes, assign it the value of "completed point", and calculate the travel time of the surrounding attributes of "far-away points" and "points to be updated". If the calculated point is a "far-away point", then after the calculation, update its attribute to "point to be updated". If the calculated point is a "point to be updated", then after the calculation, compare the original travel time of the point with the current calculated travel time, and select the smaller value as the travel time of the point.

[0055] Step 3-4: Repeat step 3-3 until the number of nodes with the attribute "points to be updated" is 0;

[0056] Furthermore, in step 4, the difference between the actual observed time and the calculated time, the uncertainty of the observed time, and the uncertainty of the initial model constitute a Kalman filter inversion equation. The inversion equation is solved to obtain the updated velocity model and the uncertainty of the updated velocity model. Specifically, the inversion equation based on the Kalman filter can be written as:

[0057]

[0058] in,

[0059] is the model velocity correction value of size N×1;

[0060] N is the total number of all grid nodes after the model is discretized;

[0061] is the modification of the current model based on the observation time of the k-th earthquake source point propagating to the i-th station, specifically including (δV1, δV2, L, δV3);

[0062] for The covariance matrix composed of the uncertainty of

[0063] is the sensitivity vector of the travel time from the kth earthquake source to the ith station with respect to velocity, specifically including

[0064]

[0065] is the uncertainty of the observation time from the kth earthquake source to the i-th station;

[0066] is the difference between the observed arrival time of the kth earthquake source to the i-th station and the calculated arrival time of the kth earthquake source to the i-th station calculated in step 3;

[0067] Furthermore, in step 4, formula (5) is solved by recursive iteration, as follows:

[0068] Step 4-1: Assume that the velocity of the jth grid node of the initial model m0 is V j ,j=1,2,…N, N is the total number of all mesh nodes after the model is discretized, and the initial model change amount is set to 0: The uncertainty corresponding to the velocity of each grid node is U j , j=1,2,…N, so the uncertainty covariance matrix Q0 of the initial model velocity can be written as:

[0069]

[0070] Step 4-2, in step 1, the observation time from the k-th shot point to the i-th receiving point is

[0071] Where sourceNum is the number of sources, receiverNum is the number of receiving points, and when it is observed The uncertainty is

[0072] Step 4-3, the calculation time from the kth shot point to the ith receiving point calculated in step 3 is

[0073] Step 4-4: Calculate the actual observed time from the kth source point to the ith receiving point When calculated The difference between

[0074]

[0075] Step 4-5: Calculate the sensitivity vector of the travel time from the kth earthquake source to the ith station with respect to velocity

[0076]

[0077] in,

[0078] L is the length of the ray passing through the mth unit of the model;

[0079] V m is the velocity distribution of the mth unit;

[0080] Step 4-6: Based on the uncertainty Q of the velocity model i-1 The first observation is used to modify the velocity model, Q0, and the uncertainty of the observation is (7) (8) Substitute into the Kalman filter inversion equation (5) to obtain the velocity model correction at this time: and its corresponding uncertainty

[0081] Step 4-7, Plus the updated velocity model based on the last observation When the first observation is made, the velocity model is modified to m0, and the updated velocity model is obtained

[0082] The velocity model is modified according to the next observation time until all observation times are reached;

[0083] By adopting the above scheme, the tomography method based on Kalman filtering disclosed in the present invention has the following advantages:

[0084] (1) The tomography method based on Kalman filtering of the present invention adopts a two-part Kalman filtering inversion strategy, one of which is a global Kalman filtering inversion strategy with time as the step, and the other is a local Kalman filtering inversion strategy with the number of iterations as the step;

[0085] The global Kalman filter inversion strategy uses the model parameters and uncertainties obtained from the previous observation data as prior information, and uses the observation data obtained at this time to calculate the likelihood function, thereby obtaining the posterior estimate of the model. The global Kalman filter inversion strategy is conducive to ensuring that a new updated model and corresponding uncertainty are obtained every time an observation data is received.

[0086] The local inversion strategy uses the model and uncertainty obtained from the previous update as prior information, and uses the new travel time residuals to calculate the likelihood function to obtain the posterior model. The local inversion strategy is conducive to ensuring the resolution of the inversion.

[0087] By combining the global Kalman filtering strategy with the local Kalman filtering strategy, the local Kalman filtering-based Kalman filtering-based tomography method of the present invention has both high temporal resolution and spatial resolution;

[0088] (2) Compared with the traditional 4D seismic travel time iterative inversion method, the tomography method based on Kalman filtering of the present invention avoids the selection of regularization parameters (inappropriate regularization parameters may lead to erroneous inversion results);

[0089] (3) The tomography method based on Kalman filtering of the present invention can provide uncertainty in the inversion results; the uncertainty in the inversion is of great significance to the evaluation of the later solution and geological interpretation;

[0090] (4) In the tomography method based on the local Kalman filter of the present invention, when the prior information is very certain about the velocity estimate of certain regions in the model, the elements of the corresponding regions in the initial covariance matrix Q0 can be set to zero to reduce the number of unknown parameters, avoid ill-conditioned problems and numerical artifacts, and greatly improve the inversion resolution;

[0091] In summary, the tomography method based on Kalman filtering disclosed in the present invention adopts a two-part Kalman filtering inversion strategy. One part is a global Kalman filtering inversion strategy with time as the step, which is conducive to ensuring that a new updated model and corresponding uncertainty are obtained every time an observation data is received; the other part is a local Kalman filtering inversion strategy with the number of iterations as the step, which is conducive to ensuring the resolution of the inversion; avoiding the selection of regularization parameters and the problem that may lead to errors in the inversion results. It can not only perform time-lapse fault imaging, but also reduce the uncertainty of the inversion results, with high accuracy and high resolution; at the same time, the results obtained can help later geological interpretation and be applied to monitoring changes in microseismic and natural seismic media, with high practicality and good applicability.

[0092] The concept, specific technical solutions and technical effects of the present invention will be further described below in conjunction with specific implementation methods to fully understand the purpose, characteristics and effects of the present invention. BRIEF DESCRIPTION OF THE DRAWINGS

[0093] Figure 1 3 is a flow chart of step 3 of a tomography method based on Kalman filtering of the present invention;

[0094] Figure 2 (a) is the present invention used to obtain the true velocity model and shot detection arrangement when the seismic wave from each shot point to the receiving point is observed;

[0095] Figure 2 (b) is the ray path diagram corresponding to a shot obtained based on the true velocity model;

[0096] Figure 2 (c) is the travel time contour map corresponding to a shot obtained based on the true velocity model;

[0097] Figure 3 This is the result diagram obtained by inversion in Example 1 of the present invention;

[0098] Figure 3 middle, Figure 3 (a) is the uncertainty of the initial velocity model (Q0), Figure 3 (b) is the Kalman filter inversion result diagram obtained based on the initial velocity model of 3km / s, the uncertainty of the initial velocity model (Q0) and the observation time. Figure 3 (c) is the uncertainty distribution diagram of the inversion results. Figure 3 (d) is the inversion result ( Figure 3 (b)) and the true model ( Figure 2 Error distribution diagram between (a));

[0099] Figure 4 is a schematic diagram of point (i, j) in step 3 of embodiment 1 of the present invention;

[0100] Figure 5 This is a Kalman filter inversion result diagram obtained by setting the initial model to a uniform model with a velocity of 3 km / s in Example 1 of the present invention;

[0101] Figure 5 middle, Figure 5 (a) is the uncertainty of the initial velocity model (Q0), Figure 5 (b) is the Kalman filter inversion result based on the initial velocity model of 3km / s. Figure 5 (c) is the uncertainty distribution diagram of the inversion results. Figure 5 (d) is the inversion result ( Figure 5(b)) and the true model ( Figure 2 Error distribution diagram between (a));

[0102] In the figure, the red inverted triangle is the receiving point, and the black dot is the shooting point;

[0103] Figure 3 In the model, the initial velocity is 3 km / s; DETAILED DESCRIPTION

[0104] The following describes several preferred embodiments of the present invention to make its technical content clearer and easier to understand. The present invention can be embodied in many different forms of embodiments, which are for illustrative purposes only and the scope of protection of the present invention is not limited to the embodiments mentioned herein.

[0105] Example 1: Using the Kalman filter-based tomography method of the present invention to perform earthquake travel-time tomography

[0106] Step 1: Obtain the observation time information of the monitoring station and its uncertainty; the observation time information includes the specific time when the seismic wave arrives at the monitoring station and the uncertainty of the observed time;

[0107] Step 2: Establish an initial velocity model and the uncertainty of the initial velocity model; the initial velocity model is: L kilometers wide and H kilometers deep, divided into dx meters * dz meters according to the velocity distribution according to the size of each grid, and the initial velocity model is described by a discrete grid of L / dx*H / dz; the initial velocity model and the uncertainty of the initial velocity model are specifically established by giving the uncertainty of the size of the initial velocity model on each discrete grid at the same time as the initial velocity model is given;

[0108] Step 3: Calculate the travel time of seismic waves and perform ray tracing based on the initial velocity model given in step 2 to obtain the arrival time and ray path of seismic waves from each shot point to each receiving point. Calculate the travel time of seismic waves and perform ray tracing based on the initial velocity model given in step 2. The travel time calculation and ray tracing method is a fast travel method, specifically:

[0109] When seismic waves propagate, they satisfy the Eikonal equation as shown in formula (1):

[0110]

[0111] in,

[0112] T(x, z) and V(x, z) are functions of the spatial coordinates (x, z), representing the travel time of the seismic wave from the source point 0 to the spatial point (x, z) and the velocity of the point, respectively;

[0113] Discretizing formula (1) yields:

[0114]

[0115] Solve by formula (2) Figure 4 The travel time of point (i, j) includes the following steps:

[0116] x-direction:

[0117] If the travel time of points (i-1, j) and (i+1, j) is known, determine T(x i-1，j )<T(x i+1，j );

[0118] If T(x i-1，j )<T(x i+1，j ) condition is established, judge T(x i-2，j )<T(x i-1，j ), if the condition is met, then If the condition is not met, then

[0119] If T(x i-1，j )<T(x i+1，j ) condition is not met, judge T(x i-1，j )>T(x i+1，j ), if the condition is met, judge T(x i+2，j )<T(x i+1，j ),

[0120] If the condition is met, then

[0121] If the condition is not met, then

[0122] If the travel time of point (i-1, j) is known, judge T(x i-2，j )<T(x i-1，j ), if the condition is met, then If the condition is not met, then

[0123] If the travel time of point (i+1, j) is known, judge T(x i+2，j )<T(x i+1，j ); if the condition is met, then

[0124]

[0125] If the condition is not met, then

[0126] If none of the above conditions are met, then

[0127] z-direction:

[0128] If the travel time of points (i, j-1) and (i, j+1) is known, determine T(x i，j-1 )<T(x i，j+1 );

[0129] If T(x i，j-1 )<T(x i，j+1 ) condition is established, judge T(x i，j-2 )<T(x i，j-1 ), if the condition is met, then If the condition is not met, then

[0130] If T(x i，j-1 )<T(x i，j+1 ) condition is not met, judge T(x i，j-1 )>T(x i，j+1 ), if the condition is met, judge T(x i，j+2 )<T(x i，j+1 ),

[0131] If the condition is met, then

[0132] If the condition is not met, then

[0133] If the travel time of point (i, j-1) is known, judge T(x i，j-2 )<T(x i，j-1 ), if the condition is met, then If the condition is not met, then

[0134]

[0135] If the travel time of point (i, j+1) is known, judge T(x i，j+2 )<T(x i，j+1 ); if the condition is met, then

[0136]

[0137] If the condition is not met, then

[0138] If none of the above conditions are met, then

[0139] The specific steps of calculating the seismic wave travel time are as follows:

[0140] Step 3-1, assign the attribute value of all grid nodes to "far away point";

[0141] Step 3-2: Assign the travel time of the earthquake source point to 0, assign its attribute to "Completed Point", calculate the travel time of the four nodes on the grid where the earthquake source is located, and assign the attributes of the four surrounding nodes to "Points to be Updated";

[0142] Step 3-3: Find the point with the smallest travel time among the "points to be updated" attributes, assign it the value of "completed point", and calculate the travel time of the surrounding attributes of "far-away points" and "points to be updated". If the calculated point is a "far-away point", then after the calculation, update its attribute to "point to be updated". If the calculated point is a "point to be updated", then after the calculation, compare the original travel time of the point with the current calculated travel time, and select the smaller value as the travel time of the point.

[0143] Step 3-4: Repeat step 3-3 until the number of nodes with the attribute "points to be updated" is 0;

[0144] Step 4: The difference between the observed arrival time obtained in step 1 and the calculated arrival time in step 3, as well as the uncertainty of the observed arrival time and the uncertainty of the initial velocity model, constitutes a Kalman filter inversion equation. The inversion equation is solved to obtain an updated velocity model and the uncertainty of the updated velocity model.

[0145] The difference between the actual observed arrival time and the calculated arrival time, the uncertainty of the observed arrival time, and the uncertainty of the initial model constitute the Kalman filter inversion equation. The inversion equation is solved to obtain the updated velocity model and the uncertainty of the updated velocity model. Specifically, the inversion equation based on the Kalman filter can be written as:

[0146]

[0147] in,

[0148] is the model velocity correction value of size N×1;

[0149] N is the total number of all grid nodes after the model is discretized;

[0150] is the modification of the current model based on the observation time of the k-th earthquake source point propagating to the i-th station, specifically including (δV1, δV2, L, δV3);

[0151] for The covariance matrix composed of the uncertainty of

[0152] is the sensitivity vector of the travel time from the kth earthquake source to the ith station with respect to velocity, specifically including

[0153]

[0154] is the uncertainty of the observation time from the kth earthquake source to the i-th station;

[0155] is the difference between the observed arrival time of the kth earthquake source to the i-th station and the calculated arrival time of the kth earthquake source to the i-th station calculated in step 3;

[0156] Formula (5) is solved by recursive iteration, as follows:

[0157] Step 4-1: Assume that the velocity of the jth grid node of the initial model m0 is V j ,j=1,2,…N, N is the total number of all mesh nodes after the model is discretized, and the initial model change amount is set to 0: The uncertainty corresponding to the velocity of each grid node is U j , j=1,2,…N, so the uncertainty covariance matrix Q0 of the initial model velocity can be written as:

[0158]

[0159] Step 4-2, in step 1, the observation time from the k-th shot point to the i-th receiving point is Where sourceNum is the number of sources, receiverNum is the number of receiving points, and when it is observed The uncertainty is

[0160] Step 4-3, the calculation time from the kth shot point to the ith receiving point calculated in step 3 is

[0161] Step 4-4: Calculate the actual observed time from the kth source point to the ith receiving point When calculated The difference between

[0162]

[0163] Step 4-5: Calculate the sensitivity vector of the travel time from the kth earthquake source to the ith station with respect to velocity

[0164]

[0165] in,

[0166] L is the length of the ray passing through the mth unit of the model;

[0167] V m is the velocity distribution of the mth unit;

[0168] Step 4-6: Based on the uncertainty Q of the velocity model i-1 The first observation is used to modify the velocity model, Q0, and the uncertainty of the observation is (7) (8) Substitute into the Kalman filter inversion equation (5) to obtain the velocity model correction at this time: and its corresponding uncertainty

[0169] Step 4-7, Plus the updated velocity model based on the last observation When the first observation is made, the velocity model is modified to m0, and the updated velocity model is obtained

[0170] The velocity model is modified according to the next observation time until all observation times are reached;

[0171] Step 5: Repeat steps 3 and 4 until the updated velocity model meets the accuracy requirements;

[0172] The results are Figure 3 As shown;

[0173] Furthermore, in step 4-1, if the initial model is only set to a uniform model with a velocity of 3 km / s, that is, there is an uncertainty of 0.3 km / s on each grid node, the inversion result is as follows: Figure 5 As shown;

[0174] Depend on Figure 5 、 Figure 3 The comparison of inversion results shows that Figure 3 The inversion of is significantly better at recovering anomalies, and the uncertainty is closer to zero. Figure 3 , Figure 5 The inversion of is more likely to produce numerical artifacts;

[0175] It is shown that when the prior information is very certain about the velocity estimates in certain regions of the model, setting the elements of the corresponding regions in the initial covariance matrix Q0 to zero can reduce the number of unknown parameters, avoid ill-conditioned problems and numerical artifacts, and greatly improve the inversion resolution;

[0176] In summary, by considering the uncertainty of the seismic wave arrival time information and the initial velocity model in the inversion equation, the present invention can not only obtain an updated velocity model, but also obtain the uncertainty of the updated velocity model, which is helpful for subsequent geological interpretation and construction design.

[0177] The above describes in detail the preferred embodiments of the present invention. It should be understood that numerous modifications and variations based on the concepts of the present invention can be made by those skilled in the art without inventive effort. Therefore, any technical solution that can be derived by those skilled in the art based on the concepts of the present invention through logical analysis, reasoning, or limited experimentation based on the existing technology should be within the scope of protection defined by the claims.

Claims

1. A tomographic imaging method based on Kalman filtering, characterized in that: The specific steps are: Step 1: Obtain the observation time information of the monitoring station and its uncertainty; Step 2: Establish an initial velocity model and the uncertainty of the initial velocity model; Step 3: Calculate the travel time of seismic waves and perform ray tracing based on the initial velocity model given in step 2 to obtain the arrival time and ray path of seismic waves from each shot point to each receiving point; Step 4: The difference between the observed arrival time obtained in step 1 and the calculated arrival time in step 3, as well as the uncertainty of the observed arrival time and the uncertainty of the initial velocity model, constitutes a Kalman filter inversion equation. The inversion equation is solved to obtain an updated velocity model and the uncertainty of the updated velocity model. Step 5: Repeat steps 3 and 4 until the updated velocity model meets the accuracy requirements.

2. The tomography method based on Kalman filtering according to claim 1, characterized in that: In the step 1, The observed time information includes the specific time when the seismic wave arrives at the monitoring station and the uncertainty of the observed time.

3. The tomography method based on Kalman filtering according to claim 1, characterized in that: In the step 2, The initial velocity model is: L kilometers wide and H kilometers deep. According to the velocity distribution, each grid size is divided into dx meters * dz meters. The initial velocity model is described by a discrete grid of L / dx*H / dz; The establishment of the initial velocity model and the uncertainty of the initial velocity model specifically include providing the uncertainty of the size of the initial velocity model on each discrete grid while providing the initial velocity model.

4. The tomography method based on Kalman filtering according to claim 1, wherein: In step 3, Based on the initial velocity model given in step 2, the seismic wave travel time calculation and ray tracing are performed. The travel time calculation and ray tracing method is a fast marching method, specifically: When seismic waves propagate, they satisfy the Eikonal equation as shown in formula (1): in, T(x,z) and V(x,z) are functions of the spatial coordinates (x,z), representing the travel time of the seismic wave from the source point O to the spatial point (x,z) and the velocity of the point, respectively; Discretizing formula (1) yields: The travel time of point (i, j) is solved by formula (2). The specific process includes the x-direction solution process and the z-direction solution process.

5. The tomography method based on Kalman filtering according to claim 4, characterized in that: The x-direction solution process includes: If the travel time of points (i-1, j) and (i+1, j) is known, determine T(x i-1,j ) <T(x i+1,j ); If T(x i-1,j ) <T(x i+1,j ) condition is established, judge T(x i-2,j ) <T(x i-1,j ) If the condition is met, then If the condition is not met, then If T(x i-1,j ) <T(x i+1,j ) condition is not met, judge Y(x i-1,j )>T(x i+1,j ), if the condition is met, judge T(x i+2,j ) <T(x i+1,j ), If the condition is met, then If the condition is not met, then If the travel time of point (i-1,j) is known, judge T(x i-2,j ) <T(x i-1,j ), if the condition is met, then If the condition is not met, then If the travel time of point (i+1,j) is known, judge T(x i+2,j ) <T(x i+1,j ); if the condition is met, then If the condition is not met, then If none of the above conditions are met, then The z-direction solution process includes: If the travel time of points (i, j-1) and (i, j+1) is known, determine T(x i,j-1 ) <T(x i,j+1 ); If T(x i,j-1 ) <T(x i,j+1 ) condition is established, judge T(x i,j-2 ) <T(x i,j-1 ) If the condition is met, then If the condition is not met, then If R(x i,j-1 ) <T(x i,j+1 ) condition is not met, judge T(x i,j-1 )>T(x i,j+1 ), if the condition is met, judge T(x i,j+2 ) <T(x i,j+1 ), If the condition is met, then If the condition is not met, then If the travel time of point (i,j-1) is known, judge T(x i,j-2 ) <T(x i,j-1 ), if the condition is met, then If the condition is not met, then If the travel time of point (i,j+1) is known, judge T(x i,j+2 ) <T(x i,j+1 ); if the condition is met, then If the condition is not met, then If none of the above conditions are met, then in, 6. The tomography method based on Kalman filtering according to claim 1, wherein: The specific steps of calculating the seismic wave travel time are as follows: Step 3-1, assign the attribute value of all grid nodes to "away point"; Step 3-2: Assign the travel time of the earthquake source point to 0 and its attribute to "Completed Point". Calculate the travel times of the four nodes on the grid where the earthquake source is located and assign the attributes of the four surrounding nodes to "Points to be Updated". Step 3-3: Find the point with the smallest travel time among the "points to be updated" attributes, assign it the "completed point" attribute, and calculate the travel times of the surrounding "faraway points" and "points to be updated" attributes. If the calculated point is a "faraway point", update its attribute to "point to be updated" after the calculation is complete. If the calculated point is a "point to be updated", compare the original travel time of the point with the currently calculated travel time after the calculation is complete, and select the smaller value as the travel time of the point. Step 3-4: Repeat step 3-3 until the number of nodes with the attribute "points to be updated" is 0.

7. The tomography method based on Kalman filtering according to claim 1, characterized in that: In the step 4, The difference between the actual observed arrival time and the calculated arrival time, as well as the uncertainty of the observed arrival time and the uncertainty of the initial model constitute the Kalman filter inversion equation. The updated velocity model and the uncertainty of the updated velocity model are obtained by solving the inversion equation. Specifically, the inversion equation based on the Kalman filter can be written as: in, is the model velocity correction value of size N×1; N is the total number of all grid nodes after the model is discretized; is the modification of the current model based on the observation time of the k-th earthquake source point propagating to the i-th station, specifically including (δV1, δV2, L, δV3); for The covariance matrix composed of the uncertainty of is the sensitivity vector of the travel time from the kth earthquake source to the ith station with respect to velocity, specifically including is the uncertainty of the observation time from the kth earthquake source to the i-th station; is the difference between the observed arrival time of the k-th earthquake source to the i-th station and the calculated arrival time of the k-th earthquake source to the i-th station calculated in step 3.

8. The tomography method based on Kalman filtering according to claim 1, wherein: In step 4, formula (5) is solved by recursive iteration, as follows: Step 4-1: Assume that the velocity of the jth grid node of the initial model m0 is V j ,j=1,2,…N, N is the total number of all mesh nodes after the model is discretized, and the initial model change amount is set to 0: The uncertainty corresponding to the velocity of each grid node is U j , j=1,2,…N, so the uncertainty covariance matrix Q0 of the initial model velocity can be written as: Step 4-2, in step 1, the observation time from the k-th shot point to the i-th receiving point is Where sourceNum is the number of sources, receiverNum is the number of receiving points, and when it is observed The uncertainty is Step 4-3, the calculation time from the kth shot point to the ith receiving point calculated in step 3 is Step 4-4: Calculate the actual observed time from the kth source point to the ith receiving point When calculated The difference between Step 4-5: Calculate the sensitivity vector of the travel time from the kth earthquake source to the ith station with respect to velocity in, L is the length of the ray passing through the mth unit of the model; V m is the velocity distribution of the mth unit; Step 4-6: Based on the uncertainty Q of the velocity model i-1 The first observation is used to modify the velocity model, Q0, and the uncertainty of the observation is (7) (8) Substitute into the Kalman filter inversion equation (5) to obtain the velocity model correction at this time: and its corresponding uncertainty Step 4-7, Plus the updated velocity model based on the last observation When the first observation is made, the velocity model is modified to m0, and the updated velocity model is obtained The velocity model is modified according to the next observation time until all observation times are reached and the cycle ends.

Images