CN108053379B - 一种基于改进的变分模态分解的dspi相位提取方法 - Google Patents

一种基于改进的变分模态分解的dspi相位提取方法 Download PDF

Info

Publication number
CN108053379B
CN108053379B CN201711330483.1A CN201711330483A CN108053379B CN 108053379 B CN108053379 B CN 108053379B CN 201711330483 A CN201711330483 A CN 201711330483A CN 108053379 B CN108053379 B CN 108053379B
Authority
CN
China
Prior art keywords
modal
component
decomposition
phase
dspi
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
CN201711330483.1A
Other languages
English (en)
Other versions
CN108053379A (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.)
Tianjin University
Original Assignee
Tianjin University
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 Tianjin University filed Critical Tianjin University
Priority to CN201711330483.1A priority Critical patent/CN108053379B/zh
Publication of CN108053379A publication Critical patent/CN108053379A/zh
Application granted granted Critical
Publication of CN108053379B publication Critical patent/CN108053379B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00Image enhancement or restoration
    • G06T5/70Denoising; Smoothing
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00Image enhancement or restoration
    • G06T5/40Image enhancement or restoration using histogram techniques

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Length Measuring Devices By Optical Means (AREA)
  • Image Processing (AREA)

Abstract

本发明公开了一种基于改进的变分模态分解的DSPI相位提取方法,包括:采集DSPI相位图,根据相位图长度和宽度估算分解后的模态分量的个数,利用正交指标提取最佳模态数量;根据最佳模态数量对DSPI相位图进行变分模态分解,获得一系列的模态函数分量;计算模态函数分量的直方图,利用直方图分析分解后的模态分量是否包含噪声,对含噪分量进行IVMD分解,分别计算二次分解后分量的直方图,根据直方图获得二次分解后的无噪声分量;将二次分解后的无噪声分量和一次分解后的无噪声分量进行重构,获得降噪后散斑图;采用希尔伯特变化法计算变形前后的相位图,对其进行相减得到DSPI相位图,对相位进行解包裹及位移提取。本方法降低噪声干扰,获取精确相位信息。

Description

一种基于改进的变分模态分解的DSPI相位提取方法
技术领域
本发明涉及光学图像处理技术领域,尤其涉及一种基于改进的变分模态分解(Improved variational mode decomposition,IVMD)的DSPI相位提取方法。
背景技术
近年来,随着航空航天的发展,各种复合材料获得广泛应用。然而,这些复合材料在加工制造和服役过程中由于变形、位移等缺陷导致构件性能显著下降,减少服役时间。为此,需要对这些复合材料进行检测及评估,为后续抢修提供方案,保障航空航天设备的健康运行。
研究表明当被相干光照射的物体表面发生形变或者位移时,物面的形变就转化成像面上散斑的相位变化,这种现象称为散斑干涉。数字散斑干涉(Digital SpecklePattern Interferometry,DSPI)是一种非接触式全场测量技术,它可以对复合材料的位移、形变、表面缺陷等进行测量。目前国内外散斑相位提取技术主要分为两类:一类是基于单幅散斑图的相位提取技术,即空间载波法,一类是以相移法为代表的基于多幅散斑图的相位提取技术。与采用多幅散斑图的相移法相比,空间载波技术是在某一时刻只采集一幅图像,受环境扰动的影响较小,更适合动态测量。采用合适的信号处理方法对散斑图片进行处理,对提高空间载波相位提取法的测量精度具有重要意义。
发明人在实现本发明的过程中,发现现有技术中至少存在以下缺点和不足:
1、传统的傅立叶变换和小波变换不具有自适应性,不能有效的获得准确相位信息;
2、由于经验模态分解过程中完全基于信号本身特性,无需人为选择基函数,因此获得广泛应用,但该方法也存在一些缺点,如:模态混叠、缺乏理论支撑等,不能有效的提取相位信息。
发明内容
本发明提供了一种基于改进的变分模态分解的DSPI相位提取方法,本发明针对DSPI测量中传统相位提取方法精度低,采用IVMD对变形前后的散斑图进行处理,无需参数设置,可以对图片进行分解,降低噪声干扰,获取精确相位信息,详见下文描述:
一种基于改进的变分模态分解的DSPI相位提取方法,所述方法包括以下步骤:
1)采集DSPI相位图,根据相位图长度和宽度估算分解后的模态分量的个数,利用正交指标提取最佳模态数量;
2)根据最佳模态数量对DSPI相位图进行变分模态分解,获得一系列的模态函数分量;
3)计算模态函数分量的直方图,利用直方图分析分解后的模态分量是否包含噪声,对含噪分量继续进行IVMD分解,分别计算二次分解后分量的直方图,根据直方图获得二次分解后的无噪声分量;
4)将二次分解后的无噪声分量和一次分解后的无噪声分量进行重构,获得降噪后的散斑图;
5)采用希尔伯特变化法计算变形前后的相位图,对其进行相减得到DSPI相位图。
6)对相位图进行解包裹并计算位移值。
上述步骤3)具体为:
计算分解后的模态分量直方图,若分解后的模态分量直方图符合正态则为含噪模态分量,反之则为无噪声分量;
对含噪模态分量继续进行IVMD分解,获得一系列分量;
计算二次分解后的直方图,根据直方图的分布提取二次分解后的无噪声分量。
上述步骤5)具体为:
采用希尔伯特法计算降噪后散斑图的相位;
利用变形前后的散斑相位进行相减,相减结果为DSPI相位图。
上述步骤1)具体为:
判断相位图的长度和宽度,两者相等的时候,利用长度或宽度可以计算出来模态分量的个数;
长度和宽度不相等的时候,首先根据尺度估算模态数量的范围,分别计算分解后的正交性,选取正交最小时对应的分量个数,即为最佳模态数量。
上述步骤2)具体为:
将模态分量定义为具有不同中心频率的有限带宽;根据不同模态的带宽之和最小原则构造约束变分方程;
针对约束变分方程,引入增广Lagrange函数将约束变分方程转化为非约束变分方程;
采用乘法算子交替方向法计算增广Lagrange函数的鞍点,即约束变分方程的最优解;
不断更新中心频率及模态分量,当模态分量满足迭代停止条件时,停止更新,输出分解后的一系列的模态函数分量。
本发明提供的技术方案的有益效果是:
1、针对变分模态分解(variational mode decomposition,VMD)分解过程中需要设置模态数量特点,本发明提出改进的变分模态分解法,它可以在变分框架内选取最佳模态数量,对变形前后散斑图进行处理,获得固有模态分量;
2、针对采集的散斑图中包含的噪声分量,本发明利用IVMD对图片进行多次分解,根据模态分量的直方图是否符合正态分布滤除噪声分量,获取滤波后的散斑图,提高信噪比,进而提高相位精确信息;
3、本发明提出基于IVMD的相位提取方法能够自适应的计算DSPI相位图,避免传统相位法繁琐的参数设置。
附图说明
图1为一种基于改进的变分模态分解的DSPI相位提取方法的流程图;
图2为数字散斑干涉测量系统的结构示意图;
图3为被测圆盘的示意图;
图4为圆盘变形前后采集的散斑图的示意图;
图5为变形相位图。
具体实施方式
为使本发明的目的、技术方案和优点更加清楚,下面对本发明实施方式作进一步地详细描述。
实施例1
一种基于改进的变分模态分解的DSPI相位提取方法,参见图1,该方法包括以下步骤:
101:采集DSPI相位图,根据相位图的长度和宽度估算分解后的模态分量的个数,利用正交指标提取最佳模态数量,根据分解后的模态分量的个数对DSPI相位图进行变分模态分解,获得一系列的模态函数分量;
102:计算模态函数分量的第一直方图,根据第一直方图获得含噪模态分量和无噪声分量,对每个含噪模态分量分别进行IVMD分解,分别计算分解后分量的第二直方图,根据第二直方图获得二次分解后的无噪声分量;
103:将二次分解后的无噪声分量和一次分解后的无噪声分量进行重构,得到降噪后的散斑图;
104:利用希尔伯特法计算降噪后散斑图的相位,对散斑变形前后的相位相减获得DSPI相位;
105:采用逐行逐列的方式对DSPI相位进行解包裹并提取位移。
其中,上述步骤101具体为:
判断相位图的长度和宽度,两者相等的时候,利用长度或宽度可以计算出来模态分量的个数;
长度和宽度不相等的时候,首先根据尺度估算模态数量的范围,分别计算分解后的正交性,选取正交最小时对应的分量个数,即为最佳模态数量。
首先将模态分量定义为具有不同中心频率的有限带宽;根据不同模态的带宽之和最小原则构造约束变分方程;
针对约束变分方程,引入增广Lagrange函数将约束变分方程转化为非约束变分方程;
采用乘法算子交替方向法计算增广Lagrange函数的鞍点,即约束变分方程的最优解;
不断更新中心频率及模态分量,当模态分量满足迭代停止条件时,停止更新,输出分解后的一系列的模态函数分量。
进一步地,上述步骤102具体为:
计算散斑图经IVMD分解后模态分量的直方图;
分析模态分量的直方图是否符合正态分布,若符合正态分布则模态分量为含噪分量,反之为无噪声分量;
利用IVMD对含噪分量继续处理,分析分解后的直方图,提取二次分解后的无噪声分量。
进一步地,上述步骤104具体为:
采用希尔伯特法分别计算降噪后散斑图的相位图;
对变形前后的相位图进行相减,获得DSPI相位图。
综上所述,本发明实施例通过上述步骤101-步骤105实现了在包含噪声的散斑图前提下,提出改进的变分模态分解对散斑图进行处理,根据直方图降低噪声干扰,采用希尔伯特法计算降噪后的相位图,提高了测量精度。
实施例2
下面结合具体的计算公式、实例、图1-图5对实施例1中的方案进行进一步地介绍,详见下文描述:
201:结合图2搭建数字散斑干涉测量系统,利用数字散斑干涉测量系统内的CCD相机采集被测圆盘变形前后的散斑图,该步骤的详细操作为:
搭建数字散斑干涉测量系统,该测量系统由CCD相机、成像透镜、激光器等组成。
其中,该测量系统的光路如图2所示,激光器出射的激光经过分光镜分成两束光,一束光照射被测物表面,另一束光经过耦合透镜沿着光纤传输作为物光,被测物的漫反射光依次经过光澜、成像透镜,与物光形成散斑干涉,通过CCD相机采集散斑图,被测圆盘(参见图3)面板的材料是铜片,采集的散斑图如图4所示。
其中,本发明实施例对铜片的尺寸不做限制,根据实际应用中的需要进行设定。
202:根据数字散斑干涉相位图的尺寸设定分解后的模态分量个数,该步骤的详细操作为:
1)对任意一维信号进行分解,其分解后的模态分量个数为:
D=log2L-1
其中,L为一维信号的长度,D为分解后的模态分量个数。
2)假设DSPI相位图的尺寸为M×N,经VMD分解后的分量个数为K;
Figure BDA0001506537740000051
利用正交性指标提取最佳模态数量,获取合适的模态分量。
VMD方法是在变分框架内通过迭代搜寻变分模型的最优解来确定分解后分量的频率中心以及带宽,从而能够自适应的对数字散斑相位图进行分解。
3)初始化
Figure BDA0001506537740000052
和n←0,对变分问题进行构造,对应的约束变分方程如下:
Figure BDA0001506537740000053
Figure BDA0001506537740000054
其中,f(x)为DSPI相位图,uk(x)为2D信号分解后的分量即本征模态函数,利用Hilbert变换得到单边频谱的2D解析信号uAS,k(x),它的数学表达式如下:
Figure BDA0001506537740000055
其中,ωk为中心频率;αk为惩罚参数;X为图片的矢量;k为分解后的数量;uk为分解后的分量;δ(<x,ωk>为狄拉克函数;δ(<x,ωk,⊥>)为ωk频带下的反傅里叶变换;⊥为反傅里变换;π<x,ωk>为参数。
4)针对约束变分问题,引入增广Lagrange函数将约束变分问题转化为非约束变分问题,增广Lagrange函数的数学表达式如下:
Figure BDA0001506537740000061
式中惩罚参数为αk;Lagrange函数乘子为λ;λ(x)为乘子函数;▽为计算范数;<.>为卷积。
5)为了解决最优解这一问题,采用乘法算子交替方向法计算增广Lagrange函数的鞍点,即约束变分方程的最优解。交替更新获得模态分量和中心频率数学表达式如下:
Figure BDA0001506537740000062
Figure BDA0001506537740000063
其中,i为参数,取值范围为1到k。
6)更新Lagrange函数乘子λ。
Figure BDA0001506537740000064
其中,
Figure BDA0001506537740000065
为乘子的频域函数;τ为系数;
Figure BDA0001506537740000066
为f(x)的频域函数,
Figure BDA0001506537740000067
Figure BDA0001506537740000068
的频域函数。
7)如果
Figure BDA0001506537740000069
结束循环,输出模态分量,否则继续循环。
其中,
Figure BDA00015065377400000610
为分解后的模态分量,ε为无穷小的正数。
203:计算散斑图经IVMD分解后分量的直方图,根据直方图判断模态分量中是否包含噪声分量,对噪声分量继续进行IVMD分解,并计算二次分解后分量的直方图,提取二次无噪声分量,该步骤的详细操作为:
1)计算分解后分量的直方图:
h=imhist(uk(x))
其中,h为分解后的模态分量的直方图,imhist为计算分解后分量的直方图函数。
2)计算归一化的直方图:
Figure BDA0001506537740000071
其中,p是归一化的直方图,numel为分解后分量的像素总数。
3)分析归一化直方图,判断其是否符合正态分布,直方图符合正态分布的模态分量为含噪分量,否则为无噪声分量。
4)对含噪的模态分量继续进行IVMD分解,对二次分解后的分量计算其直方图并分析直方图,根据直方图的分布选取无噪声分量。
204:利用希尔伯特法计算降噪后散斑图的相位,对散斑变形前后的相位进行相减获得DSPI相位图,该步骤的详细操作为:
1)首先将包含变形信息的主分量C(x,y)进行希尔伯特变换,并通过反正切获得包裹相位;
Figure BDA0001506537740000072
其中,ψ(x)为相位,fo为载频,Re{}与Im{}分别代表实部和虚部。
2)将变形前后的相位进行相减;
δ=ψ'-ψ
其中,δ为DSPI相位,参见图5,ψ'为变形后的相位图,ψ为变形前的相位图。
205:对DSPI相位采用逐行逐列的方式进行展开(即解包裹),并根据位移与相位的线性关系提取位移,该步骤的详细操作为:
1)对DSPI相位进行展开,相位展开的数学表达式:
Figure BDA0001506537740000073
2)根据相位值与位移的关系计算位移,数学表达式如下:
Figure BDA0001506537740000074
其中d为位移,λ为波长。
综上所述,本发明实施例通过上述步骤201-步骤205实现了包含大量噪声的散斑图前提下,采用IVMD对其进行自适应分解,根据直方图分布特性滤除噪声干扰,提高相位图的信噪比,利用希尔伯特变化法获取精确的相位信息。
本发明实施例对各器件的型号除做特殊说明的以外,其他器件的型号不做限制,只要能完成上述功能的器件均可。
本领域技术人员可以理解附图只是一个优选实施例的示意图,上述本发明实施例序号仅仅为了描述,不代表实施例的优劣。
以上所述仅为本发明的较佳实施例,并不用以限制本发明,凡在本发明的精神和原则之内,所作的任何修改、等同替换、改进等,均应包含在本发明的保护范围之内。

Claims (4)

1.一种基于改进的变分模态分解的DSPI相位提取方法,其特征在于,所述方法包括以下步骤:
1)采集DSPI相位图,根据相位图长度和宽度估算分解后的模态分量的个数,利用正交指标提取最佳模态数量;
2)根据最佳模态数量对DSPI相位图进行变分模态分解,获得一系列的模态函数分量;
3)计算模态函数分量的直方图,利用直方图分析分解后的模态分量是否包含噪声,对含噪分量继续进行改进的变分模态分解,分别计算二次分解后分量的直方图,根据直方图获得二次分解后的无噪声分量;
4)将二次分解后的无噪声分量和一次分解后的无噪声分量进行重构,获得降噪后的散斑图;
5)采用希尔伯特变化法计算变形前后的相位图,对其进行相减得到DSPI相位图;
6)对相位进行解包裹及位移提取;
其中,改进的变分模态分解为:
将模态分量定义为具有不同中心频率的有限带宽;根据不同模态的带宽之和最小原则构造约束变分方程;
针对约束变分方程,引入增广Lagrange函数将约束变分方程转化为非约束变分方程;
采用乘法算子交替方向法计算增广Lagrange函数的鞍点,即约束变分方程的最优解;
不断更新中心频率及模态分量,当模态分量满足迭代停止条件时,停止更新,输出分解后的一系列的模态函数分量;
计算分解后的模态分量直方图,若分解后的模态分量直方图符合正态则为含噪模态分量,反之则为无噪声分量;
对含噪模态分量继续进行模态分解,获得一系列分量;
计算二次分解后的直方图,根据直方图的分布提取二次分解后的无噪声分量。
2.根据权利要求1所述的一种基于改进的变分模态分解的DSPI相位提取方法,其特征在于,上述步骤5)具体为:
采用希尔伯特法计算降噪后散斑图的相位;
利用变形前后的散斑相位进行相减,相减结果为DSPI相位图。
3.根据权利要求1所述的一种基于改进的变分模态分解的DSPI相位提取方法,其特征在于,上述步骤1)具体为:
判断相位图的长度和宽度,两者相等的时候,利用长度或宽度可以计算出来模态分量的个数;
长度和宽度不相等的时候,首先根据尺度估算模态数量的范围,分别计算分解后的正交性,选取正交最小时对应的分量个数,即为最佳模态数量。
4.根据权利要求1所述的一种基于改进的变分模态分解的DSPI相位提取方法,其特征在于,上述步骤6)具体为:
对相减后的相位采用逐行逐列的方法进行展开;
根据展开后的相位与位移的线性关系计算位移。
CN201711330483.1A 2017-12-13 2017-12-13 一种基于改进的变分模态分解的dspi相位提取方法 Active CN108053379B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201711330483.1A CN108053379B (zh) 2017-12-13 2017-12-13 一种基于改进的变分模态分解的dspi相位提取方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201711330483.1A CN108053379B (zh) 2017-12-13 2017-12-13 一种基于改进的变分模态分解的dspi相位提取方法

Publications (2)

Publication Number Publication Date
CN108053379A CN108053379A (zh) 2018-05-18
CN108053379B true CN108053379B (zh) 2021-06-01

Family

ID=62132438

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201711330483.1A Active CN108053379B (zh) 2017-12-13 2017-12-13 一种基于改进的变分模态分解的dspi相位提取方法

Country Status (1)

Country Link
CN (1) CN108053379B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113112075B (zh) * 2021-04-14 2022-04-05 哈尔滨工程大学 一种基于vmd和narx的内燃机噪声预测方法
CN114486259A (zh) * 2022-01-05 2022-05-13 电子科技大学 优化变分模态分解的分布式光纤声传感系统的信号处理方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6753664B2 (en) * 2001-03-22 2004-06-22 Creo Products Inc. Method for linearization of an actuator via force gradient modification
CN104395175A (zh) * 2012-09-11 2015-03-04 日本精工株式会社 车载电子控制装置
CN105181300A (zh) * 2015-09-06 2015-12-23 天津大学 低相干频域干涉图的自适应干涉项提取方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20140046603A1 (en) * 2012-08-08 2014-02-13 International Business Machines Corporation Estimating losses in a smart fluid-distribution system

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6753664B2 (en) * 2001-03-22 2004-06-22 Creo Products Inc. Method for linearization of an actuator via force gradient modification
CN104395175A (zh) * 2012-09-11 2015-03-04 日本精工株式会社 车载电子控制装置
CN105181300A (zh) * 2015-09-06 2015-12-23 天津大学 低相干频域干涉图的自适应干涉项提取方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Time-Frequency Analysis Based on Improved Variational Mode Decomposition and Teager Energy Operator for Rotor System Fault Diagnosis;Shangkun Liu 等;《Mathematical Problems in Engineering》;20161231;第1-9页 *
基于改进变分模态分解的旋转机械故障时频分析方法;刘尚坤 等;《振动工程学报》;20161231;第1119-1126页 *

Also Published As

Publication number Publication date
CN108053379A (zh) 2018-05-18

Similar Documents

Publication Publication Date Title
Trusiak et al. Advanced processing of optical fringe patterns by automated selective reconstruction and enhanced fast empirical mode decomposition
Youngworth et al. An overview of power spectral density (PSD) calculations
Uss et al. Image informative maps for component-wise estimating parameters of signal-dependent noise
WO2007043036A1 (en) Method and system for object reconstruction
KR20130119910A (ko) 광학적 3차원 측량의 자기 적응 윈도우 푸리에 위상추출방법
Shukla et al. Generalized fractional filter-based algorithm for image denoising
CN107563969A (zh) 基于变分模态分解的dspi相位滤波方法
Budianto et al. Marker encoded fringe projection profilometry for efficient 3D model acquisition
CN108053379B (zh) 一种基于改进的变分模态分解的dspi相位提取方法
US8189958B2 (en) Method of fast image reconstruction
Pérez et al. Lightfield recovery from its focal stack
Dou et al. Anti-forensics of diffusion-based image inpainting
Zou et al. Color fringe-projected technique for measuring dynamic objects based on bidimensional empirical mode decomposition
CN107907542B (zh) 一种ivmd及能量估计相结合的dspi相位滤波方法
US20190320100A1 (en) Electromagnetic wave phase/amplitude generation device, electromagnetic wave phase/amplitude generation method, and electromagnetic wave phase/amplitude generation program
Fu et al. Real-time three-dimensional shape measurement based on color binary fringe projection
CN112165570B (zh) 一种基于计算鬼成像的多深度目标对焦方法
Zhong et al. Noise reduction in modulation measurement profilometry based on the wavelet transform method
Zhou et al. Fourier transform profilometry based on convolution neural network
Diaz et al. Estimating photometric properties from image collections
Florack A spatio-frequency trade-off scale for scale-space filtering
CN113011107A (zh) 基于深度卷积神经网络的一维光纤传感信号相位恢复方法
Yagnik et al. 3D shape extraction of human face in presence of facial hair: A profilometric approach
Chen et al. Fast two-step phase-shifting method for measuring the three-dimensional contour of objects
Wang et al. Suppression method for strong interference from stray light in 3D imaging system of structured light

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