Fresnel volume three-dimensional tomography method with priori knowledge correction

Abstract

The invention discloses a Fresnel volume three-dimensional tomography method with priori knowledge correction. The method comprises the following steps: preprocessing experimental data to obtain a test time propagation table; then building a grid model, and writing the initial slowness field and coordinates of the transmitting point and the receiving point; performing tomographic forward modeling, and calculating tomographic forward modeling travel time, ray paths, Fresnel body ranges and normalized weight factors of each transmitting-receiving pair; and comparing a test result with a forward modeling result, if an error is greater than a limit value, entering an inversion iteration step, carrying out first correction on the slowness field by adopting an SIRT algorithm and a back projection algorithm, carrying out second correction on the slowness field according to a priori slowness field shape parameter, and carrying out forward modeling calculation on the obtained slowness field again until the error is less than the limit value. According to the method, tomography inversion of the three-dimensional space shear wave velocity is realized, the space coverage rate of rays is improved based on the Fresnel volume method, and the tomography result is corrected in combination with prior knowledge, so that the inversion result is more reliable.

Description

Technical Field

[0001] The present invention belongs to the field of geotechnical engineering tests, and particularly relates to a three-dimensional tomographic imaging method of a Fresnel body with prior knowledge correction. Background Technique

[0002] The shear wave velocity of soil has important applications in aspects such as the classification of research site categories, the discrimination of sand liquefaction, and the conversion of soil dynamic shear modulus, and is an important parameter for characterizing the deformation and dynamic characteristics of soil. The bender element method is a commonly used method for measuring small-strain wave velocity. Its principle is that the excitation element in the bender element sensor vibrates the soil, the obtained shear wave propagates in the soil, reaches the receiving element on the opposite side, and the shear wave propagation time is obtained. For homogeneous soil, it can be converted only through the distance between the bender element and the receiving element, and the shear wave velocity is converted according to its propagation time. However, for large-scale geotechnical model tests, since the test target involves the test soil with a large spatial scale, the shear wave velocity distribution of the test soil has spatial heterogeneity after specific loading. At this time, it is necessary to use multiple shear waves emitted by multiple pairs of bender elements to perform three-dimensional scanning on the soil space, and use tomographic imaging technology to invert the shear wave velocity in the space. Therefore, a reasonable three-dimensional tomographic inversion method is required to efficiently and accurately invert the three-dimensional distribution information of the soil shear wave velocity based on the known shear wave propagation, end point, and propagation time.

[0003] The currently commonly used tomographic inversion method is a two-step solution method, that is, first, the shortest propagation path and the minimum travel time of the ray under the current assumed slowness field are obtained through forward calculation, and then the assumed slowness field is corrected through inversion iteration to make the simulated minimum travel time close to the experimental minimum travel time. Currently, most technologies use ordinary rays as the basis for inversion iteration to correct the slowness field of the grid on the ray path, and a certain difference method is used for correction in the area outside the ray, but this method has defects such as a small number of rays and a small coverage area, and it is easy to cause spatial slowness mutations and anomalies in a small part of the area during multiple inversion iterations. The Fresnel body ray tracing technology, on the basis of ray tracing, combines frequency information and considers the influence of the soil slowness field at the position of the first Fresnel zone on ray propagation, improves the spatial coverage rate of rays during inversion iteration, and reduces the generation of slowness anomaly areas during inversion iteration. It is a relatively advanced tomographic inversion method.

[0004] In actual operation, setting up too many bending element groups in the bending element test incurs high costs, and unexpected situations such as sensor damage are likely to occur during the test. Therefore, the obtained rays may be relatively sparse. This sparsity of ray information may lead to inaccurate inversion of the obtained field during the inversion iteration process, or even completely deviate from the expected physical laws. Therefore, it is necessary to estimate the inversion slowness field model based on prior knowledge of physical laws to make it satisfy certain basic laws, and then perform inversion iteration within the framework of prior knowledge, so as to obtain a satisfactory slowness field result under sparse ray information.

[0005] In summary, at present, the three-dimensional shear wave velocity tomography technologies at home and abroad have the following problems: (1) Most of the existing technologies use ordinary rays to invert the spatial slowness field, covering a small space and prone to spatial slowness mutations and anomalies in some areas; (2) The ray sparsity caused by test costs and complex condition limitations is not considered, resulting in the inversion result not conforming to the basic physical laws; (3) There is less research and application on tomography technologies in three-dimensional space. Summary of the Invention

[0006] The purpose of the present invention is to provide a Fresnel body three-dimensional tomography method with prior knowledge correction in view of the deficiencies of the existing technology. The present invention realizes the tomographic inversion of the three-dimensional shear wave velocity, improves the spatial coverage rate of rays based on the Fresnel body method, and corrects the tomographic result in combination with prior knowledge to make the inversion result more reliable.

[0007] To achieve the above purpose, a Fresnel body three-dimensional tomography method with prior knowledge correction provided by the present invention comprises the following specific steps:

[0008] 1) Pretreatment of test data: According to the waveform table scanned by the bending element array (2), after data pretreatment, an experimental time propagation table is obtained;

[0009] 2) Construction of the grid model: Create a computational grid model, write the assumed initial slowness field at the grid nodes, set the corresponding emission points m and reception points n on the grid nodes according to the positions of the excitation element m and reception element n in the test, and record the serial numbers, spatial coordinates, and node numbers of the emission points and reception points;

[0010] 3) Tomographic forward modeling: For the emission point m, calculate the minimum travel time from it to each grid node in the entire grid model, and then for each reception point n, sequentially read the tomographic forward modeling travel time of each emission-reception pair according to its node number to obtain a forward propagation time table; Subsequently, starting from the reception point n, back-calculate the ray paths of each emission-reception pair according to the grid minimum travel time, and finally calculate the Fresnel body range and normalized weight factor of each emission-reception pair to obtain a forward propagation path table;

[0011] 4) Inversion iteration: Compare the experimental propagation schedule and the forward propagation schedule. If the error between the two is greater than the limit value, enter the inversion iteration step. First, use the SIRT algorithm to calculate the slowness field correction value for the area passed by the Fresnel volume. Subsequently, use the backprojection algorithm to calculate the slowness field correction value for the area not passed by the Fresnel volume. Finally, obtain the corresponding slowness correction value for each grid node. Calculate the iteration factor according to the ray path and the slowness field correction value to update the slowness field 1. Finally, perform a scanning and reorganization calculation on the slowness field 1 according to the prior slowness field shape parameter to obtain the slowness field correction value generated by the prior knowledge correction. Consider the slowness field correction values before and after modifying the prior knowledge to obtain the final slowness field 2;

[0012] 5) Re-perform tomographic forward modeling using the slowness field 2. Compare the experimental propagation schedule and the new forward propagation schedule. If the error between the two is less than the limit value, the calculation ends. If it is still greater than the limit value, repeat step 4) until the error is less than the limit value.

[0013] Further, in step 1), the flexural element array (2) includes a seismic isolation pad (201), a flexural element bracket (202), an excitation element (203), and a receiving element (204). The flexural element array (2) is installed in the model box (1). Among them, several pairs of excitation elements (203) are arranged on the left flexural element bracket (202), and the corresponding several pairs of receiving elements (204) are arranged on the right flexural element bracket (202). The excitation element (203) emits a single vibration under the action of the excitation voltage and generates shear waves in the soil body (3). After propagating through the soil body (3), it reaches the receiving element (204) on the opposite side, thereby obtaining the received wave signal; a seismic isolation pad (201) is arranged below the flexural element bracket to prevent waves from propagating through the box body of the model box (1).

[0014] Further, in step 1), the flexural element array (2) tests the received wave waveforms obtained by combining each excitation element (203) and receiving element (204), and thus organizes and obtains the flexural element array scanning waveform table; after preprocessing each received wave by filtering, propagation time discrimination, and effective data screening, the flexural element array scanning waveform table is converted into an experimental propagation schedule, which contains the serial numbers of the excitation elements (203), the serial numbers of the receiving elements (204), and the ray propagation time obtained by preprocessing the waveforms.

[0015] Further, in step 2), the calculation grid is a cubic grid; the initial slowness field selects a slowness field that linearly decreases along the depth Z.

[0016] Further, in step 3), when back-calculating the ray paths of each transmitter-receiver pair according to the minimum travel time of the grid, it is necessary to search for the travel times of the grid points in different directions in the three layers of grids with R = 1, R = 2, and R = 3 around the starting point, calculate their travel time gradients, and find the direction with the largest descending gradient as the ray propagation path until the emission point is retrieved.

[0017] Further, in step 3), according to the formula the Fresnel volume range is determined, where S represents the emission point, R represents the reception point, P represents any point within the spatial grid, the subscript i means the i-th emission-reception pair, j represents the j-th node, and T SP is the travel time from the emission point to any point; T PR is the travel time from any point to the reception point; T SR is the travel time from the emission point to the reception point, and f is the ray frequency (unit s -1 ). Subsequently, according to the formula the weight factor ω of the influence degree of different parts of the Fresnel volume on the travel time is calculated ij . The closer to the central high-frequency part of the ray, the greater the influence weight, and the closer the weight factor is to 1; then the normalized weight factor Q ij is calculated. For an emission-reception pair i, the sum of the normalized weight factors Q ij is 1.

[0018] Further, in step 4), according to the formula it is judged whether the error between the trial propagation schedule and the forward propagation schedule is less than the limit value, where M is the total number of emission-reception pairs, ΔT is the average relative error between the trial propagation schedule and the forward propagation schedule; ΔT lim is the average relative error limit value between the trial propagation schedule and the forward propagation schedule; T i is the forward propagation time of the i-th emission-reception pair; is the trial propagation time of the i-th emission-reception team.

[0019] Further, in step 4), for the nodes through which the Fresnel volume passes, the simultaneous iterative reconstruction technique SIRT is adopted, and according to the formula the slowness correction value of the Fresnel volume formed by each emission-reception pair i for the node j is calculated, and is averaged by . N is the number of Fresnel volumes passing through this node, and s j represents the initial slowness of the j-th node. For the nodes through which no Fresnel volume passes, the formula is used to perform inverse distance interpolation on the correction values in the correction area where the Fresnel volume passes, and the slowness correction value of the node through which no Fresnel volume passes is obtained. L represents the total number of nodes within the Fresnel volume, d j represents the distance between this point and the j-th node within the Fresnel volume, represents the slowness correction amount of the j-th node within the Fresnel volume.

[0020] Further, in step 4), according to the forward ray path slowness correction value and calculate the iteration factor Enable the slowness inversion iteration to converge faster; and use to update the slowness field 1, where k represents the k-th iteration; λ k represents the iteration factor of the k-th iteration; represents the forward travel time of the i-th transmitter-receiver pair in the k-th iteration; represents the ray path length at the j-th node of the i-th transmitter-receiver pair in the k-th iteration; N i represents the number of nodes that the ray passes through for the i-th transmitter-receiver pair.

[0021] Furthermore, in step 4), a reconstruction area grid range of k*k*k is delimited in the updated slowness field 1. Within the reconstruction area, the slowness field 1 is proportionally adjusted using the prior slowness field shape parameters obtained from prior knowledge. Subsequently, switch to the adjacent next reconstruction area and repeat the proportional adjustment operation until all grid areas are scanned. After such operations, the shape of the slowness field obtained will be strongly similar to the prior slowness field shape parameters obtained from prior knowledge; the final slowness field correction value is calculated by comprehensively considering the slowness field correction value generated by the simultaneous iterative reconstruction technique SIRT and the slowness field correction value generated by prior knowledge correction according to the weight factor α, and the slowness field 2 is updated.

[0022] The beneficial effects of the present invention are: Through the Fresnel body three-dimensional tomography method with prior knowledge correction provided by the present invention, the inversion imaging of the spatial shear wave velocity is realized under the condition of knowing the coordinates of the excitation and reception points and the shear wave propagation time in the three-dimensional space, and the change of the shear wave velocity in the three-dimensional space after the test loading is obtained. Compared with the existing methods, there are mainly the following advantages:

[0023] (1) First, a tomography method based on the Fresnel body is adopted. Considering the frequency information, the influence of the soil slowness field at the position of the first Fresnel zone on the ray propagation is considered, which improves the spatial coverage rate of the rays during the inversion iteration. Compared with the ordinary ray method, the defect of the small coverage area of the ordinary ray method is effectively avoided, and after multiple inversion iterations, spatial slowness mutations and anomalies in some areas (usually the areas passed by ordinary rays) will not occur due to the small correction area.

[0024] (2) Second, prior knowledge is introduced as a correction in the inversion process. In this case, even if the input experimental ray data is relatively sparse in the three-dimensional space, because the physical laws in the prior knowledge are considered during the inversion process, the inversion result will not completely violate the expected physical laws. And the weight factor α is introduced, so that the obtained slowness field comprehensively considers the effects of both tomographic inversion and prior knowledge correction. By increasing the weight factor α, the effect of prior knowledge correction can be reduced, thereby reducing the overconstraint of prior knowledge on the inversion result.

[0025] (3)Finally, an inversion method based on three-dimensional space is presented. This method provides a grid construction method, a ray tracing method, and a tomographic inversion method based on three-dimensional space, and can also be applied to two-dimensional space through degradation, with generalizability. Description of the Drawings

[0026] Figure 1 It is a flowchart of the three-dimensional tomographic imaging method of the Fresnel body with prior knowledge correction according to the present invention;

[0027] Figure 2 It is a schematic diagram of the test instrument according to the present invention;

[0028] Figure 3 It is a scanning waveform table of the flexural element array according to the present invention;

[0029] Figure 4 It is a schematic diagram of the spatial grid according to the present invention;

[0030] Figure 5 It is a schematic diagram of the emission point and reception point according to the present invention;

[0031] Figure 6 It is a schematic diagram of the minimum travel time according to the present invention;

[0032] Figure 7 It is a schematic diagram of the ray path search according to the present invention;

[0033] Figure 8 It is a schematic diagram of the ray propagation path according to the present invention;

[0034] Figure 9 It is a schematic diagram of the Fresnel body in a uniform slowness field according to the present invention;

[0035] Figure 10 It is a schematic diagram of the Fresnel body in a non-uniform slowness field according to the present invention;

[0036] Figure 11 It is a schematic diagram of the prior slowness field shape parameter correction according to the present invention;

[0037] Figure 12 It is a flowchart of the prior slowness field shape parameter correction according to the present invention;

[0038] In the figure: 1 - model box; 2 - flexural element sensor array; 201 - vibration isolation pad; 202 - flexural element bracket; 203 - excitation element; 204 - receiving element; 3 - soil body. Detailed Embodiment

[0039] The present invention will be further described in detail below with reference to the drawings and specific embodiments.

[0040] As Figure 1As shown in the figure, a Fresnel body three-dimensional tomography method with prior knowledge correction provided by the present invention mainly includes the following steps:

[0041] 1) Pretreatment of test data: According to the waveform table scanned by the bender element array 2, after data pretreatment such as filtering each received wave, judging the propagation time, and screening valid data, a test time propagation table is obtained.

[0042] As Figure 2 shown, the bender element array 2 includes parts such as a seismic isolation pad 201, a bender element support 202, an excitation element 203, and a receiving element 204. The bender element array 2 is installed in the model box 1. A pair of excitation elements 203 are arranged on the left bender element support 202, and a pair of receiving elements 204 are arranged on the right bender element support 202. The excitation element 203 emits a single vibration under the action of the excitation voltage and generates a shear wave in the soil mass 3. After propagating through the soil mass 3, it reaches the receiving element 204 on the opposite side, thereby obtaining a received wave signal; a seismic isolation pad 201 is arranged below the bender element support to prevent waves from propagating through the model box 1 body.

[0043] As Figure 3 shown, the bender element array 2 tests the received wave waveforms obtained by combining each excitation element 203 and receiving element 204, and sorts them to obtain a bender element array scan waveform table; after pretreatment, the bender element array scan waveform table can be converted into a test propagation time table, which contains the serial numbers of the excitation elements 203, the serial numbers of the receiving elements 204, and the ray propagation time obtained by preprocessing the waveforms.

[0044] 2) Construction of the grid model: As Figure 4 and Figure 5 shown, a calculation grid is created, and an assumed initial slowness field is written at the grid nodes. According to the positions of the excitation element m and the receiving element n in the test, the corresponding emission point m and receiving point n are set at the grid nodes, and the serial numbers, spatial coordinates, and node numbers of the emission point and the receiving point are recorded.

[0045] The density of the calculation grid is determined according to requirements. Too many grids will lead to too long calculation time, while too few grids will make the calculation results unreasonable. It is necessary to gradually adjust the number of grids according to the spatial size of the soil mass in the test and the number of emission points and excitation points; the calculation grid is a cubic grid; since the initial soil layer generally has the characteristic of being soft on the top and hard on the bottom, the initial slowness field selects a slowness field that linearly decreases along the depth Z. Its selection rationality will affect the calculation efficiency and will also affect the calculation results to a certain extent;

[0046] 3) Tomographic forward modeling: As Figure 6As shown in the figure, for the emission point m, the Dijkstra algorithm is used to calculate the minimum travel time from it to each grid node on the full grid model. Then, for each receiving point n, the travel time of the forward tomographic calculation of each emission-reception pair is read in sequence according to its node number, and a forward propagation schedule is obtained, which contains information such as the emission point serial number m, the receiving point serial number n, and the ray propagation time. Subsequently, starting from the receiving point n, the gradient descent method is used to back-calculate the ray path of each emission-reception pair based on the minimum travel time of the grid, and finally, the Fresnel volume range and the normalized weight factor of each emission-reception pair are calculated to obtain the forward propagation path table.

[0047] As Figure 7 and Figure 8 shown, the gradient descent method needs to search for the travel times of grid points in different directions in three layers of grids with R = 1, R = 2, and R = 3 around the starting point, calculate their travel time gradients, and find the direction with the largest descending gradient as the ray propagation path until the emission point is retrieved.

[0048] As Figure 9 and Figure 10 shown, when calculating the Fresnel volume information, first determine the Fresnel volume range according to the formula , and then calculate the weight factor ω of the influence degree of different parts of the Fresnel volume on the travel time according to the formula . The closer the ray is to the central high-frequency part, the greater the influence weight, and the closer the weight factor is to 1. Then, use the formula ij to calculate the normalized weight factor Q . For an emission-reception pair i, the sum of the normalized weight factors Q ij is 1, and M is the total number of emission-reception pairs. ij

[0049] 4) Inversion iteration: Compare the experimental propagation schedule and the forward propagation schedule. If the error between the two is greater than the limit value, enter the inversion iteration step. First, calculate the slowness field correction value for the area passed by the Fresnel volume using the SIRT algorithm, and then calculate the slowness field correction value for the area not passed by the Fresnel volume using the back-projection algorithm. Finally, obtain Δs j(1) for each grid node, calculate the iteration factor λ according to the ray path and the slowness field correction value k , update the slowness field 1, and finally perform a scanning and reforming calculation of the slowness field according to the prior slowness field shape parameter to obtain Δs j(2) . According to Δs j(1) and Δs j(2) to obtain the final slowness field 2.

[0050] Specifically, judge whether the error between the experimental propagation schedule and the forward propagation schedule is less than the limit value according to the formula , where M is the total number of emission-reception pairs.

[0051] For the nodes through which the Fresnel body passes, the simultaneous iterative reconstruction technique SIRT is adopted. According to the formula calculate the slowness correction value of the Fresnel body formed by each emission-reception pair i for node j, and use for averaging. N is the number of Fresnel bodies passing through this node.

[0052] For the nodes through which no Fresnel body passes, use The formula performs inverse distance interpolation on the correction values in the correction area where the Fresnel body passes to obtain the slowness correction value of the node through which no Fresnel body passes.

[0053] According to the forward ray path slowness correction value and calculate the iteration factor to make the slowness inversion iteration converge faster; and use to update the slowness field 1;

[0054] As Figure 11 and 12 shown, in the updated slowness field 1, delimit the reconstruction area grid range k*k*k. Within the reconstruction area, use the prior slowness field shape parameters obtained from prior knowledge to scale the slowness field 1 to obtain Δs j(2) , then switch to the adjacent next reconstruction area, repeat the above operation until all grid areas are scanned. After such an operation, the shape of the slowness field obtained will be strongly similar to the prior slowness field shape parameters obtained from prior knowledge; finally, according to the formula Δs j =αΔs j(1) +(1-α)Δs j(2) , comprehensively consider the slowness field correction value Δs j(1) generated by the simultaneous iterative reconstruction technique SIRT and the slowness field correction value Δs j(2) generated by prior knowledge correction according to the weight factor α, and calculate the final slowness field correction value Δs j to update the slowness field 2.

[0055] 5) Use the slowness field 2 to perform tomographic forward modeling again, compare the test propagation schedule and the new forward propagation schedule. If the error between the two is less than the limit value, the calculation ends. If it is still greater than the limit value, repeat step 4) until the error is less than the limit value.

[0056] Embodiment

[0057] Build a test platform for the flexural element array 2, with 3 flexural element brackets 202 arranged on the left side. Longitudinally arranged on each flexural element bracket are 5 excitation elements 203. At the corresponding positions on the right side, 3 flexural element brackets 202 are arranged, and longitudinally arranged on each flexural element bracket are 5 receiving elements 204. The excitation elements 203 emit single vibrations in sequence under the action of the excitation voltage and generate shear waves in the soil mass 3. After propagating through the soil mass 3, they reach the receiving elements 204 on the opposite side, thereby obtaining the received wave signals. After one test, a total of 15×15 = 225 groups of wave signals are obtained, forming a flexural element array waveform table. After data preprocessing such as filtering, propagation time discrimination, and effective data screening for each received wave, a test time propagation table is obtained, which contains the numbers and propagation times of each excitation element and flexural element, totaling 225 groups.

[0058] Considering the spatial size of the soil mass and calculation efficiency in the test, build a spatial grid with a length of 35cm * 35cm * 70cm and 51 * 51 * 101 grid nodes; arrange 15 emission points and 15 receiving points at the corresponding grid nodes according to the positions of the excitation elements 203 and receiving elements 204; write a slowness field that linearly decreases with depth, with a surface slowness of 1 / 40 s / m and a bottom slowness of 1 / 80 s / m.

[0059] For emission point 1, use the Dijkstra algorithm to calculate the minimum travel time from it to each grid node on the full grid model. Then, for all receiving points, successively read the travel times of the tomographic forward calculations of each emission-reception pair according to their node numbers to obtain the forward propagation time table. Subsequently, starting from receiving point n, use the gradient descent method to back-calculate the ray paths of each emission-reception pair according to the minimum travel time of the grid, and finally calculate the Fresnel body range and normalized weight factor of each emission-reception pair composed of emission point 1 and receiving point n; then perform the above operations on emission point 2, emission point 3... emission point 15 in sequence until the forward information of 225 emission-reception pairs is obtained, forming the forward propagation time table and the forward propagation path table.

[0060] Use the formula Compare the test propagation time table and the forward propagation time table. If the error between the two is greater than the limit value, enter the inversion iteration step, which is as follows:

[0061] (1) Tomographic inversion updates the slowness field 1: Since there are 225 groups of Fresnel body rays passing through the space, most of the spatial grids have Fresnel body rays passing through once or multiple times, and there are also a small number of regions where no Fresnel body rays pass through. For the nodes with Fresnel bodies passing through, use the simultaneous iterative reconstruction technique SIRT. According to the formula Calculate the slowness correction value of the Fresnel body of each emission-reception pair i to node j, and use for averaging, where N is the number of Fresnel bodies passing through this node. For the nodes without Fresnel bodies passing through, use The formula performs inverse distance interpolation on the correction values of the Fresnel body passing through the correction area to obtain the slowness correction value of the node without the Fresnel body passing through. According to the forward ray path slowness correction value and the calculation iteration factor make the slowness inversion iteration converge faster; and use to update the slowness field 1;

[0062] (2) After the shape parameter of the prior slowness field is corrected, the slowness field 2 is obtained: In the updated slowness field 1, a reconstruction area grid range of k*k*k is delimited. Within the reconstruction area, the shape parameter of the prior slowness field obtained by prior knowledge is used to perform proportional adjustment on the slowness field 1 to obtain Δs j(2) , and then switch to the next adjacent reconstruction area, repeat the above operation until all grid areas are scanned. After such operation, the shape of the slowness field obtained will be strongly similar to the shape parameter of the prior slowness field obtained by prior knowledge; finally, according to the formula Δs j =αΔs j(1) +(1 - α)Δs j(2) , comprehensively consider the slowness field correction value Δs j(1) generated by the simultaneous iterative reconstruction technique SIRT and the slowness field correction value Δs j(2) generated by prior knowledge correction according to the weight factor α, and calculate the final slowness field correction value Δs j to update the slowness field 2.

[0063] After the iteration is completed, perform tomographic forward modeling again with the updated slowness field 2, compare the experimental propagation schedule and the new forward propagation schedule. If the error between the two is less than the limit value, the calculation ends. If it is still greater than the limit value, repeat the inversion iteration step until the error is less than the limit value.

Claims

1. A three-dimensional tomographic imaging method for Fresnel body with prior knowledge correction, characterized in that The specific steps are as follows: 1) Pretreatment of test data: According to the waveform table scanned by the bender element array (2), after data pretreatment, the test time propagation table is obtained; 2) Establishment of grid model: Create a computational grid model, write the assumed initial slowness field at the grid nodes, set the corresponding emission point m and reception point n on the grid nodes according to the positions of the excitation element m and reception element n in the test, and record the serial numbers, spatial coordinates and node numbers of the emission point and reception point; 3) Forward tomography: For the emission point m, calculate the minimum travel time from it to each grid node in the entire grid model. Then, for each reception point n, successively read the travel times calculated by forward tomography of each emission-reception pair according to its node number to obtain the forward propagation time table. Subsequently, starting from the reception point n, back-calculate the ray paths of each emission-reception pair according to the minimum travel time of the grid, and finally calculate the Fresnel volume range and normalized weight factor of each emission-reception pair to obtain the forward propagation path table; 4) Inversion iteration: Compare the test propagation time table and the forward propagation time table. If the error between the two is greater than the limit value, enter the inversion iteration step. First, calculate the slowness field correction value for the area passed by the Fresnel volume using the SIRT algorithm, and then calculate the slowness field correction value for the area not passed by the Fresnel volume using the back-projection algorithm. Finally, obtain the corresponding slowness correction value for each grid node. Calculate the iteration factor according to the ray path and the slowness field correction value to update the slowness field 1. Finally, perform a scanning and reordering calculation on the slowness field 1 according to the prior slowness field shape parameter to obtain the slowness field correction value generated by prior knowledge modification, and comprehensively consider the slowness field correction values before and after prior knowledge modification to obtain the final slowness field 2; 5) Re-perform forward tomography using the slowness field 2, compare the test propagation time table and the new forward propagation time table. If the error between the two is less than the limit value, the calculation ends. If it is still greater than the limit value, repeat step 4) until the error is less than the limit value.

2. A Fresnel body three-dimensional tomography method with prior knowledge correction according to claim 1, characterized in that In step 1), the bender element array (2) includes a seismic isolation pad (201), a bender element support (202), an excitation element (203), and a reception element (204). The bender element array (2) is installed in the model box (1). Among them, several pairs of excitation elements (203) are arranged on the left bender element support (202), and the corresponding several pairs of reception elements (204) are arranged on the right bender element support (202). The excitation element (203) emits a single vibration under the action of the excitation voltage and generates a shear wave in the soil body (3). After propagating through the soil body (3), it reaches the reception element (204) on the opposite side, thereby obtaining the received wave signal. A seismic isolation pad (201) is arranged below the bender element support to prevent the wave from propagating through the box body of the model box (1).

3. A Fresnel body three-dimensional tomography method with prior knowledge correction according to claim 1, characterized in that, In step 1), the bender element array (2) tests the received wave waveforms obtained by combining each excitation element (203) and reception element (204), and thus organizes and obtains the bender element array scanning waveform table. After preprocessing each received wave by filtering, propagation time discrimination, and effective data screening, the bender element array scanning waveform table is converted into a test propagation time table, which contains the serial number of the excitation element (203), the serial number of the reception element (204), and the ray propagation time obtained by waveform preprocessing.

4. A Fresnel body three-dimensional tomography method with prior knowledge correction according to claim 1, characterized in that In step 2), the computational grid is a cubic grid; the initial slowness field is selected as a slowness field that linearly decreases along the depth Z.

5. A Fresnel body three-dimensional tomography method with prior knowledge correction according to claim 1, characterized in that In step 3), when back-calculating the ray paths of each transmitter-receiver pair based on the minimum travel time of the grid, it is necessary to search for the travel times of grid points in different directions in the three layers of grids with R = 1, R = 2, and R = 3 around the starting point, calculate their travel time gradients, and find the direction with the largest descending gradient as the ray propagation path until the transmitter point is retrieved.

6. A Fresnel body three-dimensional tomography method with prior knowledge correction according to claim 1, characterized in that In step 3), according to the formula the Fresnel volume range is determined, where S represents the emission point, R represents the reception point, P represents any point within the spatial grid, the subscript i means the i-th emission-reception pair, j represents the j-th node, and T SP is the travel time from the emission point to any point; T PR is the travel time from any point to the reception point; T SR is the travel time from the emission point to the reception point, f is the ray frequency (unit s -1 ), and then according to the formula the weight factor ω of the influence degree of different parts of the Fresnel volume on the travel time is calculated ij . The influence weight of the ray is greater for the part closer to the central high-frequency part, and the weight factor is closer to 1; Then calculate the normalized weight factor Q ij , For a transmit-receive pair i, the sum of the normalized weight factor Q ij is 1.

7. A Fresnel body three-dimensional tomography method with prior knowledge correction according to claim 1, characterized in that In step 4), according to the formula it is determined whether the error between the test propagation schedule and the forward propagation schedule is less than the limit value, where M is the total number of transmitter-receiver pairs, and ΔT is the average relative error between the test propagation schedule and the forward propagation schedule; ΔT lim is the average relative error limit value between the test propagation schedule and the forward propagation schedule; T i is the forward propagation time of the i-th transmit-receive pair; is the trial propagation time of the i-th transmit-receive pair.

8. A Fresnel body three-dimensional tomography method with prior knowledge correction according to claim 1, characterized in that In step 4), for the nodes through which the Fresnel volume passes, the simultaneous iterative reconstruction technique SIRT is adopted. According to the formula calculate the slowness correction value of the Fresnel volume formed by each emission-reception pair i for node j, and use for averaging. N is the number of Fresnel volumes passing through this node, and s j represents the initial slowness of the j-th node. For the nodes through which no Fresnel volume passes, use The formula performs inverse distance interpolation on the correction values in the correction area where the Fresnel volume passes to obtain the slowness correction value of the nodes through which no Fresnel volume passes. L represents the total number of nodes in the Fresnel volume, and d j represents the distance between this point and the j-th node in the Fresnel volume. represents the slowness correction amount of the j-th node in the Fresnel volume.

9. A Fresnel body three-dimensional tomography method with prior knowledge correction according to claim 1, characterized in that In step 4), according to the forward ray path slowness correction value and calculate the iteration factor to make the slowness inversion iteration converge faster; and use to update the slowness field 1, where k represents the k-th iteration; λ k represents the iteration factor of the k-th iteration; represents the forward travel time of the i-th transmitter-receiver pair in the k-th iteration; represents the ray path length at the j-th node of the i-th transmitter-receiver pair in the k-th iteration; N i represents the number of nodes passed by the ray of the i-th transmitter-receiver pair.

10. A Fresnel body three-dimensional tomography method with prior knowledge correction according to claim 1, characterized in that, In step 4), a reconstruction area grid range of k*k*k is delimited in the updated slowness field 1. Within the reconstruction area, the slowness field 1 is proportionally adjusted using the prior slowness field shape parameters obtained from prior knowledge. Subsequently, switch to the next adjacent reconstruction area and repeat the proportional adjustment operation until all grid areas are scanned. After such operations, the shape of the slowness field obtained will have a strong similarity to the prior slowness field shape parameters obtained from prior knowledge; Comprehensively consider the slowness field correction value generated by the simultaneous iterative reconstruction technique SIRT and the slowness field correction value generated by prior knowledge correction according to the weight factor α, and calculate the final slowness field correction value to update the slowness field 2.

Images