CN117476039B - 基于声学信号的水轮机初生空化预警方法 - Google Patents

基于声学信号的水轮机初生空化预警方法 Download PDF

Info

Publication number
CN117476039B
CN117476039B CN202311787284.9A CN202311787284A CN117476039B CN 117476039 B CN117476039 B CN 117476039B CN 202311787284 A CN202311787284 A CN 202311787284A CN 117476039 B CN117476039 B CN 117476039B
Authority
CN
China
Prior art keywords
cavitation
acoustic
early warning
signal
acoustic signal
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202311787284.9A
Other languages
English (en)
Other versions
CN117476039A (zh
Inventor
罗兴锜
赵楠楠
冯建军
朱国俊
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xian University of Technology
Original Assignee
Xian University of Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Xian University of Technology filed Critical Xian University of Technology
Priority to CN202311787284.9A priority Critical patent/CN117476039B/zh
Publication of CN117476039A publication Critical patent/CN117476039A/zh
Application granted granted Critical
Publication of CN117476039B publication Critical patent/CN117476039B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G10MUSICAL INSTRUMENTS; ACOUSTICS
    • G10LSPEECH ANALYSIS TECHNIQUES OR SPEECH SYNTHESIS; SPEECH RECOGNITION; SPEECH OR VOICE PROCESSING TECHNIQUES; SPEECH OR AUDIO CODING OR DECODING
    • G10L25/00Speech or voice analysis techniques not restricted to a single one of groups G10L15/00 - G10L21/00
    • G10L25/48Speech or voice analysis techniques not restricted to a single one of groups G10L15/00 - G10L21/00 specially adapted for particular use
    • G10L25/51Speech or voice analysis techniques not restricted to a single one of groups G10L15/00 - G10L21/00 specially adapted for particular use for comparison or discrimination
    • FMECHANICAL ENGINEERING; LIGHTING; HEATING; WEAPONS; BLASTING
    • F03MACHINES OR ENGINES FOR LIQUIDS; WIND, SPRING, OR WEIGHT MOTORS; PRODUCING MECHANICAL POWER OR A REACTIVE PROPULSIVE THRUST, NOT OTHERWISE PROVIDED FOR
    • F03BMACHINES OR ENGINES FOR LIQUIDS
    • F03B11/00Parts or details not provided for in, or of interest apart from, the preceding groups, e.g. wear-protection couplings, between turbine and generator
    • F03B11/008Measuring or testing arrangements
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01HMEASUREMENT OF MECHANICAL VIBRATIONS OR ULTRASONIC, SONIC OR INFRASONIC WAVES
    • G01H17/00Measuring mechanical vibrations or ultrasonic, sonic or infrasonic waves, not provided for in the preceding groups
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00Pattern recognition
    • G06F18/10Pre-processing; Data cleansing
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00Pattern recognition
    • G06F18/20Analysing
    • G06F18/21Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
    • G06F18/214Generating training patterns; Bootstrap methods, e.g. bagging or boosting
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00Pattern recognition
    • G06F18/20Analysing
    • G06F18/21Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
    • G06F18/217Validation; Performance evaluation; Active pattern learning techniques
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06NCOMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00Computing arrangements based on biological models
    • G06N3/004Artificial life, i.e. computing arrangements simulating life
    • G06N3/006Artificial life, i.e. computing arrangements simulating life based on simulated virtual individual or collective life forms, e.g. social simulations or particle swarm optimisation [PSO]
    • GPHYSICS
    • G10MUSICAL INSTRUMENTS; ACOUSTICS
    • G10LSPEECH ANALYSIS TECHNIQUES OR SPEECH SYNTHESIS; SPEECH RECOGNITION; SPEECH OR VOICE PROCESSING TECHNIQUES; SPEECH OR AUDIO CODING OR DECODING
    • G10L21/00Speech or voice signal processing techniques to produce another audible or non-audible signal, e.g. visual or tactile, in order to modify its quality or its intelligibility
    • G10L21/02Speech enhancement, e.g. noise reduction or echo cancellation
    • G10L21/0208Noise filtering
    • G10L21/0216Noise filtering characterised by the method used for estimating noise
    • GPHYSICS
    • G10MUSICAL INSTRUMENTS; ACOUSTICS
    • G10LSPEECH ANALYSIS TECHNIQUES OR SPEECH SYNTHESIS; SPEECH RECOGNITION; SPEECH OR VOICE PROCESSING TECHNIQUES; SPEECH OR AUDIO CODING OR DECODING
    • G10L25/00Speech or voice analysis techniques not restricted to a single one of groups G10L15/00 - G10L21/00
    • G10L25/27Speech or voice analysis techniques not restricted to a single one of groups G10L15/00 - G10L21/00 characterised by the analysis technique
    • G10L25/30Speech or voice analysis techniques not restricted to a single one of groups G10L15/00 - G10L21/00 characterised by the analysis technique using neural networks
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2123/00Data types
    • G06F2123/02Data types in the time domain, e.g. time-series data
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/02Preprocessing
    • G06F2218/04Denoising
    • YGENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02EREDUCTION OF GREENHOUSE GAS [GHG] EMISSIONS, RELATED TO ENERGY GENERATION, TRANSMISSION OR DISTRIBUTION
    • Y02E10/00Energy generation through renewable energy sources
    • Y02E10/20Hydro energy

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • General Physics & Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Artificial Intelligence (AREA)
  • Evolutionary Computation (AREA)
  • Computational Linguistics (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Multimedia (AREA)
  • Evolutionary Biology (AREA)
  • Signal Processing (AREA)
  • Audiology, Speech & Language Pathology (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Acoustics & Sound (AREA)
  • Human Computer Interaction (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Mechanical Engineering (AREA)
  • Combustion & Propulsion (AREA)
  • Quality & Reliability (AREA)
  • Chemical & Material Sciences (AREA)
  • Biomedical Technology (AREA)
  • Biophysics (AREA)
  • General Health & Medical Sciences (AREA)
  • Molecular Biology (AREA)
  • Computing Systems (AREA)
  • Mathematical Physics (AREA)
  • Software Systems (AREA)
  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)

Abstract

本发明公开了基于声学信号的水轮机初生空化预警方法,采集各工况运行参数下不同空化状态水轮机运行时的声学信号,对声学信号去噪处理,计算低频空化声学信号熵率和高频空化声学信号的信号瞬时能量值;构建水轮机工况运行参数与初生空化声学信号预警特征值之间关系的预测模型,通过水轮机工况运行参数空化声学信号熵率和高频声学信号的信号能量值对预测模型训练;采集工况运行参数已知、空化状态未知状态的实况声学信号,计算实况声学信号的熵率和能量值,将实况声学信号工况运行参数输入训练后的预测模型中获得熵率预警阈值与能量值的预警阈值;当计算信号的熵率或能量值指标超过相应预警阈值时,发出预警,精确实现水轮机初生空化预警。

Description

基于声学信号的水轮机初生空化预警方法
技术领域
本发明属于水轮机技术领域,具体涉及基于声学信号的水轮机初生空化预警方法。
背景技术
能源已成为当今世界各国经济发展最重要的制约因素和国际竞争的焦点问题,保护生态环境,发展清洁能源,改善电网稳定性及质量是未来电力工业新的发展方向。水电是公认的绿色环保能源,对于调节电网稳定性具有重要保障作用。作为水电系统中的核心换能部件,其稳定性对整个系统乃至电网的安全运行至关重要,但因其工况条件多变,偏工况下水轮机内部流态恶化,压强降低诱发空化空蚀,直接关系到水轮机的效率以及使用寿命。因此识别机组空化状态对于保障机组运行安全具有重要意义。
由于空化问题涉及到多个学科,耦合性较强,一直是人们研究的热点和难点。一般根据空化发展程度将水轮机运行状态分为无空化、初生空化、临界空化以及完全空化状态。而由于初生空化状态时产生的空泡数较临界空化和完全空化时较少,使得难以及时识别水轮机初生空化状态。目前尚无实际应用的在线监测设备,传统方法常通过观察转轮叶片上空泡对空化状态进行判别,这种人工目测观察法对观察者能力要求较高,效率较低且因观察位置不同可能会导致误判。现有技术常通过对空化信号提取简单的时频域特征并结合机器学习方法对水轮机空化状态进行判断,并没有结合水轮机空化原理,在空化现象识别精度及效率上有所欠缺。
发明内容
本发明的目的是提供基于声学信号的水轮机初生空化预警方法,采集低频和高频声学信号,利用改进降噪方法对所采集的声学信号进行降噪处理,基于空化原理分析空化状态与声学信号之间的关联特性,提出更有效的空化特征指标,可以有效解决强环境噪声干扰下空化识别不准确的问题。
本发明所采用的技术方案是,基于声学信号的水轮机初生空化预警方法,具体按照以下步骤实施:
步骤1、通过水轮机上安装的声学传感器获取空化声学信号数据集,包含各工况运行参数下无空化及空化状态下水轮机运行时的声学信号;
步骤2、对声学信号进行去噪处理,得到去噪声学信号;具体过程为:
步骤2.1、输入所采集的声学信号对应的时间序列,确定改进自适应噪声完全集合经验模态分解的最优白噪声幅值£和加噪次数Ne
步骤2.2、根据最优白噪声幅值£和加噪次数Ne对采集的声学信号进行分解,获得K个IMF分量;
步骤2.3、分别计算K个IMF分量与分解前的声学信号之间的最大互信息系数,根据最大互信息系数对IMF分量筛选、重构,得到去噪声学信号;
步骤3、计算低频空化声学信号熵率和高频空化声学信号的信号能量值;
步骤4、构建水轮机工况运行参数与初生空化声学信号预警特征值之间关系的预测模型,将水轮机工况运行参数作为输入数据,初生空化声学信号熵率和高频声学信号的信号能量值作为输出数据,每组输入数据及对应的输出数据作为训练样本训练预测模型;
步骤5、采集工况运行参数已知、空化状态未知状态的实况声学信号,计算实况声学信号的熵率和能量值,将实况声学信号对应的工况运行参数输入训练后的预测模型中获得熵率预警阈值与能量值的预警阈值;当实况声学信号的熵率超过熵率预警阈值或实况声学信号的能量值超过能量值的预警阈值,发出初生空化预警。
本发明的特点还在于:
工况运行参数包括导叶开度α0,转速n,水头H,下游水位Hd、水轮机出力P和扭矩M。
步骤1具体过程为:在尾水管直锥段接近转轮出口处安装水听器与声发射传感器,其中水听器采样频率设置为60kHz,声发射传感器采样频率设置为2MHz;固定水轮机的运行参数,并在不同的空化状态下利用水听器与声发射传感器同步采集水轮机空化声学信号,通过数据采集卡发送到计算机并记录声学信号对应时间,直至当转轮叶片上出现空泡时,将此状态判断为初生空化状态,通过水听器与声发射传感器采集声学信号,当一个工况的空化声学信号采集完成后,分别改变工况运行参数,获得不同工况下声学信号。
步骤2.1具体过程为:
步骤2.1.1、以白噪声幅值和加噪次数为粒子坐标,设定好坐标参数搜索范围,初始化粒子群算法参数;
步骤2.1.2、在当前白噪声幅值和加噪次数下,对所采集的声学信号序列进行自适应噪声完全集合经验模态分解,计算各IMF分量的分形维数,将分形维数的和作为适应度函数;将初始粒子适应度作为个体最优值,粒子历史最优适应度值作为全局最优值;
步骤2.1.3、对粒子位置以及速度、粒子种群最优解进行更新,更新公式为:
式中,是个体历史最优位置;/>是全局历史最优位置;/>是目前粒子位置,/>是下一次迭代的粒子速度即更新后的速度;w代表惯性权重,/>,/>分别取0.9和0.4,c 1c 2为学习因子,/>,/>是(0,1)之间的随机数,下标ij表示第i个粒子和第j维,t代表迭代次数,/>表示最大进化代数,计算更新后的粒子适应度值,将当前粒子种群适应度值与个体最优值、全局最优值进行比较,更新粒子种群的个体最优以及全局最优极值;
步骤2.1.4、判断是否满足迭代条件,设置容忍度阈值、最大容忍代数,最大迭代次数,计算当前最优适应度与上一代最优适应度的相对变化量作为容忍度,比较容忍度与容忍度阈值大小,若前者小则容忍代数加1,判断容忍代数是否大于最大容忍代数或者迭代次数是否超过最大迭代次数,若满足,则输出目前全局最优解对应的白噪声幅值和加噪次数,作为最优白噪声幅值£和加噪次数Ne,若不满足,则返回步骤2.1.3。
计算各IMF分量的分形维数过程为:
设定IMF分量序列C的时间序列为x(t),t=1,2,NN表示整个时间序列的长度;构造新的以延迟时间k为参数的时间序列/>表示为:
将时间序列展开后计算时间序列矩阵的曲线长度Lm(k),表示为:
式中,N表示整个时间序列的长度,是归一化因子,[/>]表示对/>取整,上式由k=1到/>循环计算,求得一个k值下所有L m(k)的均值L(k),在双对数坐标系下利用最小二乘法拟合/>与/>得到直线斜率即是IMF分量的分形维数。
步骤2.3具体过程为:
计算IMF分量与分解前的声学信号之间的最大互信息系数:
其中,XY分别为IMF分量和分解前的声学信号,表示XY之间的互信息,表示变量间的联合概率密度,/>为变量间的最大互信息系数,B表示网格划分系数;
选择最大互信息系数前5阶IMF分量重构声学信号,得到去噪声学信号。
步骤3具体过程为:
步骤3.1、将去噪声学信号进行粗粒化,粗粒化过程中不同的起点创建不同的时间序列,重构后声学序列的粗粒化序列为:
其中,L为重构声学序列长度;τ为尺度因子,u b表示重构声学序列的第b个值;
步骤3.2、计算粗粒化时间序列中色散模式π的平均概率
其中,色散模式π为重构空间中序列的组合形式,其中为不同起始点每个可能的色散模式/>的概率;
步骤3.3、计算色散模式概率平均值的色散熵:
其中,c为类别个数,m为嵌入维数,d为时延长度;
步骤3.4、利用色散熵的前八个尺度的熵值进行最小二乘法拟合,得到的曲线斜率即为低频声学信号的熵率;
步骤3.5、由下式计算重构后的高频声学信号的Teager能量算子序列,重构公式为:
式中,表示Teager能量算子序列的第n个值,/>表示重构高频信号序列中的第n个数,/>表示重构高频信号序列中的第n+1个数,/>表示重构高频信号序列中的第n-1个数;
对信号Teager能量算子序列进行快速傅里叶变换,其主频幅值即为高频声学信号的信号能量值。
步骤3.2中不同起始点每个可能的色散模式的概率计算过程为:
对于每个粗粒化时间序列,通过下式将变量映射到{1,2,/>,c}各类别中:
式中,为序列的期望,/>为序列的标准差,round为取整函数,c为类别个数;/>表示/>通过正态累积分布函数变换后的累积概率值;
将嵌入向量映射/>种基于波动的色散模式/>中,则有:
式中m为嵌入维数,d为时延长度。
每个可能的色散模式的概率由下式计算:
式中,count是指映射至/>的模式个数,再除以嵌入维数为m对应的嵌入向量总数即为色散模式的概率/>N表示粗粒化序列长度。
步骤4中水轮机工况运行参数与初生空化声学信号预警特征值之间关系的预测模型为BiGRU深度学习模型,将水轮机每个工况运行参数作为输入数据,每个工况运行参数对应的初生空化声学信号熵率和高频声学信号的信号能量值作为输出数据,每个输入数据及对应的输出数据作为训练样本,对训练样本进行归一化处理后输入BiGRU深度学习模型,采用十折交叉验证与网格搜索优化预测模型的超参数、BiGRU层数、神经元个数以及学习率,当训练次数达到设定最大次数或模型误差小于1%时保存BiGRU深度学习模型,作为训练后的预测模型。
本发明有益效果是:
本发明基于声学信号的水轮机初生空化预警方法,基于改进降噪方法去除环境噪声干扰,通过分析声学信号随水轮机运行状态的变化规律,采用Teager能量算子为提取旋转机械空化信号中的瞬态冲击特征提供了一种有效途径,可以有效增强信号的瞬态特征,适用于处理信号中的空泡溃灭冲击特征,即为瞬时能量理论。本发明基于信号复杂度与瞬时能量理论提出了更适合用于进行水轮机空化状态预警的特征指标,并结合BiGRU深度学习模型实现水轮机初生空化在线预警,适用范围广,可以有效解决强环境噪声干扰下空化识别不准确的问题。本发明中基于水轮机空化原理将空化信号瞬时能量、随机性和复杂性与空化状态之间的规律作为判断水轮机空化状态的依据,可以有效提高空化状态识别准确率。当空化未发生时,水轮机声学信号熵率与瞬时能量值随空化系数降低而缓慢增加,当水轮机内发生空化时,声学信号随机性与复杂性迅速增加,熵率与瞬时能量值迅速增大。应该是随着水轮机空化系数的进一步降低,空泡群会淹没声学信号信息,并由于空泡之间的共振作用吸收声波能量,使得熵率与瞬时能量值有所减小;当空化程度进一步加深,这种现象消失,熵率与瞬时能量值又开始迅速增加。
附图说明
图1是本发明基于声学信号的水轮机初生空化预警方法流程图;
图2是本发明方法中低频声学信号熵率随空化系数变化趋势图;
图3是本发明方法中高频声学信号能量值随空化系数变化趋势图。
具体实施方式
下面结合附图及具体实施方式对本发明进行详细说明。
实施例1
本发明基于声学信号的水轮机初生空化预警方法,基于水轮机初生空化时声学信号随机性、复杂性以及能量的变化规律,通过采集水轮机组声学信号来实现水轮机初生空化的在线预警。如图1所示,具体按照以下步骤实施:在水轮机上安装声学传感器,采集各工况运行参数下不同空化状态水轮机运行时的声学信号;对声学信号进行去噪处理,去除背景噪声等因素的干扰,得到去噪声学信号;随后对去噪声学信号进行精细复合多尺度波动色散熵特征提取,更好地刻画水轮机系统的非线性行为以及随机性、复杂性特征;计算高频声学信号的Teager能量算子序列,对能量序列进行快速傅里叶变换,将其主频幅值作为表征高频声学信号的能量特征,以此计算低频空化声学信号熵率和高频空化声学信号的信号能量值。
如图2和图3所示,当水轮机未发生空化时,运行状态较为稳定,系统的复杂性与随机性较为稳定,所采集的低频声学信号的熵率较小,且几乎保持不变,高频声学信号的瞬时能量值同样较低。当水轮机内部压力减小,转轮边缘产生少量空泡,空泡的溃灭等随机性行为使得所采集的声学信号的随机性及复杂性开始迅速上升,空泡溃灭所产生的高频声波能量增大,使得低频声学信号的熵率值与高频声学信号的瞬时能量值开始迅速增大;当空化程度进一步加深,空泡数量增大,水轮机内部转变为云空化时,水体中的空泡改变了原有的流态物性结构,空泡对于声波的散射与空泡群中空泡共振对声波的吸收作用使得所采集的声波信号能量减小,声波传递至水听器时许多信息被淹没,当空化程度继续加深,这种现象会消失,预警指标熵率和能量的数值开始迅速上升。
随后构建水轮机工况运行参数与初生空化声学信号预警特征值之间关系的预测模型,将水轮机工况运行参数作为输入数据,初生空化声学信号熵率和高频声学信号的信号能量值作为输出数据,每组输入数据及对应的输出数据作为训练样本训练预测模型;采集工况运行参数已知、空化状态未知状态的实况声学信号,计算实况声学信号的熵率和能量值,将实况声学信号对应的工况运行参数输入训练后的预测模型中获得熵率预警阈值与能量值的预警阈值;当实况声学信号的熵率超过熵率预警阈值或实况声学信号的能量值超过能量值的预警阈值,发出初生空化预警,通过本发明方法监测水轮机组是否发生初生空化,及时发现并判断机组状态,以配合合理的维护检修措施。
实施例2
基于声学信号的水轮机初生空化预警方法,具体按照以下步骤实施:
由于不同体积空泡的特征频率不同,因此为采集信息足够丰富的空化信号,在尾水管直锥段接近转轮出口处安装水听器与声发射传感器,其中水听器采样频率设置为60kHz,声发射传感器采样频率设置为2MHz;固定水轮机的运行参数,并在不同的空化状态下利用水听器与声发射传感器同步采集水轮机空化声学信号,通过数据采集卡发送到计算机并记录声学信号对应时间,直至当转轮叶片上出现空泡时,将此状态判断为初生空化状态,通过水听器与声发射传感器采集声学信号,当一个工况的空化声学信号采集完成后,分别改变工况运行参数,获得不同工况下声学信号。工况运行参数包括导叶开度α0,转速n,水头H,下游水位Hd、水轮机出力P和扭矩M。
以白噪声幅值和加噪次数为粒子坐标,设定好坐标参数搜索范围,初始化粒子群算法参数;
在当前白噪声幅值和加噪次数下,对所采集的声学信号序列进行自适应噪声完全集合经验模态分解,计算各IMF分量的分形维数过程为:
设定IMF分量序列的时间序列为x(t),t=1,2,NN表示整个时间序列的长度;构造新的以延迟时间k为参数的时间序列/>表示为:
将时间序列展开后计算时间序列矩阵的曲线长度Lm(k),表示为:
式中,N表示整个时间序列的长度,是归一化因子,[/>]表示对/>取整,上式由k=1到/>循环计算,求得一个k值下所有L m(k)的均值L(k),在双对数坐标系下利用最小二乘法拟合/>与/>得到直线斜率即是IMF分量的分形维数。
采用分形维数作为算法参数优化中的适应度函数,计算速度更快,准确度更高,通过提取可以精确刻画信号非线性动力学特性的精细复合多尺度波动色散熵和代表瞬时能量的新型能量指标,并结合深度学习算法实现初生空化在线预警,适用范围更广,识别性能较传统方法而言精度更高,可以有效提高电站的智能化运维水平。
将分形维数的和作为适应度函数;将初始粒子适应度作为个体最优值,粒子历史最优适应度值作为全局最优值。
对粒子位置以及速度、粒子种群最优解进行更新,更新公式为:
式中,是个体历史最优位置;/>是全局历史最优位置;/>是目前粒子位置,/>是下一次迭代的粒子速度即更新后的速度;w代表惯性权重,/>,/>分别取0.9和0.4,c 1c 2为学习因子,/>,/>是(0,1)之间的随机数,下标ij表示第i个粒子和第j维,t代表迭代次数,/>表示最大进化代数,计算更新后的粒子适应度值,将当前粒子种群适应度值与个体最优值、全局最优值进行比较,更新粒子种群的个体最优以及全局最优极值;
判断是否满足迭代条件,设置容忍度阈值、最大容忍代数,最大迭代次数,计算当前最优适应度与上一代最优适应度的相对变化量作为容忍度,比较容忍度与容忍度阈值大小,若前者小则容忍代数加1,判断容忍代数是否大于最大容忍代数或者迭代次数是否超过最大迭代次数,若满足,则输出目前全局最优解对应的白噪声幅值和加噪次数,作为最优白噪声幅值£和加噪次数Ne,若不满足,则返回对粒子位置以及速度、粒子种群最优解进行更新。
根据最优白噪声幅值£和加噪次数Ne对采集的声学信号进行分解,获得K个IMF分量。
计算IMF分量与分解前的声学信号之间的最大互信息系数:
其中,XY分别为IMF分量和分解前的声学信号,表示XY之间的互信息,表示变量间的联合概率密度,/>为变量间的最大互信息系数,B表示网格划分系数;
选择最大互信息系数前5阶IMF分量重构声学信号,得到去噪声学信号。
计算低频空化声学信号熵率和高频空化声学信号的信号能量值;
构建水轮机工况运行参数与初生空化声学信号预警特征值之间关系的预测模型,将水轮机工况运行参数作为输入数据,初生空化声学信号熵率和高频声学信号的信号能量值作为输出数据,每组输入数据及对应的输出数据作为训练样本训练预测模型;
采集工况运行参数已知、空化状态未知状态的实况声学信号,计算实况声学信号的熵率和能量值,将实况声学信号对应的工况运行参数输入训练后的预测模型中获得熵率预警阈值与能量值的预警阈值;当实况声学信号的熵率超过熵率预警阈值或实况声学信号的能量值超过能量值的预警阈值,发出初生空化预警。
本实施例中,通过改进自适应噪声完全集合经验模态分解对所采集的声学信号进行去噪处理,并提出了基于分形维数的粒子群算法对改进自适应噪声完全集合经验模态分解的参数进行自适应优化,避免了人工设置模型参数的弊端,减小了模态混叠的可能性,减小了信号中所包含的环境噪声对后续的熵率及瞬时能量值计算的影响,提高了水轮机初生空化预警的准确率。
实施例3
基于声学信号的水轮机初生空化预警方法,具体按照以下步骤实施:
在水轮机上安装声学传感器,采集各工况运行参数下不同空化状态水轮机运行时的声学信号;工况运行参数包括导叶开度α0,转速n,水头H,下游水位Hd、水轮机出力P和扭矩M。
对声学信号进行去噪处理,得到去噪声学信号;
将去噪声学信号进行粗粒化,粗粒化过程中不同的起点创建不同的时间序列,重构后声学序列的粗粒化序列为:
其中,L为重构声学序列长度;τ为尺度因子,u b表示重构声学序列的第b个值;
计算粗粒化时间序列中色散模式π的平均概率
其中,色散模式π为重构空间中序列的组合形式,其中为不同起始点每个可能的色散模式/>的概率;
其中,不同起始点每个可能的色散模式的概率计算过程为:
对于每个粗粒化时间序列,通过下式将变量映射到{1,2,/>,c}各类别中:
式中,为序列的期望,/>为序列的标准差,round为取整函数,c为类别个数;/>表示/>通过正态累积分布函数变换后的累积概率值;
将嵌入向量映射到/>种基于波动的色散模式/>中,则有:
式中m为嵌入维数,d为时延长度。
每个可能的色散模式的概率由下式计算:
式中,count是指映射至/>的模式个数,再除以嵌入维数为m对应的嵌入向量总数即为色散模式的概率/>N表示粗粒化序列长度。
计算色散模式概率平均值的色散熵:
利用色散熵的前八个尺度的熵值进行最小二乘法拟合,得到的曲线斜率即为低频声学信号的熵率;
由下式计算重构后的高频声学信号的Teager能量算子序列,重构公式为:
式中,表示Teager能量算子序列的第n个值,/>表示重构高频信号序列中的第n个数,/>表示重构高频信号序列中的第n+1个数,/>表示重构高频信号序列中的第n-1个数;
对信号Teager能量算子序列进行快速傅里叶变换,其主频幅值即为高频声学信号的信号能量值。
构建水轮机工况运行参数与初生空化声学信号预警特征值之间关系的预测模型,即为BiGRU深度学习模型,将水轮机每个工况运行参数作为输入数据,每个工况运行参数对应的初生空化声学信号熵率和高频声学信号的信号能量值作为输出数据,每个输入数据及对应的输出数据作为训练样本,对训练样本进行归一化处理后输入BiGRU深度学习模型,采用十折交叉验证与网格搜索优化预测模型的超参数、BiGRU层数、神经元个数以及学习率,当训练次数达到设定最大次数或模型误差小于1%时保存BiGRU深度学习模型,作为训练后的预测模型。
采集工况运行参数已知、空化状态未知状态的实况声学信号,计算实况声学信号的熵率和能量值,将实况声学信号对应的工况运行参数输入训练后的预测模型中获得熵率预警阈值与能量值的预警阈值;当实况声学信号的熵率超过熵率预警阈值或实况声学信号的能量值超过能量值的预警阈值,发出初生空化预警。
本实施例中,基于水轮机空化原理,通过分析水轮机空化状态与声学信号之间的关联特性,基于信号复杂度与瞬时能量理论提出了更有效的空化特征指标,熵率与瞬时能量值;有效避免了传统时频域特征提取中易受噪声影响的弊端。通过Teager能量算子方法放大了空化声学信号中的微弱冲击,提高水轮机空化预警预测模型的准确度,计算方法简单且提高了空化预警的准确性。
本发明预警方法,基于不同体积的空泡释放的声学信号特征频率不同,因此通过监测水听器和声发射传感器采集空化所造成的低频声学信号和高频声学信号,利用声学信号瞬时能量、随机性、复杂性随空化系数的变化规律来判断水轮机是否发生空化。水轮机未发生空化时,低频声学信号的结构复杂性较弱,因此熵率较小;高频信号能量稳定在一个较小的值,几乎不变。水轮机由无空化到初生空化的过程中,内部流态逐渐失稳,熵率逐渐增大,在初生空化附近取得局部极大值,高频信号能量逐渐增大。当水轮机内部空化程度继续加深至临界空化前,高频信号能量继续增大,而低频信号熵率有一个短暂降低的过程,与此时水轮机特性变化趋势结论一致,当空化程度继续增大时,熵率迅速增大。由于水中空泡的数量较多,改变了水体的物性,对声波能量的传递起到缓解作用,使得声发射信号能量值降低。因此本文通过选取初生空化工况点附近的低频声学信号熵率以及高频信号的瞬时能量指标作为预警指标,当信号的预警指标值超过预警阈值时便代表此时水轮机中已出现了空化现象。

Claims (8)

1.基于声学信号的水轮机初生空化预警方法,其特征在于,具体按照以下步骤实施:
步骤1、通过水轮机上安装的声学传感器获取空化声学信号数据集,包含各工况运行参数下无空化及空化状态下水轮机运行时的声学信号;
步骤2、对声学信号进行去噪处理,得到去噪声学信号;具体过程为:
步骤2.1、输入所采集的声学信号对应的时间序列,确定改进自适应噪声完全集合经验模态分解的最优白噪声幅值£和加噪次数Ne
步骤2.2、根据最优白噪声幅值£和加噪次数Ne对采集的声学信号进行分解,获得K个IMF分量;
步骤2.3、分别计算K个IMF分量与分解前的声学信号之间的最大互信息系数,根据最大互信息系数对IMF分量筛选、重构,得到去噪声学信号;
步骤3、计算低频空化声学信号熵率和高频空化声学信号的信号能量值;具体过程为:
步骤3.1、将去噪声学信号进行粗粒化,粗粒化过程中不同的起点创建不同的时间序列,重构后声学序列的粗粒化序列为:
其中,L为重构声学序列长度;τ为尺度因子,u b表示重构声学序列的第b个值;
步骤3.2、计算粗粒化时间序列中色散模式π的平均概率
其中,色散模式π为重构空间中序列的组合形式,其中为不同起始点每个可能的色散模式/>的概率;
步骤3.3、计算色散模式概率平均值的色散熵:
其中,c为类别个数,m为嵌入维数,d为时延长度;
步骤3.4、利用色散熵的前八个尺度的熵值进行最小二乘法拟合,得到的曲线斜率即为低频声学信号的熵率;
步骤3.5、由下式计算重构后的高频声学信号的Teager能量算子序列,重构公式为:
式中,表示Teager能量算子序列的第n个值,/>表示重构高频信号序列中的第n个数,/>表示重构高频信号序列中的第n+1个数,/>表示重构高频信号序列中的第n-1个数;
对信号Teager能量算子序列进行快速傅里叶变换,其主频幅值即为高频声学信号的信号能量值;
步骤4、构建水轮机工况运行参数与初生空化声学信号预警特征值之间关系的预测模型,将水轮机工况运行参数作为输入数据,初生空化声学信号熵率和高频声学信号的信号能量值作为输出数据,每组输入数据及对应的输出数据作为训练样本训练预测模型;
步骤5、采集工况运行参数已知、空化状态未知状态的实况声学信号,计算实况声学信号的熵率和能量值,将实况声学信号对应的工况运行参数输入训练后的预测模型中获得熵率预警阈值与能量值的预警阈值;当实况声学信号的熵率超过熵率预警阈值或实况声学信号的能量值超过能量值的预警阈值,发出初生空化预警。
2.根据权利要求1所述基于声学信号的水轮机初生空化预警方法,其特征在于,所述工况运行参数包括导叶开度α0,转速n,水头H,下游水位Hd、水轮机出力P和扭矩M。
3.根据权利要求1所述基于声学信号的水轮机初生空化预警方法,其特征在于,步骤1具体过程为:在尾水管直锥段接近转轮出口上安装水听器与声发射传感器,其中水听器采样频率设置为60kHz,声发射传感器采样频率设置为2MHz;固定水轮机的运行参数,并在不同的空化状态下利用水听器与声发射传感器同步采集水轮机空化声学信号,通过数据采集卡发送到计算机并记录声学信号对应时间,直至当转轮叶片上出现空泡时,将此状态判断为初生空化状态,通过水听器与声发射传感器采集声学信号,当一个工况的空化声学信号采集完成后,分别改变工况运行参数,获得不同工况下声学信号。
4.根据权利要求1所述基于声学信号的水轮机初生空化预警方法,其特征在于,步骤2.1具体过程为:
步骤2.1.1、以白噪声幅值和加噪次数为粒子坐标,设定好坐标参数搜索范围,初始化粒子群算法参数;
步骤2.1.2、在当前白噪声幅值和加噪次数下,对所采集的声学信号序列进行自适应噪声完全集合经验模态分解,计算各IMF分量的分形维数,将分形维数的和作为适应度函数;将初始粒子适应度作为个体最优值,粒子历史最优适应度值作为全局最优值;
步骤2.1.3、对粒子位置以及速度、粒子种群最优解进行更新,更新公式为:
式中,是个体历史最优位置;/>是全局历史最优位置;/>是目前粒子位置,/>表示下一次迭代的粒子位置即更新后的粒子位置,/>是下一次迭代的粒子速度即更新后的速度;w代表惯性权重,/>,/>分别取0.9和0.4,c 1c 2为学习因子,/>,/>是(0,1)之间的随机数,下标ij表示第i个粒子和第j维,t代表迭代次数,/>表示最大进化代数,计算更新后的粒子适应度值,将当前粒子种群适应度值与个体最优值、全局最优值进行比较,更新粒子种群的个体最优以及全局最优极值;
步骤2.1.4、判断是否满足迭代条件,设置容忍度阈值、最大容忍代数,最大迭代次数,计算当前最优适应度与上一代最优适应度的相对变化量作为容忍度,比较容忍度与容忍度阈值大小,若前者小则容忍代数加1,判断容忍代数是否大于最大容忍代数或者迭代次数是否超过最大迭代次数,若满足,则输出目前全局最优解对应的白噪声幅值和加噪次数,作为最优白噪声幅值£和加噪次数Ne,若不满足,则返回步骤2.1.3。
5.根据权利要求4所述基于声学信号的水轮机初生空化预警方法,其特征在于,所述计算各IMF分量的分形维数过程为:
设定IMF分量序列的时间序列为x(t),t=1,2,NN表示整个时间序列的长度;构造新的以延迟时间k为参数的时间序列/>表示为:
将时间序列展开后计算时间序列矩阵的曲线长度L m(k),表示为:
式中,N表示整个时间序列的长度,是归一化因子,[/>]表示对/>取整,上式由k=1到/>循环计算,求得所有L m(k)的均值L(k),在双对数坐标系下利用最小二乘法拟合/>与/>得到直线斜率即是IMF分量的分形维数。
6.根据权利要求1所述基于声学信号的水轮机初生空化预警方法,其特征在于,步骤2.3具体过程为:
计算IMF分量与分解前的声学信号之间的最大互信息系数:
其中,XY分别为IMF分量和分解前的声学信号,表示XY之间的互信息,表示变量间的联合概率密度,p(x)表示x的边缘概率分布,p(y)表示y的边缘概率分布,/>为变量间的最大互信息系数,B表示网格划分系数;
选择最大互信息系数前5阶IMF分量重构声学信号,得到去噪声学信号。
7.根据权利要求1所述基于声学信号的水轮机初生空化预警方法,其特征在于,步骤3.2中所述不同起始点每个可能的色散模式的概率计算过程为:
对于每个粗粒化时间序列,通过下式将变量映射到{1,2,/>, c}各类别中:
式中,为序列的期望,/>为序列的标准差,round为取整函数,c为类别个数;/>表示/>通过正态累积分布函数变换后的累积概率值,/>表示y j经过线性变换后映射至第c类别中的集合;
将嵌入向量映射/>种基于波动的色散模式/>中,则有:
式中m为嵌入维数,d为时延长度;
每个可能的色散模式的概率由下式计算:
式中,count是指映射至/>的模式个数,再除以嵌入维数为m对应的嵌入向量总数即为色散模式的概率/>N表示粗粒化序列长度。
8.根据权利要求1所述基于声学信号的水轮机初生空化预警方法,其特征在于,步骤4中所述水轮机工况运行参数与初生空化声学信号预警特征值之间关系的预测模型为BiGRU深度学习模型,将水轮机每个工况运行参数作为输入数据,每个工况运行参数对应的初生空化声学信号熵率和高频声学信号的信号能量值作为输出数据,每个输入数据及对应的输出数据作为训练样本,对训练样本进行归一化处理后输入BiGRU深度学习模型,采用十折交叉验证与网格搜索优化预测模型的超参数、BiGRU层数、神经元个数以及学习率,当训练次数达到设定最大次数或模型误差小于1%时保存BiGRU深度学习模型,作为训练后的预测模型。
CN202311787284.9A 2023-12-25 2023-12-25 基于声学信号的水轮机初生空化预警方法 Active CN117476039B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202311787284.9A CN117476039B (zh) 2023-12-25 2023-12-25 基于声学信号的水轮机初生空化预警方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202311787284.9A CN117476039B (zh) 2023-12-25 2023-12-25 基于声学信号的水轮机初生空化预警方法

Publications (2)

Publication Number Publication Date
CN117476039A CN117476039A (zh) 2024-01-30
CN117476039B true CN117476039B (zh) 2024-03-08

Family

ID=89623885

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202311787284.9A Active CN117476039B (zh) 2023-12-25 2023-12-25 基于声学信号的水轮机初生空化预警方法

Country Status (1)

Country Link
CN (1) CN117476039B (zh)

Citations (17)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101017515A (zh) * 2007-02-12 2007-08-15 三峡大学 水电站水力学的水流精细模拟方法
CN101813512A (zh) * 2009-12-07 2010-08-25 哈尔滨电机厂有限责任公司 采用计算机程序确定模型水轮机转轮叶片初生空化的声学方法
CN110738115A (zh) * 2019-09-12 2020-01-31 浙江大学 一种基于脉冲频率特征模式识别的螺旋桨空化程度识别方法
CN110987494A (zh) * 2019-12-02 2020-04-10 吉林松江河水力发电有限责任公司 一种基于声发射对水轮机空化状态监测的方法
CN111259864A (zh) * 2020-03-04 2020-06-09 哈尔滨理工大学 一种水轮机运转状态识别方法
KR20210010404A (ko) * 2019-07-17 2021-01-27 충남대학교산학협력단 음원 출력 제어 장치 및 그 방법
CN112729836A (zh) * 2020-11-30 2021-04-30 华电电力科学研究院有限公司 一种循环改进型的水轮机空蚀初生状态判别系统及其方法
CN113255848A (zh) * 2021-07-08 2021-08-13 浙江大学 基于大数据学习的水轮机空化声信号辨识方法
CN114462448A (zh) * 2022-01-17 2022-05-10 哈尔滨理工大学 一种水轮机运转状态识别方法
CN115017445A (zh) * 2022-05-31 2022-09-06 安徽工业大学 一种基于ceemdan结合近似熵的螺旋桨状态识别方法
CN115293004A (zh) * 2022-08-19 2022-11-04 浙江理工大学 一种基于多尺度空化模型的空蚀预测方法
CN115600126A (zh) * 2022-10-10 2023-01-13 江苏大学镇江流体工程装备技术研究院(Cn) 一种基于ceemd-drsn的离心泵空化状态识别方法
CN115730241A (zh) * 2022-11-17 2023-03-03 国网新源控股有限公司 一种水轮机空化噪声识别模型的构建方法
CN115774841A (zh) * 2022-11-17 2023-03-10 国网新源控股有限公司 一种基于组合谱特征提取的水轮机空化现象智能识别方法
CN116049637A (zh) * 2023-02-28 2023-05-02 西安理工大学 基于流动噪声信号的气液两相流流型识别方法
CN116432527A (zh) * 2023-04-07 2023-07-14 中国长江三峡集团有限公司 一种空化预测方法、装置、存储介质及电子设备
CN117037841A (zh) * 2023-07-31 2023-11-10 西安电子科技大学 基于层级过渡网络的声学信号层级空化强度识别方法

Patent Citations (17)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101017515A (zh) * 2007-02-12 2007-08-15 三峡大学 水电站水力学的水流精细模拟方法
CN101813512A (zh) * 2009-12-07 2010-08-25 哈尔滨电机厂有限责任公司 采用计算机程序确定模型水轮机转轮叶片初生空化的声学方法
KR20210010404A (ko) * 2019-07-17 2021-01-27 충남대학교산학협력단 음원 출력 제어 장치 및 그 방법
CN110738115A (zh) * 2019-09-12 2020-01-31 浙江大学 一种基于脉冲频率特征模式识别的螺旋桨空化程度识别方法
CN110987494A (zh) * 2019-12-02 2020-04-10 吉林松江河水力发电有限责任公司 一种基于声发射对水轮机空化状态监测的方法
CN111259864A (zh) * 2020-03-04 2020-06-09 哈尔滨理工大学 一种水轮机运转状态识别方法
CN112729836A (zh) * 2020-11-30 2021-04-30 华电电力科学研究院有限公司 一种循环改进型的水轮机空蚀初生状态判别系统及其方法
CN113255848A (zh) * 2021-07-08 2021-08-13 浙江大学 基于大数据学习的水轮机空化声信号辨识方法
CN114462448A (zh) * 2022-01-17 2022-05-10 哈尔滨理工大学 一种水轮机运转状态识别方法
CN115017445A (zh) * 2022-05-31 2022-09-06 安徽工业大学 一种基于ceemdan结合近似熵的螺旋桨状态识别方法
CN115293004A (zh) * 2022-08-19 2022-11-04 浙江理工大学 一种基于多尺度空化模型的空蚀预测方法
CN115600126A (zh) * 2022-10-10 2023-01-13 江苏大学镇江流体工程装备技术研究院(Cn) 一种基于ceemd-drsn的离心泵空化状态识别方法
CN115730241A (zh) * 2022-11-17 2023-03-03 国网新源控股有限公司 一种水轮机空化噪声识别模型的构建方法
CN115774841A (zh) * 2022-11-17 2023-03-10 国网新源控股有限公司 一种基于组合谱特征提取的水轮机空化现象智能识别方法
CN116049637A (zh) * 2023-02-28 2023-05-02 西安理工大学 基于流动噪声信号的气液两相流流型识别方法
CN116432527A (zh) * 2023-04-07 2023-07-14 中国长江三峡集团有限公司 一种空化预测方法、装置、存储介质及电子设备
CN117037841A (zh) * 2023-07-31 2023-11-10 西安电子科技大学 基于层级过渡网络的声学信号层级空化强度识别方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
A Novel Acoustic Method for Cavitation Identification of Propeller;Yang L等;Journal of Marine Science and Engineering;20220901;第10卷(第9期);1-21 *
Fuzzy Dispersion Entropy: A Nonlinear Measure for Signal Analysis;M. Rostaghi等;IEEE Transactions on Fuzzy Systems;20211118;第30卷(第9期);3785-3796 *
水动力噪声计算方法综述;李环等;中国舰船研究;20160317;第11卷(第02期);72-89 *

Also Published As

Publication number Publication date
CN117476039A (zh) 2024-01-30

Similar Documents

Publication Publication Date Title
Wang et al. A data indicator-based deep belief networks to detect multiple faults in axial piston pumps
Lv et al. Vibration signal-based early fault prognosis: Status quo and applications
He et al. A joint adaptive wavelet filter and morphological signal processing method for weak mechanical impulse extraction
CN115688018B (zh) 一种多工况下轴承的状态监测与故障诊断方法
CN113420691A (zh) 一种基于皮尔逊相关系数的混合域特征轴承故障诊断方法
Chen et al. Fault feature extraction and diagnosis of gearbox based on EEMD and deep briefs network
CN109827777A (zh) 基于偏最小二乘法极限学习机的滚动轴承故障预测方法
Xue et al. Intelligent diagnosis method for centrifugal pump system using vibration signal and support vector machine
JP2021018818A (ja) ウェーブレットおよび主成分解析に基づくプロペラのキャビテーション状態の検出方法
CN117116291A (zh) 一种含沙水流冲击水轮机的声信号处理方法
CN114263621B (zh) 一种离心泵空化故障诊断模拟的试验方法及系统
Ding et al. Gear fault diagnosis based on VMD sample entropy and discrete hopfield neural network
CN117476039B (zh) 基于声学信号的水轮机初生空化预警方法
Yao et al. Multiband weights-induced periodic sparse representation for bearing incipient fault diagnosis
Zhang et al. Nonstationary significant wave height forecasting with a hybrid VMD-CNN model
Ren et al. Research on fault feature extraction of hydropower units based on adaptive stochastic resonance and fourier decomposition method
Zhang et al. A novel hybrid compound fault pattern identification method for gearbox based on NIC, MFDFA and WOASVM
Zhou et al. An adaptive VMD method based on improved GOA to extract early fault feature of rolling bearings
Wan et al. Adaptive asymmetric real Laplace wavelet filtering and its application on rolling bearing early fault diagnosis
CN114528869A (zh) 一种风机叶片损伤识别方法
CN114462448A (zh) 一种水轮机运转状态识别方法
CN112926504A (zh) 一种基于降噪自编码器的声发射信号去噪方法
Meng et al. Rolling bearing fault diagnosis method based on MCMF and SAIMFE
Wang et al. Research on Rolling Bearing Fault Diagnosis Based on Volterra Kernel Identification and KPCA
Xu et al. A Novel Empirical Variational Mode Decomposition for Early Fault Feature Extraction

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant