CN113887360B - 一种基于迭代扩展频散模态分解的频散波提取方法 - Google Patents

Abstract

本发明涉及一种基于迭代扩展频散模态分解的频散波提取方法，包括：获取频散波时域信号；分别进行快速傅里叶变换和离散短时傅里叶变换，获得频域信号和时频表达式；通过双向快速单脊线提取法，初步估计能量最显著的关键时频脊线；通过迭代扩展频散模态分解，分离重构关键频散模态及其群延迟；重复获取残差信号以及分离重构关键频散模态，直至满足停止迭代准则；对重构的所有关键频散模态进行快速傅里叶变换逆变换，得到频散波时域波形。与现有技术相比，本发明在无需预知频散波模态数量的前提下，通过准确提取和分析信号中的频散模态，高效地对多分量非平稳频散波进行分解，在结构健康监测、水下声学、生物医学、地球物理学等领域有重要应用。

Description

一种基于迭代扩展频散模态分解的频散波提取方法

技术领域

本发明涉及信号处理领域，尤其是涉及一种基于迭代扩展频散模态分解的频散波提取方法。

背景技术

频散信号在结构健康监测、水下声学、生物医学、地球物理学等领域中都有广泛的研究。与时域调频信号不同，频散信号表现出的是与频率相关的特性，其中最重要的是信息是群延迟(Group Delay)。群延迟定义为相位函数对频率的导数，包含众多有关信号源与系统的重要信息，因此可以利用群延迟进行信号源的定位与系统辨识。然而，当宽带信号通过非线性系统时，每个频率都会产生不同的时间延迟，从而导致群延迟在频率上发生变化，从时频图上来看，信号的时频脊线不再是与频率轴平行的直线，而是与频率轴相关的曲线，这使得窄带滤波方法不再适用。因此，准确分离这些非平稳频散波并提取它们的群延迟曲线对于实际应用至关重要。

虽然众多学者和工程人员提出了有关频散信号的提取方法，但是现有的频散信号提取方法存在着算法计算消耗大、时频分辨率低、计算时间长等问题。目前，国内有关频散信号提取技术的发明专利申请有：

授权日为2021年8月6日，授权公告号为CN110045014B的中国专利中，公开了一种名称为“基于贝叶斯学习的Lamb波频散消除方法及其系统”。其通过构造多模态频散传播字典和非频散传播字典，将频散多模态Lamb波信号进行稀疏表示，并利用贝叶斯学习方法对字典系数进行系数求解，从而得到频散波的稀疏表达并进行频散消除。然而，该方法的计算结果与字典构建的匹配性有非常强的依赖，同时字典系数的稀疏化求解需要耗费大量的计算时间和计算资源消耗。

公开日为2021年5月28日，公开号为CN112858481A的中国专利中，公开了一种名称为“一种基于时间频率域评估板材中裂纹深度的方法”。其采用时频分析方法求解Lamb波频散曲线。然而传统的非参数化时频分析方法，在群延时快速波动时，会造成分辨率降低的现象；参数化时频分析方法虽然能够提高时频分辨率，但是计算复杂度大大提高，给实际工程应用带来巨大困难。

发明内容

本发明的目的就是为了克服上述现有技术存在计算消耗大、时频分辨率低、计算时间长的缺陷而提供一种基于迭代扩展频散模态分解的频散波提取方法。

本发明的目的可以通过以下技术方案来实现：

一种基于迭代扩展频散模态分解的频散波提取方法，包括以下步骤：

步骤1：获取频散波时域信号；

步骤2：对所述频散波时域信号进行快速傅里叶变换，获得对应的频域信号；

步骤3：对所述频散波时域信号进行离散短时傅里叶变换，获得时频表达式；

步骤4：初始化预设的停止迭代准则，设定频散波信号成分数量m＝0；

步骤5：通过双向快速单脊线提取法，初步估计所述时频表达式中能量最显著的关键时频脊线；

步骤6：通过迭代扩展频散模态分解，根据所述频域信号和关键时频脊线的初步估计值，分离重构关键频散模态及其群延迟；

步骤7：根据重构的关键频散模态，从频域信号中得到残差信号，对残差信号重复依次执行步骤5至步骤7，重复获取残差信号以及进行关键频散模态的分离重构，直至满足所述停止迭代准则；

步骤8：对重构的所有关键频散模态进行快速傅里叶变换逆变换，得到频散波的时域波形。

进一步地，步骤2具体包括以下步骤：

S201：对频散波时域信号S(t)进行快速傅里叶变换，得到对应的频域信号S(f)；

S202：通过频散波时域信号的长度N_t，确定频域信号长度N_f，该频域信号长度N_f的计算表达式为：

N_f＝floor(N_t/2)+1

式中，floor(·)为向下取整符号；

S203：保留频域信号的前N_f项，第N_f+1项至第N_t项置零。

进一步地，步骤3具体包括以下步骤：

S301：设置窗长度N_win，确定傅里叶点数为2N_f-1，构造相应窗函数W(t)；

S302：通过所述窗函数W(t)对时域频散波信号S(t)进行离散短时傅里叶变换，获得时频表达式TFR(t,f)。

进一步地，所述步骤5包括以下步骤：

步骤501：进行参数初始化；

步骤502：定位所述时频表达式TFR(t,f)中能量最大值点(t_em,f_em)并沿频率方向双向搜索Δt范围内时间维度的能量极大值点；

步骤503：获得时频脊线作为关键模态群延时的初始估计值

进一步地，步骤501中，所述参数为最大时间允许波动范围Δt；

步骤502中，所述能量最大值点(t_em,f_em)的计算表达式为：

式中，t为时间，f为幅值；

所述沿频率方向双向搜索Δt范围内时间维度的能量极大值点的计算表达式为：

t_R＝t_em,t_L＝t_em

式中，N_f为频域信号长度。

进一步地，所述步骤6具体包括以下步骤：

步骤601：初始化参数，获取频域频散波信号S(f)和关键模态群延时初始值

步骤602：构造关键频散模态并建立解析幅值/>的冗余傅里叶展开式；

步骤603：将所述的冗余傅里叶展开式带入关键频散模态/>中构造矩阵表达式，通过求解最优化问题获得傅里叶参数矩阵/>并重构得/>

步骤604：利用所述重构幅值通过解卷绕求导方法估计群延时误差/>对其进行滤波光滑然后更新群延时/>

步骤605：将代入步骤603更新核矩阵K^k+1、参数矩阵/>和解析幅值矩阵重复计算达到预设的迭代次数后停止，然后求解重构关键频散模态s_key(f)；

步骤606：更新频散波信号成分数量m＝m+1。

进一步地，步骤601中，所述参数包括正则化参数λ，惩罚参数β,μ和傅里叶展开阶数L；

步骤602中，所述关键频散模态的计算表达式为：

其中，k＝0,...,end为迭代次数，为解析幅值：

其中，为初始相位；

解析幅值的冗余傅里叶展开式为：

其中，f_s为时域信号的采样频率；

步骤603中，构造的矩阵表达式为：

K^k＝E^k·T为N_f×(2L+1)的核矩阵，其中，

所述最优化问题为：

获得的所述傅里叶参数矩阵的计算表达式为：

所述重构幅值的计算表达式为：

式中，I为单位矩阵，(·)^H为矩阵的共轭转置；

步骤604中，取的相位值并进行解卷绕，通过对解卷绕相位进行求导得到群延时误差/>

式中，∠·为求复数的相位运算，unwrap(·)为相位解卷绕运算，将相位值由[-π,π]区间叠加到[0,+∞)区间；

对进行滤波光滑并更新瞬时频率/>该步骤的计算表达式为：

其中为改进二阶差分算子；

步骤605中，所述关键频散模态s_key(f)的重构表达式为：

式中，end为预设的迭代次数，K^end为更新end次后的核矩阵，为更新end次后的参数矩阵。

进一步地，所述步骤7包括以下步骤：

步骤701：将所述重构关键频散模态s_key(f)作为第m个频散模态分量s_m(f)，并将其从频散波频域信号S(f)中提取，得到残差信号S_r(f)：

S_r(f)＝S(f)-s_m(f)

步骤702：计算残差信号能量与频散波频域信号能量之比：

步骤703：当所述能量比R大于预先设定的迭代停止判别值ε时，重复依次执行步骤5、步骤6和步骤7；否则，得到频散模态分量总数M＝m，执行步骤8。

进一步地，所述步骤8包括以下步骤：

步骤801：扩展所有重构的频散模态分量s_m(f)；

步骤802：对所有扩展后的频散模态分量进行快速傅里叶变换逆变换，得到频散波的时域波形。

进一步地，步骤801中，扩展后的频散模态分量的计算表达式为：

式中，conj[·]表示求复数的共轭；

步骤802中，取快速傅里叶变换逆变换计算结果的实部作为对应的频散波时域波形该频散波时域波形/>的计算表达式为：

式中，IFFT(·)为快速傅里叶变换逆变换操作，real[·]为取复数实部操作。

本发明新提出了一种迭代扩展频散模态分解法(IEDMD)，通过双向快速单脊线提取法精确提取频散模态时频脊线，实现了精确高分辨率的非稳态频散波提取。与现有技术相比，本发明具有以下优点：

(1)本发明相比同类方法，群延迟提取精度提高；

(2)本发明无需预知频散波模态的数量，通过能量比方法自适应地确定提取数量，更有利于实际工程的应用；

(3)本发明在强噪声背景下精确重构频散波的时域波形和群延时，有效地抑制了原信号中的噪声干扰和成分混叠，信息丰富完整，有利于对其进行深入分析；

(4)本发明计算速度相比同类方法提升2倍以上，瞬时频率提取精度提高20％以上，模态重构精度提高15％以上，特征阶次误差小于1％并且计算精度提高20％以上，无需较多人为干预，有利于嵌入式开发和实际工程应用。

附图说明

图1为本发明实施例中提供一种基于迭代扩展频散模态分解的频散波提取方法的流程示意图；

图2(a)为本发明实施例中频散波仿真信号时域图；

图2(b)为本发明实施例中频散波仿真模态1时域图；

图2(c)为本发明实施例中频散波仿真模态2时域图；

图3为本发明实施例中频散波的时频图；

图4(a)为本发明实施例中第一频散模态IEDMD群延迟估计图；

图4(b)为本发明实施例中第一频散模态IEDMD群延时计算误差曲线图；

图4(c)为本发明实施例中第一频散模态IEDMD时域信号及计算误差图；

图4(d)为本发明实施例中第一频散模态GDMD时域信号及计算误差图；

图5(a)为本发明实施例中第二频散模态IEDMD群延迟估计图；

图5(b)为本发明实施例中第二频散模态IEDMD群延时计算误差曲线图；

图5(c)为本发明实施例中第二频散模态IEDMD时域信号及计算误差图；

图5(d)为本发明实施例中第二频散模态GDMD时域信号及计算误差图；

图中，Time为时间，Amplitude为振幅，Freq为频率，Frequency为频率，Dispersivemode为频散模态，Residual为残差信号，True为真实信号。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。

因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步定义和解释。

实施例1

如图1所示，本实施例提供一种基于迭代扩展频散模态分解的频散波提取方法，在无需预知频散波模态数量的前提下，通过准确提取和分析信号中的频散模态，高效地对多分量非平稳频散波进行分解，在结构健康监测、水下声学、生物医学、地球物理学等领域有重要应用。

假设采集到的频散波时域信号为S(t)，该方法具体包含以下步骤：

步骤1：根据实际应用场景，按需要安装传感器，获取相关领域的频散波时域信号S(t)；

步骤2，对所述时域信号进行快速傅里叶变换，获得频散波对应的频域信号，具体包括以下步骤：

步骤201：对时域频散波信号S(t)进行快速傅里叶变换(FFT)，得到对应的频域信号S(f)；

步骤202：通过时域信号长度N_t，确定频域信号长度N_f：

N_f＝floor(N_t/2)+1

式中，floor(·)为向下取整符号；

步骤203：保留频域信号的前N_f项，第N_f+1项至第N_t项置零。

步骤3，对所述时域信号进行离散短时傅里叶变换，获得时频表达式，具体包括以下步骤：

步骤301：设置窗长度N_win，确定傅里叶点数为2N_f-1，构造相应窗函数W(t)；

步骤302：通过所述窗函数W(t)对时域频散波信号S(t)进行离散短时傅里叶变换(STFT)，获得时频表达式TFR(t,f)。

步骤4，设定迭代停止判别值ε，设定频散波信号成分数量m＝0。

步骤5，通过双向快速单脊线提取法，初步估计所述时频表达式中能量最显著的关键时频脊线，具体包括以下步骤：

步骤501：初始化参数最大时间允许波动范围Δt；

步骤502：定位所述时频表达式TFR(t,f)中能量最大值点(t_em,f_em)：

沿频率方向双向搜索Δt范围内时间维度的能量极大值点：

t_R＝t_em,t_L＝t_em

步骤503：获得时频脊线作为关键模态群延时的初始值。

步骤6，通过迭代扩展频散模态分解，根据所述频域信号和关键时频脊线初始值，分离重构关键频散模态及其群延迟，具体包括以下步骤：

步骤601：初始化正则化参数λ，惩罚参数β,μ和傅里叶展开阶数L，输入频域频散波信号S(f)和关键模态群延时初始值

步骤602：构造关键频散模态

其中，k＝0,...,end为迭代次数，为解析幅值：

其中，为初始相位；然后建立解析幅值/>的冗余傅里叶展开式：

其中，f_s为时域信号的采样频率；

步骤603：将所述展开模型带入/>中构造矩阵表达式：

K^k＝E^k·T为N_f×(2L+1)的核矩阵：

然后求解最优化问题：

获得傅里叶参数矩阵并重构得/>

其中，I为单位矩阵，(·)^H为矩阵的共轭转置；

步骤604：取的相位值并进行解卷绕，通过对解卷绕相位进行求导得到群延时误差

式中，∠·为求复数的相位运算，unwrap(·)为相位解卷绕运算，即将相位值由[-π,π]区间叠加到[0,+∞)区间；

对进行滤波光滑并更新瞬时频率/>

其中为改进二阶差分算子；

步骤605：将所述跳转至步骤603更新核矩阵K^k+1、参数矩阵/>和解析幅值矩阵/>重复计算end次后停止，然后求解重构关键频散模态s_key(f)：

步骤606：更新频散波信号成分数量m＝m+1；

步骤7，从所述频散波频域信号中提取出重构的关键频散模态分量，对残余信号重复进行步骤5和步骤6，直到达到停止迭代准则，具体包括以下步骤：

步骤701：将所述重构关键频散模态s_key(f)作为第m个频散模态分量s_m(f)，并将其从频散波频域信号S(f)中提取，得到残差信号S_r(f)：

S_r(f)＝S(f)-s_m(f)

步骤702：计算残差信号能量与频散波频域信号能量之比：

步骤703：当所述能量比R大于预先设定的迭代停止判别值ε时，重复执行步骤5、步骤6和步骤7；否则，得到频散模态分量总数M＝m，执行步骤8。

步骤8，对所述所有重构关键频散模态进行快速傅里叶变换逆变换，得到频散波的时域波形，具体包括以下步骤：

步骤801：将所有重构频散模态分量s_m(f)在频域上进行扩展，得到扩展频散模态分量

式中，conj[·]表示求复数的共轭；

步骤802：对所述的扩展频散模态分量逐一进行快速傅里叶变换逆变换(IFFT)，并取IFFT计算结果的实部作为对应的频散波时域波形/>

式中，IFFT(·)为快速傅里叶变换逆变换操作，real[·]为取复数实部操作。

本实施例的具体实施结果如下：

如图2(a)所示，频散波仿真信号S(t)的信噪比为0dB，采样频率100Hz，采样时间15s；如图2(b)和图2(c)所示，频散波仿真信号S(t)中存在2个频散模态。如图3所示，两个频散模态在时间上和频率上都存在重叠。采用双向快速单脊线提取法，初步估计图3中能量最强的时频脊线，并利用IEDMD算法重构，可以得到精确的频散波模态分量及其群延迟。如图4(a)和图4(b)所示，IEDMD提取的第一频散模态群延迟与真实群延迟匹配度极高，与同类方法相比，群延迟的提取精度有显著提高，IEDMD方法的整体群延迟估计误差在1％以内，而GDMD方法在信号端点的群延时估计误差达到了15％；如图4(c)和图4(d)所示，IEDMD提取的第一频散模态时域信号与真实信号的误差极小，而GDMD提取的第一频散模态时域信号在信号端点产生较大误差。采用平均相对误差(RE_mean)指标，评价IEDMD对第一模态分量群延迟估计的准确性，如表1所示。与同类方法相比，第一模态分量群延迟估计的相对误差降低了36％。

表1IEDMD第一模态分量群延迟提取精度比较

将第一频散模态估计值从原始信号中提取，并对残差信号继续采用IEDMD方法提取第二模态分量。如图5(a)和图5(b)所示，IEDMD提取的第二频散模态群延迟与真实群延迟匹配度极高，与同类方法相比，群延迟的提取精度有显著提高，IEDMD方法的整体群延迟估计误差在0.2％以内，而GDMD方法在信号端点的群延时估计误差达到了50％；如图5(c)和图5(d)所示，IEDMD提取的第二频散模态时域信号与真实信号的误差极小，而GDMD提取的第二频散模态时域信号在信号端点产生较大误差。采用平均相对误差(RE_mean)指标，评价IEDMD对第二模态分量群延迟估计的准确性，如表2所示。与同类方法相比，第二模态分量群延迟估计的相对误差降低了95％。

表2IEDMD第二模态分量群延迟提取精度比较

采用均方根误差(RMSE)、信噪比(SNR)和相关系数(CC)对IEDMD方法的信号提取准确度进行评价，如表3所示。与同类方法相比，RMSE指标降低了32％，SNR指标提高了9.3％，CC指标提高了3.5％，表明IEDMD方法在大幅度降低频散波信号估计误差的同时，提取信号与原信号存在更高的相关度，并抑制了噪声。

表3IEDMD信号提取精度比较

综上所述，与同类方法相比，本方法的各项性能指标均有可观的提升。

以上详细描述了本发明的较佳具体实施例。应当理解，本领域的普通技术人员无需创造性劳动就可以根据本发明的构思做出诸多修改和变化。因此，凡本技术领域中技术人员依本发明的构思在现有技术的基础上通过逻辑分析、推理或者有限的实验可以得到的技术方案，皆应在由权利要求书所确定的保护范围内。

Claims (4)

1.一种基于迭代扩展频散模态分解的频散波提取方法，其特征在于，包括以下步骤：

步骤1：获取频散波时域信号；

步骤2：对所述频散波时域信号进行快速傅里叶变换，获得对应的频域信号；

步骤3：对所述频散波时域信号进行离散短时傅里叶变换，获得时频表达式；

步骤4：初始化预设的停止迭代准则，设定频散波信号成分数量m＝0；

步骤5：通过双向快速单脊线提取法，初步估计所述时频表达式中能量最显著的关键时频脊线；

步骤6：通过迭代扩展频散模态分解，根据所述频域信号和关键时频脊线的初步估计值，分离重构关键频散模态及其群延迟；

步骤7：根据重构的关键频散模态，从频域信号中得到残差信号，对残差信号重复依次执行步骤5至步骤7，重复获取残差信号以及进行关键频散模态的分离重构，直至满足所述停止迭代准则；

步骤8：对重构的所有关键频散模态进行快速傅里叶变换逆变换，得到频散波的时域波形；

步骤2具体包括以下步骤：

S201：对频散波时域信号S(t)进行快速傅里叶变换，得到对应的频域信号S(f)；

S202：通过频散波时域信号的长度N_t，确定频域信号长度N_f，该频域信号长度N_f的计算表达式为：

N_f＝floor(N_t/2)+1

式中，floor(·)为向下取整符号；

S203：保留频域信号的前N_f项，第N_f+1项至第N_t项置零；

步骤3具体包括以下步骤：

S301：设置窗长度N_win，确定傅里叶点数为2N_f-1，构造相应窗函数W(t)；

S302：通过所述窗函数W(t)对时域频散波信号S(t)进行离散短时傅里叶变换，获得时频表达式TFR(t,f)；

所述步骤5包括以下步骤：

步骤501：进行参数初始化；

步骤502：定位所述时频表达式TFR(t,f)中能量最大值点(t_em,f_em)并沿频率方向双向搜索Δt范围内时间维度的能量极大值点；

步骤503：获得时频脊线作为关键模态群延时的初始估计值GD_k0_ey；

步骤501中，所述参数为最大时间允许波动范围Δt；

步骤502中，所述能量最大值点(t_em,f_em)的计算表达式为：

式中，t为时间，f为幅值；

所述沿频率方向双向搜索Δt范围内时间维度的能量极大值点的计算表达式为：

t_R＝t_em,t_L＝t_em

式中，N_f为频域信号长度；

所述步骤6具体包括以下步骤：

步骤601：初始化参数，获取频域频散波信号S(f)和关键模态群延时初始值

步骤602：构造关键频散模态并建立解析幅值/>的冗余傅里叶展开式；

步骤603：将所述的冗余傅里叶展开式带入关键频散模态/>中构造矩阵表达式，通过求解最优化问题获得傅里叶参数矩阵/>并重构得/>

步骤604：利用重构幅值通过解卷绕求导方法估计群延时误差/>对其进行滤波光滑然后更新群延时/>

步骤605：将代入步骤603更新核矩阵K^k+1、参数矩阵/>和解析幅值矩阵/>重复计算达到预设的迭代次数后停止，然后求解重构关键频散模态s_key(f)；

步骤606：更新频散波信号成分数量m＝m+1；

步骤601中，所述参数包括正则化参数λ，惩罚参数β,μ和傅里叶展开阶数L；

步骤602中，所述关键频散模态的计算表达式为：

其中，k＝0,...,end为迭代次数，为解析幅值：

其中，为初始相位；

解析幅值的冗余傅里叶展开式为：

其中，f_s为时域信号的采样频率；

步骤603中，构造的矩阵表达式为：

K^k＝E^k·T为N_f×(2L+1)的核矩阵，其中，

所述最优化问题为：

获得的所述傅里叶参数矩阵的计算表达式为：

重构幅值的计算表达式为：

式中，I为单位矩阵，(·)^H为矩阵的共轭转置；

步骤604中，取的相位值并进行解卷绕，通过对解卷绕相位进行求导得到群延时误差/>

式中，∠γ为求复数的相位运算，unwrap(·)为相位解卷绕运算，将相位值由[-π,π]区间叠加到[0,+∞)区间；

对进行滤波光滑并更新瞬时频率/>该步骤的计算表达式为：

其中为改进二阶差分算子；

步骤605中，所述关键频散模态s_key(f)的重构表达式为：

式中，end为预设的迭代次数，K^end为更新end次后的核矩阵，为更新end次后的参数矩阵。

2.根据权利要求1所述的一种基于迭代扩展频散模态分解的频散波提取方法，其特征在于，所述步骤7包括以下步骤：

步骤701：将所述重构关键频散模态s_key(f)作为第m个频散模态分量s_m(f)，并将其从频散波频域信号S(f)中提取，得到残差信号S_r(f)：

S_r(f)＝S(f)-s_m(f)

步骤702：计算残差信号能量与频散波频域信号能量之比：

步骤703：当所述能量比R大于预先设定的迭代停止判别值ε时，重复依次执行步骤5、步骤6和步骤7；否则，得到频散模态分量总数M＝m，执行步骤8。

3.根据权利要求1所述的一种基于迭代扩展频散模态分解的频散波提取方法，其特征在于，所述步骤8包括以下步骤：

步骤801：扩展所有重构的频散模态分量s_m(f)；

步骤802：对所有扩展后的频散模态分量进行快速傅里叶变换逆变换，得到频散波的时域波形。

4.根据权利要求3所述的一种基于迭代扩展频散模态分解的频散波提取方法，其特征在于，步骤801中，扩展后的频散模态分量的计算表达式为：

式中，conj[·]表示求复数的共轭；

步骤802中，取快速傅里叶变换逆变换计算结果的实部作为对应的频散波时域波形该频散波时域波形/>的计算表达式为：

式中，IFFT(·)为快速傅里叶变换逆变换操作，real[·]为取复数实部操作。