\section{Introduction}
\IEEEPARstart{A}{tmospheric} turbulence mitigation methods aim at removing the image degradation in long-range imaging systems due to the atmosphere with a non-uniform refractive index. From the image quality point of view, the degradation consists of spatially varying blurring effects and geometric distortions. Such degradation causes bad visibility of the object and uncertainty in object boundary and shape, which limit the performance of long-range computer vision systems, such as long-range detection, recognition, and surveillance. In addition, the long-range images may also contain motion blur and noise due to the optics and sensors. The entanglement of these degradation makes the image restoration problem more challenging than better-known problems such as deconvolution. Hence, it is necessary to develop customized turbulence mitigation modules to process the turbulence distorted images.

\begin{figure}[ht]
  \centering
   \includegraphics[width=0.98\linewidth]{images/overall2.png}
  \caption{The overall idea of our method for turbulence mitigation. Our main contribution includes a physically accurate, diverse and fast turbulence synthesis and a specifically designed neural network called turbulence mitigation transformer (TMT) for blind multi-frame turbulence removal. The TMT is trained purely on synthetic data in a supervised framework, yet shows good generalization and robustness dealing with real world turbulence degradation.}
  \label{fig:simulation}
\end{figure}

The study of atmospheric turbulence mitigation methods has made promising progress over the past decades, yet there is still a gap between launching the algorithms in practical applications. Conventional atmospheric turbulence mitigation methods \cite{Anantrasirichai2013,mao_tci,Lau2017,Milanfar2013,Hirsch2010} are often based on the "lucky phenomenon" \cite{Fried78}. They fuse many short exposure frames to form a single sharp image. Despite the satisfactory restoration of the geometric shape of the original scene, they require a large number of input frames and their runtimes are usually very slow. Recent single frame deep learning-based methods \cite{Yasarla2021ICIP,Nair2021ICIP,Lau2021} are more flexible. They demonstrate promising performance at removing the spatially varying blur and smoothing the small pixel jittering along the object boundary. However, many of these methods are customized for specific classes of images, such as human faces. This is due to the ill-posedness of correcting large geometric errors from a single measurement. In addition, since they process the data frame by frame, the temporal consistency is less enforced.

The development of learning-based algorithms was also hindered by the challenge and cost of collecting turbulence data. Even synthetic turbulence data is not easily accessible due to the complication of the forward model. Conventional simulators \cite{Hardie2017,Lau2021_sim,Dan17Sim} are either too slow or inaccurate due to oversimplifications of the model. Recently, an efficient physics-based simulator \cite{P2S} has been proposed to enable the generation of large-scale training datasets. However, the simulator in \cite{P2S} has only be applied to images with a certain low resolution grid. This greatly limit the accuracy of the turbulence model which requires a densely sampled grid.

In this paper, we make two main contributions to the blind multi-frame atmospheric turbulence mitigation problem. Firstly, we propose a few modifications to the simulator in \cite{P2S}. With the new modifications, we can now synthesize a large-scale training dataset that allows us to train deep neural networks with more diverse turbulence conditions. Secondly, we present the \textit{Turbulence Mitigation Transformer} (TMT), a multi-frame based transformer network for atmospheric turbulence mitigation. Our TMT is built upon a new attention mechanism, called the Temporal-Channel Joint Attention (TCJA). While the common self-attention modules mostly operate on spatial dimensions, TCJA operates jointly on both image channels and temporal dimensions, which allows us to effectively select and aggregate temporal information from multiple turbulence distorted frames. Our model achieves comparable performance with the conventional methods in terms of correcting the geometric distortions using only about 1/10 of the input frames compared with conventional methods. It also outperforms the SOTA networks designed for general video restoration.

\noindent Our contributions are summarized as follows:
\begin{itemize}
    \item We analyze the statistical properties of the turbulence simulation and make modifications to the Phase-to-Space (P2S) \cite{P2S} simulator. Our changes allow the simulator to generate more diverse, realistic turbulence conditions for images without loss of speed and precision.

    \item We introduce a progressive supervision architecture inspired by the physical property of turbulence. Our model can effectively remove distortion in the first stage and then restore the clear image in the second stage.

    \item We introduce temporal-channel joint attention (TCJA), a self-attention module to efficiently fuse the temporal information of atmospheric turbulence.

    \item We demonstrate good generalization of the proposed model on real world images. Our method shows superior performance to restore images on real world scenarios across a wide range of turbulence profiles.
\end{itemize}

\section{Related Work}
\noindent {\bf{Turbulence Mitigation Methods.}} Most of the conventional turbulence mitigation methods \cite{Anantrasirichai2013,mao_tci,Lau2017,Milanfar2013,Anantrasirichai2013,Hirsch2010,Hardie2017recon} are based on a common multi-frame processing pipeline, where a large number of input frames are fused to generate a high-quality image. Before the fusion, it is often necessary to first extract a geometrically stable reference frame and align the input frames to the reference using optical flow or image registration methods. Blind deconvolution is commonly used at the end to remove any residue degradation, such as the diffraction blur. The conventional methods are often very slow, with the run time to process a $256\times256$ image ranging from a few seconds to over half an hour.

Due to the difficulties to obtain real data, most of the recent deep learning methods are trained with synthetic turbulence images. \cite{P2S,Lau2021_sim} are two common choices for efficiently simulating large-scale turbulence datasets. Notably, the \textit{P2S} method \cite{P2S} is a physically accurate simulator that retains the necessary atmospheric turbulence statistics. It is adapted from \cite{chimitt2020simulating} by replacing computational expensive steps with tiny paralleled networks. The simulation technique used in this paper is an improvement over \cite{P2S} by relaxing several assumptions made in the original paper.

\noindent {\bf{General image and video restoration.}} Besides algorithms dedicated to atmospheric turbulence, image and video restoration methods designed for general restoration tasks \cite{zamir2021multi,zamir2021restormer,liang2022vrt,chan2022basicvsrpp} can also be potentially applied to the problem after re-training the model with appropriate data. These networks have often demonstrated competitive performance across the various image and video restoration tasks, such as denoising, image or video super-resolution, deraining, etc. Among these methods, the video restoration techniques are of particular interest to us, as they often adopt various strategies to utilize the temporal information from the input frames. At a high level, our work adopts a similar design to VRT \cite{liang2022vrt}, where a chunk of input frames is stored and processed parallelly. This often leads to better temporal information aggregation and higher throughput during both training and inference stages.

\noindent {\bf{Vision transformers.}} Since ViT \cite{dosovitskiy2020image} and DETR \cite{DETR} were introduced,  transformers are gradually becoming accepted as a powerful backbone for high-level vision tasks. Recent studies \cite{zamir2021restormer,Wang2022Uformer,liang2022vrt} have shown that the transformer-based models can also produce competitive results on image and video restoration tasks. Despite the promising performance, conventional transformers require memory that is quadratic in the size of input. The Swin Transformer \cite{liu2021Swin} proposed to split images into spatially shifting patches and apply self-attention within each patch to improve the efficiency. Video Swin Transformer \cite{liu2021video} extended the original 2D windows to 3D and this design was adopted by the Video Restoration Transformer (VRT) \cite{liang2022vrt} for video restoration tasks. Just like operations in \cite{ali2021xcit}, Restormer \cite{zamir2021restormer} alleviates the memory issue by constructing its self-attention purely in channel dimension, which is typically far smaller than spatial dimension in image restoration tasks. Recently, \cite{zhou2022understanding} has shown channel attention could enhance the robustness of transformers. Our design is similar to the Restormer and we introduce a new shuffle mechanism for efficient channel fusion among different frames.

\section{Our method}
\subsection{Data synthesis with fast, physics-based simulator}
\noindent Assuming that the imaging system is incoherent, the observed image follows the equation in \cite{Hirsch2010}:
\begin{equation}\label{eq:1}
\boldsymbol{O} = \boldsymbol{h}(x,y) \circledast \boldsymbol{I} + \boldsymbol{n},
\end{equation}

\noindent where $\boldsymbol{I} \in \mathbb{R}^{H\times W}$ is the ideal clean image, $\boldsymbol{O} \in \mathbb{R}^{H\times W}$ is the degraded image,  $\circledast$ denotes convolution, and $\boldsymbol{n} \in \mathbb{R}^{H\times W}$ is the random noise. With no atmospheric distortions, $\boldsymbol{h} \in \mathbb{R}^{K\times K}$ represents the spatially-invariant diffraction blur kernel. However, in the case of atmospheric turbulence, $\boldsymbol{h}(x,y)$ is a \emph{spatially varying} blur kernel. Generally, there is a different blur kernel at each pixel location $(x, y)$, which we emphasize through our notation. The atmospheric distortions are random in nature, though can be statistically described by a set of parameters which contains physical information such as the aperture diameter of the imaging system $D$, the distance between the object and image $L$, the wavelength of light $\lambda$ and the Fried parameter $r_{0}$ \cite{Fried78}.

\noindent The atmospheric distortions described by $\boldsymbol{h}$ have two principle effects on an image: blur and geometric distortion. Since the correlation between them is small \cite{chimitt2020simulating} and geometric distortion is easier to implement, separating them is preferred for efficiency. As shown in Figure \ref{fig:simulation}, we decompose the $\boldsymbol{h}$ into spatial varying blur kernels $\boldsymbol{b}$ and a grid deformation field $\boldsymbol{\mathcal{G}}$. This leads to a new expression:
\begin{equation}\label{eq:2}
\boldsymbol{O} = \boldsymbol{b}(x,y) \circledast \boldsymbol{\mathcal{G}}(\boldsymbol{I}) + \boldsymbol{n}
\end{equation}

\noindent This is a physical grounded decomposition of degradation kernel $\boldsymbol{h}$, as wildly adopted by many optical science papers. The decomposition of $\boldsymbol{h}$ into $\boldsymbol{\mathcal{G}}$ and $\boldsymbol{b}(x,y)$ gives us the flexibility to treat the blur and geometric distortion separately. However, in image processing and computer vision community, another expression is more frequently used:
\begin{equation}\label{eq:3}
\boldsymbol{O} = \boldsymbol{\mathcal{G}}(\boldsymbol{b}(x,y) \circledast \boldsymbol{I}) + \boldsymbol{n}.
\end{equation}
Because feature alignment is crucial for multi-frame image restoration, doing distortion removal before deblurring is preferred, the inverse of this operation leads to applying blur first in data synthesis. We also used the eq. \ref{eq:3} to simulate turbulence effect and we'll see in the result that using data synthesized under this scheme, neural networks actually generalize well on real world images. Even though we made compromisation to use the "wrong" order, accurate generation of $\boldsymbol{b}(x,y)$ and $\boldsymbol{\mathcal{G}}$ must be properly considered. Previous works' assumption on the blur operation $\boldsymbol{b}(x,y)$ being spatially invariant \cite{yasarla2020learning,Lau2021} is inaccurate for general scenarios. We adopt a Zernike-based approach \cite{chimitt2020simulating, P2S} to match the turbulent statistics: we generate random samples in the Zernike space, then transform them to be in the image space, providing us $\boldsymbol{b}(x,y)$ and $\boldsymbol{\mathcal{G}}$. The Zernike polynomials are a set of basis functions which can be used to describe the phase distortions of a turbulence-distorted wave \cite{noll1976zernike}. The set of all Zernike coefficients for every pixel in an image forms a random vector field. One must randomly sample according to the statistics of this random field in order to perform a correct simulation. The random field generation in \cite{P2S} requires Cholesky decomposition of a $(HWN)\times(HWN)$ correlation matrix, where $H$ and $W$ are the width and height of the random field and $N=36$ is the number of Zernike coefficients. This imposes a serious limitation on the generation of turbulent statistics, which \cite{P2S} overcomes by interpolation from a small set of points. To overcome the speed and resolution bottleneck of the P2S simulator and generate more realistic data, we made the following modifications:

\begin{figure*}[ht]
  \centering
   \includegraphics[width=0.9\linewidth]{images/simulation.png}
  \caption{Pipeline of our physics-based simulation. The blur and geometric distortion caused by turbulence can be independently modeled.}
  \label{fig:simulation}
\end{figure*}

\noindent {\bf{Accelerate the correlation matrix generation}} The correlation for a single Zernike coefficient is a function of their displacement according to the Eq. 13 in \cite{takato1995spatial}, implying the property of WSS. However, there exist non-zero cross-correlations which do not adhere to the WSS property (i.e., not a function of difference of indices). We observe the cross-correlations of the 33 high-order Zernike coefficients are small in magnitude and decay sharply in a considerably small spatial range \cite{chimitt2020simulating}, this can be modeled approximately by a exponential decay function $\exp{(-\beta x)}$, where $x$ is the distance in pixel and $\beta$ is a coefficient. To allow for FFT generation of the random fields, we neglect the inter-mode correlations in different pixels which we justify by their relatively small statistical significance. With this, we can safely generate the random field of higher order Zernike coefficients by: (1) Applying the spatial correlation via FFT to generate a set of independent Zernike coefficient fields; (2) Performing a pixel-wising mixing by Cholesky decomposition of the covariance matrix proposed by \cite{noll1976zernike}. With these changes we only require Cholesky decomposition of a $33\times33$ matrix, furthermore we can draw random field of Zernike coefficient in the same resolution as the input image. We then use P2S transform to get a physically accurate blurry image.

With respect to the spatial correlations of the second and third Zernike coefficients, the generation of the correlation kernel for each field requires the evaluation as \cite{chimitt2020simulating}:
\begin{equation}
\mathbb{E}[a_j^{\boldsymbol{x}_1}a_j^{\boldsymbol{x}_2}] = \frac{c_2}{2^{5/3}} \left(\frac{D}{r_{0}} \right)^{5/3} \left[ I_{0}(s) +\cos(2\varphi_0) I_{2}(s)\right]
\end{equation}

\noindent where $j \in \{2,3\}$, $a_j$ denotes the $j$-th Zernike coefficient, and $c_2 = 7.7554$. Additionally, $s=|\boldsymbol{x}_1 - \boldsymbol{x}_2|/D$, $\boldsymbol{x}$ is the spatial location with $\varphi_0$ as the angle between two points $\boldsymbol{x}_1, \boldsymbol{x}_2$. The terms $I_{0}$ and $I_{2}$ are integrals of Bessel functions of the first kind:
\begin{equation}
I_{i}(s) = \int_{0}^{\infty}\zeta^{-14/3} J_{i}(2s\zeta)J_{2}^{2}(\zeta) \text{d}\zeta.
\end{equation}

This evaluation is relatively time-consuming, motivating us to store a look-up table for the potential inputs. Since the range of the input value of the integral is limited, it results in a trivial loss of precision. After the aforementioned modifications, changing the turbulence profile for a $512 \times 512$ image only costs 17 ms rather than 4 minutes in \cite{P2S}. This allows us to change the turbulence parameters online during a real-time simulation.

\noindent {\bf{Employ various blur kernel sizes}} The size of blur kernel basis in P2S simulator was set to $33 \times 33$, while in real world scenarios, the size of the point spread functions varies based on the aperture of the optical system, imaging range and environment, but empirically, the blur kernels should be smaller than $33 \times 33$ to match the real-world videos. We randomly rescaled the kernel basis of P2S by various sizes (from $9\times 9$ to $33 \times 33$) but kept their original energy unchanged in the simulation process.

\noindent {\bf{Add temporal correlation}} Real world turbulence is temporally correlated, therefore, accurately simulating a video will require the distortions to be spatio-temporally correlated. We inject temporal correlations to the simulation process by correlating the initial random seed in the simulation. We use an $AR(1)$ process to generate the initial seed at the first stage. This allows for the random seed $\boldsymbol{W}_{t}$ at time $t$, which is then used to form the distortion and blur random fields, to be related to the previous realization by
\begin{equation}
\boldsymbol{W}_{t} = \alpha \boldsymbol{W}_{t-1} + \sqrt{1-\alpha^{2}}\boldsymbol{\epsilon}_{t}
\end{equation}
The term $\alpha$ is the temporal correlation ratio and $\boldsymbol{\epsilon}_{t} \sim \mathcal{N}(0, \boldsymbol{I})$.

\noindent {\bf{Simulating color images.}} As mentioned in \cite[Chapter 8.3.1]{Goodman_StatisticalOptics}, since the visible spectrum is very narrow, the difference of turbulence impact between the three color channels is negligible. Therefore, we follow the original \cite{P2S}, where the identical distortion is applied to the RGB channels.

\subsection{Turbulence mitigation transformer}
Inspired by the forward process to generate turbulence effects, we designed a two-stage network to rectify the distortion and then remove the blur pattern on the image sequences. Our network is multi-in multi-out, the number of the input and output frames is 12, we retain full temporal resolution in every level of the network. For the distortion removal stage, we use a light-weight U-Net \cite{ronneberger2015u} with depth wise 3D convolution blocks. Recent works \cite{Li_2021_ICCV, anantrasirichai2022atmospheric} have shown that one doesn't need to use an overly complicated network to remove most of the distortion, this is also consistent with our preliminary experiments. However, those works only focused on distortion rectification and left the spatial variant blurs mostly unsolved. Following the distortion removal stage, we propose the turbulence mitigation transformer to restore the clean image from blurry inputs. The overall architecture is shown in Figure \ref{fig:pipeline}, where the U-Net network on the left is the distortion removal module to measure the tilt-free image from blurry and distorted inputs, the UNet-like network on the right is the backbone of our TMT to further rectify the distortion and restore the clean image from the blurry estimations. In the second stage, we use multi-scale inputs: the original input is downsampled by 2 and 4 times to form 2 lower resolution images, they are processed by a 3D convolutional layer to extract features in multiple scales together with the original input. We concatenate those features with features from the previous encoder stages before feed them into the second and third encoder stages.

\begin{figure*}[ht]
  \centering
   \includegraphics[width=0.9\linewidth]{images/model.png}
  \caption{The scheme of our full TMT model. The input turbulence affect image sequence was first processed by a small UNet to remove the distortion, the last output is supposed to only have spatial variant blurs, it is then restored to clear image by our TMT backbone module. The number of blocks in the proposed TMT is set to be $M_{1}=M_{2}=M_{3}=M_{4}=2$ and $N_{0}=N_{1}=2, N_{2}=N_{3}=3, N_{4}=4$}
  \label{fig:pipeline}
\end{figure*}

 \paragraph{Temporal-channel joint attention module} Just like Restormer \cite{zamir2021restormer}, our TMT uses channel-wise attention to reduce memory requirements and enjoys more flexibility on spatial dimensions to avoid the quadratic complexity with respect to spatial dimensions. Moreover, it can handle different spatial shapes of input. To align and fuse temporal information, both self-attention inside one frame and mutual attention across different frames are jointly computed via TCJA. The implementation of attention inside the transformer block is shown in Figure \ref{fig:attention}.

On the right of this figure, we present two variants of feature aggregation before the computation of attention matrix. The variant $\textbf{a}$ is a plain fully connection on temporal axis following a depth-wise convolution module. But the features in different frames don't communicate in this scheme. Directly connecting all channels in the same spatial location of different frames will require large amount of parameters and memory space, to resolve this problem, we proposed the variant $\textbf{b}$: we implemented channel shuffle and partial connection inspired by the ShuffleNet \cite{zhang2018shufflenet}. We first split all channels into $C/g$ groups, where $C$ is the number of all channels and $g$ is the volume of each group, then instead of using a huge fully connection matrix on all channel and temporal tokens, we connect the temporal axis with channels in each group. By shuffling channels in every transformer block, all channels and temporal axes can communicate across several attention blocks. The detailed implementation can be described as following:

Given an input tensor $\textbf{X}\in \mathbb{R}^{C\times T\times H \times W}$, where C, T, H, W denotes the number of channels, frames, height, and width of each frame, we first filter it with a depth-wise 2D convolution (implemented by 3D convolution) to extract the spatial features and fuse them by point-wise convolution, which can also be recognized as the channel-wise fully connection operation.

For the variant \textbf{a}, we then do a temporal projection with a fully connection layer to produce $\textbf{Q}\in \mathbb{R}^{C\times H \times W \times T}$, $\textbf{K}\in \mathbb{R}^{C\times H \times W\times T}$, $\textbf{V}\in \mathbb{R}^{C\times H \times W\times T}$.

For the variant \textbf{b}, we need to split the channel dimension into $C/g$ groups. The dimension of the feature map becomes $(g, C/g, T, H, W)$. This is followed by a \textbf{channel shuffle} operation to transpose the first and second dimensions, the feature will then be in dimension $(C/g, g, T, H, W)$. We flatten the second and third dimensions and reshape to get $\hat{\textbf{X}} \in \mathbb{R}^{C/g \times H\times W \times gT}$. The following linear projection by matrix $\textbf{W} \in \mathbb{R}^{gT \times gT}$ enables feature communication among temporal information in different channels, it produces $\textbf{Q}\in \mathbb{R}^{C/g \times H \times W \times gT}$, $\textbf{K}\in \mathbb{R}^{C/g \times H \times W \times gT}$, $\textbf{V}\in \mathbb{R}^{C/g \times H \times W \times gT}$.

Before computing the attention matrix, we need to rearrange those tensors in both variants \textbf{a} and \textbf{b} into $query$ $\hat{\textbf{Q}}\in \mathbb{R}^{CT \times HW}$, $key$ $\hat{\textbf{K}}\in \mathbb{R}^{CT\times HW}$ and $value$ $\hat{\textbf{V}}\in \mathbb{R}^{CT \times HW}$. Then a standard multi-head self-attention can be implemented, we split $\hat{\textbf{Q}}$, $\hat{\textbf{K}}$ and $\hat{\textbf{V}}$ into $h$ heads, in the $i$-th head, we compute the embedding:
\begin{equation}
    \hat{\textbf{X}}_{i} = \text{softmax}\left(\frac{\hat{\textbf{Q}}_{i}\hat{\textbf{K}}_{i}^T}{\sqrt{d_{HW}}} \tau_{i}\right) \hat{\textbf{V}}_{i}
\end{equation}
Where learnable parameter $\tau_{i}$ is the temperature of the head $i$. The complexity of the channel attention is $O(\frac{HW(CT)^{2}}{h})$, which is linear to the spatial dimension and quadratic to the channel and temporal dimension.

\begin{figure*}[h]
  \centering
   \includegraphics[width=0.9\linewidth]{images/attention.png}
  \caption{The architecture of the transformer block. It consists of the channel-temporal joint attention (CTJA) module followed by a gated feed forward network. We designed two variants (a) and (b) of the CTJA block for $\text{TMT}_{a}$ and $\text{TMT}_{b}$, respectively. The only difference between those two is the channel shuffle operation.}
  \label{fig:attention}
\end{figure*}

\paragraph{Distortion removal module} Temporal alignment in video restoration is always an important component in network design. As the scene may moving, direct temporal fusion is infeasible to ensure the features to be aggregated are highly correlated. Pixel displacement can be captured by optical flow \cite{ranjan2017optical} or offsets in deformable convolution \cite{dai2017deformable, chan2022basicvsrpp}. The displacement are usually measured between 2 frames, because it is often caused by motion and itself is not a degradation. While in turbulence caused problems, the pixel displacement is a degradation and the bi-frame displacement will just be a random field, so the real position of the feature is untractable via bi-frame alignment. Hence, we introduce an explicit, light-weight distortion removal module to measure the distortion-free, blur-only images. The network is a simple UNet which uses depth wise 3D convolution as its basic operation and ReLU as its activation. It has 4 levels, the input channel dimensions of them from shallow to deep are 64, 128, 256, 512 as the spatial scale downsampled by 2 across them. In the last 3 levels of the decoder, we output 3 sequences of warping grids to progressively rectify the input image. To enlarge the perceptual field of the network, we set the kernel sizes of the 4 encoder levels as 7, 7, 5 and 3. The downsampling is realized by maxpooling and the upsample module is a depth wise 3D transposed convolution layer.

\paragraph{Loss function}
The loss function $\boldsymbol{\mathcal{L}}$ of our work is
$$\boldsymbol{\mathcal{L}} = \sum_{l=1}^{3}\gamma_{l}\boldsymbol{\mathcal{L}}_{\text{char}}(\boldsymbol{I}^{\text{blur}}, \hat{\boldsymbol{\mathcal{G}}}_{l}(\boldsymbol{I}^{\text{in}})) +  \boldsymbol{\mathcal{L}}_{\text{char}}(\boldsymbol{I}^{\text{gt}}, \hat{\boldsymbol{O}})$$

Where $\boldsymbol{\mathcal{L}}_{\text{char}}$ denotes the Charbonnier loss \cite{charbonnier}. The $\boldsymbol{I}^{\text{blur}}$ is the synthesized blur-only, distortion-free image, it is used as ground truth in the distortion removal stage. $\hat{\boldsymbol{\mathcal{G}}}_{l}(\boldsymbol{I}^{\text{in}})$ represent the $l$-th level distortion removal process which warps the input to distortion-free image. $\gamma_{l}$ is the weight of loss in different level $l$, the value is selected empirically with $\gamma_{1} = 0.6$, $\gamma_{1} = 0.3$ and $\gamma_{1} = 0.1$. $\boldsymbol{I}^{\text{gt}}$ denotes the target clean image and \hat{\boldsymbol{O}} is the restored image.

\section{Experiment}
\subsection{Datasets}
With the proposed modification on the P2S simulator, we created 2 synthetic datasets. Since most traditional methods for turbulence mitigation can only apply on static scenes, for direct comparison, we made a static scene dataset. We also created a dynamic scene dataset for restoring videos containing moving objects. Every single sequence of frames in our datasets have a unique turbulence profile, yet they can be synthesized in only one day with a NVIDIA 2080 Ti GPU, which is impossible using any other turbulence simulators. We trained and tested deep learning based restoration methods on synthetic data. To measure the generalization of our method, we used several existing real world datasets, most of them do not have the ground truth but we present several visual comparisons.
\paragraph{Synthetic datasets}
To restore turbulence effected images in as diverse scenarios as possible, the dataset should contain various of scenes. We used the place dataset \cite{zhou2017places} for synthesis for this reason. We randomly selected 9,017 images in the original dataset as input of the simulator, we generated 50 turbulence images and their associated distortion free images for each single input, resulting in 9,017 pairs of image sequences of static scenes. We split them into 7,499 pairs and 1,518 pairs, for training and testing, respectively. Figure \ref{fig:static_example} shows an example of testing on the synthetic image.

We also synthesized a video dataset for dynamic scenes. The source datasets for our dynamic scene data are the Sport Video in the Wild (SVW) dataset \cite{safdarnejad2015sports} and all ground truth videos used in TSRWGAN \cite{Jin2021NatureMI}. We mixed those videos and got 4,684 samples, the total number of frames is 1,979,564. We genereted 4,684 pairs of full turbulence and distortion free videos, then randomly split it into 3,500 videos for training and 1,184 for testing, for each testing videos, we only kept at most 120 frames.

During data synthesis, we set parameters mainly in an empirical way based on the observed real-world cases. Overall, we distribute turbulence parameters in 3 levels, they are weak, medium, and strong. Usually, within the longer distance, the larger diameter of the aperture, and the smaller Fried parameter \cite{Fried78}, the turbulence is stronger. For a stronger turbulence level, we tend to use larger blur kernels. The correlation of pixel displacement can be computed by the parameters above, then we explicitly sample the correlation coefficient of higher-order aberration. Additionally, we empirically set temporal correlations to match the visual perception of the existing real-world data for better generalizability.

The detailed parameter sampling setting is shown in table \ref{table:parameter}. Each scene has 3 lines of settings, from top to bottom the turbulence strength becomes higher. Notice for static scenes, we prefer stronger turbulence so we set a lower probability to choose the first line which corresponds to weaker turbulence.

\paragraph{Real world data}
The ultimate goal of our work is to restore the real world image sequences. To examine for the generalization of the models, we purposely test all the models using data from different source, including the static patterns from OTIS \cite{gilles2017open} dataset, dynamic scene videos from \cite{Jin2021NatureMI} and the CLEAR \cite{anantrasirichai2022atmospheric} dataset. These datasets together contain a wide-range of turbulence conditions. We present samples of the testing result on each in Figure \ref{fig:OTIS}, \ref{fig:TSRWGAN_Data} and \ref{fig:WGAN_comparison} .

\begin{table*}[ht]
  \caption{Parameter sampling setting of the data synthesis. $\beta$ denotes the correlation coefficient of higher order aberration, $D$ denotes the aperture diameter of the imaging system, $D/r_{0}$ denotes the ratio of $D$ and the Fried parameter $r_{0}$. $U(a,b)$ denotes uniform distribution in the range $(a,b)$ and [] means random choice with equal probability.}
  \label{table:parameter}
  \centering
  \resizebox{\textwidth}{!}{\setlength\doublerulesep{0.5pt}
  \begin{tabular}{cccccccc}
    \toprule[1pt]\midrule[0.3pt]
    Dataset    & Probability & Kernel size choices & $\beta$ & $D$ (m) & $D/r_{0}$ choices & Distance (m) & Temporal correlation \\
    \midrule
    \multirow{3}{*}{Static Scene} & 0.2 & [33] & [0.05, 0.1, 0.2]  & $U(0.001, 0.005)$ & [0.5, 1, 1.2, 1.5] & $U(100, 400)$ & $U(0.2, 0.6)$\\
       & 0.4 & [33]  & [0.01, 0.02, 0.05, 0.1] & $U(0.04, 0.1)$ & [1, 1.5, 2] & $U(400, 800)$ & $U(0.2, 0.6)$\\
       &  0.4 &  [33] & [0.01, 0.02, 0.05, 0.1] & $U(0.1, 0.2)$ & [1.5, 2, 3] & $U(800, 1500)$ & $U(0.2, 0.6)$\\
    \midrule
    \multirow{3}{*}{Dynamic Scene} & $1/3$ & [9, 13, 15, 21] & [0.02, 0.05, 0.1, 0.2] & $U(0.001, 0.005)$ & [0.3, 0.6, 1, 1.2] & $U(50, 400)$ & $U(0.4, 0.8)$ \\
     & $1/3$ & [11, 17, 25, 33] & [0.02, 0.05, 0.1, 0.2] & $U(0.04, 0.1)$ & [0.3, 1, 1.5] & $U(400, 800)$ & $U(0.8, 0.95)$\\
     & $1/3$ & [15, 21, 27, 33] & [0.02, 0.05, 0.1, 0.2] & $U(0.1, 0.2)$ & [1, 1.5, 2, 2.5] & $U(800, 2000)$ & $U(0.88, 0.95)$ \\
    \midrule[0.3pt]\bottomrule[1pt]
  \end{tabular}}
\end{table*}

\subsection{Training and evaluation}
We compare the proposed TMT model with the TSRWGAN \cite{Jin2021NatureMI}, a recent multi-frame turbulence mitigation method, and two state-of-the-art general video restoration methods VRT \cite{liang2022vrt} and BasicVSR++ \cite{chan2022basicvsrpp, chan2022generalization}.

All models are trained separately for static and dynamic scene datasets. We trained our model using the Adam optimizer with the Cosin Annealing scheduler \cite{loshchilov2016sgdr}. The VRT, BasicVSR++, and TSRWGAN models are all trained under the same settings as the original paper. For VRT and BasicVSR++, we use the configurations for motion deblurring as it is the closest resemble to our task. We set the input and output as 12 frames for TMT, VRT and BasicVSR++ for fair comparison. We fine-tuned the pre-trained TSRWGAN model on our dataset to maximize the performance and hence we keep the original 15 input frames setting, which shouldn't put the method at any disadvantage.

We first trained the distortion removal module for 400K iterations, the batch size was set to 4. We used Adam optimizer \cite{kingma2014adam} and with $\beta_{1}=0.9$ and $\beta_{2}=0.99$. The learning rate was initialized as $2\times10^{-4}$ then gradually decreased to $1\times 10^{-6}$ by the Cosine Annealing schedualer \cite{loshchilov2016sgdr}. We fixed the distortion removal module later and trained the TMT backbone for 800K iterations within the scheme of online training. The same Adam optimizer was applied and the learning rate follows the Cosine Annealing scheme with initial learning rate $1 \times 10^{-4}$ and final learning rate $1 \times 10^{-6}$. The model was implemented in PyTorch with a single NVIDIA A100 GPU, it requires 10 days to finish all training.
The training of VRT model mostly follows the setting for motion deblurring in the original paper \cite{liang2022vrt}. We used the scheme of online training as well, it is feasible since both TMT and VRT only used Layer normalization \cite{ba2016layer}, which is insensitive to batch size. The initial learning rate was set to $1 \times 10^{-4}$ and the model was trained from scratch for 800K iterations as well. The training was on a single NVIDIA A100 GPU and took 40 days to finish.

We fine-tuned the TSRWGAN for 400K iterations on our datasets since it was trained for the same task as ours. We used the same training scheme as the original work. The training took 5.5 days on a single NVIDIA GeForce RTX 3090 GPU to converge.
For the BasicVSR++ \cite{chan2022basicvsrpp}, we trained it from scratch for 400K iterations with batch size 4, the other training scheme is the same as the original paper.

During training, we set the input size of both TMT and VRT as $12\times 160 \times 160$, input size of TSRWGAN as $15\times 200 \times 200$ and input size of BasicVSR++ as $12\times 256 \times 256$.

For all experiments, we used the same data augmentation strategy: random Gamma correction, random saturation, random cropping, followed by a random combination of horizontal/vertial flipping and $90^{\circ}$ rotation. Since the real turbulence images are often captured using telephoto lenses with large F-number, a high ISO is required to improve the image quality, which inevitably leads to non-negligible image noise. Therefore, we also injected Gaussian random noise from $ \mathcal{N}(0, \beta \boldsymbol{I})$, where $\beta \sim U(0,0.0004]$ during training to enhance the robustness of models.

The models are evaluated on our synthetic testing set, we measured PSNR and SSIM in both RGB space and YCbCr space for quatitative comparison. The testing results are shown in Table \ref{table:comparison}. We also provide a visual comparison in Fig. \ref{fig:static_example} to show the advance of our method.

\begin{table}
\centering
\caption{Comparison on the {\bf{STATIC}} scene dataset.}
\label{tab1}
\setlength\doublerulesep{0.5pt}
\begin{tabular}{ccccc}
\toprule[1pt]\midrule[0.3pt]
Methods & PSNR & SSIM & $\text{PSNR}_{Y}$ & $\text{SSIM}_{Y}$\\
\midrule
BasicVSR++ \cite{chan2022generalization}  & 26.5055 & 0.8121 & 27.9068  & 0.8238 \\
\midrule
TSRWGAN \cite{Jin2021NatureMI}   & 25.2888 & 0.7784 & 27.3573  & 0.8017 \\
\midrule
VRT \cite{liang2022vrt}    & 27.4556 & 0.8280 & 28.9101 & 0.8451\\
\midrule
TMT [ours]   & \textbf{27.5422} & \textbf{0.8320} & \textbf{29.0011} & \textbf{0.8487} \\
\midrule[0.3pt]\bottomrule[1pt]
\end{tabular}
\end{table}

\begin{table}
\centering
\caption{Comparison on the {\bf{DYNAMIC}} scene dataset.}
\label{tab2}
\setlength\doublerulesep{0.5pt}
\begin{tabular}{ccccc}
\toprule[1pt]\midrule[0.3pt]
Methods & PSNR & SSIM & $\text{PSNR}_{Y}$ & $\text{SSIM}_{Y}$\\
\midrule
BasicVSR++ \cite{chan2022generalization} & 27.0231 & 0.8073 & 28.3055  & 0.8314 \\
\midrule
TSRWGAN \cite{Jin2021NatureMI}   & 26.3262 & 0.7957 & 28.0491  & 0.8245 \\
\midrule
VRT \cite{liang2022vrt}    & 27.6114 & 0.8300 & 29.1128  & 0.8496\\
\midrule
TMT [ours]   & \textbf{27.7419} & \textbf{0.8323} & \textbf{29.2510}  & \textbf{0.8520}\\
\midrule[0.3pt]\bottomrule[1pt]
\end{tabular}
\end{table}

\begin{figure}
    \captionsetup[subfloat]{font=scriptsize,farskip=1pt}
    \centering
  \subfloat[Input\label{fig2:Input}]{%
       \includegraphics[width=0.48\linewidth]{images/syn_static/input.jpg}}
    \hfill
  \subfloat[Ground truth\label{1b}]{%
        \includegraphics[width=0.48\linewidth]{images/syn_static/gt.png}}
    \\
  \subfloat[{TMT [ours]} \label{1c}]{%
        \includegraphics[width=0.48\linewidth]{images/syn_static/TMT.png}}
  \hfill
  \subfloat[TSRWGAN \cite{Jin2021NatureMI}]{%
        \includegraphics[width=0.48\linewidth]{images/syn_static/WGAN.png}}
    \\
  \subfloat[VRT \cite{liang2022vrt}\label{1e}]{%
        \includegraphics[width=0.48\linewidth]{images/syn_static/VRT.png}}
    \hfill
  \subfloat[BasicVSR++ \cite{chan2022basicvsrpp} \label{1f}]{%
        \includegraphics[width=0.48\linewidth]{images/syn_static/Basic.png}}
  \caption{Example of testing results on our synthetic static scene dataset. (a). Input (b).Ground truth (c). Output of {\bf{TMT [ours]}} (d).Output of TSRWGAN \cite{Jin2021NatureMI} (e). Output of VRT \cite{liang2022vrt} (f). Output of BasicVSR++ \cite{chan2022basicvsrpp}}
  \label{fig:static_example}
\end{figure}

\subsection{Ablation study}
\begin{table*}[ht]
  \caption{Quantitative comparison on the synthetic testing set for different design choices of the TMT. $\text{TMT}_{a}$ and $\text{TMT}_{b}$ denote the TMT model with variants a and b of the attention module, respectively. Here, "warp" means the distortion removal module, "SS" and "MS" are aberrations of single-scale input and multi-scale input. FLOPs is measured with input size of $12 \times 208 \times 208$}
  \label{table:ablation}
  \centering
  \setlength\doublerulesep{0.5pt}
  \begin{tabular}{ccccccccccc}
    \toprule[1pt]\midrule[0.3pt]
    Dataset & \multicolumn{4}{c}{Static Scenes} & \multicolumn{4}{c}{Dynamic Scenes}  & \multicolumn{2}{c}{Computation Consumption}   \\
    \cmidrule(r){1-1} \cmidrule(r){2-5} \cmidrule(r){6-9} \cmidrule(r){10-11}
    Methods     & PSNR & SSIM & $\text{PSNR}_{Y}$ & $\text{SSIM}_{Y}$  & PSNR & SSIM & $\text{PSNR}_{Y}$ & $\text{SSIM}_{Y}$ & $\#$ Params (M) & FLOPs (G)\\
    \midrule
    $\text{SS-TMT}_{a}$ w.o. warp  & 27.2782 & 0.8221 & 28.7316  & 0.8398 & 27.5635 & 0.8258 & 29.0644  & 0.8459  & \textbf{22.75}  & 1490 \\
    $\text{SS-TMT}_{b}$ w.o. warp  & 27.3092 & 0.8235 & 28.7626  & 0.8411 & 27.5864 & 0.8265 & 29.0886  & 0.8465 & 22.95  & 1514 \\
    $\text{SS-TMT}_{a}$  & 27.3432 & 0.8257 & 28.8013 & 0.8431 & 27.5779 & 0.8282 & 29.0863  & 0.8487 & 24.87  & 1676 \\
    $\text{SS-TMT}_{b}$  & 27.3836 & 0.8266 & 28.8396 & 0.8440 & 27.6051 & 0.8281 & 29.1110  & 0.8485 & 25.07  & 1700 \\
    $\text{MS-TMT}_{a}$ w.o. warp  & 27.4718 & 0.8291 & 28.9243  & 0.8452 & 27.6637 & 0.8301  & 29.1685  & 0.8501  & 23.70  & \textbf{1206} \\
    $\text{MS-TMT}_{b}$ w.o. warp  & 27.5003 & 0.8301 & 28.9602  & 0.8468 & 27.6841 & 0.8310 & 29.1894  & 0.8506 & 23.92 & 1304 \\
    $\text{MS-TMT}_{a}$  & 27.5215 & 0.8307 & 28.9781 & 0.8469 & -  & -  & -   & -  & 25.82  & 1392 \\
    \midrule[0.3pt]
    $\text{MS-TMT}_{b}$ [Proposed]  & \textbf{27.5422} & \textbf{0.8320} & \textbf{29.0011} & \textbf{0.8487} & \textbf{27.7419} & \textbf{0.8323} & \textbf{29.2510}  & \textbf{0.8520} & 26.04  & 1490\\
    \midrule[0.3pt]\bottomrule[1pt]
  \end{tabular}
\end{table*}

\paragraph{Influence of number of input frames}
To estimate the clear components of an image sequence under turbulence, the network must perceive a large enough area in spatial-temporal domain to estimate the property of the degradation random field, especially for samples having high temporal correlation. If the correlation factor is 1, all frames have the same degradation, the multi-frame reconstruction will be collapsed into single-frame problem, which is extremely hard to solve as revealed by \cite{vint2020analysis}. To
understand how the number of frames would affect the restoration, we trained our $\text{TMT}_{b}$ using 4, 8, 12, 20 input frames with our synthetic static scene dataset. The testing performance is shown in figure \ref{fig:nframe}. The results shows the performance of the network can be improved with more input frames. We also note that, to make fair comparison with strong baseline such as \cite{liang2022vrt}, which requires large GPU memory space so that it is hard to be deployed with more than 12 frames, we kept 12-frame setting in this paper.

\begin{figure}[ht]
    \captionsetup[subfloat]{farskip=2pt, font=scriptsize}
    \centering
  \subfloat[static scene]{%
    \includegraphics[width=0.49\linewidth]{images/frame_static.png}}
    \hfill
  \subfloat[Dynamic scene]{%
    \includegraphics[width=0.49\linewidth]{images/frame_video.png}}
   \hfill
\caption{Number of frames and performance}
\label{fig:nframe}
\end{figure}

\paragraph{Distortion removal module}
To evaluate the value of the distortion removal module, we trained our TMT backbone to directly restore images without the preprocessing of removing the distortion. The result is shown in Table \ref{table:ablation}. As it shows, the distortion removal module can boost the overall performance by $0.02 \sim 0.08$dB. It is more essential in static scenes because the statistic of the distortion is harder to estimate the dynamic scenes due to the hybrid of motion. We also have two variants of the transformer block. As mentioned in section 3, the shuffle connection between temporal and channel axes facilitates better temporal fusion in the following attention operation.

\paragraph{Influence of channel shuffle operation}
The purpose of channel shuffle operation is to facilitate the communication of channels across different frames in the TCJA module. The comparison in Table \ref{table:ablation} demonstrates the effectiveness of the channel shuffle operation, the performance benefits $0.02 \sim 0.03$ dB from the shuffle operation compared to the plain connection.

\paragraph{Influence of multi-scale input in the TMT backbone}
For general U-Net like design \cite{ronneberger2015u}, only single input is required to feed into the encoder, there're several works suggesting feeding inputs in multiple scales/resolutions could help the network to learn more efficiently \cite{8099518, ji2022multi, zamir2021multi}. In the proposed paper, we utilize multi-scale input in a very simple way: only downsampled the raw input into lower resolution. To demonstrate the advantage of this operation, we also trained our TMT with single-scale input: the input dimension of each level in the encoder is remain the same, only the input feature is changed from the concatenation of outputs produced by the last level and the 3D convolution layer to solely from the last level. The comparison in Table \ref{table:ablation} demonstrates the effectiveness of the multi-scale input, the PSNR of networks with multi-scale input is $\sim 0.15$ dB higher than that with single-scale input, while the model consumption doesn't increase.

\subsection{Generalization on real world data}
Since our simulation is a physically grounded process, samples of real-world turbulence mostly become the interpolation of samples generated by our simulation.  Other simulation methods are either much more time-consuming or inaccurate \cite{P2S}. Those 2 recent data-driven works TSRWGAN \cite{Jin2021NatureMI} and complex CNN \cite{anantrasirichai2022atmospheric} also used synthetic data to train neural networks, however, since the limitation of their simulation methods, their generalizability is constraint. However, quantitative comparisons among models on real world images are hard because most real-world image sequences don't have ground truth. Despite this problem, direct visual comparison is still possible and the superiority of our method over others in terms of generalizability in real world image restoration.

\paragraph{Visual comparison among models trained on our datasets}
In the CVPR 2022 UG2+ Challenge, a new long-range turbulence dataset is released for benchmarking turbulence mitigation algorithms \cite{UG2}. This dataset consists of 100 image sequences of text patterns, which were captured from 300 meters away during hot weather. We tested our trained models on this dataset and got promising results, furthermore, we noticed our network performs better than other networks under heavy turbulence. Several samples are given in Fig. \ref{fig:text}.

Other than the turbulence text dataset, we also tested all models on an earlier OTIS dataset \cite{gilles2017open} and TSRWGAN's real-world test set. The OTIS has 16 sequences of static patterns in different scales and turbulence levels. The TSRWGAN's real-world test set has 27 different dynamic scenes in relatively mild turbulence strength. Since the ground truth is unavailable, we only show some visual comparisons in \ref{fig:OTIS} and \ref{fig:TSRWGAN_Data}, from which one can easily conclude that models trained on our dataset generalize well on a broad range of real-world turbulence conditions and among them, our model could restore images in better visual quality than others.

\begin{figure*}[ht]
    \captionsetup[subfloat]{farskip=2pt, font=scriptsize}
    \centering
    \includegraphics[width=0.42\linewidth]{images/text/A3/pg2_input.png}%
    \hfill
    \includegraphics[width=0.14\linewidth]{images/text/A3/basic_2.png}%
    \hfill
    \includegraphics[width=0.14\linewidth]{images/text/A3/WGAN_2.jpg}%
    \hfill
    \includegraphics[width=0.14\linewidth]{images/text/A3/VRT_2.jpg}%
    \hfill
    \includegraphics[width=0.14\linewidth]{images/text/A3/TMT_2.jpg}%
    \\
    \vspace{2pt}
    \includegraphics[width=0.42\linewidth]{images/text/A3/pg24_input.png}%
    \hfill
    \includegraphics[width=0.14\linewidth]{images/text/A3/basic_24.png}%
    \hfill
    \includegraphics[width=0.14\linewidth]{images/text/A3/WGAN_24.jpg}%
    \hfill
    \includegraphics[width=0.14\linewidth]{images/text/A3/VRT_24.jpg}%
    \hfill
    \includegraphics[width=0.14\linewidth]{images/text/A3/TMT_24.jpg}%
    \\
  \subfloat[Input sequences]{%
      \includegraphics[width=0.42\linewidth]{images/text/A3/pg96_input.png}}
    \hfill
  \subfloat[BasicVSR++ \cite{chan2022basicvsrpp}]{%
    \includegraphics[width=0.14\linewidth]{images/text/A3/basic_96.png}}
   \hfill
  \subfloat[TSRWGAN \cite{Jin2021NatureMI}]{%
        \includegraphics[width=0.14\linewidth]{images/text/A3/WGAN_96.jpg}}
    \hfill
  \subfloat[VRT \cite{liang2022vrt}]{%
        \includegraphics[width=0.14\linewidth]{images/text/A3/VRT_96.jpg}}
    \hfill
  \subfloat[{TMT [ours]}]{%
        \includegraphics[width=0.14\linewidth]{images/text/A3/TMT_96.jpg}}
  \caption{\textbf{Test on sequences of real-world text data \cite{UG2}}. The 3 lines of images from top to buttom are input and output of the 2nd, 24th, 96th sequences in that dataset. Column (a). Input sequences (b).Output of BasicVSR++ \cite{chan2022basicvsrpp} (c). Output of TSRWGAN \cite{Jin2021NatureMI} (d).VRT \cite{liang2022vrt} (e).TMT [ours]}
  \label{fig:text}
\end{figure*}

\begin{figure*}[htbp]
    \captionsetup[subfloat]{farskip=2pt, font=scriptsize}
    \centering
  \subfloat[Input]{%
    \includegraphics[width=0.33\linewidth, height=0.16\linewidth]{images/OTIS/OTIS_in.png}}
    \hfill
  \subfloat[Ground truth]{%
    \includegraphics[width=0.33\linewidth, height=0.16\linewidth]{images/OTIS/OTIS_gt.png}}
   \hfill
  \subfloat[{TMT [ours]}]{%
    \includegraphics[width=0.33\linewidth, height=0.16\linewidth]{images/OTIS/OTIS_TMT.png}}
   \\
  \subfloat[TSRWGAN \cite{Jin2021NatureMI}]{%
    \includegraphics[width=0.33\linewidth, height=0.16\linewidth]{images/OTIS/OTIS_TSRWGAN.png}}
   \hfill
  \subfloat[VRT \cite{liang2022vrt}]{%
    \includegraphics[width=0.33\linewidth, height=0.16\linewidth]{images/OTIS/OTIS_VRT.png}}
   \hfill
  \subfloat[BasicVSR++ \cite{chan2022basicvsrpp}]{%
    \includegraphics[width=0.33\linewidth, height=0.16\linewidth]{images/OTIS/OTIS_basicvsr.png}}
  \caption{Visual comparison on the OTIS dataset \cite{gilles2017open}. We show the result of Pattern 15 and 16 which contain the strongest turbulence effect.}
  \label{fig:OTIS}
\end{figure*}

\begin{figure*}[!htbp]
    \centering
    \captionsetup[subfloat]{farskip=2pt, font=scriptsize}
  \subfloat[Input]{%
    \includegraphics[width=0.33\linewidth, height=0.125\linewidth]{images/car/input_new.png}}
    \hfill
  \subfloat[Original WGAN \cite{Jin2021NatureMI}]{%
    \includegraphics[width=0.33\linewidth, height=0.125\linewidth]{images/car/WGAN_original_new.png}}
   \hfill
  \subfloat[Fine-tuned WGAN]{%
    \includegraphics[width=0.33\linewidth, height=0.125\linewidth]{images/car/WGAN_ours_new.png}}
    \\
  \subfloat[BasicVSR++ \cite{chan2022basicvsrpp}]{%
    \includegraphics[width=0.33\linewidth, height=0.125\linewidth]{images/car/Basic_new.png}}
    \hfill
  \subfloat[VRT \cite{liang2022vrt}]{%
    \includegraphics[width=0.33\linewidth, height=0.125\linewidth]{images/car/VRT_new.png}}
   \hfill
  \subfloat[{TMT [ours]}]{%
    \includegraphics[width=0.33\linewidth, height=0.125\linewidth]{images/car/TMT_new.png}}
  \caption{Visual comparison on real turbulence data from TSRWGAN's dataset \cite{Jin2021NatureMI}. Model of (b) is provided by authors of \cite{Jin2021NatureMI} and model for (c) is fine-tuned on our data. From the comparison between (b) and (c), we find the robustness to relatively strong turbulence has been enhanced using our data, despite the artifacts of railing on the window still exists. Notice (f) has the least turbulence left while preserving the most details.}
  \label{fig:TSRWGAN_Data}
\end{figure*}

\paragraph{Compare with simulation methods used by TSRWGAN}
The algorithm simulation tool TSRWGAN \cite{Jin2021NatureMI} used is \cite{repasi2011computer}. This method generates physically grounded distortion and spatially invariant blur kernels, but higher-order aberrations are not simulated. \cite{Jin2021NatureMI} generated more physically synthetic data by artificial heat sources to produce turbulence effect in a short distance. However, this physical simulation tends to generate highly correlated degradation with a weak blur. Although their synthetic dataset is reasonably big and they simulated various turbulence levels, we observe that the released model does not generalize well on CLEAR's real-world dataset, OTIS \cite{gilles2017open} and the text dataset \cite{UG2} where the turbulence effect is stronger than that in the synthetic data used by TSRWGAN. We fine-tuned TSRWAGN on our synthetic data and we show the significant improvement of the TSRWGAN model in the text dataset, the OTIS dataset and CLEAR's dataset in Fig. \ref{fig:WGAN_comparison}.

\begin{figure}[ht]
    \captionsetup[subfloat]{farskip=2pt, font=scriptsize}
    \centering
  \subfloat[]{%
    \includegraphics[width=0.45\linewidth]{images/text/WGAN_input_82.png}}
    \hfill
  \subfloat[]{%
    \includegraphics[width=0.27\linewidth]{images/text/WGAN_original_82.jpg}}
   \hfill
  \subfloat[]{%
    \includegraphics[width=0.27\linewidth]{images/text/WGAN_ours82.jpg}}
    \hfill
    \\
  \subfloat[]{%
     \includegraphics[width=0.33\linewidth]{images/OTIS/OTIS_in.png}}
    \hfill
  \subfloat[]{%
    \includegraphics[width=0.33\linewidth]{images/OTIS/Original_WGAN.jpg}}
   \hfill
  \subfloat[]{%
    \includegraphics[width=0.33\linewidth]{images/OTIS/OTIS_TSRWGAN.png}}
    \hfill
    \\
  \subfloat[]{%
      \includegraphics[width=0.33\linewidth]{images/plane_supp/input.png}}
    \hfill
  \subfloat[]{%
    \includegraphics[width=0.33\linewidth]{images/plane_supp/WGAN_original.png}}
   \hfill
  \subfloat[]{%
    \includegraphics[width=0.33\linewidth]{images/plane_supp/WGAN_ours.png}}
    \hfill
  \caption{Test TSRWAGN’s generator on real-world images. (a-c). On text dataset \cite{UG2} (d-f). On OTIS dataset (g-i). On CLEAR's dataset (a). Input sequence (d, g). Single frame sample of input (b, e, h). Restoration result of TSRWGAN trained with their original data (c, f, i). Restoration result of TSRWGAN trained with our synthetic datasets. Clearly our datasets directly facilitates the generalization on real world images}
  \label{fig:WGAN_comparison}
\end{figure}

\paragraph{Compare with the simulation method used in CLEAR} Recent learning-based extension \cite{anantrasirichai2022atmospheric} of CLEAR \cite{Anantrasirichai2013} employed 9 predefined point spread function of atmospheric turbulence introduced in \cite{Hirsch2010} to simulate the degradation. Those PSFs are not physically grounded and the diversity of the turbulence conditions is very limited. Therefore, the restoration performance of \cite{anantrasirichai2022atmospheric} on some real-world data is not ideal. Since the trained model of CLEAR for video restoration is not released, we were not able to train it on our dataset to demonstrate the potential improvement. Despite this, based on our discussions and comparisons in the previous sections, we can still conclude that using our proposed simulation method will make more positive impacts to the reconstruction.

\begin{figure*}[ht]
    \captionsetup[subfloat]{farskip=1pt, font=scriptsize}
    \centering
  \subfloat[Input: $\#12$ frame]{%
    \includegraphics[width=0.33\linewidth]{images/hillhouse/hillhouse_distorted.png}}
    \hfill
  \subfloat[Conventional CLEAR \cite{Anantrasirichai2013}]{%
    \includegraphics[width=0.33\linewidth]{images/hillhouse/hillhouse_restored.png}}
   \hfill
  \subfloat[{TMT [ours]}]{%
    \includegraphics[width=0.33\linewidth]{images/hillhouse/TMT.jpg}}
  \caption{\textbf{Comparison with a representative conventional method CLEAR \cite{Anantrasirichai2013} on the real-world image sequence}. Our model is trained on the synthetic static scene dataset. The presented input frame is the 12th in the whole sequence. {\bf{Zoom in for better view}}.}
  \label{fig:Conventional_clear}
\end{figure*}

\subsection{Pixel displacement rectification by distortion removal module}
Our distortion removal module is robust to different levels of blurs caused by turbulence, it can be used alone without the TMT backbone in some cases. For examples, there are scenarios where the field of view is large and the optical path is short, in this case, the blurry effect will be weak and the turbulence is mainly dominated by distortion. As shown in figure \ref{fig:tilt_street}, if we directly apply the distortion removal model only, we can already give a visually satisfying reconstruction.

\begin{figure*}[ht]
    \captionsetup[subfloat]{farskip=2pt, font=scriptsize}
    \centering
    \includegraphics[width=0.33\linewidth]{images/A2/input_init.png}%
    \hfill
    \includegraphics[width=0.33\linewidth]{images/A2/CLEAR_init.png}%
    \hfill
    \includegraphics[width=0.33\linewidth]{images/A2/tilt_init.png}%
    \\
  \subfloat[]{%
      \includegraphics[width=0.33\linewidth]{images/A2/input_processed.png}}
    \hfill
  \subfloat[]{%
    \includegraphics[width=0.33\linewidth]{images/A2/CLEAR_processed.png}}
   \hfill
  \subfloat[]{%
    \includegraphics[width=0.33\linewidth]{images/A2/tilt_processed.png}}
  \caption{Apply our distortion removal module on CLEAR's real-world data. Images on the top are full resolution images and below are the same patches in 5 frames starting from them. (a). Input (b). Processed by CLEAR {c}. Processed by our distortion removal module only}
  \label{fig:tilt_street}
\end{figure*}

\begin{figure*}[ht]
\captionsetup[subfloat]{farskip=2pt, font=scriptsize}
    \centering
    \includegraphics[width=0.33\linewidth]{images/A2_airport/input_init.png}%
    \hfill
    \includegraphics[width=0.33\linewidth]{images/A2_airport/CLEAR_init.png}%
    \hfill
    \includegraphics[width=0.33\linewidth]{images/A2_airport/tilt_init.png}%
    \\
  \subfloat[]{%
      \includegraphics[width=0.33\linewidth]{images/A2_airport/input_processed.png}}
    \hfill
  \subfloat[]{%
    \includegraphics[width=0.33\linewidth]{images/A2_airport/CLEAR_processed.png}}
   \hfill
  \subfloat[]{%
    \includegraphics[width=0.33\linewidth]{images/A2_airport/tilt_processed.png}}
  \caption{Apply our distortion removal module on CLEAR's real-world data. Images on the top are full resolution images and below are the same patches in 5 frames starting from them. (a). Input (b). Processed by CLEAR {c}. Processed by our distortion removal module only}
  \label{fig:tilt_airport}
\end{figure*}

\subsection{Comparison with conventional methods}
We make a comparison with two state-of-the-art conventional turbulence mitigation methods \cite{mao_tci,Anantrasirichai2013}. Due to the slow processing speed of the conventional ones, we only compare the methods on a random subset of 100 scenes from our synthetic static dataset. Conventional methods could still work because our synthetic data has similar properties compared to real-world turbulence degraded images. For example, the lucky effect can be synthesized with this simulator \cite{P2S}. It is worth mentioning that we used 50 frames as input for conventional methods and 12/15 for learning-based methods.

\begin{table}[h]
  \caption{Quantitative comparison with conventional turbulence mitigation methods on the 100 scenes dataset. $\text{PSNR}_{Y}$ and $\text{SSIM}_{Y}$ denotes the luminance channel in YCbCr color space.}
  \vspace{1ex}
  \label{table:comparison}
  \centering
  \resizebox{0.45\textwidth}{!}{\setlength\doublerulesep{0.5pt}
  \begin{tabular}{cccccc}
    \toprule[1pt]\midrule[0.3pt]
    \multicolumn{2}{c}{Methods}     & PSNR & SSIM & $\text{PSNR}_{Y}$ & $\text{SSIM}_{Y}$\\
    \midrule
    \multirow{2}{*}{Conventional} & Mao et al. \cite{mao_tci}    & 27.4154 & 0.8174 & 28.9619 & 0.8497\\
    & CLEAR \cite{Anantrasirichai2013}    & 26.9969 & 0.8295 & 28.4906 & 0.8525\\
    \midrule
    \multirow{4}{*}{Deep learning} &
    BasicVSR++ \cite{chan2022generalization}  & 28.0811 & 0.8410 & 29.5210  & 0.8588 \\
    & TSRWGAN \cite{Jin2021NatureMI}   & 25.8139 & 0.7977 & 27.9967  & 0.8451\\
    & VRT \cite{liang2022vrt}    & 28.1112 & 0.8456 & 29.5735 & 0.8632\\
    & TMT [ours]   & \textbf{28.2170} & \textbf{0.8477} & \textbf{29.6821} & \textbf{0.8642}\\
    \midrule[0.3pt]\bottomrule[1pt]
  \end{tabular}}
\end{table}

Additionally, we compare our proposed method with conventional CLEAR's method \cite{Anantrasirichai2013} on real-world images. A qualitative comparison is presented in figure \ref{fig:Conventional_clear}. CLEAR \cite{Anantrasirichai2013} used a subsample mechanism to select images in good quality among all 75 frames of the input sequence, then used a complex wavelet-based algorithm to fuse selected frames. After those operations, the output is still subjected to explicit denoising, while our method just blindly used the first 12 frames as input and produced more natural and high-contrast reconstruction without further lose of details.

\subsection{Computation consumption}
We present the inference-time computing budget of all models we trained in this section. The measurement is based on a single NVIDIA 2080 Ti GPU to match a very common deployment environment, we set the spatial size of the input image sequence to be $540 \times 960$. When the full size can't fit into the GPU, we split the input into largest possible patches. The splitting should be overlapped to reduce artifacts, for example, the maximum patches of VRT is $160 \times 160$, so we need $4\times 7$ patches to reconstruct a full resolution image sequence. From the table \ref{table:comparison}, although VRT has the closest performance in terms of PSNR and SSIM to ours, TMT is much more efficient. Because of the efficiency, we can use larger patches as input during testing and get more precise estimation on turbulence conditions.
\begin{table}[h]
  \caption{Comparison of computational consumption in inference. We show 2 settings of the distortion removal model, the first line is when it is a part of TMT, the second line shows the consumption when it's split into its maximum allowable patches on our testing environment. Note the TSRWGAN is multi-in, single-out, we need to run it 15 times to get reconstruction of 15 frames.}
  \label{table:comparison}
  \centering
  \resizebox{0.98\linewidth}{!}{\setlength\doublerulesep{0.5pt}
  \begin{tabular}{cccc}
    \toprule[1pt]\midrule[0.3pt]
     \multirow{2}{*}{Methods}     & $\#$ parameters   & Patch size & FLOPs  \\
        & M  & $\text{T}\times \text{H} \times \text{W}$ &  $\text{G} \times \#$ times
    \\
    \midrule
    BasicVSR++ \cite{chan2022generalization}  & 9.76  & $ 12\times 256 \times 256$  & $466.0 \times 12$  \\
    TSRWGAN \cite{Jin2021NatureMI}   & 46.28  & $ 15\times 540 \times 960$  & $2243 \times 15$ \\
    VRT \cite{liang2022vrt}  &  18.32 & $12\times 160 \times 160$ & $3324 \times 28$ \\
    \midrule
    $\text{TMT}_{a}$ & 24.87 & $12\times 208 \times 208$ & $1676 \times 15$ \\
    $\text{TMT}_{b}$ & 25.07  & $12\times 208 \times 208$ & $1700 \times 15$  \\
    \multirow{2}{*}{Distortion removal module} & 2.120 & $12\times 208 \times 208$ &  $186.0 \times 15$ \\
    & 2.120 & $12\times 320 \times 512$ &  $704.5 \times 4$ \\
    \midrule[0.3pt]\bottomrule[1pt]
  \end{tabular}}
\end{table}

\section{Conclusion}
In this paper, we presented a complete and generalizable solution for a multi-frame blind atmospheric turbulence mitigation problem. We refined a physics-based simulation in \cite{P2S} by approximating the spatially varying field using a wide sense stationary field. We further proposed an efficient transformer-based network to restore image sequences degraded by turbulence and our method can adapt to various scenes and turbulence conditions, our model surpassed the state of art models designed for general video restoration tasks. Compared to conventional methods on turbulence mitigation, the proposed method has superior performance in visual quality, speed, and robustness.

\section*{Acknowledgments and Disclosure of Funding}
The research is based upon work supported in part by the Intelligence Advanced Research Projects Activity (IARPA) under Contract No. 2022‐21102100004, and in part by the National Science Foundation under the grants CCSS-2030570 and IIS-2133032. The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies, either expressed or implied, of IARPA, or the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for governmental purposes notwithstanding any copyright annotation therein.

\bibliographystyle{IEEEtran}

\section{Introduction}
\IEEEPARstart{A}{tmospheric} turbulence mitigation methods aim at removing the image degradation in long-range imaging systems due to the atmosphere with a non-uniform refractive index. From the image quality point of view, the degradation consists of spatially varying blurring effects and geometric distortions. Such degradation causes bad visibility of the object and uncertainty in its boundary and shape, limiting the performance of long-range computer vision applications such as detection, recognition, and surveillance. In addition, the long-range images may also contain motion blur and noise due to the optics and sensors. The entanglement of these degradation makes the image restoration problem more challenging than better-known problems such as deconvolution. Hence, it is necessary to develop customized turbulence mitigation modules to process the turbulence distorted images.

\begin{figure}[ht]
  \centering   \includegraphics[width=0.99\linewidth]{images/intro_TMT.png}  \caption{The overall idea of our method for turbulence mitigation. Our main contribution includes a physically accurate, diverse and fast turbulence synthesis and a specifically designed neural network called turbulence mitigation transformer (TMT) for blind multi-frame turbulence removal. The TMT is trained purely on synthetic data in a supervised framework, yet shows good generalization and robustness dealing with real world turbulence degradation in both dynamic and static scenes.}
  \label{fig:intro}
\end{figure}

The study of atmospheric turbulence mitigation methods has made promising progress over the past decades, yet there remains a gap in deploying reconstruction algorithms for practical applications. Conventional atmospheric turbulence mitigation methods \cite{Anantrasirichai2013,mao_tci,Lau2017,Milanfar2013,Hirsch2010} are often based on the ``lucky phenomenon'' \cite{Fried78} where a patch in an image may be minimally degraded by turbulence due to its stochastic nature. These high-quality patches are combined to form a single sharp image. Despite the satisfactory restoration of the geometric shape of the original scene, they require a large number of input frames and their runtimes are usually very slow. Recent single frame deep learning-based methods \cite{Yasarla2021ICIP,Nair2021ICIP,Lau2021} have demonstrated promising performance at removing the spatially varying blur and smoothing the small pixel jittering along the object boundary. However, many of these methods are customized for specific classes of images, such as human faces. This is due to the ill-posed nature of correcting large geometric errors from a single measurement. In addition, frames are processed independently, disregarding any temporal consistency.

The development of learning-based algorithms is also hindered by the challenge and cost of collecting turbulence data. Even synthetic turbulence data is not easily accessible due to the complication of the forward model. Conventional simulators \cite{Hardie2017,Lau2021_sim,Dan17Sim} are either too slow or inaccurate due to oversimplifications of the model. Recently, an efficient physics-based simulator \cite{P2S} has been proposed to enable the generation of large-scale training datasets. However, the simulator in \cite{P2S} has only be applied to images with a certain low resolution grid. This greatly limits the accuracy of the turbulence model which requires a densely sampled grid.

In this paper, we make two main contributions to the blind multi-frame atmospheric turbulence mitigation problem. Firstly, we propose a few modifications to the simulator in \cite{P2S}. With the new modifications, we can now synthesize a large-scale training dataset that allows us to train deep neural networks with more diverse turbulence conditions. Secondly, we present the \textit{Turbulence Mitigation Transformer} (TMT), a multi-frame based transformer network for atmospheric turbulence mitigation. Our TMT is built upon a new attention mechanism, called the Temporal-Channel Joint Attention (TCJA). While the common self-attention modules mostly operate on spatial dimensions, TCJA operates jointly on both image channels and temporal dimensions, which allows us to effectively select and aggregate temporal information from multiple turbulence distorted frames. Our model achieves comparable performance with the conventional methods in terms of correcting the geometric distortions using only about 1/10 of the input frames compared with conventional methods. It also outperforms the SOTA networks designed for general video restoration.

\noindent Our contributions are summarized as follows:
\begin{itemize}
    \item We analyze the statistical properties of the turbulence simulation and make modifications to the Phase-to-Space (P2S) \cite{P2S} simulator. Our changes allow the simulator to generate more diverse, realistic turbulence conditions for images without loss of speed and precision.

    \item We introduce a progressive supervision architecture inspired by the physical property of turbulence. Our model can effectively remove distortion in the first stage and then restore the clear image in the second stage.

    \item We introduce temporal-channel joint attention (TCJA), a self-attention module to efficiently fuse the temporal information of atmospheric turbulence.

    \item We demonstrate good generalization of the proposed model on real world images. Our method shows superior performance to restore images on real world scenarios across a wide range of turbulence profiles.
\end{itemize}

\section{Related Work}
\noindent {\bf{Turbulence Mitigation Methods.}} Most of the conventional turbulence mitigation methods \cite{Anantrasirichai2013, mao_tci, Lau2017, Milanfar2013, Hirsch2010, Hardie2017recon} are based on a common multi-frame processing pipeline, where a large number of input frames are fused to generate a high-quality image. Before the fusion, it is often necessary to first extract a geometrically stable reference frame and align the input frames to the reference using optical flow or image registration methods. Blind deconvolution is commonly used at the end to remove any residual degradation, such as the diffraction blur. The conventional methods are often very slow, with the run time to process a $256\times256$ image ranging from a few seconds to over half an hour.

Due to the difficulties in obtaining real data, most of the recent deep learning methods are trained with synthetic turbulence images. \cite{P2S,Lau2021_sim} are two common choices for efficiently simulating large-scale turbulence datasets. Notably, the \textit{P2S} method \cite{P2S} is a physically accurate simulator that retains the necessary atmospheric turbulence statistics. It is adapted from \cite{chimitt2020simulating} by replacing computational expensive steps with tiny paralleled networks. The simulation technique used in this paper is an improvement over \cite{P2S} by relaxing several assumptions made in the original paper.

\noindent {\bf{General image and video restoration.}} Besides algorithms dedicated to atmospheric turbulence, image and video restoration methods designed for general restoration tasks \cite{zamir2021multi,zamir2021restormer,liang2022vrt,chan2022basicvsrpp} can also be potentially applied to the problem after re-training the model with appropriate data. These networks have often demonstrated competitive performance across the various image and video restoration tasks, such as denoising, image or video super-resolution, deraining, etc. Among these methods, the video restoration techniques are of particular interest to us, as they often adopt various strategies to utilize the temporal information from the input frames. At a high level, our work adopts a similar design to VRT \cite{liang2022vrt}, where a chunk of input frames is stored and processed parallelly. This often leads to better temporal information aggregation and higher throughput during both training and inference stages.

\noindent {\bf{Vision transformers.}} Since ViT \cite{dosovitskiy2020image} and DETR \cite{DETR} were introduced,  transformers are gradually becoming accepted as a powerful backbone for high-level vision tasks. Recent studies \cite{zamir2021restormer,Wang2022Uformer,liang2022vrt} have shown that the transformer-based models can also produce competitive results on image and video restoration tasks. Despite the promising performance, conventional transformers require memory that is quadratic in the size of input. The Swin Transformer \cite{liu2021Swin} proposed to split images into spatially shifting patches and apply self-attention within each patch to improve the efficiency. Video Swin Transformer \cite{liu2021video} extended the original 2D windows to 3D and this design was adopted by the Video Restoration Transformer (VRT) \cite{liang2022vrt} for video restoration tasks. Just like operations in \cite{ali2021xcit}, Restormer \cite{zamir2021restormer} alleviates the memory issue by constructing its self-attention purely in channel dimension, which is typically far smaller than the spatial dimension in image restoration tasks. Recently, \cite{zhou2022understanding} has shown channel attention could enhance the robustness of transformers. Our design is similar to the Restormer and we introduce a new shuffle mechanism for efficient channel fusion among different frames.

\section{Our method}
\subsection{Data synthesis with fast, physics-based simulator}
\noindent Assuming that the imaging system is incoherent, the observed image follows the equation in \cite{Hirsch2010}:
\begin{equation}\label{eq:1}
\boldsymbol{O} = \boldsymbol{h}(x,y) \circledast \boldsymbol{I} + \boldsymbol{n},
\end{equation}

\noindent where $\boldsymbol{I} \in \mathbb{R}^{H\times W}$ is the ideal clean image, $\boldsymbol{O} \in \mathbb{R}^{H\times W}$ is the degraded image,  $\circledast$ denotes convolution, and $\boldsymbol{n} \in \mathbb{R}^{H\times W}$ is the random noise. With no atmospheric distortions, $\boldsymbol{h} \in \mathbb{R}^{K\times K}$ represents the spatially-invariant diffraction blur kernel. However, in the case of atmospheric turbulence, $\boldsymbol{h}(x,y)$ is a \emph{spatially varying} blur kernel. Generally, there is a different blur kernel at each pixel location $(x, y)$, which we emphasize through our notation. The atmospheric distortions are random in nature, though can be statistically described by a set of parameters which contains physical information such as the aperture diameter of the imaging system $D$, the distance between the object and image $L$, the wavelength of light $\lambda$ and the Fried parameter $r_{0}$ \cite{Fried78}.

The atmospheric distortions described by $\boldsymbol{h}$ have two principal effects on an image: blur and geometric distortion. Since the correlation between them is small \cite{chimitt2020simulating} and geometric distortion is easier to implement, separating them into 2 sequential operators, the spatially varying blur kernels $\boldsymbol{b}$ and grid deformation field $\boldsymbol{\mathcal{G}}$, is preferred for efficiency. The expression can be:
\begin{equation}\label{eq:2}
\boldsymbol{O} = \boldsymbol{b}(x,y) \circledast \boldsymbol{\mathcal{G}}(\boldsymbol{I}) + \boldsymbol{n}
\end{equation}

\noindent or
\begin{equation}\label{eq:3}
\boldsymbol{O} = \boldsymbol{\mathcal{G}}(\boldsymbol{b}(x,y) \circledast \boldsymbol{I}) + \boldsymbol{n}.
\end{equation}

 The decomposition of $\boldsymbol{h}$ into $\boldsymbol{\mathcal{G}}$ and $\boldsymbol{b}(x,y)$ gives us the flexibility to treat the blur and geometric distortion separately. Comparing with equation \ref{eq:3}, equation \ref{eq:2} is a physically-based decomposition of degradation kernel $\boldsymbol{h}$, as widely adopted by many optical science papers as discussed by \cite{turb_Tutorial}.

However, solving the inverse problem according to the order of first inverting $\boldsymbol{b}(x,y)$ and then inverting $\boldsymbol{\mathcal{G}}$ in multi-frame restoration is challenging because inverting $\boldsymbol{b}(x,y)$ typically requires the features to be aligned first by inverting $\boldsymbol{\mathcal{G}}$. Therefore, when solving the inverse problem, we proceed with the approximated model by switching the order of the operators:
\begin{equation}\label{eq:4}
 \hat{\boldsymbol{I}}
 = \mathcal{B}^{-1} (\boldsymbol{\mathcal{G}}^{-1} (\boldsymbol{O}))
\end{equation}

\noindent where $\mathcal{B}$ represents the spatially varying convolution by $\boldsymbol{b}(x,y)$, and $\hat{\boldsymbol{I}}$ denotes the estimated image. This leads to using the more inaccurate equation \ref{eq:3} as widely accepted in turbulence mitigation \cite{Lau2021, yasarla2020learning}. We also use equation \ref{eq:3} in data synthesis and this forward model also enables good generalization because of the physically realistic blur kernels introduced by \cite{P2S}.

The scheme of data synthesis is shown in Figure \ref{fig:simulation}, we adopt a Zernike-based approach \cite{chimitt2020simulating, P2S} to match the turbulent statistics: we generate random samples in the Zernike space, then transform them to be in the image space, providing us $\boldsymbol{b}(x,y)$ and $\boldsymbol{\mathcal{G}}$. The Zernike polynomials are a set of basis functions which can be used to describe the phase distortions of a turbulence-distorted wave \cite{noll1976zernike}. The set of all Zernike coefficients for every pixel in an image forms a random vector field. One must randomly sample according to the statistics of this random field in order to perform a correct simulation. The random field generation in \cite{P2S} requires Cholesky decomposition of a $(HWN)\times(HWN)$ correlation matrix, where $H$ and $W$ are the width and height of the random field and $N=36$ is the number of Zernike coefficients. This imposes a serious limitation on the generation of turbulent statistics, which \cite{P2S} overcomes by interpolation from a small set of points. To overcome the speed and resolution bottleneck of the P2S simulator and generate more realistic data, we made the following modifications:

\begin{figure*}[ht]
  \centering
   \includegraphics[width=0.9\linewidth]{images/simulation.png}
  \caption{Pipeline of our physics-based simulation. The blur and geometric distortion caused by turbulence can be independently modeled.}
  \label{fig:simulation}
\end{figure*}

\noindent {\bf{Accelerate the correlation matrix generation}} The correlation for a single Zernike coefficient is a function of their displacement according to the Eq. 13 in \cite{takato1995spatial}, implying the property of WSS. However, there exist non-zero cross-correlations which do not adhere to the WSS property (i.e., not a function of difference of indices). We observe the cross-correlations of the 33 high-order Zernike coefficients are small in magnitude and decay sharply in a considerably small spatial range \cite{chimitt2020simulating}, suggesting it can be modeled approximately by a exponential decay function $\exp{(-\beta x)}$, where $x$ is the distance in pixel and $\beta$ describes the decay in spatial correlation due to the optical parameters. To allow for FFT generation of the random coefficient fields, we make an approximation to the the intermodal correlations, which we justify by their small statistical significance. The generation of the random fields for the higher order Zernike coefficients can then be described by: (1) Applying the spatial correlation via FFT to generate a set of independent Zernike coefficient fields; (2) Performing a pixel-wising mixing by Cholesky decomposition of the covariance matrix proposed by \cite{noll1976zernike}. With these changes we only require Cholesky decomposition of a $33\times33$ matrix, furthermore we can draw random field of Zernike coefficient in the same resolution as the input image. We then use P2S transform to get a physically accurate blurry image.

With respect to the spatial correlations of the second and third Zernike coefficients, the generation of the correlation kernel for each field requires the evaluation as \cite{chimitt2020simulating}:
\begin{equation}
\mathbb{E}[a_j^{\boldsymbol{x}_1}a_j^{\boldsymbol{x}_2}] = \frac{c_2}{2^{5/3}} \left(\frac{D}{r_{0}} \right)^{5/3} \left[ I_{0}(s) +\cos(2\varphi_0) I_{2}(s)\right]
\end{equation}

\noindent where $j \in \{2,3\}$, $a_j$ denotes the $j$-th Zernike coefficient, and $c_2 = 7.7554$. Additionally, $s=|\boldsymbol{x}_1 - \boldsymbol{x}_2|/D$, $\boldsymbol{x}$ is the spatial location with $\varphi_0$ as the angle between two points $\boldsymbol{x}_1, \boldsymbol{x}_2$. The terms $I_{0}$ and $I_{2}$ are integrals of Bessel functions of the first kind:
\begin{equation}
I_{i}(s) = \int_{0}^{\infty}\zeta^{-14/3} J_{i}(2s\zeta)J_{2}^{2}(\zeta) \text{d}\zeta.
\end{equation}

This evaluation is relatively time-consuming, motivating us to store a look-up table for the potential inputs. Since the range of the input value of the integral is limited, it results in a trivial loss of precision. After the aforementioned modifications, changing the turbulence profile for a $512 \times 512$ image only costs 17 ms rather than 4 minutes in \cite{P2S}. This allows us to change the turbulence parameters online during a real-time simulation.

\noindent {\bf{Employ various blur kernel sizes}} The size of blur kernel basis in P2S simulator was set to $33 \times 33$, while in real world scenarios, the size of the point spread functions varies based on the aperture of the optical system, imaging range and environment, but empirically, the blur kernels should be smaller than $33 \times 33$ to match the real-world videos. We randomly rescaled the kernel basis of P2S by various sizes (from $9\times 9$ to $33 \times 33$) but kept their original energy unchanged in the simulation process.

\noindent {\bf{Add temporal correlation}} Real world turbulence is temporally correlated, therefore, accurately simulating a video will require the distortions to be spatio-temporally correlated. We inject temporal correlations to the simulation process by correlating the initial random seed in the simulation. We use an $AR(1)$ process to generate the initial seed at the first stage. This allows for the random seed $\boldsymbol{W}_{t}$ at time $t$, which is then used to form the distortion and blur random fields, to be related to the previous realization by
\begin{equation}
\boldsymbol{W}_{t} = \alpha \boldsymbol{W}_{t-1} + \sqrt{1-\alpha^{2}}\boldsymbol{\epsilon}_{t}
\end{equation}
The term $\alpha$ is the temporal correlation ratio and $\boldsymbol{\epsilon}_{t} \sim \mathcal{N}(0, \boldsymbol{I})$.

\noindent {\bf{Simulating color images.}} As mentioned in \cite[Chapter 8.3.1]{Goodman_StatisticalOptics}, since the visible spectrum is very narrow, the difference of turbulence impact between the three color channels is negligible. Therefore, we follow the original \cite{P2S}, where the identical distortion is applied to the RGB channels.

\subsection{Turbulence mitigation transformer}
Inspired by the approximate forward process to generate turbulence effects, we designed a two-stage network to rectify the distortion and then remove the blur pattern on the image sequences. Our network is multi-in multi-out, the number of the input and output frames is 12, we retain full temporal resolution in every level of the network. For the distortion removal stage, we use a light-weight U-Net \cite{ronneberger2015u} with depth wise 3D convolution blocks. Recent works \cite{Li_2021_ICCV, anantrasirichai2022atmospheric} have shown that one doesn't need to use an overly complicated network to remove most of the distortion, this is also consistent with our preliminary experiments. However, those works only focused on distortion rectification and left the spatially variant blurs mostly unconsidered. Following the distortion removal stage, we propose the turbulence mitigation transformer to restore the clean image from blurry inputs. The overall architecture is shown in Figure \ref{fig:pipeline}, where the U-Net network on the left is the distortion removal module to measure the tilt-free image from blurry and distorted inputs, the UNet-like network on the right is the backbone of our TMT to further rectify the distortion and restore the clean image from the blurry estimations. In the second stage, we use multi-scale inputs: the original input is downsampled by 2 and 4 times to form 2 lower resolution images, they are processed by a 3D convolutional layer to extract features in multiple scales together with the original input. We concatenate those features with features from the previous encoder stages before feed them into the second and third encoder stages.

\begin{figure*}[ht]
  \centering
   \includegraphics[width=0.9\linewidth]{images/model.png}
  \caption{The scheme of our full TMT model. The input turbulence affect image sequence was first processed by a small UNet to remove the distortion, the last output is supposed to only have spatial variant blurs, it is then restored to clear image by our TMT backbone module. The number of blocks in the proposed TMT is set to be $M_{1}=M_{2}=M_{3}=M_{4}=2$ and $N_{0}=N_{1}=2, N_{2}=N_{3}=3, N_{4}=4$}
  \label{fig:pipeline}
\end{figure*}

 \paragraph{Temporal-channel joint attention module} Just like Restormer \cite{zamir2021restormer}, our TMT uses channel-wise attention to reduce memory requirements and enjoys more flexibility on spatial dimensions to avoid the quadratic complexity with respect to spatial dimensions. Moreover, it can handle different spatial shapes of input. To align and fuse temporal information, both self-attention inside one frame and mutual attention across different frames are jointly computed via TCJA. The implementation of attention inside the transformer block is shown in Figure \ref{fig:attention}.

On the right of this figure, we present two variants of feature aggregation before the computation of attention matrix. The variant $\textbf{a}$ is a plain fully connection on temporal axis following a depth-wise convolution module. But the features in different frames don't communicate in this scheme. Directly connecting all channels in the same spatial location of different frames will require large amount of parameters and memory space, to resolve this problem, we proposed the variant $\textbf{b}$: we implemented channel shuffle and partial connection inspired by the ShuffleNet \cite{zhang2018shufflenet}. We first split all channels into $C/g$ groups, where $C$ is the number of all channels and $g$ is the volume of each group, then instead of using a huge fully connection matrix on all channel and temporal tokens, we connect the temporal axis with channels in each group. By shuffling channels in every transformer block, all channels and temporal axes can communicate across several attention blocks. The detailed implementation can be described as following:

Given an input tensor $\textbf{X}\in \mathbb{R}^{C\times T\times H \times W}$, where C, T, H, W denotes the number of channels, frames, height, and width of each frame, we first filter it with a depth-wise 2D convolution (implemented by 3D convolution) to extract the spatial features and fuse them by point-wise convolution, which can also be recognized as the channel-wise fully connection operation.

For the variant \textbf{a}, we then do a temporal projection with a fully connection layer to produce $\textbf{Q}\in \mathbb{R}^{C\times H \times W \times T}$, $\textbf{K}\in \mathbb{R}^{C\times H \times W\times T}$, $\textbf{V}\in \mathbb{R}^{C\times H \times W\times T}$.

For the variant \textbf{b}, we need to split the channel dimension into $C/g$ groups. The dimension of the feature map becomes $(g, C/g, T, H, W)$. This is followed by a \textbf{channel shuffle} operation to transpose the first and second dimensions, the feature will then be in dimension $(C/g, g, T, H, W)$. We flatten the second and third dimensions and reshape to get $\hat{\textbf{X}} \in \mathbb{R}^{C/g \times H\times W \times gT}$. The following linear projection by matrix $\textbf{W} \in \mathbb{R}^{gT \times gT}$ enables feature communication among temporal information in different channels, it produces $\textbf{Q}\in \mathbb{R}^{C/g \times H \times W \times gT}$, $\textbf{K}\in \mathbb{R}^{C/g \times H \times W \times gT}$, $\textbf{V}\in \mathbb{R}^{C/g \times H \times W \times gT}$.

Before computing the attention matrix, we need to rearrange those tensors in both variants \textbf{a} and \textbf{b} into $query$ $\hat{\textbf{Q}}\in \mathbb{R}^{CT \times HW}$, $key$ $\hat{\textbf{K}}\in \mathbb{R}^{CT\times HW}$ and $value$ $\hat{\textbf{V}}\in \mathbb{R}^{CT \times HW}$. Then a standard multi-head self-attention can be implemented, we split $\hat{\textbf{Q}}$, $\hat{\textbf{K}}$ and $\hat{\textbf{V}}$ into $h$ heads, in the $i$-th head, we compute the embedding:
\begin{equation}
    \hat{\textbf{X}}_{i} = \text{softmax}\left(\frac{\hat{\textbf{Q}}_{i}\hat{\textbf{K}}_{i}^T}{\sqrt{d_{HW}}} \tau_{i}\right) \hat{\textbf{V}}_{i}
\end{equation}
Where learnable parameter $\tau_{i}$ is the temperature of the head $i$. The complexity of the channel attention is $O(\frac{HW(CT)^{2}}{h})$, which is linear to the spatial dimension and quadratic to the channel and temporal dimension.

\begin{figure*}[h]
  \centering
   \includegraphics[width=0.9\linewidth]{images/attention.png}
  \caption{The architecture of the transformer block. It consists of the channel-temporal joint attention (CTJA) module followed by a gated feed forward network. We designed two variants (a) and (b) of the CTJA block for $\text{TMT}_{a}$ and $\text{TMT}_{b}$, respectively. The only difference between those two is the channel shuffle operation.}
  \label{fig:attention}
\end{figure*}

\paragraph{Distortion removal module} Temporal alignment in video restoration is always an important component in network design. As the scene may contain moving objects, direct temporal fusion is infeasible to ensure the features to be aggregated are highly correlated. Pixel displacement can be captured by optical flow \cite{ranjan2017optical} or offsets in deformable convolution \cite{dai2017deformable, chan2022basicvsrpp}. For non-turbulent imaging, pixel displacement is typically measured between 2 frames, because it is often caused by motion and itself is not a degradation. However, in case of turbulent-distorted images, the pixel displacement may be a combination of both object motion and turbulent effects, thus making the bi-frame displacement to be a random field and the subsequent position of the ground truth object to be intractable via bi-frame alignment. Hence, we introduce an explicit, light-weight distortion removal module to measure the distortion-free, blur-only images. The network is a simple UNet which uses depth wise 3D convolution as its basic operation and ReLU as its activation. It has 4 levels, the input channel dimensions of them from shallow to deep are 64, 128, 256, 512 as the spatial scale downsampled by 2 across them. In the last 3 levels of the decoder, we output 3 sequences of warping grids to progressively rectify the input image. To enlarge the perceptual field of the network, we set the kernel sizes of the 4 encoder levels as 7, 7, 5 and 3. The downsampling is realized by maxpooling and the upsample module is a depth wise 3D transposed convolution layer.

\paragraph{Loss function}
The loss function $\boldsymbol{\mathcal{L}}$ of our work is
$$\boldsymbol{\mathcal{L}} = \sum_{l=1}^{3}\gamma_{l}\boldsymbol{\mathcal{L}}_{\text{char}}(\boldsymbol{I}^{\text{blur}}, \hat{\boldsymbol{\mathcal{G}}}_{l}(\boldsymbol{I}^{\text{in}})) +  \boldsymbol{\mathcal{L}}_{\text{char}}(\boldsymbol{I}^{\text{gt}}, \hat{\boldsymbol{O}})$$

Where $\boldsymbol{\mathcal{L}}_{\text{char}}$ denotes the Charbonnier loss \cite{charbonnier}. The $\boldsymbol{I}^{\text{blur}}$ is the synthesized blur-only, distortion-free image, it is used as ground truth in the distortion removal stage. $\hat{\boldsymbol{\mathcal{G}}}_{l}(\boldsymbol{I}^{\text{in}})$ represent the $l$-th level distortion removal process which warps the input to distortion-free image. $\gamma_{l}$ is the weight of loss in different level $l$, the value is selected empirically with $\gamma_{1} = 0.6$, $\gamma_{1} = 0.3$ and $\gamma_{1} = 0.1$. $\boldsymbol{I}^{\text{gt}}$ denotes the target clean image and $\hat{\boldsymbol{O}}$ is the restored image.

\section{Experiment}
\subsection{Datasets}
With the proposed modification on the P2S simulator, we created 2 synthetic datasets. Since most traditional methods for turbulence mitigation can only apply on static scenes, for direct comparison, we made a static scene dataset. We also created a dynamic scene dataset for restoring videos containing moving objects. Every single sequence of frames in our datasets have a unique turbulence profile, yet they can be synthesized in only one day with a NVIDIA 2080 Ti GPU, which is impossible using any other turbulence simulators. We trained and tested deep learning based restoration methods on synthetic data. To measure the generalization of our method, we used several existing real world datasets, most of them do not have the ground truth but we present several visual comparisons.
\paragraph{Synthetic datasets}
To restore turbulence effected images in as diverse scenarios as possible, the dataset should contain a variety of scenes. We used the place dataset \cite{zhou2017places} for synthesis for this reason. We randomly selected 9,017 images in the original dataset as input of the simulator, we generated 50 turbulence images and their associated distortion free images for each single input, resulting in 9,017 pairs of image sequences of static scenes. We split them into 7,499 pairs and 1,518 pairs, for training and testing, respectively. Figure \ref{fig:static_example} shows an example of testing on the synthetic image.

We also synthesized a video dataset for dynamic scenes. The source datasets for our dynamic scene data are the Sport Video in the Wild (SVW) dataset \cite{safdarnejad2015sports} and all ground truth videos used in TSRWGAN \cite{Jin2021NatureMI}. These videos are mixed, generating 4,684 samples with a total number of frames of 1,979,564. We generated 4,684 pairs of full turbulence and distortion free videos, then randomly split it into 3,500 videos for training and 1,184 for testing, keeping at most 120 frames per testing video.

During data synthesis, we set parameters mainly in an empirical way based on the observed real-world cases. Overall, we distribute turbulence parameters in 3 levels, they are weak, medium, and strong. Usually, longer distance requires a larger diameter of the aperture, and the smaller Fried parameter \cite{Fried78}, the turbulence is stronger. For a stronger turbulence level, we tend to use larger blur kernels. The correlation of pixel displacement can be computed by the parameters above, then we explicitly sample the correlation coefficient of higher-order aberration. Additionally, we empirically set temporal correlations to match the visual perception of the existing real-world data for better generalizability.

The detailed parameter sampling setting is shown in table \ref{table:parameter}. Each scene has 3 lines of settings, from top to bottom the turbulence strength becomes higher. Notice for static scenes, we prefer stronger turbulence so we set a lower probability to choose the first line which corresponds to weaker turbulence.

\paragraph{Real world data}
The ultimate goal of our work is to restore the real world image sequences. To examine for the generalization of the models, we purposely test all the models using data from different source, including the static patterns from OTIS \cite{gilles2017open} dataset, dynamic scene videos from \cite{Jin2021NatureMI} and the CLEAR \cite{anantrasirichai2022atmospheric} dataset. These datasets together contain a wide-range of turbulence conditions. We present samples of the testing result on each in Figure \ref{fig:OTIS}, \ref{fig:TSRWGAN_Data} and \ref{fig:WGAN_comparison} .

\begin{table*}[ht]
  \caption{Parameter sampling setting of the data synthesis. $\beta$ denotes the correlation coefficient of higher order aberration, $D$ denotes the aperture diameter of the imaging system, $D/r_{0}$ denotes the ratio of $D$ and the Fried parameter $r_{0}$. $U(a,b)$ denotes uniform distribution in the range $(a,b)$ and [] means random choice with equal probability.}
  \label{table:parameter}
  \centering
  \resizebox{\textwidth}{!}{\setlength\doublerulesep{0.5pt}
  \begin{tabular}{cccccccc}
    \toprule[1pt]\midrule[0.3pt]
    Dataset    & Probability & Kernel size choices & $\beta$ & $D$ (m) & $D/r_{0}$ choices & Distance (m) & Temporal correlation \\
    \midrule
    \multirow{3}{*}{Static Scene} & 0.2 & [33] & [0.05, 0.1, 0.2]  & $U(0.001, 0.005)$ & [0.5, 1, 1.2, 1.5] & $U(100, 400)$ & $U(0.2, 0.6)$\\
       & 0.4 & [33]  & [0.01, 0.02, 0.05, 0.1] & $U(0.04, 0.1)$ & [1, 1.5, 2] & $U(400, 800)$ & $U(0.2, 0.6)$\\
       &  0.4 &  [33] & [0.01, 0.02, 0.05, 0.1] & $U(0.1, 0.2)$ & [1.5, 2, 3] & $U(800, 1500)$ & $U(0.2, 0.6)$\\
    \midrule
    \multirow{3}{*}{Dynamic Scene} & $1/3$ & [9, 13, 15, 21] & [0.02, 0.05, 0.1, 0.2] & $U(0.001, 0.005)$ & [0.3, 0.6, 1, 1.2] & $U(50, 400)$ & $U(0.4, 0.8)$ \\
     & $1/3$ & [11, 17, 25, 33] & [0.02, 0.05, 0.1, 0.2] & $U(0.04, 0.1)$ & [0.3, 1, 1.5] & $U(400, 800)$ & $U(0.8, 0.95)$\\
     & $1/3$ & [15, 21, 27, 33] & [0.02, 0.05, 0.1, 0.2] & $U(0.1, 0.2)$ & [1, 1.5, 2, 2.5] & $U(800, 2000)$ & $U(0.88, 0.95)$ \\
    \midrule[0.3pt]\bottomrule[1pt]
  \end{tabular}}
\end{table*}

\subsection{Training and evaluation}
We compare the proposed TMT model with the TSRWGAN \cite{Jin2021NatureMI}, a recent multi-frame turbulence mitigation method, and two state-of-the-art general video restoration methods VRT \cite{liang2022vrt} and BasicVSR++ \cite{chan2022basicvsrpp, chan2022generalization}.

All models are trained separately for static and dynamic scene datasets. We trained our model using the Adam optimizer with the Cosine Annealing scheduler \cite{loshchilov2016sgdr}. The VRT, BasicVSR++, and TSRWGAN models are all trained under the same settings as the original paper. For VRT and BasicVSR++, we use the configurations for motion deblurring as it is the closest resemble to our task. We set the input and output as 12 frames for TMT, VRT and BasicVSR++ for fair comparison. We fine-tuned the pre-trained TSRWGAN model on our dataset to maximize the performance and hence we keep the original 15 input frames setting, which shouldn't put the method at any disadvantage.

We first trained the distortion removal module for 400K iterations, the batch size was set to 4. We used Adam optimizer \cite{kingma2014adam} and with $\beta_{1}=0.9$ and $\beta_{2}=0.99$. The learning rate was initialized as $2\times10^{-4}$ then gradually decreased to $1\times 10^{-6}$ by the Cosine Annealing schedualer \cite{loshchilov2016sgdr}. We fixed the distortion removal module later and trained the TMT backbone for 800K iterations within the scheme of online training. The same Adam optimizer was applied and the learning rate follows the Cosine Annealing scheme with initial learning rate $1 \times 10^{-4}$ and final learning rate $1 \times 10^{-6}$. The model was implemented in PyTorch with a single NVIDIA A100 GPU, it requires 10 days to finish all training.
The training of VRT model mostly follows the setting for motion deblurring in the original paper \cite{liang2022vrt}. We used the scheme of online training as well, it is feasible since both TMT and VRT only used Layer normalization \cite{ba2016layer}, which is insensitive to batch size. The initial learning rate was set to $1 \times 10^{-4}$ and the model was trained from scratch for 800K iterations as well. The training was on a single NVIDIA A100 GPU and took 40 days to finish.

We fine-tuned the TSRWGAN for 400K iterations on our datasets since it was trained for the same task as ours. We used the same training scheme as the original work. The training took 5.5 days on a single NVIDIA GeForce RTX 3090 GPU to converge.
For the BasicVSR++ \cite{chan2022basicvsrpp}, we trained it from scratch for 400K iterations with batch size 4, the other training scheme is the same as the original paper.

During training, we set the input size of both TMT and VRT as $12\times 128 \times 128$, input size of TSRWGAN as $15\times 200 \times 200$ and input size of BasicVSR++ as $12\times 256 \times 256$.

For all experiments, we used the same data augmentation strategy: random Gamma correction, random saturation, random cropping, followed by a random combination of horizontal/vertial flipping and $90^{\circ}$ rotation. Since the real turbulence images are often captured using telephoto lenses with large F-number, a high ISO is required to improve the image quality, which inevitably leads to non-negligible image noise. Therefore, we also injected Gaussian random noise from $ \mathcal{N}(0, \beta \boldsymbol{I})$, where $\beta \sim U(0,0.0004]$ during training to enhance the robustness of models.

The models are evaluated on our synthetic testing set, we measured PSNR and SSIM in both RGB space and YCbCr space for quatitative comparison. The testing results are shown in Table \ref{table:tab1} and \ref{table:tab2}. We also provide a visual comparison in Fig. \ref{fig:static_example} to show the advance of our method.

\begin{table}
\centering
\caption{Comparison on the {\bf{STATIC}} scene dataset.}
\label{table:tab1}
\setlength\doublerulesep{0.5pt}
\begin{tabular}{ccccc}
\toprule[1pt]\midrule[0.3pt]
Methods & PSNR & SSIM & $\text{PSNR}_{Y}$ & $\text{SSIM}_{Y}$\\
\midrule
BasicVSR++ \cite{chan2022generalization}  & 26.5055 & 0.8121 & 27.9068  & 0.8238 \\
\midrule
TSRWGAN \cite{Jin2021NatureMI}   & 25.2888 & 0.7784 & 27.3573  & 0.8017 \\
\midrule
VRT \cite{liang2022vrt}    & 27.4556 & 0.8280 & 28.9101 & 0.8451\\
\midrule
TMT [ours]   & \textbf{27.5422} & \textbf{0.8320} & \textbf{29.0011} & \textbf{0.8487} \\
\midrule[0.3pt]\bottomrule[1pt]
\end{tabular}
\end{table}

\begin{table}
\centering
\caption{Comparison on the {\bf{DYNAMIC}} scene dataset.}
\label{table:tab2}
\setlength\doublerulesep{0.5pt}
\begin{tabular}{ccccc}
\toprule[1pt]\midrule[0.3pt]
Methods & PSNR & SSIM & $\text{PSNR}_{Y}$ & $\text{SSIM}_{Y}$\\
\midrule
BasicVSR++ \cite{chan2022generalization} & 27.0231 & 0.8073 & 28.3055  & 0.8314 \\
\midrule
TSRWGAN \cite{Jin2021NatureMI}   & 26.3262 & 0.7957 & 28.0491  & 0.8245 \\
\midrule
VRT \cite{liang2022vrt}    & 27.6114 & 0.8300 & 29.1128  & 0.8496\\
\midrule
TMT [ours]   & \textbf{27.7419} & \textbf{0.8323} & \textbf{29.2510}  & \textbf{0.8520}\\
\midrule[0.3pt]\bottomrule[1pt]
\end{tabular}
\end{table}

\begin{figure}
    \captionsetup[subfloat]{font=scriptsize,farskip=1pt}
    \centering
  \subfloat[Input\label{fig2:Input}]{%
       \includegraphics[width=0.49\linewidth]{images/syn_static/input.jpg}}
    \hfill
  \subfloat[Ground truth\label{1b}]{%
        \includegraphics[width=0.49\linewidth]{images/syn_static/gt.png}}
    \\
  \subfloat[{TMT [ours]} \label{1c}]{%
        \includegraphics[width=0.49\linewidth]{images/syn_static/TMT.png}}
  \hfill
  \subfloat[TSRWGAN \cite{Jin2021NatureMI}]{%
        \includegraphics[width=0.49\linewidth]{images/syn_static/WGAN.png}}
    \\
  \subfloat[VRT \cite{liang2022vrt}\label{1e}]{%
        \includegraphics[width=0.49\linewidth]{images/syn_static/VRT.png}}
    \hfill
  \subfloat[BasicVSR++ \cite{chan2022basicvsrpp} \label{1f}]{%
        \includegraphics[width=0.49\linewidth]{images/syn_static/Basic.png}}
  \caption{Example of testing results on our synthetic static scene dataset. (a). Input (b).Ground truth (c). Output of {\bf{TMT [ours]}} (d).Output of TSRWGAN \cite{Jin2021NatureMI} (e). Output of VRT \cite{liang2022vrt} (f). Output of BasicVSR++ \cite{chan2022basicvsrpp}}
  \label{fig:static_example}
\end{figure}

\subsection{Ablation study}
\begin{table*}[ht]
  \caption{Quantitative comparison on the synthetic testing set for different design choices of the TMT. $\text{TMT}_{a}$ and $\text{TMT}_{b}$ denote the TMT model with variants a and b of the attention module, respectively. Here, ``warp'' means the distortion removal module, ``SS'' and ``MS'' are aberrations of single-scale input and multi-scale input. FLOPs is measured with input size of $12 \times 208 \times 208$}
  \label{table:ablation}
  \centering
  \setlength\doublerulesep{0.5pt}
  \begin{tabular}{ccccccccccc}
    \toprule[1pt]\midrule[0.3pt]
    Dataset & \multicolumn{4}{c}{Static Scenes} & \multicolumn{4}{c}{Dynamic Scenes}  & \multicolumn{2}{c}{Computation Consumption}   \\
    \cmidrule(r){1-1} \cmidrule(r){2-5} \cmidrule(r){6-9} \cmidrule(r){10-11}
    Methods     & PSNR & SSIM & $\text{PSNR}_{Y}$ & $\text{SSIM}_{Y}$  & PSNR & SSIM & $\text{PSNR}_{Y}$ & $\text{SSIM}_{Y}$ & $\#$ Params (M) & FLOPs (G)\\
    \midrule
    $\text{SS-TMT}_{a}$ w.o. warp  & 27.2782 & 0.8221 & 28.7316  & 0.8398 & 27.5635 & 0.8258 & 29.0644  & 0.8459  & \textbf{22.75}  & 1490 \\
    $\text{SS-TMT}_{b}$ w.o. warp  & 27.3092 & 0.8235 & 28.7626  & 0.8411 & 27.5864 & 0.8265 & 29.0886  & 0.8465 & 22.95  & 1514 \\
    $\text{SS-TMT}_{a}$  & 27.3432 & 0.8257 & 28.8013 & 0.8431 & 27.5779 & 0.8282 & 29.0863  & 0.8487 & 24.87  & 1676 \\
    $\text{SS-TMT}_{b}$  & 27.3836 & 0.8266 & 28.8396 & 0.8440 & 27.6051 & 0.8281 & 29.1110  & 0.8485 & 25.07  & 1700 \\
    $\text{MS-TMT}_{a}$ w.o. warp  & 27.4718 & 0.8291 & 28.9243  & 0.8452 & 27.6637 & 0.8301  & 29.1685  & 0.8501  & 23.70  & \textbf{1206} \\
    $\text{MS-TMT}_{b}$ w.o. warp  & 27.5003 & 0.8301 & 28.9602  & 0.8468 & 27.6841 & 0.8310 & 29.1894  & 0.8506 & 23.92 & 1304 \\
    $\text{MS-TMT}_{a}$  & 27.5215 & 0.8307 & 28.9781 & 0.8469 & 27.7239  & \textbf{0.8329}  & 29.2337   & \textbf{0.8526}  & 25.82  & 1392 \\
    \midrule[0.3pt]
    $\text{MS-TMT}_{b}$ [Proposed]  & \textbf{27.5422} & \textbf{0.8320} & \textbf{29.0011} & \textbf{0.8487} & \textbf{27.7419} & 0.8323 & \textbf{29.2510}  & 0.8520 & 26.04  & 1490\\
    \midrule[0.3pt]\bottomrule[1pt]
  \end{tabular}
\end{table*}

\paragraph{Influence of number of input frames}
To estimate the ground truth of an image sequence under turbulence, the network must perceive a large enough area in spatial-temporal domain to estimate the property of the degradation field, especially for samples having high temporal correlation. If the correlation factor is 1, all frames have the same degradation, the multi-frame reconstruction will be collapsed into single-frame problem, which is extremely hard to solve as revealed by \cite{vint2020analysis}. To
understand how the number of frames would affect the restoration, we trained our $\text{TMT}_{b}$ using 4, 8, 12, 20 input frames with our synthetic static scene dataset. The testing performance is shown in figure \ref{fig:nframe}. The results shows the performance of the network can be improved with more input frames. We also note that, to make fair comparison with strong baseline such as \cite{liang2022vrt}, which requires large GPU memory space so that it is hard to be deployed with more than 12 frames, we kept 12-frame setting in this paper.

\begin{figure}[ht]
    \centering
    \includegraphics[width=0.9\linewidth]{images/frame_all.png}
    \caption{Number of frames and performance}
    \label{fig:nframe}
\end{figure}

\paragraph{Distortion removal module}
To evaluate the value of the distortion removal module, we trained our TMT backbone to directly restore images without the preprocessing of removing the distortion. The result is shown in Table \ref{table:ablation}. As it shows, the distortion removal module can boost the overall performance by $0.02 \sim 0.08$dB. It is more essential in static scenes because the statistic of the distortion is harder to estimate the dynamic scenes due to the hybrid of motion. We also have two variants of the transformer block. As mentioned in section 3, the shuffle connection between temporal and channel axes facilitates better temporal fusion in the following attention operation.

\paragraph{Influence of channel shuffle operation}
The purpose of channel shuffle operation is to facilitate the communication of channels across different frames in the TCJA module. The comparison in Table \ref{table:ablation} demonstrates the effectiveness of the channel shuffle operation, the performance benefits $0.02 \sim 0.03$ dB from the shuffle operation compared to the plain connection.

\paragraph{Influence of multi-scale input in the TMT backbone}
For general U-Net like design \cite{ronneberger2015u}, only single input is required to feed into the encoder. There are several works suggesting feeding inputs in multiple scales/resolutions could help the network to learn more efficiently \cite{8099518, ji2022multi, zamir2021multi}. In the proposed paper, we utilize multi-scale input by simply downsampling the raw input into lower resolutions. To demonstrate the advantage of this operation, we also trained our TMT with single-scale input: the input dimension of each level in the encoder is remain the same, only the input feature is changed from the concatenation of outputs produced by the last level and the 3D convolution layer to solely from the last level. The comparison in Table \ref{table:ablation} demonstrates the effectiveness of the multi-scale input, the PSNR of networks with multi-scale input is $\sim 0.15$ dB higher than that with single-scale input, while the model consumption doesn't increase.

\subsection{Generalization on real world data}
With our physically-based simulation process, samples of real-world turbulence can be viewed as the interpolation of samples generated by our simulation. Other simulation methods are either much more time-consuming or inaccurate \cite{P2S}. Two recent data-driven works TSRWGAN \cite{Jin2021NatureMI} and complex CNN \cite{anantrasirichai2022atmospheric} also used synthetic data to train neural networks, however, since the limitation of their simulation methods, their generalizability is limited. However, quantitative comparisons among models on real world images are hard because most real-world image sequences don't have ground truth. Despite this problem, direct visual comparison is still possible and the superiority of our method over others in terms of generalizability in real world image restoration can be demonstrated.

\paragraph{Visual comparison among models trained on our datasets}
In the CVPR 2022 UG2+ Challenge, a new long-range turbulence dataset is released for benchmarking turbulence mitigation algorithms \cite{UG2}. This dataset consists of 100 image sequences of text patterns, which were captured from 300 meters away during hot weather. We tested our trained models on this dataset and got promising results, furthermore, we noticed our network performs better than other networks under heavy turbulence. Several samples are given in Fig. \ref{fig:text}.

Other than the turbulence text dataset, we also tested all models on an earlier OTIS dataset \cite{gilles2017open} and TSRWGAN's real-world test set. The OTIS has 16 sequences of static patterns in different scales and turbulence levels. The TSRWGAN's real-world test set has 27 different dynamic scenes in relatively mild turbulence strength. Since the ground truth is unavailable, we only show some visual comparisons in \ref{fig:OTIS} and \ref{fig:TSRWGAN_Data}, from which one can easily conclude that models trained on our dataset generalize well on a broad range of real-world turbulence conditions and among them, our model could restore images in better visual quality than others.

\begin{figure*}[ht]
    \captionsetup[subfloat]{farskip=2pt, font=scriptsize}
    \centering
    \includegraphics[width=0.42\linewidth]{images/text/A3/pg2_input.png}%
    \hfill
    \includegraphics[width=0.14\linewidth]{images/text/A3/basic_2.png}%
    \hfill
    \includegraphics[width=0.14\linewidth]{images/text/A3/WGAN_2.jpg}%
    \hfill
    \includegraphics[width=0.14\linewidth]{images/text/A3/VRT_2.jpg}%
    \hfill
    \includegraphics[width=0.14\linewidth]{images/text/A3/TMT_2.jpg}%
    \\
    \vspace{2pt}
    \includegraphics[width=0.42\linewidth]{images/text/A3/pg24_input.png}%
    \hfill
    \includegraphics[width=0.14\linewidth]{images/text/A3/basic_24.png}%
    \hfill
    \includegraphics[width=0.14\linewidth]{images/text/A3/WGAN_24.jpg}%
    \hfill
    \includegraphics[width=0.14\linewidth]{images/text/A3/VRT_24.jpg}%
    \hfill
    \includegraphics[width=0.14\linewidth]{images/text/A3/TMT_24.jpg}%
    \\
  \subfloat[Input sequences]{%
      \includegraphics[width=0.42\linewidth]{images/text/A3/pg96_input.png}}
    \hfill
  \subfloat[BasicVSR++ \cite{chan2022basicvsrpp}]{%
    \includegraphics[width=0.14\linewidth]{images/text/A3/basic_96.png}}
   \hfill
  \subfloat[TSRWGAN \cite{Jin2021NatureMI}]{%
        \includegraphics[width=0.14\linewidth]{images/text/A3/WGAN_96.jpg}}
    \hfill
  \subfloat[VRT \cite{liang2022vrt}]{%
        \includegraphics[width=0.14\linewidth]{images/text/A3/VRT_96.jpg}}
    \hfill
  \subfloat[{TMT [ours]}]{%
        \includegraphics[width=0.14\linewidth]{images/text/A3/TMT_96.jpg}}
  \caption{\textbf{Test on sequences of real-world text data \cite{UG2}}. The 3 lines of images from top to buttom are input and output of the 2nd, 24th, 96th sequences in that dataset. Column (a). Input sequences (b).Output of BasicVSR++ \cite{chan2022basicvsrpp} (c). Output of TSRWGAN \cite{Jin2021NatureMI} (d).VRT \cite{liang2022vrt} (e).TMT [ours]}
  \label{fig:text}
\end{figure*}

\begin{figure*}[htbp]
    \captionsetup[subfloat]{farskip=2pt, font=scriptsize}
    \centering
  \subfloat[Input]{%
    \includegraphics[width=0.33\linewidth, height=0.165\linewidth]{images/OTIS/OTIS_in.png}}
    \hfill
  \subfloat[Ground truth]{%
    \includegraphics[width=0.33\linewidth, height=0.165\linewidth]{images/OTIS/OTIS_gt.png}}
   \hfill
  \subfloat[{TMT [ours]}]{%
    \includegraphics[width=0.33\linewidth, height=0.165\linewidth]{images/OTIS/OTIS_TMT.png}}
   \\
  \subfloat[TSRWGAN \cite{Jin2021NatureMI}]{%
    \includegraphics[width=0.33\linewidth, height=0.165\linewidth]{images/OTIS/OTIS_TSRWGAN.png}}
   \hfill
  \subfloat[VRT \cite{liang2022vrt}]{%
    \includegraphics[width=0.33\linewidth, height=0.165\linewidth]{images/OTIS/OTIS_VRT.png}}
   \hfill
  \subfloat[BasicVSR++ \cite{chan2022basicvsrpp}]{%
    \includegraphics[width=0.33\linewidth, height=0.165\linewidth]{images/OTIS/OTIS_basicvsr.png}}
  \caption{Visual comparison on the OTIS dataset \cite{gilles2017open}. We show the result of Pattern 15 and 16 which contain the strongest turbulence effect.}
  \label{fig:OTIS}
\end{figure*}

\begin{figure*}[!htbp]
    \centering
    \captionsetup[subfloat]{farskip=2pt, font=scriptsize}
  \subfloat[Input]{%
    \includegraphics[width=0.33\linewidth, height=0.13\linewidth]{images/car/input_new.png}}
    \hfill
  \subfloat[Original WGAN \cite{Jin2021NatureMI}]{%
    \includegraphics[width=0.33\linewidth, height=0.13\linewidth]{images/car/WGAN_original_new.png}}
   \hfill
  \subfloat[Fine-tuned WGAN]{%
    \includegraphics[width=0.33\linewidth, height=0.13\linewidth]{images/car/WGAN_ours_new.png}}
    \\
  \subfloat[BasicVSR++ \cite{chan2022basicvsrpp}]{%
    \includegraphics[width=0.33\linewidth, height=0.13\linewidth]{images/car/Basic_new.png}}
    \hfill
  \subfloat[VRT \cite{liang2022vrt}]{%
    \includegraphics[width=0.33\linewidth, height=0.13\linewidth]{images/car/VRT_new.png}}
   \hfill
  \subfloat[{TMT [ours]}]{%
    \includegraphics[width=0.33\linewidth, height=0.13\linewidth]{images/car/TMT_new.png}}
  \caption{Visual comparison on real turbulence data from TSRWGAN's dataset \cite{Jin2021NatureMI}. Model of (b) is provided by authors of \cite{Jin2021NatureMI} and model for (c) is fine-tuned on our data. From the comparison between (b) and (c), we find the robustness to relatively strong turbulence has been enhanced using our data, despite the artifacts of railing on the window still exists. Notice (f) has the least turbulence left while preserving the most details.}
  \label{fig:TSRWGAN_Data}
\end{figure*}

\paragraph{Compare with simulation methods used by TSRWGAN}
The algorithm simulation tool TSRWGAN \cite{Jin2021NatureMI} used is \cite{repasi2011computer}. This method generates physically grounded distortion and spatially invariant blur kernels, but higher-order aberrations are not simulated. \cite{Jin2021NatureMI} generated more physically synthetic data by artificial heat sources to produce turbulence effect in a short distance. However, this physical simulation tends to generate highly correlated degradation with a weak blur. Although their synthetic dataset is reasonably big and they simulated various turbulence levels, we observe that the released model does not generalize well on CLEAR's real-world dataset, OTIS \cite{gilles2017open} and the text dataset \cite{UG2} where the turbulence effect is stronger than that in the synthetic data used by TSRWGAN. We fine-tuned TSRWAGN on our synthetic data and we show the significant improvement of the TSRWGAN model in the text dataset, the OTIS dataset and CLEAR's dataset in Fig. \ref{fig:WGAN_comparison}.

\begin{figure*}[htbp]
    \captionsetup[subfloat]{farskip=2pt, font=scriptsize}
    \centering
  \subfloat[Input]{%
    \includegraphics[width=0.45\linewidth]{images/text/WGAN_input_82.png}}
    \hfill
  \subfloat[Trained with TSRWGAN's data]{%
    \includegraphics[width=0.27\linewidth]{images/text/WGAN_original_82.jpg}}
   \hfill
  \subfloat[Trained with our data]{%
    \includegraphics[width=0.27\linewidth]{images/text/WGAN_ours82.jpg}}
    \hfill
    \\
  \subfloat[Input]{%
      \includegraphics[width=0.33\linewidth]{images/plane_supp/input.png}}
    \hfill
  \subfloat[Trained with TSRWGAN's data]{%
    \includegraphics[width=0.33\linewidth]{images/plane_supp/WGAN_original.png}}
   \hfill
  \subfloat[Trained with our data]{%
    \includegraphics[width=0.33\linewidth]{images/plane_supp/WGAN_ours.png}}
    \hfill
  \caption{Test TSRWAGN’s generator on real-world images. (a-c). On text dataset \cite{UG2} (d-f). On CLEAR's dataset (a). Input sequence (d). Single frame sample of input (b, e). Restoration result of TSRWGAN trained with their original data (c, f). Restoration result of TSRWGAN trained with our synthetic datasets. Clearly our datasets directly facilitates the generalization on real world images}
  \label{fig:WGAN_comparison}
\end{figure*}

\paragraph{Compare with the simulation method used in CLEAR} Recent learning-based extension \cite{anantrasirichai2022atmospheric} of CLEAR \cite{Anantrasirichai2013} employed 9 predefined point spread function of atmospheric turbulence introduced in \cite{Hirsch2010} to simulate the degradation. Those PSFs are not physically grounded and the diversity of the turbulence conditions is very limited. Therefore, the restoration performance of \cite{anantrasirichai2022atmospheric} on some real-world data is not ideal. Since the trained model of CLEAR for video restoration is not released, we were not able to train it on our dataset to demonstrate the potential improvement. Despite this, based on our discussions and comparisons in the previous sections, we can still conclude that using our proposed simulation method has a more positive impact on the reconstruction.

\begin{figure*}[htbp]
    \captionsetup[subfloat]{farskip=1pt, font=scriptsize}
    \centering
  \subfloat[Input: $\#12$ frame]{%
    \includegraphics[width=0.248\linewidth]{images/hillhouse/hillhouse_distorted2.png}}
    \hfill
  \subfloat[Conventional CLEAR \cite{Anantrasirichai2013}]{%
    \includegraphics[width=0.248\linewidth]{images/hillhouse/hillhouse_restored2.png}}
   \hfill
  \subfloat[Mao et al. \cite{mao_tci}]{%
    \includegraphics[width=0.248\linewidth]{images/hillhouse/hillhouse_tci2.png}}
   \hfill
  \subfloat[{TMT [ours]}]{%
    \includegraphics[width=0.248\linewidth]{images/hillhouse/TMT2.png}}
  \caption{Comparison with representative conventional methods on the real-world image sequence. Our model is trained on the synthetic static scene dataset. The presented input frame is the 12th in the whole sequence. {\bf{Zoom in for better view of details}}.}
  \label{fig:Conventional_clear}
\end{figure*}

\subsection{Pixel displacement rectification by distortion removal module}
Our distortion removal module is robust to different levels of blurs caused by turbulence, it can be used alone without the TMT backbone in some cases. For examples, there are scenarios where the field of view is large and the optical path is short, in this case, the blurry effect will be weak and the turbulence is mainly dominated by distortion. As shown in figure \ref{fig:tilt_street}, if we directly apply the distortion removal model only, we can already give a visually satisfying reconstruction.

\begin{figure*}[htbp]
    \captionsetup[subfloat]{farskip=2pt, font=scriptsize}
    \centering
    \includegraphics[width=0.33\linewidth]{images/A2/input_init.png}%
    \hfill
    \includegraphics[width=0.33\linewidth]{images/A2/CLEAR_init.png}%
    \hfill
    \includegraphics[width=0.33\linewidth]{images/A2/tilt_init.png}%
    \\
  \subfloat[]{%
      \includegraphics[width=0.33\linewidth, height=0.3465\linewidth]{images/A2/input_processed2.png}}
    \hfill
  \subfloat[]{%
    \includegraphics[width=0.33\linewidth, height=0.3465\linewidth]{images/A2/CLEAR_processed2.png}}
   \hfill
  \subfloat[]{%
    \includegraphics[width=0.33\linewidth, height=0.3465\linewidth]{images/A2/tilt_processed2.png}}
  \caption{Application of our distortion removal module on CLEAR's real-world data. Images on the top are full resolution images and below are the same patches in 3 frames starting from them. (a). Input (b). Processed by CLEAR (c). Processed by our distortion removal module only}
  \label{fig:tilt_street}
\end{figure*}

\subsection{Comparison with conventional methods}
We make a comparison with two state-of-the-art conventional turbulence mitigation methods \cite{Anantrasirichai2013, mao_tci}. Due to the slow processing speed of the conventional ones, we only compare the methods on a random subset of 100 scenes from our synthetic static dataset. The comparison on the 100 scenes is shown in table \ref{table:comparison} Conventional methods could still work because our synthetic data has similar properties compared to real-world turbulence degraded images. For example, the lucky effect can be synthesized with this simulator \cite{P2S}. It is worth mentioning that we used 50 frames as input for conventional methods and 12/15 for learning-based methods.

\begin{table}[ht]
  \caption{Quantitative comparison with conventional turbulence mitigation methods on the 100 scenes dataset. $\text{PSNR}_{Y}$ and $\text{SSIM}_{Y}$ denotes the luminance channel in YCbCr color space.}
  \vspace{1ex}
  \label{table:comparison}
  \centering
  \resizebox{0.45\textwidth}{!}{\setlength\doublerulesep{0.5pt}
  \begin{tabular}{cccccc}
    \toprule[1pt]\midrule[0.3pt]
    \multicolumn{2}{c}{Methods}     & PSNR & SSIM & $\text{PSNR}_{Y}$ & $\text{SSIM}_{Y}$\\
    \midrule
    \multirow{2}{*}{Conventional} & Mao et al. \cite{mao_tci}    & 27.4154 & 0.8174 & 28.9619 & 0.8497\\
    & CLEAR \cite{Anantrasirichai2013}    & 26.9969 & 0.8295 & 28.4906 & 0.8525\\
    \midrule
    \multirow{4}{*}{Deep learning} &
    BasicVSR++ \cite{chan2022generalization}  & 28.0811 & 0.8410 & 29.5210  & 0.8588 \\
    & TSRWGAN \cite{Jin2021NatureMI}   & 25.8139 & 0.7977 & 27.9967  & 0.8451\\
    & VRT \cite{liang2022vrt}    & 28.3454 & 0.8512 & 29.8080 & 0.8664\\
    & TMT [ours]   & \textbf{28.4310} & \textbf{0.8535} & \textbf{29.8968} & \textbf{0.8695}\\
    \midrule[0.3pt]\bottomrule[1pt]
  \end{tabular}}
\end{table}

Additionally, we compare our proposed method with \cite{Anantrasirichai2013} and \cite{mao_tci} on real-world images. A qualitative comparison is presented in figure \ref{fig:Conventional_clear}. CLEAR \cite{Anantrasirichai2013} used a subsample mechanism to select images in good quality among all 75 frames of the input sequence, then used a complex wavelet-based algorithm to fuse selected frames. After those operations, the output is still subjected to explicit denoising, while our method just blindly used the first 12 frames as input and the noise is mostly compressed. Mao et al. \cite{mao_tci} is based on lucky fusion algorithm, a popular framework in conventional turbulence mitigation task, it has an explicit denoising stage after the fusion so the result is more blurry. It can be seen that our method produced more natural and high-contrast reconstruction without further lose of details.

\subsection{Computation consumption}
We present the inference-time computing budget of all models we trained in this section. The measurement is based on a single NVIDIA 2080 Ti GPU to match a very common deployment environment, we set the spatial size of the input image sequence to be $540 \times 960$. When the full size can't fit into the GPU, we split the input into largest possible patches. The splitting should be overlapped to reduce artifacts, for example, the maximum patches of VRT is $160 \times 160$, so we need $4\times 7$ patches to reconstruct a full resolution image sequence. From the table \ref{table:computation_comparison}, although VRT has the closest performance in terms of PSNR and SSIM to ours, TMT is much more efficient. Because of the efficiency, we can use larger patches as input during testing and get more precise estimation on turbulence conditions.
\begin{table}[ht]
  \caption{Comparison of computational consumption in inference. We show 2 settings of the distortion removal model, the first line is when it is a part of TMT, the second line shows the consumption when it's split into its maximum allowable patches on our testing environment. Note the TSRWGAN is multi-in, single-out, we need to run it 15 times to get reconstruction of 15 frames.}
  \label{table:computation_comparison}
  \centering
  \resizebox{0.98\linewidth}{!}{\setlength\doublerulesep{0.5pt}
  \begin{tabular}{cccc}
    \toprule[1pt]\midrule[0.3pt]
     \multirow{2}{*}{Methods}     & $\#$ parameters   & Patch size & FLOPs  \\
        & M  & $\text{T}\times \text{H} \times \text{W}$ &  $\text{G} \times \#$ times
    \\
    \midrule
    BasicVSR++ \cite{chan2022generalization}  & 9.76  & $ 12\times 256 \times 256$  & $466.0 \times 12$  \\
    TSRWGAN \cite{Jin2021NatureMI}   & 46.28  & $ 15\times 540 \times 960$  & $2243 \times 15$ \\
    VRT \cite{liang2022vrt}  &  18.32 & $12\times 160 \times 160$ & $3324 \times 28$ \\
    \midrule
    $\text{TMT}_{a}$ & 25.82 & $12\times 208 \times 208$ & $1392 \times 15$ \\
    $\text{TMT}_{b}$ & 26.04  & $12\times 208 \times 208$ & $1490 \times 15$  \\
    \multirow{2}{*}{Distortion removal module} & 2.120 & $12\times 208 \times 208$ &  $186.0 \times 15$ \\
    & 2.120 & $12\times 320 \times 512$ &  $704.5 \times 4$ \\
    \midrule[0.3pt]\bottomrule[1pt]
  \end{tabular}}
\end{table}

\section{Conclusion}
In this paper, we presented a complete and generalizable solution for a multi-frame blind atmospheric turbulence mitigation problem. We refined a physics-based simulation in \cite{P2S} by approximating the spatially varying field using a wide sense stationary field. We further proposed an efficient transformer-based network to restore image sequences degraded by turbulence and our method can adapt to various scenes and turbulence conditions, our model surpassed the state of art models designed for general video restoration tasks. Compared to conventional methods on turbulence mitigation, the proposed method has superior performance in visual quality, speed, and robustness.

\section*{Acknowledgments and Disclosure of Funding}
The research is based upon work supported in part by the Intelligence Advanced Research Projects Activity (IARPA) under Contract No. 2022‐21102100004, and in part by the National Science Foundation under the grants CCSS-2030570 and IIS-2133032. The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies, either expressed or implied, of IARPA, or the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for governmental purposes notwithstanding any copyright annotation therein.

\bibliographystyle{IEEEtran}