CN115204362A - 一种机械设备剩余寿命区间预测方法 - Google Patents

一种机械设备剩余寿命区间预测方法 Download PDF

Info

Publication number
CN115204362A
CN115204362A CN202210736162.6A CN202210736162A CN115204362A CN 115204362 A CN115204362 A CN 115204362A CN 202210736162 A CN202210736162 A CN 202210736162A CN 115204362 A CN115204362 A CN 115204362A
Authority
CN
China
Prior art keywords
layer
time
feature
prediction
attention
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.)
Pending
Application number
CN202210736162.6A
Other languages
English (en)
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.)
Naval University of Engineering PLA
Original Assignee
Naval University of Engineering PLA
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 Naval University of Engineering PLA filed Critical Naval University of Engineering PLA
Priority to CN202210736162.6A priority Critical patent/CN115204362A/zh
Publication of CN115204362A publication Critical patent/CN115204362A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01D—MEASURING NOT SPECIALLY ADAPTED FOR A SPECIFIC VARIABLE; ARRANGEMENTS FOR MEASURING TWO OR MORE VARIABLES NOT COVERED IN A SINGLE OTHER SUBCLASS; TARIFF METERING APPARATUS; MEASURING OR TESTING NOT OTHERWISE PROVIDED FOR
    • G01D21/00—Measuring or testing not otherwise provided for
    • G01D21/02—Measuring two or more variables by means not covered by a single other subclass
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/04—Architecture, e.g. interconnection topology
    • G06N3/049—Temporal neural networks, e.g. delay elements, oscillating neurons or pulsed inputs
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/08—Learning methods
    • Y—GENERAL 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
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02P—CLIMATE CHANGE MITIGATION TECHNOLOGIES IN THE PRODUCTION OR PROCESSING OF GOODS
    • Y02P90/00—Enabling technologies with a potential contribution to greenhouse gas [GHG] emissions mitigation
    • Y02P90/30—Computing systems specially adapted for manufacturing

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Molecular Biology (AREA)
  • Biophysics (AREA)
  • Computational Linguistics (AREA)
  • Artificial Intelligence (AREA)
  • Evolutionary Computation (AREA)
  • General Health & Medical Sciences (AREA)
  • Biomedical Technology (AREA)
  • Computing Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Mathematical Physics (AREA)
  • Software Systems (AREA)
  • Health & Medical Sciences (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种机械设备剩余寿命区间预测方法，针对统计数据驱动的剩余寿命预测模型难以描述多维传感器数据之间的耦合关系以及数据驱动融合模型存在关键信息丢失的问题，首先，利用特征注意力机制(Featureattention,FA)从多维度、非线性和大规模的传感器信号中提取出关键特征向量；然后，采用BiGRU网络从前向和后向两个方向对注意力加权特征的时变特性进行建模学习，并通过最大似然估计损失函数来训练网络参数，获得隐含状态输出向量的概率分布，从而计算出基于LN分布的PDF，实现部件不确定性的衡量的剩余寿命区间预测方法。仿真结果表明本发明提出的基于FA‑LN‑BiGRU融合模型的剩余寿命区间预测方法能够可靠、准确地得到RUL预测区间。

Description

一种机械设备剩余寿命区间预测方法
技术领域
本发明涉及机械设备的寿命预测领域，特别涉及一种机械设备剩余寿命区间预测方法。
背景技术
机械设备(如飞机上的发动机)的剩余使用寿命(RemainingUsefulLife，简称RUL)预测通过分析机械设备的历史性能退化过程，动态感知设备健康状态的未来变化，从而预测设备衰退到失效阈值的剩余时间，目前主要包括基于机理模型、基于数据驱动和基于两者混合驱动3类预测方法。基于机理模型和混合驱动的预测方法能够精确的描述系统退化过程，通过已知的力学原理和机械原理建立明确的数学模型，进而分析设备的性能退化过程。然而以涡轮发动机为代表的大型设备机械结构复杂、任务环境恶劣和故障模式多样，建立精确的物理失效模型代价过高，难以部署在实际应用过程中。数据驱动的预测方法可以充分挖掘设备运行过程中的历史监测数据，通过提取设备性能劣化特征来构建RUL预测模型，该方法无需先验知识且规避了机械设备复杂机理的建模过程，逐渐成为当下RUL预测方法研究的热点。
数据驱动的预测模型主要分为机器学习模型和统计分析模型两大类。深度学习方法可以充分学习时间序列蕴含的信息，从而对RUL进行更为精确的预测，无需确切的物理模型和先验知识，但该方法只能提供RUL的点估计而非区间估计，难以量化RUL预测过程中不确定性的概率分布，增加后续健康管理中制定维修决策的风险。统计分析方法通过统计模型或随机过程对设备历史失效数据进行分布拟合，利用参数估计构造RUL的概率密度分布函数(probability density function,PDF)，然后对实时传感器数据进行统计分析来实现RUL的预测。统计数据驱动的方法所建立的模型参数和随机分布与设备退化失效过程紧密联系，具有较强的可解释性，同时可以准确得到RUL预测的概率分布，在度量预测的不确定性上具有天然优势。然而随着传感器和信息技术的迅猛发展，通常采用多种组网传感器收集多源数据来全面描述设备运行状态，监测大数据往往呈现出大规模、碎片化、非线性、稀疏化和相互耦合的特点。基于信息融合将多性能退化变量转换为融合变量的方式虽然可以应用于传统单变量退化建模，但是难以描述多维数据的耦合关系，存在监测数据退化信息利用不充分的问题。
鉴于此，数据驱动下属的机器学习和统计分析融合RUL预测方法得到了国内外学者的广泛关注。尽管上述方法可以综合机器学习方法和统计分析方法两者的优势，可是依然存在以下局限：①深度网络提取的退化特征难以匹配之后所采用的随机模型：深度网络的退化特征提取过程和随机模型建模过程是相互独立的，退化特征能否适用于随机模型有待商榷，且上述方法只能通过指标筛选技术从多维退化特征中选择最优指标用于随机退化建模，在非线性转换为线性的过程中存在关键信息丢失的风险；②实际部署存在困难：上述方法在实际应用中需要大量的计算资源，整个训练和预测过程被分为两步，且需要多次迭代才能得到最优置信区间，不适合实时RUL预测。
发明内容
针对上述存在的问题，本发明提出了一种机械设备剩余寿命(remain usefullife,RUL)区间预测方法，其是一种输入端直接到输出端的RUL预测方法，无需人工提取特征，使得模型对RUL变化趋势的拟合更加精准，同时基于特征注意力机制的对数正态分布模型和双向门控循环单元(feature attention-lognorm-bidirectional gated recurrentunit,FA-LN-BiGRU)实现机械设备剩余寿命区间预测。
为了实现上述目的，本发明所采用的技术方案如下：
一种机械设备剩余寿命区间预测方法，其特征在于，包括以下步骤：
步骤1：获取机械设备的传感器监测数据；
步骤2：构建FA-LN-BiGRU预测模型；
步骤3：将得到的传感器数据输入至FA-LN-BiGRU预测模型中，通过所述预测模型计算并输出航空发动机的RUL的PDF参数预测值；
步骤4：根据PDF参数预测值对发送机的剩余寿命区间进行预测。
进一步地，所述FA-LN-BiGRU预测模型包括依次连接的输入层、特征加权层、特征学习层和输出层；
输入层，用于将得到的航空发动机的历史传感器监测数据输入到特征加权层；
特征加权层，用于通过轻量级特征注意力机制提取输入的传感器监测数据的注意力加权特征向量；
特征学习层，用于通过双向BiGRU网络提取注意力加权特征向量的时序依赖关系；
输出层，用于计算航空发动机的RUL的PDF参数预测值。
进一步地，所述特征加权层包括全局平均池层、两层全连接层和sigmoid激活函数；
全局平均池化层，用于将获取到的传感器监测数据进行平均池化操作，得到每个时间样本的均值信息；
全连接层，用于量化各个特征对应的注意力权重；
sigmoid激活函数，用于注意力权重做归一化处理。
进一步地，步骤3的具体操作步骤包括：
步骤31：输入层获取不同运行工况和类型的机械设备历史传感器数据，将所述传感器数据表示为多维特征时间序列矩阵：
Figure BDA0003715808390000041
其中，
Figure BDA0003715808390000042
为在t时刻的F个特征数值序列，
Figure BDA0003715808390000043
为第f个多维特征在T个历史时刻的数值序列；
步骤32：所述特征加权层利用特征注意力机制提取多维特征时间序列与目标建模参数之间的关联关系，通过循环迭代对所有时刻各个特征分别分配注意力权重，得到输入的多维特征时间序列的注意力加权特征；
步骤33：特征学习层将输入的注意力加权特征向量分时刻输入至双向BiGRU网络的前向层和反向层，得到t时刻前向输出和后向输出，再将前向层和反向层的输出合并得到总输出；
步骤34：输出层根据得到的特征学习层的总输出，实现特征学习层输出向量对最终参数输出的映射，得到当前时刻下的对数均值μt和对数方差σt，从而计算此时RUL的PDF；
步骤35：构建最终参数与目标输出的损失函数，通过损失函数更新网络参数，再通过更新后的网络循环迭代求解下一时刻的对数均值和对数方差，最终实现对数正态分布参数的最大似然估计。
进一步地，步骤32的具体操作步骤包括：
步骤321：将输入的多维特征时间序列x通过全局平均池进行池化操作，得到每个时间样本的均值信息，将维度为T×F的多维特征时间序列降维至1×F，该计算过程记为：
Figure BDA0003715808390000051
其中，yf表示通过全局平均池化后得到的向量；g(·)为特征维度的全局平均池化操作；Xf为第f个特征所有时间点的数据；
Figure BDA0003715808390000052
为第f个特征的第t个时间点的数据；
步骤322：通过两层全连接层量化各个特征对应的注意力权重，并通过sigmoid激活函数对权重做归一化处理，得到归一化后的注意力权重αf：
αf＝σ(W2δ(W1yf+b1)+b2) (2)
根据式(1)将式(2)改写为：
Figure BDA0003715808390000053
其中，e＝[e1,e2,…,eF]∈RF×1；
步骤323：将注意力权重与相应的多维特征时间序列x相乘，得到多维特征时间序列的注意力加权特征X′：
X′＝e⊙X＝[e1x1,e2x2,…,eFxF] (8)
其中，X′可记为
Figure BDA0003715808390000061
其中F为注意力加权特征的维度，T为时间长度。
进一步地，步骤33的具体操作步骤包括：
步骤331：将1时刻得到的注意力加权特征同时输入至前向层和后向层，直至T时刻；
步骤332：分别计算从0时刻到t时刻正向计算出向前隐含层的输出hf，以及从t时刻到0时刻反向计算出向后隐含层的输出hb；
步骤333：将hf和hb通过式(5)计算得到最终特征学习层的总输出H
h＝f(hf,hb) (5)
其中，f(·)是全连接层的映射函数。
进一步地，步骤34所述计算RUL的PDF的公式为：
y′＝[μt,σt]＝woH+bo (9)
其中，y′t为第t时刻模型的预测值；μt为当前时刻的对数均值，σt为当前时刻的对数方差；wo为权重矩阵；bo为偏置系数。
进一步地，步骤35的具体操作步骤包括：
步骤351：计算t0时刻对数均值和对数方差：
Figure BDA0003715808390000062
Figure BDA0003715808390000063
其中，wμ和wσ分别代表与μ和σ相关的整体网络输出权重；xt为t0时刻传感器监测数据；
Figure BDA0003715808390000064
为t0时刻的对数均值；
Figure BDA0003715808390000065
为t0时刻的对数方差；
步骤352：定义RUL预测的统计参数模型：给定一个随机变量Y＝ln(X)，服从均值为μ、方差为σ的正态分布，所述随机变量X满足参数为(μ,σ2)的LN分布；
步骤353：通过所述统计参数模型预测每个时刻分布的参数对数均值μ和对数方差σ；
步骤354：利用最大似然估计L(μ,σ|y)训练FA-LN-BiGRU网络的参数，且最大似然估计的计算公式为：
Figure BDA0003715808390000071
其中，y是RUL的实测值；
步骤355：定义FA-LN-BiGRU网络的损失函数，并利用反向传播算法更新网络参数：
Figure BDA0003715808390000072
本发明的有益效果是：
第一，本发明通过FA机制使得端到端的预测模型更多地专注于重要特征，无需人工干预就可以实现多维特征重要性的动态加权。
第二，LN-BiGRU可以深入挖掘时间序列前后的退化信息，同时基于LN分布重新构造损失函数，衡量设备退化过程中的不确定性，能够有效提高机械设备剩余寿命点预测和区间预测的准确度和可靠性；
第三，本发明与现有技术中的统计数据驱动模型与数据驱动融合模型相比，利用FA-LN-BiGRU融合模型对运行条件复杂和故障模式多变的多维机械设备监测数据进行处理，可生成兼顾可靠性和准确性的RUL预测区间，有利于后续运行规划、维修决策等健康管理措施。
附图说明
图1为特征注意力结构图；
图2为GRU网络结构图；
图3为BiGRU网络结构图；
图4a-图4c为LN分布的PDF曲线，其中的图4a为LN分布的统计特征值、图4b为不同μ对PDF的影响、图4c为不同σ对PDF的影响；
图5为FA-LN-BiGRU模型整体框架；
图6为FD001前12个特征分布图；
图7为RUL标签设置；
图8a-图8b为RUL预测结果，图8a为FD001测试集RUL预测结果、图8b为FD001发动机重新排序预测结果；
图9a-图9f为3台发动机RUL区间预测结果，其中图9a-图9f分别为34号引擎RUL区间预测结果、34号引擎最终监测点RUL的pdf、76号引擎RUL区间预测结果、76号引擎最终监测点RUL的pdf、100号引擎RUL区间预测结果、100号引擎最终监测点RUL的pdf。
具体实施方式
为了使本领域的普通技术人员能更好的理解本发明的技术方案，下面结合附图和实施例对本发明的技术方案做进一步的描述。
本发明提出一种端对端的RUL预测方法，即基于特征注意力机制的对数正态分布模型和双向门控循环单元(feature attention-lognorm-bidirectional gatedrecurrent unit,FA-LN-BiGRU)的机械设备剩余寿命区间预测方法。首先，利用特征注意力机制(Feature attention,FA)从多维度、非线性和大规模的传感器信号中提取出关键特征向量；其次，采用BiGRU网络从前向和后向两个方向对注意力加权特征的时变特性进行建模学习，并通过最大似然估计损失函数来训练网络参数，获得隐含状态输出向量的概率分布，从而计算出基于LN分布的PDF，实现部件不确定性的衡量。
一、特征注意力机制
对于工作在复杂环境下的机械设备，往往需要综合考虑多个特征的监测数据，才能更全面地评估其健康状态和性能退化趋势，而这些特征对RUL预测的贡献有多有少。为了有效建模特征与RUL预测的关联程度，从而提高预测模型特征提取能力，本发明提出一种轻量级特征注意力机制对原始输入数据的关键特征进行选择，具体结构如图1所示，其中T表示时间序列的时间长度，F表示特征长度。
首先，通过特征维度的全局平均池化(global average pooling,GAP)得到每个时间样本的均值信息，将T×F压缩至1×F，计算过程可表示为：
Figure BDA0003715808390000091
其中，yf表示通过全局平均池化后得到的向量；g(·)为特征维度的全局平均池化操作；Xf为第f个特征所有时间点的数据；
Figure BDA0003715808390000092
为第f个特征的第t个时间点的数据；
其次，通过两层全连接层(Linear)量化各个特征对应的注意力权重，并通过sigmoid激活函数对权重做归一化处理，计算过程如式(2)所示：
αf＝σ(W2δ(W1yf+b1)+b2) (2)
其中，W1、b1和W2、b2分别为两层Linear的权重和偏置；δ为ReLU激活函数，用于解决梯度消失的问题；σ为sigmoid激活函数，用来捕捉数据的非线性因素。
将归一化的注意力权重αf与相应的输入数据x相乘以实现第f个特征向量的增强或削弱，得到加权输入x′，如式(3)所示：
x′＝αf⊙x＝[α1x1,α2x2,…,αfxf] (3)。
二、双向门限循环单元
剩余寿命预测的主要途径是从机械设备传感器数据中学习与设备性能退化的时空信息，而LSTM网络可以充分利用数据的时间相关特性对时域序列数据进行处理，然而其内部结构复杂导致训练时间较长，不利于实际工业部署。
GRU对LSTM进行优化改进，通过将遗忘门和输入门组合在一起，减少了门的数量，在保证记忆能力的同时，提升了网络的训练效率。记忆单元是GRU网络的核心，由更新门、重置门和候选隐藏层构成，其网络拓扑结构如图2所示，其中ht-1、ht和ht+1分别为t-1、t和t+1时刻的隐藏状态，表示输入序列的记忆信息；xt-1、xt和xt+1为t-1、t和t+1时刻记忆单元的输入；箭头所指方向为信息流动方向。此外，×代表矩阵的数乘计算；+代表矩阵的相加计算；1-表示该路向前传播的数据为1-zt；σ为Sigmoid激活函数，将输出变换到[0,1]区间内；tanh为双曲正切激活函数，将输出变换到[-1,1]区间内。
每个GRU记忆单元的具体计算过程如下：
Figure BDA0003715808390000101
其中，Wr和br、Wz和bz及Wn和bn分别为重置门rt、更新门zt和候选隐藏层nt的权重矩阵和偏置向量。rt负责在单元的传递过程中控制需要保留的记忆，当rt趋近于1时，nt中主要包含当前输入状态xt和上一节点输入的隐藏状态ht-1共同组成的信息，当rt趋近于0时，则nt中主要包含xt的信息；而zt控制对ht-1进行选择性遗弃的程度，当zt趋近于1时，记忆单元主要输出上一时刻的ht-1，而当zt趋近于0时，记忆单元主要输出nt，即输出当前时刻的ht；
综上所述，GRU通过上一节点的隐藏状态ht-1和当前时刻的输入xt来获取两个关键的门控状态，而重置门和更新门的输出向量决定了哪些信息能够作为GRU的最终输出。因此，GRU能够对信息流进行精准的控制，具有保存长时间序列信息的能力，同时减少了训练参数，使得训练难度大大降低。为了加深对原始时间序列特征提取的层次，进一步提高模型输出的准确性，将两个方向不同且独立的GRU叠加在一起构成BiGRU网络，具体结构如图3所示。
将xt输入至前向层，从0时刻到t时刻正向计算出向前隐含层的输出hf；输入至反向层，从t时刻到0时刻反向计算出向后隐含层的输出hb。最后，将前向层和反向层的输出结果输入至全连接层，得到最终输出h，计算过程如式(5)所示：
h＝f(hf,hb) (5)
其中，f(·)是全连接层的映射函数。
三、对数正态分布
设备的可靠性通常使用概率指标进行度量，其寿命是一个非负变量，难以通过正态分布拟合其变化趋势。对数正态分布(lognorm，LN)能够将取值差异大的寿命数据取对数后转换为正态分布处理，广泛用于描述金属材料断裂、半导体元器件老化和绝缘材料失效过程，是可靠性领域的常用寿命分布之一。假设给定一个随机变量Y＝ln(X)，服从均值为μ、方差为σ的正态分布，则随机变量X满足参数为(μ,σ2)的LN分布，其各个统计特征的计算公式如表1所示。
表1对数正态分布的统计特征
Table 1 Statistical Characteristics of Lognormal Distribution
Figure BDA0003715808390000121
LN分布的统计特征值分布如图4中(a)所示，其中众数的PDF最大，同时LN分布呈现出右偏态，样本主要集中于期望值的左侧，而预测RUL应尽量小于真实RUL以保证设备的安全性要求，因此LN分布适合描述退化过程中RUL的概率分布。图4中(b)和(c)展示了不同μ和σ对PDF的影响，μ的值越小，LN分布的偏态特性越明显，而σ的值越小，LN分布越趋近于正态分布。综上，LN分布仅根据两个参数就可以改变成不同的形状，是一种较为灵活的分布，有利于描述大致对称或向右倾斜的寿命样本。
四、机械设备剩余寿命预测模型
在机械设备剩余寿命预测中，往往需要根据多方面的状态信息来评估其健康状况，因此多维历史时间序列是性能退化的主要表现，它包含了性能退化规律的重要信息。传统信息融合方法在处理多维传感器数据时，通常是根据单调性或趋势性人为选取与性能退化相关的最优特征，这种方法难以全面表征整个设备的健康状态，增加了信息丢失的风险，并且特征提取与随机模型是孤立进行的，导致RUL的预测精度较低。而FA机制无需人工干预就可以学习输入特征的动态变化，并赋予其不同的概率权重，从多维度、非线性和大规模的传感器信号中提取出关键特征向量；基于LN分布模型和BiGRU网络融合的预测模型能够从前向和后向两个方向对注意力加权特征的时变特性进行建模学习，并通过最大似然估计实现损失函数的最小化来训练网络参数，获得隐含状态输出向量的概率分布，从而计算出基于LN分布的PDF，实现部件不确定性的衡量。因此，本发明提出的FA-LN-BiGRU机械设备RUL预测模型，能够综合统计分析模型在预测不确定性量化能力上的优势与机器学习模型在大数据处理能力上的优势，有效挖掘传感器数据中丰富的特性与规律，提高RUL预测结果的精度，同时解决不确定性难以度量的问题，为后续管理辅助制定运维策略。
1、FA-LN-BiGRU模型结构
本发明构建的FA-LN-BiGRU网络模型结构如图5所示。该模型主要分为：输入层、特征加权层、特征学习层和输出层。将机械设备传感器数据输入至特征加权层，通过GAP、两层Linear及激活操作后获得注意力权重，将之与原始数据相乘实现关键特征的提取；特征学习层从注意力加权特征中学习机械设备随时间变化的规律；最后输出层对LN分布的参数进行最大似然估计，从而对当前RUL的PDF做出精准预测。下面对各主要模块说明：
(1)输入层
输入层负责将机械设备传感器数据输入至特征加权层中。记运行工况和不同类型传感器监测数据为多维特征时间序列矩阵X＝(x1,x2,…,xF)＝(x1,x2,…,xT)T，展开可表示为：
Figure BDA0003715808390000141
式中：
Figure BDA0003715808390000142
为在t时刻的F个特征数值序列，
Figure BDA0003715808390000143
为第f个多维特征在T个历史时刻的数值序列。
(2)特征加权层
特征加权层负责提取输入时间序列的注意力权重，从而获得注意力加权特征。利用FA机制自动提取输入数据(多维特征时间序列矩阵X)与目标建模参数(W1、W2、b1、b2等一系列神经网络参数)之间的关联关系，通过不断循环迭代对所有时刻各个特征分别分配注意力权重，逐渐提取出对RUL预测贡献率较大的关键特征。根据式(1)和式(2)，注意力权重可表示为：
Figure BDA0003715808390000144
其中，e＝[e1,e2,…,eF]∈RF×1；
则注意力加权特征为：
X′＝e⊙X＝[e1x1,e2x2,…,eFxF] (8)
其中，X′可记为
Figure BDA0003715808390000145
其中F为注意力加权特征的维度，T为时间长度；
(3)特征学习层
对于机械设备剩余寿命预测任务，作为输入的传感器数据特征是一段连续的具有强相关性的时间序列，从而要求网络具有一定的“记忆功能”，通过学习前后信息的差异对当前设备的剩余寿命进行判断。特征加权层可以对多维特征赋予不同的权重，但难以学习具有明显时间相关性的序列信息。因此，在特征加权层后搭建BiGRU网络构成特征学习层，使得预测模型能够捕获数据的时序依赖关系。
因时间序列长度为T，输入特征共F个，故BiGRU层由T个顺序串联的记忆单元组成，每个GRU单元的输入维度为F。如图5所示，序列数据将分时刻输入至前向层和反向层，得到t时刻前向输出hft和后向输出hbt，最终通过式(5)计算出单元的输出ht，t∈[1,T]，而BiGRU总输出H可记为[h1,h2,…,ht,…,hT]。
(4)输出层
将特征学习层的输出H作为输出层的输入，通过全连接层计算出设备的RUL的PDF参数预测值Y′＝[μ1,μ2,…,μt,…,μT；σ1,σ2,…,σt,…,σT]T，计算公式可表示为：
y′t＝[μt,σt]＝woht+bo (9)
其中，y′t为第t时刻模型的预测值；wo为权重矩阵；bo为偏置系数。
(5)损失函数
为准确拟合设备退化过程中RUL的概率分布，首先定义RUL预测的统计参数模型，即LN分布，从而对不确定性分布进行假设；再通过模型预测每个时刻特定分布的参数(对数均值μ和对数方差σ)；之后利用最大似然估计L(μ,σ|y)来训练BiGRU网络的参数，计算公式如式11所示：
Figure BDA0003715808390000151
其中，y是RUL的实测值。
对数似然函数常用于加速模型收敛至最优值，同时可以计算单个独立统计因子的总和，有助于损失函数的反向求导。其次，深度学习方法通常会最小化损失函数，因此将损失函数定义为：
Figure BDA0003715808390000161
对数均值μt由模型输出的映射函数直接求得，如式12所示。然而对于预测不确定性的衡量，对数方差σt应该是一个正值，因此由激活函数softplus得出，计算公式如式13所示。
μt＝wμyt (12)
Figure BDA0003715808390000162
其中，wμ和wσ分别代表与μ和σ相关的输出权重；yt为t时刻RUL的实测值；μt为t时刻的对数均值；σt为t时刻的对数方差。
实施例
为了进一步验证本发明所提方法的有效性进行对比实验。
1、实验数据集
本发明采用NASA提供的航空发动机仿真状态监测公开数据集—C-MAPSS，是RUL预测领域应用最为广泛的公共数据集之一。该数据集包含了不同故障模式和工作条件下涡扇发动机从正常运行到故障失效的4组传感器监测数据(FD001-FD004)，且将21种传感器作为能够表征发动机运行状况的典型指标。每组数据由训练集、测试集和RUL标签组成，训练集包含发动机整个寿命周期的全部数据，测试集只包含测试发动机生命周期前一段的数据，RUL标签是与测试集相对应的发动机最后监测时刻的RUL。其中各组数据包含不同数量的发动机，每台发动机带有不同程度的初始磨损，因此每台发动机监测数据的序列长度也不同。4组监测数据的具体信息如表2所示。
表2航空发动机数据集
Table 2 CMAPSS dataset
Figure BDA0003715808390000171
2、数据预处理
(1)特征选择
以FD001为例说明特征选择的过程，仅以展示训练集前12个传感器数据特征的分布为例，如图6所示，其中深蓝色点表示退化数据，红色点表示最终失效点，横坐标为发动机的全寿命周期内运行的时间，纵坐标为不同传感器的测量值。由图6可知，T2、P2、P15、epr传感器的方差为0，表明这些传感器的监测数据在整个生命周期中并未产生显著变化，因此这些特征难以提供任何对RUL预测有贡献的退化信息，应予以剔除。通过对所有的特征进行分析，最终选取14个传感器数据与3种运行状态作为输入特征的数据集，名称分别为：set1、set2、set3、T24、T30、T50、P30、Nf、Nc、Ps30、Phi、NRf、NRc、BPR、htBleed、W31和W32。
(2)数据归一化及RUL标签设置
由于多维监测数据具有不同的量纲，所以在模型构建前须进行归一化预处理以消除特征之间的影响。本发明选用最小-最大标准化方法将数据统一至[0,1]范围内，具体归一化公式如下：
Figure BDA0003715808390000181
其中，
Figure BDA0003715808390000182
为归一化后的数据；
Figure BDA0003715808390000183
和
Figure BDA0003715808390000184
分别为第j个传感器监测数据的最大值和最小值。
其次，发动机在运行初期状态良好，性能退化量可忽略不计，但到了使用末期，发动机的性能会随着时间的延长而急剧下降。如果将发动机性能迅速劣化之前监测数据的标签设置为总运行周期减去当前运行周期数，会增加RUL预测结果的滞后性。因此可以认为发动机开始急速退化之前的RUL保持不变，即为训练集的RUL标签设置阈值，使之成为分段线性函数。研究表明，将训练集RUL标签突变临界值设置为第130个运行周期预测效果较好，RUL标签设置结果如图7所示。
(3)滑动窗口处理
为了利用有限的时间序列深入挖掘性能退化规律，同时将时间序列转化成BiGRU网络擅长处理的三维输入格式，采用滑动时间窗分割方法对归一化后的数据进行处理。这样既可以充分保留相邻序列之间的时间相关性，又可以提高训练集样本的数量，使模型具有更高的鲁棒性和泛化能力。记发动机原始时间序列长度为T，特征维度为N，利用窗宽为S的滑动窗沿着时间序列滑动，将每滑动一个步长所截取的时间序列叠加至第三个维度，形成(T-S,S,N)的三维张量，计算过程如下：
Figure BDA0003715808390000185
其中，X1:T-S为转换后的三维张量；xi:i+S为从第i个时间周期开始长度为S的序列；
Figure BDA0003715808390000186
表示窗口内数据在第三个维度的连接，从而形成三维张量。
3、模型超参数设置
本模型所涉及的主要超参数有：网络结构、学习率、批尺寸、迭代次数、优化器、Dropout率等，由PyTorch框架实现。模型超参数对模型性能的优劣有较大的影响，因此采用网格搜索法来使得模型预测误差最低，从而获取最优超参数组合，具体如表3所示。
表3模型最优超参数组合
Table 3 Model optimal hyperparameter combination
Figure BDA0003715808390000191
4、评价指标
本发明使用两方面的评价指标来量化预测模型的性能：点预测评价指标和区间预测评价指标。其中，点预测评价指标用以评估模型在测试集上进行确定性点预测的准确度和泛化能力，区间预测评价指标用以衡量预测区间的质量和覆盖率。
(1)点预测评价指标
①均方根误差(rootmean square error,RMSE)
RMSE用于衡量预测值同真实值之间的偏差，其计算公式为：
Figure BDA0003715808390000192
其中：yi为第i个发动机的实际RUL值；
Figure BDA0003715808390000193
为第i个发动机的RUL预测值；m为发动机总数。
②评分函数
由于发动机较高的安全性要求，不及时维修造成的灾难性后果所引起的成本远大于过度维修资源的成本，因此早期预测比滞后预测更可取。SF作为2008年PHM数据挑战赛中官方提出的权威指标，对RUL滞后预测的情况施加更高的惩罚。其计算公式如下：
Figure BDA0003715808390000201
其中：a1＝13，a2＝10，a1和a2控制SF对预测滞后性的惩罚程度，从而实现非对称预测偏好。另外，SF分数越低，表明模型的预测性能越好。
(2)区间预测评价指标
①预测区间覆盖率(predict interval coverage probability，PICP)
PICP反映实际RUL值yi落入预测值
Figure BDA0003715808390000202
波动区间
Figure BDA0003715808390000203
的概率，用于评估预测区间的可靠性，计算公式如下：
Figure BDA0003715808390000204
其中：
Figure BDA0003715808390000205
与
Figure BDA0003715808390000206
分别为
Figure BDA0003715808390000207
在额定置信水平(1-α)下预测区间的上下边界。PICP越大表示有越多的真实值落入预测区间内，模型的预测结果越可靠。其次，PICP应尽可能超过给定的置信度，否则在实际维修决策中将被判定为预测无效。
②预测区间平均带宽(predictinterval normalized average width，PINAW)
仅追求PICP会导致预测区间过宽，难以精确获得RUL预测的不确定性信息，因此还需要考虑区间宽度以提高决策价值。PINAW表示预测区间宽度的平均值，可以用来评估预测区间的清晰度，其计算公式为：
Figure BDA0003715808390000211
其中：ymax和ymin分别为实际RUL的最大值和最小值。当PICP一定时，越小的PINAW对应越窄的预测区间，模型的预测效果越好。
5、模型结果
(1)可视化分析
为了验证本方法点预测的有效性，以FD001数据集为例，将预处理后的高维时间序列输入至FA-LN-BiGRU模型以进行RUL的点预测，测试集的预测结果如图8(a)所示。为比较RUL大小对模型预测的影响，将FD001测试集中所有的发动机根据真实RUL值从大到小进行重新排序，结果如图8(b)所示。由图8(a)可知，RUL的真实值和预测值非常接近，表明了本方法点预测的有效性。当RUL值较大时，发动机处于健康状态，表明设备运行状况良好。由图8(b)可知，发动机刚开始运行时，RUL预测误差较大，同时波动较为剧烈且表现出明显的滞后性。但经过长时间运行后，预测RUL收敛在真实RUL周围，预测性能显著增强。因此，设备历史信息越充分，性能退化信息越明显，模型的预测误差越小。
为了进一步探究本方法区间预测的有效性，在数据集中随机选取3个引擎在80％～95％置信度的RUL区间预测结果与各个引擎最终监测点RUL的pdf曲线，如图9所示。其中，RUL点预测结果基于表1中位数求得。由图9可以看出FA机制使得模型专注于包含引擎关键退化信息的重要特征，即使引擎刚开始运转时历史数据较少，难以预测RUL大小，模型预测值也较为接近临界值130。其次，RUL点预测结果与真实值基本吻合，80％置信度的区间预测基本能覆盖真实RUL，另外区间宽度在引擎退化引起的点估计波动期间变宽，但随着真实RUL的变小逐渐变窄，且后期点预测精度较高，这是因为随着运行周期的增加，发动机的性能退化量逐渐积累，模型能够逐渐捕获时间序列前后的退化信息，有效提升运行后期的预测精度。
(2)对比分析
为了研究本文方法进行发动机RUL点预测和区间预测的预测性能，采用上述点预测评价指标以及区间预测评价指标来客观评估不同模型在测试集上的准确度和泛化能力。表4对比了在单一模式下FD001数据集和多种模式下FD004数据集，利用不同循环神经网络的训练效果，同时对比了不同种类的注意力机制:自注意力机制(self-attention,SA)和FA机制的特征提取效果。
由表4可知，与其余5种方法相比，本发明的RMSE和SF指标均为最小，对于FD001数据集，RMSE＝13.13，与其他方法相比平均降低了11％，SF＝238，平均降低了29％；对于FD004数据集，RMSE＝20.29，与其他方法相比平均降低了19％，SF＝3395，平均降低了36％。其次，由于LN分布的不对称性，使用中位数、期望值和众数来评估点预测精度时存在细微差异。相较于另外两个统计特征，期望值为RMSE提供了更好的结果，但考虑到滞后预测的危害及综合RMSE值的大小，本发明选取SF值更小的中位数作为点预测的结果。进一步分析可知，虽然GRU训练效果略逊于LSTM，但训练时间缩短约11％，有利于实时预测的部署；其次BiGRU能够捕捉双向长时间依赖特性，较单向网络FD001数据的RMSE和score分别下降约4％和10％，FD004数据的RMSE和score更是下降约7％和16％，达到预测精度更优的效果；FA机制的概率权重不同于SA机制的加权归一化权重，可以对不同特征进行重要性打分，而非对某一重要特征进行筛选，进一步提升了模型预测精度。因此从点预测结果纵向对比可以发现，本发明所提方法具有更优的实时性和泛化能力。
表4不同预测模型点预测对比结果
Table 4 Comparison results of point forecasts of differentforecasting models
Figure BDA0003715808390000231
为了进一步验证本文方法进行区间预测的优势，分别构建了蒙特卡洛模型、贝叶斯模型和LSTMBS模型作为对比模型，所得FD001数据集预测结果如表5所示。
表5不同预测模型预测对比结果
Table 5 Comparing the results of different forecasting models
Figure BDA0003715808390000232
如表5所示，蒙特卡洛模型在点预测的效果最好，但所使用的Dropout操作难以描述设备性能退化的不确定性，过于追求神经网络的性能而导致置信区间的平均宽度过窄，PICP难以满足对应置信度的要求。贝叶斯模型作为传统的统计分析模型，面对高维特征时稍显乏力，大数据的波动性和稀疏性导致PINAW过宽。LSTMBS模型区间覆盖率PICP与本文方法相差不大，但是平均带宽PINAW在置信度为95％、90％和80％下分别缩小19.8％、13.9和16.9％，说明本方法可以更准确地量化性能退化的不确定性，有效提高区间预测的预测精度，另外点预测精度也优于LSTMBS模型，表明本文方法具有更强大的特征提取能力和数据处理能力。综上所述，通过表5横向对比传统统计分析方法和数据驱动融合方法发现，本方法能够从高维、稀疏的大数据中提取更为准确的性能退化趋势特征，同时有效衡量设备退化过程中的不确定性，实现RUL的精准预测，其预测区间兼具可靠性高和准确性好两方面的优点，为后续运维人员制定维修决策提供有力支持。
以上显示和描述了本发明的基本原理、主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。

Claims (8)

1.一种机械设备剩余寿命区间预测方法，其特征在于，包括以下步骤：
步骤1：获取机械设备的传感器监测数据；
步骤2：构建FA-LN-BiGRU预测模型；
步骤3：将得到的传感器数据输入至FA-LN-BiGRU预测模型中，通过所述预测模型计算并输出航空发动机的RUL的PDF参数预测值；
步骤4：根据PDF参数预测值对发送机的剩余寿命区间进行预测。
2.如权利要求1所述的一种机械设备剩余寿命区间预测方法，其特征在于，所述FA-LN-BiGRU预测模型包括依次连接的输入层、特征加权层、特征学习层和输出层；
输入层，用于将得到的航空发动机的历史传感器监测数据输入到特征加权层；
特征加权层，用于通过轻量级特征注意力机制提取输入的传感器监测数据的注意力加权特征向量；
特征学习层，用于通过双向BiGRU网络提取注意力加权特征向量的时序依赖关系；
输出层，用于计算航空发动机的RUL的PDF参数预测值。
3.如权利要求2所述的一种机械设备剩余寿命区间预测方法，其特征在于，所述特征加权层包括全局平均池层、两层全连接层和sigmoid激活函数；
全局平均池化层，用于将获取到的传感器监测数据进行平均池化操作，得到每个时间样本的均值信息；
全连接层，用于量化各个特征对应的注意力权重；
sigmoid激活函数，用于注意力权重做归一化处理。
4.如权利要求3所述的一种机械设备剩余寿命区间预测方法，其特征在于，步骤3的具体操作步骤包括：
步骤31：输入层获取不同运行工况和类型的机械设备历史传感器数据，将所述传感器数据表示为多维特征时间序列矩阵：
Figure FDA0003715808380000021
其中，
Figure FDA0003715808380000022
为在t时刻的F个特征数值序列，
Figure FDA0003715808380000023
为第f个多维特征在T个历史时刻的数值序列；
步骤32：所述特征加权层利用特征注意力机制提取多维特征时间序列与目标建模参数之间的关联关系，通过循环迭代对所有时刻各个特征分别分配注意力权重，得到输入的多维特征时间序列的注意力加权特征；
步骤33：特征学习层将输入的注意力加权特征向量分时刻输入至双向BiGRU网络的前向层和反向层，得到t时刻前向输出和后向输出，再将前向层和反向层的输出合并得到总输出；
步骤34：输出层根据得到的特征学习层的总输出，实现特征学习层输出向量对最终参数输出的映射，得到当前时刻下的对数均值μt和对数方差σt，从而计算此时RUL的PDF；
步骤35：构建最终参数与目标输出的损失函数，通过损失函数更新网络参数，再通过更新后的网络循环迭代求解下一时刻的对数均值和对数方差，最终实现对数正态分布参数的最大似然估计。
5.如权利要求2所述的一种机械设备剩余寿命区间预测方法，其特征在于，步骤32的具体操作步骤包括：
步骤321：将输入的多维特征时间序列x通过全局平均池进行池化操作，得到每个时间样本的均值信息，将维度为T×F的多维特征时间序列降维至1×F，该计算过程记为：
Figure FDA0003715808380000031
其中，yf表示通过全局平均池化后得到的向量；g(·)为特征维度的全局平均池化操作；Xf为第f个特征所有时间点的数据；
Figure FDA0003715808380000032
为第f个特征的第t个时间点的数据；
步骤322：通过两层全连接层量化各个特征对应的注意力权重，并通过sigmoid激活函数对权重做归一化处理，得到归一化后的注意力权重αf：
αf＝σ(W2δ(W1yf+b1)+b2) (2)
根据式(1)将式(2)改写为：
Figure FDA0003715808380000033
其中，e＝[e1,e2,…,eF]∈RF×1；
步骤323：将注意力权重与相应的多维特征时间序列x相乘，得到多维特征时间序列的注意力加权特征X′：
X′＝e⊙X＝[e1x1,e2x2,…,eFxF] (8)
其中，X′可记为
Figure FDA0003715808380000034
其中F为注意力加权特征的维度，T为时间长度。
6.如权利要求4所述的一种机械设备剩余寿命区间预测方法，其特征在于，步骤33的具体操作步骤包括：
步骤331：将1时刻得到的注意力加权特征同时输入至前向层和后向层，直至T时刻；
步骤332：分别计算从0时刻到t时刻正向计算出向前隐含层的输出hf，以及从t时刻到0时刻反向计算出向后隐含层的输出hb；
步骤333：将hf和hb通过式(5)计算得到最终特征学习层的总输出H
h＝f(hf,hb) (5)
其中，f(·)是全连接层的映射函数。
7.如权利要求6所述的一种机械设备剩余寿命区间预测方法，其特征在于，步骤34所述计算RUL的PDF的公式为：
y′＝[μt,σt]＝woH+bo (9)
其中，y′t为第t时刻模型的预测值；μt为当前时刻的对数均值，σt为当前时刻的对数方差；wo为权重矩阵；bo为偏置系数。
8.如权利要求7所述的一种机械设备剩余寿命区间预测方法，其特征在于，步骤35的具体操作步骤包括：
步骤351：计算t0时刻对数均值和对数方差：
Figure FDA0003715808380000041
Figure FDA0003715808380000042
其中，wμ和wσ分别代表与μ和σ相关的整体网络输出权重；xt为t0时刻传感器监测数据；
Figure FDA0003715808380000043
为t0时刻的对数均值；
Figure FDA0003715808380000044
为t0时刻的对数方差；
步骤352：定义RUL预测的统计参数模型：给定一个随机变量Y＝ln(X)，服从均值为μ、方差为σ的正态分布，所述随机变量X满足参数为(μ,σ2)的LN分布；
步骤353：通过所述统计参数模型预测每个时刻分布的参数对数均值μ和对数方差σ；
步骤354：利用最大似然估计L(μ,σ|y)训练FA-LN-BiGRU网络的参数，且最大似然估计的计算公式为：
Figure FDA0003715808380000051
其中，y是RUL的实测值；
步骤355：定义FA-LN-BiGRU网络的损失函数，并利用反向传播算法更新网络参数：
Figure FDA0003715808380000052
CN202210736162.6A 2022-06-27 2022-06-27 一种机械设备剩余寿命区间预测方法 Pending CN115204362A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210736162.6A CN115204362A (zh) 2022-06-27 2022-06-27 一种机械设备剩余寿命区间预测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210736162.6A CN115204362A (zh) 2022-06-27 2022-06-27 一种机械设备剩余寿命区间预测方法

Publications (1)

Publication Number Publication Date
CN115204362A true CN115204362A (zh) 2022-10-18

Family

ID=83578693

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210736162.6A Pending CN115204362A (zh) 2022-06-27 2022-06-27 一种机械设备剩余寿命区间预测方法

Country Status (1)

Country Link
CN (1) CN115204362A (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116502544A (zh) * 2023-06-26 2023-07-28 武汉新威奇科技有限公司 一种基于数据融合的电动螺旋压力机寿命预测方法及系统
CN118332354A (zh) * 2024-06-12 2024-07-12 之江实验室 考虑设备结构信息的剩余使用寿命预测方法、装置及介质

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116502544A (zh) * 2023-06-26 2023-07-28 武汉新威奇科技有限公司 一种基于数据融合的电动螺旋压力机寿命预测方法及系统
CN116502544B (zh) * 2023-06-26 2023-09-12 武汉新威奇科技有限公司 一种基于数据融合的电动螺旋压力机寿命预测方法及系统
CN118332354A (zh) * 2024-06-12 2024-07-12 之江实验室 考虑设备结构信息的剩余使用寿命预测方法、装置及介质
CN118332354B (zh) * 2024-06-12 2024-08-20 之江实验室 考虑设备结构信息的剩余使用寿命预测方法、装置及介质

Similar Documents

Publication Publication Date Title
US11409347B2 (en) Method, system and storage medium for predicting power load probability density based on deep learning
CN109492822B (zh) 空气污染物浓度时空域关联预测方法
CN109685252B (zh) 基于循环神经网络和多任务学习模型的建筑能耗预测方法
CN112990556A (zh) 一种基于Prophet-LSTM模型的用户用电能耗预测方法
CN115204362A (zh) 一种机械设备剩余寿命区间预测方法
CN112116162B (zh) 基于ceemdan-qfoa-lstm的输电线覆冰厚度预测方法
CN111340282B (zh) 基于da-tcn的设备剩余使用寿命的估计方法及系统
CN112785092B (zh) 一种基于自适应深层特征提取的道岔剩余寿命预测方法
CN115688579A (zh) 一种基于生成对抗网络的流域多点水位预测预警方法
CN112329990A (zh) 一种基于lstm-bp神经网络的用户用电负荷预测方法
CN113128113A (zh) 一种基于深度学习和迁移学习的贫乏信息建筑负荷预测方法
CN112070322B (zh) 基于长短期记忆网络的高压电缆线路运行状态预测方法
Barzola-Monteses et al. Energy consumption of a building by using long short-term memory network: a forecasting study
CN114580545A (zh) 一种基于融合模型的风电机组齿轮箱故障预警方法
CN117154704B (zh) 一种基于多尺度时空图注意力卷积网络的光伏功率预测方法
CN117554818A (zh) 一种锂离子电池剩余使用寿命预测方法
CN114357670A (zh) 一种基于bls和自编码器的配电网用电数据异常预警方法
CN115840893A (zh) 一种多变量时间序列预测方法及装置
CN117293790A (zh) 一种考虑预测误差不确定性的短期电力负荷预测方法
CN111754033B (zh) 一种基于循环神经网络的非平稳时序数据预测方法
CN113673774A (zh) 基于自编码器和时序卷积网络的航空发动机余寿预测方法
CN116845882A (zh) 考虑温湿指数与耦合特性的综合能源负荷短期预测方法
CN117371331A (zh) 基于kddin的航空发动机剩余寿命预测方法
CN117236380A (zh) 一种电力系统故障预测方法、系统、电子设备和介质
CN117313936A (zh) 一种燃煤电厂烟气脱硫过程净烟气so2浓度预测方法

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