Multi-sampling-rate observation system optimization method based on improved simulated annealing

Abstract

The invention discloses a multi-sampling-rate observation system optimization method based on improved simulated annealing, and the method comprises the following steps: firstly, determining the total number N of shot points and detection points according to a multi-sampling-rate sampling observation system needing to be designed, and calculating mu; generating a new multi-sampling matrix according to the coordinates of the disturbance shot point or the detection point, and calculating a new mu; if the new value mu is smaller than the mu calculated by the previous iteration, assigning the new mu to the old mu; otherwise, judging whether a new solution is accepted probably; and comparing the mu with the minimum output with the mu in the whole process, taking the minimum value as the optimal value of the whole calculation process, and taking the observation system corresponding to the minimum value as the optimal multi-sampling-rate observation system at the moment. According to the method, a simulated annealing algorithm is introduced for the problems of calculation efficiency and precision, and the efficiency and precision of optimization calculation are improved by improving an annealing function and a probability function, so that a high-quality multi-sampling-rate observation system is obtained, and technical support is provided for design of observation systems in complex regions.

Description

Technical Field

[0001] The invention relates to a multi-sampling rate observation system optimization method based on improved simulated annealing, and belongs to the technical field of multi-sampling rate seismic acquisition in complex areas. Background Art

[0002] With the deepening of oil and gas exploration, seismic acquisition has become a key factor in oil and gas exploration in complex areas. However, there are still many problems to be solved in seismic acquisition in complex areas. For example, in the eastern region, a large number of living and industrial facilities seriously restrict the quality of seismic acquisition; in the western region, complex surface and underground geological conditions seriously affect the reasonable arrangement of excitation and receiving points. Obviously, whether in the eastern or western regions, the traditional regular seismic acquisition method cannot effectively avoid the acquisition difficulty area and cannot obtain high-quality seismic data. In this regard, the present invention takes the multi-sampling rate seismic acquisition method as the purpose, combines the high-precision characteristics of the curvelet transform and the simulated annealing algorithm, establishes a coherence evaluation method for the multi-sampling rate observation system, and constructs a multi-sampling rate observation system optimization design method through the Wrapping algorithm to calculate the curvelet coefficient, establish a multi-sampling rate perception matrix, solve the Gram matrix and calculate the maximum incoherent value. And by exploring the relationship between the multi-sampling rate observation system and the cooling characteristics of the annealing function, an annealing function with fast cooling in the early stage and slow cooling in the later stage is established to improve the accuracy and efficiency of the optimization calculation. Thereby, a high-quality multi-sampling rate observation system in complex areas is obtained, providing technical support for the design of observation systems in complex areas. Summary of the invention

[0003] In order to overcome the defects in the prior art, the present invention aims to provide a multi-sampling rate observation system optimization method based on improved simulated annealing.

[0004] The technical solution provided by the present invention to solve the above technical problems is: a multi-sampling rate observation system optimization method based on improved simulated annealing, comprising the following steps:

[0005] Step 1: First, according to the multi-sampling rate sampling observation system designed as needed, determine the total number of shot points and detection points N; preliminarily determine the number of rapid cooling times n, generate the initial multi-sampling rate matrix according to N, and calculate the maximum incoherent value μ;

[0006] Step 2: Generate a new multi-sampling matrix based on the coordinates of the disturbed shot point or the detection point, and calculate a new maximum incoherent value μ;

[0007] Step 3: Analyze the maximum incoherent value μ according to the improved annealing function. If the new maximum incoherent value μ is less than the maximum incoherent value μ calculated in the previous iteration, assign the new maximum incoherent value μ to the old maximum incoherent value μ; otherwise, determine whether to accept the new solution with probability according to the probability acceptance criterion.

[0008] Step 4, cool down according to the improved annealing function, and repeat the calculation process of steps 2 to 4 above; when the condition for terminating the loop is met, compare the maximum incoherent value μ with the maximum incoherent value μ of the whole process, and the minimum value between the two is the optimal value of the whole calculation process. At this time, the corresponding observation system is the optimal multi-sampling rate observation system.

[0009] A further technical solution is that the calculation process of the maximum incoherent value μ is:

[0010] A. Use the Wrapping algorithm to obtain the curvelet coefficients C1 and C2 of the multi-sampling rate sampling matrix;

[0011] B. Convert the curvelet coefficients C1 and C2 into Gram matrices G1 and G2;

[0012] C. Normalize the Gram matrices G1 and G2 to obtain normalized matrices G3 and G4 respectively;

[0013] D. Subtract the normalized matrix G3 from the normalized matrix G4 and take the absolute value. The maximum value obtained is the maximum incoherent value μ corresponding to the multi-sampling rate matrix.

[0014] A further technical solution is that the formula in step A is:

[0015] C1=fdct-wrapping(fftshift(fft(I)))

[0016] C2=Φ·C1

[0017]

[0018] Where: Φ is the multi-sampling rate matrix.

[0019] A further technical solution is that the conversion formula in step B is:

[0020]

[0021] Where: A H is the inverse of the perception matrix A, i is the number of rows in the matrix, and j is the number of columns in the matrix. A further technical solution is that the normalization formula in step C is:

[0022] G3=diag(1 / sqrt(diag(G1)))*G1*diag(1 / sqrt(diag(G1)))

[0023] G4=diag(1 / sqrt(diag(G2)))*G2*diag(1 / sqrt(diag(G2)))

[0024] Where: G1, G2 are Gram matrices; G3, G4 are normalized matrices.

[0025] A further technical solution is that the calculation formula in step D is:

[0026] μ=max(abs(G4-G3))

[0027] Where: G3 and G4 are normalized matrices.

[0028] A further technical solution is that the improved annealing function in step 3 is:

[0029]

[0030] Where: k is the number of iterations, alpha is the cooling coefficient, and T0 is the initial temperature.

[0031] A further technical solution is that the probability acceptance criterion in step 3 is:

[0032]

[0033] Where: k is the current iteration number, k max is the maximum number of iterations; F(X new ) is the new objective function value, F(X old ) is the old objective function value.

[0034] A further technical solution is that the specific process of cooling according to the improved annealing function in step 4 is:

[0035] (1) Input the initial temperature T0 and the stop temperature T min , Markov chain length iter, cooling coefficient alpha;

[0036] (2) Select the initial temperature T0, the initial solution Φ0, and calculate the objective function F(Φ);

[0037] (3) At the same temperature, generate a new neighborhood solution Φ1 based on the current solution, and calculate the difference in target value ΔF = F(Φ1) - F(Φ0);

[0038] (4) If ΔF < 0, then accept Φ1 as the new solution, otherwise accept the bad solution with probability function p = exp(-ΔF / T);

[0039] (5) If thermal equilibrium is reached, iterate iter times and go to (5) to perform the next step of calculation; otherwise, go to (2) to recalculate;

[0040] (6) Lower the temperature T k If the temperature T is lowered k<initial temperature T0, the iterative calculation ends and the optimal solution is obtained.

[0041] The present invention has the following beneficial effects: the method is based on the multi-sampling rate seismic acquisition method, establishes a multi-sampling rate observation system optimization design method, and introduces a simulated annealing algorithm to address the problems of computational efficiency and accuracy. By improving the annealing function and probability function, the efficiency and accuracy of the optimization calculation are improved, thereby obtaining a high-quality multi-sampling rate observation system, providing technical support for the design of observation systems in complex areas. BRIEF DESCRIPTION OF THE DRAWINGS

[0042] Figure 1 The cooling process of the present invention and the conventional annealing function is compared;

[0043] Figure 2 Contrast the iterative calculation convergence process of the present invention and the conventional annealing function;

[0044] Figure 3 Contrast the convergence process of iterative calculation of the present invention and conventional probability function;

[0045] Figure 4 It is a multi-sampling rate observation system obtained by calculation after the implementation of the present invention. DETAILED DESCRIPTION

[0046] The technical solution of the present invention will be described clearly and completely below in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0047] The invention discloses a multi-sampling rate observation system optimization method based on improved simulated annealing, comprising the following steps:

[0048] Step 1: First, according to the multi-sampling rate sampling observation system designed as needed, determine the total number of shot points and detection points N; preliminarily determine the number of rapid cooling times n, generate the initial multi-sampling rate matrix according to N, and calculate the maximum incoherent value μ;

[0049] The calculation process of the maximum incoherent value μ is as follows:

[0050] 1. Use the Wrapping algorithm to obtain the curve coefficients C1 and C2 of the multi-sampling rate sampling matrix;

[0051] C1=fdct-wrapping(fftshift(fft(I)))

[0052] C2=Φ·C1

[0053]

[0054] 2. Convert the curvelet coefficients C1 and C2 into Gram matrices G1 and G2;

[0055]

[0056]

[0057] Where: A H is the inverse of the perception matrix A, i is the number of rows in the matrix, and j is the number of columns in the matrix;

[0058] The construction process of the new Gram matrix is:

[0059] Multi-sampling rate theory is established.

[0060] In sparse multi-sampling theory, one of the key factors to achieve signal reconstruction is to construct the incoherence relationship between the sampling matrix and the sparse transformation basis. The mathematical model is as follows:

[0061] y=ΦX=ΦC H s=As

[0062] Where y represents the sampled data; Φ represents the multi-sampling rate matrix; X represents the complete data; C represents the curvelet transform basis; C H represents the inverse of the curvelet transform basis; s represents the sparse representation of data X in the curvelet domain; A = ΦC H represents the perception matrix;

[0063] According to the above theory, in order to successfully reconstruct incomplete data, the perception matrix A must satisfy the RIP isometry criterion (restricted isometry), that is, for any sparse vector ν, It is established. However, it is very difficult to solve the matrix that satisfies the RIP condition. Therefore, the present invention uses an equivalent condition to evaluate the incoherence between the multi-sampling rate matrix and the sparse basis, and constructs a new Gram matrix by fully considering the influence of the matrix size on the calculation results.

[0064] 3. Normalize the Gram matrices G1 and G2 to obtain normalized matrices G3 and G4 respectively;

[0065] G3=diag(1 / sqrt(diag(G1)))*G1*diag(1 / sqrt(diag(G1)))

[0066] G4=diag(1 / sqrt(diag(G2)))*G2*diag(1 / sqrt(diag(G2)))

[0067] Where: G1, G2 are Gram matrices; G3, G4 are normalized matrices;

[0068] 4. Subtract the normalized matrix G3 from the normalized matrix G4 and take the absolute value. The maximum value is the maximum incoherent value μ corresponding to the multi-sampling rate matrix;

[0069] μ=max(abs(G4-G3))

[0070] Where: G3 and G4 are normalized matrices;

[0071] Step 2: Generate a new multi-sampling matrix based on the coordinates of the disturbed shot point or the detection point, and calculate a new maximum incoherent value μ;

[0072] Step 3: Analyze the maximum incoherent value μ according to the improved annealing function. If the new maximum incoherent value μ is less than the maximum incoherent value μ calculated in the previous iteration, assign the new maximum incoherent value μ to the old maximum incoherent value μ; otherwise, determine whether to accept the new solution with probability according to the probability acceptance criterion.

[0073] Improvement of annealing function: In theory, the higher the initial temperature, the greater the probability of searching for the global optimum, but the calculation time will increase significantly. Aiming at the disadvantage that all temperatures of the conventional annealing function have the same cooling rate, the present invention proposes a segmented cooling method, each segment is assigned a different cooling rate, and rapid cooling is performed in the high temperature stage. When the temperature drops rapidly to near the optimal value (low temperature stage), slow cooling is adopted ( Figure 1 ), conduct a full search, escape the trap of local optimal solution, and improve computational efficiency and solution accuracy ( Figure 2 ).

[0074] The conventional annealing cooling function is:

[0075] T k =T0·alpha

[0076] The improved annealing function proposed by the present invention is as follows:

[0077]

[0078] Where: k is the number of iterations, alpha is the cooling coefficient, and T0 is the initial temperature.

[0079] The present invention proposes that when the number of iterations is small in the early stage of calculation, the probability of accepting a bad solution should be large to avoid falling into a local solution and enhance the global search capability; when entering the late stage of calculation, the difference in the objective function value is small, the probability of accepting a bad solution is reduced in real time, the local fine search capability is enhanced, the calculation accuracy is improved, and the global optimal solution can be obtained ( Figure 3 );

[0080] The Metropolis probability acceptance criterion of the traditional algorithm is:

[0081]

[0082] The Metropolis probability acceptance criterion proposed by the present invention is:

[0083]

[0084] Where: k is the current iteration number, k max is the maximum number of iterations; F(X new ) is the new objective function value, F(X old ) is the old objective function value;

[0085] Step 4, cool down according to the improved annealing function, and repeat the calculation process of steps 2 to 4 above; when the condition for terminating the loop is met, compare the maximum incoherent value μ with the maximum incoherent value μ of the whole process, and the minimum value of the two is the optimal value of the whole calculation process. At this time, the corresponding observation system is the optimal multi-sampling rate observation system ( Figure 4 ).

[0086] The improved annealing function is used to cool down the system. The randomly generated solution is used as the starting point. New solutions are continuously generated under the action of the temperature function. Then, the new solutions are compared with the old solutions and accepted according to the Metropolis probability criterion. The optimal solution is found in the continuous iteration, and then the temperature is reduced according to the cooling function. The above process is repeated to finally obtain the global optimal solution. The specific steps are as follows:

[0087] (1) Input the initial temperature T0 and the stop temperature T min , Markov chain length iter, cooling coefficient alpha;

[0088] (2) Select the initial temperature T0, the initial solution Φ0, and calculate the objective function F(Φ);

[0089] (3) At the same temperature, generate a new neighborhood solution Φ1 based on the current solution, and calculate the difference in target value ΔF = F(Φ1) - F(Φ0);

[0090] (4) If ΔF < 0, then accept Φ1 as the new solution, otherwise accept the bad solution with probability function p = exp(-ΔF / T);

[0091] (5) If thermal equilibrium is reached, iterate iter times and go to (5) to perform the next step of calculation; otherwise, go to (2) to recalculate;

[0092] (6) Lower the temperature T k If the temperature T is lowered k <initial temperature T0, the iterative calculation ends and the optimal solution is obtained.

[0093] The invention takes the multi-sampling rate seismic acquisition method as its purpose, utilizes the high-precision characteristics of curvelet transform and simulated annealing algorithm, and establishes a coherence evaluation optimization method for multi-sampling rate matrices.

[0094] First, the curvelet coefficient is calculated by the Wrapping algorithm; secondly, the multi-sampling rate perception matrix and Gram matrix are established, and the matrix is ​​solved. Finally, the maximum incoherent value is solved according to the column incoherence of the matrix, so as to construct a multi-sampling rate observation system design method. At the same time, in view of the low computational efficiency and accuracy of the above methods, the simulated annealing algorithm is introduced to solve the optimal maximum incoherent value. And by exploring the relationship between the multi-sampling rate observation system and the cooling characteristics of the annealing function, a new annealing function is established to form a calculation strategy of fast cooling in the early stage and slow cooling in the later stage. The calculation efficiency is improved by fast cooling in the early stage, and the trap of local optimal solution is fully searched by slow cooling in the later stage, so as to ensure the gradual convergence of the solution process of the maximum incoherent value; at the same time, a new Metropolis probability function is proposed, so that the probability function in the simulated annealing algorithm is more in line with the characteristics of the multi-sampling rate sampling matrix, thereby improving the accuracy and efficiency of the optimization calculation as a whole. In this way, a high-quality multi-sampling rate observation system in complex areas is obtained, providing technical support for the design of observation systems in complex areas.

[0095] The above description is not intended to limit the present invention in any form. Although the present invention has been disclosed through the above embodiments, it is not intended to limit the present invention. Any technician familiar with the profession can make some changes or modifications to equivalent embodiments of equivalent changes using the technical contents disclosed above without departing from the scope of the technical solution of the present invention. However, any simple modification, equivalent changes and modifications made to the above embodiments based on the technical essence of the present invention without departing from the content of the technical solution of the present invention are within the scope of the technical solution of the present invention.

Claims

1. A multi-sampling rate observation system optimization method based on improved simulated annealing, characterized in that: The following steps are involved: Step 1: First, according to the multi-sampling rate sampling observation system designed as needed, determine the total number of shot points and detection points N; preliminarily determine the number of rapid cooling times n, generate the initial multi-sampling rate matrix according to N, and calculate the maximum incoherent value μ; Step 2: Generate a new multi-sampling matrix based on the coordinates of the disturbed shot point or the detection point, and calculate a new maximum incoherent value μ; Step 3: Analyze the maximum incoherent value μ according to the improved annealing function. If the new maximum incoherent value μ is less than the maximum incoherent value μ calculated in the previous iteration, assign the new maximum incoherent value μ to the old maximum incoherent value μ; otherwise, determine whether to accept the new solution with probability according to the probability acceptance criterion. Step 4, cool down according to the improved annealing function, and repeat the calculation process of steps 2 to 4 above; when the condition for terminating the loop is met, compare the maximum incoherent value μ with the maximum incoherent value μ of the whole process, and the minimum value between the two is the optimal value of the whole calculation process. At this time, the corresponding observation system is the optimal multi-sampling rate observation system.

2. The method for optimizing a multi-sampling rate observation system based on improved simulated annealing according to claim 1, characterized in that: The calculation process of the maximum incoherent value μ is: A. Use the Wrapping algorithm to obtain the curvelet coefficients C1 and C2 of the multi-sampling rate sampling matrix; B. Convert the curvelet coefficients C1 and C2 into Gram matrices G1 and G2; C. Normalize the Gram matrices G1 and G2 to obtain normalized matrices G3 and G4 respectively; D. Subtract the normalized matrix G3 from the normalized matrix G4 and take the absolute value. The maximum value obtained is the maximum incoherent value μ corresponding to the multi-sampling rate matrix.

3. The multi-sampling rate observation system optimization method based on improved simulated annealing according to claim 2 is characterized in that: The formula in step A is: C1=fdct-wrapping(fftshift(fft(I))) C2=Φ·C1 Where: Φ is the multi-sampling rate matrix.

4. The method for optimizing a multi-sampling rate observation system based on improved simulated annealing according to claim 2, characterized in that: The conversion formula in step B is: Where: A H is the inverse of the perception matrix A, i is the number of rows in the matrix, and j is the number of columns in the matrix.

5. The method for optimizing a multi-sampling rate observation system based on improved simulated annealing according to claim 2, characterized in that: The normalization formula in step C is: G3=diag(1 / sqrt(diag(G1)))*G1*diag(1 / sqrt(diag(G1))) G4=diag(1 / sqrt(diag(G2)))*G2*diag(1 / sqrt(diag(G2))) Where: G1, G2 are Gram matrices; G3, G4 are normalized matrices.

6. The method for optimizing a multi-sampling rate observation system based on improved simulated annealing according to claim 2, characterized in that: The calculation formula in step D is: μ=max(abs(G4-G3)) Where: G3 and G4 are normalized matrices.

7. The multi-sampling rate observation system optimization method based on improved simulated annealing according to claim 1, characterized in that: The improved annealing function in step 3 is: Where: k is the number of iterations, alpha is the cooling coefficient, and T0 is the initial temperature.

8. The method for optimizing a multi-sampling rate observation system based on improved simulated annealing according to claim 1, characterized in that: The probability acceptance criterion in step 3 is: Where: k is the current iteration number, kmax is the maximum iteration number; F(X new ) is the new objective function value, F(X old ) is the old objective function value.

9. The method for optimizing a multi-sampling rate observation system based on improved simulated annealing according to claim 7, characterized in that: The specific process of cooling according to the improved annealing function in step 4 is: (1) Input the initial temperature T0 and the stop temperature T min , Markov chain length iter, cooling coefficient alpha; (2) Select the initial temperature T0, the initial solution Φ0, and calculate the objective function F(Φ); (3) At the same temperature, generate a new neighborhood solution Φ1 based on the current solution, and calculate the difference in target value ΔF = F(Φ1) - F(Φ0); (4) If ΔF < 0, then accept Φ1 as the new solution, otherwise accept the bad solution with probability function p = exp(-ΔF / T); (5) If thermal equilibrium is reached, iterate iter times and go to (5) to perform the next step of calculation; otherwise, go to (2) to recalculate; (6) Lower the temperature T k If the temperature T is lowered k <initial temperature T0, the iterative calculation ends and the optimal solution is obtained.

Images