High-precision numerical simulation method in frequency domain for three-dimensional acoustic wave equation of high-speed railway seismic source

Abstract

The present invention discloses a high-precision numerical simulation method for three-dimensional acoustic wave equation in the frequency domain of high-speed railway seismic sources, including clarifying the definite solution problem for numerical solution, including the governing equation and boundary conditions. The governing equation adopted is the three-dimensional frequency-domain acoustic wave equation, the boundary condition is set that the boundary field value is zero, and the seismic source term of the equation is the high-speed railway seismic source; the 33-point finite difference scheme is used to discretize the spatial derivative term and the mass acceleration term of the three-dimensional acoustic wave equation. The present invention realizes the high-precision numerical simulation of the high-speed railway wave field in the frequency domain, fills the gap in the forward modeling method of the high-speed railway wave field in the frequency domain, provides a forward modeling module for the full waveform inversion of the high-speed railway wave field in the frequency domain, which will provide new ideas for high-speed railway data imaging and strongly promote the practical application of near-surface imaging along and around the high-speed railway.

Description

Technical Field

[0001] The present invention relates to the technical field of seismic wave numerical simulation, and particularly to a high-precision numerical simulation method for the three-dimensional acoustic wave equation in the frequency domain of a high-speed rail seismic source. Background Art

[0002] Traditional seismic exploration usually uses seismic sources such as explosives, vibroseis, or air guns to excite seismic waves for imaging the underground space. However, such seismic sources are costly and environmentally destructive. Nowadays, with the sudden increase in transportation demand, high-speed rail has become a major choice for people to travel quickly. A moving high-speed rail, due to its nearly constant running speed, determined length, and load, can become a repeatable artificial seismic source, and the crisscross high-speed rail system broadens the distribution range of this artificial seismic source, thus enabling it to have great potential for detecting the internal structure and physical property parameters of the earth. Compared with traditional seismic sources, high-speed rail seismic sources have the advantages of strong repeatability, zero cost, non-destruction of the environment, and rich signal spectrum information. Therefore, it is of great value to use the vibration signals generated by moving high-speed rails to image the underground structure and estimate relevant underground characteristics.

[0003] Accurate high-speed rail wavefield simulation technology is the basis for using the high-speed rail wavefield to carry out near-surface structure imaging and parameter inversion. The existing high-speed rail wavefield simulation technologies mainly focus on two aspects: the first aspect, analytical and semi-analytical methods based on full space, half space, and horizontally layered models, and the second aspect, time-domain numerical simulation methods such as finite difference methods and spectral element methods. The numerical simulation schemes in the frequency domain are mainly for fixed seismic sources to synthesize wavefield signals excited by traditional seismic sources. In summary, the frequency-domain numerical simulation technology for high-speed rail wavefields is currently in a blank state.

[0004] Simulating the high-speed rail wavefield using analytical or semi-analytical methods can effectively analyze the multi-domain characteristics of the high-speed rail wavefield, but such methods can only handle several simple high-symmetry models such as homogeneous or horizontally layered models. For complex underground space models, it is difficult to find corresponding analytical solutions.

[0005] To address the above problems, the present invention provides a high-precision numerical simulation method for the three-dimensional acoustic wave equation in the frequency domain of a high-speed rail seismic source for achieving high-precision frequency-domain inversion. Summary of the Invention

[0006] The present invention provides a high-precision numerical simulation method for the three-dimensional acoustic wave equation in the frequency domain of a high-speed rail seismic source for achieving high-precision frequency-domain inversion.

[0007] The objective of the present invention is to provide a high-precision numerical simulation method for the three-dimensional acoustic wave equation in the frequency domain of a high-speed rail seismic source, and the solution includes the following steps:

[0008] Step 1: Define the definite solution problem for numerical solution, including the governing equation and boundary conditions. The governing equation used is the three-dimensional frequency-domain acoustic wave equation, and the boundary condition is set such that the boundary field value is zero. The source term of the equation is the high-speed rail source;

[0009] Step 2: Discretize the spatial derivative term and the mass acceleration term of the three-dimensional acoustic wave equation using a 33-point finite difference scheme;

[0010] Step 3: Without considering the source term, obtain the dispersion relation of the 33-point difference scheme in a homogeneous infinite medium, and determine the weighting coefficient through optimization;

[0011] Step 4: To avoid boundary reflections, apply the PML absorbing boundary condition around the target model;

[0012] Step 5: Simplify the modeling of the high-speed rail pier model, construct the high-speed rail source function, and design the loading method of the high-speed rail source;

[0013] Step 6: Obtain the high-speed rail frequency-domain single-frequency wave field.

[0014] Furthermore, the specific calculation steps of the governing equation and boundary conditions in Step 1 are as follows:

[0015] Step 1: The three-dimensional frequency-domain acoustic wave equation and boundary conditions are expressed as:

[0016]

[0017] where x = (x, y, z) represents the spatial position, p is the pressure wave field, ω is the angular frequency, v is the propagation speed, is the Laplace operator, and S is the source term;

[0018] Step 2: For a traditional fixed source, S can usually be expressed as:

[0019] S(x, y, z, ω) = S(ω)δ(x - x s )δ(y - y s )δ(z - z s ),

[0020] where S(ω) is the source function in the frequency domain, and (x s , y s , z s ) is the source spatial position;

[0021] Step 3: When the high-speed rail moves uniformly along the x direction, the source function of the m-th pier can be expressed as:

[0022]

[0023] where (x m , ys , z s ) represents the spatial position of the m-th pier, v represents the moving speed of the high-speed rail, represents the imaginary unit, and Q is the modulation factor expressed as:

[0024]

[0025] It can be seen from the expression of Q that the high-speed rail source function is modulated by the carriage parameters a and b and the moving speed v. The factor P is expressed as:

[0026]

[0027] where N is the number of carriages.

[0028] Further, the collective calculation steps of the 33-point finite difference scheme in step two are as follows:

[0029] Step one: For the discretization of the three-dimensional frequency-domain acoustic wave equation, we propose a new 33-point finite difference scheme. For the second-order spatial derivative term The 33-point discrete scheme is expressed as:

[0030]

[0031] where represents the average derivative 27-point operator, represents the classical 7-point operator.

[0032] Step two: The average derivative 27-point operator is expressed as:

[0033]

[0034] where

[0035]

[0036] where p m,l,n ≈ p(mΔx, lΔy, nΔz), Δx, Δy, and Δz respectively represent the grid intervals in the x, y, and z directions, and a IJ (I, J = 1, 2, 3) are the weighting coefficients, is expressed as:

[0037]

[0038] Step two: For the mass acceleration term, the 33-point weighting operator is expressed as:

[0039]

[0040] where

[0041]

[0042] Among them, b i (i = 0, 1, 2, 3, 4) represents the weighting coefficient. Therefore, the 33-point finite difference discretization equation can be obtained as follows:

[0043]

[0044] Step 3: When w = 1 and b 4 = 0, the 33-point method degenerates into the 27-point average derivative method.

[0045] Furthermore, the specific calculation steps of the dispersion relation in Step 3 are as follows:

[0046] Step 1: Evaluate the calculation accuracy of the 33-point difference scheme and determine its weighting coefficient, and conduct a classical dispersion analysis. Consider the classical plane wave solution in three-dimensional cases, and the specific expression is:

[0047]

[0048] Among them, p 0 is a constant, k is the wave number, θ is the propagation angle, and φ is the azimuth angle;

[0049] Step 2: Insert the discrete format in Step 1 into the 33-point finite difference discretization equation, and the normalized phase velocity can be obtained. The specific expression is:

[0050]

[0051] Among them, v ph is the phase velocity, represents the number of grid points per wavelength (PPW), Without loss of generality, Δx = max{Δx, Δy, Δz},

[0052]

[0053]

[0054] Step 3: In order to suppress the dispersion error, consider the following error function of the L 2 norm:

[0055]

[0056] Among them, θ, φ and are respectively set in the ranges of and [0, 0.36], and the constrained nonlinear optimization program fmincon in MATLAB is used for optimization.

[0057] Furthermore, the specific calculation steps for the boundary conditions in Step 4 are as follows:

[0058] Step 1: To suppress the artificial reflections generated by the boundaries introduced to limit the model size, apply the PML absorbing boundary condition around the target model. The three-dimensional acoustic wave equation in the frequency domain considering the PML absorbing boundary condition is:

[0059]

[0060] where,

[0061]

[0062] where, L x , L y and L z are the thicknesses of the PML in the x, y, and z directions respectively. The starting points in the x, y, and z directions are the outer edges of the internal model. c x = c y = c z = 180.

[0063] Furthermore, the specific calculation steps for Step 5 are as follows:

[0064] Step 1: Considering all the grid points of the model, the linear equation system is expressed as:

[0065] Ap = -s,

[0066] where, A is a sparse impedance matrix of size (N x N y N z )×(N x N y N z ), p is a wave field vector of size N x N y N z ×1, s is a source vector of size N x N y N z ×1, and N x , N y and N z are the number of grid points in the x, y, and z directions respectively;

[0067] Step 2: In the finite difference grid template, the high-speed rail source can be expressed as:

[0068]

[0069] where, P is a zero-padding sorting operator to ensure that s is N x N y N zA vector of size ×1; T represents matrix transpose, and M is the number of bridge piers.

[0070] The present invention has the following advantages: The present invention realizes the high-precision numerical simulation of the high-speed railway wave field in the frequency domain, fills the gap in the forward modeling method of the high-speed railway wave field in the frequency domain, provides a forward modeling module for the full waveform inversion of the high-speed railway wave field in the frequency domain, which will provide new ideas for high-speed railway data imaging and strongly promote the practical application of near-surface imaging along and around the high-speed railway. BRIEF DESCRIPTION OF THE DRAWINGS

[0071] Figure 1 is a flowchart of the present invention;

[0072] Figure 2 is a schematic diagram of a high-speed railway running on an elevated structure of the present invention;

[0073] Figure 3 is a high-speed railway source function diagram based on the first derivative of the E-BEB model of the present invention;

[0074] Figure 4 is a 33-point finite difference template diagram of the present invention;

[0075] Figure 5 is a dispersion curve diagram of the 27-point format and 33-point format of the present invention;

[0076] Figure 6 is a loading diagram of the high-speed railway source in the finite difference grid of the present invention;

[0077] Figure 7 is a source time function diagram of different bridge piers of the present invention;

[0078] Figure 8 is the calculation result and analytical solution of the 27-point and 33-point difference formats of the present invention;

[0079] Figure 9 is a frequency domain waveform profile calculated by the 27-point and 33-point difference formats of the present invention;

[0080] Figure 10 is a 15 Hz high-speed railway single-frequency wave field diagram calculated by using the 33-point difference format under the uniform model of the present invention;

[0081] Figure 11 is a three-dimensional SEG / EAGE overthrust velocity model diagram of the present invention;

[0082] Figure 12 is a 7 Hz high-speed railway single-frequency wave field diagram calculated by using the 33-point difference format under the three-dimensional SEG / EAGE overthrust model of the present invention. DETAILED DESCRIPTION OF THE INVENTION

[0083] The present invention provides a high-precision numerical simulation method for the three-dimensional acoustic wave equation in the frequency domain of a high-speed rail seismic source. The method includes the following steps:

[0084] Step 1: Define the definite solution problem for numerical solution, including the governing equation and boundary conditions. The governing equation used is the three-dimensional frequency-domain acoustic wave equation, the boundary condition is set that the boundary field value is zero, and the seismic source term of the equation is the high-speed rail seismic source;

[0085] Step 2: Discretize the spatial derivative term and the mass acceleration term of the three-dimensional acoustic wave equation using a 33-point finite difference scheme;

[0086] Step 3: Without considering the seismic source term, in a homogeneous infinite medium, obtain the dispersion relationship of the 33-point difference scheme and determine the weighting coefficient through optimization;

[0087] Step 4: To avoid boundary reflections, apply a PML absorbing boundary condition around the target model;

[0088] Step 5: Simplify the modeling of the high-speed rail pier model, construct the high-speed rail seismic source function, and design the loading method of the high-speed rail seismic source;

[0089] Step 6: Obtain the high-speed rail frequency-domain single-frequency wave field.

[0090] In this embodiment, the specific calculation steps of the governing equation and boundary conditions in Step 1 are as follows:

[0091] Step 1: The three-dimensional frequency-domain acoustic wave equation and boundary conditions are expressed as:

[0092]

[0093] where x = (x, y, z) represents the spatial position, p is the pressure wave field, ω is the angular frequency, v is the propagation speed, is the Laplace operator, and S is the seismic source term;

[0094] Step 2: For a traditional fixed seismic source, S can usually be expressed as:

[0095] S(x, y, z, ω) = S(ω)δ(x - x s )δ(y - y s )δ(z - z s ),

[0096] where S(ω) is the seismic source function in the frequency domain, and (x s , y s , z s ) is the seismic source spatial position;

[0097] Step 3: When the high-speed rail moves uniformly along the x direction, the seismic source function of the m-th pier can be expressed as:

[0098]

[0099] Among them, (x m , y s , z s ) represents the spatial position of the m-th pier, v represents the moving speed of the high-speed rail, represents the imaginary unit, and Q is the modulation factor expressed as:

[0100]

[0101] It can be seen from the expression of Q that the high-speed rail source function is modulated by the carriage parameters a and b and the moving speed v, and the factor P is expressed as:

[0102]

[0103] Among them, N is the number of carriages.

[0104] In this embodiment, the collective calculation steps of the 33-point finite difference scheme in step two are as follows:

[0105] Step 1: For the discretization of the three-dimensional frequency-domain acoustic wave equation, we propose a new 33-point finite difference scheme. For the second-order spatial derivative term The 33-point discretization scheme is expressed as:

[0106]

[0107] Among them, represents the 27-point operator of the average derivative, represents the classical 7-point operator.

[0108] Step 2: The 27-point operator of the average derivative is expressed as:

[0109]

[0110] Among them,

[0111]

[0112] Among them, p m,l,n ≈ p(mΔx, lΔy, nΔz), Δx, Δy, and Δz respectively represent the grid intervals in the x, y, and z directions, and a IJ (I, J = 1, 2, 3) are the weighting coefficients, is expressed as:

[0113]

[0114] Step 2: For the mass acceleration term, the 33-point weighting operator is expressed as:

[0115]

[0116] Among them,

[0117]

[0118] Among them, b i (i = 0, 1, 2, 3, 4) represents the weighting coefficient. Therefore, the 33-point finite difference discrete equation can be obtained as:

[0119]

[0120] Step 3: When w = 1 and b 4 = 0, the 33-point method degenerates into the 27-point average derivative method.

[0121] In this embodiment, the specific calculation steps of the dispersion relationship in Step 3 are as follows:

[0122] Step 1: Evaluate the calculation accuracy of the 33-point difference scheme and determine its weighting coefficient, and conduct a classical dispersion analysis. Consider the classical plane wave solution in three-dimensional cases, and the specific expression is:

[0123]

[0124] Among them, p 0 is a constant, k is the wave number, θ is the propagation angle, and φ is the azimuth angle;

[0125] Step 2: Insert the discrete format in Step 1 into the 33-point finite difference discrete equation, and the normalized phase velocity can be obtained. The specific expression is:

[0126]

[0127] Among them, v ph is the phase velocity, represents the number of grid points per wavelength (PPW), Without loss of generality, Δx = max{Δx, Δy, Δz},

[0128]

[0129] Step 3: To suppress the dispersion error, consider the following error function of the L 2 norm:

[0130]

[0131] Among them, θ, φ and are respectively set in the range of And [0, 0.36], and optimize it using the constrained nonlinear optimization program fmincon in MATLAB.

[0132] In this embodiment, the specific calculation steps of the boundary conditions in Step 4 are as follows:

[0133] Step 1: To suppress the artificial reflections generated by the boundaries introduced to limit the model size, apply the PML absorbing boundary condition to the periphery of the target model. The three-dimensional acoustic wave equation in the frequency domain considering the PML absorbing boundary condition is:

[0134]

[0135] Where,

[0136]

[0137] Where, L x , L y And L z Are the thicknesses of the PML in the x, y, and z directions respectively. The starting points in the x, y, and z directions are the outer edges of the internal model. c x = c y = c z = 180.

[0138] In this embodiment, the specific calculation steps of Step 5 are as follows:

[0139] Step 1: Considering all the grid points of the model, the linear equations are expressed as:

[0140] Ap = -s,

[0141] Where, A is a sparse impedance matrix of size (N x N y N z ) × (N x N y N z ), p is a wave field vector of size N x N y N z ×1, s is a source vector of size N x N y N z ×1, N x , N y And N z Are the number of grid points in the x, y, and z directions respectively;

[0142] Step 2: In the finite difference grid template, the high-speed rail source can be expressed as:

[0143]

[0144] where P is a zero-padding sorting operator to ensure that s is an N x N y N z vector of size N×1; T represents matrix transpose, and M is the number of bridge piers.

[0145] In this embodiment, three numerical values are analyzed theoretically: (1) calculation of the fixed-source wave field of a homogeneous model, (2) calculation of the high-speed railway wave field of a homogeneous model, and (3) calculation of the high-speed railway wave field of the SEG / EAGE overthrust model. The computing platforms for the first two examples are: a server with 128 GB of RAM and an Intel Xeon 3.3 GHz processor, calling 8 CPU cores; the computing platform for the third example is: a server with 512 GB of RAM and an Intel Xeon 3.2 GHz processor, calling 16 CPU cores. First, consider a homogeneous model with v = 3000 m / s. The number of grid points is 121×121×21, and the grid size is Δx = Δy = Δz = 72 m. A Ricker wavelet with a dominant frequency of 15 Hz is placed in the middle of the model. The analytical solution of the three-dimensional homogeneous acoustic wave model is where, The wave field is calculated using the analytical method, 27-point, and 33-point difference schemes respectively. For the 27-point difference scheme, the memory occupancy is 24 GB and the calculation time is 1426 s; for the 33-point difference scheme, the occupancy is 26 GB and the calculation time is 1449 s. As shown, as Figure 8 shown, the difference between the calculation result of the 33-point difference scheme and the analytical solution is significantly smaller. The frequency-domain waveform profiles are shown. In contrast, the waveform of the 33-point difference scheme fits better with the analytical solution. Therefore, from the above calculation results, the following conclusion can be drawn: at a similar or slightly higher computational cost, the computational accuracy of the 33-point difference scheme is significantly higher than that of the 27-point difference scheme. Further, it can be inferred that at the same computational accuracy, the computational cost of the 33-point difference scheme should be lower than that of the 27-point difference scheme.

[0146] In this embodiment, a homogeneous model with v = 3000 m / s is still considered. The number of grid points is 121×121×21, and the grid size is Δx = Δy = Δz = 50 m. It is assumed that the high-speed railway source moves from the left end to the right end of the model along the x direction (y = 3000 m, z = 0 m). The single-frequency wave field at 15 Hz is calculated. Within the model scale, the starting and ending points of the excitation by the bridge piers are (0, 3000, 0) m and (5984, 3000, 0) m respectively, with a total of 188 bridge piers. As Figure 10 shows the spatial distribution of the single-frequency wave field of the high-speed railway source.

[0147] In this embodiment, as Figure 11The figure shows a three-dimensional SEG / EAGE overthrust model. The calculation frequency is 7 Hz. The number of grid points is 197×197×55, and the grid size is Δx = Δy = Δz = 80 m. The starting and ending points of the pier excitation are (0, 7840, 0) m and (7840, 7840, 0) m respectively, with a total of 491 piers. Figure 12 It shows the single-frequency wave field of high-speed rail at 7 Hz.

[0148] Although the specific implementation manners of the present invention have been described in detail with reference to the accompanying drawings, it should not be construed as a limitation on the protection scope of this patent. Within the scope described in the claims, various modifications and deformations that can be made by those skilled in the art without creative efforts still fall within the protection scope of this patent.

Claims

1. A high-precision numerical simulation method of three-dimensional acoustic wave equation in frequency domain for high-speed railway earthquake source, characterized by: The method comprises the following steps: Step 1: Clarify the well-defined problem to be solved numerically, including the control equation and boundary conditions. The control equation used is the three-dimensional frequency domain acoustic wave equation. The boundary condition is set to zero boundary field value, and the source term of the equation is the high-speed rail source. Step 2: discretize the spatial derivative term and mass acceleration term of the three-dimensional acoustic wave equation using a 33-point finite difference format; The specific calculation steps of the 33-point finite difference format in step 2 are as follows: Step 21: For the discretization of the three-dimensional frequency domain acoustic wave equation, for the second-order spatial derivative term The 33-point discrete format is represented as: in, represents the average derivative 27-point operator, represents the classic 7-point operator, w is and The weighting coefficient of Step 22, the average derivative 27 point operator is expressed as: in, Among them, p m,l,n ≈p(mΔx,lΔy,nΔz), where Δx, Δy and Δz represent the grid spacing in the x, y and z directions respectively, and a IJ is the weighting coefficient, I,J=1,2,3, It is expressed as: Step 23: For the mass acceleration term, the 33-point weighted operator is expressed as: in, Among them, b i It means that the weight coefficient i=0,1,2,3,4, so the 33-point finite difference discrete equation is obtained as follows: Step 24, when w=1 and b4=0, the 33-point method degenerates into the 27-point average derivative method; Step 3, without considering the source term, obtain the dispersion relation of the 33-point difference format in a uniform infinite medium, and determine the weighting coefficient by optimization; Step 4: To avoid boundary reflection, apply PML absorbing boundary conditions around the target model; Step 5: Simplify the high-speed railway pier model, construct the high-speed railway source function, and design the loading method of the high-speed railway source; Step 6: Obtain the high-speed rail frequency domain single-frequency wave field.

2. The high-precision numerical simulation method of the three-dimensional acoustic wave equation frequency domain of the high-speed railway earthquake source according to claim 1 is characterized by: The specific calculation steps of the control equations and boundary conditions in step 1 are as follows: Step 11: The three-dimensional frequency domain acoustic wave equation and boundary conditions are expressed as: Where X = (x, y, z) represents the spatial position, p is the pressure wave field, ω is the angular frequency, and v is the propagation velocity. is the Laplace operator, S is the source term; Step 12: For a traditional fixed source, S is expressed as: S(x,y,z,ω)=S(ω)δ(x-x s )δ(y-y s )δ(z-z s ), Among them, S(ω) is the source function in the frequency domain, (x s ,y s ,z s ) is the spatial location of the earthquake source; Step 13: When the high-speed railway moves at a constant speed along the x direction, the source function of the mth bridge pier is expressed as: Among them, (x m ,y s ,z s ) represents the spatial position of the mth bridge pier, v represents the moving speed of the high-speed rail, represents the imaginary unit, Q is the modulation factor expressed as: From the expression of Q, it can be seen that the source function of the high-speed rail is modulated by the carriage parameters a and b and the moving speed v, and the factor P is expressed as: Among them, N is the number of carriages and L is the length of a single carriage.

3. The high-precision numerical simulation method of the three-dimensional acoustic wave equation frequency domain of the high-speed railway earthquake source according to claim 1 is characterized by: The specific calculation steps of the dispersion relation in step 3 are as follows: Step 31, evaluate the calculation accuracy of the 33-point difference format and determine its weighting coefficient, perform classical dispersion analysis, consider the classical plane wave solution in three-dimensional case, the specific expression is: Where p0 is a constant, k is the wave number, θ is the propagation angle, and φ is the azimuth angle; Step 32: insert the discrete format of step 21 into the 33-point finite difference discrete equation to obtain the normalized phase velocity. The specific expression is: Among them, v ph is the phase velocity, represents the number of grid points per wavelength (PPW), Without loss of generality, Δx=max{Δx,Δy,Δz}, Step 33: In order to suppress the dispersion error, consider the following L2 norm error function: in, θ, φ and The ranges are set to and [0,0.36], and were optimized using the constrained nonlinear optimization program fmincon in MATLAB.

4. The high-precision numerical simulation method of the three-dimensional acoustic wave equation frequency domain of the high-speed railway earthquake source according to claim 1 is characterized by: The specific calculation steps of the boundary conditions in step 4 are as follows: In order to suppress the artificial reflection generated by the boundary introduced to limit the model size, a PML absorption boundary condition is applied to the periphery of the target model. The frequency domain three-dimensional acoustic wave equation considering the PML absorption boundary condition is: in, Among them, L x , L y and L z are the thickness of the PML in the x, y and z directions respectively. The starting point in the x, y and z directions is the outer edge of the inner model. x =c y =c z =180.

5. The high-precision numerical simulation method of the three-dimensional acoustic wave equation frequency domain of the high-speed railway earthquake source according to claim 1 is characterized by: The specific calculation steps of step 5 are as follows: Step 51, considering all grid points of the model, the linear equation system is expressed as: Ap=-s, Where A is a number of size (N x N y N z )×(N x N y N z ), p is a sparse impedance matrix of size N x N y N z ×1 wave field vector, s is a wave field vector of size N x N y N z ×1 source vector, N x , N y and N z are the number of grid points in the x, y and z directions respectively; Step 52: In the finite difference grid template, the high-speed rail source is represented as: Among them, P is the zero-filling sorting operator, ensuring that s is N x N y N z ×1 vector; T represents matrix transpose, M is the number of piers, and d is the interval between adjacent piers.

Images