A method and system for seismic wave full waveform inversion using the RBF-FD fusion differential method

Abstract

This application relates to the field of seismic wave full waveform inversion technology, specifically to a seismic wave full waveform inversion method and system based on the RBF-FD fusion differential method. The method includes: setting up artificial seismic sources and detection points and triggering source signals to acquire measured seismic waveform data; scanning the tunnel of the coal seam to be measured to generate point cloud data; constructing a three-dimensional geometric model of the tunnel based on the point cloud data, dividing it into regular surface regions, irregular surface regions, and transition zones; performing partitioned discrete solution to obtain partitioned discrete solution results; pre-setting the medium parameter field of the coal seam to be measured; performing forward seismic wave simulation based on the source signal, the current medium parameter field of the coal seam to be measured, and the partitioned discrete solution results to generate a forward simulated wavefield and extract the forward simulated waveform; comparing the residuals of the forward simulated waveform and the measured seismic waveform, calculating the residual norm, and updating the medium parameter field through an optimization algorithm until the accuracy requirements are met. This application can achieve high-fidelity seismic wave full waveform inversion.

Description

Technical Field

[0001] This application relates to the field of seismic wave full waveform inversion technology, specifically to a seismic wave full waveform inversion method and system using the RBF-FD fusion differential method. Background Technology

[0002] Safe and efficient coal mining relies heavily on precise geological structure detection technology. Seismic wave full-waveform inversion, as a core detection method, is crucial for identifying hidden hazards by comparing measured and simulated wavefield data to retrieve underground medium parameters.

[0003] Current techniques primarily employ the fixed-grid finite difference method for wavefield simulation and inversion calculations. This method offers high computational efficiency in regular regions and is easy to implement. For complex boundary problems, some studies have introduced meshless methods such as radial basis function generating finite difference (RBF-FD), leveraging their flexible node distribution to adapt to irregular geometries and improve boundary matching accuracy.

[0004] However, existing methods still have significant limitations: First, the fixed-grid finite difference method requires approximating irregular surfaces with stepped grids, leading to interference from artificial scattered waves and reducing the accuracy of forward modeling; Second, although the global application of the RBF-FD method can accurately match the boundary, it is difficult to meet the efficiency and stability requirements of actual engineering due to the large amount of computation, the susceptibility of the interpolation matrix to ill-conditioned phenomena, and the numerical instability caused by the co-location of variables; Third, existing tunnel modeling often simplifies irregular surfaces into regular geometric shapes, ignoring real geological disturbances and micro-undulations, which limits the fidelity of the inversion model. Summary of the Invention

[0005] To address the technical problems in existing seismic wave full waveform inversion techniques, such as the use of stepped grids to approximate irregular surfaces leading to artificial scattered wave errors, the global application of the RBF-FD method causing low computational efficiency and numerical instability, and the simplification of tunnel models ignoring real boundary morphology and reducing inversion accuracy, this application provides a seismic wave full waveform inversion method and system that integrates RBF-FD with the finite difference method. This method constructs a three-dimensional model from point cloud data and divides the region. The region is discretized using staggered-node RBF-FD and the finite difference method, combined with a dynamically coupled transition zone to eliminate artificial scattered waves and improve computational efficiency and stability, achieving high-fidelity medium parameter inversion.

[0006] In a first aspect, this application provides a seismic wave full waveform inversion method based on the RBF-FD fusion differential method, comprising the following steps:

[0007] S1. Set up artificial seismic sources and detection points in the coal seam to be tested and trigger the seismic source signal. Collect the measured seismic waveform data of the detection points through the seismic detectors set at the detection points;

[0008] S2. A laser scanning device is deployed in the tunnel of the coal seam to be measured. The laser scanning device emits pulsed lasers and receives reflected signals to generate point cloud data of the tunnel.

[0009] S3. Construct a three-dimensional geometric model of the tunnel based on point cloud data, divide the tunnel into regular surface regions and irregular surface regions in the three-dimensional geometric model, and set a transition zone at the junction of the regular surface regions and irregular surface regions;

[0010] S4. Perform partitioned discretization to obtain the partitioned discretization results, including:

[0011] In irregular surface regions, the staggered node RBF-FD method is used for discretization to obtain the RBF-FD weight coefficient matrix of the irregular regions.

[0012] The staggered mesh finite difference method is used to discretize and solve the problem in the regular surface region to obtain the finite difference coefficient matrix of the regular region.

[0013] The transition zone is solved by a dynamic coupling discretization method based on staggered node layout to obtain the mixed weight coefficient matrix;

[0014] S5. Preset the medium parameter field of the coal seam to be tested;

[0015] S6. Based on the source signal, the medium parameter field of the current coal seam to be measured, and the partitioned discrete solution results, perform forward modeling of seismic waves to generate a forward modeling wave field, and extract the forward modeling waveform of the detection point from the forward modeling wave field;

[0016] S7. Compare the residuals of the forward modeling waveforms at the detection points with the measured seismic waveforms, and calculate the residual norm;

[0017] If the residual norm is greater than the preset threshold, the medium parameter field is updated by the optimization algorithm and then the process returns to step S6.

[0018] If the residual norm is less than or equal to the preset threshold, the spatial distribution model of the current medium parameters is output.

[0019] It should be further noted that in step S2, the step of the laser scanning device emitting pulsed laser and receiving reflected signals to generate point cloud data of the tunnel includes:

[0020] S201. Configuration Each sampling point is targeted by a laser scanning device. Emit a pulsed laser beam and record the vertical direction angle of the pulsed laser beam. and horizontal direction angle ;

[0021] S202. Receive the laser signal reflected from the tunnel surface and record each sampling point. slant distance ;

[0022] S203. Calculate the coordinate information of each sampling point, including the three-dimensional coordinates of the point. Calculated based on the principle of laser ranging:

[0023]

[0024] S204. Summarize the coordinate information of all sampled points to generate the original point cloud dataset;

[0025] S205. Denoise the initial point cloud dataset to generate a valid point cloud dataset.

[0026] It should be further noted that step S202 also includes receiving the laser signal reflected from the tunnel surface and recording the reflection intensity at each sampling point. ;

[0027] The coordinate information in step S203 also includes the reflection intensity of the sampling point. .

[0028] It should be further explained that, in step S3, the step of dividing the tunnel's three-dimensional geometric model into regular surface regions and irregular surface regions, and setting a transition zone at the boundary between the regular surface regions and irregular surface regions, includes:

[0029] S301. Based on the effective tunnel surface point cloud data, calculate the value of each sampling point. normal vector ;

[0030] S302. For each sampling point Calculate its neighborhood point set rate of change of normal vector :

[0031] ;

[0032] S303. Set the threshold for the first normal vector change rate. The threshold of the rate of change of the second normal vector ,satisfy ;

[0033] Will and , Compare, among which:

[0034] like Then mark It belongs to an irregular surface area;

[0035] like Then mark It belongs to a regular surface area;

[0036] like Then mark It belongs to the transition zone;

[0037] S304. Perform morphological closing operations on the marking results to output continuously closed irregular surface regions. Regular surface areas and transition zone :

[0038]

[0039]

[0040] .

[0041] It should be further explained that in step S4, the specific steps for discretizing the irregular surface region using the staggered node RBF-FD method to obtain the RBF-FD weight coefficient matrix of the irregular region include:

[0042] S401. Pressure nodes and velocity nodes are arranged alternately in irregular surface areas;

[0043] The density of pressure nodes and velocity nodes is dynamically adjusted according to the local curvature. The higher the local curvature, the higher the density of pressure nodes and velocity nodes.

[0044] S402. For each pressure node as the central pressure node, select the M nearest neighbor pressure nodes to construct the support domain of the central pressure node;

[0045] Each velocity node is taken as the central velocity node, and the N nearest neighbor velocity nodes are selected to construct the support domain of the central velocity node.

[0046] S403. Set the radial basis functions:

[0047]

[0048]

[0049]

[0050]

[0051] in, This represents the Euclidean distance between the p-th pressure node and the central pressure node within the support domain of this central pressure node;

[0052] This represents the position vector of the central pressure node within the support domain of this central pressure node;

[0053] This represents the position vector of the p-th pressure node within the pressure node support domain of this center;

[0054] This represents the Euclidean distance between the q-th velocity node and the central velocity node within the support domain of this central velocity node;

[0055] This represents the position vector of the central velocity node within the support domain of this central velocity node;

[0056] This represents the position vector of the q-th velocity node in the support domain of this central velocity node;

[0057] Represents the calculation of Euclidean norm;

[0058] S404. Setting the pressure approximation function based on radial basis functions and velocity approximation function :

[0059]

[0060]

[0061] in, This represents the pressure RBF weight coefficient of the p-th pressure node in the pressure node support domain of this center;

[0062] This represents the velocity RBF weight coefficient of the q-th velocity node in the velocity node support domain of this center;

[0063] S405. For each central pressure node, the derivative of the pressure approximation function is used to obtain the expression for the spatial pressure derivative:

[0064]

[0065] in, For the p-th pressure node in the support domain of this center pressure node, the derivative weight coefficient is...

[0066] For spatial differential operators;

[0067] For each central velocity node, the derivative of the velocity approximation function is used to obtain the expression for the spatial velocity derivative:

[0068]

[0069] in, The derivative weight coefficient of the q-th velocity node in the velocity node support domain of this center. ;

[0070] S406. Solve for the spatial pressure derivative expression and the spatial velocity derivative expression to obtain the pressure derivative weighting coefficient between each central pressure node and all pressure nodes in its support domain, and the velocity derivative weighting coefficient between each central velocity node and all velocity nodes in its support domain.

[0071] S407. Summarize the pressure derivative weight coefficients of all central pressure nodes to form a pressure weight submatrix. :

[0072] The row index corresponds to the global number of the central pressure node;

[0073] The column index corresponds to the global number of all pressure nodes in the entire irregular region;

[0074] The element value is the pressure derivative weight coefficient between the center pressure node corresponding to the row index and the pressure node corresponding to the column index. For each row, a non-zero value is stored only in the column position corresponding to the pressure node within its support domain.

[0075] The velocity derivative weight coefficients of all central velocity nodes are summarized to form a velocity weight submatrix. :

[0076] The row index corresponds to the global number of the central velocity node;

[0077] The column index corresponds to the global number of all velocity nodes in the entire irregular region;

[0078] The element value is the velocity derivative weight coefficient between the center velocity node corresponding to the row index and the velocity node corresponding to the column index. For each row, a non-zero value is stored only in the column position corresponding to the velocity node within its supporting domain.

[0079] S408. Construct the RBF-FD weight coefficient matrix :

[0080] .

[0081] It should be further explained that in step S4, the specific steps for discretizing and solving the finite difference coefficient matrix of the regular surface region using the staggered mesh finite difference method include:

[0082] S411. Divide the regular surface area into a three-dimensional uniform grid, with the grid step size set according to the preset accuracy requirements;

[0083] S412. Arrange pressure nodes at each grid vertex and velocity nodes at the center of each grid cell;

[0084] S413. Based on the Taylor expansion formula, the wave field function is expanded into a series at each pressure node, and the pressure finite difference coefficients of each order spatial derivative are solved by fitting the least squares method.

[0085] Based on the Taylor expansion formula, the wave field function is expanded into a series at each velocity node, and the velocity finite difference coefficients of each order spatial derivative are solved by fitting the least squares method.

[0086] S414. Arrange the differential coefficients of each pressure node or velocity node according to the grid index to form a finite pressure differential coefficient submatrix. and finite velocity difference coefficient submatrix ;

[0087] S415. Constructing the finite difference coefficient matrix :

[0088] .

[0089] It should be further explained that in step S4, the step of using a dynamic coupling discretization method based on staggered node layout to discretize the transition region and obtain the mixed weight coefficient matrix includes:

[0090] S421. Arrange alternating pressure nodes and velocity nodes in the transition zone;

[0091] S422. Using each transition zone pressure node as the central transition zone pressure node, select the closest one. Each transition zone pressure node constructs a pressure hybrid support domain for the central transition zone pressure node;

[0092] Each transition zone velocity node is used as the central transition zone velocity node, and the closest one is selected. Each transition zone velocity node constructs the velocity hybrid support domain of the central transition zone velocity node;

[0093] S423. Calculate the pressure mixing weighting coefficient for each pressure node in the central transition zone and the velocity mixing weighting coefficient for each velocity node in the central transition zone:

[0094]

[0095]

[0096]

[0097]

[0098] in, This indicates that in this pressure-mixed support domain, the first... Pressure mixing weighting coefficient for each mixed pressure node;

[0099] This indicates that in this pressure-mixed support domain, the first... The pressure derivative weighting coefficients of each mixed pressure node are calculated using the staggered node RBF-FD method;

[0100] This indicates that in this pressure-mixed support domain, the first... The finite pressure difference coefficients of the mixed pressure nodes are calculated using the staggered mesh finite difference method;

[0101] This represents the Euclidean distance between the pressure node in the central transition zone of this pressure hybrid support domain and the nearest pressure node in the regular region;

[0102] This represents the Euclidean distance between the pressure node in the central transition zone of this pressure hybrid support domain and the nearest pressure node in the irregular region.

[0103] This indicates that in this velocity hybrid support domain, the first... The velocity mixing weight coefficient of each mixed velocity node;

[0104] This indicates that in this velocity hybrid support domain, the first... The velocity derivative weighting coefficients of each mixed velocity node are calculated using the staggered node RBF-FD method;

[0105] This indicates that in this velocity hybrid support domain, the first... The finite velocity difference coefficients of the mixed velocity nodes are calculated using the staggered mesh finite difference method;

[0106] This represents the Euclidean distance between the velocity node in the central transition zone of this velocity hybrid support domain and the nearest regular region velocity node;

[0107] This represents the Euclidean distance between the velocity node in the central transition zone of this velocity hybrid support domain and the nearest irregular region velocity node;

[0108] S424. Summarize the pressure mixing weight coefficients of all pressure nodes in the central transition zone to form a pressure mixing weight coefficient submatrix. :

[0109] The row index corresponds to the global number of the pressure node in the central transition zone;

[0110] The column index corresponds to the global number of all pressure nodes in the entire transition zone;

[0111] The element value is the pressure mixing weight coefficient between the pressure node in the central transition zone corresponding to the row index and the pressure node in the transition zone corresponding to the column index. For each row, a non-zero value is stored only in the column position corresponding to the pressure node in its support domain.

[0112] The velocity derivative weight coefficients of all central velocity nodes are summarized to form a velocity weight submatrix. :

[0113] The row index corresponds to the global number of the central velocity node;

[0114] The column index corresponds to the global number of all velocity nodes in the entire irregular region;

[0115] The element value is the velocity derivative weight coefficient between the center velocity node corresponding to the row index and the velocity node corresponding to the column index. For each row, a non-zero value is stored only in the column position corresponding to the velocity node within its supporting domain.

[0116] S425. Constructing the finite difference coefficient matrix :

[0117] .

[0118] It should be further noted that the medium parameter field includes the wave velocity field and the density field.

[0119] It should be further explained that in step S6, the steps of performing seismic wave forward modeling based on the source signal, the medium parameter field of the current coal seam to be measured, and the partitioned discrete solution results to generate the forward modeling wavefield include:

[0120] S601. Construct the governing equations describing the propagation of seismic waves. The governing equations include medium parameter fields, wave field variables, and spatial derivative terms.

[0121] S602. Discretize the control equations to obtain discretized control equations, including:

[0122] In irregular surface regions, the spatial derivative term is calculated using the RBF-FD weighting coefficient matrix;

[0123] In regular surface regions, the spatial derivative term is calculated using a finite difference coefficient matrix;

[0124] For the transition region, the spatial derivative term is calculated using the mixed weight coefficient matrix;

[0125] S603. Based on the source signal and medium parameter field, forward modeling is performed using discretized control equations to output the forward modeling wave field.

[0126] It should be further noted that the optimization algorithm used in step S6 is gradient descent.

[0127] Secondly, this application provides a seismic wave full waveform inversion system based on the RBF-FD fusion differential method, used to implement the aforementioned seismic wave full waveform inversion method, including:

[0128] The seismic data acquisition module is used to set up artificial seismic sources and detection points in the coal seam to be measured and trigger the seismic source signal. The measured seismic waveform data of the detection points is acquired by the seismic detectors set at the detection points.

[0129] The point cloud data acquisition and processing module is used to deploy laser scanning equipment in the tunnel of the coal seam to be measured. The laser scanning equipment emits pulsed laser and receives reflected signals to generate point cloud data of the tunnel.

[0130] The 3D model construction and region division module is used to construct a 3D geometric model of the tunnel based on point cloud data, divide the 3D geometric model of the tunnel into regular surface regions and irregular surface regions, and set a transition zone at the boundary between regular surface regions and irregular surface regions.

[0131] The partitioned discretization module is used to perform partitioned discretization and obtain the partitioned discretization results, including:

[0132] In irregular surface regions, the staggered node RBF-FD method is used for discretization to obtain the RBF-FD weight coefficient matrix of the irregular regions.

[0133] The staggered mesh finite difference method is used to discretize and solve the problem in the regular surface region to obtain the finite difference coefficient matrix of the regular region.

[0134] The transition zone is solved by a dynamic coupling discretization method based on staggered node layout to obtain the mixed weight coefficient matrix;

[0135] The medium parameter field initialization module is used to preset the medium parameter field of the coal seam to be tested.

[0136] The seismic wave forward modeling module is used to perform seismic wave forward modeling based on the source signal, the medium parameter field of the current coal seam to be measured, and the partitioned discrete solution results, generate the forward modeling wave field, and extract the forward modeling waveform of the detection point in the forward modeling wave field;

[0137] The residual calculation and parameter optimization module is used to compare the residuals of the forward modeling waveforms and the measured seismic waveforms at the detection points, calculate the residual norm, and update the medium parameter field through optimization algorithms until the accuracy requirements are met.

[0138] Thirdly, this application provides an electronic device, including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement the steps of the above-described seismic wave full waveform inversion method.

[0139] Fourthly, this application provides a storage medium storing a computer program, which, when executed by a processor, implements the steps of the above-described seismic wave full waveform inversion method.

[0140] As can be seen from the above technical solutions, this application has the following advantages:

[0141] 1. This application generates tunnel point cloud data and constructs a three-dimensional geometric model by laser scanning, dividing the surface into regular, irregular, and transitional regions. This solves the problem of traditional simplified models ignoring the true boundary morphology, and achieves high-fidelity tunnel modeling, which can restore the influence of microscopic undulation features on the wave field.

[0142] 2. This application solves the problems of low computational efficiency and numerical instability of the global RBF-FD method by using staggered node RBF-FD discretization for irregular surface regions, finite difference method discretization for regular surface regions, and dynamic coupling discretization for transition regions. It achieves synergistic optimization of computational efficiency and stability, while avoiding ill-conditioned interpolation matrix.

[0143] 3. This application solves the problem of wavefield distortion caused by the introduction of artificial scattered waves by accurately matching the boundary in irregular surface regions using the staggered node RBF-FD method, replacing the traditional stepped mesh approximation. This achieves a substantial improvement in the accuracy of wavefield simulation in complex boundary regions and effectively suppresses the interference of false waveforms on inversion.

[0144] 4. This application solves the limitation of boundary modeling error accumulation in traditional inversion by performing seismic wave forward modeling based on partitioned discretization results and medium parameter fields, and by combining residual comparison and iterative updating of medium parameters. It achieves high-resolution medium parameter field inversion and provides reliable imaging support for the identification of hidden disasters. Attached Figure Description

[0145] To more clearly illustrate the technical solution of this application, the accompanying drawings used in the description will be briefly introduced below. Obviously, the accompanying drawings described below are only some embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0146] Figure 1 This is a flowchart of a seismic wave full waveform inversion method using the RBF-FD fusion differential method in one embodiment of this application.

[0147] Figure 2 This is a schematic block diagram of a seismic wave full waveform inversion system using the RBF-FD fusion differential method in one embodiment of this application.

[0148] Figure 3 This is a schematic diagram of the hardware structure of an electronic device in one embodiment of this application. Detailed Implementation

[0149] To make the purpose, features, and advantages of this application more apparent and understandable, specific embodiments and accompanying drawings will be used to clearly and completely describe the technical solution protected by this application. Obviously, the embodiments described below are only some embodiments of this application, and not all embodiments. Based on the embodiments in this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0150] The following will describe in detail the seismic wave full waveform inversion method involved in this application. Specific details such as particular system structures and techniques are presented for illustrative purposes and not for limitation, in order to provide a thorough understanding of the embodiments of this application. However, those skilled in the art will understand that this application can also be implemented in other embodiments without these specific details.

[0151] In the seismic wave full waveform inversion method involved in this application, the term "comprising" indicates the presence of the described feature, whole, step, operation, element, and / or component, but does not exclude the presence or addition of one or more other features, wholes, steps, operations, elements, components, and / or sets thereof. The terms "comprising," "including," "having," and variations thereof all mean "including but not limited to," unless otherwise specifically emphasized.

[0152] To facilitate a clear description of the technical solutions of this application, the terms "first" and "second" are used to distinguish identical or similar items with essentially the same function and effect. Those skilled in the art will understand that the terms "first" and "second" do not limit the quantity or execution order, and that the terms "first" and "second" do not necessarily imply that they are different.

[0153] The terms "one embodiment" or "some embodiments" used in this application mean that one or more embodiments of this application include the specific features, structures, or characteristics described in that embodiment. Therefore, the terms "in one embodiment," "in some embodiments," "in other embodiments," "in still other embodiments," etc., appearing in different parts of this application do not necessarily refer to the same embodiment, but rather mean "one or more, but not all, embodiments," unless otherwise specifically emphasized.

[0154] The technical solutions in the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings.

[0155] The seismic wave full waveform inversion method provided in this application embodiment is executed by a computer device. Correspondingly, the seismic wave full waveform inversion system of the RBF-FD fusion differential method runs in the computer device.

[0156] Figure 1 This is a flowchart of a seismic wave full waveform inversion method based on the RBF-FD fusion differential method according to an embodiment of this application. Figure 1 The executing entity can be a seismic wave full waveform inversion system. Depending on different requirements, the order of the steps in this flowchart can be changed, and some can be omitted.

[0157] like Figure 1 As shown, the RBF-FD fusion differential method for seismic wave full waveform inversion includes:

[0158] Step S1: Set up an artificial seismic source and detection point in the coal seam to be tested and trigger the seismic source signal. Collect the measured seismic waveform data of the detection point through the seismic detector set at the detection point.

[0159] Step S2: A laser scanning device is deployed in the tunnel of the coal seam to be tested. The laser scanning device emits pulsed lasers and receives reflected signals to generate point cloud data of the tunnel.

[0160] By setting up artificial seismic sources and detection points to collect measured seismic waveform data, and deploying laser scanning equipment in the tunnel to emit pulsed lasers and receive reflected signals to generate point cloud data, comprehensive data acquisition of the seismic wave field and tunnel geometry was achieved. This provides a high-precision input foundation for subsequent model construction and inversion, ensuring the reliability and integrity of the data source.

[0161] In some specific embodiments, the step of generating point cloud data of the tunnel by emitting pulsed laser light and receiving reflected signals through a laser scanning device includes:

[0162] S201. Configuration Each sampling point is targeted by a laser scanning device. Emit a pulsed laser beam and record the vertical direction angle of the pulsed laser beam. and horizontal direction angle ;

[0163] S202. Receive the laser signal reflected from the tunnel surface and record each sampling point. slant distance ;

[0164] S203. Calculate the coordinate information of each sampling point, including the three-dimensional coordinates of the point. Calculated based on the principle of laser ranging:

[0165]

[0166] S204. Summarize the coordinate information of all sampled points to generate the original point cloud dataset;

[0167] S205. Denoise the initial point cloud dataset to generate a valid point cloud dataset.

[0168] By limiting the generation steps of point cloud data, a high-precision, noise-free point cloud dataset was generated, providing reliable input for the construction of the tunnel's 3D geometric model, avoiding model distortion caused by noise interference, and ensuring the accuracy and integrity of the subsequent inversion base data.

[0169] In some specific embodiments, step S202 further includes receiving the laser signal reflected from the tunnel surface and recording the reflection intensity at each sampling point. ;

[0170] The coordinate information in step S203 also includes the reflection intensity of the sampling point. .

[0171] By receiving and recording the reflection intensity information of each sampling point and incorporating it into the point cloud coordinate information, multidimensional enhancement of point cloud data is achieved, providing additional surface reflection characteristic parameters, assisting in tunnel surface division and medium parameter inversion analysis, improving data utilization and the physical authenticity of model construction, and avoiding the limitations brought by single coordinate information.

[0172] Step S3: Construct a three-dimensional geometric model of the tunnel based on the point cloud data, divide the tunnel into regular surface areas and irregular surface areas in the three-dimensional geometric model, and set a transition zone at the junction of the regular surface areas and irregular surface areas.

[0173] By constructing a three-dimensional geometric model of the tunnel based on point cloud data, dividing the surface into regular and irregular regions, and setting transition zones at the boundaries, the structured partitioning of the tunnel surface is realized. This provides a clear framework for discrete solution of the partitions, ensuring that different regions are adapted to the corresponding solution methods and avoiding inconsistencies in model processing.

[0174] In some specific embodiments, the step of dividing the tunnel's three-dimensional geometric model into regular surface regions and irregular surface regions, and setting a transition zone at the boundary between the regular surface regions and irregular surface regions, includes:

[0175] S301. Based on the effective tunnel surface point cloud data, calculate the value of each sampling point. normal vector ;

[0176] S302. For each sampling point Calculate its neighborhood point set rate of change of normal vector :

[0177] ;

[0178] S303. Set the threshold for the first normal vector change rate. The threshold of the rate of change of the second normal vector ,satisfy ;

[0179] Will and , Compare, among which:

[0180] like Then mark It belongs to an irregular surface area;

[0181] like Then mark It belongs to a regular surface area;

[0182] like Then mark It belongs to the transition zone;

[0183] S304. Perform morphological closing operations on the marking results to output continuously closed irregular surface regions. Regular surface areas and transition zone :

[0184]

[0185]

[0186] .

[0187] By calculating the normal vector of each sampling point, setting a threshold based on the rate of change of the normal vector of the neighborhood point set to divide the surface into regular, irregular, and transition regions, and applying morphological closing operations to output continuous closed partitions, the automatic and objective partitioning of the tunnel surface is realized. This ensures the scientific division of regular and irregular regions, reduces subjective judgment errors, and provides a structured basis for subsequent partitioning and discrete solution.

[0188] Step S4: Perform partitioned discretization to obtain the partitioned discretization results, including:

[0189] In irregular surface regions, the staggered node RBF-FD method is used for discretization to obtain the RBF-FD weight coefficient matrix of the irregular regions.

[0190] The staggered mesh finite difference method is used to discretize and solve the problem in the regular surface region to obtain the finite difference coefficient matrix of the regular region.

[0191] The transition zone is solved by a dynamic coupling discretization method based on staggered node layout, and the mixed weight coefficient matrix is ​​obtained.

[0192] By employing the staggered-node RBF-FD method in irregular surface regions, the staggered-grid finite difference method in regular surface regions, and the dynamic coupling discretization method in transition regions for discretization, and generating corresponding weight coefficient matrices, efficient discretization processing for different surface characteristics is achieved, improving overall solution accuracy and computational efficiency, and ensuring the robustness of the inversion basic steps.

[0193] In some specific embodiments, the steps for discretizing and solving the irregular surface region using the staggered node RBF-FD method to obtain the RBF-FD weight coefficient matrix of the irregular region include:

[0194] S401. Pressure nodes and velocity nodes are arranged alternately in irregular surface areas;

[0195] The density of pressure nodes and velocity nodes is dynamically adjusted according to the local curvature. The higher the local curvature, the higher the density of pressure nodes and velocity nodes.

[0196] S402. For each pressure node as the central pressure node, select the M nearest neighbor pressure nodes to construct the support domain of the central pressure node;

[0197] Each velocity node is taken as the central velocity node, and the N nearest neighbor velocity nodes are selected to construct the support domain of the central velocity node.

[0198] S403. Set the radial basis functions:

[0199]

[0200]

[0201]

[0202]

[0203] in, This represents the Euclidean distance between the p-th pressure node and the central pressure node within the support domain of this central pressure node;

[0204] This represents the position vector of the central pressure node within the support domain of this central pressure node;

[0205] This represents the position vector of the p-th pressure node within the pressure node support domain of this center;

[0206] This represents the Euclidean distance between the q-th velocity node and the central velocity node within the support domain of this central velocity node;

[0207] This represents the position vector of the central velocity node within the support domain of this central velocity node;

[0208] This represents the position vector of the q-th velocity node in the support domain of this central velocity node;

[0209] Represents the calculation of Euclidean norm;

[0210] S404. Setting the pressure approximation function based on radial basis functions and velocity approximation function :

[0211]

[0212]

[0213] in, This represents the pressure RBF weight coefficient of the p-th pressure node in the pressure node support domain of this center;

[0214] This represents the velocity RBF weight coefficient of the q-th velocity node in the velocity node support domain of this center;

[0215] S405. For each central pressure node, the derivative of the pressure approximation function is used to obtain the expression for the spatial pressure derivative:

[0216]

[0217] in, For the p-th pressure node in the support domain of this center pressure node, the derivative weight coefficient is...

[0218] For spatial differential operators;

[0219] For each central velocity node, differentiating the velocity approximation function yields the expression for the spatial velocity derivative:

[0220]

[0221] in, The derivative weight coefficient of the q-th velocity node in the velocity node support domain of this center. ;

[0222] S406. Solve for the spatial pressure derivative expression and the spatial velocity derivative expression to obtain the pressure derivative weighting coefficient between each central pressure node and all pressure nodes in its support domain, and the velocity derivative weighting coefficient between each central velocity node and all velocity nodes in its support domain.

[0223] S407. Summarize the pressure derivative weight coefficients of all central pressure nodes to form a pressure weight submatrix. :

[0224] The row index corresponds to the global number of the central pressure node;

[0225] The column index corresponds to the global number of all pressure nodes in the entire irregular region;

[0226] The element value is the pressure derivative weight coefficient between the center pressure node corresponding to the row index and the pressure node corresponding to the column index. For each row, a non-zero value is stored only in the column position corresponding to the pressure node within its support domain.

[0227] The velocity derivative weight coefficients of all central velocity nodes are summarized to form a velocity weight submatrix. :

[0228] The row index corresponds to the global number of the central velocity node;

[0229] The column index corresponds to the global number of all velocity nodes in the entire irregular region;

[0230] The element value is the velocity derivative weight coefficient between the center velocity node corresponding to the row index and the velocity node corresponding to the column index. For each row, a non-zero value is stored only in the column position corresponding to the velocity node within its supporting domain.

[0231] S408. Construct the RBF-FD weight coefficient matrix :

[0232] .

[0233] By staggering pressure and velocity nodes in irregular surface regions, dynamically adjusting node density based on local curvature, constructing radial basis function approximations by selecting support domains, solving spatial derivative weight coefficients, and forming an RBF-FD weight coefficient matrix, high-precision discrete solutions for irregular regions are achieved, improving the stability and adaptability of spatial derivative calculations and avoiding the numerical divergence problem of traditional methods in complex curvature regions.

[0234] In some specific embodiments, the steps for discretizing and solving the finite difference coefficient matrix of the regular surface region using the staggered mesh finite difference method include:

[0235] S411. Divide the regular surface area into a three-dimensional uniform grid, with the grid step size set according to the preset accuracy requirements;

[0236] S412. Arrange pressure nodes at each grid vertex and velocity nodes at the center of each grid cell;

[0237] S413. Based on the Taylor expansion formula, the wave field function is expanded into a series at each pressure node, and the pressure finite difference coefficients of each order spatial derivative are solved by fitting the least squares method.

[0238] Based on the Taylor expansion formula, the wave field function is expanded into a series at each velocity node, and the velocity finite difference coefficients of each order spatial derivative are solved by fitting the least squares method.

[0239] S414. Arrange the differential coefficients of each pressure node or velocity node according to the grid index to form a finite pressure differential coefficient submatrix. and finite velocity difference coefficient submatrix ;

[0240] S415. Constructing the finite difference coefficient matrix :

[0241] .

[0242] By dividing the regular surface region into three-dimensional uniform meshes, arranging pressure nodes and velocity nodes in an alternating manner, and solving the difference coefficients based on the Taylor expansion formula and the least squares method to form a finite difference coefficient matrix, efficient discretization of the regular region is achieved, optimizing the utilization of computing resources, ensuring the rapid convergence and numerical stability of spatial derivative calculations, and improving the overall inversion efficiency.

[0243] In some specific embodiments, the steps of using a dynamic coupling discretization method based on staggered node layout to discretize the transition region and obtain the hybrid weight coefficient matrix include:

[0244] S421. Arrange alternating pressure nodes and velocity nodes in the transition zone;

[0245] S422. Using each transition zone pressure node as the central transition zone pressure node, select the closest one. Each transition zone pressure node constructs a pressure hybrid support domain for the central transition zone pressure node;

[0246] Each transition zone velocity node is used as the central transition zone velocity node, and the closest one is selected. Each transition zone velocity node constructs the velocity hybrid support domain of the central transition zone velocity node;

[0247] S423. Calculate the pressure mixing weighting coefficient for each pressure node in the central transition zone and the velocity mixing weighting coefficient for each velocity node in the central transition zone:

[0248]

[0249]

[0250]

[0251]

[0252] in, This indicates that in this pressure-mixed support domain, the first... Pressure mixing weighting coefficient for each mixed pressure node;

[0253] This indicates that in this pressure-mixed support domain, the first... The pressure derivative weighting coefficients of each mixed pressure node are calculated using the staggered node RBF-FD method;

[0254] This indicates that in this pressure-mixed support domain, the first... The finite pressure difference coefficients of the mixed pressure nodes are calculated using the staggered mesh finite difference method;

[0255] This represents the Euclidean distance between the pressure node in the central transition zone of this pressure hybrid support domain and the nearest pressure node in the regular region;

[0256] This represents the Euclidean distance between the pressure node in the central transition zone of this pressure hybrid support domain and the nearest pressure node in the irregular region.

[0257] This indicates that in this velocity hybrid support domain, the first... The velocity mixing weight coefficient of each mixed velocity node;

[0258] This indicates that in this velocity hybrid support domain, the first... The velocity derivative weighting coefficients of each mixed velocity node are calculated using the staggered node RBF-FD method;

[0259] This indicates that in this velocity hybrid support domain, the first... The finite velocity difference coefficients of the mixed velocity nodes are calculated using the staggered mesh finite difference method;

[0260] This represents the Euclidean distance between the velocity node in the central transition zone of this velocity hybrid support domain and the nearest regular region velocity node;

[0261] This represents the Euclidean distance between the velocity node in the central transition zone of this velocity hybrid support domain and the nearest irregular region velocity node;

[0262] S424. Summarize the pressure mixing weight coefficients of all pressure nodes in the central transition zone to form a pressure mixing weight coefficient submatrix. :

[0263] The row index corresponds to the global number of the pressure node in the central transition zone;

[0264] The column index corresponds to the global number of all pressure nodes in the entire transition zone;

[0265] The element value is the pressure mixing weight coefficient between the pressure node in the central transition zone corresponding to the row index and the pressure node in the transition zone corresponding to the column index. For each row, a non-zero value is stored only in the column position corresponding to the pressure node in its support domain.

[0266] The velocity derivative weight coefficients of all central velocity nodes are summarized to form a velocity weight submatrix. :

[0267] The row index corresponds to the global number of the central velocity node;

[0268] The column index corresponds to the global number of all velocity nodes in the entire irregular region;

[0269] The element value is the velocity derivative weight coefficient between the center velocity node corresponding to the row index and the velocity node corresponding to the column index. For each row, a non-zero value is stored only in the column position corresponding to the velocity node within its supporting domain.

[0270] S425. Constructing the finite difference coefficient matrix :

[0271] .

[0272] By arranging staggered nodes in the transition zone, constructing a hybrid support domain, and dynamically calculating the pressure hybrid weight coefficient and velocity hybrid weight coefficient based on distance weight to form a hybrid weight coefficient matrix, smooth coupling discretization of the transition zone is achieved, ensuring the continuity of the solution between regular and irregular regions, avoiding boundary numerical instability, and improving the consistency and accuracy of the overall inversion.

[0273] Step S5: Preset the medium parameter field of the coal seam to be tested.

[0274] By pre-setting the medium parameter field of the coal seam to be tested, the initial parameter conditions for the forward modeling of seismic waves are provided, ensuring the initiation of the inversion process, avoiding simulation interruption due to missing parameters, and laying the foundation for iterative optimization.

[0275] In some specific embodiments, the medium parameter field includes the wave velocity field and the density field.

[0276] By clearly defining the medium parameter fields, including the wave velocity field and the density field, a comprehensive coverage of the physical properties of coal seams is achieved, providing a complete parameter framework for seismic wave forward modeling, ensuring the physical authenticity of the wave field simulation and the comprehensiveness of the inversion results, and avoiding model bias caused by missing parameters.

[0277] Step S6: Based on the source signal, the medium parameter field of the current coal seam to be measured, and the partitioned discrete solution results, perform forward modeling of seismic waves to generate a forward modeling wave field, and extract the forward modeling waveform of the detection point from the forward modeling wave field.

[0278] Based on the partitioned discrete solution results and the current medium parameter field, a forward modeling simulation of seismic waves is performed to generate a forward modeling wavefield and extract the waveforms of the detection points. This realizes the numerical reconstruction of the actual wavefield, provides reference data for residual comparison, and supports parameter update decisions.

[0279] In some specific embodiments, the steps of performing forward seismic wave modeling based on the source signal, the medium parameter field of the current coal seam to be measured, and the results of partitioned discretization to generate the forward modeling wavefield include:

[0280] S601. Construct the governing equations describing the propagation of seismic waves. The governing equations include medium parameter fields, wave field variables, and spatial derivative terms.

[0281] S602. Discretize the control equations to obtain discretized control equations, including:

[0282] In irregular surface regions, the spatial derivative term is calculated using the RBF-FD weighting coefficient matrix;

[0283] In regular surface regions, the spatial derivative term is calculated using a finite difference coefficient matrix;

[0284] For the transition region, the spatial derivative term is calculated using the mixed weight coefficient matrix;

[0285] S603. Based on the source signal and medium parameter field, forward modeling is performed using discretized control equations to output the forward modeling wave field.

[0286] By constructing the seismic wave propagation control equations, calculating the spatial derivative terms in different regions using corresponding weight coefficient matrices based on the partitioned discrete solution results, and executing forward modeling to output the wavefield, an efficient and accurate forward modeling process was achieved. This ensures the comparability of the simulated wavefield with the actual wavefield, provides a reliable benchmark for residual comparison, and supports inversion iterative optimization.

[0287] Step S7: Compare the residuals of the forward modeling waveform at the detection point with the measured seismic waveform, and calculate the residual norm.

[0288] If the residual norm is greater than the preset threshold, the medium parameter field is updated by the optimization algorithm and then the process returns to step S6.

[0289] If the residual norm is less than or equal to the preset threshold, the spatial distribution model of the current medium parameters is output.

[0290] By comparing the residuals of the forward-modeled waveforms at the detection points with the measured seismic waveforms to calculate the residual norm, and updating the medium parameter field or outputting the spatial distribution model based on the results, iterative optimization and convergence control of the parameter field are achieved, ensuring that the inversion results reach the preset accuracy threshold, and finally outputting a reliable spatial distribution model of medium parameters.

[0291] In some specific embodiments, the optimization algorithm used is gradient descent.

[0292] By employing gradient descent as the optimization algorithm to update the medium parameter field, rapid convergence of parameter iteration is achieved, reducing computation time and resource consumption, ensuring efficient and stable inversion process, and avoiding the computational burden caused by complex algorithms.

[0293] In one specific embodiment, the steps of the RBF-FD fusion difference method for seismic wave full waveform inversion include:

[0294] Step S1: Set up an artificial seismic source and detection point in the coal seam to be tested and trigger the seismic source signal. Collect the measured seismic waveform data of the detection point through the seismic detector set at the detection point.

[0295] Step S2 involves deploying a laser scanning device in the tunnel of the coal seam to be measured. The laser scanning device emits pulsed laser light and receives reflected signals to generate point cloud data of the tunnel. The steps include:

[0296] S201. Configuration Each sampling point is targeted by a laser scanning device. Emit a pulsed laser beam and record the vertical direction angle of the pulsed laser beam. and horizontal direction angle ;

[0297] S202. Receive the laser signal reflected from the tunnel surface and record each sampling point. slant distance and reflection intensity ;

[0298] S203. Calculate the coordinate information of each sampling point, including the three-dimensional coordinates of the point. Calculated based on the principle of laser ranging:

[0299]

[0300] It also includes the reflection intensity of the sampling point. ;

[0301] S204. Summarize the coordinate information of all sampled points to generate the original point cloud dataset;

[0302] S205. Denoise the initial point cloud dataset to generate a valid point cloud dataset.

[0303] Step S3: Construct a 3D geometric model of the tunnel based on the point cloud data. Divide the 3D geometric model into regular surface regions and irregular surface regions, and set transition zones at the boundaries between the regular and irregular surface regions. The steps include:

[0304] S301. Based on the effective tunnel surface point cloud data, calculate the value of each sampling point. normal vector ;

[0305] S302. For each sampling point Calculate its neighborhood point set rate of change of normal vector :

[0306] ;

[0307] S303. Set the threshold for the first normal vector change rate. The threshold of the rate of change of the second normal vector ,satisfy ;

[0308] Will and , Compare, among which:

[0309] like Then mark It belongs to an irregular surface area;

[0310] like Then mark It belongs to a regular surface area;

[0311] like Then mark It belongs to the transition zone;

[0312] S304. Perform morphological closing operations on the marking results to output continuously closed irregular surface regions. Regular surface areas and transition zone :

[0313]

[0314]

[0315] .

[0316] Step S4: Perform partitioned discretization to obtain the partitioned discretization results, including:

[0317] In irregular surface regions, the staggered-node RBF-FD method is used for discretization to obtain the RBF-FD weight coefficient matrix of the irregular regions. The steps include:

[0318] S401. Pressure nodes and velocity nodes are arranged alternately in irregular surface areas;

[0319] The density of pressure nodes and velocity nodes is dynamically adjusted according to the local curvature. The higher the local curvature, the higher the density of pressure nodes and velocity nodes.

[0320] S402. For each pressure node as the central pressure node, select the M nearest neighbor pressure nodes to construct the support domain of the central pressure node;

[0321] Each velocity node is taken as the central velocity node, and the N nearest neighbor velocity nodes are selected to construct the support domain of the central velocity node.

[0322] S403. Set the radial basis functions:

[0323]

[0324]

[0325]

[0326]

[0327] in, This represents the Euclidean distance between the p-th pressure node and the central pressure node within the support domain of this central pressure node;

[0328] This represents the position vector of the central pressure node within the support domain of this central pressure node;

[0329] This represents the position vector of the p-th pressure node within the pressure node support domain of this center;

[0330] This represents the Euclidean distance between the q-th velocity node and the central velocity node within the support domain of this central velocity node;

[0331] This represents the position vector of the central velocity node within the support domain of this central velocity node;

[0332] This represents the position vector of the q-th velocity node in the support domain of this central velocity node;

[0333] Represents the calculation of Euclidean norm;

[0334] S404. Setting the pressure approximation function based on radial basis functions and velocity approximation function :

[0335]

[0336]

[0337] in, This represents the pressure RBF weight coefficient of the p-th pressure node in the pressure node support domain of this center;

[0338] This represents the velocity RBF weight coefficient of the q-th velocity node in the velocity node support domain of this center;

[0339] S405. For each central pressure node, the derivative of the pressure approximation function is used to obtain the expression for the spatial pressure derivative:

[0340]

[0341] in, For the p-th pressure node in the support domain of this center pressure node, the derivative weight coefficient is...

[0342] For spatial differential operators;

[0343] For each central velocity node, differentiating the velocity approximation function yields the expression for the spatial velocity derivative:

[0344]

[0345] in, The derivative weight coefficient of the q-th velocity node in the velocity node support domain of this center. ;

[0346] S406. Solve for the spatial pressure derivative expression and the spatial velocity derivative expression to obtain the pressure derivative weighting coefficient between each central pressure node and all pressure nodes in its support domain, and the velocity derivative weighting coefficient between each central velocity node and all velocity nodes in its support domain.

[0347] S407. Summarize the pressure derivative weight coefficients of all central pressure nodes to form a pressure weight submatrix. :

[0348] The row index corresponds to the global number of the central pressure node;

[0349] The column index corresponds to the global number of all pressure nodes in the entire irregular region;

[0350] The element value is the pressure derivative weight coefficient between the center pressure node corresponding to the row index and the pressure node corresponding to the column index. For each row, a non-zero value is stored only in the column position corresponding to the pressure node within its support domain.

[0351] The velocity derivative weight coefficients of all central velocity nodes are summarized to form a velocity weight submatrix. :

[0352] The row index corresponds to the global number of the central velocity node;

[0353] The column index corresponds to the global number of all velocity nodes in the entire irregular region;

[0354] The element value is the velocity derivative weight coefficient between the center velocity node corresponding to the row index and the velocity node corresponding to the column index. For each row, a non-zero value is stored only in the column position corresponding to the velocity node within its supporting domain.

[0355] S408. Construct the RBF-FD weight coefficient matrix :

[0356] ;

[0357] The staggered mesh finite difference method is used to discretize and solve the problem in the regular surface region to obtain the finite difference coefficient matrix of the regular region. The specific steps include:

[0358] S411. Divide the regular surface area into a three-dimensional uniform grid, with the grid step size set according to the preset accuracy requirements;

[0359] S412. Arrange pressure nodes at each grid vertex and velocity nodes at the center of each grid cell;

[0360] S413. Based on the Taylor expansion formula, the wave field function is expanded into a series at each pressure node, and the pressure finite difference coefficients of each order spatial derivative are solved by fitting the least squares method.

[0361] Based on the Taylor expansion formula, the wave field function is expanded into a series at each velocity node, and the velocity finite difference coefficients of each order spatial derivative are solved by fitting the least squares method.

[0362] S414. Arrange the differential coefficients of each pressure node or velocity node according to the grid index to form a finite pressure differential coefficient submatrix. and finite velocity difference coefficient submatrix ;

[0363] S415. Constructing the finite difference coefficient matrix :

[0364] ;

[0365] The transition region is discretized using a dynamic coupling discretization method based on staggered node layout to obtain the mixed weight coefficient matrix. The steps include:

[0366] S421. Arrange alternating pressure nodes and velocity nodes in the transition zone;

[0367] S422. Using each transition zone pressure node as the central transition zone pressure node, select the closest one. Each transition zone pressure node constructs a pressure hybrid support domain for the central transition zone pressure node;

[0368] Each transition zone velocity node is used as the central transition zone velocity node, and the closest one is selected. Each transition zone velocity node constructs the velocity hybrid support domain of the central transition zone velocity node;

[0369] S423. Calculate the pressure mixing weighting coefficient for each pressure node in the central transition zone and the velocity mixing weighting coefficient for each velocity node in the central transition zone:

[0370]

[0371]

[0372]

[0373]

[0374] in, This indicates that in this pressure-mixed support domain, the first... Pressure mixing weighting coefficient for each mixed pressure node;

[0375] This indicates that in this pressure-mixed support domain, the first... The pressure derivative weighting coefficients of each mixed pressure node are calculated using the staggered node RBF-FD method;

[0376] This indicates that in this pressure-mixed support domain, the first... The finite pressure difference coefficients of the mixed pressure nodes are calculated using the staggered mesh finite difference method;

[0377] This represents the Euclidean distance between the pressure node in the central transition zone of this pressure hybrid support domain and the nearest pressure node in the regular region;

[0378] This represents the Euclidean distance between the pressure node in the central transition zone of this pressure hybrid support domain and the nearest pressure node in the irregular region.

[0379] This indicates that in this velocity hybrid support domain, the first... The velocity mixing weight coefficient of each mixed velocity node;

[0380] This indicates that in this velocity hybrid support domain, the first... The velocity derivative weighting coefficients of each mixed velocity node are calculated using the staggered node RBF-FD method;

[0381] This indicates that in this velocity hybrid support domain, the first... The finite velocity difference coefficients of the mixed velocity nodes are calculated using the staggered mesh finite difference method;

[0382] This represents the Euclidean distance between the velocity node in the central transition zone of this velocity hybrid support domain and the nearest regular region velocity node;

[0383] This represents the Euclidean distance between the velocity node in the central transition zone of this velocity hybrid support domain and the nearest irregular region velocity node;

[0384] S424. Summarize the pressure mixing weight coefficients of all pressure nodes in the central transition zone to form a pressure mixing weight coefficient submatrix. :

[0385] The row index corresponds to the global number of the pressure node in the central transition zone;

[0386] The column index corresponds to the global number of all pressure nodes in the entire transition zone;

[0387] The element value is the pressure mixing weight coefficient between the pressure node in the central transition zone corresponding to the row index and the pressure node in the transition zone corresponding to the column index. For each row, a non-zero value is stored only in the column position corresponding to the pressure node in its support domain.

[0388] The velocity derivative weight coefficients of all central velocity nodes are summarized to form a velocity weight submatrix. :

[0389] The row index corresponds to the global number of the central velocity node;

[0390] The column index corresponds to the global number of all velocity nodes in the entire irregular region;

[0391] The element value is the velocity derivative weight coefficient between the center velocity node corresponding to the row index and the velocity node corresponding to the column index. For each row, a non-zero value is stored only in the column position corresponding to the velocity node within its supporting domain.

[0392] S425. Constructing the finite difference coefficient matrix :

[0393] .

[0394] Step S5: Preset the medium parameter field of the coal seam to be tested, including wave velocity field and density field.

[0395] Step S6: Based on the source signal, the medium parameter field of the current coal seam to be measured, and the results of the partitioned discretization solution, perform forward seismic wave modeling to generate the forward modeling wave field. The steps include:

[0396] S601. Construct the governing equations describing the propagation of seismic waves. The governing equations include medium parameter fields, wave field variables, and spatial derivative terms.

[0397] S602. Discretize the control equations to obtain discretized control equations, including:

[0398] In irregular surface regions, the spatial derivative term is calculated using the RBF-FD weighting coefficient matrix;

[0399] In regular surface regions, the spatial derivative term is calculated using a finite difference coefficient matrix;

[0400] For the transition region, the spatial derivative term is calculated using the mixed weight coefficient matrix;

[0401] S603. Based on the source signal and medium parameter field, forward modeling is performed using discretized control equations to output the forward modeling wave field;

[0402] Extract the forward-modeled waveform of the detection point in the forward-modeled wave field.

[0403] Step S7: Compare the residuals of the forward modeling waveform at the detection point with the measured seismic waveform, and calculate the residual norm.

[0404] If the residual norm is greater than the preset threshold, the medium parameter field is updated by gradient descent and then the process returns to step S6.

[0405] If the residual norm is less than or equal to the preset threshold, the spatial distribution model of the current medium parameters is output.

[0406] The following are embodiments of the RBF-FD fusion differential method seismic wave full waveform inversion system provided in this application. This RBF-FD fusion differential method seismic wave full waveform inversion system belongs to the same inventive concept as the seismic wave full waveform inversion methods in the above embodiments. For details not described in detail in the embodiments of the seismic wave full waveform inversion system, please refer to the embodiments of the RBF-FD fusion differential method seismic wave full waveform inversion method described above.

[0407] like Figure 2 As shown, the RBF-FD fusion difference method seismic wave full waveform inversion system includes:

[0408] The point cloud data acquisition and processing module is used to deploy laser scanning equipment in the tunnel of the coal seam to be measured. The laser scanning equipment emits pulsed laser and receives reflected signals to generate point cloud data of the tunnel.

[0409] The seismic data acquisition module is used to set up artificial seismic sources and detection points in the coal seam to be measured and trigger the seismic source signal. The measured seismic waveform data of the detection points is acquired by the seismic detectors set at the detection points.

[0410] The 3D model construction and region division module is used to construct a 3D geometric model of the tunnel based on point cloud data, divide the 3D geometric model of the tunnel into regular surface regions and irregular surface regions, and set a transition zone at the boundary between regular surface regions and irregular surface regions.

[0411] The partitioned discretization module is used to perform partitioned discretization and obtain the partitioned discretization results, including:

[0412] In irregular surface regions, the staggered node RBF-FD method is used for discretization to obtain the RBF-FD weight coefficient matrix of the irregular regions.

[0413] The staggered mesh finite difference method is used to discretize and solve the problem in the regular surface region to obtain the finite difference coefficient matrix of the regular region.

[0414] The transition zone is solved by a dynamic coupling discretization method based on staggered node layout to obtain the mixed weight coefficient matrix;

[0415] The medium parameter field initialization module is used to preset the medium parameter field of the coal seam to be tested.

[0416] The seismic wave forward modeling module is used to perform seismic wave forward modeling based on the source signal, the medium parameter field of the current coal seam to be measured, and the partitioned discrete solution results, generate the forward modeling wave field, and extract the forward modeling waveform of the detection point in the forward modeling wave field;

[0417] The residual calculation and parameter optimization module is used to compare the residuals of the forward modeling waveforms and the measured seismic waveforms at the detection points, calculate the residual norm, and update the medium parameter field through optimization algorithms until the accuracy requirements are met.

[0418] The seismic wave full waveform inversion system in this embodiment is used to realize the seismic wave full waveform inversion method of RBF-FD fusion difference method.

[0419] This application also provides an electronic device for implementing the various embodiments of this application. Figure 3 To illustrate the hardware structure of an electronic device according to various embodiments of this application, as shown in the following diagram... Figure 3 As shown, the electronic device includes a memory, a processor, and a computer program stored in the memory and capable of running on the processor.

[0420] Those skilled in the art will understand that the electronic device structure involved in the embodiments of this application does not constitute a limitation on the electronic device. The electronic device may include more or fewer components than shown in the figure, or combine certain components, or have different component arrangements.

[0421] In embodiments of this application, electronic devices include, but are not limited to, laptop computers, desktop computers, workstations, personal digital assistants, servers, blade servers, mainframe computers, and other suitable computers. Electronic devices may also represent various forms of mobile devices and other similar computing devices. The components shown herein, their connections and relationships, and their functions are merely examples and are not intended to limit the implementation of the embodiments of this application described and / or claimed herein.

[0422] In this application embodiment, the processor can be implemented using at least one of an Application-Specific Integrated Circuit (ASIC), a Digital Signal Processor (DSP), a Digital Signal Processing Device (DSPD), a processor, a controller, a microcontroller, a microprocessor, or an electronic unit designed to perform the functions described herein. In some cases, such implementations can be implemented within a controller. For software implementations, implementations such as processes or functions can be implemented with separate software modules that allow the performance of at least one function or operation. The software code can be implemented by a software application (or program) written in any suitable programming language, and the software code can be stored in memory and executed by the controller.

[0423] In addition, the electronic device includes some functional modules not shown, which will not be described in detail here.

[0424] Those skilled in the art will understand that the various aspects of the electronic device provided in this application can be implemented as a system, method, or program product. Therefore, the various aspects of this application can be specifically implemented in the following forms: a completely hardware implementation, a completely software implementation (including firmware, microcode, etc.), or a combination of hardware and software aspects, collectively referred to herein as a "circuit," "module," or "system."

[0425] This application also provides a storage medium storing a program product capable of implementing the RBF-FD fusion differential method for seismic wave full waveform inversion. In some possible embodiments, various aspects of this application can also be implemented as a program product comprising program code that, when run on a terminal device, causes the terminal device to perform the steps described in the foregoing "Exemplary Methods" section of this specification according to the various exemplary embodiments of this application.

[0426] The storage medium may be any combination of one or more readable media. A readable medium may be a readable signal medium or a readable storage medium. A readable storage medium may be, for example,, but not limited to, an electrical, magnetic, optical, electromagnetic, infrared, or semiconductor system, apparatus, or device, or any combination thereof. More specific examples (a non-exhaustive list) of readable storage media include: electrical connections having one or more wires, portable disks, hard disks, random access memory (RAM), read-only memory (ROM), erasable programmable read-only memory (EPROM or flash memory), optical fiber, portable compact disk read-only memory (CD-ROM), optical storage devices, magnetic storage devices, or any suitable combination thereof.

[0427] The above description of the disclosed embodiments enables those skilled in the art to make or use this application. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of this application. Therefore, this application is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A method for seismic wave full waveform inversion using the RBF-FD fusion differential method, characterized in that, include: S1. Set up artificial seismic sources and detection points in the coal seam to be tested and trigger the seismic source signal. Collect the measured seismic waveform data of the detection points through the seismic detectors set at the detection points; S2. A laser scanning device is deployed in the tunnel of the coal seam to be measured. The laser scanning device emits pulsed lasers and receives reflected signals to generate point cloud data of the tunnel. S3. Construct a three-dimensional geometric model of the tunnel based on point cloud data, divide the tunnel into regular surface regions and irregular surface regions in the three-dimensional geometric model, and set a transition zone at the junction of the regular surface regions and irregular surface regions; S4. Perform partitioned discretization to obtain the partitioned discretization results, including: In irregular surface regions, the staggered node RBF-FD method is used for discretization to obtain the RBF-FD weight coefficient matrix of the irregular regions. The staggered mesh finite difference method is used to discretize and solve the problem in the regular surface region to obtain the finite difference coefficient matrix of the regular region. The transition zone is solved by a dynamic coupling discretization method based on staggered node layout to obtain the mixed weight coefficient matrix; S5. Preset the medium parameter field of the coal seam to be tested; S6. Based on the source signal, the medium parameter field of the current coal seam to be measured, and the partitioned discrete solution results, perform forward modeling of seismic waves to generate a forward modeling wave field, and extract the forward modeling waveform of the detection point from the forward modeling wave field; S7. Compare the residuals of the forward modeling waveforms at the detection points with the measured seismic waveforms, and calculate the residual norm; If the residual norm is greater than the preset threshold, the medium parameter field is updated by the optimization algorithm and then the process returns to step S6. If the residual norm is less than or equal to the preset threshold, the spatial distribution model of the current medium parameters is output.

2. The seismic wave full waveform inversion method as described in claim 1, characterized in that, In step S2, the laser scanning device emits pulsed laser light and receives reflected signals to generate point cloud data of the tunnel. S201. Configuration Each sampling point is targeted by a laser scanning device. Emit a pulsed laser beam and record the vertical direction angle of the pulsed laser beam. and horizontal direction angle ; S202. Receive the laser signal reflected from the tunnel surface and record each sampling point. slant distance ; S203. Calculate the coordinate information of each sampling point, including the three-dimensional coordinates of the point. Calculated based on the principle of laser ranging: S204. Summarize the coordinate information of all sampled points to generate the original point cloud dataset; S205. Denoise the initial point cloud dataset to generate a valid point cloud dataset.

3. The seismic wave full waveform inversion method as described in claim 1, characterized in that, Step S3, which involves dividing the tunnel's three-dimensional geometric model into regular and irregular surface regions and setting a transition zone at the boundary between these regions, includes the following steps: S301. Based on the effective tunnel surface point cloud data, calculate the value of each sampling point. normal vector ; S302. For each sampling point Calculate its neighborhood point set rate of change of normal vector : ； S303. Set the threshold for the first normal vector change rate. The threshold of the rate of change of the second normal vector ,satisfy ; Will and , Compare, among which: like Then mark It belongs to an irregular surface area; like Then mark It belongs to a regular surface area; like Then mark It belongs to the transition zone; S304. Perform morphological closing operations on the marking results to output continuously closed irregular surface regions. Regular surface areas and transition zone : 。 4. The seismic wave full waveform inversion method as described in claim 1, characterized in that, Step S4 involves using the staggered-node RBF-FD method to discretize and solve the irregular surface region, obtaining the RBF-FD weight coefficient matrix for the irregular region. The specific steps include: S401. Pressure nodes and velocity nodes are arranged alternately in irregular surface areas; The density of pressure nodes and velocity nodes is dynamically adjusted according to the local curvature. The higher the local curvature, the higher the density of pressure nodes and velocity nodes. S402. For each pressure node as the central pressure node, select the M nearest neighbor pressure nodes to construct the support domain of the central pressure node; Each velocity node is taken as the central velocity node, and the N nearest neighbor velocity nodes are selected to construct the support domain of the central velocity node. S403. Set the radial basis functions: in, This represents the Euclidean distance between the p-th pressure node and the central pressure node within the support domain of this central pressure node; This represents the position vector of the central pressure node within the support domain of this central pressure node; This represents the position vector of the p-th pressure node within the pressure node support domain of this center; This represents the Euclidean distance between the q-th velocity node and the central velocity node within the support domain of this central velocity node; This represents the position vector of the central velocity node within the support domain of this central velocity node; This represents the position vector of the q-th velocity node in the support domain of this central velocity node; Represents the calculation of Euclidean norm; S404. Setting the pressure approximation function based on radial basis functions and velocity approximation function : in, This represents the pressure RBF weight coefficient of the p-th pressure node in the pressure node support domain of this center; This represents the velocity RBF weight coefficient of the q-th velocity node in the velocity node support domain of this center; S405. For each central pressure node, the derivative of the pressure approximation function is used to obtain the expression for the spatial pressure derivative: in, For the p-th pressure node in the support domain of this center pressure node, the derivative weight coefficient is... For spatial differential operators; For each central velocity node, the derivative of the velocity approximation function is used to obtain the expression for the spatial velocity derivative: in, The derivative weight coefficient of the q-th velocity node in the velocity node support domain of this center. ; S406. Solve for the spatial pressure derivative expression and the spatial velocity derivative expression to obtain the pressure derivative weighting coefficient between each central pressure node and all pressure nodes in its support domain, and the velocity derivative weighting coefficient between each central velocity node and all velocity nodes in its support domain. S407. Summarize the pressure derivative weight coefficients of all central pressure nodes to form a pressure weight submatrix. : The row index corresponds to the global number of the central pressure node; The column index corresponds to the global number of all pressure nodes in the entire irregular region; The element value is the pressure derivative weight coefficient between the center pressure node corresponding to the row index and the pressure node corresponding to the column index. For each row, a non-zero value is stored only in the column position corresponding to the pressure node within its support domain. The velocity derivative weight coefficients of all central velocity nodes are summarized to form a velocity weight submatrix. : The row index corresponds to the global number of the central velocity node; The column index corresponds to the global number of all velocity nodes in the entire irregular region; The element value is the velocity derivative weight coefficient between the center velocity node corresponding to the row index and the velocity node corresponding to the column index. For each row, a non-zero value is stored only in the column position corresponding to the velocity node within its supporting domain. S408. Construct the RBF-FD weight coefficient matrix : 。 5. The seismic wave full waveform inversion method as described in claim 1, characterized in that, In step S4, the specific steps for discretizing and solving the finite difference coefficient matrix of the regular surface region using the staggered mesh finite difference method include: S411. Divide the regular surface area into a three-dimensional uniform grid, with the grid step size set according to the preset accuracy requirements; S412. Arrange pressure nodes at each grid vertex and velocity nodes at the center of each grid cell; S413. Based on the Taylor expansion formula, the wave field function is expanded into a series at each pressure node, and the pressure finite difference coefficients of each order spatial derivative are solved by fitting the least squares method. Based on the Taylor expansion formula, the wave field function is expanded into a series at each velocity node, and the velocity finite difference coefficients of each order spatial derivative are solved by fitting the least squares method. S414. Arrange the differential coefficients of each pressure node or velocity node according to the grid index to form a finite pressure differential coefficient submatrix. and finite velocity difference coefficient submatrix ; S415. Constructing the finite difference coefficient matrix : 。 6. The seismic wave full waveform inversion method as described in claim 1, characterized in that, Step S4, which involves using a dynamic coupling discretization method based on staggered node layout to discretize the transition region and obtain the mixed weight coefficient matrix, includes the following steps: S421. Arrange alternating pressure nodes and velocity nodes in the transition zone; S422. Using each transition zone pressure node as the central transition zone pressure node, select the closest one. Each transition zone pressure node constructs a pressure hybrid support domain for the central transition zone pressure node; Each transition zone velocity node is used as the central transition zone velocity node, and the closest one is selected. Each transition zone velocity node constructs the velocity hybrid support domain of the central transition zone velocity node; S423. Calculate the pressure mixing weighting coefficient for each pressure node in the central transition zone and the velocity mixing weighting coefficient for each velocity node in the central transition zone: in, This indicates that in this pressure-mixed support domain, the first... Pressure mixing weighting coefficient for each mixed pressure node; This indicates that in this pressure-mixed support domain, the first... The pressure derivative weighting coefficients of each mixed pressure node are calculated using the staggered node RBF-FD method; This indicates that in this pressure-mixed support domain, the first... The finite pressure difference coefficients of the mixed pressure nodes are calculated using the staggered mesh finite difference method; This represents the Euclidean distance between the pressure node in the central transition zone of this pressure hybrid support domain and the nearest pressure node in the regular region; This represents the Euclidean distance between the pressure node in the central transition zone of this pressure hybrid support domain and the nearest pressure node in the irregular region. This indicates that in this velocity hybrid support domain, the first... The velocity mixing weight coefficient of each mixed velocity node; This indicates that in this velocity hybrid support domain, the first... The velocity derivative weighting coefficients of each mixed velocity node are calculated using the staggered node RBF-FD method; This indicates that in this velocity hybrid support domain, the first... The finite velocity difference coefficients of the mixed velocity nodes are calculated using the staggered mesh finite difference method; This represents the Euclidean distance between the velocity node in the central transition zone of this velocity hybrid support domain and the nearest regular region velocity node; This represents the Euclidean distance between the velocity node in the central transition zone of this velocity hybrid support domain and the nearest irregular region velocity node; S424. Summarize the pressure mixing weight coefficients of all pressure nodes in the central transition zone to form a pressure mixing weight coefficient submatrix. : The row index corresponds to the global number of the pressure node in the central transition zone; The column index corresponds to the global number of all pressure nodes in the entire transition zone; The element value is the pressure mixing weight coefficient between the pressure node in the central transition zone corresponding to the row index and the pressure node in the transition zone corresponding to the column index. For each row, a non-zero value is stored only in the column position corresponding to the pressure node in its support domain. The velocity derivative weight coefficients of all central velocity nodes are summarized to form a velocity weight submatrix. : The row index corresponds to the global number of the central velocity node; The column index corresponds to the global number of all velocity nodes in the entire irregular region; The element value is the velocity derivative weight coefficient between the center velocity node corresponding to the row index and the velocity node corresponding to the column index. For each row, a non-zero value is stored only in the column position corresponding to the velocity node within its supporting domain. S425. Constructing the finite difference coefficient matrix : 。 7. The seismic wave full waveform inversion method as described in claim 1, characterized in that, The medium parameter field includes the wave velocity field and the density field.

8. The seismic wave full waveform inversion method as described in claim 1, characterized in that, In step S6, the steps of performing forward seismic wave modeling based on the source signal, the medium parameter field of the current coal seam to be measured, and the partitioned discrete solution results to generate the forward modeling wavefield include: S601. Construct the governing equations describing the propagation of seismic waves. The governing equations include medium parameter fields, wave field variables, and spatial derivative terms. S602. Discretize the control equations to obtain discretized control equations, including: In irregular surface regions, the spatial derivative term is calculated using the RBF-FD weighting coefficient matrix; In regular surface regions, the spatial derivative term is calculated using a finite difference coefficient matrix; For the transition region, the spatial derivative term is calculated using the mixed weight coefficient matrix; S603. Based on the source signal and medium parameter field, forward modeling is performed using discretized control equations to output the forward modeling wave field.

9. The seismic wave full waveform inversion method as described in claim 1, characterized in that, The optimization algorithm used in step S6 is gradient descent.

10. A seismic wave full waveform inversion system based on the RBF-FD fusion differential method, characterized in that, To implement the seismic wave full waveform inversion method as described in any one of claims 1-9, the method includes: The seismic data acquisition module is used to set up artificial seismic sources and detection points in the coal seam to be measured and trigger the seismic source signal. The measured seismic waveform data of the detection points is acquired by the seismic detectors set at the detection points. The point cloud data acquisition and processing module is used to deploy laser scanning equipment in the tunnel of the coal seam to be measured. The laser scanning equipment emits pulsed laser and receives reflected signals to generate point cloud data of the tunnel. The 3D model construction and region division module is used to construct a 3D geometric model of the tunnel based on point cloud data, divide the 3D geometric model of the tunnel into regular surface regions and irregular surface regions, and set a transition zone at the boundary between regular surface regions and irregular surface regions. The partitioned discretization module is used to perform partitioned discretization and obtain the partitioned discretization results, including: In irregular surface regions, the staggered node RBF-FD method is used for discretization to obtain the RBF-FD weight coefficient matrix of the irregular regions. The staggered mesh finite difference method is used to discretize and solve the regular surface region to obtain the finite difference coefficient matrix of the regular region; The transition zone is solved by a dynamic coupling discretization method based on staggered node layout, and the mixed weight coefficient matrix is ​​obtained. The medium parameter field initialization module is used to preset the medium parameter field of the coal seam to be tested. The seismic wave forward modeling module is used to perform seismic wave forward modeling based on the source signal, the medium parameter field of the current coal seam to be measured, and the partitioned discrete solution results, generate the forward modeling wave field, and extract the forward modeling waveform of the detection point in the forward modeling wave field; The residual calculation and parameter optimization module is used to compare the residuals of the forward modeling waveforms and the measured seismic waveforms at the detection points, calculate the residual norm, and update the medium parameter field through optimization algorithms until the accuracy requirements are met.

Images