1.   引言

海表面温度对于研究海洋动力过程、海气相互作用、海洋资源开发利用、地球气候变化检测等问题^[1,2]具有重要意义。浮标和卫星遥感是获取海表面温度数据的重要方法。但是，浮标在进行海表面温度监测时，由于数量有限，难以覆盖整个海洋区域，且由于传感器故障等原因，使得观测海域的海表面温度数据难免存在数据缺失问题；此外，由于海洋浮标放置的位置受到海水流动影响，导致收集到的浮标数据呈现空间上的不规则分布；而利用卫星遥感数据观测海表面温度时，会受到云层、降水和气溶胶等因素的影响，导致获取的数据呈现空间上的不规则分布以及一定的缺失。因此，研究存在缺失的不规则海表面温度数据重构问题具有重要应用价值。现有对于此类不完整数据进行重构的方法，可以分为基于统计的经验模型^[3]和插值方法^[4]两类。对于前者，需要特定区域大量的历史数据.，而对于缺少历史数据的区域，重构精度有待提高；对于插值方法，存在缺失的测量数据不足以确定插值方法的回归模型，如果预先选择的基函数不能反映数据之间的空间相关性，则会导致重构性能退化。上述方法都没有利用海表面温度数据的时空特性，难以有效处理非规则海表面温度数据重构问题。本文引入图信号处理方法^[5]来重构这类不规则海洋观测数据，将海表面温度数据建模为时变图信号，通过图的结构来捕获时空数据的关联关系，建立时变图信号重构模型用于缺失数据重构，从而提高重构的准确性。

当前，已提出了多种时变图信号重构算法。为降低算法的计算复杂度，文献[6]提出了一种采用截断Jacobi迭代的广义牛顿法求解图信号矩阵补全问题；文献[7]提出一种基于笛卡尔乘积图上Sobolev平滑的分布式批量重构算法，有效提高大规模数据重构效率。为降低重构误差，学者们提出了一系列图信号重构算法，该类算法主要是基于数据的相关性，然后采用不同的先验信息进行重构。文献[8]提出了基于1阶时差信号平滑性时变图信号重构算法，显著降低了重构误差。在文献[8]的基础上，文献[9]加入了低秩先验信息，提出了基于低秩和差分平滑的重构算法(Low Rank and Differential Smoothness, LRDS) ，从而降低了重构误差；文献[10]为降低重构误差和提高重构效率，提出了基于Sobolev差分平滑的重构算法(Graph Time-varying-signal Reconstruction via Sobolev Smoothness, GraphTRSS)；文献[11]考虑了时变图信号的时域-图域联合平滑性，提出基于联合平滑的重构算法 (Joint Smooth Reconstruction of Time-Varying Graph Signals, JSR-TVGS)，进一步提高了重构精度。较文献[9]而言，文献[12]考虑信号的空间结构和时间相关性，建立了一个1阶马尔可夫模型用于时变图信号重构，降低了重构误差；文献[13]基于2阶时差信号平滑与低秩性提出了重构算法，提高了重构精度。此外，文献[14–16]从图信号处理的角度重建海洋数据，结果表明利用图信号处理的方法较已有方法而言更具有优越性。

受上述相关研究的启发，为了以较低误差重构非规则海表面温度，本文以图信号处理为基础，引入2阶时差信号，提出了基于低秩和联合平滑性(Low Rank and Joint Smoothness, LRJS)的时变图信号重构模型，采用交替方向乘子法求解，并进一步分析了所提算法的计算复杂度和重构误差的极限。最后在两个公开的真实数据集上进行实验，实验结果表明，本文所提方法较已有的方法而言，重构误差更小。

2.   图信号处理基础理论

在图信号处理^[5]中，图可以表示为由节点和边组成的图结构，其中节点表示数据样本，边表示节点之间的关联关系。本文将海表面温度数据抽象为时变图信号，建立一个无向图$G = (V,E,{\boldsymbol{W}})$来表示海域面积，其中$V = \{ {v_1},{v_{2,}} \cdots ,{v_{{N}}}\}$是顶点集，$N$为顶点数；$E = \{ {e_{ij}}\}$是边集，${e_{ij}}$表示顶点${v_i}$和${v_j}$有边相连；${\boldsymbol{W}} \in {R^{N \times N}}$是图的权矩阵，${{\boldsymbol{W}}_{ij}}$表示顶点${v_i}$和${v_j}$边的权重。若海表面温度数据的采样时刻总数为$T$，对于图中第$i$个顶点${v_i}$，假设该点在$t$时刻的图信号为${x_{it}}$，那么整个图上的$t$时刻的静态图信号可以表示为${{\boldsymbol{x}}_{{t}}} = {[{{{x}}_{1{{t}}}},{{{x}}_{2{{t}}}}, \cdots ,{{{x}}_{{{Nt}}}}]^{{\mathrm{T}}} } \in {{{R}}^{{N}}}$，对于时变图信号${\boldsymbol{X}}$，可以将其表示为形如${\boldsymbol{X}} = [{{\boldsymbol{x}}_1},{{\boldsymbol{x}}_2}, \cdots ,{{\boldsymbol{x}}_T}] \in {R^{N \times T}}$的矩阵形式。

图拉普拉斯矩阵^[5]定义为${{\boldsymbol{L}}_{\rm{G}}}{\boldsymbol{ = D - W}}$，式中矩阵${\boldsymbol{D}} = {{\mathrm{diag}}} ([{d_1},{d_2}, \cdots ,{d_{{N}}}])$表示图$G$的度矩阵；拉普拉斯矩阵${{\boldsymbol{L}}_{\rm{G}}}$满足${{\boldsymbol{L}}_{\rm{G}}}{\boldsymbol{ = U}}{{\boldsymbol{\varLambda }}_{\rm{G}}}{{\boldsymbol{U}}^{\text{T}}}$，其中${\boldsymbol{U}}$为${{\boldsymbol{L}}_{\rm{G}}}$的正交特征向量矩阵，${{\boldsymbol{\varLambda }}_{\rm{G}}}$为${{\boldsymbol{L}}_{\rm{G}}}$的非负特征值对角矩阵。图傅里叶变换^[5]定义为$\widehat {\boldsymbol{X}}{\boldsymbol{ = }}{{\boldsymbol{U}}^{\text{T}}}{\boldsymbol{X}}$ 。

3.   基于低秩和联合平滑性的时变图信号重构方法

3.1   联合平滑性

海表面温度数据存在两种相关性。一是空间上相邻节点的图信号存在相关性；二是时间上相邻采样时刻获得的图信号存在相关性。这种相关性可通过信号的变差，或者平滑性来体现。

常用图拉普拉斯2次形^[5]来表示$t$时刻的图信号${{\boldsymbol{x}}_{{t}}} \in {{{R}}^{{N}}}$的全局平滑性

上述2次形值越小，表示图上顶点与其领域内顶点上的信号值越相似，即图信号越平滑。

对于时变图信号矩阵${\boldsymbol{X}} \in {R^{N \times T}}$，其平滑性可定义为

其中${{\mathrm{tr}}} ( \cdot )$表示矩阵的迹。

已有研究表明，时差信号比信号本身更具有平滑性^[8]。Liu等人^[13]基于Hodrick-Prescott滤波，引入了2阶时差算子${{\boldsymbol{D}}_{{\mathrm{h}}} }$来捕获数据在连续时间上的变化特性，具体形式为

那么，基于2阶时差算子的2阶时差信号${\boldsymbol{X}}{{\boldsymbol{D}}_{\mathrm{h}}}$可定义为

基于上述分析，2阶时差信号${\boldsymbol{X}}{{\boldsymbol{D}}_{\mathrm{h}}}$在图域上的平滑性可定义为

另外，信号的傅里叶变换可以定量反映该信号变化快慢特性。对时变图信号${\boldsymbol{X}}$应用图傅里叶变换(Graph Fourier Transform, GFT)，即${{\mathrm{GFT}}} ({\boldsymbol{X}}) = {{\boldsymbol{U}}^{\text{T}}}{\boldsymbol{X}}$，可以分析其图域上的变化特性，该操作相当于对矩阵${\boldsymbol{X}}$每列独立地进行图傅里叶变换，但忽视了信号的时间维度。由于每个节点上的海洋温度值都是时变信号，所以对${\boldsymbol{X}}$的每一行应用离散傅里叶变换(Discrete Fourier Transform, DFT)^[17,18]，可分析其时间上的变化特性，即${{\mathrm{DFT}}} ({\boldsymbol{X}}) = {\boldsymbol{X}}{{\boldsymbol{V}}^*}$，其中${\boldsymbol{V}}$为归一化DFT矩阵，${( \cdot )^*}$表示复共轭矩阵。

为了沿时间域和图域同时捕获时变图信号${\boldsymbol{X}}$的变化特性，对${\boldsymbol{X}}$同时应用GFT和DFT，称为时域-图域联合图傅里叶变换(Joint Fourier Transform, JFT)^[17,18]，表示为

联合图拉普拉斯矩阵${{\boldsymbol{L}}_{\mathrm{J}}}$定义为^[17,18]

其中${( \cdot )^{\text{H}}}$为Hermitian矩阵，${{\mathrm{vec}}} ( \cdot )$表示矢量化，${{\boldsymbol{I}}_{{N}}}$和${{\boldsymbol{I}}_{{T}}}$为单位矩阵，$\mathop \otimes \limits^{}$为克罗内克积，$\mathop \oplus \limits^{}$为直积，${{\boldsymbol{\varLambda }}_{{T}}}$为${{\boldsymbol{L}}_{{T}}}$的非负特征值对角矩阵，时域图拉普拉斯矩阵${{\boldsymbol{L}}_{{T}}}$定义为^[17,18]

根据式(7)和式(8)，时变图信号${\boldsymbol{X}}$在图域-时域的联合平滑性可定义为

由于2阶时差信号在图域上比时变图信号本身更具有平滑性^[13]，所以本文将2阶时差信号${\boldsymbol{X}}{{\boldsymbol{D}}_{\mathrm{h}}}$在图域与时变图信号${\boldsymbol{X}}$在时域的联合平滑性构建为${{\mathrm{tr}}} ({({\boldsymbol{X}}{{\boldsymbol{D}}_{\mathrm{h}}})^{\text{T}}}{{\boldsymbol{L}}_{\rm{G}}}({\boldsymbol{X}}{{\boldsymbol{D}}_{\mathrm{h}}})) + {{\mathrm{tr}}} ({\boldsymbol{X}}{{\boldsymbol{L}}_{{T}}}{{\boldsymbol{X}}^{\text{T}}})$。

3.2   基于低秩和联合平滑性的海表面温度数据重构模型

由于海表面温度数据具有显著的时空相关性，所以具有低秩性，故由海表面温度数据建模而成的时变图信号通常近似为低秩的。

利用时变图信号${\boldsymbol{X}}$在图域和时域的联合平滑性，图上每个顶点(即${\boldsymbol{X}}$的每一行)上的时间序列的2阶时差平滑，以及低秩性，存在缺失的时变图信号重构问题可以描述为式(10)的优化问题

其中，${\lambda _1},{\lambda _2} \gt 0$为正则化参数，$\varepsilon \ge 0$为噪声水平，$\circ$为哈达玛积，${{\mathrm{rank}}} ( \cdot )$表示矩阵的秩，${\boldsymbol{J}} \in {\{ 0,1\} ^{N \times T}}$为采样算子，定义为

其中，$v$表示顶点，${{{S}}_t}$表示在$t$时刻被采样的顶点的集合。

因为矩阵秩最小化是NP-hard问题，为了有效地求解式(10)，将其进行凸松弛，从而优化问题转化为

其中，核范数项$||{\boldsymbol{X}}|{|_*}$是指矩阵${\boldsymbol{X}}$中非零奇异值之和，$\mu ,\alpha ,\beta \gt 0$为正则化参数。

上述目标函数中的4项在重构问题中作用如下所示：

(1)第1项$||{\boldsymbol{X}}|{|_*}$为估计矩阵X的低秩性质，保证了时空信号的全局相关性；

(2)第2项$||{\boldsymbol{J}} \circ {\boldsymbol{X}} - {\boldsymbol{Y}}||_{{\mathrm{F}}} ^2$使估计值与观测值保持一致，从而保证了重建结果的可靠性；

(3)第3项${{\mathrm{tr}}} ({({\boldsymbol{X}}{{\boldsymbol{D}}_{\mathrm{h}}})^{\text{T}}}{{\boldsymbol{L}}_{\rm{G}}}({\boldsymbol{X}}{{\boldsymbol{D}}_{\mathrm{h}}}))$保证了2阶时差信号在图域上的平滑性；

(4)第4项${{\mathrm{tr}}} ({\boldsymbol{X}}{{\boldsymbol{L}}_{{T}}}{{\boldsymbol{X}}^{\text{T}}})$保证了时变图信号在时域上的平滑性。

3.3   交替方向乘子法求解

通过观察上述优化问题模型式(12)，可以注意到上述式子的第1项是可逼近的，后3项是可微的，因此选择应用交替方向乘子法^[9]来求解该优化问题。

首先，引入辅助变量${\boldsymbol{Z}}$将上述优化模型式(12)的等价表达表示为

然后，将上述等价模型的增广拉格朗日函数写为

其中，${\boldsymbol{P}}$表示拉格朗日乘子矩阵，$\rho$表示罚参数，$\left\langle , \right\rangle$表示两个矩阵的内积。

交替方向乘子法通过以下迭代方案找到了鞍点

上述子问题式(15)相当于式(18)的形式

采用算法1的共轭梯度法来迭代求解式(18)，将式(18)所示目标函数记为，其梯度为

表 1  共轭梯度法

输入：${\boldsymbol{Y}},{\boldsymbol{J}},{{\boldsymbol{L}}_{\rm{G}}},{{\boldsymbol{L}}_T},{{\boldsymbol{Z}}^k},{{\boldsymbol{P}}^k},{\boldsymbol{D}},\alpha ,\beta ,\rho$，终止迭代阈值$\varepsilon$，最 　大迭代次数$K$
输出：重构信号时变图信号${{\boldsymbol{X}}_{^i}}$
初始化设置：${{\boldsymbol{X}}_i} = 0,\Delta {{\boldsymbol{X}}_i} = 0,i = 0$
迭代：
(1) 确定步长： 　　　 $\tau = - \dfrac{{\left\langle {\Delta {{\boldsymbol{X}}_i},\nabla f({{\boldsymbol{X}}_i})} \right\rangle }}{{\left\langle {\Delta {{\boldsymbol{X}}_i},\nabla f({{\boldsymbol{X}}_i}) + {\boldsymbol{Y}} + \rho {{\boldsymbol{Z}}_i} - {{\boldsymbol{P}}_i}} \right\rangle }}$
其中，$\Delta {{\boldsymbol{X}}_i}$为第$i$步的搜索方向，$\tau$为第$i$步的最优步长，其由 　线性最小化步长准则$\mathop {\min }\limits_\tau f({{\boldsymbol{X}}_i} + \tau \Delta {{\boldsymbol{X}}_i})$决定。
确定步长：
(2) 更新搜索方向：
$\begin{aligned} & {{\boldsymbol{X}}_{i + 1}} = {{\boldsymbol{X}}_i} + \tau \Delta {{\boldsymbol{X}}_i}; \xi {\text{ = }}\frac{{||\nabla f({{\boldsymbol{X}}_{i + 1}})||_{{\mathrm{F}}} ^2}}{{||\nabla f({{\boldsymbol{X}}_{i + 1}})||_{{\mathrm{F}}} ^2}};\\& \Delta {{\boldsymbol{X}}_{i + 1}} = - \nabla f({{\boldsymbol{X}}_{i + 1}}) + \xi \Delta {{\boldsymbol{X}}_i}\end{aligned}$
终止条件：如果$i = K$或者${\text{||}}\Delta {{\boldsymbol{X}}_i}|{|_{\mathrm{F}}} \le \varepsilon$，则停止迭代；否则令$i = i + 1$，继续迭代，直至满足终止条件。

下载: 导出CSV 

| 显示表格

上述子问题式(16)相当于

式(20)的闭式解表示为

其中，SVT为奇异值阈值算子，定义为

其中，${\boldsymbol{U}},\;{\boldsymbol{V}},\;{\boldsymbol{\varSigma }}$来自于${\boldsymbol{X}}$的奇异值分解，即${\boldsymbol{X}} = {\boldsymbol{U\varSigma }}{{\boldsymbol{V}}^{\text{T}}}$, ${{{{\boldsymbol{\varLambda}} }}_\tau }({\boldsymbol{x}}): = {{\mathrm{sign}}} ({\boldsymbol{x}})\max (|{\boldsymbol{x}}| - \tau ,0)$是$\tau \in {R^ + }$定义的软阈值算子。

通过上述分析，本文提出了一种基于低秩和联合平滑性(LRJS)的时变图信号重构算法，算法具体实现步骤如算法2所示。

表 2  LRJS算法求解步骤

输入：${\boldsymbol{Y}},{\boldsymbol{J}},{{\boldsymbol{L}}_{{\mathrm{G}}} },{{\boldsymbol{L}}_{{{T}}} },{{\boldsymbol{Z}}^k},{{\boldsymbol{P}}^k},{\boldsymbol{D}},\alpha ,\beta ,\rho ,\mu$，终止迭代阈值$\varepsilon$， 　最大迭代次数$K$
输出：重构海表面温度数据${{\boldsymbol{X}}^k}$
初始化设置：${{\boldsymbol{X}}^0} = {{\boldsymbol{Z}}^0} = {\boldsymbol{Y}},{{\boldsymbol{P}}^0} = 0,k = 0$
迭代：
(1) 更新${{\boldsymbol{X}}^{k + 1}}$：利用算法1的共轭梯度法；
(2) 更新${{\boldsymbol{Z}}^{k + 1}}$：利用式(21)；
(3) 更新${{\boldsymbol{P}}^{k + 1}}$：利用式(17)；
终止条件：如果$k = K$或者${\text{||}}\Delta {{\boldsymbol{X}}^k}|{|_{{\mathrm{F}}} } \le \varepsilon$，则停止迭代；否则 　令$k = k + 1$，继续迭代，直至满足终止条件。

下载: 导出CSV 

| 显示表格

3.4   算法复杂度分析

为了更好地理解LRJS算法的可拓展性，现在对LRJS算法的计算复杂度进行分析。

在每次迭代当中，矩阵${\boldsymbol{X}}$通过算法1的共轭梯度法进行更新，其中主要计算为式(19)的梯度计算中的矩阵乘法${{\boldsymbol{L}}_{\rm{G}}}{\boldsymbol{X}}{{\boldsymbol{D}}_{\mathrm{h}}}{{\boldsymbol{D}}_{\mathrm{h}}}^{\text{T}}$，其计算消耗为$O({N^2}T + N{T^2})$。在利用共轭梯度法进行迭代时，需要的迭代次数为$O(1/{\ln ((\sqrt {\kappa ({\boldsymbol{A}})} + 1)} / {(\sqrt {\kappa ({\boldsymbol{A}})} - 1)})$，其中${\boldsymbol{A}} = \widehat {\boldsymbol{J}} + \alpha {\widehat {\boldsymbol{L}}_{{\mathrm{G}}} }\widehat {\boldsymbol{D}} + \beta {\widehat {\boldsymbol{L}}_{{{T}}} } + \rho \widehat {\boldsymbol{I}}$, $\widehat {\boldsymbol{J}} = {{\mathrm{diag}}} ({{\mathrm{vec}}} ({\boldsymbol{J}}))$, ${\widehat {\boldsymbol{L}}_{\rm{G}}} = ({{\boldsymbol{I}}_{{T}}} \otimes {{\boldsymbol{L}}_{\rm{G}}})$, $\widehat {\boldsymbol{I}} = {{\boldsymbol{I}}_{{\mathrm{NT}}}}$, $\widehat {\boldsymbol{D}} = ({{\boldsymbol{D}}_{\mathrm{h}}}{{\boldsymbol{D}}_{\mathrm{h}}}^{\text{T}} \otimes {{\boldsymbol{I}}_{{N}}})$, ${\widehat {\boldsymbol{L}}_{{T}}} = ({{\boldsymbol{L}}_{{T}}} \otimes {{\boldsymbol{I}}_{{{N}}} })$^[19]。

上述$\kappa ({\boldsymbol{A}}) = {\sigma _{\max }}({\boldsymbol{A}})/{\sigma _{\min }}({\boldsymbol{A}})$为矩阵${\boldsymbol{A}}$的条件数，${\sigma }_{\mathrm{max}}({\boldsymbol{A}})，{\sigma }_{\mathrm{min}}({\boldsymbol{A}})$分别为矩阵${\boldsymbol{A}}$的最大特征值和最小特征值，所以矩阵${\boldsymbol{X}}$更新的计算复杂度为

下面分析$\kappa ({\boldsymbol{A}})$的上界：

命题1　若${d_{\max }}$是图$G = (V,E,{\boldsymbol{W}})$的所有顶点$N$的度集合$\{ {d_1},{d_2}, \cdots ,{d_{{N}}}\}$的最大度，则它满足

证明　估算矩阵${\boldsymbol{A}}$的最大特征值和最小特征值，其中矩阵$\widehat {\boldsymbol{J}}$, ${\widehat {\boldsymbol{L}}_{\rm{G}}}$, $\widehat {\boldsymbol{I}}$, $\widehat {\boldsymbol{D}}$, ${\widehat {\boldsymbol{L}}_T}$和${\widehat {\boldsymbol{L}}_T}$都为半正定矩阵，并且$\alpha ,\beta \gt 0$。现在有${\sigma _{\min }}({\boldsymbol{A}}) \ge \rho$，对于矩阵${\boldsymbol{A}}$最大特征值的范围，有

因为$\widehat {\boldsymbol{J}}$为一个主对角线上元素只有0或者1的对角矩阵，所以${\sigma _{\max }}(\widehat {\boldsymbol{J}}) \le 1$；此外，由式(3)对矩阵${{\boldsymbol{D}}_{\mathrm{h}}}$的定义，可知矩阵${{\boldsymbol{D}}_{\mathrm{h}}}{{\boldsymbol{D}}_{\mathrm{h}}}^{\text{T}}$为一个方阵，利用Gershgorin circle定理^[20]，可计算得到${\sigma _{\max }}(\widehat {\boldsymbol{D}}) \le 16$；同理，对于图域拉普拉斯矩阵${{\boldsymbol{L}}_{\rm{G}}}$，可计算得到${\sigma _{\max }}({\widehat {\boldsymbol{L}}_{\rm{G}}}) \le 2{d_{\max }}$；对于式(8)中的时域拉普拉斯矩阵${{\boldsymbol{L}}_T}$，可计算得到${\sigma _{\max }}({\widehat {\boldsymbol{L}}_{{T}}}) \le 4$；将${\sigma _{\max }}(\widehat {\boldsymbol{J}}), {\sigma _{\max }}({\widehat {\boldsymbol{L}}_{\rm{G}}}), {\sigma _{\max }}(\widehat {\boldsymbol{D}}),{\sigma _{\max }}({\widehat {\boldsymbol{L}}_{{T}}})$的上界代入式(26)，命题1得证。

应用基本不等式$\ln (1 + x) \ge x/(1 + x), \forall x \ge 0$和命题1，则有

令$K = 1 + \dfrac{{\sqrt {1 + 32\alpha {d_{\max }} + 4\beta } }}{{2\sqrt \rho }}$，则矩阵${\boldsymbol{X}}$的更新需消耗$O(K \cdot ({N^2}T + N{T^2}))$；

矩阵${\boldsymbol{Z}}$的更新中，SVT计算占主要地位，式(22)中的矩阵乘积计算复杂度为$O(\min ({N^2}T,N{T^2}))$，而矩阵进行奇异值分解(Singular Value Decomposition,SVD)计算的复杂度也是$O(\min ({N^2}T,N{T^2}))$，则矩阵${\boldsymbol{Z}}$的更新计算复杂度为$O(\min ({N^2}T,N{T^2}))$；对于矩阵${\boldsymbol{P}}$，其每次更新只涉及到矩阵求和与矩阵标量积，因此其计算复杂度为$O(NT)$。

综上所述，本文所提LRJS算法的计算消耗由矩阵${\boldsymbol{X}}$的更新所主导，总的计算复杂度为$O(K \cdot ({N^2}T + N{T^2}))$。

3.5   算法误差上界分析

令${\boldsymbol{X}}$是真实的时变图信号，$\widehat {\boldsymbol{X}}$是重构的时变图信号，由下述命题2可以得到LRJS算法的误差上界。

命题2　若对于任意矩阵${\boldsymbol{M}} \in {R^{N \times T}}$，存在常数$\gamma \in [0,1)$，有

则LRJS算法有误差上界

证明　对于误差上界，有

首先，为了处理A项，讨论收敛点$\widehat {\boldsymbol{X}}$，利用Karush-Kuhn-Tucker条件，有

共同应用上述3个条件，则有

通过式(31)，有

根据文献[21]的定理3，对于任何$\tau \in {R^ + }$，有${{{\mathrm{SVT}}} _\tau } = {({\boldsymbol{I}} + {\partial _{|| \cdot |{|_*}}})^{ - 1}}$，在式(32)中，通过令$\tau = 1/16$，有

根据文献[21]引理1中关于近端算子的性质，即对于任意两个矩阵${{\boldsymbol{Y}}_1},{{\boldsymbol{Y}}_2}$，以及任意$\tau \in {R^ + }$，有

现在将式(33)代入式(29)中的A项，并应用式(34)，有

进一步考虑式(27)，并且令${\boldsymbol{M}} = \widehat {\boldsymbol{X}} - {\boldsymbol{X}}$有

接下来，对于B项

其中，${\boldsymbol{\varSigma }}$是以矩阵${\boldsymbol{X}} - \dfrac{\alpha }{{16}}{{\boldsymbol{L}}_{\rm{G}}}{\boldsymbol{X}}{{\boldsymbol{D}}_{\mathrm{h}}}{{\boldsymbol{D}}_{\mathrm{h}}}^{\text{T}} - \dfrac{\beta }{{16}}{\boldsymbol{X}}{{\boldsymbol{L}}_{{T}}}$的奇异值为对角线元素的对角矩阵。

最后将式(36)和式(37)的上界代入式(29)，则可证得命题2。

4.   仿真实验

为测试本文所提LRJS方法的重构性能， 将本文所提方法与GraphTRSS算法^[8] 、LRDS算法^[13]、JSR-TVGS算法^[14]、近似SVD分解的固定点迭代法 (Fixed Point Continuation with Approximate SVD, FPCA)^[21]、自然邻点插值法(Natural Neighbor Interpolation, NNI)^[22]5种方法进行比较。

仿真实验中，利用图信号处理工具箱GSPBOX^[23]构建图$G = (V,E,{\boldsymbol{W}})$。利用k最近邻法构造一个无向加权图来表征观测点的地理邻接关系，为了得到一个连接的连通图，取$k{\text{ = }}10$。设置随机采样比例为10%～90%，对应数据缺失率为90%～10%。在每次测试中，采样数据的节点均是随机的，并在不同采样率情况下分别进行100次蒙特卡洛测试，最大迭代次数为$K{\text{ = }}2 \times {10^3}$，终止迭代阈值$\varepsilon {\text{ = }}{10^{{{ - }}6}}$。

采用均方根误差(Root Mean Square Error, RMSE)和平均绝对误差(Mean Absolute Error, MAE)两个指标来评价方法的重构性能，计算公式为

其中，$\widehat {\boldsymbol{x}}$表示重构信号，${\boldsymbol{x}}$为真实信号，${N_{\boldsymbol{x}}}$是信号长度。MAE和RMSE越小，表示重构数据与真实数据的匹配度越高，即重构算法的重构精度越高。

本文所有仿真实验的硬件环境均为配置2.50 GHz Intel i5-7200U处理器、8 GB内存的计算机，软件环境为Matlab R2022a。

4.1   实验结果

(1)太平洋海海平面温度数据集。该实验数据集来自美国NOAA地球系统研究实验室发布的海面温度数据^[24]，该数据集覆盖从1870年～2014年采集的月平均海表面温度数据，空间分辨率为1°×1°全球网格。在太平洋海域上随机选择100个位置采集的数据进行实验，选取观测时间为200个月，因此数据大小为100×200。

算法中设置的参数为$\mu = 0.5,\alpha = 0.01,\beta = 0.001, \rho = 1$。利用GSPBOX构建的太平洋海域海表面温度数据采集连接图如图1所示。对于该数据集，不同方法在不同采样率条件下海表面温度重构的MAE, RMSE和运行时间结果如图2所示。

图 1  太平洋海域海表面温度数据采集连接图

下载: 全尺寸图片 幻灯片

图 2  所有方法在太平洋海域海表面温度数据集上的性能

下载: 全尺寸图片 幻灯片

(2) 南海海平面温度数据集。该实验使用ARGO南海海平面温度数据集^[25]作为实验数据集。该数据集的空间分辨率为1°×1°全球网格，实验中，在南海海域随机选择66个地点进行模拟，选取的观测时间为2010年到2019年共计120个月，因此数据的大小为66×120。

算法中设置的参数为$\mu = 0.5,\alpha = 0.01,\beta = 0.001, \rho = 0.5$。太平洋海域海表面温度数据采集连接图如图3所示。对于该数据集，不同方法在不同采样率条件下海表面温度重构的MAE,RMSE和运行时间结果如图4所示。

图 3  南海海域海表温度数据采集连接图

下载: 全尺寸图片 幻灯片

图 4  所有方法在南海海域海表温度数据集上的性能

下载: 全尺寸图片 幻灯片

4.2   实验结果分析

上述仿真实验主要从MAE, RMSE以及算法运行时间3个方面比较了各算法的性能，具体分析如下：

(1)从图2和图4可以看出，随着采样率的提高，LRJS算法的MAE和RMSE逐渐降低，这是因为随着采样率的增加提供了更充分的图信号先验信息，从而减小了重构误差，这验证了本文算法的合理性。在相同的采样率条件下，与其他5种算法相比，本文所提LRJS算法表现出更小的MAE和RMSE，表明其重构精度更高，这是因为LRJS算法充分利用了更多的先验知识。

(2)从图2(a)和图2(b)可知，在太平洋海域海平面温度数据集的实验中，随着采样率的增加，本文所提LRJS算法相较于LRDS算法而言，MAE下降了6%～34%，RMSE下降了4%～32%，相较于JSR-TVGS算法而言，MAE下降了37%～45%，RMSE下降了39%～51%。

(3)从图4(a)和图4(b)可知，随着采样率的增加，在南海海平面温度数据集的实验中，本文所提LRJS算法相较于LRDS算法而言，MAE下降了5%～19%，RMSE下降了4%～17%；本文所提LRJS算法相较于JSR-TVGS算法而言，MAE下降了20%～31%，RMSE下降了17%～29%。

(4)但是从图2(c)和图4(c)可知，在算法重构效率方面，在相同的采样率下，LRJS算法的运行时间仅仅低于FPCA算法，这是因为在每次迭代过程中，核范数都需要进行SVD分解和利用共轭梯度法进行子问题的更新，导致计算成本较高。

上述仿真实验结果表明，本文所提LRJS算法尽管在重构效率方面，未能表现出太大优势，但是提高了海表温度的重构精度，综合考虑所有进行比较的性能指标，本文所提LRJS算法仍占有优势。

5.   结束语

本文研究了基于图信号处理的海表温度数据重构问题。通过将海表面温度数据建模为时变图信号，并将时变图信号重构问题转化为凸优化问题，利用海表面温度的低秩性，以及该信号在图域和时域联合平滑性，提出了LRJS算法，并采用交替方向乘子法来进行求解，最终实现时变图信号重构，解决海表温度数据重构问题。在两种不同的真实数据集上进行了仿真实验，结果表明，与其他相关的时变图信号重构算法相比，本文所提LRJS算法能够获得更高的重构精度，表明了LRJS算法在时变图信号重构方面的有效性。下一步的研究将探索在确保较小重构误差的基础上，改进优化算法的重构效率，以更好地满足实际应用需求。