CN110260812B - 一种基于欠定盲源分离双通道光学三维干涉方法及系统 - Google Patents

一种基于欠定盲源分离双通道光学三维干涉方法及系统 Download PDF

Info

Publication number
CN110260812B
CN110260812B CN201910400871.5A CN201910400871A CN110260812B CN 110260812 B CN110260812 B CN 110260812B CN 201910400871 A CN201910400871 A CN 201910400871A CN 110260812 B CN110260812 B CN 110260812B
Authority
CN
China
Prior art keywords
interference
signal
signals
frequency
source
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.)
Expired - Fee Related
Application number
CN201910400871.5A
Other languages
English (en)
Other versions
CN110260812A (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.)
Guangdong University of Technology
Original Assignee
Guangdong 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 Guangdong University of Technology filed Critical Guangdong University of Technology
Priority to CN201910400871.5A priority Critical patent/CN110260812B/zh
Publication of CN110260812A publication Critical patent/CN110260812A/zh
Priority to US16/875,017 priority patent/US11060849B2/en
Application granted granted Critical
Publication of CN110260812B publication Critical patent/CN110260812B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
    • G01B9/00—Measuring instruments characterised by the use of optical techniques
    • G01B9/02—Interferometers
    • G01B9/02001—Interferometers characterised by controlling or generating intrinsic radiation properties
    • G01B9/02002—Interferometers characterised by controlling or generating intrinsic radiation properties using two or more frequencies
    • G01B9/02004—Interferometers characterised by controlling or generating intrinsic radiation properties using two or more frequencies using frequency scans
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
    • G01B11/00—Measuring arrangements characterised by the use of optical techniques
    • G01B11/24—Measuring arrangements characterised by the use of optical techniques for measuring contours or curvatures
    • G01B11/2441—Measuring arrangements characterised by the use of optical techniques for measuring contours or curvatures using interferometry
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
    • G01B9/00—Measuring instruments characterised by the use of optical techniques
    • G01B9/02—Interferometers
    • G01B9/02015—Interferometers characterised by the beam path configuration
    • G01B9/02017—Interferometers characterised by the beam path configuration with multiple interactions between the target object and light beams, e.g. beam reflections occurring from different locations
    • G01B9/02021—Interferometers characterised by the beam path configuration with multiple interactions between the target object and light beams, e.g. beam reflections occurring from different locations contacting different faces of object, e.g. opposite faces
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
    • G01B9/00—Measuring instruments characterised by the use of optical techniques
    • G01B9/02—Interferometers
    • G01B9/02015—Interferometers characterised by the beam path configuration
    • G01B9/02027—Two or more interferometric channels or interferometers
    • G01B9/02028—Two or more reference or object arms in one interferometer
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
    • G01B9/00—Measuring instruments characterised by the use of optical techniques
    • G01B9/02—Interferometers
    • G01B9/02015—Interferometers characterised by the beam path configuration
    • G01B9/02032—Interferometers characterised by the beam path configuration generating a spatial carrier frequency, e.g. by creating lateral or angular offset between reference and object beam
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
    • G01B9/00—Measuring instruments characterised by the use of optical techniques
    • G01B9/02—Interferometers
    • G01B9/02041—Interferometers characterised by particular imaging or detection techniques
    • G01B9/02044—Imaging in the frequency domain, e.g. by using a spectrometer
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
    • G01B9/00—Measuring instruments characterised by the use of optical techniques
    • G01B9/02—Interferometers
    • G01B9/02055—Reduction or prevention of errors; Testing; Calibration
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
    • G01B9/00—Measuring instruments characterised by the use of optical techniques
    • G01B9/02—Interferometers
    • G01B9/02083—Interferometers characterised by particular signal processing and presentation
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
    • G01B9/00—Measuring instruments characterised by the use of optical techniques
    • G01B9/02—Interferometers
    • G01B9/02083—Interferometers characterised by particular signal processing and presentation
    • G01B9/02084—Processing in the Fourier or frequency domain when not imaged in the frequency domain
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01J—MEASUREMENT OF INTENSITY, VELOCITY, SPECTRAL CONTENT, POLARISATION, PHASE OR PULSE CHARACTERISTICS OF INFRARED, VISIBLE OR ULTRAVIOLET LIGHT; COLORIMETRY; RADIATION PYROMETRY
    • G01J3/00—Spectrometry; Spectrophotometry; Monochromators; Measuring colours
    • G01J3/28—Investigating the spectrum
    • G01J3/2823—Imaging spectrometer

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Engineering & Computer Science (AREA)
  • Signal Processing (AREA)
  • Mathematical Physics (AREA)
  • Instruments For Measurement Of Length By Optical Means (AREA)
  • Investigating Or Analysing Materials By Optical Means (AREA)

Abstract

本发明公开了一种基于欠定盲源分离双通道光学三维干涉方法，通过CCD相机采集到的干涉数据，盲分离出载玻片被测件表面之间的干涉信号，求解干涉信号参数，包括：干涉信号幅频和干涉信号相频；基于双通道光学三维迈克尔逊式干涉实验，通过K均值聚类算法得到混合矩阵的估计及通过L1范数最短路径法实现源信号的恢复；最终实现激光波数扫描通过CCD相机采集的光强值准确地盲分离出四表面干涉信号，实现四表面干涉信号的盲分离；相对于普通的干涉测量方法无法准确地识别出各个峰值位置，本发明利用双通道检测和采样，通过盲源分离，使得信号即使面对较为强烈的噪声，也可以检测出各个干涉信号所在位置。

Description

一种基于欠定盲源分离双通道光学三维干涉方法及系统
技术领域
本发明涉及双通道激光波数扫描三维迈克尔逊式干涉技术领域，具体涉及一种基于欠定盲源分离双通道光学三维干涉方法及系统。
背景技术
干涉指满足一定条件的两列相干波相遇叠加，在叠加区域某些点的振动始终加强，某些点的振动始终减弱，即在干涉区域内振动强度有稳定的空间分布。干涉测量法在各种参数的测量中，具有很高的测试灵敏度和准确度，是一种高精度的测量方法。
盲源分离是指在源信号以及传输通道过程未知的条件下，仅仅根据观测到的混合信号来恢复出源信号，当观测信号个数小于源信号个数时，为欠定盲源分离。
现有技术中，由于观测信号个数少于源信号个数，无法采用求逆矩阵的形式求解源信号；如果噪声波动过于强烈，而干涉信号的幅值又较小，则很容易使两个较为接近的干涉信号混叠为一个信号。
发明内容
本发明的目的在于克服现有技术的缺点与不足，提供一种基于欠定盲源分离双通道光学三维干涉方法，该方法基于迈克尔逊式干涉仪使得干涉参考光路与测量光路相互独立，通过在干涉参考光路中设置液晶片，起到压制参考光光强作用。
本发明的另一目的在于提供一种基于欠定盲源分离双通道光学三维干涉系统。
本发明的目的通过下述技术方案实现：
一种基于欠定盲源分离双通道光学三维干涉方法，包括下述步骤：
步骤一，计算机发出控制指令，并通过激光控制器和温控模块对半导体激光器波数输出进行线性调频；
其中所述线性调频公式为：
Figure BDA0002059705130000021
其中k(t)为激光波数，λc为半导体激光器工作时的波长，Δλ为激光输出波长的变化范围，T为激光输出波长对应的结束时刻；
步骤二，激光输出光经透镜L1准直为平行光后，被50:50的正方体分光棱镜分成两束光，一束光经液晶片照射到载玻片的两表面S3和S4，另一束光照射到光楔前后两表面S1和S2上，载玻片两表面S3和S4的散射光束和光楔前后两表面S1和S2的散射光束再次经过50:50的正方体分光棱镜后，形成一条返回光路，返回的光路通过透镜L2聚焦到CCD相机里的数据采集卡，并在数据采集卡上相互叠加形成干涉信号；实验中激光器发出的相干光照射到所有被测表面时，通道上表面反射光强相互叠加，形成M(M-1)/2个峰值的干涉信号；其干涉信号因液晶片所产生的光强不同，在通道中产生的幅值高度就不同，故被分为两路干涉信号I1和I2；
其中，I1的干涉信号公式为：
Figure BDA0002059705130000031
其中Λ1pq(x,y)＝Z1pq(x,y)τ1pq(x,y)；
上式中(x,y)为空间坐标；下标p,q分别代表被测表面SP和Sq；I1为反射光强；Λ1pq为被测表面SP和Sq之间的光程差；Z1pq为被测表面SP和Sq之间沿z方向上的位置差；τ1pq为被测表面SP和Sq之间介质的光学折射率；
Figure BDA0002059705130000032
为被测表面SP和Sq之间的初始相位差；
对上述公式I1作傅里叶变换得：
Figure BDA0002059705130000033
其中，I2的干涉信号公式为：
Figure BDA0002059705130000034
其中Λ2pq(x,y)＝Z2pq(x,y)τ2pq(x,y)；
上式中(x,y)为空间坐标；下标p,q分别代表被测表面SP和Sq；I2为反射光强；Λ2pq为被测表面SP和Sq之间的光程差；Z2pq为被测表面SP和Sq之间沿z方向上的位置差；τ2pq为被测表面SP和Sq之间介质的光学折射率；
Figure BDA0002059705130000041
为被测表面SP和Sq之间的初始相位差；
对上述公式I2作傅里叶变换得：
Figure BDA0002059705130000042
基于I1干涉信号能采集到的幅频特性有M(M-1)/2个峰值，每个峰值对应载玻片被测件表面的干涉信号Spq，记该峰值处的幅频和相频为f1pq和Φ1pq，则得到以下公式：
Figure BDA0002059705130000043
Figure BDA0002059705130000044
基于I2干涉信号能采集到的幅频特性有M(M-1)/2个峰值，每个峰值对应载玻片被测件表面的干涉信号Spq，记该峰值处的幅频和相频为f2pq和Φ2pq，则得到以下公式：
Figure BDA0002059705130000045
Figure BDA0002059705130000046
干涉幅频fpq和干涉相频Φpq均包含载玻片被测表面深度轮廓信息，可以通过它们解调出Λpq；
为作简述方便，省略I1(x,y,k)和I2(x,y,k)中的(x,y)，因(x,y)在整个过程中保持一致，故只对k做划分；上述k均是在时域中，信号经傅里叶变换变换到频域中用f代替，以作区分；
步骤三，由于欠定盲源分离要求信号为稀疏信号，而时域光强信号不满足稀疏性，所以对信号进行快时傅里叶变换FFT变换至频域，形成稀疏信号；
基于观测信号个数I＝2小于源信号个数S＝6，故以K均值聚类方法估计混合矩阵A；
设置聚类个数为源信号个数，即K＝S，K均值聚类方法估计混合矩阵A求解步骤如下：
(1)从M(M-1)/2个峰值中选出k个初始聚类中心，将全部峰值随机分成k类；
(2)计算每个峰值信号到各聚类中心的欧氏距离；所述的欧氏距离公式为：
Figure BDA0002059705130000051
(3)误差函数不在变化，则聚类结束；所述的误差函数为：
Figure BDA0002059705130000052
其中Ψ(cj)代表以cj为中心的列向量集合，d(xi(t),cj)为每个采样信号点xi(t)与聚类中心cj的欧氏距离，当且仅当dcj(xi(t),cj)＝min{d(xi(t),cp)，p＝1,…,k}时，函数获得最优解；
步骤四，基于上述源信号在时域上不满足稀疏性的要求，所以对源信号进行快时傅里叶变换FFT，使源信号在变换后的频域上满足稀疏性，从而在变换后的频域上进行盲源分离，变换后的I1(k)和I2(k)满足频域：
Figure BDA0002059705130000061
其中
Figure BDA0002059705130000062
为频域中I1(f)和I2(f)的观测实验结果；
Figure BDA0002059705130000063
为I1(f)、I2(f)干涉光强中幅频的峰值，即S个数为6；故可以从两路观测信号中分离出六路源信号，则观测信号个数小于源信号个数；A为未知的2×6混合矩阵；A和S(f)均未知且混合矩阵A不可逆；
故所述观测信号向量
Figure BDA0002059705130000064
在幅频中展开可得矩阵：
Figure BDA0002059705130000065
根据稀疏理论，所以在某一频率中只有一个信号取值幅频相对较大；假在f1＝1,...,1024，只有S1(f1)的取值较大，则根据上述公式有：
(1)I1(f1)＝a11S1(f1)，I2(f2)＝a21S1(f1)，将前边两式相除得
I2(f1)/I1(f1)＝a21/a11，由此可得，在观测信号的幅频图上，S1(f1)取值的点均聚集在斜率为a21/a11的直线周围；
(2)同理当S2(fi)，S3(fi)，S4(fi)，S5(fi)，S6(fi)幅值较大时点分别均聚集在斜率为a22/a12，a23/a13，a24/a14，a25/a15，a26/a16的直线周围；
所以根据这6条直线的斜率tanθ的值求得sinθ和cosθ，对观测信号数据点进行聚类划分就能得到混合矩阵的各列向量，即
Figure BDA0002059705130000071
为求解的矩阵A中的一列向量；
步骤五，基于K均值聚类法估计出混合矩阵A后，通过最小化L1范数最短路径法来实现源信号的恢复；操作原理为通过对观测信号进行分解，寻找到最靠近观测信号向量的线性组合来作为对源信号的估计，提取部分源的基矢量来最小干扰地恢复出源信号；该方法是先根据观测信号所在方向，计算其与混合矩阵中列向量所代表方向的差，设定一定的阔值，从中挑选出差值较小的多个代表性方向作为观测信号分解的潜在方向；
根据对源信号的先验性假设，当源信号满足稀疏分布时，绝大部分时刻都存在一个源信号从原点到
Figure BDA0002059705130000072
的一条路径；
根据上述欠定盲源分离模型的矩阵形式，可知观测信号向量
Figure BDA0002059705130000073
其中
Figure BDA0002059705130000074
包括I1(f)和I2(f)，可由基向量aij(i＝1,2)(j＝1,2,3,4,5,6)线性组合，而组合系数为各源信号S1(f)，S2(f)，...，S6(f)；所述
Figure BDA0002059705130000075
和aijsk(f)(K＝1,2,3,4,5,6)首尾相连组成一个封闭几何图形，所有向量aijsk(f)的长度之和即为其系数的绝对值之和
Figure BDA0002059705130000081
因此在所有的可行解中，最小值
Figure BDA0002059705130000082
的解是
Figure BDA0002059705130000083
的最短路径；故在混合矩阵估计出来后，由最大后验法，基于任给的
Figure BDA0002059705130000084
都有
Figure BDA0002059705130000085
此时，源信号恢复就转为L1范数最小路径，从而估计出源信号S；
基于最小二乘的参数估计思路是通过使测量误差的平方和最小，寻找待匹配函数的最优参数值；
步骤六，由于I1、I2的方程
Figure BDA0002059705130000086
中干涉频谱的实部和虚部均包含被测表面光程差信息，因此分别对实部和虚部作最小二乘的方法，故在频域内对干涉测量频谱构建误差方程E(x,y)如下：
Figure BDA0002059705130000087
其中L为傅里叶变换频率点数；
Figure BDA0002059705130000088
为干涉测量频谱，可通过CCD相机拍摄的干涉信号作傅里叶变换得到；
干涉相位和幅值可求解为：
Figure BDA0002059705130000089
Figure BDA00020597051300000810
此时表面SP，Sq之间的干涉信号Spq的波数域光强表达式为：
Figure BDA0002059705130000091
即为干涉信号复数域线性最小二乘算法；复数域线性最小二乘算法仅通过CCD相机采集的干涉信号便可自动分离出各个表面的干涉信号，有助于载玻片干涉信号的提取；另一方面，干涉信号成功地盲分离可以进一步提高激光波数扫描干涉检测的深度轮廓分辨率和减少相位波动误差。
优选地，步骤三中所述稀疏信号是指信号个数在1024个频率中的幅频取值接近零，同时又有六个明显的波峰幅频较大；双通道峰值除了六个明显的波峰外，绝大多数地方平滑下降，知源信号之间是相互统计独立的，满足稀疏特性的要求，所以绝大部分时刻最多只有一个源信号取值占优。
一种基于欠定盲源分离双通道光学三维干涉系统，包括半导体激光器、激光控制器、温控模块、光楔、CCD相机、分光棱镜、计算机、液晶片、透镜L1、透镜L2和载玻片；
其中，计算机发出控制指令通过激光控制器和温控模块对半导体激光器波数输出进行线性调频；激光输出光经透镜L1准直为平行光后，被50:50的正方体分光棱镜分成两束光，一束光经液晶片照射到载玻片两表面S3和S4，另一束照射到光楔前后两表面S1和S2上，载玻片两表面S3和S4的散射光束和光楔前后两表面S1和S2的散射光束再次经过50:50的正方体分光棱镜后，形成一条返回光路，返回光路通过透镜L2聚焦到CCD相机里的数据采集卡，并在数据采集卡上相互叠加形成干涉信号，该信号通过CCD相机里的处理器形成干涉图像，并传入计算机最后以图片形式展示进行分析。
本发明与现有技术相比具有以下的有益效果：
(1)相对于普通的干涉测量方法无法准确地识别出各个峰值位置，本发明利用双通道检测和采样，通过盲源分离，使得信号即使面对较为强烈的噪声，也可以检测出各个干涉信号所在位置；
(2)本发明引入K均值聚类算法以及L1范数最小路径算法分别对欠定矩阵进行求解以实现源信号的分离；
(3)由于欠定盲源分离要求信号为稀疏信号，而时域光强信号不满足稀疏性，本发明对信号进行快时傅里叶变换FFT变换至频域，形成稀疏信号；
(4)本发明观测信号I＝2，即为多通道欠定盲源分离提供了有利的条件；
(5)本发明为解决欠定盲源分离，首先估计混合矩阵，再由混合矩阵恢复出源信号。
附图说明
图1为本发明双通道激光波数扫描三维干涉测量系统示意图；
图2为本发明载玻片(S3，S4)和光楔(S1，S2)盲分离散点图；
图3为本发明载玻片和光楔前后两表面干涉盲分离散点图；
图4为本发明6个峰值盲分离散点图；
图5为本发明载玻片(S3，S4)和光楔(S1，S2)干涉幅频图；
图6为本发明载玻片和光楔前后两表面干涉幅频图；
图7为本发明6个峰值干涉幅频图；
图8为本发明干涉条纹图。
图中附图标记为：1半导体激光器；2、激光控制器；3、温控模块；4、光楔；5、CCD相机；6、分光棱镜；7、计算机；8、液晶片；9、透镜L2；10、透镜L1；11、载玻片。
具体实施方式
下面结合实施例及附图对本发明作进一步详细的描述，但本发明的实施方式不限于此。
本实施例所述是一种双通道激光波数扫描三维干涉测量系统，如图1所示，该系统主要包括半导体激光器，激光控制器，温控模块，光楔，CCD相机，分光棱镜，计算机，液晶片，透镜L1，透镜L2和载玻片。
本实施例成像光路的关键器为CCD相机，要求CCD相机动态范围大，分辨率高以及响应速度快。基于载玻片型被测件受环境和温度影响因素大，实验过程中应保持室内温度恒温且稳定。
本实施例提供一种载玻片干涉实验，确定好载玻片干涉实验的条件后，载玻片干涉实验可分为两个阶段：
第一阶段为数据采集阶段，第二阶段为幅频，相位图片切割阶段。对于第一阶段，打开激光控制器，通过计算机选择恒激光电流模式，设置温控范围：28℃-15℃；温控电流：0.6A；激光电流：100mA；曝光时间：70000μs；采样间隔：50ms；样本数量：500张，从而根据上述实验参数设置，对实验进行数据采集，采集过程中我们可以观测到干涉条纹图，如图2所示；
对于第二阶段，基于第一阶段采样数据完成后，进入MATLAB程序进行横向图片切割，其选择存放形为数据包；数据类型为短整型；进行CZT变换和傅里叶变换，其中傅里叶变换不加窗，所述的傅里叶变换为插值法傅里叶变换。
根据上述实施例进入MATLAB程序进行横向图片切割，其原理是把多张图片的同一行或同一列放在一个mat文件中，方便以后进行图片分析处理。因为激光波束序列的变化是与时间t有关系的，从而可把拍摄张数等同于时间t。
本发明的另一核心是激光波数扫描干涉检测。激光波数扫描干涉检测的光路是基于双通道分光路的迈克尔逊式干涉仪。
本实施例基于迈克尔逊式干涉仪，如图1所示。
本实施例的激光控制器的目的是对温度进行线性调制，从而实现激光波数线性扫描。要求温度控制精度极高。
本实施例基于激光波数扫描对载玻片型被测件进行光学干涉检测的实验表明，仅通过CCD相机拍摄的干涉条纹便可以实现被测件表面之间干涉信号的盲分离。
本实施例的干涉信号盲分离问题可以简化为6个干涉信号的叠加，当深度z方向上存在4个表面时，包括：基于载玻片、光楔两者上，将会有6个干涉信号叠加，其中一个干涉为一个叠加。准确地提取干涉参考面与被测载玻片之间的干涉信号。根据所述设计被测表面之间合适的距离，实现有用的干涉信号分离。
本实施例的实验原理为：
计算机发出控制指令通过激光控制器和温控模块对半导体激光器波数输出进行线性调频；激光输出光经透镜L1准直为平行光后，被50:50的正方体分光棱镜分成两束光，一束光经液晶片照射到载玻片两表面S3和S4，另一束照射到光楔前后两表面S1和S2上，载玻片两表面S3和S4的散射光束和光楔前后两表面S1和S2的散射光束再次经过50:50的正方体分光棱镜后，形成一条返回光路，返回光路通过透镜L2聚焦到CCD相机里的数据采集卡，并在数据采集卡上相互叠加形成干涉信号，该信号通过CCD相机里的处理器形成干涉图像，并传入计算机最后以图片形式展示进行分析。
实验中激光器发出的相干光照射到被测表面S1,S2,S3,S4,时，通道上表面反射光强相互叠加，形成6个峰值的干涉信号。
基于干涉信号能采集到的幅频特性有6个峰值，每个峰值对应载玻片被测件表面的干涉信号Spq，记该峰值处的幅频和相频为fpq和Φpq，则得到以下公式：
Figure BDA0002059705130000131
Figure BDA0002059705130000132
本实施例观测信号个数I＝2小于源信号个数S＝6，故以K-均值聚类方法估计混合矩阵A，成功从载玻片被测件表面的干涉实验中，分离出干涉信号。
本实施例所述的K均值聚类算法是设置聚类个数等于源信号个数，在本实验中即为6，故A为2x6的矩阵。对信号进行快时傅里叶变换(FFT)，使源信号稀疏化，在变换后的频域上进行分离，得到公式：
Figure BDA0002059705130000141
本实施例中K均值聚类算法实验过程如下：
(1)在1024个散点的分布数据中选出6类初始聚类中心，把距离这6类点最接近的其他点归为一类，取当前类的所有点的均值，作为中心点；参见散点图4。
(2)分别计算6个峰值信号到各聚类中心的欧氏距离；所述的欧氏距离公式：
Figure BDA0002059705130000142
其中y1i为通道1上的y坐标，y2i为通道2上的y坐标。
其数据结果参见双通道干涉幅频图7。因欧氏距离是到原点的距离，故y0＝0，以下统一省略，则：
当y11＝39.44,y21＝36.76时，
Figure BDA0002059705130000143
当y12＝22.33,y22＝13.72时，
Figure BDA0002059705130000144
当y13＝24.61,y23＝16.46时，
Figure BDA0002059705130000145
当y14＝27.89,y24＝15.94时，
Figure BDA0002059705130000146
当y15＝21.70,y25＝13.33时，
Figure BDA0002059705130000147
当y16＝22.76,y26＝15.00时，
Figure BDA0002059705130000151
(3)误差函数不在变化，则聚类结束。所述的误差函数：
Figure BDA0002059705130000152
基于K均值聚类法估计出混合矩阵A后，通过最小化L1范数最小路径法来实现源信号的恢复。操作原理为通过对观测信号进行分解，寻找到最靠近观测信号向量的线性组合来作为对源信号的估计，提取部分源的基矢量来最小干扰地恢复出源信号。该方法是先根据观测信号所在方向，计算其与混合矩阵中列向量所代表方向的差，设定一定的阔值，从中挑选出差值较小的多个代表性方向作为观测信号分解的潜在方向。
根据对源信号的先验性假设，当源信号满足稀疏分布时，绝大部分时刻都存在一个源信号从原点到
Figure BDA0002059705130000153
的一条路径。
基于
Figure BDA0002059705130000154
可知观测信号向量
Figure BDA0002059705130000155
可由基向量aij(i＝1,2)(j＝1,2,3,4,5,6)线性组合，而组合系数为各源信号
Figure BDA0002059705130000156
和aijsK(f)(K＝1,2,3,4,5,6)首尾相连组成一个封闭几何图形，所有向量aijsK(f)的长度之和即为其系数的绝对值之和
Figure BDA0002059705130000157
因此在所有的可行解中，最小值
Figure BDA0002059705130000158
的解是
Figure BDA0002059705130000159
的最短路径。故在混合矩阵估计出来后，由最大后验法，基于任给的
Figure BDA00020597051300001510
都有
Figure BDA00020597051300001511
此时源信号恢复就转为L1范数最小路径。从而估计出源信号S，实现盲源分离。
基于干涉频谱的实部和虚部均包含被测表面光程差信息，为简化，以下频率数均转换为0至1.5的范围内，故由4表面所产生的6个峰值可知：
频率f1对应载玻片前后两表面S3和S4，频率f2对应光楔前后两表面S1和S2，频率f3对应光楔前表面S1和载玻片前表面S3，频率f4对应光楔前表面S1和载玻片后表面S4，频率f5对应光楔后表面S2和载玻片前表面S3，频率f6对应光楔后表面S2和载玻片后表面S4。
其在通道1上实验：
f1＝0.470Hz,f2＝2.704Hz,f3＝5.523Hz,f4＝6.014Hz,f5＝8.241Hz,f6＝8.718Hz，
可求解出光程差：Λ1pq(x,y)＝π·f1pq(x,y)
从而S3和S4光程差Λ34＝3.14*0.470＝1.476mm
S1和S2光程差Λ12＝3.14*2.704＝8.491mm
S1和S3光程差Λ13＝3.14*5.523＝17.342mm
S1和S4光程差Λ14＝3.14*6.014＝18.884mm
S2和S3光程差Λ23＝3.14*8.241＝25.877mm
S2和S4光程差Λ24＝3.14*8.718＝27.375mm
在通道2上实验：
f1＝0.470Hz,f2＝2.704Hz,f3＝5.525Hz,f4＝6.012Hz,f5＝8.239Hz,f6＝8.719Hz，
可求解出光程差：Λ2pq(x,y)＝π·f2pq(x,y)mm
从而S3和S4光程差Λ34＝3.14*0.470＝1.476mm
S1和S2光程差Λ12＝3.14*2.704＝8.491mm
S1和S3光程差Λ13＝3.14*5.525＝17.348mm
S1和S4光程差Λ14＝3.14*6.012＝18.878mm
S2和S3光程差Λ23＝3.14*8.239＝25.870mm
S2和S4光程差Λ24＝3.14*8.719＝27.378mm
本发明迈克尔逊式干涉的原理是由激光波数扫描双通道分光路的干涉测量系统通过不可见激光经透镜准直为平行光后由分光棱镜分为两束相干光，一束光直射到光楔上，一束光被分光棱镜照射到被测件载玻片的两表面，从而发生振幅分割干涉。故本发明通过K均值聚类算法估计出混合矩阵，通过L1范数最短路径法恢复源信号。在欠定条件下，混叠矩阵未知且不可逆时，其混合矩阵的估计是基于K均值聚类分析。
相对于普通的干涉测量方法无法准确地识别出各个峰值位置，本发明利用双通道检测和采样，通过盲源分离，使得信号即使面对较为强烈的噪声，也可以检测出各个干涉信号所在位置；引入K均值聚类算法以及L1范数最小路径算法分别对欠定矩阵进行求解以实现源信号的分离；由于欠定盲源分离要求信号为稀疏信号，而时域光强信号不满足稀疏性，本发明对信号进行快时傅里叶变换FFT变换至频域，形成稀疏信号；观测信号I＝2，即为多通道欠定盲源分离提供了有利的条件；为解决欠定盲源分离，首先估计混合矩阵，再由混合矩阵恢复出源信号。
上述为本发明较佳的实施方式，但本发明的实施方式并不受上述内容的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。

Claims (2)

1.一种基于欠定盲源分离双通道光学三维干涉方法，其特征在于，包括下述步骤：
步骤一，计算机发出控制指令，并通过激光控制器和温控模块对半导体激光器波数输出进行线性调频；
其中所述线性调频公式为：
Figure FDA0002751911470000011
其中k(t)为激光波数，λc为半导体激光器工作时的波长，Δλ为激光输出波长的变化范围，T为激光输出波长对应的结束时刻；
步骤二，激光输出光经透镜L1准直为平行光后，被50:50的正方体分光棱镜分成两束光，一束光经液晶片照射到载玻片的两表面S3和S4，另一束光照射到光楔前后两表面S1和S2上，载玻片两表面S3和S4的散射光束和光楔前后两表面S1和S2的散射光束再次经过50:50的正方体分光棱镜后，形成一条返回光路，返回的光路通过透镜L2聚焦到CCD相机里的数据采集卡，并在数据采集卡上相互叠加形成干涉信号；实验中激光器发出的相干光照射到所有被测表面时，通道上表面反射光强相互叠加，形成M(M-1)/2个峰值的干涉信号；其干涉信号因液晶片所产生的光强不同，在通道中产生的幅值高度就不同，故被分为两路干涉信号I1和I2；
其中，I1的干涉信号公式为：
Figure FDA0002751911470000021
,
其中Λ1pq(x,y)＝Z1pq(x,y)τ1pq(x,y)；
上式中(x,y)为空间坐标；下标p,q分别代表被测表面SP和Sq；I1为反射光强；Λ1pq为被测表面SP和Sq之间的光程差；Z1pq为被测表面SP和Sq之间沿z方向上的位置差；τ1pq为被测表面SP和Sq之间介质的光学折射率；
Figure FDA0002751911470000023
为被测表面SP和Sq之间的初始相位差；
对上述公式I1作傅里叶变换得：
Figure FDA0002751911470000022
其中，I2的干涉信号公式为：
Figure FDA0002751911470000031
,
其中Λ2pq(x,y)＝Z2pq(x,y)τ2pq(x,y)；
上式中(x,y)为空间坐标；下标p,q分别代表被测表面SP和Sq；I2为反射光强；Λ2pq为被测表面SP和Sq之间的光程差；Z2pq为被测表面SP和Sq之间沿z方向上的位置差；τ2pq为被测表面SP和Sq之间介质的光学折射率；
Figure FDA0002751911470000035
为被测表面SP和Sq之间的初始相位差；
对上述公式I2作傅里叶变换得：
Figure FDA0002751911470000032
基于I1干涉信号能采集到的幅频特性有M(M-1)/2个峰值，每个峰值对应载玻片被测件表面的干涉信号Spq，记该峰值处的幅频和相频为f1pq和Φ1pq，则得到以下公式：
Figure FDA0002751911470000033
Figure FDA0002751911470000034
基于I2干涉信号能采集到的幅频特性有M(M-1)/2个峰值，每个峰值对应载玻片被测件表面的干涉信号Spq，记该峰值处的幅频和相频为f2pq和Φ2pq，则得到以下公式：
Figure FDA0002751911470000041
Figure FDA0002751911470000042
干涉幅频fpq和干涉相频Φpq均包含载玻片被测表面深度轮廓信息，可以通过它们解调出Λpq；
为作简述方便，省略I1(x,y,k)和I2(x,y,k)中的(x,y)，因(x,y)在整个过程中保持一致，故只对k做划分；上述k均是在时域中，信号经傅里叶变换变换到频域中用f代替，以作区分；
步骤三，由于欠定盲源分离要求信号为稀疏信号，而时域光强信号不满足稀疏性，所以对信号进行快时傅里叶变换FFT变换至频域，形成稀疏信号；
基于观测信号个数I＝2小于源信号个数S＝6，故以K均值聚类方法估计混合矩阵A；
设置聚类个数为源信号个数，即K＝S，K均值聚类方法估计混合矩阵A求解步骤如下：
(1)从M(M-1)/2个峰值中选出K个初始聚类中心，将全部峰值随机分成K类；
(2)计算每个峰值信号到各聚类中心的欧氏距离；所述的欧氏距离公式为：
Figure FDA0002751911470000051
(3)误差函数不在变化，则聚类结束；所述的误差函数为：
Figure FDA0002751911470000052
其中Ψ(cj)代表以cj为中心的列向量集合，d(xi(t),cj)为每个采样信号点xi(t)与聚类中心cj的欧氏距离，当且仅当dcj(xi(t),cj)＝min{d(xi(t),cp)，p＝1,…,k}时，函数获得最优解；
步骤四，基于上述源信号在时域上不满足稀疏性的要求，所以对源信号进行快时傅里叶变换FFT，使源信号在变换后的频域上满足稀疏性，从而在变换后的频域上进行盲源分离，变换后的I1(k)和I2(k)满足频域：
Figure FDA0002751911470000053
其中
Figure FDA0002751911470000054
为频域中I1(f)和I2(f)的观测实验结果；
Figure FDA0002751911470000055
为I1(f)、I2(f)干涉光强中幅频的峰值，即S个数为6；故可以从两路观测信号中分离出六路源信号，则观测信号个数小于源信号个数；A为未知的2×6混合矩阵；A和S(f)均未知且混合矩阵A不可逆；
故所述观测信号向量
Figure FDA0002751911470000063
在幅频中展开可得矩阵：
Figure FDA0002751911470000061
根据稀疏理论，所以在某一频率中只有一个信号取值幅频相对较大；假在f1＝1,...,1024，只有S1(f1)的取值较大，则根据上述公式有：
(1)I1(f1)＝a11S1(f1)，I2(f2)＝a21S1(f1)，将前边两式相除得I2(f2)/I1(f1)＝a21/a11，由此可得，在观测信号的幅频图上，S1(f1)取值的点均聚集在斜率为a21/a11的直线周围；
(2)同理当S2(fi)，S3(fi)，S4(fi)，S5(fi)，S6(fi)幅值较大时点分别均聚集在斜率为a22/a12，a23/a13，a24/a14，a25/a15，a26/a16的直线周围；
所以根据这6条直线的斜率tanθ的值求得sinθ和cosθ，对观测信号数据点进行聚类划分就能得到混合矩阵的各列向量，即
Figure FDA0002751911470000062
为求解的矩阵A中的一列向量；
步骤五，基于K均值聚类法估计出混合矩阵A后，通过最小化L1范数最短路径法来实现源信号的恢复；操作原理为通过对观测信号进行分解，寻找到最靠近观测信号向量的线性组合来作为对源信号的估计，提取部分源的基矢量来最小干扰地恢复出源信号；该方法是先根据观测信号所在方向，计算其与混合矩阵中列向量所代表方向的差，设定一定的阈值，从中挑选出差值较小的多个代表性方向作为观测信号分解的潜在方向；
根据对源信号的先验性假设，当源信号满足稀疏分布时，绝大部分时刻都存在一个源信号从原点到
Figure FDA0002751911470000071
的一条路径；
根据上述欠定盲源分离模型的矩阵形式，可知观测信号向量
Figure FDA0002751911470000072
其中
Figure FDA0002751911470000073
包括I1(f)和I2(f)，可由基向量aij(i＝1,2)(j＝1,2,3,4,5,6)线性组合，而组合系数为各源信号S1(f)，S2(f)，...，S6(f)；所述
Figure FDA0002751911470000074
和aijsk(f)(k＝1,2,3,4,5,6)首尾相连组成一个封闭几何图形，所有向量aijsk(f)的长度之和即为其系数的绝对值之和
Figure FDA0002751911470000075
因此在所有的可行解中，最小值
Figure FDA0002751911470000076
的解是
Figure FDA0002751911470000077
的最短路径；故在混合矩阵估计出来后，由最大后验法，基于任给的
Figure FDA0002751911470000078
都有
Figure FDA0002751911470000081
此时，源信号恢复就转为L1范数最小路径，从而估计出源信号S；
基于最小二乘的参数估计思路是通过使测量误差的平方和最小，寻找待匹配函数的最优参数值；
步骤六，由于I1、I2的方程
Figure FDA0002751911470000082
中干涉频谱的实部和虚部均包含被测表面光程差信息，因此分别对实部和虚部作最小二乘的方法，故在频域内对干涉测量频谱构建误差方程E(x,y)如下：
Figure FDA0002751911470000083
其中L为傅里叶变换频率点数；
Figure FDA0002751911470000084
为干涉测量频谱，可通过CCD相机拍摄的干涉信号作傅里叶变换得到；
干涉相位和幅值可求解为：
Figure FDA0002751911470000085
Figure FDA0002751911470000086
此时表面SP，Sq之间的干涉信号Spq的波数域光强表达式为：
Figure FDA0002751911470000091
即为干涉信号复数域线性最小二乘算法；复数域线性最小二乘算法仅通过CCD相机采集的干涉信号便可自动分离出各个表面的干涉信号，有助于载玻片干涉信号的提取；另一方面，干涉信号成功地盲分离可以进一步提高激光波数扫描干涉检测的深度轮廓分辨率和减少相位波动误差。
2.根据权利要求1所述的基于欠定盲源分离双通道光学三维干涉方法，其特征在于，步骤三中所述稀疏信号是指信号个数在1024个频率中的幅频取值接近零，同时又有六个明显的波峰幅频较大；双通道峰值除了六个明显的波峰外，绝大多数地方平滑下降，知源信号之间是相互统计独立的，满足稀疏特性的要求，所以绝大部分时刻最多只有一个源信号取值占优。
CN201910400871.5A 2019-05-15 2019-05-15 一种基于欠定盲源分离双通道光学三维干涉方法及系统 Expired - Fee Related CN110260812B (zh)

Priority Applications (2)

Application Number Priority Date Filing Date Title
CN201910400871.5A CN110260812B (zh) 2019-05-15 2019-05-15 一种基于欠定盲源分离双通道光学三维干涉方法及系统
US16/875,017 US11060849B2 (en) 2019-05-15 2020-05-15 Dual-channel optical three-dimensional interference method and system based on underdetermined blind source separation

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910400871.5A CN110260812B (zh) 2019-05-15 2019-05-15 一种基于欠定盲源分离双通道光学三维干涉方法及系统

Publications (2)

Publication Number Publication Date
CN110260812A CN110260812A (zh) 2019-09-20
CN110260812B true CN110260812B (zh) 2021-02-09

Family

ID=67913154

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910400871.5A Expired - Fee Related CN110260812B (zh) 2019-05-15 2019-05-15 一种基于欠定盲源分离双通道光学三维干涉方法及系统

Country Status (2)

Country Link
US (1) US11060849B2 (zh)
CN (1) CN110260812B (zh)

Families Citing this family (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111129916B (zh) * 2019-12-29 2021-01-15 中国科学院西安光学精密机械研究所 一种具有调试补偿功能的时间整形系统
CN111829954B (zh) * 2020-09-09 2023-07-25 广东工业大学 一种提高全场扫频光学相干层析测量量程的系统及方法
CN112288762B (zh) * 2020-10-15 2023-05-09 西北工业大学 一种有限角度ct扫描的离散迭代重建方法
CN112432590B (zh) * 2020-12-14 2022-07-05 西安邮电大学 一种基于约束欠定方程的三波长数字全息成像光路及方法
CN113009417B (zh) * 2021-02-05 2022-09-20 中国人民解放军国防科技大学 利用声场干涉特性的海底声学阵列阵形估计方法
US20220364848A1 (en) * 2021-05-13 2022-11-17 Industrial Technology Research Institute Depth measurement apparatus and depth measurement method
CN113504182A (zh) * 2021-05-26 2021-10-15 杭州久益机械股份有限公司 基于激光声表面波的物体表面裂纹在线检测方法
CN114965367B (zh) * 2022-06-14 2024-06-25 广东工业大学 一种用于光学层析测量的混叠正弦波信号分离方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101727908A (zh) * 2009-11-24 2010-06-09 哈尔滨工业大学 基于混合信号局部峰值方差检测的盲源分离方法
US7853031B2 (en) * 2005-07-11 2010-12-14 Siemens Audiologische Technik Gmbh Hearing apparatus and a method for own-voice detection
CN105354594A (zh) * 2015-10-30 2016-02-24 哈尔滨工程大学 一种针对欠定盲源分离的混合矩阵估计方法
CN105890538A (zh) * 2014-12-30 2016-08-24 广东工业大学 三表面干涉式高精度曲面轮廓测量系统及方法
CN108009584A (zh) * 2017-12-01 2018-05-08 西安交通大学 基于单源点检测的欠定盲源分离方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7853031B2 (en) * 2005-07-11 2010-12-14 Siemens Audiologische Technik Gmbh Hearing apparatus and a method for own-voice detection
CN101727908A (zh) * 2009-11-24 2010-06-09 哈尔滨工业大学 基于混合信号局部峰值方差检测的盲源分离方法
CN105890538A (zh) * 2014-12-30 2016-08-24 广东工业大学 三表面干涉式高精度曲面轮廓测量系统及方法
CN105354594A (zh) * 2015-10-30 2016-02-24 哈尔滨工程大学 一种针对欠定盲源分离的混合矩阵估计方法
CN108009584A (zh) * 2017-12-01 2018-05-08 西安交通大学 基于单源点检测的欠定盲源分离方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
基于频谱分离最小二乘法的三维轮廓测量系统;黄启权;《中国学位论文全文数据库》;20181219;第11-12、27-30、36页 *

Also Published As

Publication number Publication date
CN110260812A (zh) 2019-09-20
US11060849B2 (en) 2021-07-13
US20200378756A1 (en) 2020-12-03

Similar Documents

Publication Publication Date Title
CN110260812B (zh) 一种基于欠定盲源分离双通道光学三维干涉方法及系统
CN105066908B (zh) 一种基于多波长和多偏振态的数字全息三维形貌检测装置
CN106767489B (zh) 数字散斑干涉面内微小动态形变测量系统及测量方法
CN105606217B (zh) 一种图像、光谱、偏振态一体化获取装置及方法
US7847954B2 (en) Measuring the shape and thickness variation of a wafer with high slopes
CN102889853B (zh) 分光同步移相共光路干涉显微检测装置及检测方法
JP6256879B2 (ja) 偏光感受型光画像計測システム及び該システムに搭載されたプログラム
CN109000781B (zh) 一种结构微振动线域测量装置及方法
US11131539B2 (en) Multimodal image data acquisition system and method
CN110361099B (zh) 一种谱域低相干光干涉光程差解调方法
CN113418469B (zh) 光谱共焦扫描共光路数字全息测量系统及测量方法
CN111538027A (zh) 一种用于高分辨率测量的激光测距装置及方法
US20220113674A1 (en) Differential holography
Wu et al. Universal optical setup for phase-shifting and spatial-carrier digital speckle pattern interferometry
CN106091974B (zh) 一种物体形变测量仪器、方法和设备
CN108982510A (zh) 利用90°光学混波器数字剪切散斑动态检测系统及方法
US11892801B2 (en) Systems and methods for simultaneous multi-channel off-axis holography
CN109341519A (zh) 用于确定结构中的兴趣区域的参数的方法和系统
CN110132846A (zh) 基于马赫曾德光路的多方向剪切散斑干涉系统及测量方法
CN107894204B (zh) 干涉仪及其成像方法
KR20210157969A (ko) 복합막 단층 측정을 위한 영상 시스템 및 그 방법
CN208833364U (zh) 一种结构微振动线域测量装置
CN105716521A (zh) 增大频域低相干光干涉测量范围的装置及其方法
KR20170037800A (ko) 투명 매질의 미세 결함 검출방법 및 시스템
CN114965367B (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
CF01 Termination of patent right due to non-payment of annual fee
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20210209

Termination date: 20210515