，获得一系列固有模态分量（Intrinsic Mode Function, IMF）及残余分量。因经过EMD分解后的IMF分量无法确定有用信号和含噪分量的信噪分界点，所以引入自相关函数确定含噪分量，对含噪分量进行小波阈值去噪，得到新的IMF分量。最后将IMF分量与残余分量重构得到降噪后的时间参数序列
。摘要 针对铁路继电器参数的噪声问题，为提取其有效信息，该文建立一种基于改进波形匹配延拓法优化的经验模态分解（EMD）算法，利用自相关函数分离含噪信号，小波阈值去除噪声的混合降噪模型。首先根据继电器参数特点对波形匹配方法进行改进，重新定义匹配误差度公式，并引入匹配精度误差系数，采用改进波形匹配延拓法优化EMD分解过程产生的端点效应，得到有效的固有模态分量（IMF）和余项；然后求解其自相关函数，并根据自相关函数图像结合噪声信号特征分离出含噪分量；最后对含噪分量进行小波阈值去噪，去噪后与剩余分量和余项结合，得到重构后的参数序列。同时，提出利用结构相似性（SSIM）评价指标，结合信噪比（SNR）、方均误差（MSE）指数对模型可靠度评判。通过结果分析，并与EMD分解后重构和小波阈值去噪方法作对比，证明该模型可优化铁路继电器参数的降噪效果。
关键词：端点效应 波形匹配 端点延拓法 自相关函数 结构相似性（SSIM）指数
铁路继电器是实现控制目标、信号传递、电路隔离的关键设备，在铁路系统中发挥着重要的作用。它的可靠性将直接影响整个系统的稳定和安全，因此对铁路继电器参数进行有效分析是十分必要的[1]。继电器参数是反映其工作特性及可靠性的一项重要指标[2]，对继电器参数的有效分析，不仅有助于分析继电器失效的机理，还可为其性能分析和寿命预测提供重要数据支持。由于参数数量达百万之多，因此对继电器参数进行预处理是十分必要的，而滤波降噪方法的选择正是前期处理的重中之重。
小波阈值去噪是目前常用的继电器参数降噪方法，针对单一方法滤波降噪的不足，其他领域的学者提出了此算法与经验模态分解（Empirical Mode Decomposition, EMD）相结合的降噪方法。文献[3]提出基于EMD分解的小波阈值降噪算法，应用于轴承、齿轮等零件的振动信号降噪分析，解决了离心泵振动信号复杂等问题。文献[4]在电力系统低频振荡模态参数的提取中，结合EMD和小波的优点，提高了抗模态混叠能力。文献[5]将总体经验模态分解（Ensemble Empirical Mode Decomposition, EEMD）与新型小波阈值去噪相结合用于谐波检测。文献[6]提出了一种采用经验模态分解与小波分析相结合的方法，能较好地分析原功率的特征，更好地还原原有功率值，从而获得更为准确的风电并网功率的参考值。研究结果表明，将EMD分解与小波阈值去噪相结合，能够解决小波阈值去噪算法对含噪分量分解不细致的问题，并且在EMD分解的基础上利用小波去噪进行二次去噪，能够有效抑制高频白噪声，使联合方法优于单一方法。但EMD分解产生的端点效应问题还需进一步改善。
目前，广大学者已提出一些抑制EMD端点效应的方法。文献[7]提出镜像连续法，克服了由于镜放置位置不是极值点而产生的连续误差。这种方法直观性强，操作也简单，但是对信号本身的对称性要求高，在处理对称性不强的信号时，容易产生误差。文献[8]提出一种基于改进极值波延拓与对称中点插值相结合的端点效应抑制方法。此方法在对比趋势时，对于内部波形与边缘波形的联系考虑较少。文献[9]提出一种基于时间尺度的端点延拓方法，应用到原始风速信号的数据处理中。此方法对于数据量较大的参数序列运行时间较长，不适合本文所用超百万次的试验数据。文献[10]提出了一种基于自适应序贯相似性检测波形匹配延拓的端点效应抑制方法，提高了子波截取的合理性与匹配的准确性。此方法将端点处数据与内部数据波形相关联，合理匹配波形，抑制端点效应。
本文总结上述方法的优缺点，以铁路继电器的两个电参数和3个时间参数为研究对象，利用改进波形匹配的端点延拓方法优化EMD分解算法，引入自相关函数确定信噪分界点。采用小波阈值去噪对含噪分量进行降噪处理，将处理后的分量与其余分量和余项进行重构得到降噪序列，通过降噪前后的折线图对比及信噪比（Signal Noise Ratio, SNR）、方均误差（Mean Squared Error, MSE）和结构相似性（Structural Similarity Index Measurement, SSIM）指数的求解，相比单独利用EMD分解后重构以及小波阈值去噪方法进行降噪处理，本文所用模型使降噪后的序列保留了原始序列的基本趋势，优化了降噪效果。
本小节提出一种基于改进波形匹配延拓法优化的EMD算法，利用自相关函数分离含噪信号，小波阈值去除噪声的混合降噪模型。改进的EMD_小波阈值去噪模型结构如图1所示，采用改进波形匹配延拓法修正EMD分解的端点效应，将改进EMD分解方法应用于铁路继电器的时间参数序列，获得一系列固有模态分量（Intrinsic Mode Function, IMF）及残余分量。因经过EMD分解后的IMF分量无法确定有用信号和含噪分量的信噪分界点，所以引入自相关函数确定含噪分量，对含噪分量进行小波阈值去噪，得到新的IMF分量。最后将IMF分量与残余分量重构得到降噪后的时间参数序列。
图1 改进的EMD_小波阈值去噪模型
Fig.1 Improved EMD_ wavelet threshold denoising model
EMD分解是一种不需要先验基底的自适应分解方法，可以将信号分解为从高频到低频排列的若干个IMF分量[3]。设时间参数序列为，经EMD分解后的原始信号可表示为
（1）
式中，
为第i个固有模态分量；R为残余分量。
EMD算法在分解继电器时间参数时，因信号的两个端点不一定是极值点，从而导致分解后得到的分量曲线在序列的两端出现发散，这被称为端点效应[11]。端点效应会随着分解的进行逐渐向曲线内部发展，污染原始数据，导致结果出现较大偏差，甚至使分解的结果失去意义。
波形匹配端点延拓法通过在序列内部寻找与边缘波形匹配度较高的波形进行延拓，同时能够考虑序列的边缘特性进行延拓[12]。但基本的波形匹配方法一般直接采用匹配子波和待匹配子波幅值差的绝对值，并不考虑匹配子波的幅值，这样很难直观反映两者的匹配精度[10]。本文结合继电器参数的数据特点及分解效果，采用改进波形匹配的端点延拓方法对继电器参数序列的端点振荡问题进行改善。以左端延拓为例，对此方法进行详细描述：
（1）对于参数序列，取其最左端的第
个极值点分别为
和
（假设为极大值，为极小值），选择更靠近右端的极值点位置，记从序列初始值到此极值点的这段波形为
，此时参数动作次数为
。
（2）假设步骤（1）中取点为极小值点，设序列的极小值点集合为
，以中除之外的右端其余各点
为参考点，取为右端点向左延拓动作次数次的波形
，计算和的匹配度
，有
（2）
式中，
为
序列的第i点；
为
序列的第i点。
（3）为比较各段波形匹配度，得到最合适的匹配子波，定义匹配精度误差
为
（3）
使最小，得到波形匹配成功的序列，此时波形匹配度为
，若此序列仅包含一个极值点，则取
左端包含此极小值点的子波作为序列左端的延拓，延拓完成，转步骤（4）。
（4）右端序列的延拓与左端类似。假设右端序列延拓了
个数据点，对延拓后的序列重新进行EMD分解，得到以
为数据长度的IMF分量，将各分量左端减少个数据，右端减少个数据，得到原始数据被分解后的分量。
以进行极小值点延拓过程为例，需要说明的是，一般波形匹配的待匹配子波选取序列初始点到最左端的极小值点，在进行波形匹配时，由于本文所分解的数据量较大，为适应大数据达到更好的分解效果，故根据不同参数分解效果选择经过的极值点个数。
随机噪声在各个时刻具有随机性，关联性很弱，其自相关函数在自变量的零点上达到最大后，会迅速衰减到0左右[13]。但对具有物理意义的信号，时间点之间有一定的关联性，因此其自相关函数值不会迅速衰减。
现定义序列的自相关函数为IMF分量在两个不同动作次数时的相关程度
，有
（4）
式中，
为序列X(t)中不同元素乘积的均值；
为第n个不同元素的乘积。
利用得到的序列画出各个分量的自相关函数图像，根据其自相关函数图像的特点寻找符合噪声信号的分量，确定信噪分界点。
一般而言，经EMD分解后的含噪分量多为高频分量，小波阈值去噪算法对高频分量的抑制效果较强[14]。小波阈值去噪流程如下：
（1）根据具体含噪信号选择适合的小波基以及分解层数N，然后将含噪信号进行N层小波分解，得到相应的小波分解系数。
（2）对分解得到的小波序列选用合适的阈值函数进行阈值处理，对细节量进行去噪。
（3）对阈值处理过后的小波量进行重构，得到去噪信号。
本文选用db5小波作为小波基函数，选择软阈值函数作为进行阈值的函数，其分解层数则根据不同分量的不同特点进行选取，小波去噪可最大程度保留原始序列真实性并去除序列中的噪声信号。
在得到降噪参数序列的情况下，可根据不同的评判指标对降噪效果进行评价。本文引用SSIM、SNR、MSE这三个指数对模型进行评判。
图像的结构携带了图像最主要的信息，SSIM从图像组成的角度解释结构信息，图像的结构信息由亮度、对比度和结构度3个要素组成[15]。SSIM指数在[0, 1]范围内，指数越大的两个图像越相似。SSIM可建模为
（5）
其中
（6）
式中，w、y分别为原始数据折线图和滤波后数据折线图；
为亮度比较函数；
为对比度比较函数；
为结构对比函数，为常数；
、
为原始数据折线图的亮度均值和标准差；
、
为滤波后数据折线图的亮度均值和标准差；
为原始数据折线图和滤波后数据折线图的相关度；c1、c2和c3为常数，避免分母为0带来的系统错误；
、
和
为大于0的常数，实际工程计算中一般设定
。
需要说明的是，在进行SSIM指数求解时，主要对降噪前后参数的折线图进行结构对比，考虑到图例和坐标轴等内容的影响，输入的两张图像只保留折线进行SSIM指数的求解。
信噪比和方均误差是衡量去噪效果的两个重要指标，铁路继电器参数的信噪比和方均误差公式定义[16]为
（7）
（8）
式中，为铁路继电器参数的原始序列；为降噪后的参数序列。信噪比越大，方均误差越小，说明降噪效果越好。
铁路相关应用对继电器的特殊要求包括动作必须可靠准确，有稳定的电气特性和时间特性，在周围温度和湿度变化很大的情况下，保持很高的电气绝缘强度，铁路继电器的具体电气参数包括电参数、时间参数和机械参数等[17]。为此开展铁路继电器可靠性试验并测试铁路继电器的各类电气参数。
目前，针对铁路继电器在低电平运行测试设备控制系统运行的试验条件及参数见表1。
表1 铁路继电器电寿命试验条件及参数
Tab.1 Conditions and parameters of electrical life test for railway relay
条件参 数 试验环境温度: 40℃湿度: 65% 动作试验平台低电平运行测试设备控制系统 试验条件线圈电压: DC 24V触点负载: DC 24V/1A阻性 试验样品铁路继电器 样品数量8台继电器同时运行 检测参数电参数: 线圈稳流值、线圈稳压值等时间参数: 吸合时间、释放时间、动断超程时间、动合超程时间等
为便于观察检测参数的波动趋势，选取线圈稳流值、线圈稳压值、吸合时间、动断超程时间、动合超程时间这5个参数作图分析，各参数序列原始波形如图2所示。
图2 各参数原始序列波形
Fig.2 Original sequence waveforms of each parameter
通过观察图2中不同参数的变化趋势，发现吸合时间序列的波动趋势较为显著，而其他参数没有显著的下降或上升的趋势，故以继电器的吸合时间参数为例进行降噪模型介绍，原始参数序列为，对进行EMD分解，部分分解结果如图3所示。
图3所示为吸合时间序列经EMD分解后出现端点效应的部分分量，其中不难发现，
都有较为明显的端点振荡现象，这种现象即为端点效应。这会使接下来的去噪和重构出现极大的误差，导致降噪后的序列波形与原始序列结构差距较大，整个降噪模型将失去意义。
为此使用第1.1节中提出的改进波形匹配的端点延拓法对EMD算法进行改进，在进行左端端点延拓时，经验所得在经历10个极值点时分解效果最佳，第10个极大值点所在位置为38，第10个极小值点所在位置为37，故
，待匹配子波序列为
，计算此序列与各匹配序列的匹配精度误差，使最小，得到波形匹配成功的序列为
，将匹配子波延拓到原始序列左端，右端延拓与左端类似，延拓的序列长度
，最终得到的匹配子波序列为
。延拓后重新对吸合时间序列进行EMD分解，得到的部分分解结果如图4所示。
图3 吸合时间序列EMD分解的部分分量
Fig.3 The partial components of EMD decomposition of accretion time series
对比图3和图4的分解结果，经过改进波形匹配延拓法优化的EMD分解使端点效应得到了明显的抑制，纵坐标范围的大幅减小正体现了其端点振荡的改善，解决了端点处数值与波形内部数值相比高出较多的问题。
经分解共得到20个固有模态函数分量及其残余分量，按第1.2节中介绍的算法计算其自相关函数，得到吸合时间序列的IMF分量自相关函数图像，如图5所示。
由图5可知，噪声主导的高频分量的自相关函数呈现高频振荡，继电器动作次数越多，振荡幅度越大，随着动作次数的增加，迅速减小到0位置上下，这是典型的噪声自相关函数图像。而到了
之后的8个分量则出现了明显的放缓，说明其中的相关性越来越强，也就是信号越来越强。由此可知，前12个分量符合随机噪声特征，说明前12个分量信息中主要含有随机噪声。
图4 吸合时间序列经改进EMD分解后的部分分量
Fig.4 The partial components of the time series after EMD decomposition
图5 吸合时间序列的IMF分量自相关函数图像
Fig.5 Image of autocorrelation function of IMF component in time series
选用db5小波作为小波基函数，阈值函数选择软阈值函数，分解层数根据不同分量的不同特点进行选取，对含有噪声信号的12个分量进行小波阈值去噪，去噪后与剩余的8个分量和余项R进行重构。吸合时间序列的降噪后参数序列图像与原始序列图像对比如图6所示。
利用本文所述模型依次对继电器的线圈稳流值、线圈稳压值、动断超程时间、动合超程时间参数序列进行降噪处理，得到的降噪结果对比如图7所示。
观察图7中继电器不同参数的降噪前后对比，发现降噪后的序列不仅能保持其随动作次数增大的变化趋势，而且能够使其带宽减小，达到降噪效果。下面采用第2节中的评价指标表征具体的降噪效果。利用不同去噪方法（EMD分解与重构、小波阈值去噪、本文所用方法）得到的SNR、MSE和SSIM指数求解结果如图8所示。
图6 吸合时间序列降噪前后对比
Fig.6 Comparison chart of suction time series before and after noise reduction
由图8针对SSIM、MSE和SNR指数的对比，对于吸合时间、线圈稳流值、线圈稳压值、动断超程时间以及动合超程时间这5个参数序列来说，其信噪比和结构相似度指数都得到提升，方均误差降低，可见本文所用方法能将EMD分解和小波阈值去噪方法的优点相结合，有效提高降噪效果。
图7 继电器各参数序列降噪前后对比
Fig.7 Comparison of relay parameters before and after noise reduction
本文采用铁路继电器参数作为试验数据，建立了一种基于改进波形匹配延拓法优化的EMD算法，利用自相关函数分离含噪信号，小波阈值去除噪声的混合降噪模型。并得出以下结论：
1）在对铁路继电器参数序列进行EMD分解时，根据分解后的IMF分量两端出现的端点振荡问题，提出基于波形匹配的端点延拓法，并在波形匹配时，将匹配子波的幅值考虑进去，重新定义了匹配误差度公式，并引入匹配精度误差系数，让单纯的匹配误差转化为具有统一标准的匹配精度误差，优化后的EMD分解有效抑制了端点效应引发的波形失真。
2）利用自相关函数确定信噪分界点，对比自相关函数图像和噪声信号特点分离出含噪分量，在一定程度上解决了仅依靠IMF分量图像来确定噪声信号的片面性。
图8 不同模型去噪效果对比
Fig.8 Comparison of different models' denoising effects
3）在一般的降噪评价指标SNR和MSE基础上，引入SSIM指数评价降噪效果，从各指标论证本文所用方法能够保留原始参数序列的结构，弥补EMD分解与小波阈值去噪的不足，提高降噪效果。
在后续工作中，将利用本文降噪后的参数序列进行继电器性能退化过程及寿命预测的研究。
参考文献
[1] 苏伟龙, 许志红. 高压直流继电器磁吹系统的建模与设计[J]. 电工技术学报, 2022, 37(6): 1583-1594.
Su Weilong, Xu Zhihong. Modeling and design of magnetic blowing system for high voltage direct current relay[J]. Transactions of China Electro- technical Society, 2022, 37(6): 1583-1594.
[2] 李奎, 张国盼, 郑淑梅, 等. 基于实时服役参数的交流接触器电寿命最大化控制策略[J]. 电工技术学报, 2021, 36(9): 1976-1985.
Li Kui, Zhang Guopan, Zheng Shumei, et al. A control strategy for maximizing the electrical life of AC contactors based on real-time operating parame- ters[J]. Transactions of China Electrotechnical Society, 2021, 36(9): 1976-1985.
[3] Sun Biao, Zhou Shaoping, Wang Congyi. Application of wavelet soft threshold denoising algorithm based on EMD decomposition in vibration signals[C]//2019 6th International Conference on Systems and Infor- matics (ICSAI), Shanghai, China, 2019: 7-11.
[4] 陈琳, 陈静, 王惠民, 等. 基于小波包能量熵的电池剩余寿命预测[J]. 电工技术学报, 2020, 35(8): 1827-1835.
Chen Lin, Chen Jing, Wang Huimin, et al. Prediction of battery remaining useful life based on wavelet packet energy entropy[J]. Transactions of China Electrotechnical Society, 2020, 35(8): 1827-1835.
[5] 孙曙光, 庞毅, 王景芹, 等. 一种基于新型小波阈值去噪预处理的EEMD谐波检测方法[J]. 电力系统保护与控制, 2016, 44(2): 42-48.
Sun Shuguang, Pang Yi, Wang Jingqin, et al. EEMD harmonic detection method based on the new wavelet threshold denoising pretreatment[J]. Power System Protection and Control, 2016, 44(2): 42-48.
[6] 史林军, 周佳佳, 温荣超, 等. 基于经验模态分解与小波分析相结合的风电功率平滑控制[J]. 电力系统保护与控制, 2016, 44(24): 9-16.
Shi Linjun, Zhou Jiajia, Wen Rongchao, et al. Power smoothing control of wind power based on combination of empirical mode decomposition and wavelet analysis[J]. Power System Protection and Control, 2016, 44(24): 9-16.
[7] Liu Zhengkun, Zhang Ze. The improved algorithm of the EMD endpoint effect based on the mirror con- tinuation[C]//2016 Eighth International Conference on Measuring Technology and Mechatronics Auto- mation (ICMTMA), Macau, 2016: 792-795.
[8] 宿文才, 张树团, 刘陵顺. 基于改进极值波延拓的极点对称模态分解端点效应抑制方法[J]. 电工技术学报, 2020, 35(增刊1): 294-301.
Su Wencai, Zhang Shutuan, Liu Lingshun. Extreme- point symmetric mode decomposition endpoint reduction method based on improved extremal wave extension[J]. Transactions of China Electrotechnical Society, 2020, 35(S1): 294-301.
[9] 程启明, 陈路, 程尹曼, 等. 基于EEMD和LS-SVM模型的风电功率短期预测方法[J]. 电力自动化设备, 2018, 38(5): 27-35.
Cheng Qiming, Chen Lu, Cheng Yinman, et al. Short- term wind power forecasting method based on EEMD and LS-SVM model[J]. Electric Power Automation Equipment, 2018, 38(5): 27-35.
[10] 杨剑锋, 石戈戈, 周天奇, 等. 基于自适应序贯相似性检测波形匹配延拓的EMD端点效应抑制[J]. 振动与冲击, 2018, 37(18): 121-125.
Yang Jianfeng, Shi Gege, Zhou Tianqi, et al. Wave- form extension method based on similarity sequential detection for the end effects reduction of EMD[J]. Journal of Vibration and Shock, 2018, 37(18): 121- 125.
[11] 胡明, 郭健鹏, 李富强, 等. 基于自适应CEEMD方法的电能质量扰动检测与分析[J]. 电力系统保护与控制, 2018, 46(21): 103-110.
Hu Ming, Guo Jianpeng, Li Fuqiang, et al. Power quality disturbance detection and analysis based on adaptively complementary ensemble empirical mode decomposition method[J]. Power System Protection and Control, 2018, 46(21): 103-110.
[12] 徐可, 陈宗海, 张陈斌, 等. 基于经验模态分解和支持向量机的滚动轴承故障诊断[J]. 控制理论与应用, 2019, 36(6): 915-922.
Xu Ke, Chen Zonghai, Zhang Chenbin, et al. Rolling bearing fault diagnosis based on empirical mode decomposition and support vector machine[J]. Control Theory & Applications, 2019, 36(6): 915-922.
[13] 李泽, 王辉, 钱勇, 等. 基于加速鲁棒特征的含噪局部放电模式识别[J]. 电工技术学报, 2022, 37(3): 775-785.
Li Ze, Wang Hui, Qian Yong, et al. Pattern recogni- tion of partial discharge in the presence of noise based on speeded up robust features[J]. Transactions of China Electrotechnical Society, 2022, 37(3): 775- 785.
[14] 马星河, 张登奎. 基于改进经验小波变换的高压电缆局部放电噪声抑制研究[J]. 电工技术学报, 2021, 36(增刊1): 353-361.
Ma Xinghe, Zhang Dengkui. Research on suppression of partial discharge noise of high voltage cable based on improved empirical wavelet transform[J]. Transa- ctions of China Electrotechnical Society, 2021, 36(S1): 353-361.
[15] Zhou Yufeng, Yu Mei, Ma Hualin, et al. Weighted-to- spherically-uniform SSIM objective quality evalu- ation for panoramic video[C]//2018 14th IEEE International Conference on Signal Processing (ICSP), Beijing, China, 2018: 54-57.
[16] 肖贤贵, 李开成, 蔡得龙, 等. 一种电能质量扰动信号的联合去噪算法[J]. 电工技术学报, 2021, 36(21): 4418-4428.
Xiao Xiangui, Li Kaicheng, Cai Delong, et al. A combined de-noising method for power quality distur- bances events[J]. Transactions of China Electro- technical Society, 2021, 36(21): 4418-4428.
[17] 中华人民共和国国家质量监督检验检疫总局, 中国国家标准化管理委员会. GB/T 6902—2010 铁路信号继电器试验方法[S]. 北京: 中国标准出版社, 2011.
Application of Empirical Mode Decomposition Algorithm Based on Waveform Matching Endpoint Continuation Method in Noise Reduction of Railway Relay Parameters
Abstract Aiming at the noise problem of railway relay parameters,this paper establishes a hybrid noise reduction model based on empirical mode decomposition algorithm optimized by the improved waveform matching continuation method to extract the effective information.This model uses autocorrelation function to separate noisy signals and wavelet threshold to remove noise. Firstly, according to the characteristics of relay parameters, the waveform matching method is improved.The matching error degree formula is redefined, and the matching accuracy error coefficient is introduced. The end effect produced by EMD decomposition process is optimized by the improved waveform matching continuation method, and the effective IMF component and residual term are obtained. Then, the autocorrelation function is solved, and the noise signal feature is separated according to the autocorrelation function image. Finally, the noisy component is denoised by wavelet threshold. After denoising, it is combined with the residual component and remainder to obtainthe reconstructed parameter sequence. Moreover, the SSIM index, SNR and MSE index are used to evaluate the reliability of the model. Compared with EMD reconstruction and the wavelet threshold de-noising method, the noise reduction effect of the model is proved.
keywords：Endpoint effect, waveform matching, endpoint continuation method, autocorrelation function, structural similarity index measurement (SSIM) index
DOI: 10.19595/j.cnki.1000-6753.tces.201491
中图分类号：TM581
李文华 男，1973年生，教授，研究方向为电器可靠性与寿命预测、储能技术及其可靠性。E-mail: liwenhua@hebut.edu.cn
潘如政 男，1977年生，副教授，研究方向为电器可靠性、高电压技术、放电等离子体应用等。E-mail: prz368@126.com（通信作者）
河北省自然科学基金项目（E2020202221）和河北省自然科学基金创新群体项目（E2020202142）资助。
收稿日期 2020-11-14
改稿日期 2021-01-12
（编辑 崔文静）