Exact 3D seismic data reconstruction using Tubal-Alt-Min algorithm
Feng Qian1 , Quan Chen2 , Ming-Jun Su3 , Guang-Min Hu2 , Xiao-Yang Liu∗
1 Center of Information Geoscience, University of Electronic Science and Technology of China,
2 School of Rresources and Environment, University of Electronic Science and Technology of China,
3 PetroChina Research Institute of Petroleum Exploration and Development (RIPED)-Northwest, Lanzhou, China,
∗
Department of Electrical Engineering, Columbia University, NY, USA

arXiv:1704.02445v1 [cs.CE] 8 Apr 2017

SUMMARY
Data missing is an common issue in seismic data, and many
methods have been proposed to solve it. In this paper, we
present the low-tubal-rank tensor model and a novel tensor
completion algorithm to recover 3D seismic data. This is a
fast iterative algorithm, called Tubal-Alt-Min which completes
our 3D seismic data by exploiting the low-tubal-rank property
expressed as the product of two much smaller tensors. TubalAlt-Min alternates between estimating those two tensor using
least squares minimization. We evaluate its reconstruction performance both on synthetic seismic data and land data survey. The experimental results show that compared with the
tensor nuclear norm minimization algorithm, Tubal-Alt-Min
improves the reconstruction error by orders of magnitude.

INTRODUCTION
Seismic data quality is vital to various geophysical processing. However, due to the financial and physical constraints,
the real seismic survey data are usually incomplete. Seismic
data reconstruction is a complex problem, and a large number
of researchers have devoted themselves to the research of this
field. Consequently, many approaches have been proposed to
handle this problem. From the point of view of data organization, these methods can be divided into two categories.
The low dimensional based methods, such as transform-based
methods which utilize the properties of the seismic trace in an
auxiliary domain (Hennenfent et al., 2010). Liner prediction
theory use the predictability of the signal in the f-x or t-x domain (Naghizadeh and Sacchi, 2007). And methods which exploit the low-rank nature of seismic data embedded in Hankel
matrices (Oropeza and Sacchi, 2011).
Those low dimensional based methods usually ignores the spatial structure of the seismic traces. While the spatial structure coherence is very important for the seismic data completion, hence recent developments in high dimensional tensor
completion approaches exploit various tensor decomposition
model are widely used in seismic data reconstruction (Kreimer
et al., 2013; Ely et al., 2015a). Those high-order tensor decomposition approaches have became a trend for seismic data
completion, and the exist different definitions for tensor decomposition that lead to different tensor completion model,
i.e. the higher-order singular value decomposition (HOSVD)
(Kreimer and Sacchi, 2011), the tuker decomposition (Da Silva
and Herrmann, 2013), and tensor SVD (tSVD) decomposition
(Ely et al., 2015b).

In this paper, we focus on a new seismic data reconstruction
algorithm which based on low-tubal-rank tensor decomposition model possesses extremely high precision seismic data
recovery performance. Because of the high redundancy or coherence between one seismic trace to the others (Ely et al.,
2015a), we assume that the fulled sampled seismic volume has
low-tubal-rank property in the tSVD domain. Therefore, we
can solve the seismic tensor completion problem through an
alternating minimization algorithm for low-tubal-rank tensor
completion (Tubal-Alt-Min) (Liu et al., 2016a,b). We have
evaluated the performance of this approach on both synthetic
and field seismic data.

NOTATIONS
The data in seismic survey is a natural high-dimensional tensor, such as the 3D poststack seismic data which consists of
one time or frequency dimension and two spatial dimensions
corresponding to xline and inline directions.
Throughout the paper, we denote those 3D seismic tensor in
time domain by uppercase calligraphic letter, T ∈ Rm×n×k ,
and denote the frequency domain 3D seismic tensor by Te ∈
Rm×n×k correspondingly. Uppercase letter A ∈ Rm×n denotes
matrix, and lowercase boldface letter x ∈ Rn denotes vector.
Let [n] denotes the set {1, 2, . . . , n}. In addition, we introduce
an important tensor operator t-product (Kilmer et al., 2013).
t-product. The tensor-product T = X ∗ Y of X ∈ Rn1 ×n2 ×k
n2 ×n3 ×k is a tensor of size n × n × k, T (i, j, :) =
and
1
2
Pn2Y ∈ R
s=1 X (i, s, :) ∗ Y(s, j, :), for i ∈ [n1 ] and j ∈ [n3 ].
PROBLEM SETUP
From what we have stated above, we know that seismic data
comprises many traces that provide a spatio-temporal sampling of the reflected wavefield. However, caused by various
factors, such lost many important informations. Such as the
exist of reservoir, residential or any other obstacle in the seismic data acquisition areas will lead to under-sampled seismic
record. The missing traces will complicate certain data processing steps such as the prediction accuracy of underground
reservoirs. Hence, the completion step in seismic data processing is of grate significance.
In this paper, we explored the relationship between low-tubalrank and under-sampled rate firstly, and found that low-tubal
property is positively correlated with the sampling rate. Figure
1 shows the detail experiment result. Base on this work, we as-

Empirical CDF

100

SOLUTION
90

In the Tubal-Alt-Min algorithm, the target 3D seismic volume
Tb ∈ Rm×n×k can be decomposed as Tb = X ∗ Y, X ∈ Rm×r×k ,
Y ∈ Rr×n×k , and r is the target tubal-rank. With this decomposition, the problem (1) reduces to

80
70

CDF

60
50

Tb =

40
20% sampling rate
40% sampling rate
60% sampling rate
80% sampling rate
100% sampling rate

30
20
10

10

20

30

40

50

60

70

80

90

kPΩ (T ) − PΩ (X ∗ Y)k2F .

(2)

This cost function can be solved by the alternating minimization algorithm for low tubal rank tensor completion designed
by Liu. The main algorithm steps are showing as Alg 1.

0
0

argmin
X ∈Rm×r×k ,Y∈Rr×n×k

100

Top portion of singular value %

Figure 1: The empirical CDF of singular value.

Alg. 1 Tubal Alternating Minimization: AM(PΩ (T ), Ω, L, r)
Input: Observation set Ω ∈ [n] × [m] × [k] and the corresponding elements PΩ (T ) , number of iterations L, target tubal-rank
r.
1: X0 ← Initialize(PΩ (T ), Ω, r),
2: Xe0 ← fft(X0 , [], 3), TeΩ ← fft(TΩ , [], 3),
eΩ ← fft(PΩ , [], 3),
P
3: for l = 1 to L do
eΩ , Xel−1 , r ),
el ← LS Y( TeΩ , P
Y
eΩ , Y
el , r),
Xel ← LS X( TeΩ , P
4: end for
eL , [], 3),
5: X ← ifft( XeL , [], 3); Y ← ifft( Y

The key problem of Alg. 1 is the tensor least square minimization, which was solved by the providing methods in Liu’s
paper. The main ideal is to decompose (2) into n separate standard least squares minimization problem in the frequency domain. Then, we just need to solve a least square problem like
the following form each step:
b
x = argmin kb − A1 A2 xk2F .

(3)

x∈Rrk×1

PERFORMANCE EVALUATION
To evaluate the performance of the algorithm we adopted, two
commonly used evaluation criteria in seismic data completion
filed have been used for comparison - the reconstruction error
and the convergence speed.
• Reconstruction error: here we adopted the relative square error as a scale standard which defined as RSE =
kTb − T kF /kT kF .

Output: Pair of tensors ( X , Y ).
sume that the full sampled seismic data volume has low-tubalrank property and under-sampled traces will increase the tubal
rank. Therefore, the poststack 3D seismic data reconstruction can be tackled with tensor completion tools that using
the low tubal rank property. The statement above transform
to mathematical representation is that T ∈ Rm×n×k is a 3D
seismic data with tubal-rank equal to r. Then, by observing a
set Ω ⊂ Rm×n×k of T ’s elements, we get the under-sampled
seismic data TΩ . Then, our aim is to recover T from TΩ . The
tensor reconstruction problem can be formulated as following

• Convergence speed: we measure decreasing rate of the
RSE across iterations by linearly fitting the measured
RSEs.
In order to have an intuitive performance comparison, we compared with two other seismic data volume completion algorithm. The Parallel matrix factorization algorithm (PMF) (Gao
et al., 2015) and the tensor singular value decomposition based
algorithm, also called tensor nuclear norm algorithm (TNN)
(Ely et al., 2015b) for seismic reconstruction. We applied those
algorithm both on synthetic and real seismic data, the following subsections will demonstrate the detail performance comparison.

optimization function:
Tb = argmin rank(X ),
X ∈Rm×n×k

(1)

s.t. PΩ (X ) = PΩ (T ).
Here, PΩ (·) denote the projection of a tensor onto the observed
set Ω, TΩ = PΩ (T ). The Tubal-Alt-Min algorithm proposed
by Liu recently, which can complete the low tubal rank tensor with very high currency in several iterations is a perfect
approach to solve this problem.

Synthetic data example
We use a Ricker wavelet with central frequency of 40 Hz to
generate a simple 3D seismic model with two dipping planes.
The seismic data corresponds to a spatial tensor of size 64 ×
64 × 256, 256 time samples with the time sampling rate of 1
ms and 64 corresponding to inline and xline direction. Then,
through tSVD decomposition we get U, S, and V. According
to the decomposition result, we choose the first 2 tubes of S,
and make other tubes elements equate to zeros, form a new
b Then we get the tensor Tb = U ∗ Sb ∗ V which tubaltensor S.
rank is 2.

Reconstruction error

RSE in log-scale

100

PMF
TNN
Tuabl-Alt-Min

10-5

10-10

10-15

(a)
0

10

20

30

40

50

60

70

80

90

(b)

100

Sampling rate %

Trace number
0

20

40

60

0

Figure 2: Reconstruction error for the synthetic data set as a
function of the downsampling rate.
Time(ms)

50

Convergence speed

100

PMF
TNN
Tubal-Alt-Min

150

200

10-1

RSE in log-scale

100

250

10-2

(c)

Figure 4: (a) Full sampled synthetic data. (b) Under-sampled
measured data for the case when 40% of the traces in (a) was
removed. (c) The recovered traces of (b), it’s RSE< 1e − 4. (d)
The relative residual between (c) and (a), when sampling rate
equal to 40%, there almost no residual between original data
and reconstruction data.

10-3

10-4

(d)

0

10

20

30

40

50

60

70

80

Iterations

Figure 3: Convergence speed for the three algorithms on synthetic data.

Firstly, we apply Algorithm 1 to decompose the low-tuabl-rank
tensor of our data set at different sampling rate vary from 5%
to 95% and set the tubal-rank equate to 2. Using these decompositions, we generate our reconstructed data and compare the
error between the low-tubal-rank reconstruction and full sampled data. Applying two other algorithm we stated above to
reconstruct the data at the same condition. Because of all of
the three algorithms include random sampling operator, we
averaged those algorithms’ performance under 20 above experiments. Figure 2 shows the relative square error (RSE) of
the three algorithms. From these curves, we see that the performance of our algorithm is very outstanding. Even there are
only 30% sampling traces, it can also get a relative perfect reconstruction data. As the increase of sampling rate, the reconstruction error of Tubal-Alt-Min algorithm decrease rapidly.
Comparing with TNN, when the sampling rate over 50%, it
improves the recovery error by orders of magnitude. It’s recovery error also better than the other algorithm almost in the
whole sampling rate range.
Secondly, to evaluate the reconstruction error, we fixed the sa-

mpling rate to be 40% and set the RSE tolerance of the three algorithms equal to 1e-4. Then evaluate relative square error and
convergence speed for the three algorithms. Figure 3 shows
that Tubal-Alt-Min convergence at 9th iteration, TNN-ADMM
terminates at the 60-th iteration, PMF used 77-th iterations to
reach the preset threshold. The convergence speed of our algorithm is obviously better than other algorithms.
Then, we fixed the sampling rate at 40% to evaluate the completion performance on seismic data. The portion of reconstruction result shows in figure 4. From figure 4 (b), 4 (c) and 4 (d),
we observe that the missing traces are effectively recovered.
Filed data example
We also test the performance of the reconstruction method on
a land data set that was acquired to monitor a heavy oil field
in Anyue mountain, China. Our test data is a 25 × 200 × 600
tensor cut from the full-sampled area of the filed data. On this
basis, we randomly sample 50% traces in the seismic tensor,
then reconstruct it. To reconstruct the traces, we evaluated the
low-tubal-rank of the data manually. Determined by the tSVD,
we found it’s first 16 eigentubes’s l2 norm is much larger than
the rest in almost 40 times. Consequently, we set the parameter
of Alg. 1 r = 16. Figure 4 shows portion of the result cut along
inline direction. From the comparison of figure 5 (a), 5 (c) and

Columns
100

Columns
150

200

50

100

100

200

200
Time(ms)

Time(ms)

50

300

100

150

200

is that it can achieve very high recovery accuracy in several it-

300

400

400

500

500
(a)

600

600
(b)

Columns

100

150

200

50

100

100

200

200
Time(ms)

Time(ms)

50

(a)

Columns

300

150

200

300

400

400

500

500

600

100

(b)

600
(c)

(d)

Figure 6: The reconstruction result of filed data. (a) The
horizontal-slice of filed data from Anyue survey areas. (b)
Randomly sampled 40% trances from (a). (c) The data recovered by Tubal-Alt-Min algorithm, it’s RSE reached 1e-7. (d)
The data recovered by TNN whose RSE converged to 1e-3.
(c)

5 (d), the reconstruction performance is such amazing that almost no residual between 5 (a) and 5 (c) Compared with the
reconstruction result of TNN, it seems no difference. However,
the RSE gap between them almost have a few orders of magnitude. Figure 5 shows the difference between the specific details
of the seismic data traces restored by the two algorithms.

Figure 7: The reconstruction traces of filed data. (a) The
lateral-slices of filed data from Anyue survey areas. (b) Traces
recovered by TNN. (c) Traces recovered by our algorithm.
From the recovered traces in the blue rectangles, our algorithm
has an obvious performance advantages, which can recover
more details of the seismic traces.

CONCLUSION
In this article, we formulate the problem of poststack seismic
data reconstruction as an low-tubal-rank decomposition problem. The adopted method is based on the alternating minimization approach for low-tubal rank tensor completion. The
unknown low-tubal-rank tensor is parameterized as the product of two much smaller tensors with the low-tubal-rank property being automatically incorporated, and Tubal-Alt-Min alternates between estimating those two tensors using tensor least
squares minimization. The biggest advantage of our algorithm

erations. From the experimental result, it potently proved that
compared with contrast algorithms our algorithm improves the
recovery error by orders of magnitude with much better convergence speed for higher sampling rates.

ACKNOWLEDGMENTS
We would like to acknowledge financial support from the National Natural Science Foundation of China (Grant No.U1562218).

REFERENCES
Da Silva, C., and F. Herrmann, 2013, Hierarchical tucker tensor optimization-applications to 4d seismic data interpolation: Presented at the 75th EAGE Conference & Exhibition incorporating SPE EUROPEC 2013.
Ely, G., S. Aeron, N. Hao, and M. E. Kilmer, 2015a, 5d and 4d pre-stack seismic data completion using tensor nuclear norm (tnn),
in SEG Technical Program Expanded Abstracts 2013: Society of Exploration Geophysicists, 3639–3644.
Ely, G. T., S. Aeron, N. Hao, and M. E. Kilmer, 2015b, 5d seismic data completion and denoising using a novel class of tensor
decompositions: Geophysics, 80, V83–V95.
Gao, J., A. Stanton, and M. D. Sacchi, 2015, Parallel matrix factorization algorithm and its application to 5d seismic reconstruction
and denoising: Geophysics, 80, V173–V187.
Hennenfent, G., L. Fenelon, and F. J. Herrmann, 2010, Nonequispaced curvelet transform for seismic data reconstruction: A
sparsity-promoting approach: Geophysics, 75, WB203–WB210.
Kilmer, M. E., K. Braman, N. Hao, and R. C. Hoover, 2013, Third-order tensors as operators on matrices: A theoretical and
computational framework with applications in imaging: SIAM Journal on Matrix Analysis and Applications, 34, 148–172.
Kreimer, N., and M. D. Sacchi, 2011, A tensor higher-order singular value decomposition (hosvd) for pre-stack simultaneous
noise-reduction and interpolation: Seg Technical Program Expanded Abstracts, 3069–3074.
Kreimer, N., A. Stanton, and M. D. Sacchi, 2013, Tensor completion based on nuclear norm minimization for 5d seismic data
reconstruction: Geophysics, 78, V273–V284.
Liu, X.-Y., S. Aeron, V. Aggarwal, and X. Wang, 2016a, Low-tubal-rank tensor completion using alternating minimization: SPIE
Defense+ Security, International Society for Optics and Photonics, 984809–984809.
——–, 2016b, Low-tubal-rank tensor completion using alternating minimization: 9848, 984809.
Naghizadeh, M., and M. D. Sacchi, 2007, Multistep autoregressive reconstruction of seismic records: Geophysics, 72, V111–V118.
Oropeza, V., and M. Sacchi, 2011, Simultaneous seismic data denoising and reconstruction via multichannel singular spectrum
analysis: Geophysics: Geophysics, 76, V25–V32.

