叶文海, 林红波

(1. 吉林省科维交通工程有限公司 工程部, 长春 130021； 2. 吉林大学 通信工程学院, 长春 130012)

0 引 言

地震勘探数据中存在大量的随机噪声, 受复杂的地表条件和自然环境影响, 地震随机噪声通常表现为非高斯、 非平稳的特性。这些复杂的随机噪声存在严重干扰反射地震信号的辨识和提取, 掩盖了地震数据中蕴含的地下结构信息。此外, 恶劣的采集条件及复杂的地形地貌易导致地震检波器工作异常, 某些地区无法布置检波器, 这些因素使地震数据中存在地震道缺失问题, 影响利用地震勘探数据对地下成像的精度。因此, 压制地震勘探随机噪声并重构地震缺失道是地震勘探数据处理的重要且基础环节, 能提高地震数据质量, 有助于改善地下结构成像的精度[1-3]。

目前地震勘探学者已经提出了许多地震信号恢复方法, 并在实际应用中取得了显著效果。从低信噪比地震数据中重构地震信号是一类病态的反问题, 需要利用地震信号的先验作为约束以提高降噪和地震信号恢复的效果。边缘导向的滤波方法能提供平滑的地震信号, 例如偏微分方程[4-5]、 形态学滤波[6], 然而在缺失数据的情况下信号恢复效果不理想。稀疏先验通过从学习的过完备字典中选择合适的原子对地震信号进行稀疏表示, 但其计算稀疏系数和字典更新的计算量大[7]。全变分是一种常用且有效的信号先验[8], 全变分先验利用梯度刻画地震信号结构特征, 且优化问题是凸的, 能在压制地震随机噪声的同时保持复杂信号结构, 然而同相轴边缘易被平滑[9]。各向异性全变分及其高阶版本挖掘地震信号的形态结构特征, 从而调整滤波强度和方向实现沿同相轴滤波, 增强对信号结构的保持能力[10]。针对变分模型仅学习当前数据低层次特征的问题, Liang等[11]将变分模型和可学习模型相结合, 充分利用了期望对数似然块先验表征地震数据的精细结构, 从而改善对地震信号恢复能力。近年来, 结合地震信号低秩特性的降噪算法在地震噪声压制方面取得了显著效果[12], 常用于恢复高斯噪声下的地震信号, Hankel矩阵的引入则进一步提高了低秩降噪算法的有效性[13-14]。上述信号恢复模型具有灵活的降噪框架, 可以嵌入不同层次的目标数据特征提升滤波效果, 对非均匀、 不平稳、 多类型的地震噪声具有一定的鲁棒性, 但在低信噪比情况下获得的先验易受到随机噪声的干扰而降低地震信号恢复质量。

流形学习是一种数据驱动的信号分析方法, 其核心假设为高维空间获得的观测数据本质上位于一个更低维的流形空间, 可通过寻找信号的低维表示解决信号处理和图像处理中的信号恢复问题[15]。将信号的低维作为正则项的低维流形模型(LDMM: Low Dimensional Manifold Model)是一种有效的信号恢复方法[16-17], 该模型将数据块作为光滑低维流形的均匀采样, 利用图拉普拉斯刻画数据块光滑流形的块间相似或相关结构, 则最小化流形维度可通过最小化图拉普拉斯正则项实现。地震信号恢复问题中, 地震同相轴具有横向相干特性, 使数据块之间具有相似特性, 基于块流形的假设, 地震数据块可看作低维光滑流形的采样。Yu等[18]验证了地震数据块流形分布的低维特性, 通过寻找使维度最低的块流形实现了地震随机噪声压制和缺失道恢复。块流形学习在低维流形上利用局部数据块间的内在关联, 同时对相似的块进行去噪, 解决传统方法不能有效处理具有复杂结构地震数据的不足。然而, 现有方法未考虑地震数据的非局部特征而不能很好地恢复地震信号的全局形态。

笔者将框架变换和低维流形模型相结合, 利用框架变换的局部基和非局部基全面分析地震块间的内在关联, 借助框架变换和低维流形的图拉普拉斯矩阵之间的内在联系, 将块流形的低维约束问题转化为框架变换系数的能量结构化稀疏问题, 提出基于卷积框架变换的低维流形模型消除地震噪声并恢复地震信号。基于卷积框架变换的低维流形模型既能利用框架变换提供地震信号的结构化稀疏特征, 又能充分利用块流形的局部和非局部特性, 实现复杂地震噪声消减与缺失道重构相结合的地震有效信号恢复。

1 相关方法

笔者利用框架变换系数的能量集中模式与低维流形模型(LDMM)的低维正则项之间的关联, 将低维流形学习的非局部基与数据块驱动的局部基相卷积构建紧框架小波, 从局部和非局部联合地刻画非平稳地震信号特征, 重构位于高维空间中光滑的低维流形以恢复地震信号。

1.1 卷积框架变换

卷积框架变换是一种能以局部和全局联合刻画信号的表示方法[19]。对N点长实值信号y∈RN, 令信号块第i行Yi=(y(i),…,y(i+l-1))∈Rl, 0≤i≤N-1为在i点开始长为l的行矢量, 按信号块在信号中的顺序垂直堆叠构建信号的块矩阵

(1)

其中上角标T为转置操作。块矩阵是Hankel矩阵, 包含了冗余的信号信息。块矩阵的第i列矢量Yi由信号的循环移位构造, 包含了信号的全部信息, 而块矩阵的行矢量Yi为信号的块, 包含了信号的局部信息。若将各行看作信号块构成的点, 块矩阵则为信号块的点云函数。

令Φ和V为两个正交矩阵, 其维度分别为N×N和l×l,Φ的列矢量φi和V的列矢量vj的外积构造了所有N×l矩阵在RN×l空间的正交基, 表示为

(2)

利用正交基将块矩阵Y改写为

(3)

其中tr(·)为矩阵的迹。由式(3)可知, 两个不同的基Φ和V在块矩阵的行空间和列空间联合地表示块矩阵。由于信号y可通过块矩阵的反对角平均重构且块矩阵Y与矢量vj的乘积为矢量y与矢量vj翻褶的卷积[19], 所以块矩阵Y可表示为

(4)

(5)

则{ψij|1≤i≤N, 1≤j≤l}为定义在RN上函数的紧框架。

信号y的卷积框架变换定义为

Cij=〈y,φi*vj〉

(6)

其中Cij为框架变换系数。利用紧框架变换, 信号y可分解为

(7)

由式(7)可知, 利用基Φ的列矢量φi和V的列矢量vj的卷积框架小波可近似地表示信号y。基V的列矢量对应块矩阵的行, 刻画了信号的局部特征, 在信号近似中的作用类似于对信号进行局部加窗变换。因此, 将构造的基V称作卷积框架变换的局部基, 可由块矩阵的右奇异矢量构造[20]。相应地, 基Φ的列矢量对应块矩阵的列, 由于块矩阵的每列均为信号的循环移位, 对块矩阵的列表示是一种对信号的非局部表示。考虑到块矩阵沿列方向可看作数据块的点云, 对列的非局部表示等效于在卷积框架变换下基Φ对块矩阵的非局部近似。因此, 基Φ称作卷积框架变换的非局部基。

L=D1/2(I-D-1W)D-1/2

(8)

其中W为加权邻接矩阵, 尺寸为N×N, 利用高斯函数度量第i个数据点Yi与第j个数据点Yj之间的非相似度, 公式如下

Wij=exp(-‖Yi-Yj‖2/ρ)

(9)

对归一化图拉普拉斯特征分解：L=ΦΛΦT,Λ的对角元素为降序排列的特征值λ1,…,λN。基于扩散映射原理,ΦΛ1/2的列是正交的, 可作为块集谱嵌入的坐标函数。由于干净数据是Rl的p维嵌入, 则可选择与L的最小p个特征值对应的ΦΛ1/2的p列作为非局部基函数, 使原数据所有点对的距离与嵌入空间对应点的距离近似。

卷积框架变换利用局部基和非局部基对信号稀疏表示, 能使地震信号的框架变换系数具有能量集中模式, 这意味着地震信号或数据的系数能量集中在少数低频分量。另一方面, 从非线性降维中获得的非局部基被视为定义数据集嵌入的坐标函数, 用相对较少的基函数捕捉尽可能多的数据集内变化。这种基函数构建方案能通过耦合非局部(局部)基而放大局部(非局部)基的能量集中效应。因此, 笔者首先给定非局部基, 利用最大特征值对应的列给出块矩阵行的近似等距嵌入, 进而利用线性重构代价函数求取局部基, 使卷积框架变换系数的能量集中在系数矩阵的左上部分。

1.2 低维流形模型

在低维流形模型中, 基于LDMM的降噪算法利用块流形的维数作为正则项, 使块流形的维度尽可能小以恢复原始信号。

流形是位于RN的子空间点的集合, 流形的维度即为这些点的自由度。地震勘探记录

y=Sf+n

(10)

其中n为噪声,S为观测矩阵, 式(10)表示记录过程中对地震信号f的降质操作, 例如数据的丢失, 下采样等。在地震数据降噪任务中, 观测矩阵为一个单位阵。所有地震数据块集P(y)定义为

P(y)={R[y(i)]∈Rl:i∈Ω}

(11)

其中R[y(i)]为在索引i从地震数据y提取的地震数据块, 块尺寸为q×q, 令l=q2。Ω为地震记录中所有采样点的集合。在低维流形模型中, 所有地震信号块的集合位于光滑的低维流形上[18], 等距地嵌入在欧几里得空间Rl。块流形的一个重要特征是低维, 因此, 地震数据块的低维特性可用作正则项, 使块流形的维度尽可能小, 以从含噪的观测数据中恢复地震信号f。为此, 基于LDMM的降噪模型可表示为优化目标函数

(12)

其中μ为正则化参数, 第1项d(M(y))为块流形M(y)的维度, 对一个等距嵌入在Rl中光滑子流形M, 其维度可利用嵌入空间中在块流形M(y)上的梯度近似计算, 得到

(13)

其中αj为流形M(y)上的第j个坐标函数, 即x=(α1(x),…,αl(x)), ∀x∈M(f)⊂Rl, 对应块矩阵的第j列。M为在黎曼流形M上的梯度算子。将块流形维度带入式(12), 优化问题改写为

(14)

利用分裂Bregman迭代方案求解上述优化问题。首先固定流形M, 更新数据f； 然后固定数据f, 对每个坐标函数优化维正则项, 更新流形。重复上述步骤, 优化问题会收敛到一个最优解, 得到滤波后的地震信号f。

在第m次迭代, 更新流形可设置为

(15)

根据获得的流形, 更新坐标函数和数据。对维正则化优化子问题解耦, 并利用点积分方法将其离散化为

(16)

2 基于卷积框架变换的低维流形降噪模型

信号的块流形位于光滑的低维流形上是基于低维流形模型的降噪算法的一个重要假设。由此, 地震数据去噪问题可转化为在嵌入空间中求取维度最低的光滑流形的优化问题。将嵌入空间信号流形的低维特性作为正则项, 约束信号重构。笔者进一步考虑黎曼流形上低维特性与卷积框架变换系数结构化稀疏的关系, 将基于数据驱动的非局部基和局部基的卷积框架变换系数能量集中模式等效为流形维正则项, 利用地震信号块流形特性对地震数据降噪, 提高地震信号的恢复能力。

2.1 卷积框架变换与流形维的关系

在LDMM降噪算法中, 低维正则项优化子问题通过最小化流形维二次型更新块矩阵的列矢量。从流形的角度出发, 块矩阵的列矢量可理解为Rl的坐标函数, 块矩阵是Rl中点云的表示函数。如果点云是Rl的低维子空间的采样, 则可将点云嵌入到一个低维欧几里得空间, 且每两点之间的距离与原空间相比没有明显的畸变。在卷积框架变换下, 如果块矩阵存在一个好的低维嵌入Φ, 选择与Φ适应的局部基V, 则系数矩阵的能量集中在左上部三角区。上述分析表明, 更低的维度对应系数矩阵左上部更小的区域, 即系数能量集中。因此, LDMM中的低维正则项等效于卷积框架变换系数在左上部的区域集中, 将框架小波系数能量作为正则项的卷积框架低维流形降噪, 等价于LDMM降噪。下面具体分析卷积框架变换与流形维度的关系。

(17)

(18)

其中第1项为流形维正则项, 第2项是数据保真项, ‖·‖F,D1/2=‖D1/2·‖F表示D1/2加权的F范数, 为表述方便省略了式(16)中的角标(m-1), 在流形维正则项中,RL=D1/2((I-L)-1-I)D-1/2。为分析二者之间的联系, 将流形维正则项改写为与框架小波系数二次型类似的形式, 则

(19)

上述分析表明, 卷积框架变换系数能量与流形维正则化项的二次型是等效的, 最小化流形维等效于卷积框架小波变换系数能量的集中模式。

2.2 基于卷积框架小波的低维流形模型

基于卷积框架小波系数能量集中和流形低维的等价关系, 笔者将LDMM模型中的流形低维正则项改写为更简单的卷积框架变换系数能量的正则项, 称为基于卷积框架小波的低维流形模型, 定义为

(20)

通过最小化目标函数对地震数据进行降噪。求解过程与LDMM相同, 其中正则项优化子问题是求解的关键, 笔者利用点积分方法, 将欧拉拉格朗日公式转换为线性系统

(21)

基于卷积框架变换的低维流形模型(CFR-LDMM: Convolutional Framelet Regularization Based Low Dimensional Manifold Model)总结如算法1。

算法1 基于卷积框架小波正则的低维流形模型(CFR-LDMM)。

输入： 地震数据y

初始化迭代次数为M,Y(0)=y,

当m

1) 构造Hankel矩阵Y(m)；

2) 由式(9)计算Y(m)的邻接矩阵W(m)和度矩阵D(m)；

3) 对块流形Hankel矩阵部分SVD(Singular Value Decomposition)分解构建局部基V(m)；

4) 令E=Y(m)-d(m)；

5) 计算处理后的Hankel矩阵在局部基下的系数矩阵Z(m)；

6) 对Z(m)进行局部基反变换得到Hankel矩阵Y(m+1)；

7) 对Y(m+1)进行反对角平均得到y(m+1)；

8) 令m=m+1。

结束循环

3 实验和结果

笔者利用合成地震勘探记录和实际地震勘探记录测试所提出的基于卷积框架变换的低维流形地震信号恢复模型, 并与降维算法主分量分析(PCA： Principal Components Analysis), 以及稀疏变换方法曲波变换的降噪效果进行比较。同时进一步计算滤波后地震数据的信噪比和均方误差, 定量评估几种去噪方法的降噪性能。

3.1 合成地震勘探记录

首先将提出的信号恢复模型应用于合成地震勘探记录以测试去噪性能。采用的合成地震记录(见图1a)包含4条同相轴, 由雷克子波模拟地震信号, 主频分别为20 Hz,18 Hz,18 Hz,20 Hz, 采样时间为2 ms。在合成地震记录中加入高斯白噪声, 生成含噪合成记录, 如图1b所示, 其信噪比为-2.35 dB。由于随机噪声的存在, 含噪记录中的反射地震信号变得不清晰, 难于辨识和恢复同相轴。

利用PCA降噪方法、 曲波降噪方法和所提出的CFR-LDMM方法消减含噪合成地震记录中的随机噪声。在这些方法中, PCA利用前10个最大主成分进行信号重构, 曲波降噪算法在曲波域对系数进行硬阈值处理, CFR-LDMM利用50×10的数据块构建块矩阵, 利用SVD分解的右奇异矢量构建局部基, 得到的滤波结果如图2a～图2c所示。对比滤波结果可见, PCA方法虽然有效压制了随机噪声, 但有效反射信号也同时被消减。曲波变换方法在信号恢复方面虽然有所改善, 但信号恢复效果仍不太理想, 且信号边缘被模糊。与前两种方法相比, 笔者降噪方法给出了更干净的背景和更清晰的信号, 有效反射信号得到很好的恢复。

为进一步评估上述3种算法的滤波结果, 分别计算了3种滤波方法在不同噪声强度下降噪结果的信噪比和均方误差(MSE： Mean Squared Error), 计算公式如下

(22)

(23)

表1 不同方法滤波前后信噪比和均方误差比较

图1 合成记录Fig.1 Synthetic seismic data

图2 合成记录信号恢复效果比较Fig.2 Comparison of signal recovery of synthetic seismic data with PCA, curvelet and our method

为验证对实际地震噪声的有效性, 笔者进一步在合成地震记录中加入实际随机噪声进行测试。加入的实际随机噪声为中国塔里木某矿区采集的实际噪声记录, 得到的含噪记录如图3a所示。从图3a中可见, 实际地震随机噪声具有慢变的波形, 与模拟的反射信号的波形相似, 噪声强度分布不均匀, 遮掩甚至截段了地震同相轴。利用PCA、 曲波降噪和笔者方法处理含实际噪声的合成地震记录, 滤波结果如图3b～图3d所示。对比滤波结果可见, 3种方法均能较好地压制实际随机噪声并恢复有效信号。然而, PCA滤波结果中存在信号损失问题, 具有复杂结构的同相轴没有得到恢复。曲波和笔者方法均能较好地恢复地震信号, 相比之下, 笔者方法恢复的信号同相轴更连续, 背景更干净。其降噪效果在滤波结果与干净数据的差图像中得到进一步验证(见图4)。与PCA和曲波方法相比, 笔者方法得到的差图像中有更少的信号残留。

上述结果表明, 笔者提出的基于卷积框架变换正则化低维流形地震信号恢复模型能从低品质地震勘探数据中有效恢复地震信号, 消减地震随机噪声, 地震信号结构得以良好恢复, 地震数据的信噪比得到显著提升。

图3 实际噪声污染的合成记录信号恢复效果比较Fig.3 Comparison of signal recovery of the synthetic seismic data contaminated by field random noise with PCA, curvelet and our method

图4 合成记录去噪前后的差图比较Fig.4 Comparison of the difference before and after denoising

3.2 实际地震勘探记录处理

笔者将基于卷积框架变换的低维流形模型应用于实际沙漠地震勘探数据处理, 采用的实际沙漠地震勘探记录来源于中国塔里木盆地沙漠地区, 该记录包含100个地震道, 每道由2 000个采样点组成, 采样频率为500 Hz, 对此实际数据随机抽取, 地震道缺失率为25%, 如图5a所示。从含噪地震记录可以看到, 随机噪声波形变化缓慢,噪声较强, 且强度随时间和空间的变化而变化。有效信号受到了噪声的严重干扰,部分区域的同相轴结构被沙漠噪声破坏导致其难以识别。笔者分别采用PCA滤波、 曲波变换和基于卷积框架变换的低维流形模型对该地震记录进行处理。3种模型的结构和参数与合成数据分析的参数相同。

图5b～图5d分别给出了3种方法的信号恢复结果, 可看出PCA滤波对沙漠噪声的抑制能力有限, 滤波结果中仍然存在大量的沙漠噪声, 缺失的地震道没有得到很好的恢复。 曲波降噪方法实现了对大部分背景噪声的有效压制, 缺失道得到恢复。但许多有效信号被过度地平滑。相比之下, 基于卷积框架变换的低维流形模型能更有效地恢复反射地震信号, 缺失道恢复得更好, 信号更为光滑、 连续。 此外, 在矩形框中被强干扰掩盖的信号得以更好地恢复。

图5 实际地震记录信号恢复效果比较Fig.5 Comparison of signal recovery of the field seismic data with PCA, curvelet and our method

4 结 语

笔者提出了一种基于卷积框架变换的低维流形地震信号恢复模型, 借助卷积框架变换系数的结构稀疏约束实现块流形的低维约束, 既能利用低维块流形的自相似特性提升地震信号恢复能力, 又可避免显示定义流形坐标, 使低维流形模型的应用拓展到地震信号恢复。理论分析和沙漠地震勘探数据处理结果验证了基于卷积框架变换正则化的低维流形模型对同相轴的恢复能力, 在极低信噪比条件下实现了低频随机噪声的压制, 同时实现了缺失道的有效恢复。

参考文献：

[1]DONG X, ZHONG T, LI Y. A Deep-Learning-Based Denoising Method for Multiarea Surface Seismic Data [J]. IEEE Geoscience and Remote Sensing Letters, 2021, 18(5): 925-929.

[2]SANG W, YUAN S, YONG X, et al. DCNNs-Based Denoising with a Novel Data Generation for Multidimensional Geological Structures Learning [J]. IEEE Geoscience and Remote Sensing Letters, 2021, 18(10): 1861-1865.

[3]FANG Z, LIN H,XU X. Mask-Guided Model for Seismic Data Denoising [J]. Geoscience and Remote Sensing Letters, 2022, 19: 8026705.

[4]ZHANG Y, LIN H, LI Y, et al. Low-Frequency Seismic Noise Reduction Based on Deep Complex Reaction-Diffusion Model [J]. IEEE Transactions on Geoscience and Remote Sensing, 2022, 60： 5903814.

[5]林红波, 张丹丹. 自适应结构方向复扩散压制沙漠地震噪声 [J]. 吉林大学学报(信息科学版), 2019, 37(5): 463-469.

LIN Hongbo, ZHANG Dandan. Adaptive Structure-Oriented Complex Diffusion Filtering for Desert Seismic Random Noise Attenuation [J]. Journal of Jilin University (Information Science Edition), 2019, 37(5): 463-469.

[6]HUANG W, WANG R, WANG S, et al. F-XY Domain Morphological Filtering for Three-Dimensional Seismic Erratic Noise Suppression [J]. Geophysical Journal International, 2019, 216: 81-102.

[7]张良, 韩立国, 方金伟, 等. 双稀疏字典和FISTA的地震数据去噪 [J]. 地球物理学报, 2019, 62(7): 2671-2683.

ZHANG Liang, HAN Liguo, FANG Jinwei, et al.Seismic Data Denoising via Double Sparsity Dictionary and Fast Iterative Shrinkage-Thresholding Algorithm [J]. Chinese J Geophys, 2019, 62(7): 2671-2683.

[8]CHAN S H, KHOSHABEH R, GIBSON K B, et al. An Augmented Lagrangian Method for Total Variation Video Restoration [J]. IEEE Transactions on Image Processing, 2011, 20(11): 3097-3111.

[9]ZHANG Y, LIN H, LI Y, et al. Seismic Signal Enhancement and Noise Suppression Using Structure-Adaptive Nonlinear Complex Diffusion [J]. IEEE Transactions on Geoscience and Remote Sensing, 2020, 58(3)： 2198-2211.

[10]GEMECHU D, YUAN H, MA J. Random Noise Attenuation Using an Improved Anisotropic Total Variation Regularization [J]. Journal of Applied Geophysics, 2017, 144: 173-187.

[11]LIANG Chao, LIN Hongbo, MA Haitao. Low-Frequency Seismic Random Noise Attenuation Based on EPLL-TV under Double Prior Constraints [J]. Journal of Applied Geophysics, 2022, 203: 104689.

[12]FENG J, LIU X, LI X, et al. Low-Rank Tensor Minimization Method for Seismic Denoising Based on Variational Mode Decomposition [J]. IEEE Geoscience and Remote Sensing Letters, 2022, 19: 7504805.

[13]MA H, WANG Y, LI Y, et al. Desert Seismic Low-Frequency Noise Attenuation Using Low-Rank Decomposition-Based Denoising Convolutional Neural Network [J]. IEEE Transactions on Geoscience and Remote Sensing, 2022, 60: 5900809.

[14]李佳, 马海涛, 李月. 自适应秩收敛低秩算法压制沙漠地震随机噪声 [J]. 吉林大学学报(信息科学版), 2021, 39(3): 237-245.

LI Jia, MA Haitao, LI Yue. Low-Rank Algorithm Based on Adaptive Rank Convergence for Desert Seismic Random Noise Attenuation [J]. Journal of Jilin University (Information Science Edition), 2021, 39(3): 237-245.

[15]ZENG J, CHEUNG G, NG M, et al. 3D Point Cloud Denoising Using Graph Laplacian Regularization of a Low Dimensional Manifold Model [J]. IEEE Transactions on Image Processing, 2020, 29: 3474-3489.

[16]OSHER S, SHI Z, ZHU W. Low Dimensional Manifold Model for Image Processing [J]. SIAM Journal on Imaging Sciences, 2017, 10(4): 1669-1690.

[17]CHEN X, YE Y, WU Q, at al, Fast Manifold Ranking with Local Bipartite Graph [J]. IEEE Transactions on Image Processing, 2021, 30: 6744-6756.

[18]YU S, OSHER S, MA J, et al. Noise Attenuation in a Low-Dimensional Manifold [J]. Geophysics, 2017, 82(5): V321-V334.

[19]YIN R, GAO T, LU Y M, et al. A Tale of Two Bases: Local-Nonlocal Regularization on Image Patches with Convolution Framelets [J]. SIAM Journal on Imaging Sciences, 2017, 10(2): 711-750.

[20]陶顺, 郭傲, 刘云博, 等. 基于矩阵束和奇异值分解的间谐波检测算法 [J]. 电力系统保护与控制, 2021, 49(2)： 57-64.

TAO Shun, GAO Ao, LIU Yunbo, et al. Interharmonic Detection Algorithm Based on a Matrix Pencil and Singular Value Decomposition [J]. Power System Protection and Control, 2021, 49(2)： 57-64.

[21]PANG J, CHEUNG G. Graph Laplacian Regularization for Image Denoising: Analysis in the Continuous Domain [J]. IEEE Transactions on Image Processing, 2017, 26(4): 1770-1785.

[22]AXELSSON O, LINDSKOG G. On the Rate of Convergence of the Preconditioned Conjugate Gradient Method [J]. Numerische Mathematik, 1986, 48(5): 499-523.