CN111010193B - 基于压电导纳的数据压缩重构及结构损伤识别方法、设备 - Google Patents

基于压电导纳的数据压缩重构及结构损伤识别方法、设备 Download PDF

Info

Publication number
CN111010193B
CN111010193B CN201911375347.3A CN201911375347A CN111010193B CN 111010193 B CN111010193 B CN 111010193B CN 201911375347 A CN201911375347 A CN 201911375347A CN 111010193 B CN111010193 B CN 111010193B
Authority
CN
China
Prior art keywords
sparsity
data
matrix
reconstruction
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
CN201911375347.3A
Other languages
English (en)
Other versions
CN111010193A (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.)
Huazhong University of Science and Technology
Original Assignee
Huazhong University of Science and 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 Huazhong University of Science and Technology filed Critical Huazhong University of Science and Technology
Priority to CN201911375347.3A priority Critical patent/CN111010193B/zh
Publication of CN111010193A publication Critical patent/CN111010193A/zh
Application granted granted Critical
Publication of CN111010193B publication Critical patent/CN111010193B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • H—ELECTRICITY
    • H03—ELECTRONIC CIRCUITRY
    • H03M—CODING; DECODING; CODE CONVERSION IN GENERAL
    • H03M7/00—Conversion of a code where information is represented by a given sequence or number of digits to a code where the same, similar or subset of information is represented by a different sequence or number of digits
    • H03M7/30—Compression; Expansion; Suppression of unnecessary data, e.g. redundancy reduction

Abstract

本发明属于结构健康监测领域，涉及基于压电导纳的数据压缩重构及结构损伤识别方法、设备。该方法通过采集原始电导数据，并对其作最大最小归一化处理；利用压缩感知理论对处理后的采样数据进行压缩，通过随机矩阵对数据实施线性变换，将恢复原始信号的问题转化为残差迭代问题；基于优化稀疏度参数分析的自适应正交匹配追踪算法对电导数据进行重构；利用重构数据得到的均方根偏差和相关系数对结构进行损伤诊断。该方法能够有效地将基于压电导纳技术的损伤识别与基于正交匹配追踪算法的数据压缩算法结合，能够有效识别结构不同类型的损伤，并且在传感器带伤作业的情况下依然有效。

Description

基于压电导纳的数据压缩重构及结构损伤识别方法、设备
技术领域
本发明属于结构健康监测领域，涉及一种基于压电导纳的数据压缩重构及结构损伤识别方法、设备，更具体地，涉及一种结合运用数据压缩和重构理论与压电导纳技术的损伤识别方法。
背景技术
由于土木工程结构处于复杂多变的环境下作业，受自身材料退化、环境腐蚀、外力干扰(如风、地震)等耦合因素影响，造成结构或构件的损伤。其中绝大部分损伤属于一个渐变的过程，期间引发的微小损伤通过肉眼难以识别，因而结构早期的微小损伤逐渐积累而诱发大的隐患至今都是土木工程领域面对的难题。基于传统压电阻抗技术的损伤识别系通过将PZT传感器粘结于结构表面或植入结构内部，设置局部高频激励，利用PZT正逆压电效应，以PZT自身同时作为驱动器和传感器，获取结构局部性能(质量，刚度，阻尼等)相关信息，一般以初始(未损)情形为基线状态，再根据同样方法获取待测(未知)状态的相关信息，以比较其与基线状态的不同而识别结构局部性能变化的信息，从而实现局部微小损伤的识别。
然而其缺陷在于(并不限于)：监测过程中所需采样的数据量极为庞大，数据采集系统储存负荷过重，不利于试验/工程人员进行实时数据收集和处理分析；且采样过程易受到环境、人为干扰因素强烈的影响，导致识别效果大幅降低。
发明内容
针对传统技术的缺陷或改进需求，本发明提供了一种基于压电导纳的数据压缩重构及结构损伤识别方法、设备，其目的在于，对于压电导纳数据中的电导信号，通过将基于正交匹配追踪算法的压缩感知理论引入到基于压电导纳技术的信号压缩重构与损伤识别过程中，可以对复杂庞大的信号进行压缩采样并且重构后不失真，从而降低数据采样、存储的成本，并提高结构损伤识别的效率。
为了实现上述目的，按照本发明的一个方面，提供了一种基于压电导纳的数据压缩重构方法，包括如下步骤：
(1)基于压电导纳技术，采集结构初始状态时的导纳数据，将原始电导信号作如下最大最小归一化处理：
Figure GDA0003118647020000021
其中，
Figure GDA0003118647020000022
为原始电导信号，
Figure GDA0003118647020000023
为
Figure GDA0003118647020000024
的最小值，
Figure GDA0003118647020000025
为
Figure GDA0003118647020000026
的最大值，X为处理后的原始电导信号，长度为m；
(2)设定稀疏度区间范围及步长[Sstart：Sstep：Send]，其中，Sstart为所选稀疏度最小值，Sstep为所选稀疏度步长，Send为所选稀疏度最大值，从而得到一系列稀疏度Si，下标i表示稀疏度序号；
(3)输入当前稀疏度Si，下标i表示第i个稀疏度，设置相对应的测量数n，
n＝cst×ceil(Si/log(m/Si))
设置转换基Ψ(m×m)和测量矩阵Φ(n×m)，下标m×m和n×m表示矩阵行列数，基于公式Y＝ΦX＝ΦΨα＝Θα，对信号X进行基于压缩感知理论与正交匹配追踪算法的数据压缩与重构，包括：
即将传感矩阵Θ作为完备字典库，对观测向量Y进行分解，从传感矩阵Θ中选择一个与信号Y最匹配的即列向量，构建稀疏系数逼近Y，并求出信号残差，然后继续选择与信号残差最匹配的列向量，再次进行残差分解，反复迭代，最终信号Y可以被传感矩阵Θ中的所有列向量表达，这些列向量也在迭代过程中被重新排列，按照每次迭代获得的稀疏系数的顺序，将所有稀疏系数排列组成稀疏系数矩阵α，稀疏系数矩阵α即为压缩数据，对应的重构数据为
Figure GDA0003118647020000031
(4)计算当前稀疏度Si对应的重构数据
Figure GDA0003118647020000032
与原始数据X比对下得到的总误差与压缩率；
(5)回到步骤(3)，输入稀疏度Si+1，重复步骤(3)和(4)，直到遍历稀疏度区间范围内的所有稀疏度；Si+1＝Si+Sstep；
(6)根据步骤(5)获得的所有稀疏度对应的总误差及压缩率，绘制稀疏度与总误差的关系曲线、稀疏度与压缩率的关系曲线，基于上述关系曲线，找出保证总误差低于期望值的同时尽量使得压缩率最高的稀疏度，作为最优稀疏度，并得到当前最优稀疏度对应的稀疏系数矩阵和重构数据，完成数据压缩与重构。
进一步地，步骤(2)设定稀疏度区间范围及步长[Sstart：Sstep：Send]的方法如下：
从一长度不大于信号X长度的几何级数序列中，以信号X长度的十分之一为基，截取一系列几何级数作为稀疏度区间进行试算，试算方法如下：
根据试错法，根据截取出的稀疏度区间生成随机测量矩阵
Figure GDA0003118647020000033
代入
Figure GDA0003118647020000034
并计算，挑选出一个使得上式值尽量小的稀疏度以确定稀疏度的区间范围[Sstart：Sstep：Send]；
其中，
Figure GDA0003118647020000035
为试算时设定的随机测量矩阵，Y为观测向量，|| ||2表示l2范数，Y＝ΦX，Φ为不同于
Figure GDA0003118647020000036
的随机测量矩阵。
进一步地，步骤(3)包括如下子步骤：
(3.1)基于压缩感知理论将X进行线性映射，包括：
引入一个m×m的转换基矩阵Ψ＝[ψ1，ψ2，…，ψm]，其中ψ1～ψm均为列向量，具有m列，将信号(离散序列)X通过正交基Ψ变换为基系数矩阵α，有：X＝Ψα；
又由于Y＝ΦX，因此有Y＝ΦX＝ΦΨα＝Θα，其中Θ为传感矩阵；
(3.2)基于正交匹配追踪理论的稀疏度自适应正交匹配追踪算法，对步骤(2)所得到的稀疏度区间[Sstart：Sstep：Send]中每一个稀疏度，遍历传感矩阵Θ，从而求解每个稀疏度对应的α以及每个稀疏度对应的重构数据。
进一步地，步骤(3.2)包括如下子步骤：
(3.2.1)输入测量矩阵Φ(n×m)、转换基矩阵Ψ(m×m)、观测向量Y(n×1)和当前稀疏度Si；
(3.2.2)设定k＝0，初始残差R0＝Y，D0＝{ }，I0＝{1，…，m}，I0为传感矩阵Θ中列向量的序号构成的初始序号矩阵；
(3.2.3)k≥1时，Y的第k步分解如下：
Figure GDA0003118647020000041
其中，p＝0，1，…，k-1，第k阶模型为
Figure GDA0003118647020000042
Figure GDA0003118647020000043
为分解系数；Rk为第k步的当前残差，满足：
Figure GDA0003118647020000044
Ik-1＝Ik-2\{rk-2}
Figure GDA0003118647020000045
即，第k步分解时，遍历传感矩阵Θ中序号属于集合Ik-1的列向量，从中找出一个列向量dj，使得当前残差Rk和该列向量dj能够取得最大内积，将该列向量dj作为
Figure GDA0003118647020000046
加入到矩阵Dk中，然后将
Figure GDA0003118647020000047
在传感矩阵Θ中对应的列向量的序号从集合Ik-1中剔除；
(3.2.4)令k＝k+1，重复步骤(3.2.3)，直至k＝m时，传感矩阵Θ中的所有列向量遍历完毕，终止迭代，得到Dk，Dk为传感矩阵Θ重排列后的结果；
(3.2.5)根据公式(2)，在进行Y的第k步分解时，取k阶残差迭代模型为：
Figure GDA0003118647020000051
则k+1阶残差迭代模型减去k阶残差迭代模型得：
Figure GDA0003118647020000052
引入辅助模型
Figure GDA0003118647020000053
Figure GDA0003118647020000054
Figure GDA0003118647020000055
其中，
Figure GDA0003118647020000056
为第k-1阶分解系数，q＝0，1，…，k-1，γk为残差项；
由于
Figure GDA0003118647020000057
为了计算
Figure GDA0003118647020000058
将式(5)重写为一个线性系统：
uk＝Ckβk (7)
其中：
Figure GDA0003118647020000059
Figure GDA00031186470200000510
Figure GDA00031186470200000511
Ck满足非奇异性，因而式(7)可改写为
Figure GDA00031186470200000512
其中，
Figure GDA00031186470200000513
为Ck的逆矩阵，
由于Ck+1可以改写为：
Figure GDA00031186470200000514
其中*表示共轭转置，因而可以得到：
Figure GDA0003118647020000061
其中
Figure GDA0003118647020000062
因而可以通过递归计算
Figure GDA0003118647020000063
得到βk+1；
(3.2.6)将式(4)代入式(3)，得到：
Figure GDA0003118647020000064
假如式(15)和式(16)成立，则式(14)成立，可得：
Figure GDA0003118647020000065
Figure GDA0003118647020000066
令
Figure GDA0003118647020000067
由(7)左右两边对
Figure GDA0003118647020000068
作内积，得：
Figure GDA0003118647020000069
其中，|| ||2表示向量作差之后各分量的平方和；
Figure GDA00031186470200000610
(3.2.7)令k＝k+1，重复步骤(3.2.5)～(3.2.6)，直至k＝m时，求解出所有的
Figure GDA00031186470200000611
将所有
Figure GDA00031186470200000612
按照对应的
Figure GDA00031186470200000613
在Dk中的顺序，组成稀疏系数矩阵α，α即为
Figure GDA00031186470200000614
和X的压缩数据；
(3.2.8)输出当前稀疏度Si对应的稀疏系数矩阵α和重构数据
Figure GDA00031186470200000615
进一步地，步骤(4)中，总误差TE与压缩率CR计算如下：
CR＝Sorig/scomp
Figure GDA00031186470200000616
其中，Sorig表示初始数据
Figure GDA00031186470200000617
所占字节，Scomp表示压缩数据α所占字节。进一步地，步骤(6)中，根据稀疏度与总误差、稀疏度与压缩率的对应关系曲线，寻找总误差最低的3～4个值所对应的稀疏度，再比较上述稀疏度对应的压缩率大小，压缩率不小于2，从而筛选出最优稀疏度。
为了实现上述目的，本发明还提供了一种基于压电导纳的数据压缩重构设备，包括压电导纳信号采集模块、处理器和数据压缩重构程序模块；
所述数据压缩重构程序模块在被所述处理器调用时，按照如前任意一项所述的数据压缩重构方法对压电导纳信号采集模块采集的电导信号，进行压缩重构处理。
为了实现上述目的，本发明还提供了一种结构损伤识别方法，在待测结构上布置多个采样点，基于压电导纳技术采集各个采样点位置的电导纳信号，并对其中的电导信号部分，按照如前任意一项所述的数据重构方法筛选出最优稀疏度，并将最优稀疏度所对应的重构数据
Figure GDA0003118647020000071
与基线电导数据进行均方根偏差和相关系数计算，得出工况指标RMSD和CC：
Figure GDA0003118647020000072
Figure GDA0003118647020000073
式中，
Figure GDA0003118647020000074
表示第i个采样点下的基线电导信号值，
Figure GDA0003118647020000075
表示相应采样点的待测电导信号值，σ0和σ1分别表示基线电导信号和待测电导信号的标准差，Cov(*)为协方差运算符，
Figure GDA0003118647020000076
RMSD指标越大，则结构发生损伤的可能性越大，CC指标越接近1，则结构健康状况越好。
为了实现上述目的，本发明还提供了一种结构损伤识别设备，包括压电导纳信号采集模块、处理器、数据压缩重构程序模块以及结构损伤识别程序模块；
所述数据压缩重构程序模块在被所述处理器调用时，按照如前任意一项所述的数据压缩重构方法对压电导纳信号采集模块采集的电导信号，进行压缩重构处理，得到最优稀疏度及其对应的重构数据
Figure GDA0003118647020000077
所述结构损伤识别程序模块在被所述处理器调用时，按照如前所述的结构损伤识别方法，计算出工况指标RMSD和CC。
总体而言，通过本发明所构思的以上技术方案与现有技术相比，能够取得下列有益效果：
(1)本发明通过最大最小归一化处理简化原始电导数据，然后基于压缩感知理论及正交匹配追踪算法对简化后的原始电导信号进行稀疏度选取和压缩重构，可以大幅减少采样数据，重构出的信号还原度高，仍可准确携带原始信号的典型特征信息，从而能够用于损伤状态的精确识别。
(2)通过试错法先进行稀疏度区间及步长试算，可以大幅提高稀疏度的选取效率，提升数据压缩速度。
(3)通过对原始信号的线性映射、基于正交匹配的残差分解，对传感矩阵进行重排得到最逼近于观测向量的传感矩阵原子顺序，从而得到按照相应顺序排列的观测向量分解系数，组成稀疏度系数矩阵，获得当前稀疏度下的压缩及重构数据，可以在保证数据重构误差和识别精度的基础上，有效地减少数据采样所需的成本，为实时数据监测、处理与识别提供了有效的保障；
(4)本发明提供的压缩重构压电导纳数据的方法可以与其它技术如PZT无线传感技术有机结合，将本发明所提出的算法流程嵌入模块中，能够使得数据采集和处理系统更加方便，并且应用在实际工程中。经测试，本发明提出的方法及设备，对于金属结构的健康监测有着良好的识别效果，并且在处理传感器自身存在轻微损伤的情形下，也能够良好地获得有效的损伤评估指标值。
附图说明
图1是钢梁试验传感器布置和工况设置示意图；
图2(a)～图2(i)是试验钢梁的原始电导信号对比图；
图3(a)～图3(i)是预设不同稀疏度时对应重构数据与原始数据之间的总误差关系图；
图4(a)～图4(c)是各个工况RMSD指标图；
图5(a)～图5(c)是各个工况CC指标图。
图6是本发明的组成模块示意图；
图7是图6中数据采集模块的处理流程图；
图8是图6中数据压缩与重构模块的处理流程图；
图9是图6中损伤识别模块的处理流程图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。
本发明提供的方法，系一种将数据压缩重构算法与压电导纳技术相结合以进行结构健康监测的方法。首先，在建立的结构健康监测系统中采集与主结构相连的PZT传感器原始导纳信号(导纳信号包括电导信号和电纳信号，本发明在进行数据处理的过程中采用的是电导信号)，对原始电导数据作最大最小归一化处理；然后利用压缩感知理论对采样数据进行压缩，通过随机矩阵对数据实施线性变换，将恢复原始信号的欠定多项式问题转换成贪婪算法来求解；基于压缩感知理论，利用一种稀疏度自适应的正交匹配追踪算法将数据重构，并且通过在数据处理阶段进行关于稀疏度的参数分析，提高了其数据重构效果和采样数据的压缩率；利用均方根误差和相关系数作为统计指标对数据向量的特性进行识别，从而识别结构损伤。
如图6～9所示，具体地，该方法包括如下步骤：
(1)基于压电导纳技术，采集结构初始状态时的导纳数据，将原始电导信号作最大最小归一化处理，作为基于正交匹配追踪算法的数据压缩重构子模块的输入；
优选地，在所述的最大最小归一化处理中，根据以下公式对原始导纳数据进行变换：
Figure GDA0003118647020000101
其中，
Figure GDA0003118647020000102
为原始电导信号，
Figure GDA0003118647020000103
为
Figure GDA0003118647020000104
的最小值，
Figure GDA0003118647020000105
为
Figure GDA0003118647020000106
的最大值，X为处理后的原始电导信号，长度为m。
(2)根据经验与调试选择信号X的稀疏度范围和区间步长[Sstart：Sstep：Send]，其中，Sstart为所选稀疏度最小值，Sstep为所选稀疏度步长，Send为所选稀疏度最大值，从而得到一系列稀疏度Si，下标i表示稀疏度序号。
优选地，步骤(2)所述的根据经验与调试选择稀疏度范围和区间步长的步骤，包括如下子步骤：
从一长度不大于信号X长度的几何级数序列中，以信号X长度的十分之一为基，截取一系列几何级数作为稀疏度区间进行试算，试算方法如下：
根据试错法，根据截取出的稀疏度区间生成随机测量矩阵
Figure GDA0003118647020000107
代入
Figure GDA0003118647020000108
并计算，挑选出一个使得上式值尽量小的稀疏度以确定稀疏度的区间范围[Sstart：Sstep：Send]；其中，
Figure GDA0003118647020000109
为试算的测量矩阵，X为原始数据，Y为观测向量，|| ||2表示l2范数，Y＝ΦX，Φ为随机测量矩阵，优选高斯随机测量矩阵。试错法并不一定要求找出绝对最好的结果，通常只要找出尽可能合理的结果即可，具体试错至什么程度可以由操作者视情况或者依经验自由设定。
(3)输入当前稀疏度Si，下标i表示第i个稀疏度，设置相对应的测量数n，
n＝cst×ceil(Si/log(m/Si))
设置转换基矩阵Ψ和测量矩阵Φ，并对信号X进行基于压缩感知理论与正交匹配追踪算法的数据压缩与重构步骤，由Y＝ΦX＝ΦΨα＝Θα，得到对应的传感矩阵Θ、稀疏系数矩阵α和重构数据
Figure GDA0003118647020000111
转换基矩阵Ψ只与信号
Figure GDA0003118647020000112
的长度有关，对于任意确定长度的信号，其转换基矩阵唯一。
所谓匹配追踪算法的基本思想，即将传感矩阵Θ作为完备字典库，对观测向量Y进行分解，从传感矩阵Θ中选择一个与信号Y最匹配(内积最大)的原子(即列向量)，构建系数逼近，并求出信号残差，然后继续选择与信号残差最匹配的原子(即列向量)，再次进行残差分解，反复迭代，最终信号Y可以被传感矩阵Θ中的所有原子线性表达，这些原子也在迭代过程中被重新排列，得到重构后的传感矩阵Θ。
优选地，步骤(3)包括如下子步骤：
(3.1)基于压缩感知理论将X线性映射，包括：
(3.1.1)压缩感知理论建立在原始信号自身拥有稀疏性或者在某个变换上呈现出稀疏性的基础上。然而，现实中的时域/频域信号大都不具备稀疏性，因而需要利用变换基对原始信号进行变换；首先引入一m×m的转换基矩阵Ψ＝[ψ1，ψ2，…，ψm]，其中ψi是列向量，具有m列，将信号X(也是离散序列)通过正交基Ψ变换为基系数矩阵α：X＝Ψα；
(3.1.2)通过将X线性映射为Y：Y＝ΦX，其中Φ为测量矩阵(本发明设置为高斯随机测量矩阵)，然后得到Y＝ΦX＝ΦΨα＝Θα，其中Θ为传感矩阵；
(3.2)基于正交匹配追踪理论的稀疏度自适应正交匹配追踪算法来对Θ进行重排列，并求解α，包括如下子步骤：
(3.2.1)输入测量矩阵Φ(n×m)、转换基矩阵Ψ(m×m)、观测向量Y(n×1)和当前稀疏度Si等参数；本发明中，下标(n×m)、(m×m)、(n×1)仅用于对各矩阵、向量的行数和列数进行标记，例如，Φ(n×m)表示Φ为n×m的矩阵，并非不同于Φ的新参数。
(3.2.2)设定k＝0，初始残差R0＝Y，D0＝{ }，I0＝{1，…，m}，I0为传感矩阵Θ中列向量的序号构成的初始序号矩阵。
(3.2.3)k≥1时，Y的第k步分解如下：
Figure GDA0003118647020000121
其中，p＝0，1，…，k-1，第k阶模型为
Figure GDA0003118647020000122
Figure GDA0003118647020000123
为分解系数；Rk为第k步的当前残差，满足：
Figure GDA0003118647020000124
Ik-1＝Ik-2\{rk-2}
Figure GDA0003118647020000125
即，第k步分解时，遍历传感矩阵Θ中序号属于集合Ik-1的列向量，从中找出一个列向量dj，使得当前残差Rk和该列向量dj能够取得最大内积，将该列向量dj作为
Figure GDA0003118647020000126
加入到矩阵Dk中，然后将
Figure GDA0003118647020000127
在传感矩阵Θ中对应的列向量的序号从集合Ik-1中剔除。
以k＝1和k＝2为例：
第k＝1步分解的目标是寻找出一个列向量
Figure GDA0003118647020000128
使当前残差R1和该列向量
Figure GDA0003118647020000129
能够取得最大内积，满足
Figure GDA00031186470200001210
其中sup表示上确界，令I0＝{1，…，m}，并且更新选择列向量集合
Figure GDA00031186470200001211
因此第k＝1步分解得到
Figure GDA00031186470200001212
满足
Figure GDA00031186470200001213
Figure GDA00031186470200001214
且
Figure GDA00031186470200001215
第k＝2步，对R1再次进行分解得到
Figure GDA00031186470200001216
满足
Figure GDA00031186470200001217
且
Figure GDA00031186470200001218
以此类推即可得到式(2)。
(3.2.4)令k＝k+1，重复步骤(3.2.3)，直至k＝m时，传感矩阵Θ中的所有列向量遍历完毕，终止迭代，得到Dk。
(3.2.5)根据公式(2)，在进行Y的第k步分解时，取k阶残差迭代模型为：
Figure GDA0003118647020000131
则k+1阶残差迭代模型减去k阶残差迭代模型得：
Figure GDA0003118647020000132
引入辅助模型
Figure GDA0003118647020000133
Figure GDA0003118647020000134
其中，
Figure GDA0003118647020000135
为第k-1阶分解系数，γk为
Figure GDA0003118647020000136
的残差项，q＝0，1，…，k-1；
由于
Figure GDA0003118647020000137
其中，t＝0，1，…，k-1；
为了计算
Figure GDA0003118647020000138
将式(5)重写为一个线性系统
uk＝Ckβk (7)
Figure GDA0003118647020000139
Figure GDA00031186470200001310
Figure GDA00031186470200001311
由(3.2.4)可知，Ck满足非奇异性，因而式(7)可改写为
Figure GDA00031186470200001312
其中，
Figure GDA00031186470200001313
为Ck的逆矩阵，
由于Ck+1可以改写为
Figure GDA00031186470200001314
其中*表示共轭转置，因而可以得到
Figure GDA0003118647020000141
其中，
Figure GDA0003118647020000142
因而可以通过递归计算
Figure GDA0003118647020000143
以得到βk+1；
(3.2.6)将式(4)代入(3)，
Figure GDA0003118647020000144
假如式(15)和式(16)成立，则式(14)成立，可得：
Figure GDA0003118647020000145
Figure GDA0003118647020000146
令
Figure GDA0003118647020000147
由(7)左右两边对
Figure GDA0003118647020000148
作内积，因为Rk+1与前面所有每个分量即
Figure GDA0003118647020000149
均正交，因而可得
Figure GDA00031186470200001410
其中，|| ||2表示向量作差之后各分量的平方和。
Figure GDA00031186470200001411
(3.2.7)更新模型k＝k+1：
Figure GDA00031186470200001412
Rk+1＝Y-Yk+1 (20)
Dk+1＝Dk∪{dk} (21)
Ik＝Ik-1\{rk-1} (22)
重复步骤(3.2.5)～(3.2.6)，直至k＝m时，求解出所有的
Figure GDA00031186470200001413
按照每次迭代时
Figure GDA00031186470200001414
对应的
Figure GDA00031186470200001415
在Dk中的顺序组成稀疏系数矩阵α，即得到原始信号的压缩数据；
(3.2.8)输出当前稀疏度Si对应的稀疏系数矩阵α和重构数据
Figure GDA00031186470200001416
(4)根据式(22)和式(23)计算当前稀疏度下重构数据所得到的总误差TE与压缩率CR；
将(3.2.8)中的重构数据代入下式计算压缩率CR和总误差TE；
CR＝Sorig/scomp (23)
Figure GDA0003118647020000151
式中，Sorig表示初始数据
Figure GDA0003118647020000152
所占字节，Scomp表示压缩数据(即稀疏系数矩阵α)所占字节，
Figure GDA0003118647020000153
表示重构数据；
(5)回到步骤(3)，输入稀疏度Si+1，重复步骤(3)和(4)，直到遍历稀疏度区间范围内的所有稀疏度；Si+1＝Si+Sstep；
(6)根据最终的稀疏度与总误差的关系、稀疏度与压缩率的关系进行参数分析，具体地：通过稀疏度和总误差的对应关系寻找平衡，即在保证总误差低于期望值的同时，尽量使得压缩率最高，从而得出每一次电导数据处理时的最优稀疏度，存储当前最优稀疏度对应的重构数据，准备进入损伤识别模块。优选地，根据稀疏度与总误差、稀疏度与压缩率的对应关系曲线，寻找总误差最低的3～4个值所对应的稀疏度，再比较上述稀疏度对应的压缩率大小，压缩率不小于2，从而筛选出最优稀疏度。
(7)损伤识别：
通过将最优稀疏度对应的重构数据与原始数据进行统计指标计算对比，即利用均方根偏差和相关系数计算，分析比较并得出损伤诊断结果，包括如下子步骤：
(7.1)以未损工况下的电导数据为基线，以均方根偏差和相关系数为损伤因子计算各工况指标RMSD(单位：％)：
Figure GDA0003118647020000154
CC＝Cov(X0，X1)/σ0σ1 (26)
式中，
Figure GDA0003118647020000161
表示第i个采样点下的基线导纳信号值，
Figure GDA0003118647020000162
表示相应采样点的待测电导信号值，σ0和σ1分别表示基线电导信号和待测电导信号的标准差，Cov(*)为协方差运算符，
Figure GDA0003118647020000163
优选地，还可以进行如下验证步骤：
(7.2)将步骤(7.1)计算获得的RMSD和CC与根据原始数据
Figure GDA0003118647020000164
计算所得到的损伤指标进行对比，验证该方法的有效性。
下面，结合一个具体应用实例对本发明的效果进行介绍。
在本实例中，试件是一根500mm×35mm×5mm(长宽高)的矩形截面Q235钢梁，如图1所示，简支梁两端用橡胶块固定在桌子上，整个小型健康监测系统共设置6个PZT传感器，传感器布置情况如图1所示。如图1，损伤识别试验共设置4种工况，PZT#4作为参考传感器，与梁粘结情况始终为部分粘结(具体系用高强环氧树脂仅粘贴PZT片面积的3/4)。第一种工况中，梁未损且传感器状况良好；第二种工况中，梁保持原样不变，但将PZT#5用刀片划伤，将PZT#6损坏；第三种工况中，在钢梁跨中放置8枚直径19mm、厚度1.67mm、质量为3.20g的硬币，以增加集中质量；第四种工况中，移去梁跨中硬币，用短锯将钢梁割出一道凹槽。整个试验过程在实验室环境下进行，温度大致20度，利用笔记本和4294A阻抗分析仪(1V激励)完成所有导纳数据的测量。
当稀疏度取为1300时，原始电导信号X与重构信号Y的对比如图2(a)～图2(i)所示，工况1下PZT#1、PZT#2、PZT#3、PZT#5和PZT#6传感器可视为基线数据，而PZT#4在整个过程中均作为一个针对传感器粘结问题的参考对象。由图可知，虽然在曲线在上升段中有一部分出现了一定的振荡现象，但总体而言，重构数据基本与相应的原始信号一致。每一个工况，所有6个传感器数据的压缩和重构过程耗时仅5260s，对于实施结构的实时监测来说是可以接受的。图中结果表明不同种类的工况(结构和传感器)不会影响到本发明压缩重构算法的可行性。
通过统计上述不同工况下稀疏度与压缩率的对应关系可知，当稀疏度为400时，压缩率可达到最高约10左右，而当稀疏度取为2100时，压缩率低至1.7上下。压缩率曲线呈现非线性，且随着稀疏度的增加而单调递减。
图3(a)～图3(i)依次是PZT#2(工况1)、PZT#4(工况1)、PZT#5(工况2)、PZT#6(工况2)、PZT#2(工况3)、PZT#5(工况3)、PZT#6(工况3)、PZT#5(工况4)、PZT#6(工况4)不同稀疏度下的重构数据对应误差。图3(a)～图3(i)表明，随着稀疏度的增加，总误差曲线先是大幅度下降，然后到达一个平稳阶段，最后又缓慢上升。大部分的最优误差在稀疏度为1400-1600之间时取得，且最优误差取值大部分在1.67％-2.79％范围内，但是PZT#4的最优误差相对较大，达到了6.84％(工况1)、7.24％(工况3)、6.69％(工况4)，原因可能系传感器有粘结缺陷导致了重构数据的部分畸变。
由图4(a)可知，PZT#2原始数据工况3与基线工况1的RMSD指标为0.071，而对应的重构数据为0.076(识别误差为7.183％)。而PZT#3重构的效果更佳，原始数据和重构数据的RMSD指标相差误差仅为1.499％。除此之外，由图4(a)、4(b)可知，对于梁跨中增加质量块的工况与凹槽损伤的工况，能够清晰的辨别两者RMSD指标的不同，从前者到后者RMSD值呈现一个大幅度增长。图4(c)描述了PZT#5和PZT#6的RMSD指标情况，梁跨中受损工况对比质量增加工况均有上涨，其中工况4中PZT#6的RMSD值为0.442，8倍于其在工况3的RMSD值，而PZT#5的值相对于工况3仅增加32.899％，表明从增加质量块到结构受损虽然属于损伤进阶，但由于PZT#6传感器损坏较为严重，而导致部分数据失真、RMSD相差过大。并且工况3中PZT#6数据重构误差达到17.082％，表明该传感器受损且不能提供有效用的数据反馈，相比较而言，仅有划痕的传感器PZT#5的重构误差均在6％以内，表明在传感器仅受到轻微损伤时数据依然有效。
图5(a)表示PZT#2和#3的重构数据得到的CC指标与原始数据吻合较好，并且由图5(c)可以得到PZT#5和PZT#6有着类似的结果。如图5(a)～图5(c)所示，由于CC指标与结构损伤程度存在一定的负相关，PZT#2、PZT#3、PZT#4导纳信号得到的CC值，从工况3到工况4均呈现一定的降低趋势，表明凹槽带来的损伤程度大于质量增加的情形。另外，工况2中PZT#5(0.984)和PZT#6(0.734)的CC值的下降可归结于传感器的划痕和损伤，且PZT#5从工况3的0.995至工况4的0.989、以及PZT#6从工况3的0.997至工况4的0.772，均表明了结构损伤程度的加重。尽管表现方式相反，相关系数所呈现的识别结果与均方根偏差是一致的。除了PZT#4重构与原始数据的识别误差达到了0.964％外(原因应该系传感器部分连接导致重构的导纳信号存在过多震荡峰值)，其它误差均在0.2％以内，表明了本发明的重构数据利用相关系数作为损伤因子识别的高效性。
本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (9)

1.一种基于压电导纳的数据压缩重构方法，其特征在于，包括如下步骤：
(1)基于压电导纳技术，采集结构初始状态时的导纳数据，将原始电导信号作如下最大最小归一化处理：
Figure FDA0003146912780000011
其中，
Figure FDA0003146912780000012
为原始电导信号，
Figure FDA0003146912780000013
为
Figure FDA0003146912780000014
的最小值，
Figure FDA0003146912780000015
为
Figure FDA0003146912780000016
的最大值，X为处理后的原始电导信号，长度为m；
(2)设定稀疏度区间范围及步长[Sstart:Sstep:Send]，其中，Sstart为所选稀疏度最小值，Sstep为所选稀疏度步长，Send为所选稀疏度最大值，从而得到一系列稀疏度Si，下标i表示稀疏度序号；设定稀疏度区间范围及步长[Sstart:Sstep:Send]的方法如下：
从一长度不大于信号X长度的几何级数序列中，截取一系列几何级数作为稀疏度区间进行试算，试算方法如下：
根据试错法，根据截取出的稀疏度区间生成随机测量矩阵
Figure FDA0003146912780000017
代入
Figure FDA0003146912780000018
并计算，挑选出一个使得上式值尽量小的稀疏度以确定稀疏度的区间范围[Sstart:Sstep:Send]；
其中，
Figure FDA0003146912780000019
为试算时设定的随机测量矩阵，Y为观测向量，‖‖2表示l2范数，Y＝ΦX，Φ为不同于
Figure FDA00031469127800000110
的随机测量矩阵；
(3)输入当前稀疏度Si，下标i表示第i个稀疏度，设置相对应的测量数n，
Figure FDA00031469127800000111
设置m×m的转换基Ψ和n×m的测量矩阵Φ，m×m和n×m表示矩阵行列数，基于公式Y＝ΦX＝ΦΨα＝Θα，对信号X进行基于压缩感知理论与正交匹配追踪算法的数据压缩与重构，包括：
将传感矩阵Θ作为完备字典库，对观测向量Y进行分解，即从传感矩阵Θ中选择一个与信号Y最匹配的列向量，构建稀疏系数逼近Y，并求出信号残差，然后继续选择与信号残差最匹配的列向量，再次进行残差分解，反复迭代，最终信号Y被传感矩阵Θ中的所有列向量表达，这些列向量也在迭代过程中被重新排列，按照每次迭代获得的稀疏系数的顺序，将所有稀疏系数排列组成稀疏系数矩阵α，稀疏系数矩阵α即为压缩数据，对应的重构数据为
Figure FDA0003146912780000021
(4)计算当前稀疏度Si对应的重构数据
Figure FDA0003146912780000022
与原始数据X比对下得到的总误差与压缩率；
(5)回到步骤(3)，输入稀疏度Si+1，重复步骤(3)和(4)，直到遍历稀疏度区间范围内的所有稀疏度；Si+1＝Si+Sstep；
(6)根据步骤(5)获得的所有稀疏度对应的总误差及压缩率，绘制稀疏度与总误差的关系曲线、稀疏度与压缩率的关系曲线，基于上述关系曲线，找出保证总误差低于期望值的同时尽量使得压缩率最高的稀疏度，作为最优稀疏度，并得到当前最优稀疏度对应的稀疏系数矩阵和重构数据，完成数据压缩与重构。
2.如权利要求1所述的一种基于压电导纳的数据压缩重构方法，其特征在于，步骤(2)中，从一长度不大于信号X长度的几何级数序列中，以信号X长度的十分之一为基，截取一系列几何级数作为稀疏度区间进行试算。
3.如权利要求2所述的一种基于压电导纳的数据压缩重构方法，其特征在于，步骤(3)包括如下子步骤：
(3.1)基于压缩感知理论将X进行线性映射，包括：
引入一个m×m的转换基矩阵Ψ＝[ψ1,ψ2,…,ψm]，其中ψ1～ψm均为列向量，具有m列，将信号(离散序列)X通过正交基Ψ变换为基系数矩阵α，有：X＝Ψα；
又由于Y＝ΦX，因此有Y＝ΦX＝ΦΨα＝Θα，其中Θ为传感矩阵；
(3.2)基于正交匹配追踪理论的稀疏度自适应正交匹配追踪算法，对步骤(2)所得到的稀疏度区间[Sstart:Sstep:Send]中每一个稀疏度，遍历传感矩阵Θ，从而求解每个稀疏度对应的α以及每个稀疏度对应的重构数据。
4.如权利要求3所述的一种基于压电导纳的数据压缩重构方法，其特征在于，步骤(3.2)包括如下子步骤：
(3.2.1)输入测量矩阵Φ(n×m)、转换基矩阵Ψ(m×m)、观测向量Y(n×1)和当前稀疏度Si；
(3.2.2)设定k＝0，初始残差R0＝Y，D0＝{}，I0＝{1,…,m}，I0为传感矩阵Θ中列向量的序号构成的初始序号矩阵；
(3.2.3)k≥1时，Y的第k步分解如下：
Figure FDA0003146912780000031
其中，p＝0,1,…,k-1，第k阶模型为
Figure FDA0003146912780000032
Figure FDA0003146912780000033
为分解系数；Rk为第k步的当前残差，满足：
Figure FDA0003146912780000034
Ik-1＝Ik-2\{rk-2}
Figure FDA0003146912780000035
即，第k步分解时，遍历传感矩阵Θ中序号属于集合Ik-1的列向量，从中找出一个列向量dj，使得当前残差Rk和该列向量dj能够取得最大内积，将该列向量dj作为
Figure FDA0003146912780000036
加入到矩阵Dk中，然后将
Figure FDA0003146912780000037
在传感矩阵Θ中对应的列向量的序号从集合Ik-1中剔除；
(3.2.4)令k＝k+1，重复步骤(3.2.3)，直至k＝m时，传感矩阵Θ中的所有列向量遍历完毕，终止迭代，得到Dk，Dk为传感矩阵Θ重排列后的结果；
(3.2.5)根据公式(2)，在进行Y的第k步分解时，取k阶残差迭代模型为：
Figure FDA0003146912780000041
则k+1阶残差迭代模型减去k阶残差迭代模型得：
Figure FDA0003146912780000042
引入辅助模型
Figure FDA0003146912780000043
Figure FDA0003146912780000044
Figure FDA0003146912780000045
其中，
Figure FDA0003146912780000046
为第k-1阶分解系数，q＝0,1,…,k-1，γk为残差项；
由于
Figure FDA0003146912780000047
为了计算
Figure FDA0003146912780000048
将式(5)重写为一个线性系统：
uk＝Ckβk (7)
其中：
Figure FDA0003146912780000049
Figure FDA00031469127800000410
Figure FDA00031469127800000411
Ck满足非奇异性，因而式(7)可改写为
Figure FDA00031469127800000412
其中，
Figure FDA00031469127800000413
为Ck的逆矩阵，
由于Ck+1可以改写为：
Figure FDA0003146912780000051
其中*表示共轭转置，因而可以得到：
Figure FDA0003146912780000052
其中
Figure FDA0003146912780000053
因而可以通过递归计算
Figure FDA0003146912780000054
得到βk+1；
(3.2.6)将式(4)代入式(3)，得到：
Figure FDA0003146912780000055
假如式(15)和式(16)成立，则式(14)成立，可得：
Figure FDA0003146912780000056
Figure FDA0003146912780000057
令
Figure FDA0003146912780000058
由(7)左右两边对
Figure FDA0003146912780000059
作内积，得：
Figure FDA00031469127800000510
其中，‖‖2表示向量作差之后各分量的平方和；
Figure FDA00031469127800000511
(3.2.7)令k＝k+1，重复步骤(3.2.5)～(3.2.6)，直至k＝m时，求解出所有的
Figure FDA00031469127800000512
将所有
Figure FDA00031469127800000513
按照对应的
Figure FDA00031469127800000514
在Dk中的顺序，组成稀疏系数矩阵α，α即为
Figure FDA00031469127800000515
和X的压缩数据；
(3.2.8)输出当前稀疏度Si对应的稀疏系数矩阵α和重构数据
Figure FDA00031469127800000516
5.如权利要求1～4任意一项所述的一种基于压电导纳的数据压缩重构方法，其特征在于，步骤(4)中，总误差TE与压缩率CR计算如下：
CR＝Sorig/Scomp
Figure FDA0003146912780000061
其中，Sorig表示初始数据
Figure FDA0003146912780000062
所占字节，Scomp表示压缩数据α所占字节。
6.如权利要求5所述的一种基于压电导纳的数据压缩重构方法，其特征在于，步骤(6)中，根据稀疏度与总误差、稀疏度与压缩率的对应关系曲线，寻找总误差最低的3～4个值所对应的稀疏度，再比较上述稀疏度对应的压缩率大小，压缩率不小于2，从而筛选出最优稀疏度。
7.一种基于压电导纳的数据压缩重构设备，其特征在于，包括压电导纳信号采集模块、处理器和数据压缩重构程序模块；
所述数据压缩重构程序模块在被所述处理器调用时，按照如权利要求1～6任意一项所述的数据压缩重构方法对压电导纳信号采集模块采集的电导信号，进行压缩重构处理。
8.一种结构损伤识别方法，其特征在于，在待测结构上布置多个采样点，基于压电导纳技术采集各个采样点位置的电导纳信号，并对其中的电导信号部分，按照权利要求1～6任意一项所述的数据压缩重构方法筛选出最优稀疏度，并将最优稀疏度所对应的重构数据
Figure FDA0003146912780000063
与基线电导数据进行均方根偏差和相关系数计算，得出工况指标RMSD和CC：
Figure FDA0003146912780000064
Figure FDA0003146912780000065
式中，
Figure FDA0003146912780000066
表示第i个采样点下的基线电导信号值，
Figure FDA0003146912780000067
表示相应采样点的待测电导信号值，σ0和σ1分别表示基线电导信号和待测电导信号的标准差，Cov(*)为协方差运算符，
Figure FDA0003146912780000068
RMSD指标越大，则结构发生损伤的可能性越大，CC指标越接近1，则结构健康状况越好。
9.一种结构损伤识别设备，其特征在于，包括压电导纳信号采集模块、处理器、数据压缩重构程序模块以及结构损伤识别程序模块；
所述数据压缩重构程序模块在被所述处理器调用时，按照如权利要求1～6任意一项所述的数据压缩重构方法对压电导纳信号采集模块采集的电导信号，进行压缩重构处理，得到最优稀疏度及其对应的重构数据
Figure FDA0003146912780000071
所述结构损伤识别程序模块在被所述处理器调用时，按照如权利要求8所述的结构损伤识别方法，计算出工况指标RMSD和CC。
CN201911375347.3A 2019-12-27 2019-12-27 基于压电导纳的数据压缩重构及结构损伤识别方法、设备 Active CN111010193B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911375347.3A CN111010193B (zh) 2019-12-27 2019-12-27 基于压电导纳的数据压缩重构及结构损伤识别方法、设备

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911375347.3A CN111010193B (zh) 2019-12-27 2019-12-27 基于压电导纳的数据压缩重构及结构损伤识别方法、设备

Publications (2)

Publication Number Publication Date
CN111010193A CN111010193A (zh) 2020-04-14
CN111010193B true CN111010193B (zh) 2021-08-20

Family

ID=70119188

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911375347.3A Active CN111010193B (zh) 2019-12-27 2019-12-27 基于压电导纳的数据压缩重构及结构损伤识别方法、设备

Country Status (1)

Country Link
CN (1) CN111010193B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112378958A (zh) * 2020-10-23 2021-02-19 华中科技大学 一种结构损伤识别方法及系统

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103944579A (zh) * 2014-04-10 2014-07-23 东华大学 一种压缩感知重构的编码解码系统
CN108832934A (zh) * 2018-05-31 2018-11-16 安徽大学 一种基于奇异值分解的二维正交匹配追踪优化算法
CN110061744A (zh) * 2019-04-11 2019-07-26 池州学院 一种基于压缩感知理论的电能质量信号自适应重构方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103944579A (zh) * 2014-04-10 2014-07-23 东华大学 一种压缩感知重构的编码解码系统
CN108832934A (zh) * 2018-05-31 2018-11-16 安徽大学 一种基于奇异值分解的二维正交匹配追踪优化算法
CN110061744A (zh) * 2019-04-11 2019-07-26 池州学院 一种基于压缩感知理论的电能质量信号自适应重构方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
基于压电传感机械阻抗的结构损伤识别方法研究;艾德米;《中国博士学位论文全文数据库 工程科技Ⅱ辑》;20190415(第4期);全文 *

Also Published As

Publication number Publication date
CN111010193A (zh) 2020-04-14

Similar Documents

Publication Publication Date Title
CN107666322B (zh) 一种基于字典学习的自适应微震数据压缩感知方法
CN109212018B (zh) 一种基于压缩感知的管道漏磁数据分级压缩方法
CN106841402A (zh) 一种基于贪婪算法的相控阵超声信号重构优化方法
CN106443183B (zh) 一种从电磁频谱中提取频率源数量及基频值的方法
CN111010193B (zh) 基于压电导纳的数据压缩重构及结构损伤识别方法、设备
Torkamani et al. Detection of system changes due to damage using a tuned hyperchaotic probe
CN109815940B (zh) 小波包能量谱法损伤识别方法
CN111611654A (zh) 一种铆接结构的疲劳预测方法、装置及设备、存储介质
CN109738392A (zh) 面向tdlas在线检测瓶内氧气浓度的压缩感知重构方法及系统
CN116150594B (zh) 一种频谱测试数据中开关类要素特征的识别方法
Chen et al. Application of EMD-AR and MTS for hydraulic pump fault diagnosis
CN117434337A (zh) 一种宽频电压信号测量方法及系统
CN112183297B (zh) 一种超声相控阵信号稀疏特征提取方法
CN111678991A (zh) 一种用于混凝土结构无损检测损伤识别的方法
CN116522080A (zh) 局部放电信号降噪方法
CN114186624A (zh) 基于二维卷积神经网络的损伤识别模型建立方法及设备
CN110261493A (zh) 一种基于滤波信号压缩感知的航空材料破损检测方法
CN113052018B (zh) 一种钣金拉伸过程声发射信号识别方法
CN117743752A (zh) 基于ceemd和压缩感知的大坝安全监测数据降噪方法及系统
CN113592010B (zh) 基于局部omp的快速卷积稀疏字典学习方法
CN114626248B (zh) 一种基于多应力加速退化数据的螺旋弹簧可靠性评估方法
US10840942B2 (en) Method and apparatus for compressing data
CN116818574A (zh) 基于太赫兹谱特征参数表征的材料力学性能测定方法
Lundin et al. On the estimation of quantizer reconstruction levels
CN115982548A (zh) 一种用于线谱估计的多测量向量深度原子范数去噪网络及其构建方法

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