Seismic wave field simulation method based on convolution kernel decomposition-Winograd fast convolution algorithm

Abstract

A seismic wave field simulation method based on convolution kernel decomposition-Winograd fast convolution algorithm is provided to solve the problem of low efficiency of convolution operation in large-scale seismic wave field simulation. Efficient seismic wave field simulation method: 1. Establish the intrinsic connection between spatial difference and convolution; 2. Establish the Winograd acceleration form of F(3,2) convolution; 3. Construct a seismic wave field. For any length of velocity field or stress field p and convolution kernel q, first decompose p into several sub-blocks so that the length of each sub-block can form a convolution process of size F(3,j) with convolution kernel q; 4. Divide the convolution kernel into blocks and apply the F(3,2) convolution acceleration form; 5. Merge all sub-block results of velocity or stress field p to obtain the differential field of velocity or stress field. The present invention derives the acceleration form of F(3,2) convolution process, which can decompose large-size convolution kernels into F(3,2) basic units for acceleration.

Description

Technical Field

[0001] The invention belongs to the field of seismology, and in particular relates to an efficient simulation method for seismic wave fields based on convolution kernel decomposition-Winograd fast convolution algorithm. Background Art

[0002] Seismic wavefield simulation is the process of simulating the propagation of seismic waves in the earth's internal medium on a computer. It has important applications and values ​​in many fields such as geological exploration, inversion of the earth's internal structure, and earthquake disaster assessment. At present, the most commonly used numerical method for seismic wavefield simulation is the finite difference method, which has the advantages of simple principle, easy implementation, and fast operation speed. Its basic principle is to approximate the differential terms in the wave equation through time-space differences, and for the spatial difference with the largest computational complexity, its calculation is equivalent to the convolution process. With the continuous advancement of urbanization in my country, when conducting seismic design or evaluation of megacities and urban agglomerations, we often face the need for large-scale or even super-large-scale seismic wavefield simulation. Therefore, there is an urgent need to improve the computational efficiency of convolution and develop a fast convolution algorithm.

[0003] The most commonly used fast convolution algorithm is Winograd fast convolution, the basic principle of which is to reduce the number of multiplications required for convolution operations by constructing an accelerated form, significantly reducing the amount of calculation. However, the Winograd algorithm currently only supports convolution kernels of odd lengths such as 3, that is, the Winograd basic unit of F(2,3) is used for calculation. This method does not match the convolution kernel of even order (4,8...) and cannot be directly applied. Therefore, the present invention proposes an efficient simulation method of seismic wave fields based on convolution kernel decomposition-Winograd fast convolution algorithm. Summary of the invention

[0004] The purpose of the present invention is to solve the problem of low efficiency of convolution operation in large-scale seismic wave field simulation and to provide a seismic wave field efficient simulation method based on convolution kernel decomposition-Winograd fast convolution algorithm.

[0005] The present invention provides a method for efficiently simulating seismic wave fields based on convolution kernel decomposition-Winograd fast convolution algorithm according to the following steps:

[0006] Step 1: Establish the intrinsic connection between spatial difference and convolution;

[0007] In the staggered grid finite difference scheme for the wave equation, the spatial difference of velocity or stress at any location is expressed as:

[0008]

[0009] where f n′ represents the velocity or stress differential value at n points in the staggered grid, c -i represents the convolution kernel, that is, the differential coefficient, f n-i represents the velocity or stress at position ni in the staggered grid, and M represents half of the difference order;

[0010] Step 2: Establish the Winograd acceleration form of F(3,2) convolution:

[0011] Let F(3,2) convolution refer to the convolution process with a convolution kernel length of 2 and a convolution output length of 3;

[0012] The convolution calculation process of the signal d(d0,d1,d2,d3) and the convolution kernel g(g0,g1) is d(d0,d1,d2,d3)*g(g0,g1)=r(r0,r1,r2), where r0=g0d0+d1g1, r1=g0d1+g1d2, r2=g0d2+g1d3;

[0013] Suppose vector m(m1,m2,m3,m4), if m1=(d0-d2)g0,m2=(d1+d2)(g0+g1) / 2,m3=(d2-d1)(g0-g1) / 2,m4=(d3-d1)g1, then r0=m1+m2+m3,r1=m2-m3,r3=m2+m3+m4, the Winograd acceleration form of F(3,2) convolution is established;

[0014] Step 3: Velocity or stress field segmentation:

[0015] Construct seismic wave fields. For any length of velocity field or stress field p(p0,p1,...,p i-1 ) and the convolution kernel q(q0,q1,...,q j-1 ), first decompose the velocity field or stress field p into several sub-blocks, so that the length of each sub-block can form a convolution process of size F(3,j) with the convolution kernel q; the length of each sub-block after p is divided into blocks is 3+j-1, the step size is 3, and there is overlap between each sub-block, and the overlap length is j-1;

[0016] Step 4: Divide the convolution kernel into blocks and apply the Winograd acceleration form of F(3,2) convolution:

[0017] For the convolution kernel q(q0,q1,...,q j-1 ), decompose the convolution kernel q into several sub-kernels, so that the length of each sub-kernel can form a convolution process of size F(3,2) with the sub-block, apply the Winograd accelerated form of F(3,2) convolution in step 2 to any sub-block and all the corresponding sub-kernels, and finally add and merge the results of different sub-kernels;

[0018] The calculation formula for convolution kernel decomposition is as follows:

[0019]

[0020] Step 5: Merge sub-blocks:

[0021] All sub-block results of the velocity or stress field p are merged to obtain the differential field of the velocity or stress field, thereby completing the differential calculation in the seismic wave field simulation.

[0022] The present invention establishes the intrinsic connection between spatial difference and convolution, and based on the Winograd fast convolution principle, derives the accelerated form of the F(3,2) convolution process, and proposes a convolution kernel decomposition strategy, which can decompose large-size convolution kernels into F(3,2) basic units for acceleration, making up for the deficiency that the current Winograd fast convolution algorithm cannot be applied to even-order convolution kernels, greatly improving the computational efficiency of convolution, thereby realizing efficient simulation of seismic wave fields.

[0023] The present invention proposes an efficient simulation method for seismic wave fields based on convolution kernel decomposition-Winograd fast convolution algorithm, and applies it to the field of seismology. BRIEF DESCRIPTION OF THE DRAWINGS

[0024] Figure 1 The overall framework flow chart of the Winograd fast convolution algorithm based on the F(3,2) kernel decomposition strategy of the embodiment;

[0025] Figure 2 This is a comparison diagram of the x-direction velocity at (2km, 2km, 2km) when the present invention and ordinary convolution are used for wave field simulation;

[0026] Figure 3 This is a comparison diagram of the velocity in the y direction at (2km, 2km, 2km) when the present invention and ordinary convolution are used for wave field simulation;

[0027] Figure 4 This is a comparison diagram of the velocity in the z direction at (2km, 2km, 2km) when the present invention and ordinary convolution are used for wave field simulation;

[0028] Figure 5 This is a comparison diagram of the x-direction velocity at (3km, 3km, 3km) when the present invention and ordinary convolution are used for wave field simulation;

[0029] Figure 6 This is a comparison diagram of the velocity in the y direction at (3km, 3km, 3km) when the present invention and ordinary convolution are used for wave field simulation;

[0030] Figure 7This is a comparison diagram of the velocity in the z direction at (3km, 3km, 3km) when the wave field simulation is performed using the present invention and ordinary convolution. DETAILED DESCRIPTION

[0031] Specific implementation method 1: The efficient simulation method of seismic wave field based on convolution kernel decomposition-Winograd fast convolution algorithm in this implementation method is implemented according to the following steps:

[0032] Step 1: Establish the intrinsic connection between spatial difference and convolution;

[0033] In the staggered grid finite difference scheme for the wave equation, the spatial difference of velocity or stress at any location is expressed as:

[0034]

[0035] where f n ′ represents the velocity or stress differential value at n points in the staggered grid, c -i represents the convolution kernel, that is, the differential coefficient, f n-i represents the velocity or stress at position ni in the staggered grid, and M represents half of the difference order;

[0036] Step 2: Establish the Winograd acceleration form of F(3,2) convolution:

[0037] Let F(3,2) convolution refer to the convolution process with a convolution kernel length of 2 and a convolution output length of 3;

[0038] The convolution calculation process of the signal d(d0,d1,d2,d3) and the convolution kernel g(g0,g1) is d(d0,d1,d2,d3)*g(g0,g1)=r(r0,r1,r2), where r0=g0d0+d1g1, r1=g0d1+g1d2, r2=g0d2+g1d3;

[0039] Suppose vector m(m1,m2,m3,m4), if m1=(d0-d2)g0,m2=(d1+d2)(g0+g1) / 2,m3=(d2-d1)(g0-g1) / 2,m4=(d3-d1)g1, then r0=m1+m2+m3,r1=m2-m3,r3=m2+m3+m4, the Winograd acceleration form of F(3,2) convolution is established;

[0040] Step 3: Velocity or stress field segmentation:

[0041] Construct seismic wave fields. For any length of velocity field or stress field p(p0,p1,...,p i-1 ) and the convolution kernel q(q0,q1,...,q j-1), first decompose the velocity field or stress field p into several sub-blocks, so that the length of each sub-block can form a convolution process of size F(3,j) with the convolution kernel q; the length of each sub-block after p is divided into blocks is 3+j-1, the step size is 3, and there is overlap between each sub-block, and the overlap length is j-1;

[0042] Step 4: Divide the convolution kernel into blocks and apply the Winograd acceleration form of F(3,2) convolution:

[0043] For the convolution kernel q(q0,q1,...,q j-1 ), decompose the convolution kernel q into several sub-kernels, so that the length of each sub-kernel can form a convolution process of size F(3,2) with the sub-block, apply the Winograd accelerated form of F(3,2) convolution in step 2 to any sub-block and all the corresponding sub-kernels, and finally add and merge the results of different sub-kernels;

[0044] The calculation formula for convolution kernel decomposition is as follows:

[0045]

[0046] Step 5: Merge sub-blocks:

[0047] All sub-block results of the velocity or stress field p are merged to obtain the differential field of the velocity or stress field, thereby completing the differential calculation in the seismic wave field simulation.

[0048] In step 1 of this embodiment, in the staggered grid finite difference format of the wave equation, taking the fourth order as an example, the spatial difference of velocity or stress at any position (0, 1, 2... represents the ordinal number of the spatial point) can be expressed as:

[0049] f′0=c -2 f -2 +c -1 f -1 +c1f1+c2f2

[0050] f′1=c -2 f -1 +c -1 f0+c1f2+c2f3

[0051] f′2=c -2 f0+c -1 f1+c1f3+c2f4 ...

[0053] where f i The value of velocity or stress at point i, f i ′ is the differential value of velocity or stress at point i, c -2-c2 is the differential coefficient, so it is equivalent to a convolution process f′=c*f. -1 and f -2 represents the two points on the left side of point 0, and f1 and f2 represent the two points on the right side of point 0.

[0054] Specific implementation method 2: The difference between this implementation method and specific implementation method 1 is that the calculation formula of the velocity field or stress field p and the convolution kernel q after the block division in step 3 is as follows:

[0055]

[0056] Specific implementation method three: The difference between this implementation method and specific implementation method one or two is that the length of the convolution kernel q in step three is 4 to 10.

[0057] Specific implementation method 4: The difference between this implementation method and specific implementation method 3 is that the length of the convolution kernel q in step 3 is 6 to 8.

[0058] Specific implementation method five: This implementation method is different from specific implementation methods one to four in that the seismic wave field constructed in step three is a full-space wave field.

[0059] Specific implementation method 6: This implementation method is different from specific implementation method 5 in that the grid size of the seismic wave field in step 3 is 50-200 m.

[0060] Specific implementation method seven: The difference between this implementation method and any one of specific implementation methods one to six is ​​that in step four, when the convolution kernel q is decomposed into several sub-kernels, there is no overlapping part, the step length is 2, and the length of each sub-kernel after decomposition is 2.

[0061] Specific implementation eight: The difference between this implementation and any one of specific implementations one to seven is that the merging method when merging sub-blocks in step five is to merge them by matrix connection.

[0062] The method of merging sub-blocks in this implementation is not to add the elements in a certain position, but to connect the elements in the sub-blocks end to end and merge them into a total matrix.

[0063] Embodiment: The efficient simulation method of seismic wave field based on convolution kernel decomposition-Winograd fast convolution algorithm of this embodiment is implemented according to the following steps:

[0064] Step 1: Establish the intrinsic connection between spatial difference and convolution;

[0065] In the staggered grid finite difference scheme for the wave equation, the spatial difference of velocity or stress at any location is expressed as:

[0066]

[0067] where f n′ represents the velocity or stress differential value at n points in the staggered grid, c i represents the convolution kernel, that is, the differential coefficient, f n-i represents the velocity or stress at position ni in the staggered grid, and M represents half of the difference order;

[0068] Step 2: Establish the Winograd acceleration form of F(3,2) convolution:

[0069] Let F(3,2) convolution refer to the convolution process with a convolution kernel length of 2 and a convolution output length of 3;

[0070] The convolution calculation process of the signal d(d0,d1,d2,d3) and the convolution kernel g(g0,g1) is d(d0,d1,d2,d3)*g(g0,g1)=r(r0,r1,r2), where r0=g0d0+d1g1, r1=g0d1+g1d2, r2=g0d2+g1d3, and a total of 6 multiplications are required;

[0071] For vector m(m1,m2,m3,m4), if m1=(d0-d2)g0,m2=(d1+d2)(g0+g1) / 2,m3=(d2-d1)(g0-g1) / 2,m4=(d3-d1)g1,then r0=m1+m2+m3,r1=m2-m3,r3=m2+m3+m4. Since the convolution kernel is known during the convolution process, (g0+g1) / 2 in m2 and (g0-g1) / 2 in m3 can be quickly calculated, and the Winograd acceleration form of F(3,2) convolution is established, which only requires 4 multiplications.

[0072] Step 3: Velocity field segmentation:

[0073] To construct the seismic wave field, the three-dimensional grid is configured as 60×60×60, and the grid size is 100m. For any length of velocity field p(p0,p1,...,p i-1 ) and the convolution kernel q(q0,q1,...,q j-1 ), in this embodiment, j=8, first decompose p into several sub-blocks, so that the length of each sub-block can form a convolution process of size F(3,j) with the convolution kernel q; the length of each sub-block after the velocity field p is divided into blocks is 3+j-1, the step size during the division is 3, and there is overlap between each sub-block, and the overlap length is j-1; the calculation formula of the velocity field p and the convolution kernel q after the division is as follows:

[0074]

[0075] Step 4: Divide the convolution kernel into blocks and apply the Winograd acceleration form of F(3,2) convolution:

[0076] For the convolution kernel q(q0,q1,...,q j-1 ), decompose the convolution kernel q into several sub-kernels, so that the length of each sub-kernel can form a convolution process of size F(3,2) with the sub-block, and apply the Winograd acceleration form of F(3,2) convolution to any sub-block and all the corresponding sub-kernels, that is, the velocity field p is equivalent to the vector m in step 2. After the convolution kernel q is divided into blocks, each sub-block is equivalent to the convolution kernel g in step 2, and finally the results of different sub-kernels are added and merged; there is no overlapping part when the convolution kernel q is decomposed into several sub-kernels, the step size is 2, and the length of each sub-kernel after decomposition is 2;

[0077] The calculation formula for convolution kernel decomposition is as follows:

[0078]

[0079] Step 5: Merge sub-blocks:

[0080] The convolution results of all sub-blocks of the velocity field p are combined through matrix connection operations to obtain the differential field of the velocity field, thereby completing the differential calculation in the seismic wave field simulation.

[0081] Step 6: Result verification:

[0082] The results of the present invention are compared with the performance of the common convolution process. A full-space wave field simulation is designed to compare with the common convolution process, with a convolution kernel length of 8 (eighth-order difference), a three-dimensional grid number configuration of 60×60×60, and a grid size of 100m.

[0083] The fast convolution algorithm proposed in this invention has obvious advantages over ordinary convolution. Figure 3-7 The speed comparison diagram of three directions at different positions when the present invention and ordinary convolution are used for wave field simulation is shown. It can be seen that compared with ordinary convolution, the present invention does not produce any accuracy loss. Table 1 is a speed comparison when the fast convolution proposed by the present invention and ordinary convolution are used for wave field simulation. It can be seen that the time taken by the present invention is only 1 / 2.55 of that of ordinary convolution, and the efficiency is significantly improved.

[0084] Table 1 Comparison of the time taken to perform wave field simulation using the present invention and using ordinary convolution

[0085] .

Claims

1. A seismic wave field simulation method based on convolution kernel decomposition-Winograd fast convolution algorithm, characterized by The efficient seismic wave field simulation method is implemented according to the following steps: Step 1: Establish the intrinsic connection between spatial difference and convolution; In the staggered grid finite difference scheme for the wave equation, the spatial difference of velocity or stress at any location is expressed as: where f n ′ represents the velocity or stress differential value at n points in the staggered grid, c -i represents the convolution kernel, that is, the differential coefficient, f n-i represents the velocity or stress at position ni in the staggered grid, and M represents half of the difference order; Step 2: Establish the Winograd acceleration form of F(3,2) convolution: Let F(3,2) convolution refer to the convolution process with a convolution kernel length of 2 and a convolution output length of 3; The convolution calculation process of the signal d(d0,d1,d2,d3) and the convolution kernel g(g0,g1) is d(d0,d1,d2,d3)*g(g0,g1)=r(r0,r1,r2), where r0=g0d0+d1g1, r1=g0d1+g1d2, r2=g0d2+g1d3; Suppose vector m(m1,m2,m3,m4), if m1=(d0-d2)g0,m2=(d1+d2)(g0+g1) / 2,m3=(d2-d1)(g0-g1) / 2,m4=(d3-d1)g1, then r0=m1+m2+m3,r1=m2-m3,r3=m2+m3+m4, the Winograd acceleration form of F(3,2) convolution is established; Step 3: Velocity or stress field segmentation: Construct seismic wave fields. For any length of velocity field or stress field p(p0,p1,...,p i-1 ) and the convolution kernel q(q0,q1,...,q j-1 ), first decompose the velocity field or stress field p into several sub-blocks, so that the length of each sub-block can form a convolution process of size F(3,j) with the convolution kernel q; the length of each sub-block after p is divided into blocks is 3+j-1, the step size is 3, and there is overlap between each sub-block, and the overlap length is j-1; Step 4: Divide the convolution kernel into blocks and apply the Winograd acceleration form of F(3,2) convolution: For the convolution kernel q(q0,q1,...,q j-1 ), decompose the convolution kernel q into several sub-kernels, so that the length of each sub-kernel can form a convolution process of size F(3,2) with the sub-block, apply the Winograd accelerated form of F(3,2) convolution in step 2 to any sub-block and all the corresponding sub-kernels, and finally add and merge the results of different sub-kernels; The calculation formula for convolution kernel decomposition is as follows: Step 5: Merge sub-blocks: All sub-block results of the velocity or stress field p are merged to obtain the differential field of the velocity or stress field, thereby completing the differential calculation in the seismic wave field simulation.

2. The seismic wave field simulation method based on convolution kernel decomposition-Winograd fast convolution algorithm according to claim 1 is characterized in that The calculation formulas of the velocity field or stress field p and the convolution kernel q after segmentation in step 3 are as follows:

3. The seismic wave field simulation method based on convolution kernel decomposition-Winograd fast convolution algorithm according to claim 1 is characterized in that The length of the convolution kernel q in step 3 is 4 to 10.

4. The seismic wave field simulation method based on convolution kernel decomposition-Winograd fast convolution algorithm according to claim 3 is characterized in that The length of the convolution kernel q in step 3 is 6 to 8.

5. The seismic wave field simulation method based on convolution kernel decomposition-Winograd fast convolution algorithm according to claim 1 is characterized in that The seismic wave field constructed in step three is a full-space wave field.

6. The seismic wave field simulation method based on convolution kernel decomposition-Winograd fast convolution algorithm according to claim 5 is characterized in that The grid size of the seismic wave field in step 3 is 50 to 200 m.

7. The seismic wave field simulation method based on convolution kernel decomposition-Winograd fast convolution algorithm according to claim 1 is characterized in that In step 4, when the convolution kernel q is decomposed into several sub-kernels, there is no overlapping part, the step size is 2, and the length of each sub-kernel after decomposition is 2.

8. The seismic wave field simulation method based on convolution kernel decomposition-Winograd fast convolution algorithm according to claim 1, characterized in that The merging method when merging sub-blocks in step 5 is to merge them through matrix connection.

Images