# Radar Camera Fusion via Representation Learning in Autonomous Driving

Xu Dong Binnan Zhuang Yunxiang Mao Langechuan Liu $\dagger$
XSense.ai
11010 Roselle Street, San Diego, CA 92121

###### Abstract

Radars and cameras are mature, cost-effective, and robust sensors and have
been widely used in the perception stack of mass-produced autonomous driving
systems. Due to their complementary properties, outputs from radar detection
(radar pins) and camera perception (2D bounding boxes) are usually fused to
generate the best perception results. The key to successful radar-camera
fusion is accurate data association. The challenges in radar-camera
association can be attributed to the complexity of driving scenes, the noisy
and sparse nature of radar measurements, and the depth ambiguity from 2D
bounding boxes. Traditional rule-based association methods are susceptible to
performance degradation in challenging scenarios and failure in corner cases.
In this study, we propose to address rad-cam association via deep
representation learning, to explore feature-level interaction and global
reasoning. Concretely, we design a loss sampling mechanism and an innovative
ordinal loss to overcome the difficulty of imperfect labeling and to enforce
critical human reasoning. Despite being trained with noisy labels generated by
a rule-based algorithm, our proposed method achieves a performance of 92.2% F1
score, which is 11.6% higher than the rule-based teacher. Moreover, this data-
driven method also lends itself to continuous improvement via corner case
mining.

00footnotetext: $\dagger$ indicates corresponding author<EMAIL_ADDRESS>
##  1 Introduction

LiDAR, radar, and camera are the three main sensory modalities employed by the
perception system of an autonomous driving vehicle. Though LiDAR-based 3D
object detection is very popular in high-level autonomy, its wide adoption is
still limited by some unsolved issues. First, LiDAR is prone to adversarial
conditions (e.g. rainy weather); second, current LiDAR systems still exhibit
prohibitively high maintenance need and cost; third, the mass-production of
LiDAR is not ready to meet the growing demand.

Figure 1: An illustration of the associations between radar detections (radar
pins) and camera detections (2D bounding boxes). The context of the scene is
illustrated in the top picture, with the image captured by the camera along
with the detected bounding boxes and the projected radar pins (shown as
numbered blue circles). The bottom picture highlights the association
relationships (shown in red lines ) between radar pins and bounding boxes. The
orange line in the middle denotes uncertain association relationship, which
will be explained later.

An automotive millimeter-wave radar can also provide a certain level of
geometrical information with relatively precise range and speed estimates.
Moreover, as a widely-adopted sensor in automobiles for decades, radar is
relatively robust, low-cost, and low-maintenance. The fusion between radar and
camera combines radar’s geometrical information and camera’s appearance and
semantic information, which is still the mainstream perception solution in
many practical autonomous driving and assisted driving systems.

Traditionally, the radar-camera fusion is achieved by the combination of rule-
based association algorithms and kinematic model based tracking. The key is
data association between radar and camera detections. The noisy and sparse
nature of radar detection and the depth ambiguity from a mono camera makes
such association problem very challenging. Traditionally, the association
process is hand-crafted based on minimizing certain distance metrics along
with some heuristic rules. It not only requires a large amount of engineering
and tuning but is also hard to adapt to ever-growing data.

An emerging solution is to use learning-based methods to replace the rule-
based radar-camera fusion. The latest advances focus on direct 3D object
detection with the combined radar and camera data as the input [16, 25, 28].
These approaches all rely on LiDAR-based ground-truth to build the link
between radar and camera. This is feasible on most public datasets such as
nuScenes [4], Waymo [29] etc. However, it cannot be applied to a large fleet
of commercial autonomous vehicles, often equipped with only radars and
cameras. In this study, we propose a scalable learning-based framework to
associate radar and camera information without the costly LiDAR-based ground-
truth.

Our goal is to find representations of radar and camera detection results,
such that matched pairs are close and unmatched ones are far. We convert the
detection results into image channels and combine them with the original image
to feed into a convolutional neural networks (CNN), namely, AssociationNet.
Training is performed based on imperfect labels obtained from a traditional
rule-based association method. A loss sampling mechanism is introduced to
mitigate false labels. To further boost the performance, we guide the
reasoning logic of AssociationNet by adding a novel ordinal loss. The proposed
AssociationNet significantly outperforms the rule-based method through global
reasoning about the actual scene.

Our main contributions are summarized as follows:

  * •

We proposed a scalable learning-based radar-camera fusion framework without
using ground-truth labels from LiDAR, which is suitable for building a low-
cost, production-ready perception system in real autonomous driving
applications.

  * •

We designed a loss sampling mechanism to alleviate the impact of the label
noise, and also invented an ordinal loss to enforce critical association logic
into the model for performance enhancement.

  * •

We developed a robust model via representation learning, which is capable of
handling various challenging scenarios, and also outperforms the traditional
rule-based algorithm by 11.6% regarding the F1 score.

##  2 Related Work

###  2.1 Sensor Fusion

Traditionally, different sensory modules process their data separately. A
downstream sensor fusion module augments the sensory outputs (typically
detected objects) to form a more comprehensive understanding of the
surroundings. Such an object-level fusion method is the mainstream approach
[9, 17, 19, 12, 31] and is still widely used on many Advanced Driver
Assistance Systems (ADAS). In object-level fusion, object detection is
independently performed on each sensor, and the fusion algorithm combines such
object detection results to create so-called global tracks for kinematic
tracking [1].

Data association is the most critical and challenging task in object-level
fusion. Precise association can easily lead to 3D object detection and
multiple objects tracking solutions [1, 5]. Traditional approaches tend to
manually craft various distance metrics to represent the similarities between
different sensory outputs. Distance minimization [9] and other heuristic rules
are applied to find the associations. To handle the complexity and
uncertainty, probabilistic models are also sometimes adopted in the
association process [2].

###  2.2 Learning-Based Radar-Camera Fusion

The learning-based radar-camera fusion algorithms can be primarily categorized
into three groups, data-level fusion, feature-level fusion, and object-level
fusion. The data-level fusion and feature-level fusion combine the radar and
camera information at the early stage [28, 14] and the middle stage [16, 7,
26], respectively, but both directly perform 3D object detection. Hence, they
rely on LiDAR to provide ground-truth labels during training, which prohibits
their usage to autonomous vehicles without LiDAR.

The learning-based object-level fusion remains under-explored due to the
lesser amount of information contained in the detected results. In this study,
our proposed method belongs to this category in that we focus on associating
radar and camera detection results. Thus, our method is more compatible with
traditional sensor fusion pipeline. On the other hand, our method also
directly takes the raw camera and radar data for further performance
enhancement.

###  2.3 CNN for Heterogeneous Data

The tremendous success of CNN on structured image data inspires its
application to many other types of heterogeneous data, such as sensor
parameters, point clouds, and association relationships between two groups of
data [27]. In order to get compatible with CNN, a popular approach is to adapt
the heterogeneous data into a form of pseudo-images. Examples include encoding
camera intrinsic into images with normalized coordinates and field of view
maps [11], projecting radar data into the image plane to form new image
channels [6, 28], and the various forms of projection-based LiDAR point-cloud
representations [30, 23]. We adopted a similar approach in this study to
handle the heterogeneous radar and camera outputs.

###  2.4 Representation Learning

Representation learning has been considered as the key to understanding
complex environments and problems [3, 20, 18]. Representation learning has
been widely used in many natural language processing tasks such as word
embedding [24], and many computer vision tasks, such as image classification
[8], object detection [13], and keypoint matching [10]. In this study, we aim
at learning a vector in the high-dimensional feature space as the
representation for each object in the scene, in order to establish the
interactions between objects as well as enable global reasoning about the
scene.

##  3 Problem Formulation

We use a front-facing camera and a millimeter-wave mid-range front radar for
the proposed radar-camera fusion, while the approach can be easily generalized
to 360 perception with proper hardware setups. The camera intrinsic and the
extrinsics of both sensors are obtained through offline calibration. The radar
and camera operate asynchronously at 20Hz and 10Hz, respectively. Both the
radar and the camera are front-facing and have the field-of-views (FOVs) of
120 degrees and 52 degrees, respectively. The output of the camera sensor at
each frame is an RGB image with the size of 1828 pixels (width) by 948 pixels
(height), whereas the output of the radar sensor at each frame is a list of
processed points with many attributes (conventionally referred to as radar
pins). Since the radar used here performs internal clustering, each output
radar pin is on the object level.111The proposed fusion technique also applies
to lower level detection, e.g., radar locations. There are several tens of
radar pins per-frame depending on the actual scene and traffic. The attributes
of each radar pin are listed in Table 1. There are two noteworthy
characteristics of the radar pins. First, we only consume the 2D position
information in the Bird’s-Eye View (BEV) without the elevation angle, due to
poor resolution and large measurement noise in the elevation dimension.
Second, each radar pin either corresponds to a movable object (cars, cyclists,
pedestrians etc.) or an interfering static structure such as a traffic sign, a
street light, or a bridge.

In this study, we focus on associating 2D bounding boxes detected from a
camera image to radar pins detected in the corresponding radar frame. With
precise associations, many subsequent tasks like 3D object detection and
tracking become much easier if not trivial.

Table 1: The Features of Each Radar Pin Feature | Explanation
---|---
object id | the id of the radar pin
obstacle prob |  the probability of the existence of an obstacle being detected by the radar pin
position x |  the x coordinate of the position of the detected obstacle in radar frame
position y |  the y coordinate of the position of the detected obstacle in radar frame
velocity x |  the velocity of the detected obstacle along the x coordinate in radar frame
velocity y |  the velocity of the detected obstacle along the y coordinate in radar frame
Table 2: The Features of Each 2D Bounding Box Feature | Explanation
---|---
center x |  the x coordinate in the image plane of the center of the bounding box
center y |  the y coordinate in the image plane of the center of the bounding box
height |  the height of the bounding box in the image plane
width |  the width of the bounding box in the image plane
category |  the category of the detected moving object, including sedan, suv, truck, bus, bicycle, tricycle, motorcycle, person, and unknown

##  4 Methods

Our proposed method mainly consists of a preprocessing step to align radar and
camera data, a CNN-based deep representation learning network, AssociationNet,
and a postprocessing step to extract representations and make associations. An
overview of the method is shown in Fig. 2 and details are explained the
following sections.

Figure 2: An overview of AssociationNet. Process a illustrates how the radar
pins and 2D bounding boxes are first projected into the camera image plane and
then produce a pseudo-image. Process b illustrates how the final pseudo-image
is composed by concatenating all the features of radar pins, bounding boxes,
and the original RGB camera image. The pseudo-image is then fed into a neural
network to learn high-level semantic representations. Process c illustrates
how the learned representation vectors for objects are finally extracted from
the feature-map generated in the last layer of the neural network.

###  4.1 Radar and Camera Data Preprocessing

Temporal and spatial alignment is performed in the preprocessing stage. For
each camera frame, we look for the nearest radar frame to perform data
alignment. We align the nearest radar frame to the time instant of the camera
frame, by moving the radar pin locations forward/backward along the time axis
under a constant velocity assumption. After the temporal alignment, the radar
pins are further transformed from the radar coordinate to the camera
coordinate using the known extrinsics. All the attributes of the aligned radar
pins will be used in AssociationNet.

Each camera frame is first fed into a 2D object detection network to produce a
list of 2D bounding boxes corresponding to the movable objects in the scene.
The output attributes for each detected 2D bounding box are displayed in Table
2. Though the network used in this study is an anchor-based RetinaNet [22]
network, any 2D object detector will serve the purpose. After preprocessing, a
list of temporally and spatially aligned radar pins and bounding boxes will be
ready for association.

###  4.2 Deep Association by Representation Learning

We employ AssociationNet to learn a semantic representation (or a descriptor)
of each radar pin and each bounding box. Under such representation, a pair of
matched radar pin and bounding box will look similar, in the sense that the
distance between the learned representations are small. An overview of the
general process is shown in Fig. 2.

To leverage the powerful CNN architecture, we project each radar pin and 2D
bounding box into the image plane to generate a pseudo-image, with each
attribute occupying an independent channel. Specifically, each bounding box is
assigned to the pixel location of its center. Each radar pin is assigned to
the pixel location which is obtained through projecting its 3D location into
the image plane using the camera intrinsic. The process is illustrated in
Process a of the Fig. 2. Next, we concatenate the raw RGB camera image with
the corresponding pseudo-image to incorporate the rich pixel-level
information. AssociationNet is then applied to perform representation
learning.

As shown in Fig. 3, the network consists of a ResNet-50 [15] as the backbone,
a Feature Pyramid Network [21] for feature-map decoding, along with two extra
layers to restore the size of the output feature-map to the original input
size. The output feature-map contains the high-level semantic representations
of radar pins and bounding boxes. As each radar pin or bounding box has a
unique pixel location in the feature-map, we extract the representation vector
of each on the output feature-map at its corresponding pixel location. The
process is illustrated in Process c of the Fig. 2.

The input pseudo-image contains seven radar pin channels, four bounding box
channels, and three raw camera image RGB channels. The radar pin channels
include object-id, obstacle-prob, position-x, position-y, velocity-x,
velocity-y222Positions and velocities used here are under camera coordinate,
as it is after the spatial alignment step in the preprocessing., and a heatmap
to indicate the projected pixel location. The bounding box channels include
height, width, category, and also a heatmap to indicate the pixel location.
The output feature-map contains 128 channels, resulting in the dimension of
the representation vector to be 64 for each radar pin and bounding box.

Figure 3: The architecture of the neural network. It consists of a ResNet-50
as the backbone, a feature pyramid for feature decoding, and two extra layers
to restore the feature-map size. The feature-map in the last layer will be
used to extract the representation vectors as shown in the Process c of the
Fig. 2 Figure 4: An overview of the process of obtaining final associations
from the learned representation vectors.

The obtained representation vector captures the semantic meaning of each radar
pin and each bounding box in a high dimension space. If a radar pin and a
bounding box come from the same object in real world, we treat the pair of
radar pin and bounding box as a positive sample, otherwise it is considered as
a negative sample. We try to minimize the distance between the representation
vectors of any positive sample and maximize the distance between the
representation vectors of any negative sample. Based on such logic, we design
loss functions according to the association ground-truth labels. We pull
together the representation vectors of positive samples with the following
pull loss:

$L_{pull}=\frac{1}{n_{pos}}\sum_{(i_{1},i_{2})\in\mathbb{POS}}\max(0,\lVert
h_{i_{1}}-h_{i_{2}}\rVert-m_{1})$ (1)

And we push apart the representation vectors of negative samples with the
following push loss:

$L_{push}=\frac{1}{n_{neg}}\sum_{(i_{1},i_{2})\in\mathbb{NEG}}\max(0,m_{2}-\lVert
h_{i_{1}}-h_{i_{2}}\rVert)$ (2)

Here $\mathbb{POS}$ and $\mathbb{NGE}$ are the set of positive samples and the
set of negative samples, respectively in each frame; $n_{pos}$ and $n_{neg}$
are the total number of associations in $\mathbb{POS}$ and $\mathbb{NGE}$
respectively; $(i_{1},i_{2})$ denotes the $i^{th}$ association pair consisting
of radar pin $i_{1}$ and bounding box $i_{2}$; $h_{i_{1}}$ and $h_{i_{2}}$
denotes the learned representation vectors; and $m_{1}$ and $m_{2}$ are the
thresholds for the desired distances of representations among positive
associations and negative associations, which were preset to be 2.0 and 8.0 in
our experiments.

During inference, we calculate the Euclidean distance between the
representation vectors of all possible radar-pin-bounding-box pairs. If the
distance falls below a certain threshold, the radar pin and the bounding box
will be considered as a successful association. More details of the inference
process will be explained later.

####  4.2.1 Loss Sampling

The association labels used for supervising the learning process is ultimately
from the traditional rule-based method, which are far from 100% accurate and
contain some noise. To mitigate the impact of the inaccurate labels, we first
purify the labels by applying some simple filters to remove low-confidence
associations, which increases the precision in the remaining association
labels at the cost of the undermined recall. During the push loss calculation
in the training of AssociationNet, instead of exhausting all negative pairs (a
pair of a radar pin and a bounding box that is not present in the association
labels), we only sample a fraction of those to be used for push loss
calculation to alleviate pushing apart positive pairs by mistake. The number
of sampled negative pairs is set to be equal to the number of positive ones at
each frame.

####  4.2.2 Ordinal Loss

One particular kind of error made by AssociationNet is that it could violate
the simple ordinal rule, i.e., a farther radar pin associates to a closer
bounding box and a closer radar pin associates to a farther bounding box. To
solve this issue, an ordinal loss is introduced.

Denote the $y$ coordinate of bounding box $i$’s bottom edge as
$y_{\text{max}}^{i}$ and the corresponding depth in 3D world as $d^{i}$. For
any two bounding boxes on the same image we have the property:

$y_{\text{max}}^{i}>y_{\text{max}}^{j}\iff d^{i}>d^{j}$ (3)

The ordering of the objects in the 3D world can be interpreted as the relative
vertical ordering of the bottom edges of the corresponding bounding boxes.
Hence, we design an additional ordinal loss to enforce self-consistency
according to the ordinal rule, which is written as:

$L_{ord}=\frac{2}{\widehat{n_{pos}}\cdot(\widehat{n_{pos}}-1)}\cdot\\\
\sum_{\begin{subarray}{c}(i_{1},i_{2})\in\widehat{\mathbb{POS}}\\\
(j_{1},j_{2})\in\widehat{\mathbb{POS}}\end{subarray}}\sigma(-(d^{i_{1}}-d^{j_{1}})\cdot(y_{\text{max}}^{i_{2}}-y_{\text{max}}^{j_{2}})),$
(4)

where $\widehat{\mathbb{POS}}$ denotes the set of predicted positive
associations and $\widehat{n_{pos}}$ is the size of the set; $i_{1}$ and
$i_{2}$ represent the associated radar pin and bounding box of the $i^{th}$
predicted association in $\widehat{\mathbb{POS}}$; similarly, $j_{1}$ and
$j_{2}$ represent the associated radar pin and bounding box of the $j^{th}$
predicted association in $\widehat{\mathbb{POS}}$; $d^{*}$ represents the
depth of a radar pin in camera coordinate, and $y_{\text{max}}^{*}$ represents
the $y$ coordinate of a bounding box’s bottom edge; $\sigma$ is the sigmoid
function to smooth the loss values.

Finally, the total loss is calculated as:

$L_{tot}=L_{pull}+L_{push}+w_{ord}\cdot L_{ord},$ (5)

where the $w_{ord}$ is the adjustable weight to balance losses.

###  4.3 Training and Inference

The AssociationNet was trained with a batch size of 48 frames at four NVIDIA
GeForce RTX 2080 Ti GPUs. The SGD optimizer was used for training at a total
of 10K iterations. The learning rate was set to be $10^{-4}$ initially, and
then was decreased by a factor of 10 at the end of 8K iterations and 9K
iterations, respectively.

At the inference time, the representation vectors for all radar pins and
bounding boxes are first predicted using the trained model. An affinity matrix
is then calculated, where each matrix element corresponds to the distance
between the representations of a radar pin and a bounding box. In reality,
each bounding box may be associated with multiple radar pin333This is usually
the case where the corresponding vehicles are of large sizes, such as trailer
trucks and buses., while each radar pin can only match to at most one bounding
box. As a result, we devise each radar pin to get associated to the bounding
box with the smallest distance in the affinity matrix. Lastly, the improbable
associations with a distance larger than a threshold are filtered out, which
usually consists of radar pins from interfering static objects. The whole
inference process is summarized in Fig. 4.

###  4.4 Evaluation

The predicted associations are compared against human-annotated ground-truth
associations in the test dataset. We use precision, recall, and F1 score as
the metrics for evaluating the performance.

In some very complicated scenes, correctly associating all radar pins and
bounding boxes is very challenging even for human annotators. Therefore, we
mark those plausible but dubious associations as uncertain in the evaluation
process. An example is shown in Fig. 5. For those uncertain associations, they
are counted as neither positive nor negative associations, which will be
excluded from both true and false positive predictions.

Figure 5: An illustration of radar pins, bounding boxes, and their association
relationships under BEV perspective. This BEV image corresponds to the same
scene as displayed in Fig. 1. Each grid in the image represents a 10-meter-
by-10-meter square in physical space. The bounding boxes are represented as
solid cycles in this image. The location of a bounding box is estimated by the
Inverse Projective Mapping (IPM) method from the bounding box’s center, to
provide a rough reference for its real 3D location. A truncated frustum
accompanying each bounding box is also plotted, for better assisting human
curators to determine the association relationships555The frustum is
calculated also by the IPM method, from the two side-edges of each bounding
box. According to projective geometry, the real object detected by a bounding
box has to be within the bounding box’s frustum, and hence the possibly
matched radar pins as well. We truncated the frustums for the ease of
visualizing. The widths of each frustum at the truncated positions are set to
be one meter and five meters, respectively. As the physical width of a vehicle
is most likely to be within the range, the possibly matched radar pins also
tend to lie within the truncated frustum..

Figure 6: Examples of AssociationNet predictions. Here, the red solid lines
represent the true-positive associations; and the pink solid lines represent
predicted positive associations but labeled as uncertain in the ground-truth.
In the second example, the added green lines represent the false-positive
predictions; and the added black lines represent the false-negative
predictions. Also note that each bounding box on the left corresponds to a
solid circle with the same color on the right.

##  5 Experiments and Discussion

###  5.1 Dataset

The AssociationNet was trained and evaluated on an in-house dataset with 12
driving sequences collected by a testing fleet, which consists of 14.8 hours
of driving in various driving scenarios, including highway, urban, and city
roads. The radar and camera were synchronized at 10 Hz initially and further
downsampled to 2 Hz, in order to reduce the temporal correlation among
adjacent frames. Eleven sequences out of the twelve were used for training
with the other one left for the test. Therefore, there are 104,314
synchronized radar and camera frames in the training dataset, and 2,714 in the
test dataset. For the training data, the association labels were generated by
a traditional rule-based algorithm with additional filtering to increase the
precision. For the test data, we manually curated the labels with human
annotators to obtain high-quality ground-truth.

###  5.2 Effect of Loss Sampling

We studied the effect of loss sampling on the AssociationNet’s performance.
Experiments were conducted with no sampling (meaning that all the negative
pairs present in the label are used for push loss calculation), and loss
sampling with different sampling ratios. The sampling ratio is defined as the
ratio between the number of positive pairs and the number of negative pairs at
each frame. The result is shown in Table 3. We can see that the best sampling
ratio is 1:1 with the loss sampling mechanism, which boosts the performance by
1.1% in terms of the F1 score.

Table 3: The Effect of Loss Sampling Sample Ratio | Performance
---|---
| Precision / Recall / F1
no sampling | 0.896 / 0.925 / 0.911
1:2 | 0.901 / 0.931 / 0.916
1:1 | 0.906 / 0.939 / 0.922
2:1 | 0.899 / 0.933 / 0.915
3:1 | 0.899 / 0.929 / 0.914

###  5.3 Effect of Ordinal Loss

The effect of the ordinal loss is shown in Table 4. The ordinal loss can
facilitate both precision and recall to some degree. With the optimal loss
weight, the performance is boosted by 1.8% in terms of the F1 score.

Table 4: The Effect of Ordinal Loss Loss Weight | Performance
---|---
$w_{ord}$ | Precision / Recall / F1
0.0 | 0.897 / 0.912 / 0.904
0.5 | 0.897 / 0.923 / 0.910
1.0 | 0.899 / 0.931 / 0.915
2.0 | 0.906 / 0.939 / 0.922
5.0 | 0.889 / 0.918 / 0.903

###  5.4 Comparison with Rule-Based Algorithm

We compared the performance of AssociationNet with the traditional rule-based
algorithm, as shown in Table 5. Notably, though the traditional rule-based
algorithm was used to generate association labels to supervise the training of
AssociationNet, AssociationNet significantly outperforms the rule-based
alternative. This demonstrates the inherent robustness of learning-based
algorithms in handling complex scenarios.

Table 5: Comparison with Rule-based Algorithm Algorithm | Performance
---|---
| Precision / Recall / F1
Rule-based | 0.890 / 0.736 / 0.806
Learning-based | 0.906 / 0.939 / 0.922

###  5.5 Visualization

Examples of the predicted associations are shown in Fig. 6. Despite multiple
big trucks present in both examples, AssociationNet correctly predicted their
associations, which demonstrates the robustness of the algorithm. On the other
hand, in the second example, there are two bounding boxes incorrectly
associated, with one bounding box having no predicted associations and the
other associated to a wrong radar pin. The two bounding boxes correspond to
vehicles at the very far range. The mistakes are largely due to the small
sizes of the objects in the camera image and also the heavy occlusions.

##  6 Conclusion

In this work, we developed a scalable learning-based radar-camera fusion
algorithm, without using LiDAR for ground-truth labels generation. Such a
solution has many practical merits at the current technological stage,
including low cost, low maintenance, high reliability, and more importantly,
readiness for mass-production. We employed deep representation learning to
tackle the challenging association problem, with the benefits of enabled
feature-level interaction and global reasoning. We also designed a loss
sampling mechanism and a novel ordinal loss to mitigate the impact of label
noise and enforce critical human logic into the learning process. Although
imperfect labels generated by a traditional rule-based algorithm were used to
train the network, our proposed algorithm outperforms the rule-based teacher
by 11.6% in terms of the F1 score.

## References

  * [1] Michael Aeberhard, Stefan Schlichtharle, Nico Kaempchen, and Torsten Bertram.  Track-to-track fusion with asynchronous sensors using information matrix fusion for surround environment perception.  IEEE Transactions on Intelligent Transportation Systems, 13(4):1717–1726, 2012.
  * [2] Yaakov Bar-Shalom, Fred Daum, and Jim Huang.  The probabilistic data association filter.  IEEE Control Systems Magazine, 29(6):82–100, 2009.
  * [3] Yoshua Bengio, Aaron Courville, and Pascal Vincent.  Representation learning: A review and new perspectives.  IEEE transactions on pattern analysis and machine intelligence, 35(8):1798–1828, 2013.
  * [4] Holger Caesar, Varun Bankiti, Alex H. Lang, Sourabh Vora, Venice Erin Liong, Qiang Xu, Anush Krishnan, Yu Pan, Giancarlo Baldan, and Oscar Beijbom.  nuscenes: A multimodal dataset for autonomous driving.  arXiv preprint arXiv:1903.11027, 2019.
  * [5] Josip Ćesić, Ivan Marković, Igor Cvišić, and Ivan Petrović.  Radar and stereo vision fusion for multitarget tracking on the special euclidean group.  Robotics and Autonomous Systems, 83:338–348, 2016.
  * [6] Simon Chadwick, Will Maddern, and Paul Newman.  Distant vehicle detection using radar and vision.  In 2019 International Conference on Robotics and Automation (ICRA), pages 8311–8317. IEEE, 2019.
  * [7] Shuo Chang, Yifan Zhang, Fan Zhang, Xiaotong Zhao, Sai Huang, Zhiyong Feng, and Zhiqing Wei.  Spatial attention fusion for obstacle detection using mmwave radar and vision sensor.  Sensors, 20(4):956, 2020.
  * [8] Ting Chen, Simon Kornblith, Mohammad Norouzi, and Geoffrey Hinton.  A simple framework for contrastive learning of visual representations.  In International conference on machine learning, pages 1597–1607. PMLR, 2020.
  * [9] Hyunggi Cho, Young-Woo Seo, BVK Vijaya Kumar, and Ragunathan Raj Rajkumar.  A multi-sensor fusion system for moving object detection and tracking in urban driving environments.  In 2014 IEEE International Conference on Robotics and Automation (ICRA), pages 1836–1843. IEEE, 2014.
  * [10] Daniel DeTone, Tomasz Malisiewicz, and Andrew Rabinovich.  Superpoint: Self-supervised interest point detection and description.  In Proceedings of the IEEE conference on computer vision and pattern recognition workshops, pages 224–236, 2018.
  * [11] Jose M. Facil, Benjamin Ummenhofer, Huizhong Zhou, Luis Montesano, Thomas Brox, and Javier Civera.  Cam-convs: Camera-aware multi-scale convolutions for single-view depth.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), June 2019.
  * [12] Fernando Garcia, Pietro Cerri, Alberto Broggi, Arturo de la Escalera, and José María Armingol.  Data fusion for overtaking vehicle detection based on radar and optical flow.  In 2012 IEEE Intelligent Vehicles Symposium, pages 494–499. IEEE, 2012.
  * [13] Ross Girshick, Jeff Donahue, Trevor Darrell, and Jitendra Malik.  Rich feature hierarchies for accurate object detection and semantic segmentation.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR), June 2014.
  * [14] Xiao-peng Guo, Jin-song Du, Jie Gao, and Wei Wang.  Pedestrian detection based on fusion of millimeter wave radar and vision.  In Proceedings of the 2018 International Conference on Artificial Intelligence and Pattern Recognition, pages 38–42, 2018.
  * [15] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 770–778, 2016.
  * [16] Vijay John and Seiichi Mita.  Rvnet: deep sensor fusion of monocular camera and radar for image-based obstacle detection in challenging environments.  In Pacific-Rim Symposium on Image and Video Technology, pages 351–364. Springer, 2019.
  * [17] Naoki Kawasaki and Uwe Kiencke.  Standard platform for sensor fusion on advanced driver assistance system using bayesian network.  In IEEE Intelligent Vehicles Symposium, 2004, pages 250–255. IEEE, 2004.
  * [18] Alexander Kolesnikov, Xiaohua Zhai, and Lucas Beyer.  Revisiting self-supervised visual representation learning.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 1920–1929, 2019.
  * [19] Dirk Langer and Todd Jochem.  Fusing radar and vision for detecting, classifying and avoiding roadway obstacles.  In Proceedings of Conference on Intelligent Vehicles, pages 333–338. IEEE, 1996.
  * [20] Yann LeCun, Yoshua Bengio, and Geoffrey Hinton.  Deep learning.  nature, 521(7553):436–444, 2015.
  * [21] Tsung-Yi Lin, Piotr Dollár, Ross Girshick, Kaiming He, Bharath Hariharan, and Serge Belongie.  Feature pyramid networks for object detection.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 2117–2125, 2017.
  * [22] Tsung-Yi Lin, Priya Goyal, Ross Girshick, Kaiming He, and Piotr Dollár.  Focal loss for dense object detection.  In Proceedings of the IEEE international conference on computer vision, pages 2980–2988, 2017.
  * [23] Gregory P. Meyer, Ankit Laddha, Eric Kee, Carlos Vallespi-Gonzalez, and Carl K. Wellington.  Lasernet: An efficient probabilistic 3d object detector for autonomous driving.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), June 2019.
  * [24] Tomas Mikolov, Kai Chen, Greg Corrado, and Jeffrey Dean.  Efficient estimation of word representations in vector space.  arXiv preprint arXiv:1301.3781, 2013.
  * [25] Ramin Nabati and Hairong Qi.  Rrpn: Radar region proposal network for object detection in autonomous vehicles.  In 2019 IEEE International Conference on Image Processing (ICIP), pages 3093–3097. IEEE, 2019.
  * [26] Ramin Nabati and Hairong Qi.  Centerfusion: Center-based radar and camera fusion for 3d object detection.  In Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision, pages 1527–1536, 2021.
  * [27] Alejandro Newell and Jia Deng.  Pixels to graphs by associative embedding.  arXiv preprint arXiv:1706.07365, 2017.
  * [28] Felix Nobis, Maximilian Geisslinger, Markus Weber, Johannes Betz, and Markus Lienkamp.  A deep learning-based radar and camera sensor fusion architecture for object detection.  In 2019 Sensor Data Fusion: Trends, Solutions, Applications (SDF), pages 1–7. IEEE, 2019.
  * [29] Pei Sun, Henrik Kretzschmar, Xerxes Dotiwalla, Aurelien Chouard, Vijaysai Patnaik, Paul Tsui, James Guo, Yin Zhou, Yuning Chai, Benjamin Caine, Vijay Vasudevan, Wei Han, Jiquan Ngiam, Hang Zhao, Aleksei Timofeev, Scott Ettinger, Maxim Krivokon, Amy Gao, Aditya Joshi, Yu Zhang, Jonathon Shlens, Zhifeng Chen, and Dragomir Anguelov.  Scalability in perception for autonomous driving: Waymo open dataset, 2019\.
  * [30] Yue Wang, Alireza Fathi, Abhijit Kundu, David Ross, Caroline Pantofaru, Tom Funkhouser, and Justin Solomon.  Pillar-based object detection for autonomous driving.  arXiv preprint arXiv:2007.10323, 2020.
  * [31] Ziguo Zhong, Stanley Liu, Manu Mathew, and Aish Dubey.  Camera radar fusion for increased reliability in adas applications.  Electronic Imaging, 2018(17):258–1, 2018.


# Hyperspectral Image Denoising via Multi-modal and Double-weighted Tensor
Nuclear Norm

Sheng Liu, Xiaozhen Xie, and Wenfeng Kong Sheng Liu, Xiaozhen Xie, and Wenfeng
Kong are with College of Science, Northwest A&F University, Yangling 712100,
China (e-mail<EMAIL_ADDRESS>e-mail<EMAIL_ADDRESS>(Corresponding author); e-mail: kksama@nwafu.edu.cn).

###### Abstract

Hyperspectral images (HSIs) usually suffer from different types of pollution.
This severely reduces the quality of HSIs and limits the accuracy of
subsequent processing tasks. HSI denoising can be modeled as a low-rank tensor
denoising problem. Tensor nuclear norm (TNN) induced by tensor singular value
decomposition plays an important role in this problem. In this letter, we
first reconsider three inconspicuous but crucial phenomenons in TNN. In the
Fourier transform domain of HSIs, different frequency slices (FS) contain
different information; different singular values (SVs) of each FS also
represent different information. The two physical phenomenons lie not only in
the spectral mode but also in the spatial modes. Then based on them, we
propose a multi-modal and double-weighted TNN. It can adaptively shrink the FS
and SVs according to their physical meanings in all modes of HSIs. In the
framework of the alternating direction method of multipliers, we design an
effective alternating iterative strategy to optimize our proposed model.
Denoised experiments on both synthetic and real HSI datasets demonstrate their
superiority against related methods.

###### Index Terms:

Hyperspectral image, tensor nuclear norm, double weighting, frequency slices,
multi-modal.

##  I Introduction

Hyperspectral image (HSI) has been widely used in many fields [1] due to its
wealthy spatial and spectral information of a real scene. However, the
observed HSIs are usually corrupted by different noises, e.g., Gaussian noise,
impulse noise, deadlines, stripes and their mixtures. This seriously affects
the subsequent applications of HSI, such as unmixing, target detection, and so
on. Therefore, HSI denoising , as a preprocessing step to remove mixed noise
for various subsequent applications, is a valuable and active research topic.

In HSIs, different spectral bands are images from the same scene under
different wavelengths. It means the global correlation or low-rank prior lies
in HSIs [2]. Based on this prior, how to measure the HSI low-rankness becomes
the key to denoising tasks. As HSIs can be treated as 3rd-order tensor, this
problem is turned into a 3rd-order tensor low-rank problem. Due to the
nonunique definitions of the tensor rank, different tensor decompositions and
their corresponding tensor ranks are proposed, such as the Tucker
decomposition [3, 4], PARAFAC decomposition [5], and tensor singular value
decomposition (t-SVD) [6, 7], to exploit the low-rankness of HSIs.

Among them, the tensor tubal rank induced by t-SVD can characterize the low-
rank structure of HSIs very well.Its convex relaxation is the tensor nuclear
norm (TNN) [8]. TNN is effective to keep the intrinsic structure of HSIs.
Hence, TNN has attracted extensive attention for HSI denoising problems in
recent years [7, 9, 10]. However, during the definition of TNN, there are
three kinds of prior knowledge that are underutilized for further exploiting
the low-rankness in HSIs. Firstly, in the Fourier transform domain of HSIs,
the low-frequency slices carry the profile information of HSIs, while the
high-frequency slices mainly carry the noise information of HSIs. Secondly, in
each frequency slices, bigger singular values mainly contain information on
clean data and smaller singular values mainly contain information on noise.
Thirdly, low-rankness not only exists in the spectral dimension but also lies
in the spatial dimensions [3]. The classical TNN only takes the Fourier
transform to connect the spatial dimensions with the spectral dimension and
lacks flexibility for handling different correlations along with different
modes of HSIs [7].

In this letter, to take full advantage of the above prior knowledge and
improve the capability and flexibility of TNN, we propose a multi-modal and
double-weighted TNN for HSI denoising tasks. The merits of our model are four-
fold. First, according to information types in different frequency slices in
the Fourier transform domain, we adaptively assign bigger weights to slices
that mainly contain noise information and smaller weights to slices that
mainly contain profile information, which can depress noise more and
simultaneously preserve the profile information of clean HSIs better. Second,
in each frequency slice, we use the partial sum of singular values to only
shrink small singular values, which can better protect the clean data
information contained in big singular values. Third, we apply the double-
weighted TNN in all modes of HSIs, which can achieve a more flexible and
accurate characterization of HSI low-rankness. Finally, we develop an
alternating direction method of multiplier based algorithm to efficiently
solve the proposed model, Compared with various competing HSI denoising
methods, the best-denoised performances are obtained by our method synthetic
and real HSI datasets.

##  II Preliminaries

###  II-A Notations

In this letter, matrix and tensor are denoted as bold upper-case letter
$\mathbf{X}$ and calligraphic letter $\mathcal{X}$, respectively. For a 3rd-
order tensor $\mathcal{X}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$, its
$(i,j,k)$-th component is represented as $\mathcal{X}(i,j,k)$. For
$\mathcal{X}$, $\mathcal{Y}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$,
their inner product is defined as
$<\mathcal{X},\mathcal{Y}>=\sum_{i=1}^{n_{1}}\sum_{j=1}^{n_{2}}\sum_{k=1}^{n_{3}}x_{ijk}y_{ijk}$.
Then the Frobenius norm of a tensor $\mathcal{X}$ is defined as
$\|\mathcal{X}\|_{F}=\sqrt{<\mathcal{X},\mathcal{X}>}$. The $k$-th frontal
slice of $\mathcal{X}$ is denoted as $\mathbf{X}^{(k)}=\mathcal{X}(:,:,k)$.
The fast Fourier transform along the third mode of $\mathcal{X}$ is
represented as $\bar{\mathcal{X}}=\texttt{fft}(\mathcal{X},[],3)$ and its
inverse operation is $\mathcal{X}=\texttt{ifft}(\bar{\mathcal{X}},[],3)$. The
mode-$p$ permutation of $\mathcal{X}$ is defined as
$\mathcal{X}_{p}=\texttt{permute}(\mathcal{X},p)$, $p=1,2,3$, where the $m$-th
mode-3 slice of $\mathcal{X}_{p}$ is the $m$-th mode-$p$ slice of
$\mathcal{X}$, i.e.,
$\mathcal{X}(i,j,k)=$$\mathcal{X}_{1}(j,k,i)=$$\mathcal{X}_{2}(k,i,j)=$$\mathcal{X}_{3}(i,j,k)$.
Also, its inverse operation is
$\mathcal{X}=\texttt{ipermute}(\mathcal{X}_{p},p)$.

###  II-B Problem Formulation

An ideal HSI can be viewed as a 3rd-order tensor
$\mathcal{X}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$ and usually is
assumed to be low-rank. Corrupted by mixed noise, its observed version can be
modeled as

$\mathcal{Y}=\mathcal{X}+\mathcal{S}+\mathcal{N},$ (1)

where $\mathcal{Y},\mathcal{S},\mathcal{N}\in\mathbb{R}^{n_{1}\times
n_{2}\times n_{3}}$; $\mathcal{S}$ denotes the sparse noise; $\mathcal{N}$
denotes the Gaussian white noise.

HSI denoising aims to recover the ideal HSI $\mathcal{X}$ from the observed
HSI $\mathcal{Y}$ in (1). Under the framework of regularization theory, it can
briefly be formulated as

$\displaystyle\arg\min_{\mathcal{X},\mathcal{S},\mathcal{N}}$
$\displaystyle\texttt{Rank}(\mathcal{X})+\lambda\|\mathcal{S}\|_{1}+\tau\|\mathcal{N}\|_{F}^{2},$
(2) $\displaystyle s.t.$
$\displaystyle\mathcal{Y}=\mathcal{X}+\mathcal{S}+\mathcal{N},$

where $\|\cdot\|_{1}$ is $L_{1}$ norm to detect the sparse noise;
$\|\cdot\|_{F}$ describes the Gaussian noise; $\texttt{Rank}(\cdot)$
represents the rank of unknown ideal HSI; $\lambda$ and $\tau$ are non-
negative parameters.

In model (2), regularization term Rank is approximated by different
relaxations. As mentioned above, TNN is widely used convex relaxation, which
can be defined as

$\|\mathcal{X}\|_{*}:=\frac{1}{n_{3}}\sum_{k=1}^{n_{3}}\|\bar{\textbf{X}}^{(k)}\|_{*}.$
(3)

##  III The Weighted TNN

###  III-A Frequency-Weighted TNN

In (3), the frontal slice of $\bar{\mathcal{X}}$ corresponds to the frequency
component of $\mathcal{X}$ [11]. Specifically, for $\mathcal{X}$, its profile
information is contained in the low-frequency frontal slices, while its
detailed information is contained in the high-frequency ones. When
$\mathcal{X}$ is distorted by outliers, the effects on high-frequency frontal
slices are more severe. However, different frequency slices of
$\bar{\mathcal{X}}$ have the same impact on TNN in (3), which is obviously
inconsistent with the physics meaning of frequency components. Therefore, we
improve TNN in (3) by assigning different weights for different frequency
slices, and propose the frequency-weighted TNN as follows:

$\|\mathcal{X}\|_{w*}:=\sum_{k=1}^{n_{3}}w_{k}(\bar{\textbf{X}}^{(k)})\|\bar{\textbf{X}}^{(k)}\|_{*},$
(4)

where $w_{k}(\bar{\textbf{X}}^{(k)})$ is the $k$-th weight parameter. For HSI
denoising problems, the lower the frequencies are, the less the corresponding
frequency slices should be punished. By amounts of data simulations, the
weights $w_{k}$ approximatively consist with the frequencies and are inversely
proportionate to $\|\bar{\textbf{X}}^{(k)}\|_{F}$. We let

$w_{k}(\bar{\textbf{X}}^{(k)})=\frac{c_{1}}{\log(\|\bar{\textbf{X}}^{(k)}\|_{F}^{2})+\varepsilon}+c_{2},$
(5)

where $\varepsilon=10^{-10}$ is to avoid dividing by zero; $c_{1}$ is the
scaling factor after frequency normalization; $c_{2}$ is a constant. The
proposed FWTNN is different from the frequency-filtered tensor nuclear norm
(FTNN) [11]. Our weights are data dependent, but FTNN’s weights are pre-
weight.

###  III-B Double-Weighted TNN

For $\bar{\textbf{X}}^{(k)}$ in (3), the matrix nuclear norm is used as the
tightest convex surrogate for rank. However, it has limitation in the accuracy
of approximation due to its convexity. Recently, a series of improvement
methods are proposed for better approximation [12, 13]. To differently treat
singular values of $\bar{\textbf{X}}^{(k)}$, we choose partial sum of singular
values (PSSV) to only punish the smaller singular values which mainly contain
the noise information of HSIs. Then, a double-weighted TNN is proposed by
replacing the matrix nuclear norm in (4) with the PSSV of
$\bar{\textbf{X}}^{(k)}$, which is defined as

$\|\mathcal{X}\|_{dw*}:=\frac{1}{n_{3}}\sum_{k=1}^{n_{3}}w_{k}(\bar{\textbf{X}}^{(k)})\|\bar{\textbf{X}}^{(k)}\|_{\textrm{PSSV}},$
(6)

where
$\|\bar{\textbf{X}}^{(k)}\|_{\textrm{PSSV}}=\sum_{r=R+1}^{\min\\{n_{1},n_{2}\\}}\sigma_{r}(\bar{\textbf{X}}^{(k)})$;
$\sigma_{r}(\bar{\textbf{X}}^{(k)})$ is the $r$-th biggest singular value of
matrix $\bar{\textbf{X}}^{(k)}$; $R$ is a parameter indicating the number of
main singular values. The double-weighted TNN minimization problem can be
solved by following theorem.

######  Theorem 1.

Assuming that $\tau>0$, $\mathcal{X},\mathcal{Y}\in$ $\mathbb{R}^{n_{1}\times
n_{2}\times n_{3}}$, for the minimization problem

$\mathcal{X}^{*}=\arg\min_{\mathcal{X}}\tau\lVert\mathcal{X}\rVert_{dw*}+\frac{1}{2}\lVert\mathcal{X}-\mathcal{Y}\rVert_{F}^{2},$
(7)

its solution is

$\begin{array}[]{rl}\mathcal{X}^{*}=\mathcal{D}\mathcal{W}^{w,R,\tau}\left(\mathcal{Y}\right)=\texttt{ifft}(\bar{\mathcal{U}}\cdot\bar{\mathcal{S}}_{dw*}\cdot\bar{\mathcal{V}}^{T},[],3),\end{array}$
(8)

where
$\bar{\mathcal{Y}}=\bar{\mathcal{U}}\cdot\bar{\mathcal{S}}\cdot\bar{\mathcal{V}}^{T}$;
$\bar{\mathcal{S}}_{dw*}(r,r,k)=\max(\bar{\mathcal{S}}(r,r,k)-\tau
w_{r}w_{k},0);$, $R=TW(k)$
$w_{r}=\left\\{\begin{array}[]{l}\begin{matrix}0,&r\leq R\\\ \end{matrix}\\\
\begin{matrix}1,&r>R\\\ \end{matrix}\\\ \end{array}\right.$;
$w_{k}=\frac{c_{1}}{\log(\lVert\boldsymbol{\bar{\textbf{X}}}^{(k)}\rVert_{F}^{2}+\varepsilon)}+c_{2}$.

######  Lemma 1.

(PSVT [12]). Let $\tau>0$, $l=\min(m,n)$ and
$\mathbf{X},\mathbf{Y}\in\mathbb{R}^{m\times n}$ which can be decomposed by
SVD. $\mathbf{Y}$ can be considered as the sum of two matrices,
$\mathbf{Y}=\mathbf{Y}_{1}+\mathbf{Y}_{2}=\mathbf{U}_{Y1}\mathbf{D}_{Y1}\mathbf{V}_{Y1}^{\top}+\mathbf{U}_{Y2}\mathbf{D}_{Y2}\mathbf{V}_{Y2}^{\top}$
, where $\mathbf{U}_{Y1},\mathbf{V}_{Y1}$ are the singular vector matrices
corresponding to the $R$ largest singular values by S V D , and
$\mathbf{U}_{Y2}$, $\mathbf{V}_{Y2}$ from the $(R+1)$ -th to the last singular
values. Define a minimization problem for the PSSV as

$\underset{\mathbf{X}}{\arg\min}\frac{1}{2}\|\mathbf{X}-\mathbf{Y}\|_{F}^{2}+\tau\|\mathbf{X}\|_{\textrm{PSSV}}$
(9)

Then, the optimal solution of Eq.(9) can be expressed by the PSVT operator
defined as:

$\displaystyle\mathbb{P}_{R,\tau}[\mathbf{Y}]$
$\displaystyle=\mathbf{U}_{Y}\left(\mathbf{D}_{Y1}+\mathcal{S}_{\tau}\left[\mathbf{D}_{Y2}\right]\right)\mathbf{V}_{Y}^{\top}$
(10)
$\displaystyle=\mathbf{Y}_{1}+\mathbf{U}_{Y2}\mathcal{S}_{\tau}\left[\mathbf{D}_{Y2}\right]\mathbf{V}_{Y2}^{\top}$

where

$\begin{array}[]{l}\mathbf{D}_{Y1}=\operatorname{diag}\left(\sigma_{1},\cdots,\sigma_{R},0,\cdots,0\right)\\\
\mathbf{D}_{Y2}=\operatorname{diag}\left(0,\cdots,0,\sigma_{R+1},\cdots,\sigma_{l}\right)\end{array}$
(11)

and $\mathcal{S}_{\tau}[x]=\max(|x|-\tau,0)$ is the soft-thresholding operator
[14].

###### Proof.

Let $\tau>0$, $\mathcal{X},\mathcal{Y}\in$ $\mathbb{R}^{n_{1}\times
n_{2}\times n_{3}}$, $MR\in\mathbb{R}^{n_{3}}$ is multi-rank of $\mathcal{X}$,
and $c_{1},c_{2}$ are given constants. According to the definition of DWTNN
(6) and the properties of tensors in the Fourier domain, the Eq. (7) is
equivalent to

$\displaystyle\mathcal{X}^{*}=\arg\min_{\mathcal{X}}\tau\|X\|_{dw*}+\frac{1}{2}\|\mathcal{X}-\mathcal{Y}\|_{F}^{2},$
(12)
$\displaystyle=\arg\min_{\mathcal{X}}\tau\frac{1}{n_{3}}\sum_{k=1}^{n_{3}}w_{k}\left(\overline{\mathbf{X}}^{(k)}\right)\left\|\overline{\mathbf{X}}^{(k)}\right\|_{\textrm{PSSV}}+\frac{1}{2n_{3}}\|\overline{\mathcal{X}}-\overline{\mathcal{Y}}\|_{F}^{2}$
$\displaystyle=\arg\min_{\mathcal{X}}\frac{1}{n_{3}}\sum_{k=1}^{n_{3}}\tau
w_{k}\left(\overline{\mathbf{X}}^{(k)}\right)\left\|\overline{\mathbf{X}}^{(k)}\right\|_{\textrm{PSSV}}+\frac{1}{2}\left\|\overline{\mathbf{X}}^{(k)}-\overline{\mathbf{Y}}^{(k)}\right\|_{F}^{2}.$

Therefore, the Eq. (12) can be divided into $n_{3}$ subproblems as follows:

$\arg\min_{\mathcal{X}}\tau
w_{k}\left(\overline{\mathbf{X}}^{(k)}\right)\|\overline{\mathbf{X}}\|_{p=R}+\frac{1}{2}\|\overline{\mathbf{X}}-\overline{\mathbf{Y}}\|_{F}^{2},$
(13)

the Eq.(13) can be solved by Lemma 1. Let

$w_{r}=\left\\{\begin{array}[]{ll}0,&i\leqslant R\\\ 1,&i>R\end{array}\right.$
(14)

Base on Eq.(14), Eq.(10) in Lemma 1 is equivalent to

$\displaystyle\mathbb{P}_{R,\tau,w_{k}}\left[\overline{\mathbf{Y}}\right]=\overline{\mathbf{U}}_{Y}\left(\mathcal{S}_{R,\tau,w_{k}}\left[\overline{\mathbf{D}}\right]\right)\overline{\mathbf{V}}_{Y}^{\top}$
(15)

where

$\displaystyle\mathcal{S}_{R,\tau,w_{k}}\left[\overline{\mathbf{D}}_{i}\right]$
$\displaystyle=\begin{cases}\sigma_{i}\left(\overline{\mathbf{Y}}\right)&\text{if}\
i<R+1\\\ \max\left(\sigma_{i}\left(\overline{\mathbf{Y}}\right)-\tau
w_{k},0\right)&\text{otherwise}\\\ \end{cases}$ (16)
$\displaystyle=\max\left(\sigma_{i}\left(\overline{\mathbf{Y}}\right)-\tau
w_{k}w_{r},0\right)$

Therefore, the solution of Eq. (7) is

$\begin{array}[]{rl}\mathcal{X}^{*}=\mathcal{D}\mathcal{W}^{w,R,\tau}\left(\mathcal{Y}\right)=\texttt{ifft}(\bar{\mathcal{U}}\cdot\bar{\mathcal{S}}_{dw*}\cdot\bar{\mathcal{V}}^{T},[],3),\end{array}$
(17)

where
$\bar{\mathcal{Y}}=\bar{\mathcal{U}}\cdot\bar{\mathcal{S}}\cdot\bar{\mathcal{V}}^{T}$;
$\bar{\mathcal{S}}_{dw*}(r,r,k)=\max(\bar{\mathcal{S}}(r,r,k)-\tau
w_{r}w_{k},0);$, $R=TW(k)$
$w_{r}=\left\\{\begin{array}[]{l}\begin{matrix}0,&r\leq R\\\ \end{matrix}\\\
\begin{matrix}1,&r>R\\\ \end{matrix}\\\ \end{array}\right.$;
$w_{k}=\frac{c_{1}}{\log(\lVert\boldsymbol{\bar{\textbf{X}}}^{(k)}\rVert_{F}^{2}+\varepsilon)}+c_{2}$.
∎

######  Remark 1.

In theorem 1, there are two key weights that need to be given in advance.
Frequency weights $w_{k}$ can be calculated by Eq. (5). Truncation weight
$w_{r}$ can be calculated by Algorithm 1. It refers to the way multi-rank is
calculated in reference [15]. In Algorithm 1, the ratio of the maximum value
of the singular values is chosen, which can also be replaced by the ratio of
the sum of the singular values.

Algorithm 1 Calculate truncation weight $TW\in\mathbb{R}^{n_{3}}$.

Input: A tensor $\mathcal{Y}$; ratio $\eta$.
Output:Truncation weight $TW\in\mathbb{R}^{n_{3}}$.
1: Compute $\overline{\mathcal{Y}}=\operatorname{fft}(\boldsymbol{Y},[],3)$ .
2: $\text{for}\ i=1,\cdots,n_{3}\ \text{do}$
$[\overline{\mathbf{U}}^{(i)},\overline{\mathbf{S}}^{(i)},\overline{\mathbf{V}}^{(i)}]$
$s=\text{diag}(\overline{\mathbf{S}}^{(i)})$
$TW(i)=\text{length}(s(s>max(s)\times\eta))$.
end for

###  III-C Multi-modal and Double-Weighted TNN

TNN in (3) only approximates the correlations connected by mode-3 Fourier
transform in the spatial dimensions with the spectral dimension. It lacks of
flexibility for describing low-rankness in all modes of HSIs. To connect the
$p$-th mode with other two modes, we can define the double-weighted TNN for
each mode-$p$ permutation of HSIs, i.e., $\|\mathcal{X}_{p}\|_{dw*},p=1,2,3.$
As $\|\mathcal{X}_{p}\|_{dw*}$ are different according to different modes, we
use the weighted average of double-weighted TNNs along all modes to
approximate the tensor rank of HSIs. Finally, the multi-modal and double-
weighted TNN (MDWTNN) is proposed as follows:

$\displaystyle\|\mathcal{X}\|_{mdw*}:=\sum_{p=1}^{3}\alpha_{p}\|\mathcal{X}_{p}\|_{dw*}=\sum_{p=1}^{3}\sum_{k=1}^{n_{p}}\alpha_{p}w_{k}^{p}\|\bar{\textbf{X}}^{(k)}_{p}\|_{\textrm{PSSV}},$
(18)

where $\bar{\mathcal{X}}_{p}=\texttt{fft}(\mathcal{X}_{p},[],3)$;
$\bar{\textbf{X}}^{(k)}_{p}$ is the $k$-th frontal slice of
$\bar{\mathcal{X}}_{p}$ and its assigned weight is $w_{k}^{p}$; $\alpha_{p}>0$
and $\sum_{p=1}^{3}\alpha_{p}=1$. Fig. 1 shows the schematic diagram of
MDWTNN.

Figure 1:  MDWTNN Schematic diagram. (I) Multi-modal permutations and Fourier
transforms. (II) Double-weighted TNN. (a) mode-3 Fourier transform of
$\mathcal{X}_{p},p=1,2,3$. Due to the repeatability of the operation, II only
represents one of $\mathcal{X}_{1}$, $\mathcal{X}_{2}$, $\mathcal{X}_{3}$, and
the other two only need to repeat the operation represented by II. (b)
different frequency components. (c) PSSV relaxation on each frequency slice
$\\{\bar{\mathcal{X}}_{i}\\}$. (d) The singular values of frequency slices.
(e) Frequency weights. The top is the sum of singular values in each frequency
slice, and the bottom is the weight of each frequency slice. (f) Doubel-
weighted TNN minimization with respect to $\mathcal{X}_{p}$. (III) Synthesis
from multi-modal denoised results.

##  IV HSI denoising via MDWTNN Minimization

MDWTNN in (18) takes full advantage of physical meanings in frequency
components, singular values, and modes of HSIs, which can provide a better
approximation to the tensor rank. Then we use MDWTNN to replace the
regularization term Rank in (2) and propose the HSI denoising model as
follows:

$\begin{array}[]{rl}\arg\min_{\mathcal{X},\mathcal{S},\mathcal{N}}&\|\mathcal{X}\|_{mdw*}+\lambda\|\mathcal{S}\|_{1}+\tau\|\mathcal{N}\|_{F}^{2},\\\
s.t.&\mathcal{Y}=\mathcal{X}+\mathcal{S}+\mathcal{N}.\end{array}$ (19)

Introducing auxiliary variables, model (19) is equivalent to

$\begin{array}[]{rl}\displaystyle\arg\min_{\mathcal{X},\mathcal{S},\mathcal{N}}\sum_{p=1}^{3}\alpha_{p}\|\mathcal{Z}_{p}\|_{dw*}+\lambda\lVert\mathcal{S}\rVert_{1}+\tau\lVert\mathcal{N}\rVert_{F}^{2},\\\
\displaystyle s.t.\ \mathcal{Y}=\mathcal{X}+\mathcal{S}+\mathcal{N},\
\mathcal{Z}_{p}=\mathcal{X}_{p},\ p=1,2,3.\\\ \end{array}$ (20)

By augmented Lagrangian multiplier method, the Lagrangian function of model
(20) can be written as

$\begin{array}[]{l}L_{\mu_{p},\beta}\left(\mathcal{X},\mathcal{Z}_{p},\mathcal{N},\mathcal{S},\Gamma_{p},\Lambda\right)=\lambda\lVert\mathcal{S}\rVert_{1}+\tau\lVert\mathcal{N}\rVert_{F}^{2}\\\
+<\mathcal{Y}-\left(\mathcal{X}+\mathcal{S}+\mathcal{N}\right),\Lambda>+\frac{\beta}{2}\lVert\mathcal{Y}-\left(\mathcal{X}+\mathcal{S}+\mathcal{N}\right)\rVert_{F}^{2},\\\
\displaystyle+\sum_{p=1}^{3}{\left\\{\alpha_{p}\lVert\mathcal{X}_{p}\rVert_{dw*}+<\mathcal{X}_{p}-\mathcal{Z}_{p},\Gamma_{p}>+\frac{\mu_{p}}{2}\lVert\mathcal{X}_{p}-\mathcal{Z}_{p}\rVert_{F}^{2}\right\\}},\\\
\end{array}$

where $\Lambda$ and $\Gamma_{p}$ are the Lagrangian multipliers; $\beta$ and
$\mu_{p}$ are the Lagrange penalty parameters. Its minimization problem can be
efficiently solved in the framework of ADMM [16]. At the $(n+1)$-th iteration,
each variable in the Lagrangian function can be updated by solving its
corresponding subproblem respectively when other variables are fixed at the
$n$-th iteration.

For $\mathcal{Z}_{p}$, $p=1,2,3$, their corresponding subproblems can be
written as

$\arg\min_{\mathcal{Z}_{p}}\alpha_{p}\lVert\mathcal{Z}_{p}\rVert_{dw*}+\frac{\mu_{p}}{2}\left\|\mathcal{Z}_{p}-\left(\mathcal{X}_{p}^{n}+\frac{\Gamma_{p}^{n}}{\mu_{p}}\right)\right\|_{F}^{2}.$
(21)

The closed-form solution of (21) obtained from theorem 1 are as follows:

$\begin{array}[]{rl}\mathcal{Z}_{p}^{n+1}=\mathcal{D}\mathcal{W}^{w\left(\mathcal{X}_{p}^{n}\right),R,\frac{\alpha_{p}}{\mu_{p}}}\left(\mathcal{X}_{p}^{n}+\frac{\Gamma_{p}^{n}}{\mu_{p}}\right).\\\
\end{array}$ (22)

For $\mathcal{X}$, its corresponding subproblem can be written as

$\displaystyle\mathcal{X}^{n+1}=$
$\displaystyle\arg\min_{\mathcal{X}}\sum_{p=1}^{3}{\frac{\mu_{p}}{2}}\lVert\mathcal{X}-\mathcal{Z}_{p}^{n+1}+\frac{\Gamma_{p}^{n}}{\mu_{p}}\rVert_{F}^{2}$
(23)
$\displaystyle+\frac{\beta}{2}\lVert\mathcal{Y}-\left(\mathcal{X}+\mathcal{S}^{n}+\mathcal{N}^{n}\right)+\frac{\Lambda^{n}}{\beta}\rVert_{F}^{2}.$

It has the closed-form solution as follows:

$\begin{array}[]{rl}\mathcal{X}^{n+1}=\frac{\sum_{p=1}^{3}{\mu_{p}}\left(\mathcal{Z}_{p}^{n+1}-\frac{\Gamma_{p}^{n}}{\mu_{p}}\right)+\beta\left(\mathcal{Y}-\mathcal{S}^{n}-\mathcal{N}^{n}+\frac{\Lambda^{n}}{\beta}\right)}{1+\beta}.\\\
\end{array}$ (24)

For $\mathcal{S}$, its corresponding subproblem can be written as

$\arg\min_{\mathcal{S}}\lambda\lVert\mathcal{S}\rVert_{1}+\frac{\beta}{2}\lVert\mathcal{Y}-\left(\mathcal{X}^{n+1}+\mathcal{S}+\mathcal{N}^{n}\right)+\frac{\Lambda^{n}}{\beta}\rVert_{F}^{2}.$
(25)

It can be solved by the soft-thresholding operator [14] as:

$\begin{array}[]{rl}\mathcal{S}^{n+1}=\texttt{shrink}\left(\mathcal{Y}-\mathcal{X}^{n+1}-\mathcal{N}^{n+1}+\frac{\Lambda^{n}}{\beta},\frac{\lambda}{\beta}\right).\end{array}$
(26)

For $\mathcal{N}$, its corresponding subproblem can be written as

$\arg\min_{\mathcal{N}}\tau\|\mathcal{N}\|_{F}^{2}+\frac{\beta}{2}\|\mathcal{Y}-\left(\mathcal{X}^{n+1}+\mathcal{S}^{n+1}+\mathcal{N}\right)+\frac{\Lambda^{n}}{\beta}\|_{F}^{2}.$
(27)

It has the closed-form solution as follows :

$\begin{array}[]{rl}\mathcal{N}^{n+1}=\frac{\beta\left(\mathcal{Y}-\mathcal{X}^{n+1}-\mathcal{S}^{n}+\frac{\Lambda^{n}}{\beta}\right)}{2\tau+\beta}.\end{array}$
(28)

For multipliers $\Gamma_{p}$ and $\Lambda$, they can be updated as follows:

$\left\\{\begin{array}[]{l}\Gamma_{p}^{n+1}=\Gamma_{p}^{n}+\mu_{p}\left(\mathcal{Z}_{p}^{n+1}-\mathcal{X}^{n+1}\right),p=1,2,3\\\
\Lambda^{n+1}=\Lambda^{n}+\beta\left(\mathcal{Y}-\mathcal{X}^{n+1}-\mathcal{S}^{n+1}-\mathcal{N}^{n+1}\right).\\\
\end{array}\right.$ (29)

The proposed algorithm for our HSI denoising model is summarized in Algorithm
2.

Algorithm 2 HSI denoising via the MDWTNN minimization

Input: The observed tensor $\mathcal{Y}$; weight parameters $c_{1}$, $c_{2}$,
ratio $\eta$; regularization parameters $\lambda$, $\tau$; and stopping
criterion $\epsilon$.
Output: Denoised image $\mathcal{X}$.
1: Initialize:
$\mathcal{Y}$=$\mathcal{X}$=$\mathcal{S}$=$\mathcal{N}$=$\mathcal{Z}_{p}$;
$\Gamma_{p}=\Lambda=0$; $\mu_{p}$=$\beta$=$10^{-3}$;
$p=1,2,3$; $\mu_{max}=10^{10}$; $\rho=1.2$ and $n=0$.
2: Repeat until convergence:
3\. Update
$\mathcal{Z}_{p}$,$\mathcal{X}$,$\mathcal{S}$,$\mathcal{N}$,$\mathcal{X}$
$\Gamma_{p},\Lambda$ by (22), (24), (26), (28), (29)
Update $\mu_{p}=\rho\mu_{p}$, $\beta=\rho\beta$, $w_{k}$ by (5)
4: Check the convergence condition.

##  V Experiments

To verify the effectiveness of our MDWTNN based HSI denoising model, various
experiments were performed on two challenging simulated datasets and two real
HSI datasets. For comparison, four state-of-the-art HSI denoising methodes
were employed as the benchmark in the experiments, i.e., BM4D [17], LRMR [18],
LRTDTV [3] and 3DTNN [7]. Since the BM4D method was only suitable to remove
Gaussian noise, we implemented it on HSIs which were preprocessed by the RPCA
denoising method [6].

###  V-A Simulated Data Experiments

In the simulation experiments, we selected two datasets. From the Washington
DC Mall dataset111http://lesun.weebly.com/hyperspectral-data-set.html, we
chose a sub-block with the size of $256\times 256\times 191$ as a simulation
dataset. From the Pavia City Center
dataset222http://www.ehu.eus/ccwintco/index.php/, we chose a sub-block with
the size of $200\times 200\times 80$ as a simulation dataset. The hybrids of
white Gaussian and impulse noises with 5 different intensity levels were added
to the simulation datasets band by band. Let $G$ and $P$ denoted the variance
of Gaussian white noise and percentage of impulse noise, respectively. In
noise case 1-3, the same intensity noise was added to all the bands. In noise
case 1, $G$=0.1 and $P$=0.2; In noise case 2, $G$=0.2 and $P$=0.2; In noise
case 3, $G$=0.1 and $P$=0.4; In noise case 4 and 5, the noise intensities were
different for different bands. In noise case 4, $G$ was randomly selected from
0.1 to 0.2 and $P$=0.2; In noise case 5, $G$=0.1 and $P$ was randomly selected
from 0.2 to 0.4.

For quantitatively evaluating the denoised results of all the test methods,
the CPU times and the means of PSNR, SSIM and SAM in each band, i.e., MPSNR,
MSSIM and MSAM, were listed in Table I. Although the CPU times of our model
were not the shortest, one could update $\mathcal{Z}_{p}$ by (22) in parallel
to further shorten the CPU times of our model. For visual evaluation in Fig.
2, we showed the denoised results of the Pavia City Center dataset in Case 1.

TABLE I: Quantitative comparison and time of all competing methods under
different levels of noises on simulate dataset.

Dataset | Noise case | Index | Noise | BM4D | LRMR | LRTDTV | 3DTNN | Our
---|---|---|---|---|---|---|---|---
Washington DC Mall | Case1 | PSNR | 11.068 | 31.014 | 31.567 | 32.800 | 34.270 | 36.095
SSIM | 0.085 | 0.893 | 0.867 | 0.896 | 0.936 | 0.950
MSAM | 43.139 | 4.576 | 5.042 | 4.327 | 3.481 | 2.937
time | - | 547.005 | 378.898 | 538.160 | 270.211 | 334.656
Case2 | PSNR | 10.216 | 27.192 | 27.642 | 29.489 | 29.055 | 32.525
SSIM | 0.061 | 0.791 | 0.743 | 0.807 | 0.801 | 0.894
MSAM | 45.297 | 6.830 | 7.859 | 6.382 | 6.726 | 4.393
time | - | 529.243 | 395.461 | 584.461 | 309.977 | 378.728
Case3 | PSNR | 8.305 | 29.691 | 29.183 | 30.270 | 29.572 | 34.185
SSIM | 0.037 | 0.866 | 0.798 | 0.844 | 0.784 | 0.928
MSAM | 50.092 | 5.264 | 6.579 | 6.182 | 6.542 | 3.604
time | - | 528.440 | 394.272 | 582.318 | 316.547 | 375.719
Case4 | PSNR | 10.648 | 28.970 | 29.518 | 31.073 | 31.852 | 34.379
SSIM | 0.073 | 0.848 | 0.810 | 0.859 | 0.887 | 0.930
MSAM | 44.265 | 5.695 | 6.439 | 5.468 | 4.873 | 3.579
time | - | 541.960 | 396.478 | 537.374 | 273.466 | 340.113
Case5 | PSNR | 9.669 | 30.419 | 30.412 | 31.560 | 32.783 | 35.180
SSIM | 0.060 | 0.883 | 0.836 | 0.874 | 0.902 | 0.942
MSAM | 47.270 | 4.865 | 5.761 | 5.405 | 4.311 | 3.220
time | - | 546.378 | 397.474 | 541.124 | 277.123 | 340.162
Pavia City Center | Case1 | MPSNR | 11.122 | 29.701 | 31.259 | 32.297 | 31.696 | 33.951
MSSIM | 0.105 | 0.920 | 0.905 | 0.914 | 0.924 | 0.942
MSAM | 45.712 | 5.84 | 6.824 | 4.93 | 4.844 | 4.426
time | - | 112.235 | 113.723 | 131.91 | 54.461 | 72.224
Case2 | MPSNR | 10.265 | 25.619 | 27.321 | 28.605 | 27.009 | 30.124
MSSIM | 0.074 | 0.835 | 0.791 | 0.821 | 0.799 | 0.879
MSAM | 46.978 | 7.170 | 8.385 | 6.685 | 6.923 | 5.666
time | - | 113.246 | 113.115 | 129.155 | 56.610 | 71.322
Case3 | MPSNR | 8.384 | 27.745 | 28.937 | 29.741 | 27.696 | 31.924
MSSIM | 0.043 | 0.892 | 0.848 | 0.873 | 0.790 | 0.916
MSAM | 48.580 | 6.745 | 7.74 | 6.007 | 9.875 | 5.095
time | - | 115.908 | 116.961 | 127.981 | 54.448 | 71.435
Case4 | MPSNR | 10.659 | 27.361 | 29.094 | 30.251 | 28.934 | 32.041
MSSIM | 0.088 | 0.877 | 0.853 | 0.870 | 0.860 | 0.915
MSAM | 46.457 | 6.6 | 7.796 | 5.877 | 6.389 | 5.048
time | - | 115.684 | 133.241 | 129.088 | 51.669 | 68.251
Case5 | MPSNR | 9.565 | 28.674 | 30.044 | 30.977 | 30.055 | 33.038
MSSIM | 0.070 | 0.907 | 0.878 | 0.893 | 0.879 | 0.932
MSAM | 47.885 | 6.293 | 7.333 | 5.493 | 6.795 | 4.682
time | - | 133.811 | 81.989 | 120.418 | 61.742 | 77.520

(a) Original image (b) Noise image (c) BM4D(24.53dB) (d) LRMR(27.84dB)

(e)  LRTDTV (27.13dB) (f) 3DTNN(27.55dB) (g) Our(28.08dB)

Figure 2: The denoised results of the Pavia City Center dataset in Case 1.
Pseudocolor image with bands (78, 58, 14).

###  V-B Real Data Experiments

In the real experiments, we selected two datasets. From the Indian Pines
dataset333https://engineering.purdue.edu/∼biehl/MultiSpec/hyperspectral, we
chose a sub-block with the size of $145\times 145\times 224$. From the
Australian dataset444http://remote-
sensing.nci.org.au/u39/public/html/index.shtml, we chose a sub-block with the
size of $200\times 200\times 150$. In Fig. 3, we listed all denoised results
of the Indian Pines and Australian datasets. For BM4D and LRTDTV, although
they could remove more noise, they also lost more details. This made the
denoised result too smooth. For LRMR and 3DTNN, they mainly used the low-rank
information of HSI. Although they retained more details, they also retained
more noise. Compared with them, our proposed model could remove more noise
while retaining more details. Fig. 4 showed that the vertical mean profiles of
band 218 before and after denoising. Here, one could see that the curve of the
proposed MDWTNN method was most stable.

(a) Real band (b) BM4D (c) LRMR (d) LRTDTV

(e) 3DTNN (f) Our

(g) Real band (h) BM4D (i) LRMR (j) LRTDTV

(k) 3DTNN (l) Our

Figure 3:  All denoised results for the Indian Pines and Australian dataset.
(a)-(f) are the denoised results of the 150 band of the Indian Pines dataset.
(g)-(l) are the denoised results of the 48 band of the Australian dataset.

(a) Real dataset (b) BM4D (c) LRMR (d) LRTDTV

(e) 3DTNN (f) Our

Figure 4:  Vertical mean profiles of band 218 by all denoised results for the
Indian Pines.

##  VI Conclusion

In this letter, we propose a multi-modal and double-weighted TNN for HSI
denoising tasks. The proposed TNN can efficiently characterize the physical
meanings of the frequency components, singular values, and orientations
ignored by the standard TNN. And the weight parameters also can be obtained
adaptively. They powerfully improve capability and flexibility for describing
low-rankness in HSIs. The experiments conducted with both simulate and real
HSI datasets show that our MDWTNN based HSI denoising model is a competitive
method to remove the hybrid noise. Besides, our proposed MDWTNN regularization
term can also be applied to other low-rankness based tasks, i.e.,
hyperspectral imagery classification, tensor completion, MRI reconstruction.

## References

  * [1] J. M. Bioucas-Dias, A. Plaza, G. Camps-Valls, P. Scheunders, N. Nasrabadi, and J. Chanussot, “Hyperspectral remote sensing data analysis and future challenges,” _IEEE Geoscience and Remote Sensing Magazine_ , vol. 1, no. 2, pp. 6–36, 2013.
  * [2] F. Xu, Y. Chen, C. Peng, Y. Wang, X. Liu, and G. He, “Denoising of hyperspectral image using low-rank matrix factorization,” _IEEE Geoscience and Remote Sensing Letters_ , vol. 14, no. 7, pp. 1141–1145, 2017.
  * [3] Y. Wang, J. Peng, Q. Zhao, Y. Leung, X.-L. Zhao, and D. Meng, “Hyperspectral image restoration via total variation regularized low-rank tensor decomposition,” _IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing_ , vol. 11, no. 4, pp. 1227–1243, 2017.
  * [4] S. Meng, L.-T. Huang, and W.-Q. Wang, “Tensor decomposition and pca jointed algorithm for hyperspectral image denoising,” _IEEE Geoscience and Remote Sensing Letters_ , vol. 13, no. 7, pp. 897–901, 2016.
  * [5] X. Liu, S. Bourennane, and C. Fossati, “Denoising of hyperspectral images using the parafac model and statistical performance analysis,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 50, no. 10, pp. 3717–3724, 2012.
  * [6] C. Lu, J. Feng, Y. Chen, W. Liu, Z. Lin, and S. Yan, “Tensor robust principal component analysis with a new tensor nuclear norm,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 42, no. 4, pp. 925–938, 2019\.
  * [7] Y.-B. Zheng, T.-Z. Huang, X.-L. Zhao, T.-X. Jiang, T.-H. Ma, and T.-Y. Ji, “Mixed noise removal in hyperspectral image via low-fibered-rank regularization,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 58, no. 1, pp. 734–749, 2019.
  * [8] M. E. Kilmer, K. Braman, N. Hao, and R. C. Hoover, “Third-order tensors as operators on matrices: A theoretical and computational framework with applications in imaging,” _SIAM Journal on Matrix Analysis and Applications_ , vol. 34, no. 1, pp. 148–172, 2013.
  * [9] H. Zeng, X. Xie, H. Cui, H. Yin, and J. Ning, “Hyperspectral image restoration via global $l_{1-2}$ spatial-spectral total variation regularized local low-rank tensor recovery,” _IEEE Transactions on Geoscience and Remote Sensing_ , 2020.
  * [10] H. Zeng, X. Xie, and J. Ning, “Hyperspectral image denoising via global spatial-spectral total variation regularized nonconvex local low-rank tensor approximation,” _Signal Processing_ , vol. 178, p. 107805, 2021.
  * [11] S. Wang, Y. Liu, L. Feng, and C. Zhu, “Frequency-weighted robust tensor principal component analysis,” _arXiv preprint arXiv:2004.10068_ , 2020.
  * [12] T.-H. Oh, Y.-W. Tai, J.-C. Bazin, H. Kim, and I. S. Kweon, “Partial sum minimization of singular values in robust pca: Algorithm and applications,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 38, no. 4, pp. 744–758, 2015.
  * [13] S. Gu, L. Zhang, W. Zuo, and X. Feng, “Weighted nuclear norm minimization with application to image denoising,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2014, pp. 2862–2869.
  * [14] E. T. Hale, W. Yin, and Y. Zhang, “Fixed-point continuation for $\backslash$ell_1-minimization: Methodology and convergence,” _SIAM Journal on Optimization_ , vol. 19, no. 3, pp. 1107–1130, 2008.
  * [15] T.-X. Jiang, T.-Z. Huang, X.-L. Zhao, and L.-J. Deng, “Multi-dimensional imaging data recovery via minimizing the partial sum of tubal nuclear norm,” _Journal of Computational and Applied Mathematics_ , vol. 372, p. 112680, 2020\.
  * [16] S. Boyd, N. Parikh, and E. Chu, _Distributed optimization and statistical learning via the alternating direction method of multipliers_. Now Publishers Inc, 2011.
  * [17] M. Maggioni and A. Foi, “Nonlocal transform-domain denoising of volumetric data with groupwise adaptive variance estimation,” in _Computational Imaging X_ , vol. 8296. International Society for Optics and Photonics, 2012, p. 82960O.
  * [18] H. Zhang, W. He, L. Zhang, H. Shen, and Q. Yuan, “Hyperspectral image restoration using low-rank matrix recovery,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 52, no. 8, pp. 4729–4743, 2013.