CN110780604A - 一种基于空时平滑性和时间相关性的空时信号恢复方法 - Google Patents

Abstract

本发明公开了一种基于空时平滑性和时间相关性的空时信号恢复方法该方法为一种从不完整并且有噪声的样本中恢复空时信号的算法，在这个算法中，将具有网络结构的空时信号视为时变图信号，其中图的拓扑结构不仅由空间结构决定，而且还受时变数据的时间相关性的影响。利用空时平滑性和时间相关性进行建模，把信号恢复问题从数学上建模为一个优化问题，然后把这个非凸的优化问题分解为两个子优化问题，这两个子问题分别迭代的恢复空时信号和时间相关矩阵。应用交替方向乘子法来解决第一个子问题以恢复空时信号，接着用CVX优化工具包解决第二个子问题以更新时间相关矩阵，这两步交替迭代，最终恢复出来完整的空时信号。

Description

一种基于空时平滑性和时间相关性的空时信号恢复方法

【技术领域】

本发明属于信号处理技术领域，具体涉及一种基于空时平滑性和时间相关性的空时信号恢复方法。

【背景技术】

在过去的15年里，图表已经渗透到维数据的复杂的社会、生物和技术系统的研究中，或网络的动态过程中，如流行病的传播。在网络科学中，研究者感兴趣的问题不是通过分析数据来解决的，而是通过分析图表的结构来解决。因此关于分析和处理图上信号的理论日渐成熟。这些理论主要分为两类，第一类是图信号处理。由于已经存在关于离散信号处理的成熟理论，研究者首先在图上引入了带限信号和相应的采样定理。然后建立了一个全面的图信号框架，并将传统理论推广到图形信号的研究中，给出了图信号相关概念的一些重要定义，如滤波，卷积，脉冲响应等。在具有理论框架之后，出现了许多采样和恢复图信号的方法。此外，研究者认为图形信号会随着时间的推移而发生变化，因此提出了恢复时变图信号的相关方法。首先，基于时变图信号在每个时刻都受到严格限制的假设，图信号处理的发展提供了许多有效的方法来解决这个问题。第二类是空时信号处理，而在空时信号处理中，最重要的一个问题就是如何从残缺的空时信号出发，恢复出完整的空时信号。当恢复时空信号时，总会利用一些空时信号固有的特性，也就是相关性，即空时信号在空间和时间上具有多种不同的相似性。利用不同的相似性，衍生出了许多不同的空时信号恢复方法。而目前的空时信号恢复方法，都没有充分的利用空时信号的相关性，导致恢复出来的空时信号存在误差。因此为了降低空时信号恢复的误差，有必要在充分的考虑各种相关性的基础和前提下，提出一种新的空时信号恢复方法。

另一方面，当今社会的网络结构十分庞大且复杂，从这些网络采集数据时由于各种客观条件会导致采集到的数据不全。如图1所示，以传感器网络为例，在传感器节点采集数据并把数据传输到基站的过程中，传感器节点损坏或者通信故障都会导致数据缺失，但接下来如果要进行数据分析的话很有可能需要完整的数据，在这种情况下，就需要从已知的数据出发来恢复完整的数据。数据本身有着高的相关性，这种相关性分为两个方面，一方面是时间上的相关性，即随着时间推移信号变化较为平缓；另一方面是空间上的相关性，即地理位置相近的数据相似度较高。而这些数据本身具有的相关性质都可以被用来改善恢复数据的性能。

【发明内容】

本发明的目的在于克服上述现有技术的缺点，提供一种基于空时平滑性和时间相关性的空时信号恢复方法；该方法用于解决恢复的空时信号存在误差，导致恢复的空时信号不完整的问题。

为达到上述目的，本发明采用以下技术方案予以实现：

一种基于空时平滑性和时间相关性的空时信号恢复方法，包括以下步骤：

步骤1，将采样空时信号恢复至原始空时信号的过程转化为优化问题，对应公式为：

其中，J为采样算子，X为原始空时信号，Y为采样空时信号，R为时间相关矩阵，tr(·)为矩阵的迹，α,β,γ,μ为正的正则化参数；L为拉普拉斯矩阵；

为矩阵的哈达马积；

步骤2，通过交替方向乘子法将采样空时信号恢复为原始空时信号，在每一次交替方向乘子法恢复空时信号的过程中，同时通过CVX优化工具包恢复时间相关矩阵R；每一次空时信号的恢复均基于上一次恢复的空时信号和时间相关矩阵R；直至交替方向乘子法迭代完成设定次数，完成空时信号恢复。

本发明的进一步改进在于：

优选的，步骤1中，所述原始空时信号X为无向加权图G的空时信号，L为无向加权图G的拉普拉斯矩阵，L的计算式为：

L＝D-W (2)

其中，D为无向加权图G的度矩阵，W为无向加权图G的权重矩阵。

优选的，步骤1中，采样空时信号Y的计算公式为：

其中，V为小方差的高斯噪声。

优选的，公式(10)中，tr(D(X)^TLD(X))表示采样空时信号Y的差分平滑性，其中D(X)为修正后的时间差分信号，其计算公式为：

D(X)＝X-RXB＝[a₁,a₂-Ra₁,…,a_M-Ra_M-1] (8)

其中，a_i＝[a_1,i,a_2,i,…,a_N,i]，为在第i个时刻采集到的信号；B为差分矩阵。

优选的，公式(10)中，tr(X^TLX)为采样空时信号Y的单一时刻的平滑性。

优选的，步骤2中，恢复原始空时信号X过程中，迭代恢复每一次X^k+1对应的公式为：

将公式(11)转化为等效分裂形式：

公式(12)对应的增广拉格朗日函数为：

其中，A为一个优化变量，P为拉格朗日乘子，ρ为惩罚项系数，<·,·>表示两个矩阵的内积；

针对式(13)对原始空时信号X通过迭代k次恢复，具体包括以下步骤：

步骤2.1.初始化X⁰＝Y,Z⁰＝Y,P⁰＝0,R⁰＝I,k＝0；其中I为单位矩阵；

步骤2.2.计算k＝1时的X¹；

步骤2.3.通过k＝1时的X¹和通过CVX工具包计算出的R¹计算出k＝2时的X²；

步骤2.4.通过k＝2时的X²和通过CVX工具包计算出的R²计算出k＝3时的X³，依次类推，直至k为设定次数，计算出最终的X^E，为恢复的最终空时信号X^E。

优选的，在步骤2.2-步骤2.4的任意一次迭代中，通过交替乘子迭代算法恢复原始空时信号X，具体来说，每一次对三个变量X^k、Z^k和P^k进行迭代更新，三个变量的更新公式分别为：

P^k+1＝P^k+1+ρ(X^k+1-Z^k+1) (16)

其中，SVT为奇异值阈值算子；C为一个优化变量；

通过共轭梯度方法求解(14)，计算当前k值对应的X^k+1，通过X^k+1结合式(15)计算出Z^k+1，通过X^k+1和Z^k+1，结合式(16)计算出P^k+1；

将X^k+1、Z^k+1和P^k+1输入至CVX工具包计算出R^k+1。

优选的，所述求解式(14)的共轭梯度方法包括以下步骤：

(1)确定搜索步长τ：

(2)更新变量X的值：

X_i+1 ^k＝X_i ^k+τΔX_i ^k (19)

(3)更新搜索方向λ：

ΔX_i+1 ^k＝-▽f(X_i+1 ^k)+λΔX_i ^k (21)。

优选的，所述奇异值阈值算子SVT的定义如下：

其中，U,V,Σ为X的奇异值分解，Λ_τ(x)＝sign(x)max(|x|-τ,0)为软阈值算子。

与现有技术相比，本发明具有以下有益效果：

本发明公开了一种基于空时平滑性和时间相关性的空时信号恢复方法；该方法为一种从不完整并且有噪声的样本中恢复空时信号的算法，在这个算法中，将具有网络结构的空时信号视为时变图信号，其中图的拓扑结构不仅由空间结构决定，而且还受时变数据的时间相关性的影响。利用空时平滑性和时间相关性进行建模，把信号恢复问题从数学上建模为一个优化问题，然后把这个非凸的优化问题分解为两个子优化问题，这两个子问题分别迭代的恢复空时信号和时间相关矩阵。应用交替方向乘子法来解决第一个子问题以恢复空时信号，接着用CVX优化工具包解决第二个子问题以更新时间相关矩阵，通过两步交替迭代，最终恢复出来完整的空时信号。

【附图说明】

图1为空时信号模型及其拓扑示例图；

图2为海平面温度数据库中采样率是0.4时各种方法各个时刻的信号平均绝对误差比较；

图3为海平面温度数据库中不同采样率下各种方法的100个时刻的总信号平均绝对误差比较；

图4为加州每日平均PM2.5数据库中采样率是0.3时各种方法各个时刻的信号平均绝对误差比较；

图5为加州每日平均PM2.5数据库中不同采样率下各种方法的100个时刻的总信号平均绝对误差比较。

【具体实施方式】

下面结合附图和具体过程对本发明做进一步详细描述：

以一个无向加权图G(V,ε,W)作为待恢复的图表，其中V代表N个节点的集合，ε代表所有边的集合，对称非负矩阵W代表该无向加权图的权重矩阵。矩阵W中第i行第j列的元素w_ij定义为连接第i个节点与第j个节点的边的权重。当第i个节点与第j个节点之间有边时，w_ij不为0；否则w_ij为0。该无向加权图的度矩阵定义为：

D＝diag(d₁,…,d_N) (1)

其中，

定义为第i个节点的度，且该图的拉普拉斯矩阵定义为：

L＝D-W (2)

步骤1，建立原始空时信号与采样空时信号之间的关系

针对无向加权图G，一个总共在M个时刻内被采集到的空时信号被定义为矩阵X＝[a₀,a₁,…,a_M-1],其中a_i＝[a_1,i,a_2,i,…,a_N,i]为在第i个时刻采集到的信号。

空时信号的恢复即从它的采样中恢复出完整的信号，信号采样模型为：

其中，X为原始空时信号，Y为采样空时信号，V为小方差的高斯噪声，J定义为采样算子，

为矩阵的哈达马积。

步骤2，建立空时信号矩阵的相关性数学表示

恢复空时信号时，可以利用空时信号本身所具有的相关性。相关性分为两个部分，第一部分为全局相关性，即每次观测的信号都来自有限的模式，在数学上表示为整个原始空时信号X是一个低秩矩阵，相关性的第二部分为局部相关性，即平滑性。平滑性指出，基于给定拓扑，空时信号在时间和空间上都应该是缓慢变化的，即平滑的。平滑是一个定性的表示，需要找到平滑性的数学表示来对信号恢复问题进行数学建模。针对单一时刻采样的信号，一种典型的平滑性的数学表示为：

S(a)越小，该信号就越平滑。上述是针对单一时刻采样的信号求平滑性，因此需要将其推广到所有时刻的信号，即针对无向加权图G，把所有单一时刻信号的平滑性进行求和，得到如下表达式：

其中tr(·)定义为矩阵的迹。

步骤3，建立空时信号的差分平滑性

信号不仅自身具有平滑性，而且信号的差分也具有平滑性，即信号变化量的改变也是缓慢的，该性质称之为差分平滑性。为了应用这个性质，定义差分矩阵B，这个矩阵的具体形式为：

因此，对于原始空时信号X的时间差分信号为：

X-RXB＝[a₁,a₂-a₁,…,a_M-a_M-1] (7)

因为图中的每个节点在时间上的相关程度不同，加入时间相关矩阵R修正原始空时信号X的时间差分信号，修正过后的时间差分信号为：

D(X)＝X-RXB＝[a₁,a₂-Ra₁,…,a_M-Ra_M-1] (8)

其中R为时间相关矩阵，该矩阵是一个对角矩阵，且它的第i个对角元素r_i定义为第i个节点的相关系数。把该修正过的时间差分信号代入式(5)，即表示为原始空时信号X的差分平滑性；因此信号的差分平滑性可表示为：

S(D(X))＝tr(D(X)^TLD(X)) (9)

步骤4，将空时信号的恢复问题进行优化

利用空时信号的空时相关性，将空时信号的恢复问题建模为如下的优化问题：

其中α,β,γ,μ为正的正则化参数，优化问题中的第一项控制了原信号与恢复信号的误差；第二项应用了式(9)所示的差分平滑性；在第三项中，为了应用信号的低秩特性即全局相关性，通过约束矩阵的核范数来约束矩阵的秩；第四项反映了信号的空时平滑性，即式(5)，最后一项控制了时间相关矩阵R的对角线元素大小。此外，优化问题中的第一个限制条件控制时间相关矩阵R的对角线元素，第二个限制条件控制了不同节点之间时间相关系数的差异程度。

通过上述公式(10)恢复原始空时信号X的基本思路为通过不断的迭代，将设定的初始X⁰通过反复的迭代得到最终恢复的X^E，在每一次迭代恢复空时信号的过程中，同时通过CVX优化工具包恢复时间相关矩阵R；因此优化后的问题(10)有两个优化变量X^k和R^k，然而该问题对这两个优化变量并不是联合凸的。因此为了解决该问题，在本发明中，采用块坐标下降的方案，即在迭代的每步中，保持其中一个优化变量不变，对另一个优化变量进行更新。因此，采用这种思想，一个复杂的优化问题被分解为两个分别针对单一优化变量的子问题，接下来分别解决这两个子问题的方法：

1)恢复原始空时信号X的方法：

在恢复原始空时信号X的过程中，每一次迭代恢复X^k+1，如下所示：

公式(11)是针对恢复空时信号X^k+1的子问题，A为一个优化变量；而该优化问题对于X^k+1来说是一个凸优化问题，因此采用交替乘子迭代算法来解决该问题，采用该算法的原因是该算法具有可分解性和优越的收敛性能。首先，将公式(11)转化为能够应用在交替乘子迭代算法上的等效分裂形式：

优化问题(12)的增广拉格朗日函数为：

其中，P为拉格朗日乘子，ρ为惩罚项系数，<·,·>表示两个矩阵的内积。

当运用交替乘子迭代算法时，需要对三个变量X^k、Z^k和P^k进行迭代更新，这个三个变量的更新规则分别为：

P^k+1＝P^k+ρ(X^k+1-Z^k+1) (16)

其中，C为一个优化变量；

其中，SVT定义为奇异值阈值算子，具体定义如下：

SVT_τ(X)＝UΛ_τ(∑)V^T (17)

其中，U,V,Σ为X的奇异值分解，Λ_τ(x)＝sign(x)max(|x|-τ,0)为软阈值算子。

为了找到子问题(14)的解，可以直接令目标函数的导数为0。但这种解法会带来十分高的计算复杂度，因为它涉及大矩阵求逆的操作。因此在本发明中采用共轭梯度方法来求解子问题(14)以降低计算复杂度，具体步骤为：

(4)先根据如下准则来确定搜索步长τ

(5)再根据如下准则更新变量X的值

X_i+1 ^k＝X_i ^k+τΔX_i ^k (19)

(6)最后根据如下准则更新搜索方向λ：

ΔX_i+1 ^k＝-▽f(X_i+1 ^k)+λΔX_i ^k (21)

通过上述共轭梯度方法计算出X^k，将X^k分别带入公式(15)和(16)中，计算出Z^k和P^k。

2)恢复时间相关矩阵R的方法

针对恢复时间相关矩阵R的子问题也是一个凸优化问题，因此采用凸优化工具包CVX来解决这个问题，将上面计算出的X^k、Z^k和P^k输入至CVX后，恢复了时间相关矩阵R^k。

通过上述的公式(14)-(16)求解两个子问题，具体的步骤为：

(1)计算k＝1时的X¹,Z¹,P¹和R¹；

(2)计算k＝2时的X²,Z²,P²和R²；

(3)计算k＝3…的X³,Z³,P³和R³，直至k为最大迭代次数，k的迭代次数为设定值。

上述迭代结束后，完成了空时信号X和时间相关矩阵R的恢复，得到最终恢复的X^E和R^E，公式(10)的优化问题完成。

上述整个完整的优化算法，如表1所示。

表1恢复空时信号X和时间相关矩阵R的优化算法

直到k达到最大迭代次数，迭代次数为根据情况设定。

实施例

为了验证本发明所提出算法的信号恢复性能，在两个现实数据库上进行仿真，这两个数据库分别是海平面温度数据库和加州每日平均PM2.5数据库。在仿真中，对每个数据库，首先构建一个无向加权图，接着对数据在不同采样率下进行采样，采样策略为简单随机抽样。由于数据总量太大，在仿真中只分析并恢复前100个时刻的数据。图2～图5表明在这两种数据库下，本发明所提算法在恢复准确度方面更具优势。

从图2中可以看出，对于海平面温度数据库来说，在每一时刻，当不考虑相关矩阵时运用本发明所提方法进行数据恢复，效果已经明显好于现存的其他方案，当进一步考虑相关矩阵后，恢复性能又进一步提升。从图3中可以看出，对于不同采样率，总的空时信号的恢复性能也会有着明显提升。而从图4、图5中也能看出，对于加州PM2.5数据库来说，也有着相似的结论。但恢复性能的提升并没有之前在海平面温度数据库中那么明显，这可能是因为海平面温度数据库更贴近本发明中的模型假设，而PM2.5数据库会和本发明的模型稍有偏差。但尽管是在这种情况下，本发明提出的方案依然有恢复性能的提升。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (9)

1.一种基于空时平滑性和时间相关性的空时信号恢复方法，其特征在于，包括以下步骤：

步骤1，将采样空时信号恢复至原始空时信号的过程转化为优化问题，对应公式为：

其中，J为采样算子，X为原始空时信号，Y为采样空时信号，R为时间相关矩阵，tr(·)为矩阵的迹，α,β,γ,μ为正的正则化参数；L为拉普拉斯矩阵；

为矩阵的哈达马积；

步骤2，通过交替方向乘子法将采样空时信号恢复为原始空时信号，在每一次交替方向乘子法恢复空时信号的过程中，同时通过CVX优化工具包恢复时间相关矩阵R；每一次空时信号的恢复均基于上一次恢复的空时信号和时间相关矩阵R；直至交替方向乘子法迭代完成设定次数，完成空时信号恢复。

2.根据权利要求1所述的一种基于空时平滑性和时间相关性的空时信号恢复方法，其特征在于，步骤1中，所述原始空时信号X为无向加权图G的空时信号，L为无向加权图G的拉普拉斯矩阵，L的计算式为：

L＝D-W (2)

其中，D为无向加权图G的度矩阵，W为无向加权图G的权重矩阵。

3.根据权利要求1所述的一种基于空时平滑性和时间相关性的空时信号恢复方法，其特征在于，步骤1中，采样空时信号Y的计算公式为：

其中，V为小方差的高斯噪声。

4.根据权利要求1所述的一种基于空时平滑性和时间相关性的空时信号恢复方法，其特征在于，公式(10)中，tr(D(X)^TLD(X))表示采样空时信号Y的差分平滑性，其中D(X)为修正后的时间差分信号，其计算公式为：

D(X)＝X-RXB＝[a₁,a₂-Ra₁,…,a_M-Ra_M-1] (8)

其中，a_i＝[a_1,i,a_2,i,…,a_N,i]，为在第i个时刻采集到的信号；B为差分矩阵。

5.根据权利要求1所述的一种基于空时平滑性和时间相关性的空时信号恢复方法，其特征在于，公式(10)中，tr(X^TLX)为采样空时信号Y的单一时刻的平滑性。

6.根据权利要求1所述的一种基于空时平滑性和时间相关性的空时信号恢复方法，其特征在于，步骤2中，恢复原始空时信号X过程中，迭代恢复每一次X^k+1对应的公式为：

将公式(11)转化为等效分裂形式：

公式(12)对应的增广拉格朗日函数为：

其中，A为一个优化变量，P为拉格朗日乘子，ρ为惩罚项系数，<·,·>表示两个矩阵的内积；

针对式(13)对原始空时信号X通过迭代k次恢复，具体包括以下步骤：

步骤2.1.初始化X⁰＝Y,Z⁰＝Y,P⁰＝0,R⁰＝I,k＝0；其中I为单位矩阵；

步骤2.2.计算k＝1时的X¹；

步骤2.3.通过k＝1时的X¹和通过CVX工具包计算出的R¹计算出k＝2时的X²；

步骤2.4.通过k＝2时的X²和通过CVX工具包计算出的R²计算出k＝3时的X³，依次类推，直至k为设定次数，计算出最终的X^E，为恢复的最终空时信号X^E。

7.根据权利要求6所述的一种基于空时平滑性和时间相关性的空时信号恢复方法，其特征在于，在步骤2.2-步骤2.4的任意一次迭代中，通过交替乘子迭代算法恢复原始空时信号X，具体来说，每一次对三个变量X^k、Z^k和P^k进行迭代更新，三个变量的更新公式分别为：

P^k+1＝P^k+1+ρ(X^k+1-Z^k+1) (16)

其中，SVT为奇异值阈值算子；C为一个优化变量；

通过共轭梯度方法求解(14)，计算当前k值对应的X^k+1，通过X^k+1结合式(15)计算出Z^k+1，通过X^k+1和Z^k+1，结合式(16)计算出P^k+1；

将X^k+1、Z^k+1和P^k+1输入至CVX工具包计算出R^k+1。

8.根据权利要求7所述的一种基于空时平滑性和时间相关性的空时信号恢复方法，其特征在于，所述求解式(14)的共轭梯度方法包括以下步骤：

(1)确定搜索步长τ：

(2)更新变量X的值：

X_i+1 ^k＝X_i ^k+τΔX_i ^k (19)

(3)更新搜索方向λ：

9.根据权利要求7所述的一种基于空时平滑性和时间相关性的空时信号恢复方法，其特征在于，所述奇异值阈值算子SVT的定义如下：

其中，U,V,Σ为X的奇异值分解，Λ_τ(x)＝sign(x)max(|x|-τ,0)为软阈值算子。

Images