Application of an improved particle filter for random seismic noise suppression

Random noise is inevitable during seismic prospecting. Seismic signals, which are variable in time and space, are damaged by conventional random noise suppression methods, and this limits the accuracy in seismic data imaging. In this paper, an improved particle filtering strategy based on the firefly algorithm is proposed to suppress seismic noise. To address particle degradation problems during the particle filter resampling process, this method introduces a firefly algorithm that moves the particles distributed at the tail of the probability to the high-likelihood area, thereby improving the particle quality and performance of the algorithm. Finally, this method allows the particles to carry adequate seismic information, thereby enhancing the accuracy of the estimation. Synthetic and field experiments indicate that this method can effectively suppress random seismic noise.

The particle filter (PF) is a nonlinear filtering method that combines Monte Carlo and recursive Bayesian estimation. Unlike the Kalman filter, the PF can be well adapted to the state of a nonlinear Gaussian system. A common problem that occurs because of the traditional PF algorithm is particle degradation. When dealing with seismic data, most particles in the low-weight state are discarded during the resampling process. Only a small portion of the particles contain effective seismic information and they cannot fully reflect an effective signal.
In this study, particle filtering was combined with the firefly algorithm to reduce the impact of particle degradation on noise suppression. Using this algorithm, particles in the low-weight state are iteratively optimised and driven to the high-likelihood area. This improves the overall quality of the particle sample so that a more accurate state estimation can be obtained, consequently clarifying the seismic information better.

Particle filter
Particle filtering is based on the Bayesian theoretical framework (Zhu & Xu 2014;Tian et al. 2015). Compared with the Kalman filter, the PF is not limited to a linear Gaussian distribution; therefore, more distributions can be implemented. During seismic data noise reduction processing, the measured seismic data are regarded as observational data superimposed with noise. Therefore, there is a need to establish a suitable state-space model and perform filtering processing using a PF algorithm.
Suppose the state-space function model describing the nonlinear dynamic system is: where x k indicates the effective signal of the system at time k, y k is the observed noisy signal at time k; f k : process noise of the system and u k is the observational noise. The primary idea of particle filtering is to express the actual state of the system through a set of random weighted samples {x i k , w i k } N i = 1 according to the following workflow: (i) Initialisation: assuming that the initial probability density of the system is p(x), a set of random initial samples are generated from it, (x i 0 , i = 1, 2, 3, … … , N) (ii) Importance sampling (Tian et al. 2015): N particles are extracted from the importance sampling density function x i k ∼ q(x k |x i k−m:k−1 , y k ) to obtain a new particle set, iii) Weight calculation and normalisation processing: (iv) State estimation:x During the filtering process, the state and weight of the particles are updated according to the latest observations. After several iterations, the variance in the particle weights may become very large. Except for a few particles, most particles have minimal or even zero weights, causing information loss. Substantial computing resources will thus be consumed (and wasted) on these insignificant particles if iteration is continued.

Firefly optimisation algorithm
The emergence of the firefly algorithm stems from observing and summarising the behavior of firefly populations in nature (Yang 2010). The algorithm uses the brightness and attraction differences of individual fireflies to update the position for iterative optimisation. The specific optimisation mechanism is as follows: The relative brightness and attractiveness of fireflies are defined as: where I 0 is the maximum brightness of the firefly, 0 is the maximum attraction, is the light intensity absorption coefficient and r ij is the spatial distance between firefly i and firefly j. The updated equation for the position where firefly i is attracted by firefly j and moves to it is where x i and x j are the spatial positions of fireflies i and j, respectively, [0, 1] is the step size factor and rand is a random factor that obeys a uniform distribution [0,1].

PF optimisation (FAPF)
The swarm intelligence optimisation algorithm is the current development direction for improving the PF (Fang et al. 2007;Qiao 2014;Tian et al. 2015;Fang et al. 2019). The FAPF algorithm replaces the standard PF mechanism of copying high-weight particles with the firefly-optimised particle spatial distribution, which fundamentally solves the particle degradation problem and improves the diversity of particles.
(i) Defining the brightness formula for a firefly In the firefly algorithm, individuals with high brightness attract other fireflies. Here, the fluorescence brightness of the firefly is defined as the normalised weight of the particles as follows: where x i t represents particle x i at time t, andw(x i t ) represents the normalised weight of particle x i at time t. The normalised weights of the particles reflect the characteristics of the individual firefly fluorescence. During the optimisation process, 944 the particle swarm can be driven to an area with a significant weight to solve the particle degradation problem.
(ii) Defining the location update formula In the standard firefly algorithm, the step factor is typically fixed. When the algorithm iteration reaches a certain level, all particles gather near the best point and oscillate repeatedly. This phenomenon results in a wastage of computing resources. Therefore, in this study, a new adaptive adjustment strategy for the step factor was defined as: where k is the parameter that controls the rate of step reduction,w(x t ) is the weight of particle i set at time t and w max (x t ) is the maximum weight of the particles set at time t, which makes and the normalised weight of the particle inversely proportional. A significant weight indicates that the distance is close to the optimal solution, so a small step size is used. Conversely, a small weight indicates that the distance is very different from that of the optimal solution, and a large step size is used. The position update equation of the revised firefly algorithm is: where g represents the current number of iterations, x high t (g) represents the average value of particles in a high-likelihood region at time t and the step size factor is shown in equation (9). In addition, by setting a threshold, when the minimum weight of all particlesw min (x t ) ≤ Threshold, the iteration is stopped; otherwise, further iterations are performed until the maximum number of times, and resampling is conducted. Choosing an appropriate number of iterations or setting a threshold can improve the coverage of particles near large weights and increase their diversity. The filtering process is shown in figure 1 and is described as follows: (i) Seismic data input. (ii) Modeling: the state-space model of the system is established, including system equations, measurement equations and data length. (iii) Initialisation: N particles are sampled from the initial probability density; the weight of each particle is set as 1/N to obtain the particle set, cent brightness of the fireflies is considered to be the normalised weight of the particles. (vi) A threshold and a maximum number of iterations are set. The adaptive step size firefly algorithm is used to move particles with small weights (low brightness) to high-likelihood regions (high brightness), the normalised weights of the updated particles are calculated, and the particle brightness is updated simultaneously. (vii) During the iterations, the high-likelihood particles are continuously calculated and updated. The iteration is stopped when the most negligible weightw min (x t ) > Threshold, otherwise the iteration is continued to the maximum number of iterations and this portion of the particles is updated.

Simulation
To verify the effectiveness of the FAPF algorithm proposed in this study, a commonly used nonlinear system was selected as the model for comparison experiments. The state equation and observation equation were as follows. State equation: Observation equation: The process noise u and the observation noise v are independent of each other, but both obey the Gaussian distribution. Different particle numbers, N = 20, 50 and 100, were tested for comparison. The results are shown in figure 2. Figure 2 parts a, c and e show that the overall estimation accuracy can be maintained through both methods. Compared with the standard PF algorithm, the FAPF algorithm has a higher estimation accuracy because the iterative optimisation of the particles is more reasonable in terms of the 946 spatial distribution and state update. As shown in figure 2 parts b, d and e, the error level of the FAPF algorithm is lower than that of the standard PF algorithm. In addition, by comparing the average root mean square error of 50 experiments, as shown in Table 1, it can be seen that the FAPF algorithm obtained a better estimation accuracy with a smaller number of particles.

Numerical simulation
In this study, a single-channel convolution record of the Ricker wavelet was used for numerical simulation. The wave function is where f 0 represents the center frequency and is equal to 30 Hz. It can be seen from figure 3a and b that the use of PF and FAPF to denoise the noisy seismic convolution records has an excellent overall effect and improves the signal-to-noise ratio. The statistical characteristics were calculated using the average absolute error, variance and root mean square error (RMSE). Tables 2 and 3 present the results.
It can be seen from the tables that the mean absolute errorē, variance s 2 and RMSE of the FAPF were slightly lower than those of the PF algorithm. The signal-to-noise ratio became higher after FAPF filtering. As particle filtering depends on the model, the result of one filtering was modeled and filtered again, and the concept of multistage filtering was proposed. Figure 3 parts c and d show the experimental results. It can be seen that the FAPF method restored the waveform of the seismic record better and had a better noise reduction performance. The result of FAPF filtering was smoother, and the time-varying characteristics of the signal were better restored. As shown in figure 3e and f, random noise was nearly flat across the entire system frequency band. High-frequency noise remained in the spectrum after particle filtering. In contrast, in the spectrum after FAPF filtering, low-frequency signals were highlighted and restored, and high-frequency noise was more thoroughly suppressed. It can be seen in figure 3d that as the filter stages increased, the increasing trend of the signal-to-noise ratio gradually eased. It can be concluded that the selection of an appropriate filter stage is necessary to obtain the best filtering effect depending on the actual situation.
To further verify the effectiveness of the method used in this study, a large amount of random noise was added to the standard poststack seismic data model (sigmoid model). The model included a unidirectional sloping stratum at the top, a syncline structure in the middle, an anticline structure at the bottom, an angular unconformity contact between the strata and a normal fault. The time depth was set at 1000 sampling points, and the sampling interval was 4 ms. The line length was set at 200 tracks and the track spacing was 0.01 km. The experimental results are presented next. Figure 4c shows the results of the EMD method. Although the profile information was highlighted to a certain extent, random noise was not completely removed from the filtered profile. As the traditional EMD method decomposes the signal into multiple IMF components, it only suppresses random noise by discarding the first IMF components. Although this effectively suppresses high-frequency noise, it cannot separate low-frequency information from noise. The difference profile (figure 4d) shows that the noise suppressed by this method is relatively low. Figure 4e shows the results of the FXDECON method. The filtered profile still contained a small amount of random noise. The difference profile in figure 4f shows that this method eliminated part of the structure information while suppressing noise. The protection of structural information was poor, especially for protecting curved event shafts and faults. When data are processed by the FXDECON method, the data are often divided into small windows to meet the assumption of in-phase axis linearity. As such, this method adaptively selects the processing window size. Figure 4g shows the results of the PF method. It can be seen from the figure that the noise was effectively suppressed, but figure 4h shows that this method has a poor ability to protect effective signals. Figure 4i shows the results of the FAPF method proposed in this paper. The profile was clearer, and random noise was substantially suppressed. From the difference profile in figure 4j, it can be seen that this method protected structural information well.

Actual data processing
The proposed method was applied to actual singleshot records and poststack profile data for a theoretical analysis.
The results of EMD, FXDECON and FAPF processing are shown in figure 5. It can be seen from the figure that both the FXDECON and FAPF methods suppressed a large amount of random noise, while the EMD method only removed a small amount of noise. From the difference profile (figure 5c and g) it can be seen that after EMD and FXDE-CON processing there was more signal loss compared to that after FAPF processing.
A fault is shown in figure 6. The processing results of the EMD method are shown in figure 6b. It can be seen that a portion of the random noise was suppressed, but a lot of noise remained and the overall clarity of the profile increased only slightly. From the difference profile (figure 6c), it can be seen that this method discarded the IMF components with high-frequency information and then reconstructed the signal, which had little effect on low-frequency information. The processing results of the FXDECON and FAPF methods are shown in figure 6d and ). It can be seen that both methods effectively suppressed the noise and enhanced the clarity of the in-phase axis. As shown in figure 6e and g, the FXDE-CON method had a poor protection effect on valid information, and part of the fault information was lost (as shown in the red box). The FAPF method protected structural information better.

Conclusion
This study proposed a nonlinear denoising method based on an improved firefly algorithm that optimised particle filtering and applied it to seismic random noise suppression. The firefly algorithm could move the particles distributed at the tail of a probability distribution to the high-likelihood area, thereby improving the particle quality and the performance of the algorithm. In addition, by solving the degradation problems, the particle diversity was increased and more particles were able to carry adequate seismic structure information. Simulation and application results demonstrated that this method can be effectively applied for noise suppression, can highlight effective signals and can improve the signal-to-noise ratio.