Missing seismic traces are inevitable in marine seismic data acquisition due to various objective factors. Cable towing issues, for instance, can arise from sudden shifts in ocean currents that pull the cables off their intended paths, or from contact with uneven seabed features that cause temporary disruptions in data transmission. Ocean currents themselves, with their varying velocities and directional changes, exert continuous forces on the acquisition equipment, leading to intermittent loss of signal reception. Weather conditions such as storms or strong winds not only affect the stability of the survey vessel but also create background noise that interferes with seismic wave detection, resulting in missing data points. Additionally, restricted mining areas, where access is limited by environmental regulations, territorial boundaries, or existing infrastructure, prohibit the placement of sensors in critical locations, leaving gaps in the data coverage. These missing traces can significantly compromise subsequent data analysis: reservoir characterization may misjudge the size and distribution of potential hydrocarbon reservoirs, and structural interpretation might fail to identify key geological features like faults or folds. Existing reconstruction methods, mostly based on single-algorithm frameworks, struggle to address these issues effectively. Some rely solely on wavelet-based algorithms, which lack the ability to capture directional features of seismic events, while others use Fourier transform methods that perform poorly with non-stationary data. These single-algorithm approaches often suffer from drawbacks such as slow convergence, where hundreds of iterations may be needed to achieve a stable result, or high computational costs that make them impractical for processing the massive datasets generated by modern high-resolution seismic surveys.

To address this issue, this study proposes a joint reconstruction method that combines the accelerated linear Bregman method （ALBM） and the iterative shrinkage-thresholding algorithm （ISTA）, with Curvelet transform employed as the sparse basis. The choice of Curvelet transform is rooted in its unique set of properties that make it ideal for seismic data processing. Its multi-scale nature allows it to decompose seismic signals into different frequency bands, each corresponding to geological features of varying sizes: from large-scale tectonic structures to small-scale reservoir heterogeneities. Directionality enables it to capture seismic events propagating in different directions, such as dipping reflectors that indicate subsurface stratigraphic layers. Sparsity ensures that complex seismic signals can be represented with a relatively small number of non-zero coefficients, reducing computational load, while anisotropy makes it highly effective at capturing curve and edge features, such as the boundaries of faults or reservoir edges, with far fewer coefficients than alternative transforms like the wavelet transform.

The joint algorithm is designed to leverage the complementary strengths of ALBM and ISTA in a structured workflow. In the early iterations, ALBM takes the lead by retaining a larger proportion of unthresholded Curvelet coefficients. These coefficients, which have not been filtered by thresholding operations, contain most of the energy from the main seismic events. By preserving them, ALBM accelerates the convergence process, allowing the algorithm to quickly approximate the general structure of the missing data. As the iterations progress, the algorithm shifts to relying more on ISTA. In the later stages, the remaining signals to be recovered are often weak valid signals, such as reflections from deep reservoirs, which were overshadowed by stronger events in the early iterations. ISTA excels here by applying a series of shrinkage and thresholding operations to the Curvelet coefficients, gradually isolating and enhancing these weak signals, thereby improving the overall reconstruction accuracy.

A key innovation of this method is the development of a new acceleration factor that dynamically adjusts from one to two as the iterations proceed. This factor directly controls the proportion of unthresholded Curvelet coefficients retained at each step. In the early iterations, when the factor is closer to one, a higher percentage of unthresholded coefficients is kept, ensuring that the main seismic events are captured quickly. As iterations continue and the factor increases toward two, the proportion of unthresholded coefficients is reduced, allowing ISTA to focus on refining the signal by removing noise and recovering weak features. This dynamic adjustment stands in contrast to traditional linear acceleration factors, which change at a fixed rate and often retain too few coefficients in early stages, and typical acceleration factors that lack flexibility, leading to either slow convergence or loss of weak signals. By enhancing the retention of effective signals in early iterations, the new acceleration factor significantly boosts both convergence speed and reconstruction precision.

Theoretical simulations were conducted under three representative scenarios to validate the method’s superiority: non-uniform undersampling, uniform undersampling, and noisy conditions. In the non-uniform undersampling test, where 40% of traces were randomly missing, a common scenario in real acquisitions, the proposed method achieved a signal-to-noise ratio （SNR） of 21.5 dB in just 20 iterations. This outperformed traditional linear acceleration factor methods, which required 35 iterations to reach 21.2 dB, and typical acceleration factor methods that took 30 iterations to achieve 21.0 dB. It also showed better alias suppression, with a 30% reduction in alias energy compared to linear factors, as measured by spectral analysis, and lower reconstruction errors （root mean square error of 0.021 versus 0.035 for linear factors）. Under uniform undersampling with 30% missing traces, simulating systematic equipment failures, the method reached an SNR of 19.8 dB in 18 iterations, preserving the continuity of seismic events more effectively than linear （28 iterations, 19.5 dB） and typical （24 iterations, 19.3 dB） factors. In noisy conditions, where Gaussian noise with a variance of 0.01 was added to 40% undersampling data, the method maintained an SNR of 18.2 dB, demonstrating strong noise resistance compared to linear （16.5 dB） and typical （17.1 dB） factors.

Application to real seismic data from deep-sea block A yielded promising results. Block A is located in a geologically complex region with water depth ranging from 1 500 to 3 000 meters. The seismic data from this area is characterized by strong multiples, which owns the unwanted reflections from the sea surface and seabed weak reflections from deep reservoirs located over 4000 meters below the seabed, and a network of complex faults that complicate signal interpretation. The dataset used in the application had a sampling rate of 4 milliseconds, a trace interval of 25 meters, and a total of 2 000 traces, covering an area of approximately 50 square kilometers. For 50% randomly undersampling data with 16 consecutive missing traces, which is a challenging scenario that often breaks the continuity of seismic events, the method reconstructed the missing information effectively, restoring the continuity of reflectors across the gaps.

Compared with methods using linear and conventional acceleration factors, the proposed method showed clear advantages. After 15 iterations, it achieved an SNR of 13.95 dB, significantly higher than the 12.30 dB of the linear acceleration factor method and 12.85 dB of the typical acceleration factor method. Energy preservation analysis revealed that it retained over 90% of the energy in the main frequency band （10-60 Hz）, which is critical for identifying subsurface features, compared to 82% for the linear method and 85% for the typical method. It also preserved weak deep reservoir reflections more effectively, with visual inspections showing clearer and more continuous weak signals in the reconstructed sections. In terms of computational efficiency, the method reached an SNR of 13.5 dB in just 73 seconds, saving nearly two-thirds of the time compared to the typical acceleration factor method, which required 220 seconds, and outperforming the 185 seconds of the linear factor method.

In summary, the proposed method, with its new acceleration factor, achieves an optimal balance between reconstruction speed and precision. It has demonstrated its ability to adapt to complex scenarios such as uniform missing, non-uniform missing, and noise interference that all common in marine seismic acquisition. By meeting the dual requirements of high efficiency and high precision, it provides a practical solution for processing large-scale marine seismic datasets. Future work will focus on optimizing the acceleration factor’s adjustment strategy to handle even more complex geological conditions, such as areas with severe faulting, and to explore its application in 3D seismic data reconstruction, where the data volume is exponentially larger and the challenges are more pronounced.