CN111207780B - 一种基于cft的牛顿环参数估计方法 - Google Patents

一种基于cft的牛顿环参数估计方法 Download PDF

Info

Publication number
CN111207780B
CN111207780B CN201911363780.5A CN201911363780A CN111207780B CN 111207780 B CN111207780 B CN 111207780B CN 201911363780 A CN201911363780 A CN 201911363780A CN 111207780 B CN111207780 B CN 111207780B
Authority
CN
China
Prior art keywords
newton
ring
cft
calculating
parameter estimation
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
CN201911363780.5A
Other languages
English (en)
Other versions
CN111207780A (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.)
Inner Mongolia University of Science and Technology
Original Assignee
Inner Mongolia 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 Inner Mongolia University of Science and Technology filed Critical Inner Mongolia University of Science and Technology
Priority to CN201911363780.5A priority Critical patent/CN111207780B/zh
Publication of CN111207780A publication Critical patent/CN111207780A/zh
Application granted granted Critical
Publication of CN111207780B publication Critical patent/CN111207780B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

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
    • G01D5/00—Mechanical means for transferring the output of a sensing member; Means for converting the output of a sensing member to another variable where the form or nature of the sensing member does not constrain the means for converting; Transducers not specially adapted for a specific variable
    • G01D5/26—Mechanical means for transferring the output of a sensing member; Means for converting the output of a sensing member to another variable where the form or nature of the sensing member does not constrain the means for converting; Transducers not specially adapted for a specific variable characterised by optical transfer means, i.e. using infrared, visible, or ultraviolet light
    • G01D5/266—Mechanical means for transferring the output of a sensing member; Means for converting the output of a sensing member to another variable where the form or nature of the sensing member does not constrain the means for converting; Transducers not specially adapted for a specific variable characterised by optical transfer means, i.e. using infrared, visible, or ultraviolet light by interferometric means
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Algebra (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Databases & Information Systems (AREA)
  • Length Measuring Devices By Optical Means (AREA)
  • Testing Of Optical Devices Or Fibers (AREA)

Abstract

本发明公开了一种基于CFT的牛顿环参数估计方法，属于光学测量领域，该方法包括：将牛顿环的灰度值变换到[0,1]，计算相邻两行的差，从上述结果中随机选取K行，应用离散时间chirp‑Fourier变换计算其幅度谱，将此幅度谱作为ALO的代价函数，优化得到K组最优值，利用离群点检测算法去掉其中的离群点，并基于剩余结果计算其平均值得到最优值点，根据牛顿环在CFT域的冲激特性，估计相位参数，进而得出曲率半径的测量值，因此，本发明比传统的牛顿环测量方法具有更好的噪声和遮挡鲁棒性，且具有计算复杂度低、精度高、速度快等优点。

Description

一种基于CFT的牛顿环参数估计方法
技术领域
本发明涉及一种基于CFT的牛顿环参数估计方法，属于光学测量领域。
背景技术
干涉测量是通过干涉条纹获得物理参数的光学测量技术，在工程计量、光学计量、地震、光谱学等领域具有重要的应用价值。牛顿环作为经典的二次相位干涉条纹之一，在干涉测量中有着广泛的应用，如判断透镜表面的凹凸、测量透镜表面曲率半径和液体折射率等。
随着光电技术的出现，现代信号处理理论可以应用到干涉测量中，形成了信号处理与传统干涉测量交叉的新的研究领域。这一新的交叉研究领域应用了一些现代信号处理工具来估计干涉条纹的参数，如傅里叶变换、加窗傅里叶变换、调频傅里叶变换、S变换、Kalmana滤波器、小波变换、极坐标变换和分数傅里叶变换。这些方法改善了传统干涉测量中的噪声和障碍物鲁棒性，但这些方法仍然存在一些缺点。例如，在精细确定条纹参数时，基于傅里叶变换的方法需要三次迭代来分离一个边带，迭代过程耗时较长；基于小波变换的方法可以提取任意局部范围内干涉条纹的频率特性，但仍存在阈值选择等问题，窗口大小的选择和耗时；基于小波变换的方法具有较好的时频局部分析能力，但小波函数的选取仍是一个难点；基于FRFT的方法采用二维干涉条纹和遍历搜索来获得匹配角，速度和精度不能同时得到保证。
作为傅立叶变换的推广，Chirp-Fourier变换 (CFT)将时间信号变换到频率和调频率平面上的二元函数，可以用来匹配具有多个分量的chirp型信号，这些良好的性质使得CFT成为信号处理和光学测量领域中的一个重要工具。
发明内容
本发明的目的是提供一种基于CFT的牛顿环参数估计方法；该方法不需要相位解缠绕，可以直接准确地估计牛顿环的相位导数，进而可以应用到球面曲率半径的测量；同时本发明比传统的牛顿环测量方法具有更好的噪声和遮挡鲁棒性，且具有计算复杂度低、精度高、速度快等优点。
本发明的目的是通过下述技术方案实现的：
一种基于CFT的牛顿环参数估计方法，包含如下步骤：
步骤一、 将
Figure 747734DEST_PATH_IMAGE002
牛顿环的灰度值变换到
Figure 100002_DEST_PATH_IMAGE004
, 记为
Figure 100002_DEST_PATH_IMAGE006
, 其中
Figure 100002_DEST_PATH_IMAGE008
；
步骤二、 计算
Figure 100002_DEST_PATH_IMAGE010
中相邻两行的差，记为
Figure 100002_DEST_PATH_IMAGE012
，其中
Figure 100002_DEST_PATH_IMAGE014
；
步骤三、 从
Figure 100002_DEST_PATH_IMAGE016
中随机选取K行，分别利用DTCFT计算其幅度谱, 即
Figure 100002_DEST_PATH_IMAGE018
其中
Figure DEST_PATH_IMAGE020
表示DTCFT算子；
步骤四、 令
Figure 100002_DEST_PATH_IMAGE022
为数值优化算法ALO的代价函数, 优化变量为
Figure 100002_DEST_PATH_IMAGE024
和
Figure 100002_DEST_PATH_IMAGE026
,设置ALO的初始参数并运行得到全局最优值
Figure 100002_DEST_PATH_IMAGE028
, 即
Figure 100002_DEST_PATH_IMAGE030
步骤五、利用离群点检测算法去掉其中的离群点，并基于剩余结果计算其平均值得到
Figure 100002_DEST_PATH_IMAGE032
；
步骤六、牛顿环的相位导数可由下述方程给出，即
Figure DEST_PATH_IMAGE034
其中
Figure DEST_PATH_IMAGE036
表示牛顿环的相位；
步骤七、利用相位导数中所蕴含的信息，计算球面的曲率半径，即
Figure DEST_PATH_IMAGE038
其中
Figure DEST_PATH_IMAGE040
表示入射光的波长，
Figure DEST_PATH_IMAGE042
表示牛顿环的物理尺寸，
Figure DEST_PATH_IMAGE044
表示球面曲率半径的估计值。
进一步，牛顿环的每一行或列可以看作一个包含3个分量的chirp信号，其中常数分量对参数估计是无用的，但是它的存在会影响估计的稳定性。因此，本发明提出了一种简单实用的常数分量去除办法，通过计算牛顿环相邻两行或两列的差，有效降低了常数分量对估计结果的影响。另外，由于相邻两行或两列的差仍是chirp信号，且保持初始频率和调频率不变，同样可用于参数估计。
进一步，本发明使用了异常点检测算法，仅需使用牛顿环的某些行来估计参数，大大降低了计算量。
本发明的有益效果在于：
1. 本发明使用离散时间CFT来估计牛顿环的参数，突破了离散CFT对时频分辨率的限制，提高了参数估计的精度。
2. 由于异常点检测算法的应用，该方法仅使用牛顿环的某些行来估计参数，大大降低了计算量。
3. 建立了牛顿环参数估计的优化模型，并进一步应用ALO算法搜索最优参数，有效地改善了现有一些方法的不足。
与传统的牛顿环干涉法相比，基于CFT的测量方法具有更好的噪声和遮挡鲁棒性。与基于FRFT的方法相比，该测量方法具有计算复杂度低、精度高、速度快等优点。
附图说明
图1—512x512的仿真牛顿环。
图2—仿真牛顿环的相位导数及其误差。
图3—720x720的实际牛顿环。
图4—实际牛顿环的相位导数及其误差。
图5—图1的
Figure DEST_PATH_IMAGE046
和相应的CFT幅度谱。
图6—高斯白噪声鲁棒性分析。
图7—遮挡鲁棒性分析。
图8—实际牛顿环。
图9—基于720x720实际牛顿环的性能分析。
图10—2幅中心不可见的实际牛顿环及其CFT幅度谱。
具体实施方式
下面结合具体的实例和附图对本发明做详细说明：
为了验证所提方法在参数估计中的性能，基于仿真和实际牛顿环进行了实验分析。数值优化算法ALO的初始化参数为
Figure DEST_PATH_IMAGE048
其中
Figure DEST_PATH_IMAGE050
表示种群数目，
Figure DEST_PATH_IMAGE052
表示最大迭代次数，LB和UB表示优化变量的下界和上界，
Figure DEST_PATH_IMAGE054
表示优化变量的个数。
对于如图1中所示的仿真牛顿环，其相位参数为，估计得到的相位参数为
Figure DEST_PATH_IMAGE056
。为了更直观地显示估计结果，估计的相位
Figure DEST_PATH_IMAGE058
导数及其误差如图2 (a)和(b)所示。对于实际的牛顿环（见图3），实际牛顿环的相位参数是
Figure DEST_PATH_IMAGE060
，对应的估计相位参数为
Figure DEST_PATH_IMAGE062
。估计的相位导数及其误差如图4 (a)和(b)所示。这些数据表明，仿真和实际牛顿环的相位导数可以在不使用解缠绕方法的情况下直接准确地估计出来。
基于牛顿环参数所包含的信息，基于CFT的牛顿环参数估计方法可以进一步应用于测量球面的曲率半径。对于图1中的纯净仿真牛顿环，其
Figure DEST_PATH_IMAGE064
和相应的CFT幅度谱如图5所示，从图中可以看出明显的冲激特性。该方法得到的曲率半径估计值分别为0.8606m和0.8595m，相对误差分别为0.0668%和0.0591%，处理时间约为3.8s，基于FRFT方法得到的曲率半径相对误差为0.22%，处理时间约为990s。这些数据说明了，基于CFT的测量方法的效率远高于基于FRFT的方法。
由于在实际操作中，光学仪器和人工操作不可避免地会造成污染，为了更好地模拟实际情况，有必要在牛顿环受到污染的情况下对该方法的性能进行分析。主要从以下两点分析：1）污染源为高斯白噪声，信噪比在-10db到10db之间变化；2）污染源为较大的噪声点，由半径为100像素、中心位置不同的圆形遮挡物模拟。对于以上两点，所得结果列于表1。通过分析这些数据，可以得到以下结论：
表1 基于污染牛顿环的性能分析
Figure DEST_PATH_IMAGE066
1）该方法对噪声具有较强的鲁棒性，当
Figure DEST_PATH_IMAGE068
时，噪声对测量结果的影响较小。为了将上述结果可视化，将图1中的牛顿环添加了
Figure DEST_PATH_IMAGE070
的高斯白噪声(见图6 (a))，未添加噪声和添加噪声之后的
Figure DEST_PATH_IMAGE072
如图6(b)所示。应用CFT到添加噪声之后的
Figure DEST_PATH_IMAGE074
，得到的幅度谱及其投影如图6(c)和(d)所示。从图中可得，虽然牛顿环几乎被噪声湮没，但在图6(d)中可以看到明显的冲激。因此，基于CFT的测量方法对高斯白噪声具有较好的鲁棒性。
2）该方法对遮挡具有很强的鲁棒性。与纯净牛顿环的结果相比，遮挡物对估计结果的影响很小。为了将上述结果可视化，牛顿环被圆形障碍物遮挡，其中心位于（256, 256）像素(见图7 (a))。图7(b)中展示了被遮挡牛顿环的
Figure DEST_PATH_IMAGE076
和
Figure DEST_PATH_IMAGE078
。由此图可以看出，由于障碍物的存在，
Figure DEST_PATH_IMAGE080
不再是chirp信号，而
Figure DEST_PATH_IMAGE082
仍然是chirp信号。因此，使用
Figure DEST_PATH_IMAGE084
来测量球面的曲率半径，其CFT谱和投影分别如图7(c)和(d)所示，图中可见明显的冲激。该方法利用牛顿环的某些行来测量曲率半径，由此我们可以选择那些没有被遮挡的行，这也说明基于CFT的方法对遮挡具有很强的鲁棒性。
3）与基于FRFT的方法相比，二者都具有较强的遮挡鲁棒性，但基于CFT的方法具有较好的噪声鲁棒性和较高的精度，特别是测量时间大大缩短了，有利于曲率半径的实时快速测量。
进而，使用如图8所示的实际牛顿环来验证所提测量方法的性能，分辨率分别为720×720像素（像素尺寸为5.9
Figure DEST_PATH_IMAGE086
×5.9
Figure 484395DEST_PATH_IMAGE086
）、1080×1080像素（像素尺寸为3.9
Figure 410763DEST_PATH_IMAGE086
×3.9
Figure 490714DEST_PATH_IMAGE086
）和1200×1200像素（像素尺寸为6.4
Figure 672949DEST_PATH_IMAGE086
×6.4
Figure 471140DEST_PATH_IMAGE086
）。下面以分辨率为720×720的牛顿环为例，对分析结果进行展示。在实际牛顿环的采集过程中，由于环境、人工操作等因素的影响，在实际牛顿环中存在许多噪声点。如图9（b）所示，图9(a)的第一行和第一列不再是标准chirp信号。但是其CFT幅度谱和投影（见图9（c）和（d））中仍然可以看到明显的冲激，这进一步说明了基于CFT的方法具有更好的噪声和障碍鲁棒性。对于图9（a-c），基于CFT的测量方法得到的曲率半径误差约为0.4964%、0.2675%和0.1328%。对于基于FRFT的方法，曲率半径误差大约为1.27%、0.9%和0.34%。
为了分析所提方法对中心不可见的实际牛顿环的性能，从图8(c)中裁剪了两幅512×512的实际牛顿环。图10(a)和(b)展示了中心不可见的牛顿环及其CFT幅度谱，从其幅度谱中可观察到明显的冲激。因此，对于中心不可见的实际牛顿环，基于CFT的方法仍然具有测量曲率半径的能力，得到的曲率半径分别为0.8635m和0.8625m，相对误差分别为0.4091%和0.2854%。这进一步说明，由于基于CFT的曲率半径测量方法使用牛顿环的本质特征，而不是传统牛顿环测量法所使用的物理特征(明环或暗环半径)，所以在球面曲率半径测量中具有更广泛的适用性。

Claims (2)

1.一种基于CFT的牛顿环参数估计方法，其特征在于：
步骤一、 将
Figure DEST_PATH_IMAGE001
牛顿环的灰度值变换到
Figure DEST_PATH_IMAGE002
, 记为
Figure DEST_PATH_IMAGE003
, 其中
Figure DEST_PATH_IMAGE004
；
步骤二、 计算
Figure DEST_PATH_IMAGE005
中相邻两行的差，记为
Figure DEST_PATH_IMAGE006
，其中
Figure DEST_PATH_IMAGE007
；
步骤三、 从
Figure DEST_PATH_IMAGE008
中随机选取K行，分别利用DTCFT计算其幅度谱, 即
Figure DEST_PATH_IMAGE009
其中
Figure DEST_PATH_IMAGE010
表示DTCFT算子，
Figure DEST_PATH_IMAGE012
和
Figure DEST_PATH_IMAGE014
分别表示频率和调频率；
步骤四、 令
Figure DEST_PATH_IMAGE015
为数值优化算法ALO的代价函数, 优化变量为
Figure DEST_PATH_IMAGE016
和
Figure DEST_PATH_IMAGE017
, 设置ALO的初始参数并运行得到全局最优值
Figure DEST_PATH_IMAGE018
, 即
Figure DEST_PATH_IMAGE019
其中
Figure DEST_PATH_IMAGE021
表示牛顿环的物理尺寸；
步骤五、利用离群点检测算法去掉其中的离群点，并基于剩余结果计算其平均值得到
Figure DEST_PATH_IMAGE022
；
步骤六、牛顿环的相位导数可由下述方程给出，即
Figure DEST_PATH_IMAGE023
其中
Figure DEST_PATH_IMAGE024
表示牛顿环的相位，
Figure DEST_PATH_IMAGE026
和
Figure DEST_PATH_IMAGE028
分别表示初始频率和调频率；
步骤七、利用相位导数中所蕴含的曲率半径信息，计算球面的曲率半径，得
Figure DEST_PATH_IMAGE030
其中
Figure DEST_PATH_IMAGE031
表示入射光的波长，
Figure DEST_PATH_IMAGE032
表示牛顿环的物理尺寸，
Figure DEST_PATH_IMAGE033
表示球面曲率半径的估计值。
2.如权利要求1所述的基于CFT的牛顿环参数估计方法，其特征在于：牛顿环的每一行或列可以看作一个包含3个分量的chirp信号；由于相邻两行或两列的差仍是chirp信号，且保持初始频率和调频率不变，同样可用于参数估计。
CN201911363780.5A 2019-12-26 2019-12-26 一种基于cft的牛顿环参数估计方法 Active CN111207780B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911363780.5A CN111207780B (zh) 2019-12-26 2019-12-26 一种基于cft的牛顿环参数估计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911363780.5A CN111207780B (zh) 2019-12-26 2019-12-26 一种基于cft的牛顿环参数估计方法

Publications (2)

Publication Number Publication Date
CN111207780A CN111207780A (zh) 2020-05-29
CN111207780B true CN111207780B (zh) 2021-07-16

Family

ID=70788391

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911363780.5A Active CN111207780B (zh) 2019-12-26 2019-12-26 一种基于cft的牛顿环参数估计方法

Country Status (1)

Country Link
CN (1) CN111207780B (zh)

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2509772B2 (ja) * 1991-10-25 1996-06-26 株式会社三協精機製作所 光コネクタの端面検査装置
CN101504717B (zh) * 2008-07-28 2012-07-11 上海高德威智能交通系统有限公司 特征区域的定位方法、车身深浅色与车身颜色的识别方法
CN101881821B (zh) * 2010-06-28 2012-11-28 北京理工大学 一种分数阶傅里叶域信道化接收方法
CN104931982B (zh) * 2015-05-29 2017-03-22 西安电子科技大学 高动态弱信号下基于dcft的块补零码捕获方法
CN106092158B (zh) * 2016-08-19 2018-08-21 北京理工大学 物理参数估计方法、装置和电子设备
CN107703626B (zh) * 2017-10-19 2020-04-17 北京理工大学 一种基于变换材料的光学分数阶傅里叶变换的装置和设计方法

Also Published As

Publication number Publication date
CN111207780A (zh) 2020-05-29

Similar Documents

Publication Publication Date Title
CN104006765B (zh) 单幅载频干涉条纹相位提取方法及检测装置
JP6752212B2 (ja) 画像システムとその使用方法
Salvi et al. A state of the art in structured light patterns for surface profilometry
KR102518289B1 (ko) 간섭계의 광학 성능을 최적화하기 위한 방법 및 장치
CN114993452B (zh) 基于宽带相位运动放大的结构微小振动测量方法与系统
CN110726372B (zh) 一种准确处理单缝衍射图像方法
CN111207780B (zh) 一种基于cft的牛顿环参数估计方法
CN117031684B (zh) 一种数字全息成像自动聚焦方法及装置
CN112381731A (zh) 一种基于图像去噪的单帧条纹图像相位分析方法及系统
CN103383353B (zh) 一种基于光学涡旋的动态散斑测试方法
Gladines et al. A phase correlation based peak detection method for accurate shape from focus measurements
Chen et al. A Voronoi-Diagram-based method for centerline extraction in 3D industrial line-laser reconstruction using a graph-centrality-based pruning algorithm
Zhong et al. Noise reduction in modulation measurement profilometry based on the wavelet transform method
Zhou et al. Speckle-correlation-based non-line-of-sight imaging under white-light illumination
TWI637166B (zh) 微分相位對比顯微系統與方法
Zhong et al. Elimination of nonlinearity in modulation measurement profilometry by wavelet transform
CN110608677B (zh) 一种基于圆光栅径向剪切干涉仪的三维位移测量方法
CN116243242A (zh) 基于ф-otdr的自适应滤波一体化降噪高精度3d定位方法
CN103471528B (zh) 一种干涉条纹倾斜角测量方法
Guo et al. Chirp-Fourier transform for quadratic phase interference fringe analysis: Principles, method and application
Osten Digital processing and evaluation of fringe patterns in optical metrology and non-destructive testing
Silvetti et al. Quadratic self-correlation: An improved method for computing local fractal dimension in remote sensing imagery
TWI729534B (zh) 量測物體形貌的方法
He et al. Enhanced accuracy in 3D structured illumination microscopy through binary encoding with accelerated speed using sampling Moiré
Jacques et al. Refractive index map reconstruction in optical deflectometry using total-variation regularization

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