# Unsupervised Seismic Footprint Removal With Physical Prior Augmented Deep
Autoencoder

Feng Qian, Yuehua Yue, Yu He, Hongtao Yu, Yingjie Zhou,
Jinliang Tang, and Guangmin Hu Feng Qian is with the School of Information and
Communication Engineering and the Center for Information Geoscience,
University of Electronic Science and Technology of China, Chengdu, 611731
China (e-mail: fengqian@uestc.edu.cn).Yuehua Yue, Yu He, and Guangmin Hu are
all with the School of Resource and Environment and the Center for Information
Geoscience, University of Electronic Science and Technology of China, Chengdu,
611731 China (e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>hgm@uestc.edu.cn).Yingjie Zhou is with the College of Computer Science,
Sichuan University (SCU), Chengdu 610065, China (e-mail:
yjzhou@scu.edu.cn).Jinliang Tang is with the Sinopec Geophysical Research
Institute (SGRI), Institute of Seismic Imaging Technology, Nanjing 211103,
China (e-mail: tangjl.swty@sinopec.com).

###### Abstract

Seismic acquisition footprints appear as stably faint and dim structures and
emerge fully spatially coherent, causing inevitable damage to useful signals
during the suppression process. Various footprint removal methods, including
filtering and sparse representation (SR), have been reported to attain
promising results for surmounting this challenge. However, these methods,
e.g., SR, rely solely on the handcrafted image priors of useful signals, which
is sometimes an unreasonable demand if complex geological structures are
contained in the given seismic data. As an alternative, this article proposes
a footprint removal network (dubbed FR-Net) for the unsupervised suppression
of acquired footprints without any assumptions regarding valuable signals. The
key to the FR-Net is to design a unidirectional total variation (UTV) model
for footprint acquisition according to the intrinsically directional property
of noise. By strongly regularizing a deep convolutional autoencoder (DCAE)
using the UTV model, our FR-Net transforms the DCAE from an entirely data-
driven model to a prior-augmented approach, inheriting the superiority of the
DCAE and our footprint model. Subsequently, the complete separation of the
footprint noise and useful signals is projected in an unsupervised manner,
specifically by optimizing the FR-Net via the backpropagation (BP) algorithm.
We provide qualitative and quantitative evaluations conducted on three
synthetic and field datasets, demonstrating that our FR-Net surpasses the
previous state-of-the-art (SOTA) methods.

###### Index Terms:

Seismic acquisition footprint, unidirectional total variation (UTV), deep
convolutional autoencoder (DCAE), footprint removal network (FR-Net), prior-
augmented.

##  I Introduction

Acquisition footprints are extensively perceived as anomalous and periodic
amplitude stripes that are related to the acquisition geometry rather than
geology, always appearing to cross the time slices or horizon slices of 3-D
seismic data [1]. The leading causes of footprint noise are diverse and occur
throughout the data acquisition and processing procedures, but two crucial
aspects, including acquisition design and accentuation during processing, may
play a role [2]–[4]. The presence of footprint noise makes seismic images
unclear and may cause confusion when identifying and analyzing the true
subsurface response, which undoubtedly boosts uncertainty in seismic image
interpretation [5], [6]. As a result, footprint removal is a tremendously
supreme yet inseparable step in seismic data processing [7], [8] and is also
the focus of this article.

Great efforts have been made by the geophysics community to boost the
development of footprint noise removal approaches in recent decades [1],
[9]–[26]. Notably, a surplus of vital approaches has been advanced in the
seismic literature, and they are roughly grouped into two categories:
filtering [1], [9]–[17] and sparse representation (SR) methods [18]–[26].
Filtering methods are the most extensively used techniques and can be divided
into two classes: wavenumber/frequency-domain and time-domain filtering. The
first class, initially used by Gulunay [9], handles footprints as spatially
periodic noise, resulting in suppression attempts in the wavenumber space;
examples include classic wavenumber filtering [9], [10], adaptive wavenumber
filtering [11], spectral notch filtering [12], two-pass frequency-wavenumber
filtering [1], and wavenumber notch filtering [13]. The second class is
primarily appropriate for random noise abatement but can also be applied to
footprint removal. Typical methods include singular value decomposition
(SVD)-based [14]–[16] and structure-oriented filtering [17]. However, in case
studies, conflict arises from the fact that filtering technology attempts to
achieve significant footprint removal while ensuring that the true amplitude
is not influenced.

Compared to the wavenumber/frequency-domain approaches, SR has been suggested
to achieve more compact representations of both clear seismic images and
footprints, which allows them to be better separated in the sparsifying domain
[18]–[26]. Depending on how correctly the utilized dictionaries are formed, SR
methods can be divided into analytical sparsifying transforms and data-driven
sparse coding (SC) approaches. In the analytical sparsifying transform case,
it is assumed that noise suppression is based on the calculation of sparse
coefficients for footprints and clear images with different distribution laws
in a given analytical transform basis; representative techniques include the
wavelet transform [7], the 3-D complex wavelet transform [18], the stationary
wavelet transform (SWT) [19], [20], basis pursuit [21], and the curvelet
transform [22]. Concretely, based on the analysis of the energy distribution
laws of coefficients, the complete separation of footprints and clear seismic
images is achieved by performing threshold or filtering operations on the
sparse coefficients. Entirely unlike analytical sparsifying transforms,
considering that the atoms of the learned dictionary are found to represent
clear images and footprint patterns [23], [24], SC-based denoising methods
naturally deal with these learned atoms and ignore the role of sparse
coefficients. In particular, SC methods consider a footprint suppression
problem as a binary pattern classification [23]–[25] or filtering task [26]
for all learned atoms, which enables the separation of footprints and clear
images at the atom level. However, although the existing methods have achieved
great progress in footprint removal, their performance relies on handcrafted
priors and may be degraded if the input seismic image has unexpected
properties beyond the model assumptions.

Inspired by the overwhelming success of deep learning (DL) as a data
processing technique, DL-based methods have emerged for seismic image
denoising in recent years [27]–[51]. In contrast with conventional noise
cancellation methods, a DL-based denoising method permanently seeks to
intuitively learn an underlying mapping between a noisy seismic image and the
corresponding clean image without any handcrafted assumptions [27]–[38]. In
this case, DL-based denoising methods tend to concentrate on random noise and
aim to minimize the mean squared error (MSE) induced over the training set
using various DL models. Typical DL networks include convolutional neural
networks (CNNs) [27]–[30], 3-D denoising CNNs (3D-DnCNNs) [31], [32], deep
convolutional autoencoders (DCAEs) [33], [34], 3-D DCAEs [35], [36], and
tensor DL (TDL) methods [37], [38]. Other areas of open research include
utilizing DL for the suppression of erratic noise [39], coupled noise
[41]–[44], swell noise [45], [46], ground-roll noise [47]–[49], and desert
noise [50], [51]. However, no DL methods are currently deployed in the
geophysics community to assist with suppressing acquired footprints, since the
collection of footprints and their clean counterparts is a nontrivial task.

To ameliorate this issue, we propose a footprint removal network (FR-Net) that
is directly driven by footprint’s physical prior and related model, thus
allows a more complete unsupervised decoupling of footprints and reflected
signals compared to prevailing removal methods. The superior performance of
our suggested FR-Net can be attributed to the following characteristics. 1)
Physical prior augmented DL: Pure DL methods do not work for footprint noise;
the well-equipped MSE loss pertains to the Gaussian errors incurred during
error propagation, causing it to fail to remove footprints. Based on the
footprint’s physical nature and its modeling, our FR-Net transforms DL from
fully data-driven to physical prior augmented, enabling it to inherit the
superiority of both DL and our footprint model and thus to achieve outstanding
performance. 2) Explainable: Efforts are made in this article to produce
explainable deep neural network (xDNN) or explainable machine learning (XML)
techniques, with the goal of providing a stronger descriptive approach for the
FR-Net as well as additional prior information for users, hence enhancing the
insights obtained regarding footprint noise. 3) Ground truth-free: A somewhat
surprising feature of our FR-Net model is that, despite the use of a DL model,
the FR-Net can effectively remove footprint noise in an unsupervised manner,
which does not require the acquisition of footprints and their clean
counterparts.

Relying on the advantages listed above, it is clarified that the FR-Net is a
promising method for conducting footprint removal on noisy 3-D seismic data,
which is the primary purpose of this article. In contrast with current
footprint suppression methods, three substantial evolutions make our FR-Net
model technically feasible, and they are summarized below as our main
contributions.

  * 1)

Footprint modeling: The key to the FR-Net is to design a footprint separation
model (i.e., a unidirectional total variation (UTV) model) according to the
intrinsically directional property of footprint noise (x-line axis or y-line
axis). From the above literature survey, it can be seen that our footprint
model is quite distinct from all current removal methods.

  * 2)

Physical prior augmented DL: Inspired by [33], [34], [39], we seek to develop
a DL model that can automatically describe the ideal image behavior hidden in
contaminated seismic data. Then, strongly regularizing the DL approach using
the UTV model produces an ultimate FR-Net for thoroughly decoupling footprints
and reflected signals. To the best of our knowledge, this article is the first
work that uses DL for unsupervised footprint removal.

  * 3)

Solution: The complete separation of footprint noise and useful signals is
achieved when the FR-Net can be trained by the derived backpropagation (BP)
algorithm, which is based on the classic BP algorithm. Additionally, this BP
deduction process may supply insight into the physical prior augmented nature
of the derivative calculations during the training of the FR-Net.
Nevertheless, in the specific implementation, the derivatives of this new cost
function are processed automatically by TensorFlow [52].

  * 4)

Geophysics community: As footprints are typical coherent noise, the present
physical prior augmented idea is unprecedented with respect to the
aforementioned seismic data denoising process, particularly coherent noise
suppression. It will undoubtedly pertain to our future research direction and
dramatically enrich the theory of seismic noise abatement itself.

Extensive experimental results obtained on three public datasets demonstrate
the supremacy and efficacy of the proposed FR-Net in comparison with existing
state-of-the-art (SOTA) methods.

The rest of this paper is organized as follows. Section II gives some
background on the DCAE and total variation (TV) model used throughout this
article. The problem formulation and modeling-related issues are discussed in
Section III. Section IV completely describes the proposed work. Section V
present the experimental results and an evaluation to prove the efficacy of
the proposed FR-Net. Finally, Section VI concludes our work.

Figure 1: (a) Seismic image with footprint noise. (b) Horizontal derivative.
(c) Vertical derivative.

##  II Background

For the convenience of the reader’s understanding of this work, this part
briefly reviews two essential topics to provide the necessary preliminaries
for presenting the proposed FR-Net. Specifically, some rudimentary knowledge
concerning DCAEs is revisited in Section II-A. Then, TV is succinctly
described in Section II-B.

###  II-A DCAE

A DCAE is an end-to-end model that consists of two components: an encoder with
multiple layers that first encodes an input seismic image to lower-dimensional
latent features and a decoder that decodes the latent representation back into
a seismic image [53]. Suppose that the input $\\{\mathbf{X}_{n}\\}_{n=1}^{N}$
is known; the hidden representations for each layer can be obtained by the
following:

${\color[rgb]{0,0,0}\begin{array}[]{l}\mathbf{X}_{n,\ell}=\mathbf{h}\left(\mathbf{X}_{n,\ell-1}\right)=\sigma\left(\mathbf{W}_{\ell-1}\mathbf{X}_{n,\ell-1}+\mathbf{b}_{\ell}\right)\\\
{\mathbf{g}}\left(\mathbf{X}_{n,\ell-1}\right)=\sigma\left(\mathbf{W}_{\ell-1}^{\top}\mathbf{h}\left(\mathbf{X}_{n,\ell-1}\right)+\mathbf{b}_{\ell-1}\right),\quad\ell\in[K],\end{array}}$
(1)

where $(\cdot)^{\top}$ denotes the transpose of a given matrix and $N$
signifies the number of input sample couples. Since the decoder is designed as
the reverse of the encoder, one can feed $\mathbf{X}_{n,\ell-1}$ into the
decoder and compare the output $\mathbf{X}_{n,\ell}$ of the latter with the
original images by minimizing the subsequent reconstruction error cost.

${\color[rgb]{0,0,0}\mathcal{L}_{\mathrm{MSE}}(\mathbf{Y}_{n},\mathbf{X}_{n,\ell-1};\theta)=\arg\min_{\theta}\frac{1}{N}\sum_{n=1}^{N}\frac{1}{2}\left\|\mathbf{Y}_{n}-\mathbf{X}_{n,\ell}\right\|_{\mathrm{F}}^{2},}$
(2)

where $\theta=\\{\mathbf{W}_{\ell},\mathbf{b}_{\ell}\\}_{\ell=1}^{K}$ denotes
the parameters of the AE model.

###  II-B TV

The TV model was initially suggested by Rudin et al. [43] and was formed to
tackle grayscale image denoising issues due to its ability to retain edge
information and efficiently facilitate piecewise smoothness. Later, the
extrusion of the image denoising technique by the TV model became increasingly
critical in view of its superior performance, making it a natural choice for
extending seismic image noise lessening.

Formally, for a matrix $\mathbf{X}\in\mathbb{R}^{n_{1}\times n_{2}}$, its
derivative matrices $\nabla_{x}\mathbf{X}\in\mathbb{R}^{(n_{1}-1)\times
n_{2}}$ and $\nabla_{y}\mathbf{X}\in\mathbb{R}^{n_{1}\times(n_{2}-1)}$ are
defined as

$\left\\{\begin{aligned}
\nabla_{x}\mathbf{X}(i,j)&=\mathbf{X}(i+1,j)-\mathbf{X}(i,j)\\\
\nabla_{y}\mathbf{X}(i,j)&=\mathbf{X}(i,j+1)-\mathbf{X}(i,j),\end{aligned}\right.$
(3)

where $\mathbf{X}(i,j)$ denotes the $(i,j)$th element of $\mathbf{X}$, and
$\nabla_{x}$ and $\nabla_{y}$ denote the first-order vertical and horizontal
derivative operators, respectively. The TV of $\mathbf{X}$ is written as

$\|\mathbf{X}\|_{\mathrm{TV}}=\left\|\nabla_{x}\mathbf{X}\right\|_{{1}}+\left\|\nabla_{y}\mathbf{X}\right\|_{{1}}.$
(4)

After determining the TV norm of image denoising, we can treat it as a prior
and incorporate it to constrain the maximum a posteriori probability (MAP)
estimation, which is defined as follows:

$\hat{\mathbf{X}}=\underset{\mathbf{X}\in\mathbb{R}^{m\times
n}}{\arg\min}~{}\|\mathbf{Y}-\mathbf{X}\|_{F}^{2}+\tau\|\mathbf{X}\|_{\mathrm{TV}},$
(5)

where $\tau$ denotes the noise level of the Gaussian random noise.

Figure 2: Schematic of our FR-Net model concept.

##  III Problem Formulation and Modeling

This section proposes a physical prior augmented DL approach for footprint
removal. To elaborate, we begin by stating and formulating the footprint
removal problem in Section III-A. Then, the key to such a physical prior
augmented approach is establishing a UTV model for footprints in Section
III-B. Finally, strongly regularizing the DL approach using the UTV model
produces an ultimate FR-Net for thoroughly removing footprints in Section
III-C.

###  III-A Problem Formulation

As noted in [1], acquired footprints refer to linear spatial grid patterns
that are evident on 3-D seismic time or horizon slices. Therefore, as depicted
in Fig. 1(a), the noisy data $\mathbf{Y}$ are treated as a mixture of the
ground-truth image $\mathbf{X}$, footprint noise $\mathbf{F}$, and random
noise $\mathbf{N}$. Furthermore, given that acquired footprint noise and
random noise are additive [26], the deterioration process caused by footprints
and random noise can be stated as follows.

$\mathbf{Y}=\mathbf{X}+\mathbf{F}+\mathbf{N}.$ (6)

Formally, our objective is to infer the latent clear seismic image
$\mathbf{X}$ from the noisy image $\mathbf{Y}$ by minimizing the following
equation.

$\min_{\mathbf{X},\mathbf{F}}~{}\frac{1}{2}\|\mathbf{X}+\mathbf{F}-\mathbf{Y}\|_{F}^{2}.$
(7)

Mathematically, since the task is a typical ill-posed inverse problem, the key
is to design appropriate regularization terms to simultaneously remove both
footprints $\mathbf{F}$ and random noise $\mathbf{N}$ as follows:

$\min_{\mathbf{X},\mathbf{F}}~{}\frac{1}{2}\|\mathbf{X}+\mathbf{F}-\mathbf{Y}\|_{F}^{2}+\tau
P(\mathbf{X})+\lambda P(\mathbf{F}),$ (8)

where $P(\mathbf{X})$ and $P(\mathbf{F})$ denote the priors concerning the
clear seismic image and footprint components, respectively, and both $\tau$
and $\lambda$ denote regularization parameters. To date, the two most
typically used techniques, explicit filtering [1], [9]–[17] and SR [18]–[26],
have been advanced for footprint suppression. From these methods, it can be
seen that the key to footprint noise reduction is to build an appropriate
model for clean images and footprints, which significantly facilitates the
separation of the two components.

###  III-B Footprint Modeling

The most important aspect of footprint removal is the modeling of footprint
noise, and the main existing methods treat footprints as spatially periodic
noise that is sparse in the wavenumber or time-frequency space [9]–[13].
Similar to the existing methods, we blatantly take advantage of another
property: the spatial pattern of a footprint is regular linear stripes or grid
patterns in terms of time slices [9]–[13]. Making use of this linear feature,
a TV technique is naturally chosen for the footprint removal task because of
its edge preservation property or the spatial local smoothness [43]. More
importantly, TV can easily be imposed on the loss of the DL model in the form
of a regular term [39], [54].

However, as depicted in Fig. 1(a), it is highly visible that a footprint has a
clear directional signature. For this, we argue that classic TV is
inappropriate for footprint removal due to its noncompliance with the
directional characteristics of footprints. As an illustration, Figs. 1(b) and
(c) depict the horizontal and vertical unidirectional gradients of the
degraded seismic image, respectively. In Figs. 1(b) and (c), the footprints
have an enormous impact on the vertical gradient across the footprints, but
they have no effect whatsoever on the gradient that is detected along the
footprints. This fact compels us to limit the gradient between footprints
while retaining the gradient along footprints. To accomplish this, one can
utilize the UTV model initially developed by Bouali and Ladjal [40] that
exploits the directional characteristics of footprints. The formalization of
the UTV model is expressed as

Figure 3: Synthetic example. (a and d) Clean seismic image. (b and e)
Footprint. (c and f) Noisy seismic image.

$\mathcal{L}_{\mathrm{UTV}}(\mathbf{X})=\lambda_{1}\left\|\nabla_{x}\mathbf{X}\right\|_{1}+\lambda_{2}\left\|\nabla_{y}(\mathbf{Y}-\mathbf{X})\right\|_{1},$
(9)

where $\lambda_{1}$ and $\lambda_{2}$ denote the tradeoff parameters.
$\nabla_{x}$ and $\nabla_{y}$ represent the linear first-order difference
operator in the x-line direction and y-line direction, respectively. In this
UTV model, the first term penalizes the $\ell_{1}$-norm of the gradient along
the perpendicular direction of the footprints to remove footprints noise,
while the second term enforces the $\ell_{1}$ norm constraints on the
difference between the gradient along the direction of the footprints to
reveal or separate underlying footprints. Substituting equation (9) into (8),
one can obtain a robust stripe removal method by minimizing an objective
function that includes unidirectional variations as follows:

$\min_{\mathcal{X}}~{}\frac{1}{2}\|\mathbf{X}-\mathbf{Y}\|_{F}^{2}+\tau
P(\mathbf{X})+\lambda_{1}\left\|\nabla_{x}\mathbf{X}\right\|_{1}+\lambda_{2}\left\|\nabla_{y}(\mathbf{Y}-\mathbf{X})\right\|_{1}.$
(10)

To solve equation (10) with respect to $\mathbf{X}$, it is necessary to
construct a model for $\mathbf{X}$ that imposes additional constraints on the
highly ill-posed system equation.

###  III-C Physical Prior Augmented DL Modeling

Based on the literature survey in Section I, it is found that the current
research on modeling $\mathbf{X}$ relies mainly on handcrafted priors.
However, these prior assumptions do not always hold true for $\mathbf{X}$, and
it is necessary to use DL to perform modeling in a data-driven way. Typically,
a clear seismic image $\widehat{\mathbf{X}}$ can be modeled as a nonlinear
autoregression in this situation:

$\widehat{\mathbf{X}}=\bm{h}_{\Theta}(\mathbf{Y}).$ (11)

By conducting unsupervised training for a DNN on $\mathbf{X}$, one may
approximate the denoising model $\bm{h}$ (parameterized by $\Theta$). Using
this relation in equation (11), we may rewrite equation (2) as:

$\min_{\Theta}~{}\frac{1}{2}\|\bm{h}_{\Theta}(\mathbf{Y})-\mathbf{Y}\|^{2}.$
(12)

Afterward, by substituting equation (12) into (10), one can achieve a final
cost function for the incorrect footprint removal result, incorporating both
the physical model of $\mathbf{F}$ and the data-driven model of $\mathbf{X}$,
as follows.

$\displaystyle\min_{\Theta}{}$
$\displaystyle\frac{1}{2}\left\|\bm{h}_{\Theta}(\mathbf{Y})-\mathbf{Y}\right\|_{2}^{2}+\lambda_{1}\left\|\nabla_{y}\bm{h}_{\Theta}(\mathbf{Y})\right\|_{1}$
(13)
$\displaystyle+\lambda_{2}\left\|\nabla_{x}\left(\bm{h}_{\Theta}(\mathbf{Y})-\mathbf{Y}\right)\right\|_{1}.$

From equation (13), it is surprising to find that the optimal parameters are
changed from $\mathbf{X}$ to $\Theta$ to successfully derive a new DL model.
An additional challenge involves accounting for the exclusive DL network
structure that should be associated with the use of unsupervised footprint
removal.

Figure 4: Results obtained on the synthetic dataset in the single-trace
comparison.

##  IV Methodology

This section is separated into three distinct parts. Section IV-A gives the
proposed FR-Net model for footprint removal. Section IV-B provides the
corresponding BP algorithm for determining the new model parameters. Section
IV-C presents the implementation details for the FR-Net.

###  IV-A UTV Regularization for DCAEs

As proven by [34], [35], a DCAE can successfully denoise seismic images in an
unsupervised manner, and it outperforms other model-based learning methods. As
in our previous works [37], [39], one can use a DCAE to model clear image
patches $\widehat{\mathbf{X}}$ and train a DL network to restore a clean image
from a corrupted seismic image in a convenient form.

$\widehat{\mathbf{X}}=\bm{h}_{\Theta}(\mathbf{Y})=\bm{D}_{\Theta}\left(\bm{E}_{\Theta}(\mathbf{Y})\right),$
(14)

where $\bm{D}_{\Theta}(\cdot)$ represents the encoder with $\Theta$. According
to equation (1), $\bm{D}_{\Theta}(\cdot)$ is defined in detail as follows.

$\mathbf{X}_{n,\ell}=\bm{E}_{\Theta}\left(\mathbf{X}_{n,\ell-1}\right)=\sigma\left(\mathbf{W}_{\ell-1}\mathbf{X}_{n,\ell-1}+\overrightarrow{\mathbf{B}}_{\ell}\right),\quad\ell\in[K].$
(15)

In contrast, $\bm{E}_{\Theta}(\cdot)$ represents the decoder that upsamples
and maps the low-resolution code to a reconstruction of the original input,
which is defined as

$\bm{D}_{\Theta}\left(\mathbf{X}_{n,\ell-1}\right)=\sigma\left(\mathbf{W}_{\ell-1}^{\top}\mathbf{h}\left(\mathbf{X}_{n,\ell-1}\right)+\overrightarrow{\mathbf{B}}_{\ell-1}\right),\quad\ell\in[K].$
(16)

Afterward, substituting equation (14) into (12), we obtain

$\displaystyle\mathcal{L}=\min_{\Theta}{}$
$\displaystyle\frac{1}{2}\left\|\bm{D}_{\Theta}\left(\bm{E}_{\Theta}(\mathbf{Y})\right)-\mathcal{Y}\right\|_{2}^{2}+\lambda_{1}\left\|\nabla_{y}\bm{h}_{\Theta}(\mathbf{Y})\right\|_{1}$
(17)
$\displaystyle+\lambda_{2}\left\|\nabla_{x}\left(\bm{h}_{\Theta}(\mathbf{Y})-\mathbf{Y}\right)\right\|_{1}.$

The footprint removal problem can be modeled with a UTV-regularized AE and
solved by the BP algorithm, as in any conventional AE. However, for the
$\mathcal{L}_{\mathrm{UTV}}$ regularity problem, it is more appropriate to
solve the problem with this approach than the normal
$\mathcal{L}_{\mathrm{MSE}}$, as their cost functions are quite diverse. Thus,
it is essential to deduce a dedicated BP method to resolve the parameters
$\Theta$ of equation (17) by minimizing the new cost function (17); its
detailed derivation is described in Section IV-B.

Figure 5: Results obtained on the synthetic dataset in the x-line section
comparison. Footprint removal results are obtained using (a) the DRR method,
(b) the CBDNet method, (c) the SGRDL method, and (d) the proposed FR-Net
method. The corresponding residuals are obtained using (e) the DRR method, (f)
the CBDNet method, (g) the SGRDL method, and (h) the proposed FR-Net method.
Figure 6: Results obtained on the synthetic dataset in the constant time slice
comparison. Footprint removal results are obtained using (a) the DRR method,
(b) the CBDNet method, (c) the SGRDL method, and (d) the proposed FR-Net
method. The corresponding residuals are obtained using (e) the DRR method, (f)
the CBDNet method, (g) the SGRDL method, and (h) the proposed FR-Net method.

###  IV-B Optimization of the FR-Net

To explore the essential role of UTV regularization in footprint removal, we
provide a straightforward derivation of the BP algorithm for the FR-Net in
this section. Concretely, as an iterative optimization approach, gradient
descent (GD) is used to gradually adjust the model parameters $\Theta$ by
minimizing $\mathcal{L}\left(\mathbf{Y}_{n},\Theta\right)$ in (17). To follow
the derivation, one can obtain a clearer idea about how GD modifies
$\mathbf{W}_{\ell-1}$ and $\mathbf{B}_{\ell}$:

$\begin{array}[]{r}\mathbf{W}_{\ell-1}\leftarrow\mathbf{W}_{\ell-1}-\alpha\frac{\partial\mathcal{L}(\mathbf{W},\mathbf{B})}{\partial\mathbf{W}_{\ell-1}},\end{array}$
(18)
$\begin{array}[]{r}\mathbf{B}_{\ell}\leftarrow\mathbf{B}_{\ell}-\alpha\frac{\partial\mathcal{L}(\mathbf{W},\mathbf{B})}{\partial\mathbf{B}_{\ell}},\end{array}$
(19)

where $\alpha$ denotes the learning rate, which is a hyperparameter that
controls the step size per iteration. Typical GD algorithms proceed with
alternate updates of $\mathbf{W}_{\ell}$ and $\mathbf{B}_{\ell}$. In
particular, given the current estimate of $\mathbf{B}_{\ell}$,
$\mathbf{W}_{\ell}$ can be readily updated by computing
$\partial\mathcal{L}/\partial\mathbf{W}_{\ell-1}$ in equation (18), and
$\partial\mathcal{L}/\partial\mathbf{W}_{\ell-1}$ can be performed with
standard multiplicative rules as follows.

$\delta\mathbf{W}_{\ell-1}=\frac{\partial\mathcal{L}}{\partial\mathbf{W}_{\ell-1}}=\frac{\partial\mathcal{L}_{\mathrm{MSE}}}{\partial\widehat{\mathbf{X}}_{n,\ell}}\frac{\partial\widehat{\mathbf{X}}_{n,\ell}}{\partial\mathbf{W}_{\ell-1}}+\frac{\partial\mathcal{L}_{\mathrm{UTV}}}{\partial\widehat{\mathbf{X}}_{n,\ell}}\frac{\partial\widehat{\mathbf{X}}_{n,\ell}}{\partial\mathbf{W}_{\ell-1}}.$
(20)

First,
$\partial\mathcal{L}_{\mathrm{Welsch}}/\partial\widehat{\mathbf{X}}_{n,\ell}$
is found by taking the derivative of (2), which is:

$\frac{\partial\mathcal{L}_{\mathrm{MSE}}}{\partial\widehat{\mathbf{X}}_{n,\ell}}=\mathbf{Z}.$
(21)

To simplify the presentation, we introduce the following definition:
$\mathbf{Z}=\mathbf{Y}_{n}-\mathbf{X}_{n,\ell}$. Next,
$\partial\mathcal{L}_{\mathrm{UTV}}/\partial\widehat{\mathbf{X}}_{n,\ell}$ is
found by taking the derivative of (9):

$\frac{\partial\mathcal{L}_{\mathrm{UTV}}}{\partial\widehat{\mathbf{X}}_{n,\ell}}=\lambda_{1}\operatorname{div}\left(\frac{\nabla_{x}\widehat{\mathbf{X}}_{n,\ell}}{\|\nabla_{x}\widehat{\mathbf{X}}_{n,\ell}\|}\right)+\lambda_{2}\operatorname{div}\left(\frac{\nabla_{y}\mathbf{Z}}{\|\nabla_{y}\mathbf{Z}\|}\right),$
(22)

where $\operatorname{div}(\cdot)$ is referred to as the divergence of a given
matrix at a given point.

The second essential step is to focus on the calculation of
$\partial\mathcal{L}/\partial\mathbf{B}_{\ell}$, which is computed in a
similar way to $\delta\mathcal{W}_{\ell}$ and is defined in detail as follows:

$\delta\mathbf{B}_{\ell}=\frac{\partial\mathcal{L}}{\partial\mathbf{B}_{\ell}}=\frac{\partial\mathcal{L}_{\mathrm{MSE}}}{\partial\widehat{\mathbf{X}}_{n,\ell}}\frac{\partial\widehat{\mathbf{X}}_{n,\ell}}{\partial\mathbf{B}_{\ell}}+\frac{\partial\mathcal{L}_{\mathrm{UTV}}}{\partial\widehat{\mathbf{X}}_{n,\ell}}\frac{\partial\widehat{\mathbf{X}}_{n,\ell}}{\partial\mathbf{B}_{\ell}}.$
(23)

Now, we must smoothly compute the remaining derivatives
${\partial\widehat{\mathbf{X}}_{n,\ell}}/{\partial\mathbf{W}_{\ell-1}}$ and
${\partial\widehat{\mathbf{X}}_{n,\ell}}/{\partial\mathbf{B}_{\ell}}$; this is
specifically done by differentiating (1) with respect to $\mathbf{W}_{\ell-1}$
to compute
${\partial\widehat{\mathbf{X}}_{n,\ell}}/{\partial\mathbf{W}_{\ell-1}}$ as
follows.

$\frac{\partial\widehat{\mathbf{X}}_{n,\ell}}{\partial\mathbf{W}_{\ell-1}}=\sigma^{\prime}\left(\mathbf{W}_{\ell-1}\mathbf{X}_{n,\ell-1}+\mathbf{b}_{\ell-1}\right)\mathbf{X}_{n,\ell-1},$
(24)

where $\sigma^{\prime}$ denotes the derivative of the activation function
$\sigma$. Next, one can differentiate (1) with respect to the bias
$\mathbf{B}_{\ell}$ and compute
${\partial\widehat{\mathbf{X}}_{n,\ell}}/{\partial\mathbf{B}_{\ell-1}}$ as
follows:

$\frac{\partial\widehat{\mathbf{X}}_{n,\ell}}{\partial\mathbf{B}_{\ell-1}}=\sigma^{\prime}\left(\mathbf{W}_{\ell-1}\mathbf{X}_{n,\ell-1}+\mathbf{b}_{\ell-1}\right).$
(25)

In conclusion, the above derivation is helpful for comprehending how equations
(20) and (23) are partitioned into new UTV and MSE terms to achieve footprint
removal during the alternating updates of
$\partial\mathbf{X}_{n,\ell}/\partial\mathbf{W}_{\ell-1}$ and
$\partial\mathbf{X}_{n,\ell}/\partial\mathbf{B}_{\ell}$. Certainly, an in-
depth mathematical analysis of the above problem exceeds this article’s scope.
In addition, it should be noted that $\lambda_{1}$ and $\lambda_{2}$ are the
most relevant parameters for the FR-Net model, and their specific optimization
choices are discussed in Section V-D.

###  IV-C Implementation Details

This section highlights and treats two issues related to the details of the
FR-Net implementation: the network architecture and the construction of
unsupervised training samples. Specifically, for the $\bm{h}_{\Theta}(\cdot)$
function in equation (17), one can adopt a convolutional U-Net AE as the FR-
Net, which takes $\mathbf{Y}_{n,\ell}$ as its input to remove footprint noise.
In this case, the nineteen-layer architecture details of U-Net are shown in
Fig. 2. A contracting path and an expanding path are included in this
architecture; they are joined by symmetric skip connections, stride
convolutions, and max pooling. Specifically, the contracting path consists of
two $3\times 3$ convolutions followed by a rectified linear unit (ReLU) and a
$2\times 2$ max pooling operation with a stride of 2 for downsampling. The
expanding path performs the reverse of the contracting path process.

The second issue concerns the development of an unsupervised training sample
generation process for input 3-D seismic data. Considering that footprints
appear to cross the time slices, it is natural to flatten the 3D seismic data
into many time slice images. Then, we crop each time slice image $\mathbf{Y}$
into overlapping or nonoverlapping patches $\\{\mathbf{Y}_{n}\\}_{n=1}^{N}$ to
guarantee the adequacy of the training samples. Finally, once the trained FR-
Net completes the footprint removal procedure, the inverse process is executed
to reorganize the images, and the final result is obtained.

Figure 7: Results obtained on the Penobscot-3D dataset in the constant time
slice comparison. (a) The constant time slice and footprint removal results
are obtained using (b) the proposed FR-Net, (c) the DRR method, (d) the SGRDL
method, and (e) the CBDNet method. Figure 8: Results obtained on the
Penobscot-3D dataset in the constant time slice comparison. The removed
footprints are obtained using (a) the DRR method, (b) the proposed FR-Net
method, (c) the CBDNet method, and (d) the SGRDL method. Figure 9: Results
obtained on the Penobscot-3D dataset in the inline section comparison. (a)
Original 200th inline section. The footprint removal results are obtained
using (b) the proposed FR-Net, (d) the DRR method, (f) the CBDNet method, and
(h) the SGRDL method. The corresponding residuals are obtained using (c) the
proposed FR-Net, (e) the DRR method, (g) the CBDNet method, and (i) the SGRDL
method.  Figure 10: Results obtained on the Penobscot-3D dataset in the x-line
section comparison. (a) Original 40th x-line section. The footprint removal
results are obtained using (b) the proposed FR-Net, (d) the DRR method, (f)
the CBDNet method, and (h) the SGRDL method. The corresponding residuals are
obtained using (c) the proposed FR-Net, (e) the DRR method, (g) the CBDNet
method, and (i) the SGRDL method.  Figure 11: Results obtained on the Kerry-3D
dataset in the constant time slice comparison. (a) The constant time slice and
footprint removal results are obtained using (b) the proposed FR-Net, (c) the
DRR method, (d) the SGRDL method, and (e) the CBDNet method. Figure 12:
Results obtained on the Kerry-3D dataset in the constant time slice
comparison. The removed footprints are obtained using (a) the DRR method, (b)
the proposed FR-Net method, (c) the CBDNet method, and (d) the SGRDL method.
Figure 13: Results obtained on the Kerry-3D dataset in the inline section
comparison. (a) Original 163rd inline section. The footprint removal results
are obtained using (b) the proposed FR-Net, (d) the DRR method, (f) the CBDNet
method, and (h) the SGRDL method. The corresponding residuals are obtained
using (c) the proposed FR-Net, (e) the DRR method, (g) the CBDNet method, and
(i) the SGRDL method.  Figure 14: Results obtained on the Kerry-3D dataset in
the x-line section comparison. (a) Original 70th x-line section. The footprint
removal results are obtained using (b) the proposed FR-Net, (d) the DRR
method, (f) the CBDNet method, and (h) the SGRDL method. The corresponding
residuals are obtained using (c) the proposed FR-Net, (e) the DRR method, (g)
the CBDNet method, and (i) the SGRDL method.

##  V Experimental Results

Once our FR-Net model build is implemented, it is time to ascertain the
efficacy of such a model through a series of comparison experiments conducted
on both synthetic and real field data. First, the experimental configuration
consists of the utilized baselines, datasets, performance metrics, and
parameter settings, which are presented in Sections V-A–V-D, respectively. All
experiments are divided into two parts: one part uses synthetic data in
Section V-E, and the other part uses real field data in Section V-F. For a
related discussion concerning the experiments, see Section V-G.

###  V-A Baseline Algorithms

Since our proposed FR-Net has both classic DCAE and UTV components, we compare
our proposed method against the classic TV-regularized DCAE model [54] to
verify the benefits provided by the use of UTV. According to the seismic
literature review in Section I, we also compare our approach with the recent
footprint removal model [23] and SVD-based filtering [16], which learns how to
efficiently eliminate noise from contaminated data, as done in [23]. These
SOTA methods are detailed as follows.

  1. 1.

Convolutional blind denoising network (CBDNet)
[54]111https://github.com/GuoShi28/CBDNet: For a fair comparison, both our FR-
Net and the CBDNet model share the same network architecture but use different
TV regularity terms; i.e., the CBDNet includes conventional TV, while the FR-
Net includes UTV. This comparison helps to highlight the benefits of our UTV
penalty term for footprint removal.

  2. 2.

Statistics-guided residual dictionary learning (SGRDL)
[23]222https://github.com/chenyk1990/reproducible_research/tree/master/foot:
SGRDL generates an algorithmic framework for learning the dictionary atoms of
signal waveforms and isolating the footprint characteristics from the learned
atoms. In this instance, SGRDL utilizes an approach based on statistical
analysis for classifying dictionary atoms as footprint-affected or footprint-
free.

  3. 3.

Damped rank reduction (DRR) [16]2: By introducing a damping operator into the
conventional truncated SVD model, DRR tends to partition the vector space of
the Hankel matrix of the noisy signals into a useful signal subspace and a
noise subspace.

In the experiment, DRR and SGRDL are implemented in MATLAB, while the CBDNet
and FR-Net are realized in TensorFlow333https://www.tensorflow.org/ on the
Google Colab platform444https://colab.research.google.com/; keep in mind that
Colab assigns GPUs to users at random on the platform. In addition, for all of
the experiments, the model parameters of DRR and SGRDL are carefully adjusted
to reach performance comparable to that reported in [16] and [23].

###  V-B Testing Datasets

This section seeks to exemplify the superior performance of the proposed FR-
Net in three comparisons, one of which uses synthetic field data while the
others use actual field data. Each piece of data used in this article is
freely accessible to the public and is described in detail below.

  1. 1.

Synthetic data2: This 3-D synthetic dataset, used by [23], simply consists of
three linear seismic events, each of which has a different wavelet. The clear
data deliberately add a severe footprint to produce the final 3D observations.

  2. 2.

Penobscot-3D555https://wiki.seg.org/wiki/Penobscot_3D: Kington’s Penobscot
subset gives researchers access to a publicly available dataset with footprint
patterns that are visible on shallow time slices. Each footprint pattern
consists of a succession of stripes with different lengths that follow the
direction of a crossline.

  3. 3.

Kerry-3D666https://wiki.seg.org/wiki/Kerry-3D: The upper portion of the
Kerry-3D dataset is extensively diluted by footprint noise resulting from
horizontal bands with varied widths and lengths; this dataset is far more
polluted than Penobscot.

These datasets are sufficient as benchmarks for comparing the footprint
removal algorithms and for comprehensively illuminating the features of the
FR-Net because they vary in difficulty from easy to tough.

###  V-C Performance Metric

In the synthetic data experiments, the ground-truth data are available and can
be compared with the footprint removal results. Apart from visual comparisons,
it is quite natural to use quantitative evaluation metrics to compare the
denoising performance of the FR-Net and SOTA methods. Inspired by [23], we
also use the signal-to-noise ratio (SNR) as a quantitative evaluation metric.

$\mathrm{SNR[dB]}=10\log_{10}\frac{\|\mathbf{X}\|_{2}^{2}}{\|\mathbf{X}-\widehat{\mathbf{X}}\|_{2}^{2}},$
(26)

where $\mathbf{X}$ and $\widehat{\mathbf{X}}$ denote the clean image and the
estimated image, respectively.

###  V-D Parameter Settings

After the FR-Net architecture in Section IV-C is arbitrarily chosen, the next
task is to preset two categories of parameters: the parameters of the cost
function and the training algorithm parameters. The cost function of the FR-
Net involves two explicit tradeoff parameters that must be chosen:
$\lambda_{1}$ and $\lambda_{2}$. Here, we set $\lambda_{1}=6\times 10^{4}$ and
$\lambda_{2}=4\times 10^{5}$ for the synthetic data, $\lambda_{1}=4.5\times
10^{4}$ and $\lambda=3\times 10^{5}$ for the Penobscot data, and
$\lambda_{1}=5\times 10^{4}$ and $\lambda_{2}=3\times 10^{5}$ for the Kerry-3D
data. Once the parameters of the cost function are known, hyperparameter
optimization (i.e., with the Adam optimizer) is performed to select the
optimal training algorithm parameters, such as the batch size and learning
rate. Specifically, the final learning rate is set to $1.0\times 10^{-5}$ for
both the synthetic and Kerry-3D data and to $1.0\times 10^{-4}$ for the
Penobscot data, and the batch size is set as 5 for all datasets.

###  V-E Synthetic Data Experiments

To validate the FR-Net model, we test it on a 3-D synthetic dataset created by
the open source code in [23] with a deliberately added vertical footprint that
decays with the inline, as plotted in Fig. 3. Notably, the FR-Net does not see
clean seismic images at all during the footprint removal process because it
exists only to calculate SNR values as it did earlier. Actually, our
preprocessing stage must split the noisy seismic images into $50\times 50$
patches, where all patches can be fed to the FR-Net.

Figs. 5 and 6 show the results of two constant slice comparisons conducted on
the synthetic datasets in terms of the obtained recovery effects and
residuals. From Figs. 5 and 6, it is also confirmed that the FR-Net leaves
almost no residual footprint, while the SOTA methods cannot completely
eliminate all footprint noise. Furthermore, the footprint residuals of the
SOTA and FR-Net methods can be proven more intuitively by a single-trace
comparison, as shown in Fig. 4. It is worth noting that the FR-Net results in
Fig. 4 are nearly superb and keep the useful signal completely separate from
footprint-contaminated seismic data, thanks to the fact that all simulation
footprints attempt to strictly satisfy the UTV directional assumptions. Apart
from visual comparisons, a quantitative SNR comparison between them is
examined to characterize the residual footprint levels, as shown in Table I.
From Table I, our FR-Net model produces the smallest residual footprint (36.00
$\mathrm{dB}$) relative to those of the DRR, CBDNet, and SGRDL techniques,
which achieve values of 10.29, 20.69, and 20.10 $\mathrm{dB}$, respectively.
As observed, it is apparent that the proposed method seamlessly decouples the
valid reflected signal from noisy seismic images without any loss of
information.

###  V-F Real Field Data Experiments

To verify the superiority and effectiveness of the proposed method in
practice, the FR-Net is applied to two distinct real field datasets:
Penobscot-3D and Kerry-3D. As stated in Section V-B, these real field data
cover an extensive range of footprint noise levels, which assists in
ascertaining the adaptability of the proposed FR-Net to a variety of noise
intensities. Thus, depending on the selected footprint noise level, the entire
real field data test currently consists of the following two subexperiments.

In the first example, we demonstrate the efficacy and accuracy of the FR-Net
on the Penobscot-3D dataset. The Kerry-3D dataset possesses a large
acquisition footprint, which can involve significant engineering challenges
for the Penobscot data, as shown in Fig. 7(a). In detail, Figs. 7(b)-(e) and
Figs. 8(a)-(d) show the noise cancellation outcomes and corresponding
footprint removal effects achieved by employing the proposed FR-Net method and
the SOTA methods. From Fig. 7 and Fig. 8, it is clear that the FR-Net can
effectively remove the footprint without degrading the seismic image
resolution, while the SOTA algorithms tend to either produce oversmoothed
image edges or leave residual footprints in the resulting image. To further
compare the footprint removal effects, two lines in Fig. 7(a) are selected for
section inspections, as shown in Figs. 9 and 10. Undoubtedly, all results in
Figs. 9 and 10 replicate the initial findings and therefore corroborate the
conclusion that the FR-Net method achieves significant footprint attenuation
where no signal leakage occurs, and the edges are effectively preserved.

In the second example, we use Kerry-3-D field data with serious footprint
noise to further test our proposed FR-Net. As before, it is natural to first
compare the four methods on a time slice, as shown in Figs. 11 and 12.
Regarding visual quality, the proposed FR-Net produces the closest result to
the original seismic image except for the removal of footprint noise, making
it superior to the other SOTA methods. Similarly, Figs. 13 and 14 show the two
lines selected from Fig. 11(a) to further investigate the footprint noise
removal performance of the compared methods. From Figs. 13 and 14, it can be
observed that the FR-Net deletes fewer meaningful signals in the footprint
residual, while the SOTA approaches always cause too much damage to the useful
signals. This leads us to conclude that the proposed FR-Net approach can
greatly improve upon the footprint removal results of the SOTA methods for
Kerry-3-D data with heavy footprint noise.

###  V-G Discussion

As in [23], [24], [26], it is assumed that useful signals, footprints, and
random noise are additive. In this case, we train an FR-Net to separate useful
signals from footprints and concentrate them into the corresponding
independent components. From Figs. 5, 6, 7, and 11, it is quite evident that
the footprints are removed very cleanly, but inevitably, much random noise
remains in the useful signals, which is the principal drawback of our FR-Net.
The root cause is the problem that the pure U-Net fails to effectively
eliminate random noise in an unsupervised learning manner [54]. Notably, for
the classic noise reduction task, various methods have been introduced to
compensate for this deficiency, such as structure-oriented filtering [17].

In terms of computational complexity, our FR-Net is almost the same as the
standard U-Net because we only add a regular UTV term to its cost function and
do not change the network architecture or incur any inference overhead. As
typical DL networks, both U-Net and the FR-Net also suffer from high
computational complexity due to the recursive nature of the training
procedure, which results in a common issue for any DL model.

As with the standard U-Net, the FR-Net is also parameter-sensitive, and it
would require exhaustive experiments to obtain the optimal parameters. In
addition to these network parameters, such as the training rate and batch size
used by the standard U-Net, it is essential for the FR-Net to adjust the
tradeoff values $\lambda_{1}$ and $\lambda_{2}$ of equation (13), which
significantly affect the consequences of footprint removal. Specifically, we
must set the tradeoff parameters $\lambda_{1}$ and $\lambda_{2}$ to
appropriate values that are neither too low nor too high. If $\lambda_{1}$ and
$\lambda_{2}$ are too large, the FR-Net will remove too much footprint noise
and distort the underlying signal. If the tradeoff parameters $\lambda_{1}$
and $\lambda_{2}$ are too small, the FR-Net will be inefficient.

TABLE I: Quantitative Comparison Among the Footprint Removal Results of DRR, the CBDNet, SGRDL, and the Proposed FR-Net Method on Synthetic Datasets Method | Synthetic Dataset
---|---
DRR | $10.29~{}\mathrm{dB}$
CBDNet | ${20.69}~{}\mathrm{dB}$
SGRDL | $\text{20.10}~{}\mathrm{dB}$
Proposed FR-Net | $\textbf{36.00}~{}\mathrm{dB}$

##  VI Conclusion

This article proposes a physical prior augmented footprint removal network,
namely, the FR-Net, to suppress acquired footprints in an unsupervised way.
Specifically, to capture the intrinsically directional properties of
footprints, we design a UTV model, which is quite distinct from all current
removal methods. Then, strongly regularizing the DL method using the UTV model
produces an ultimate FR-Net, which transforms DL from an entirely data-driven
approach to a physical prior augmented technique. By jointly optimizing the
FR-Net model in an unsupervised manner, the complete separation of footprint
noise and useful signals is projected by the BP algorithm. Experimental
results obtained on synthetic and real field datasets demonstrate the
substantial gains achieved by our method over the SOTA methods in terms of
their SNR valuations and visual quality.

Future work in this area will continue to extend the proposed physical prior
augmented DL framework to other types of noise attenuation, e.g., swell noise,
ground-roll noise, and desert noise. We will also explore the combination of
UTV and other models, for example, low-rank approximation (LRA), sparse
representation (SR), and dictionary learning. In this case, we not only
successfully decoupled the footprint with the help of UTV but also nurtured
some unforeseen properties, such as the ability to compensate for FR-Net
defects that fail to suppress the random noise.

## Acknowledgments

The authors would like to thank the associate editor, and the two anonymous
reviewers, whose constructive suggestions helped to greatly improve and
clarify this article.

## References

  * [1] S. Chopra and G. Larsen, “Acquisition footprint its detection and removal,” _CSEG Recorder_ , vol. 25, no. 8, pp. 16–20, Oct. 2000.
  * [2] S. Hill, M. Shultz, and J. Brewer, “Acquisition footprint and fold-of-stack plots,” _The Leading Edge_ , vol. 18, no. 6, pp. 686–695, Jun. 1999.
  * [3] J. E. Savage and J. C. Mathewson, “Prediction of 3-d seismic footprint from existing 2-d data,” _The Leading Edge_ , vol. 20, no. 5, pp. 464–473, May 2001.
  * [4] J. Li, J. Gao, J. Mei, and X. Chen, “The footprint simulation and analysis for offshore seismic tow-cable acquisition,” in _Proc. 78th Annu. Int. Meeting, SEG Tech. Program Expanded Abstr._ , Dec. 2008, pp.120-124.
  * [5] K. J. Marfurt, R. M. Scheet, J. A. Sharp, and M. G. Harper, “Suppression of the acquisition footprint for seismic sequence attribute mapping,” _Geophys._ , vol. 63, no. 3, pp. 1024–1035, Feb. 1998.
  * [6] J. K. Cooper, “Seismic acquisition footprint: modelling and mitigation,” Master’s thesis, Graduate Studies, 2010.
  * [7] A. Alali, G. Machado, and K. J. Marfurt, “Attribute-assisted footprint suppression using a 2d continuous wavelet transform,” _Interpretation_ , vol. 6, no. 2, pp. T457–T470, May 2018.
  * [8] S. K. Sahai and K. A. Soofi, “Use of simple 2-d filters to reduce footprint noise in seismic data,” _Geohorizons_ , vol. 7, pp. 14–17, Jul. 2006.
  * [9] N. Gulunay, F. Martin, and R. Martinez, “3d data acquisition artifacts removal-spot editing in the spatial-temporal frequency domain,” in _Proc. 56th Eur. Assoc. Geoscientists Engineers (EAEG) Meeting_ , Jun. 1994, pp. cp-47-00184.
  * [10] R. Soubaras, “Attenuation of acquisition footprint for non-orthogonal 3d geometries,” in _Proc. 72nd Annu. Int. Meeting, SEG Tech. Program Expanded Abstr._ , Jan. 2002, pp. 2142–2145.
  * [11] S. Falconer and K. J. Marfurt, “Attribute-driven footprint suppression,” in _Proc. 78th Annu. Int. Meeting, SEG Tech. Program Expanded Abstr._ , Dec. 2008, pp. 2667–2671.
  * [12] N. Gülünay, “Acquisition geometry footprints removal,” in _Proc. 69th Annu. Int. Meeting, SEG Tech. Program Expanded Abstr._ , Jan. 1999, pp. 637–640.
  * [13] J. M. Drummond, A. J. Budd, and J. W. Ryan, “Adapting to noisy 3d data-attenuating the acquisition footprint,” in _Proc. 70th Annu. Int. Meeting, SEG Tech. Program Expanded Abstr._ , Jan. 2000, pp. 9–12.
  * [14] M. S. Al-Bannagi, K. Fang, P. G. Kelamis, and G. S. Douglass, “Acquisition footprint suppression via the truncated svd technique: Case studies from saudi arabia,” _The Leading Edge_ , vol. 24, no. 8, pp. 832–834, Jan. 2005.
  * [15] M. Al-Bannagi, K. Fang, P. Kelamis, and G. Douglass, “Acquisition footprint suppression via the truncated svd technique,” in _Proc. 74th Annu. Int. Meeting, SEG Tech. Program Expanded Abstr._ , Dec. 2004, pp. 1957–1960.
  * [16] Y. Chen, D. Zhang, Z. Jin, X. Chen, S. Zu, W. Huang, and S. Gan, “Simultaneous denoising and reconstruction of 5-d seismic data via damped rank-reduction method,” _Geophysical J. Int._ , vol. 206, no. 3, pp. 1695–1717, Sep. 2016.
  * [17] G. C. Fehmers and C. F. Höcker, “Fast structural interpretation with structure-oriented filtering,” _Geophysics_ , vol. 68, no. 4, pp. 1286–1293, Jul. 2003.
  * [18] Z. Yu, R. Abma, J. Etgen, and C. Sullivan, “Attenuation of noise and simultaneous source interference using wavelet denoising,” _Geophys._ , vol. 82, no. 3, pp. V179–V190, May 2017.
  * [19] X. Chen, W. Yang, Z. He, and W. Zhong, “Adaptive acquisition footprint suppression based on a 3d stationary wavelet transform: A case study from china,” _J. Appl. Geophys._ , vol. 77, pp. 1–6, Feb. 2012.
  * [20] M. Cvetkovic, S. Falconer, K. J. Marfurt, and S. Chávez-Pérez, “2d stationary-wavelet transform-based acquisition footprint suppression,” in _Proc. 77th Annu. Int. Meeting, SEG Tech. Program Expanded Abstr._ , Sep. 2007, pp. 2590–2594.
  * [21] R. Zhang, “Footprint suppression with basis pursuit denoising,” in _Proc. 79th Annu. Int. Meeting, SEG Tech. Program Expanded Abstr._ , Oct. 2009, pp. 3361–3366.
  * [22] A. Górszczyk, M. Malinowski, and G. Bellefleur, “Enhancing 3d post-stack seismic data acquired in hardrock environment using 2d curvelet transform,” _Geophysical Prospecting_ , vol. 63, no. 4, pp. 903–918, Mar. 2015.
  * [23] W. Chen, O. M. Saad, H. Wang, and Y. Chen, “Statistics-guided residual dictionary learning for footprint noise removal,” _IEEE Trans. Geosci. Remote Sens._ , vol. 60, pp. 1–11, 2022, Art no. 4502011.
  * [24] D. Liu, L. Gao, X. Wang, and W. Chen, “A dictionary learning method with atom splitting for seismic footprint suppression,” _Geophys._ , vol. 86, no. 6, pp. V509–V523, Nov. 2021.
  * [25] L. Gao, W. Xiaokai, W. Chen, and W. Chao, “A local morphological-component analysis-based acquisition footprint suppression method,” in _Proc. 86th Annu. Int. Meeting, SEG Tech. Program Expanded Abstr._ , Sep. 2016, pp. 4710-4714.
  * [26] J. L. Gómez and D. R. Velis, “Footprint removal from seismic data with residual dictionary learning,” _Geophys._ , vol. 85, no. 4, pp. V355–V365, Jul. 2020.
  * [27] S. Yu, J. Ma, and W. Wang, “Deep learning for denoisingdeep learning for denoising,” _Geophys._ , vol. 84, no. 6, pp. V333–V350, Nov. 2019.
  * [28] E. Wang and J. Nealon, “Applying machine learning to 3d seismic image denoising and enhancement,” _Interpretation_ , vol. 7, no. 3, pp. SE131–SE139, Aug. 2019.
  * [29] W. Li, H. Liu, and J. Wang, “A deep learning method for denoising based on a fast and flexible convolutional neural network,” _IEEE Trans. Geosci. Remote Sens._ , vol. 60, pp. 1–13, 2022, Art no. 5902813.
  * [30] L. Yang, W. Chen, H. Wang, and Y. Chen, “Deep learning seismic random noise attenuation via improved residual convolutional neural network,” _IEEE Trans. Geosci. Remote Sens._ , vol. 59, no. 9, pp. 7968–7981, Sept. 2021.
  * [31] D. Liu, W. Wang, W. Chen, X. Wang, Y. Zhou, and Z. Shi, “Random noise suppression in seismic data: What can deep learning do?” in _Proc. 88th Annu. Int. Meeting, SEG Tech. Program Expanded Abstr._ , Aug. 2018, pp. 2016–2020.
  * [32] D. Liu, W. Wang, X. Wang, C. Wang, J. Pei, and W. Chen, “Poststack seismic data denoising based on 3-d convolutional neural network,” _IEEE Trans. Geosci. Remote Sens._ , vol. 58, no. 3, pp. 1598–1629, Mar. 2020.
  * [33] M. Zhang, Y. Liu, and Y. Chen, “Unsupervised seismic random noise attenuation based on deep convolutional neural network,” _IEEE Access_ , vol. 7, pp. 179 810–179 822, Dec. 2019.
  * [34] O. M. Saad and Y. Chen, “Deep denoising autoencoder for seismic random noise attenuation,” _Geophys._ , vol. 85, no. 4, pp. V367–V376, Jul. 2020.
  * [35] L. Yang, S. Wang, X. Chen, O. M. Saad, W. Chen, Y. A. S. I. Oboué, and Y. Chen, “Unsupervised 3-d random noise attenuation using deep skip autoencoder,” _IEEE Trans. Geosci. Remote Sens._ , vol. 60, pp. 1–16, 2022, Art no. 5905416.
  * [36] O. M. Saad and Y. Chen, “A fully unsupervised and highly generalized deep learning approach for random noise suppression,” _Geophysical Prospecting_ , vol. 69, no. 4, pp. 709–726, Apr. 2021.
  * [37] F. Qian, Z. Liu, Y. Wang, Y. Zhou, and G. Hu, “Ground truth-free 3-d seismic random noise attenuation via deep tensor convolutional neural networks in the time-frequency domain,” _IEEE Trans. Geosci. Remote Sens._ , vol. 60, pp. 1-17, 2022, Art no. 5911317.
  * [38] F. Qian, Z. Liu, Y. Wang, S. Liao, S. Pan, and G. Hu, “Dtae: Deep tensor autoencoder for 3-d seismic data interpolation,” _IEEE Trans. Geosci. Remote Sens._ , vol. 60, pp. 1-19, 2022, Art no. 5904219.
  * [39] F. Qian, W. Guo, Z. Liu, H. Yu, G. Zhang, and G. Hu, “Unsupervised erratic seismic noise attenuation with robust deep convolutional autoencoders,” _IEEE Trans. Geosci. Remote Sens._ , vol. 60, pp. 1-16, 2022, Art no. 5913016.
  * [40] M. Bouali and S. Ladjal, “Toward optimal destriping of modis data using a unidirectional variational model,” _IEEE Trans. Geosci. Remote Sens._ , vol. 49, no. 8, pp. 2924–2935, Jun. 2011.
  * [41] Y. Zhao, Y. Li, and N. Wu, “Distributed acoustic Sens. vertical seismic profile data denoiser based on convolutional neural network,” _IEEE Trans. Geosci. Remote Sens._ , vol. 60, pp. 1-11, 2022, Art no. 5900511.
  * [42] Q. Feng and Y. Li, “Denoising deep learning network based on singular spectrum analysis das seismic data denoising with multichannel svddcnn,” _IEEE Trans. Geosci. Remote Sens._ , , vol. 60, pp. 1–11, 2022, Art no. 5902911.
  * [43] L. I. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” _Physica D: nonlinear phenomena_ , vol. 60, no. 1-4, pp. 259–268, Nov. 1992.
  * [44] Y. Zhao, Y. Li, and N. Wu, “Coupled noise reduction in distributed acoustic Sens. seismic data based on convolutional neural network,” _IEEE Geosci. Remote Sens. Lett._ , vol. 19, pp. 1-5, 2022, Art no. 8025605.
  * [45] X. Zhao, P. Lu, Y. Zhang, J. Chen, and X. Li, “Swell-noise attenuation: A deep learning approach,” _The Leading Edge_ , vol. 38, no. 12, pp. 934–942, Dec. 2019.
  * [46] B. Farmani and M. W. Pedersen, “Application of a convolutional neural network to classification of swell noise attenuation,” in _Proc. 90th Annu. Int. Meeting, SEG Tech. Program Expanded Abstr._ , Sep. 2020, pp. 2868–2872.
  * [47] H. Li, W. Yang, and X. Yong, “Deep learning for ground-roll noise attenuation,” in _Proc. 88th Annu. Int. Meeting, SEG Tech. Program Expanded Abstr._ , Aug. 2018, pp. 1981–1985.
  * [48] N. Pham and W. Li, “Physics-constrained deep learning for ground roll attenuation,” _Geophys._ , vol. 87, no. 1, pp. V15–V27, Jan. 2022.
  * [49] Y. Yuan, X. Si, and Y. Zheng, “Ground-roll attenuation using generative adversarial networksground-roll attenuation using gans,” _Geophys._ , vol. 85, no. 4, pp. WA255–WA267, Jul. 2020.
  * [50] H. Ma, H. Yao, Y. Li, and H. Wang, “Deep residual encoder–decoder networks for desert seismic noise suppression,” _IEEE Geosci. Remote Sens. Lett._ , vol. 17, no. 3, pp. 529–533, Mar. 2020.
  * [51] Y. Zhao, Y. Li, X. Dong, and B. Yang, “Low-frequency noise suppression method based on improved dncnn in desert seismic data,” _IEEE Geosci. Remote Sens. Lett._ , vol. 16, no. 5, pp. 811–815, May 2019.
  * [52] M. Abadi, P. Barham, J. Chen, Z. Chen, A. Davis, J. Dean, M. Devin, S. Ghemawat, G. Irving, M. Isard _et al._ , “TensorFlow: a system for Large-Scale machine learning,” in _Proc. 12th USENIX Symp. Operating Syst. Des. Implementation (OSDI)_ , Nov. 2016, pp. 265–283.
  * [53] G. E. Hinton and R. Zemel, “Autoencoders, minimum description length and helmholtz free energy,” in _Proc. Adv. Neural Inf. Process. Syst. (NeurIPS)_ , Dec. 1993, vol. 6, pp. 3–10.
  * [54] S. Guo, Z. Yan, K. Zhang, W. Zuo, and L. Zhang, “Toward convolutional blind denoising of real photographs,” in _Proc. IEEE/CVF Conf. Comput. Vis. Pattern Recognit. (CVPR)_ , Jun. 2019, pp. 1712–1722.
  * [55] Y. Chen, W. Huang, D. Zhang, and W. Chen, “An open-source matlab code package for improved rank-reduction 3d seismic data denoising and reconstruction,” _Comput. Geosci._ , vol. 95, pp. 59–66, Oct. 2016.


#  Lagrangian formalism and classical statistical ensemble

Sikarin Yoo-Kong
_The Institute for Fundamental Study (IF),_
_Naresuan University (NU), Phitsanulok-Nakhon Sawan,_
_99 Moo 9, Tha Pho, Mueang Phitsanulok, 65000 Phitsanulok, Thailand._
_e-mail<EMAIL_ADDRESS>
###### Abstract

The Lagrangian formulation in the classical statistical mechanics is
introduced. A key important point is that one requires to replace the standard
real time with the imaginary time through the Wick’s rotation. The area of a
constant energy-shell in the tangent bundle is preserved under the time
evolution. Consequently, a definition of the statistical ensemble can be
defined.

Keywords: Lagrangian, statistical ensemble, imaginary time.

##  1 Introduction

Classical statistical mechanics is concerned with describing the behavior of
systems consisting of a large number of particles (e.g., atoms or molecules)
where quantum effects are negligible. Traditionally, it’s based on Hamiltonian
formulation in classical mechanics, which assumes that particles have well-
defined positions and momenta [1]. In this context, the cotangent bundle
(phase space) of a system, where each point represents a microstate, is a
mathematical space that combines all possible positions and momenta of the
particles. Liouville’s theorem states that, in a Hamiltonian system as a
energy function, the hyper-volume of a constant-energy shell in contangent
bundle is conserved over time. Consequently, one can define a set of different
microstates, but share certain macroscopic properties (like energy, volume,
and particle number), called an ensemble [2].

Here comes to a question: “Is there a Lagrangian approach in the classical
statistical mechanics ?”. According to a limited knowledge of the author, it
seems to be many attempts to make a connection between Lagrangian mechanics
and classical statistical mechanics, but not in the similar fashion with those
in the Hamiltonian set up. Naively, if one shall try to construct the
classical statistical ensemble on the tangent bundle with a given Lagrangian
$L(\dot{q},q)=\dot{q}^{2}/2-V(q)$, where $(\dot{q},q)$ is a set of coordinates
in tangent bungle and $V(q)$ is a potential energy, it can be immediately
noticed that the Lagrangian cannot be used as the energy function
$L(\dot{q},q)\neq E$. Therefore, this problem prevents us to proceed further
on formulating the ensemble in tangent bundle. One interesting fact is that
one can go the quantum level by employing Hamiltonian approach or Lagrangian
approach (Feynman path integration with an imaginary time) [3]. Then, this
leads to an in-completed big picture, see figure 1, since we do not have a
proper Lagrangian formulation in the classical level of the statistical
mechanics.

Figure 1: Hamiltonian and Lagrangian approaches in the classical and quantum
statistical mechanics.

In this work, we shall provide a way to formulate the classical statistic
ensemble in the tangent bundle. To achieve this goal, the imaginary-time
(Wick’s rotation) must be applied. In section 2, the imaginary-time Lagrangian
mechanics will be discussed and the classical statistical ensemble in the
tangent bundle will be given. A simple example, harmonic oscillator with one
degree of freedom, will be used to illustrate the consistent result on
computing physical quantities both Hamiltonian and imaginary-time Lagrangian.
In section 3, the summary together with some remarks will be mentioned.

##  2 Imaginary-time Lagrangian mechanics and classical statistical ensemble

In this section, we shall provide a way to construct the classical ensemble
directly from the Lagrangian mechanics. For simplicity, we shall consider a
system with one degree of freedom and the Lagrangian is given by

$L(\dot{q},q;t)=\frac{\dot{q}^{2}}{2}-V(q)\;.$ (2.1)

Under the Wick rotation on time variable: $t\rightarrow i\tau$, the Lagrangian
becomes

$-L(\tilde{q},q)=\frac{\tilde{q}^{2}}{2}+V(q)\equiv E.$ (2.2)

where $\tilde{q}=dq/d\tau$. Interestingly, we see that a negative Lagrangian
is now a total energy of the system. With (2.2), it is not difficult to see
that

$\tilde{q}=-\frac{\partial L}{\partial\tilde{q}}\;,$ (2.3)

and the Euler-Lagrange equation is

$\frac{\partial L}{\partial q}+\frac{d}{d\tau}\frac{\partial
L}{\partial\tilde{q}}=0\;.$ (2.4)

Inserting (2.3) into (2.4), one obtains

$\tilde{\tilde{q}}=\frac{\partial L}{\partial q}\;.$ (2.5)

We next would like to show that the area of the tangent bundle preserves under
the time evolution. Suppose there is an element area: $d\tilde{q}dq|_{\tau=0}$
and a later time $\tau>0$ area is given by $d\tilde{q}dq|_{\tau>0}$. Using the
fact that $q_{\tau}=q_{0}+\tilde{q}d\tau$ and
$\tilde{q}_{\tau}=\tilde{q}_{0}+\tilde{\tilde{q}}d\tau$, we have

$\displaystyle d\tilde{q}dq|_{\tau>0}$ $\displaystyle=$
$\displaystyle\left(d\tilde{q}_{0}+d\tilde{\tilde{q}}d\tau\right)\left(dq_{0}+d\tilde{q}d\tau\right)\;.$
(2.6) $\displaystyle=$
$\displaystyle\left(d\tilde{q}_{0}+\frac{\partial\tilde{\tilde{q}}}{\partial\tilde{q}}d\tau\right)\left(dq_{0}+\frac{\partial\tilde{q}}{\partial
q}d\tau\right)$ $\displaystyle\approx$ $\displaystyle
d\tilde{q}dq|_{\tau=0}+\left(\frac{\partial\tilde{\tilde{q}}}{\partial\tilde{q}}+\frac{\partial\tilde{q}}{\partial
q}\right)d\tau\;.$

Using (2.3) and (2.5), equation (2.6) becomes

$\displaystyle d\tilde{q}dq|_{\tau>0}$ $\displaystyle=$ $\displaystyle
d\tilde{q}dq|_{\tau=0}+\left(\frac{\partial}{\partial\tilde{q}}\frac{\partial
L}{\partial q}-\frac{\partial}{\partial q}\frac{\partial
L}{\partial\tilde{q}}\right)d\tau$ $\displaystyle d\tilde{q}dq|_{\tau>0}$
$\displaystyle=$ $\displaystyle d\tilde{q}dq|_{\tau=0}\;.$ (2.7)

Therefore, the area on the tangent bundle is preserved under the time
evolution and we shall treat is feature as a Lagrangian version of the
Liouville’s theorem, see figure 2.

Figure 2: Area under the time evolution on the tangent bundle.

It now prompts to define a density function of the states on the tangent
bundle: $\rho(\tilde{q},q;\tau)$ and it is not difficult to see that the
density function satisfies the equation

$0=\frac{\partial\rho}{\partial\tau}+\frac{\partial\rho}{\partial
q}\frac{\partial
q}{\partial\tau}+\frac{\partial\rho}{\partial\tilde{q}}\frac{\partial\tilde{q}}{\partial\tau}=\frac{\partial\rho}{\partial\tau}+\\{\rho,L\\}\;,$
(2.8)

where $\\{*,L\\}$ is treated as a Lagrange bracket111We note that this
definition of the bracket is not the same with the Lagrangian bracket.

$\\{*,L\\}=\left(\frac{\partial*}{\partial\tilde{q}}\frac{\partial L}{\partial
q}-\frac{\partial*}{\partial q}\frac{\partial L}{\partial\tilde{q}}\right)\;$
(2.9)

with a property $\\{*,L\\}=-\\{L,*\\}$. Here, $*$ is a function defined on the
tangent bundle. Interesting fact is that the Lagrange bracket (2.9) provides a
set of dynamic equations as follows

$\displaystyle\tilde{q}$ $\displaystyle=$ $\displaystyle-\frac{\partial
L}{\partial\tilde{q}}\;,$ (2.10) $\displaystyle\tilde{\tilde{q}}$
$\displaystyle=$ $\displaystyle\frac{\partial L}{\partial q}\;,$ (2.11)

which are (2.3) and (2.5), respectively. One can note that a set of (2.10) and
(2.11) can be treated as a Lagrange version of the Hamilton’s eqauations.

Let $dN(\tilde{q},q;\tau)$ be a number of points in an element area
$d\tilde{q}dq$ around a point $(\tilde{q},q)$. Then probability density is
given by

$\rho(\tilde{q},q;\tau)d\tilde{q}dq=\lim_{N\rightarrow\infty}\frac{dN}{N}\;,$
(2.12)

which satisfies the normalisation condition

$\int_{\Gamma}d\Gamma\rho(\tilde{q},q;\tau)=1\;,\;\;\;d\Gamma=d\tilde{q}dq\;.$
(2.13)

Any macroscopic quantity can be computed through

$\langle\mathcal{Q}(\tilde{q},q)\rangle=\int
d\Gamma\rho(\tilde{q},q;\tau)\mathcal{Q}(\tilde{q},q)\;,$ (2.14)

which is treated as an ensemble average. Moreover, the time evolution of the
ensemble average is given by

$\displaystyle\frac{d}{d\tau}\langle\mathcal{Q}(\tilde{q},q)\rangle=\int
d\Gamma\frac{\partial\rho(\tilde{q},q;\tau)}{\partial\tau}\mathcal{Q}(\tilde{q},q)=\int
d\Gamma\mathcal{Q}(\tilde{q},q)\left(\frac{\partial\rho}{\partial\tilde{q}}\frac{\partial
L}{\partial q}-\frac{\partial\rho}{\partial q}\frac{\partial
L}{\partial\tilde{q}}\right)\;.$ (2.15)

Applying integrating by parts, we obatin

$\displaystyle\frac{d}{d\tau}\langle\mathcal{Q}(\tilde{q},q)\rangle$
$\displaystyle=$ $\displaystyle-\int
d\Gamma\rho\left[\left(\frac{\partial\mathcal{Q}}{\partial\tilde{q}}\frac{\partial
L}{\partial q}-\frac{\partial\mathcal{Q}}{\partial q}\frac{\partial
L}{\partial\tilde{q}}\right)+\mathcal{Q}\left(\frac{\partial}{\partial\tilde{q}}\frac{\partial
L}{\partial q}-\frac{\partial}{\partial q}\frac{\partial
L}{\partial\tilde{q}}\right)\right]\;.$ (2.16) $\displaystyle=$
$\displaystyle-\int
d\Gamma\rho\\{L,\mathcal{Q}\\}=\langle\\{\mathcal{Q},L\\}\rangle.$

For an equilibrium macroscopic state, the ensemble average does not explicitly
depend on time: $\frac{\partial\rho}{\partial\tau}=0$, demanding

$\\{\rho_{eq},L\\}=0\;.$ (2.17)

This means that a possible solution is for $\rho_{eq}$ to be a function of the
Lagrangian: $\rho_{eq}(L(\tilde{q},q))$. Consequently, one finds that

$\\{\rho_{eq}(L),L\\}=\rho_{eq}^{\prime}(L)\\{L,L\\}=0\;.$ (2.18)

This means that $\rho$ is constant on the energy surface $E=-L$, in tangent
bundle.

Microcanonical ensemble: For a system with fixed internal energy $U=\langle
E\rangle$, volume $V$ and a number of particles $n$, the density function is
given by

$\rho(\tilde{q},q)=\frac{1}{\Sigma}\delta(L(\tilde{q},q)+U)\;.$ (2.19)

With the $\int\rho(\tilde{q},q)d\Gamma=1$, it demands

$\Sigma(U)=\int d\Gamma\delta(L(\tilde{q},q)+U)\;.$ (2.20)

Then the ensemble defined on tangent bundle with energy less than or equal $U$
is given by

$\Omega(U)=\int_{-L(\tilde{q},q)\leq U}d\Gamma=\int_{0}^{U}dE\Sigma(E)\;.$
(2.21)

With a given ensemble in (2.21), the classical statistical entropy is given by

$S(U,V,n)=k_{B}\ln\Omega(U)\;.$ (2.22)

For two separated systems with $(U_{1},V_{1},n_{1})$ and
$(U_{2},V_{2},n_{2})$, the total ensemble is given by
$\Omega_{12}=\Omega_{1}\Omega_{2}$ resulting in

$S_{12}=S_{1}+S_{2}\;,$ (2.23)

which is known as the additive property of the entropy.

Now, we are ready to make a connection between microscopic and macroscopic
worlds. Let’s start with a definition of the temperature

$\frac{1}{T}=\left(\frac{\partial S}{\partial U}\right)\Big{|}_{V,n}\;.$
(2.24) $dS_{12}=dS_{1}+dS_{2}\;.$ (2.25)

At equilibrium, $dS_{12}=0$ and $d(U_{1}+U_{2})=0$, one obtains

$\frac{1}{T_{1}}=\left(\frac{\partial S_{1}}{\partial
U_{1}}\right)\Big{|}_{V_{1},n_{1}}=\left(\frac{\partial S_{2}}{\partial
U_{2}}\right)\Big{|}_{V_{2},n_{2}}=\frac{1}{T_{2}}\;,$ (2.26)

which gives a thermal equilibrium condition between two systems or the zeroth
law of thermodynamics.

Example: We shall now consider the imaginary-time Lagrangian for the harmonic
oscillator

$-L(\tilde{q},q;\tau)=\left(\tilde{q}^{2}+q^{2}\right)=E\;.$ (2.27)

This new form of the Lagrangian allows us to construct the classical ensemble
on tangent bundle. In the energy range: $E\rightarrow E+\delta E$, see figure
3, the number of microsates is give by

$\delta\Omega=2\pi\delta E\;,$ (2.28)

and, obviously, in the energy range $0\rightarrow E$, the number of
microstates is given by

$\Omega(E)=\frac{1}{h}\int_{0}^{E}d\Omega=\frac{E}{\hbar}\;.$ (2.29)

Applying the relation $U=\langle E\rangle=k_{B}T$ and the first law $\delta
Q=\delta U$, one finds

$\frac{\delta\Omega}{\Omega}=\frac{\delta U}{U}=\frac{\delta
Q}{k_{B}T}\rightarrow dS=\frac{\delta Q}{T}\sim\frac{\delta\Omega}{\Omega}$
(2.30)

Finally, we find that the entropy is proportional to the total number of
microstates: $S=k_{B}\ln\Omega(E)$.

Figure 3: Ensemble in the energy range $E\rightarrow E+\delta E$ on the
tangent bundle.

Canonical ensemble: A system 1 characterised by $(U_{1},V_{1},n_{1})$ is in
thermal equilibrium at temperature $T$ with a heat bath, labeled as 2,
characterised by $(E_{2},V_{2},n_{2})$ with conditions $E_{2}\gg E_{1}$ and
$n_{2}\gg n_{1}$ and $E_{12}=E_{1}+E_{2}$. The condition here is that the
energy is allowed to exchange but not of particles. The Lagrangian of the
total system is given by

$L_{12}(\tilde{q}_{1},\tilde{q}_{2},q_{1},q_{2})=L_{1}(\tilde{q}_{1},q_{1})+L(\tilde{q}_{2},q_{2})\;.$
(2.31)

Since the total system is isolated, the density function of the total system
is given by

$\rho_{12}(\tilde{q},q)=\frac{1}{\Sigma_{12}}\delta(L_{12}+E_{12})\;,$ (2.32)

where $\tilde{q}=(\tilde{q}_{1},\tilde{q}_{2})$ and $q=(q_{1},q_{2})$ and

$\Sigma_{12}=\int d\Gamma_{12}\delta(L_{12}+E_{12})\;.$ (2.33)

Then the classical ensemble is given by

$\Omega_{12}(E_{12})=\int d\Gamma_{12}\delta(L_{12}+E_{12})\;.$ (2.34)

Actually, we are interested in study the property the system 1. Then one has
to trace out 2

$\displaystyle\rho_{1}(\tilde{q}_{1},q_{1})$ $\displaystyle=$ $\displaystyle
Tr_{2}\rho_{12}(\tilde{q},q)$ (2.35) $\displaystyle=$ $\displaystyle\frac{\int
d\tilde{q}_{2}\int dq_{2}\delta(L_{1}+L_{2}+E_{12})}{\Omega_{12}(E_{12})}$
$\displaystyle=$
$\displaystyle\frac{\Omega_{2}(E_{12}+L_{1})}{\Omega_{12}(E_{12})}\;,$

where $\Omega_{2}(E_{2})=\Omega_{12}(E_{12}+L_{1})$ is a classical ensemble
for the system 2. With a condition $E_{1}\ll E_{12}$, one can expand
$\ln\Omega_{2}(E_{2})$ around the $E_{2}=E_{12}$ resulting in

$\ln\Omega_{2}(E_{2})\approx\ln\Omega_{2}(E_{12})+\frac{\partial\ln\Omega_{2}}{\partial
E_{2}}\Big{|}_{E_{2}=E_{12}}L_{1}\;.$ (2.36)

Since the system 1 and system 2 are in thermal equilibrium with temperature
$T$, then we have

$\Omega_{2}(E_{12}+L_{1})=\Omega_{2}(E_{12})e^{\frac{L_{1}}{K_{B}T}}$ (2.37)

Inserting (2.37) into (2.35), one gets

$\rho_{1}(\tilde{q}_{1},q_{1})=\frac{\Omega_{2}(E_{12})}{\Omega_{12}(E_{12})}e^{\frac{L_{1}}{K_{B}T}}=\frac{e^{\beta
L_{1}}}{\int d\Gamma_{1}e^{\beta L_{1}}}\propto e^{\beta L_{1}}\;,$ (2.38)

where $e^{\beta L_{1}}$ will be treated as a Lagrangian version of the
Boltzmann factor and $\beta=\frac{1}{k_{B}T}$.

What we have now is that, for any system with thermal equilibrium with
surrounding, the density function is given by

$\rho(\tilde{q},q)=\frac{e^{\beta L(\tilde{q},q)}}{\int d\Gamma e^{\beta
L(\tilde{q},q)}}\;,$ (2.39)

and the canonical partition function $Z$ is defined as

$Z\equiv-\frac{1}{h}\int d\Gamma e^{\beta L(\tilde{q},q)}\;.$ (2.40)

Next, we consider

$\displaystyle\frac{\partial}{\partial\beta}\ln Z$ $\displaystyle=$
$\displaystyle-\frac{\int d\Gamma(L)e^{\beta L}}{\int d\Gamma e^{\beta L}}$
(2.41) $\displaystyle=$ $\displaystyle\langle-L\rangle=U\;.$

Then we employ the relation $U=\frac{\partial}{\partial\beta}(\beta F)$, where
$F$ is the Helmholtz free energy. The final relation is

$Z=e^{\beta F}\;,\;\;\text{or}\;\;F(T,V,n)=k_{B}T\ln Z(T)\;.$ (2.42)

Example: We shall work out with the harmonic oscillator with one degree of
freedom again. The partition is given by

$Z=-\frac{1}{h}\int_{-\infty}^{+\infty}d\tilde{q}\int_{-\infty}^{+\infty}dqe^{\beta(\tilde{q}^{2}+q^{2})}=\frac{k_{B}T}{\hbar}\;,$
(2.43)

which is identical with the Hamiltonian approach.

##  3 Concluding summary

We successfully provide a systematical way on constructing the classical
statistical ensemble in the tangent bundle. A magic trick used here is the
Wick’s rotation from the real time to the imaginary time on the Lagrangian.
This transformation gives a couple of first order differential equations which
can be viewed as the Lagrangian analogue of the Hamilton’s equations. With
this new structure on tangent bundle, the area under the consideration does
not change under the time evolution and, of course, this feature could be
considered as the Lagrangian version of the Liouville’s theorem in tangent
bundle. Throughout these processes, one can naturally construct the
statistical ensemble in tangent bundle and an important quantity in this
context known as the statistical entropy or Boltzmann entropy is constructed.
With a well known example, harmonic oscillator with one degree freedom, one
can show that both approaches, Hamiltonian and imaginary-time Lagrangian, give
identical result on computing physical quantities, i.e. entropy and canonical
partition function. We hope very much that, with this preliminary work, the
missing piece in figure 1 is filled and completes a whole picture. Moreover,
with this alternative approach on the Lagrangian side, new mathematical
features might be possibly explored leading a new play ground on studying
physics.

Acknowledgements

> S. Yoo-Kong would like to express a depth of gratitude to colleagues for
> valuable discussions.

## References

  * [1] Flamm D, 1998, _History and outlook of statistical physics_ , arXiv:physics/9803005.
  * [2] Reif F, 2009, _Fundamentals of Statistical and Thermal Physics_ , Waveland Press.
  * [3] Feymann R P, 1972, _Statistical Mechanics. A set of lectures (Frontiers in physics)_ , London, United Kingdom: Perseus Books.