CN112097678B - 基于频率盲估计的多表面面形测量方法 - Google Patents

Abstract

本发明公开了一种基于频率盲估计的多表面面形测量方法，使用波长移相干涉仪对干涉图进行采集，通过对多表面干涉图中的各表面干涉频率的盲估计，对所对应的信号频率特性进行精确求解和匹配。将求解得到的各信号的干涉频率带入到所构建的新型最小二乘求解方程求取被测件各干涉信号的幅值和初始相位，完成透明平行平板的多表面干涉测量。所设计的测量方法成本较低，避免了被测件平均厚度的测量，并且能够对被测件各表面的对比度进行精确求解。

Description

基于频率盲估计的多表面面形测量方法

技术领域

本发明涉及一种基于干涉频率盲估计的多表面面形测量方法，尤其是通过盲分离方法求取各个信号的干涉频率，进而通过新型最小二乘方法求取初始相位进行面形重建的多表面测量方法。本发明通过各表面信号频率和相位的精确求解，从而实现干涉信号的盲分离和重建。

背景技术

高精度多表面平行平板在光学***的构建过程中具有重要作用。其中多表面透镜的精确表面测量和重建，对于光学***成像和传输质量的提高和评估具有实用意义与研究价值，是领域内人员研究的热点之一。

传统的透镜测量多是使用硬件干涉仪或者表面接触的逐点扫描的方法进行测量。但是这些测量方法的局限性在于：对于硬件干涉仪，由于使用压电陶瓷驱动改变局部的干涉腔长，从而进行移相操作，这种情况下各个信号的移相频率是相同的，无法进行分离，因此在实际应用时多使用逐次测量和表面涂抹消光材料的方法进行测量。但是这种方法会对被测件的高精度表面造成损伤，以及在涂抹、清洗和被测件的旋转过程中纳入多种误差和残余干涉信号，无法实现被测件的一次性无损测量。对于逐点扫描的方法而言，测量成本较高，并且容易对高精度表面造成损伤，无法达到非接触测量的目的。

近年来，使用波长移相干涉仪及其相关算法对多表面透镜进行测量的方法得到了发展。但是现存的算法多是基于对干涉频率的精确获取的前提下进行测量，例如通过被测件的腔长值与被测件的光学厚度的比值的精确计算，去获取各个信号的干涉相位。但是该种信号获取的方法的局限性在于，被测件的腔长值难以被精确地测量，因为干涉仪的参考镜被夹持在夹具上，精确的获取腔长值，即参考镜表面到达被测件前表面的距离值，这在实际测量中是不易实现的，甚至在获取上述测量值的过程中容易对参考镜以及被测件的表面造成划伤，不能满足高精度测量的要求。与此同时，传统的干涉测量方法是对于干涉图上的所有点全部纳入计算，即：以单个点随着时间或者干涉图帧序数的变化的信号作为整张干涉图的代表，进行对全部的移相后的一组干涉图进行处理。这种方法的缺点在于：如果所选择的点噪声较大或者表面存在缺陷等信息，则会致使整组干涉图均不会得到精确结果；反之亦然，若所选取的代表性的点不足以代表整个像素点分布的信息变化，仍然也不会取得精确计算结果。此外，多表面透镜的测量难点在于：在对被测件进行干涉图的采集过程中，由于各个平面距离参考镜的距离是不同的，因此在使用波长移相干涉仪进行移相操作时，对各干涉信号的的相位改变量是不同的，即各干涉信号之间具有不同的移相频率，这是进行信号分离的基础之一。因此在干涉图中，各个信号均会混叠在采集的数据中，使得信号混叠的干涉图无法被直接地进行求解。

发明内容

为了解决现有技术问题，本发明的目的在于克服已有技术存在的不足，提供一种基于频率盲估计的多表面面形测量方法，特别是通过被测件各表面的干涉频率的盲估计和所构建的新型最小二乘方法求解初始相位分布的多表面透明透镜的测量方法。本发明所针对的对象是高精度透明多表面平行平板，通过对移相后采集的到的混叠的干涉图通过本发明所提出的测量方法进行处理，能一次性得到被测件前表面、后表面的面型和厚度变化。

为达到上述发明创造目的，本发明采用如下发明构思：

本发明的原理如下：

对于所采集的移相后的干涉图(每两帧干涉图之间具有相位差)而言，因光的传播和振动特性，其在时域上可将图上某一点的连续移相看作是余弦信号随着时间或干涉图帧序数的变化，因此包含前后表面和厚度变化的叠加干涉图可当作是三组叠加的余弦信号，并且三个信号之间具有不同的信号频率，即干涉频率，这也是进行多表面测量的基础。通过信号盲估计的方法求解各个信号的精确的相移频率，再带入到本发明所提出的构建的新型最小二乘求解方法中，进行各个信号的初始相位的求解，即可对被测件的面形进行恢复和重建，达到一次性多表面被测件的干涉测量的目的。本发明实现多表面被测件的非接触式无损测量，并且一次性得到多表面干涉信号所包含的表面形貌分布，使用盲信号频率精确估计进行干涉频率求解，进而通过一种新型构造的最小二乘法求取各个被测表面面形的测量。

根据上述发明构思，本发明采用如下技术方案：

一种基于频率盲估计的多表面面形测量方法，通过对采集得到的干涉图中的各点干涉频率的盲估计和幅值、相位的求解，完成具有多表面的被测件的一次性无接触测量；使用波长移相干涉仪对被测件进行干涉图的采集，所述基于频率盲估计的多表面面形测量方法包括如下步骤：

a.对叠加的余弦信号的光强变化条纹图进行基于各信号干涉频率独立分布的盲频率估计，从而确定各个干涉信号的精确频率分布；

b.通过所构建的新型最小二乘初始相位和幅值求解方法进行计算，从而恢复各个干涉信号中所包含的面形信息，对透明的平行平板进行多表面分离和重建，并对被测件的腔长值和平均厚度进行求取，从而完成被测件的多表面干涉测量。

优选地，在进各个信号的干涉频率的求解时，考虑由前后表面和厚度变化的叠加的干涉图，根据激光传播特性，采集的到的条纹图表示为三组余弦信号相叠加的形式，将

表示为由K个余弦信号及噪声e(n)叠加组成的信号，写为以下表达式：

式中，实采数据的序数为n＝1，...，N，{A_k，ω_k：k＝1，...K}是未知实量，θ_k：k＝1，...K是在-π～π之间的独立随机变量，具有均匀分布特性，并且对应各个信号的初始相位；ω_k是各个信号的频率；A_k是各信号的对比度；所述盲频率估计是通过余弦检测从给定的N个实采数据

中精确计算各未知数；假设e(n)的均值为零、方差为σ²的噪声分量；采用消除噪声的方法为：对每一个像素点计算各帧干涉图下的平均值，将采集到的干涉中各点减去此值，消除背景光强和噪声；移相值设置为π/4，为了保证测量精度，总采集帧数选择为200帧；

对因此将无噪声谐波y(n)及其组成表示为：

改写为下式：

式中，

其中e是自然底数，j是虚数单位；因此当频率分量ω_k求出后，以上N个方程组是b_k的线性方程，因此求出各个b_k是具有可行性的，并由实采序列

的各组成信号的频率进行计算和分析；进而，任意频率的K个无噪声谐波的实数余弦信号可用线性预测模型中的预测系数a_k进行表达，并且其表达式具有对称性；由实采序列可求出预测系数a_k的估计值，因此对预测模型表达式做Z变换得到下式：

式中，Ψ(z)是无噪声谐波表达式通过线性预测模型进行表示的形式，进而使用Z变换后得到的原式的表达式，其2K个根均为z平面上单位圆的共轭复根；可由z的表达式得到余弦信号频率ω(n)可由下式求得：

优选地，在进行各个信号的幅值和相位求解时，考虑由3组余弦信号组成的混叠的干涉信号，此时引入一个N×N维的正交矩阵Q_H，为6个初等反射变换阵的乘积，是总的初等反射变换矩阵，使得：

因为||Q_H||＝I，所以有如下公式：

||Ax-b||₂＝||Q_HAx-Q_Hb||₂＝||Tx-Q_Hb||₂

令

且取

的前6个分量，有如下公式：

式中，m＝6；选择采集200帧干涉图，即N＝200；上式的右边的第2项不含有待估计参量，为使上式达到最小值，取

由此可得待估参量的最小二乘估计式，用符号x表示为

然后进一步得到各个信号的幅值和初始相位的估计值，其中角标为1，2，3分别对应后表面、前表面和厚度的信息，参见如下公式：

求解得到的初始相位；在得到各表面的初始相位以后，因为所得到的是包裹结果，因此需要进行一次解包裹操作；解包裹以后，得到前表面、后表面、厚度变化三个信号所对应的初始相位为θ_f、θ_r、θ_f-r，通过如下操作得到三个信号的对应的面形分布公式：

其中n1代表被测件折射率，h代表干涉仪参考镜的后表面到达被测件外表面的距离，W_f、W_r、W_f-r分别代表被测件的前表面、后表面和厚度变化信号的面形分布。

优选地，在进行被测件的腔长值和平均厚度的求解时，通过求取的各信号的频率分布，通过光程差的计算，反算出被测件的腔长值以及其平均厚度，参见如下公式：

其中Δλ为单步移相的波长移相干涉仪中激光器的调谐量，人为进行设置；λ₀为起始调谐波长；n₁是被测件的折射率，此值与材料有关，是已知量；后表面的光程h+n₁×d是与腔长h和被测件光学厚度n₁×d的相加值，在计算腔长和被测件的平均厚度d时，将前表面信号的干涉频率带入到上式计算腔长值，而后将后表面信号计算的腔长与被测件光学厚度的和的值减去通过前表面干涉频率计算的腔长值，剩下的部分与厚度变化信号干涉频率所估算的光学厚度值进行平均值的计算，从而进一步减少估算误差，得到更为精确的被测件的光学厚度，从而与折射率做差求得被测件的平均厚度d。

优选地，根据在多表面干涉过程中的光程差对于干涉频率的影响，当求解出频率的数据组以后，在进行频率匹配时，应当把该数据组进行降序排列，此时各表面信号所对应的频率分量为后表面、前表面、和厚度变化信号。

优选地，进行初始相位和幅值的求解，将所构造的最小二乘系数矩阵经过N次初等反射变换后得到新元素，N为系数矩阵A的列数；经过若干次初等反射变换以后，所得到的新的矩阵的范数与进行变换之前是不改变的，但是形式已经化为部分三角矩阵，在新组成的矩阵中，非零元素的列数不变，但是得到的经过初等反射变换以后的矩阵的非零行数M＜m。

优选地，所述新型最小二乘初始相位和幅值求解方法的目标在于寻求a_k的最佳值使得下式中的扰动矩阵E和扰动项r的总误差

最小，这里

为矩阵的F-范数：

(A+E)a＝(b+r)

上式中扰动矩阵E和扰动项r的总体误差

最小，A，b为数据矩阵；将上式表示为增广形式：

并且：

C＝[A，b]∈C^{(N-2K)×(K+1)}

L＝[E，r]

因在实际测量过程中不可避免地存在噪声，因此C是满秩的；记C₁＝C+L为扰动矩阵的增广形式；d＝[a，-1]^T为未知向量；上述增广矩阵又可表示为：

(C+L)d＝C₁d＝0

根据上式，使用所构造的新型最小二乘方法求解未知量d的最佳值的问题又表述为：计算或搜寻最小F-范数矩阵L使C₁＝C+L为非满秩形式；与此同时，应当使C₁的零空间向量d的最后一个元素非0，从而通过d的最后一个元素按照

的求解方法归一化地求解a的值；

在矩阵

中，假设(N-2K)≥(K+1)，因此矩阵C₁的秩最大为(K+1)，将矩阵C₁最后一列全部化为0，这样矩阵C₁便为非满秩矩阵，从而解出增广形式矩阵的最佳值；

上述最佳值的求解问题由矩阵C的奇异值分解进行计算，将C写成奇异值分解的内积表示：

上述表示形式有最小F-范数，可使C降秩的矩阵L便是上式中的最后一项的相反数

C的零空间向量d＝v_K+1与

都是正交的，k＝1，...，K；因此，将v_K+1的最后一个元素对v_k进行归一化处理，即可求出上述方程组的最小二乘的最佳值a。

优选地，采用所构建的新型最小二乘方法进行叠加的余弦信号中的盲频率估计的过程包括以下步骤：

第一步：设置移相值为π/4，使用波长移相干涉仪采集干涉图200帧；

第二步：计算每个像素点的平均值，并且每个像素点减去该值以消除噪声和背景光强，得到纯谐波形式的混叠的干涉光强；

第三步：将纯谐波形式的数据矩阵A，b构造其增广表达式C＝[A，b]∈C^{(N-2K)×(K+1)}，其元素C(i，j)为：

第四步：对上述C矩阵进行奇异值分解C＝U∑V^T，并对第(K+1)个向量的最后一个元素

进行归一化操作得到a_k；

第五步：按式所构造的ψ(z)求其根z_k，最后由式

求出各信号的频率，带入到通过对采集得到的干涉图中的各点干涉频率的盲估计和幅值、相位的求解过程中，根据在多表面干涉过程中的光程差对于干涉频率的影响，当求解出频率的数据组以后，在进行频率匹配时，将该数据组进行降序排列，此时各表面信号所对应的频率分量为：后表面、前表面、厚度变化信号。

优选地，所述被测件的夹具可以选用自定心透镜安装架4SCML-8。

优选地，所述干涉图的采集通过选用INF300-LP干涉仪进行采集。

因此，在干涉腔长、被测件厚度和折射率未知的情况下，根据当前的所提出的信号盲频率估计方法和相位求解方法即可实现多表面透明被测件的一次性测量。

本发明与现有技术相比较，具有如下显而易见的突出实质性特点和显著优点：

1.本发明无需被测件的腔长、平均厚度值和折射率的精确获取，节省测量成本，直接通过信号频率的盲估计达到干涉频率求解的目的；

2.本发明通过干涉频率的精确求解，可使得干涉信号之间分离的较为纯净，从而为各个信号的初始相位的求解提供基础；

3.本发明通过所构建的新型最小二乘法求取被测件各个表面的初始相位，可实现一次性的多表面被测件的干涉初始相位的求解，从而获取被测件的前表面、后表面的面形和厚度变化；可达到非接触式无损测量的目的；

4.能够对被测目标的各表面的对比度进行求解，获得精度较高的对比度分布图，可对被测件各表面缺陷点等部分进行快速定位等操作提供基础；

5.本发明所使用的测量方法是对一组干涉图中的各像素点进行逐点计算，各点计算均独立进行，即使在噪声较大的像素点或者缺陷点处的计算不够精确，仍然能够满足其余的点能够保证求解精度，与传统算法中使用一个点进行整组干涉图的频率估计的方法相比，求解可靠性和精度大大提高，即使干涉信号之间频率相近时，仍然能求解出精确的频率分布，求解精度远高于传统快速傅里叶变换方法；

6.本发明方法不需要过多的先验知识便可以进行计算，不仅如此，只需要将被测件的折射率输入测量算法，则本发明算法还能实现对于腔长值和被测件的平均厚度的测量，大大降低了对于多表面干涉测量的测量需求和测量难度。

附图说明

图1是本发明方法的流程图。

图2是本发明方法采集的初始干涉图。

图3是通过本发明方法求得的各表面初始包裹相位。

图4是本发明方法的面形求解误差。

图5是本发明方法的对比度求解误差。

图6是本发明方法求解得到的对比度分布。

图7是本发明方法求解得到的各表面面形分布。

具体实施方式

以下结合具体的实施例子对上述方案做进一步说明，本发明的优选实施例详述如下：

实施例一：

在本实施例中，参见图1，一种基于频率盲估计的多表面面形测量方法，通过对采集得到的干涉图中的各点干涉频率的盲估计和幅值、相位的求解，完成具有多表面的被测件的一次性无接触测量；使用波长移相干涉仪对被测件进行干涉图的采集，所述基于频率盲估计的多表面面形测量方法包括如下步骤：

a.对叠加的余弦信号的光强变化条纹图进行基于各信号干涉频率独立分布的盲频率估计，从而确定各个干涉信号的精确频率分布；

b.通过所构建的新型最小二乘初始相位和幅值求解方法进行计算，从而恢复各个干涉信号中所包含的面形信息，对透明的平行平板进行多表面分离和重建，并对被测件的腔长值和平均厚度进行求取，从而完成被测件的多表面干涉测量。

本实施例基于频率盲估计的多表面面形测量方法，特别是通过被测件各表面的干涉频率的盲估计和所构建的新型最小二乘方法求解初始相位分布的多表面透明透镜的测量方法。本发明所针对的对象是高精度透明多表面平行平板，通过对移相后采集的到的混叠的干涉图通过本发明测量方法进行处理，能一次性得到被测件前表面、后表面的面型和厚度变化。

实施例二：

本实施例与实施例一基本相同，特别之处在于：

在本实施例中，参见图1，在进各个信号的干涉频率的求解时，考虑由前后表面和厚度变化的叠加的干涉图，根据激光传播特性，采集的到的条纹图表示为三组余弦信号相叠加的形式，将

表示为由K个余弦信号及噪声e(n)叠加组成的信号，写为以下表达式：

式中，实采数据的序数为n＝1，...，N，{A_k，ω_k：k＝1，...K}是未知实量，θ_k中的k＝1，...，K是在-π～π之间的独立随机变量，具有均匀分布特性，并且对应各个信号的初始相位；ω_k是各个信号的频率；A_k是各信号的对比度；所述盲频率估计是通过余弦检测从给定的N个实采数据

中精确计算各未知数；假设e(n)的均值为零、方差为σ²的噪声分量；采用消除噪声的方法为：对每一个像素点计算各帧干涉图下的平均值，将采集到的干涉中各点减去此值，消除背景光强和噪声；移相值设置为π/4，为了保证测量精度，总采集帧数选择为200帧；

对因此将无噪声谐波y(n)及其组成表示为：

式中：n＝1，...，k；

改写为下式：

式中，

k＝1，...，K其中e是自然底数，j是虚数单位；因此当频率分量ω_k求出后，以上N个方程组是b_k的线性方程，因此求出各个b_k是具有可行性的，并由实采序列

的各组成信号的频率进行计算和分析；进而，任意频率的K个无噪声谐波的实数余弦信号可用线性预测模型中的预测系数a_k进行表达，并且其表达式具有对称性；由实采序列可求出预测系数a_k的估计值，因此对预测模型表达式做Z变换得到下式：

式中，Ψ(z)是无噪声谐波表达式通过线性预测模型进行表示的形式，进而使用Z变换后得到的原式的表达式，其2K个根均为z平面上单位圆的共轭复根；可由z的表达式得到余弦信号频率ω(n)可由下式求得：

在本实施例中，在进行各个信号的幅值和相位求解时，考虑由3组余弦信号组成的混叠的干涉信号，此时引入一个N×N维的正交矩阵Q_H，为6个初等反射变换阵的乘积，是总的初等反射变换矩阵，使得：

因为||Q_H||＝I，所以有如下公式：

||Ax-b||₂＝||Q_HAx-Q_Hb||₂＝||Tx-Q_Hb||₂

令

且取

的前6个分量，有如下公式：

式中，m＝6，N＞m；选择采集200帧干涉图，即n＝200；上式的右边的第2项不含有待估计参量x，为使上式达到最小值，取

由此可得x的最小二乘估计式为

然后进一步得到各个信号的幅值和初始相位的估计值，其中角标为1，2，3分别对应后表面、前表面和厚度的信息，参见如下公式：

求解得到的初始相位；在得到各表面的初始相位以后，因为所得到的是包裹结果，因此需要进行一次解包裹操作；解包裹以后，得到前表面、后表面、厚度变化三个信号所对应的初始相位为θ_f、θ_r、θ_f-r，通过如下操作得到三个信号的对应的面形分布公式：

其中n1代表被测件折射率，h代表干涉仪参考镜的后表面到达被测件外表面的距离，W_f、W_r、W_f-r分别代表被测件的前表面、后表面和厚度变化信号的面形分布。

在本实施例中，在进行被测件的腔长值和平均厚度的求解时，通过求取的各信号的频率分布，通过光程差的计算，反算出被测件的腔长值以及其平均厚度，参见如下公式：

其中Δλ为单步移相的波长移相干涉仪中激光器的调谐量，人为进行设置；λ₀为起始调谐波长；n₁是被测件的折射率，此值与材料有关，是已知量；后表面的光程h+n₁×d是与腔长h和被测件光学厚度n₁×d的相加值，在计算腔长和被测件的平均厚度d时，将前表面信号的干涉频率带入到上式计算腔长值，而后将后表面信号计算的腔长与被测件光学厚度的和的值减去通过前表面干涉频率计算的腔长值，剩下的部分与厚度变化信号干涉频率所估算的光学厚度值进行平均值的计算，从而进一步减少估算误差，得到更为精确的被测件的光学厚度，从而与折射率做差求得被测件的平均厚度d。

在本实施例中，根据在多表面干涉过程中的光程差对于干涉频率的影响，当求解出频率的数据组以后，在进行频率匹配时，应当把该数据组进行降序排列，此时各表面信号所对应的频率分量为后表面、前表面、和厚度变化信号。

在本实施例中，进行初始相位和幅值的求解，将所构造的最小二乘系数矩阵经过N次初等反射变换后得到新元素，n为系数矩阵A的列数；经过若干次初等反射变换以后，所得到的新的矩阵的范数与进行变换之前是不改变的，但是形式已经化为部分三角矩阵，在新组成的矩阵中，非零元素的列数不变，但是得到的经过初等反射变换以后的矩阵的非零行数M＜m。

在本实施例中，所述新型最小二乘初始相位和幅值求解方法的目标在于寻求a_k的最佳值使得下式中的扰动矩阵E和扰动项r的总误差

最小，这里

为矩阵的F-范数：

(A+E)a＝(b+r)

上式中扰动矩阵E和扰动项r的总体误差

最小，这里

为矩阵的Fresenius范数，A，b为数据矩阵；将上式表示为增广形式：

并且：

C＝[A，b]∈C^{(N-2K)×(K+1)}

L＝[E，r]

因在实际测量过程中不可避免地存在噪声，因此C是满秩的；记C₁＝C+L为扰动矩阵的增广形式；d＝[a，-1]^T为未知向量；上述增广矩阵又可表示为：

(C+L)d＝C₁d＝0

根据上式，使用所构造的新型最小二乘方法求解未知量d的最佳值的问题又表述为：计算或搜寻最小F-范数矩阵L使C₁＝C+L为非满秩形式；与此同时，应当使C₁的零空间向量d的最后一个元素非0，从而通过d的最后一个元素按照

的求解方法归一化地求解a的值；

在矩阵

中，假设(N-2K)≥(K+1)，因此矩阵C₁的秩最大为(K+1)，将矩阵C₁最后一列全部化为0，这样矩阵C₁便为非满秩矩阵，从而解出增广形式矩阵的最佳值；

上述最佳值的求解问题由矩阵C的奇异值分解进行计算，将C写成奇异值分解的内积表示：

上述表示形式有最小F-范数，可使C降秩的矩阵L便是上式中的最后一项的相反数

C的零空间向量d＝v_K+1与

都是正交的，k＝1，...，K；因此，将v_K+1的最后一个元素对v_k进行归一化处理，即可求出上述方程组的最小二乘的最佳值a。

在本实施例中，采用所构建的新型最小二乘方法进行叠加的余弦信号中的盲频率估计的过程包括以下步骤：

第一步：设置移相值为π/4，使用波长移相干涉仪采集干涉图200帧；

第二步：计算每个像素点的平均值，并且每个像素点减去该值以消除噪声和背景光强，得到纯谐波形式的混叠的干涉光强；

第三步：将纯谐波形式的数据矩阵A，b构造其增广表达式C＝[A，b]∈C^{(N-2K)×(K+1)}，其元素C(i，j)为：

第四步：对上述C矩阵进行奇异值分解C＝U∑V^T，并对第(K+1)个向量的最后一个元素

进行归一化操作得到a_k；

第五步：按式所构造的ψ(z)求其根z_k，最后由式

求出各信号的频率，带入到通过对采集得到的干涉图中的各点干涉频率的盲估计和幅值、相位的求解过程中，根据在多表面干涉过程中的光程差对于干涉频率的影响，当求解出频率的数据组以后，在进行频率匹配时，将该数据组进行降序排列，此时各表面信号所对应的频率分量为：后表面、前表面、厚度变化信号。

本实施例方法使用波长移相干涉仪对干涉图进行采集，通过对多表面干涉图中的各表面干涉频率的盲估计，对所对应的信号频率特性进行精确求解和匹配。将求解得到的各信号的干涉频率带入到所构建的新型最小二乘求解方程求取被测件各干涉信号的幅值和初始相位，完成透明平行平板的多表面干涉测量。所设计的测量方法成本较低，避免了被测件平均厚度的测量，并且能够对被测件各表面的对比度进行精确求解。

实施例三：

本实施例与实施例二基本相同，特别之处在于：

在本实施例中，参见图1-7，考虑到在波长移相时，不同的表面对应的干涉信号移相频率不同，可从时域信号频率考虑，所采集到的干涉条纹图相当于多组余弦信号叠加。因此基于不同的移相频率可以对多组信号进行分离，从而提取出不同的面形信息。这样能消除被测件放置时无法满足腔长的影响。本实施例所提出的方法主要包括两个部分：

(1)对叠加的余弦信号的光强变化条纹图进行基于各信号干涉频率独立分布的盲频率估计，从而确定各个干涉信号的精确频率分布；

(2)通过所构建的新型最小二乘初始相位和幅值求解方法进行计算，从而恢复各个干涉信号中所包含的面形信息，对透明的平行平板进行多表面分离和重建，以及对被测件的腔长值和平均厚度进行求取。所述的腔长在数值上等于被测件的前表面到达干涉仪参考镜之间的实际距离，在物理意义上，因为空气的折射率近似为1，因此也等于前表面信号的光程，等价于前表面的反射信号与激光到达参考镜后表面的反射信号之间的光程差。所述的平均厚度也就是被测件的物理厚度，数值上等于前表面各点到达后表面各点的距离的平均值。总体而言，本实施例所采用的测量方法及流程图如图1所示。

若考虑由前后表面和厚度变化的叠加干涉图，如图2所示，根据激光传播特性，则采集的到的条纹图可表示为三组余弦信号相叠加的形式，因此可以将

表示为由K个频率不同的余弦信号及噪声e(n)叠加组成的信号，写为以下表达式：

式中，实采数据的序数为n＝1，...，N，{A_k，ω_k：k＝1，...K}是未知实量，θ_k：k＝1，...K，k代表不同的谐波频次，对于本实施例所针对的透明平行平板被测件而言，最大的谐波频次值为3。是在-π～π之间的独立随机变量，具有均匀分布特性，并且对应各个信号的初始相位。ω_k是各个信号的频率。A_k是各信号的对比度。本实施例的盲频率估计是通过余弦检测从给定的N个实采数据

中精确计算各未知数。这里假设e(n)的均值为零、方差为σ²的噪声分量。本实施例所选择的消除噪声的方法为对每一个像素点计算各帧干涉图下的平均值，将采集到的干涉中各点减去此值，即可消除背景光强和噪声。移相值设置为π/4，为了保证测量精度，总采集帧数选择为200帧，即N＝200。

对因此将无噪声谐波y(n)及其组成表示为：

改写为下式：

式中，

其中e是自然底数，j是虚数单位。因此当频率分量ω_k求出后，以上N个方程组是b_k的线性方程，因此求出各个b_k便是具有可行性的。所以此部分对由实采序列

的各组成信号的频率进行计算和分析。进而，任意频率的K个无噪声谐波的实数余弦信号可用线性预测模型中的预测系数a_k进行表达，并且其表达式具有对称性。由实采序列可求出a_k的估计值，因此对预测模型表达式做Z变换可得到下式：

式中，Ψ(z)是无噪声谐波表达式通过线性预测模型进行表示的形式，进而使用Z变换后得到的原式的表达式，其2K个根均为z平面上单位圆的共轭复根。因此可由Z域变量z的表达式得到余弦信号频率ωk可由下式求得：

利用系数a_k的对称性，无噪声谐波叠加式所表示的N-2K个方程中只只有K个待求解系数{a₁，...a_k}。写成矩阵形式：

并记为：

Aa＝b

需要注意上述矩阵方程中，左右两边的量均可确定，唯一不确定的量就只有系数a_k。假设系数矩阵A的行数N-2K大于列数K，则以上无噪声干扰时的超定方程组可以采用最小二乘法求解。但是在实际应用中，由于实采数据{x(n)]存在测量误差，因此Aa≈b中A矩阵及向量b均受到噪声的干扰，因此应当采用一种新型最小二乘法求解。本实施例方法所使用的新构造的最小二乘方法的目标在于寻求a_k的最佳值使得下式中的扰动矩阵E和扰动项r的总误差

最小，这里

为矩阵的F-范数：

(A+E)a＝(b+r)

上式中扰动矩阵E和r的总体误差

最小，A，b为数据矩阵。将上式表示为增广形式：

并且：

C＝[A，b]∈C^{(N-2K)×(K+1)}

L＝[E，r]

因在实际测量过程中不可避免地存在噪声，因此C通常是满秩的；记C₁＝C+L为扰动矩阵的增广形式；d＝[a，-1]^T为未知向量。因此，上述增广矩阵又可表示为：

(C+L)d＝C₁d＝0

根据上式，使用本实施例所构造的新型最小二乘方法求解未知量d的最佳值的问题又可表述为：计算或搜寻最小F-范数矩阵L使C₁＝C+L为非满秩形式；与此同时，应当使C₁的零空间向量d的最后一个元素非0，从而通过d的最后一个元素按照

的求解方法归一化地求解a的值。

应当注意：因为在矩阵

中，假设的(N-2K)≥(K+1)，因此矩阵C₁的秩最大应为(K+1)，而要解出增广形式矩阵的最佳值，可将矩阵C₁最后一列可以全部可以化为0，这样，矩阵C₁便为非满秩矩阵。

上述最佳值的求解问题可由矩阵C的奇异值分解进行计算。将C写成奇异值分解表示：

显然上述表示形式有最小F-范数，可使C降秩的矩阵L便是上式中的最后一项的相反数

C的零空间向量d＝v_K+1与

部是正交的，k＝1，...，K。因此，将v_K+1的最后一个元素对v_k进行归一化处理，即可求出上述方程组的最小二乘的最佳值a。

根据上述分析，本实施采用所构建的新型最小二乘方法进行叠加的余弦信号中的盲频率估计的过程可以归结为以下步骤：

第一步：设置移相值为π/4，使用波长移相干涉仪采集干涉图200帧；

第二步：计算每个像素点的平均值，并且每个像素点减去该值以消除噪声和背景光强，得到纯谐波形式的混叠的干涉光强；

第三步：将纯谐波形式的数据矩阵A，b构造其增广表达式C＝[A，b]∈C^{(N-2k)×(K+1)}，其元素C(i，j)为：

第四步：对上述C矩阵进行奇异值分解C＝U∑V^T，并对第(K+1)个向量的最后一个元素

进行归一化操作得到a_k；

第五步：按式所构造的ψ(z)求其根z_k，最后由式

求出各信号的频率，带入到下面的幅值和相位求解过程中。根据在多表面干涉过程中的光程差对于干涉频率的影响，当求解出频率的数据组以后，在进行频率匹配时，应当把该数据组进行降序排列，此时各表面信号所对应的频率分量为：后表面、前表面、厚度变化信号。

下面计算叠加的余弦信号的对比度和初始相位。首先，若式中有K个频率相近的信号，对于透明平行平板而言，K＝3时即可得到被测件的前表面、后表面和厚度变化信号的初始相位，进而可以得到面形分布，t＝t₁，t₂，t₃……t_N表示采集到的各帧干涉图的时刻，相邻时刻相差的间距相同，在其含义上可以理解为干涉图在不同时刻下的帧序数，在数值上与上面提及的n相同，为区分两个求解过程此处用t表示。本实施例方法取200帧干涉图，即N＝200。其表达式可写为下式形式：

A₁cos(ω₁t+θ₁)+A₂cos(ω₂t+θ₂)+...+A_mcos(ω_mt+θ_m)＝A₁cosω₁tcosθ₁-A₁sinω₁tsinθ₁+A₂cosω₂tcosθ₂-A₂sinω₂tsinθ₂+...+A_mcosω_mtcosθ_m-A_msinω_mtsinθ_m

设x₁＝A₁cosθ₁，x₂＝A₁sinθ₁，x₃＝A₂cosθ₂，x₄＝A₂sinθ₂...，x_2m＝A_msinθ_m，这样待估计的参数成为x₁，x₂，x₃，…x_m，各次的实采值为b₁，b₂，b₃……b_N，则：

在对实际的测量操作过程中，为了能够获得较为准确的求解结果，完整地对各原始组成信号进行重建，此处建议在保证采样定理的前提下，一般应选取实采点数Z远大于待求解系数。本实施例所采取的构造的新型最小二乘求解方法，即是通过使2-范数：

a_ij是系数矩阵

中的元素。使用上式趋向于最小的形式便可以构造最小二乘误差，再引入初等反射变换法来对上式进行处理。在解线性最小二乘问题时，多数情况下需要对正规方程进行求解，但是如果正规方程的条件数较大，则求解精度不能保证。如果采用初等反射变换，则可避免正规方程的计算。假设U是m维向量，U＝[u₁，u₂...，u_m]^T。用U构成的下列m×N方阵Q称为初等反射矩阵：

因此，经过N次初等反射变换后得到新元素，N为系数矩阵A的列数。应当注意：经过若干次初等反射变换以后，所得到的新的矩阵的范数与进行变换之前是不改变的，但是形式已经化为部分三角矩阵，目的在于使得计算更为简单方便。在新组成的矩阵中，非零元素的列数不变，但是得到的经过上述初等反射变换以后的矩阵的非零行数M＜m。

信号相位与幅值提取。首先，若考虑由3组余弦信号组成的混叠的干涉信号，此时引入一个N×N维的正交矩阵Q_H，为6个初等反射变换阵的乘积，是总的初等反射变换矩阵，使得：

因为||Q_H||＝I所以有：

||Ax-b||₂＝||Q_HAx-Q_Hb||₂＝||Tx-Q_Hb||₂

令

且取

的前6个分量，有：

式中，m＝6，如前所述本实施例选择采集200帧干涉图，即N＝200。右边的第2项不含有待估计参量x，为使上式达到最小值，取

由此可得x的最小二乘估计式为

然后进一步得到各个信号的幅值和初始相位的估计值，其中角标为1，2，3分别对应后表面、前表面和厚度的信息：

求解得到的初始相位如图3所示。在得到各表面的初始相位以后，因为所得到的是包裹结果，因此需要进行一次解包裹操作。解包裹以后，得到前表面、后表面、厚度变化三个信号所对应的初始相位为θ_f、θ_r、θ_f-r通过如下操作即可得到三个信号的对应的面形分布。

其中，n1代表被测件折射率，h代表干涉仪参考镜的后表面到达被测件外表面的距离，W_f、W_r、W_f-r分别代表被测件的前表面、后表面和厚度变化信号的面形分布。

本实施例方法的另一功用在于，能通过求取的各信号的频率分布，通过光程差的计算，反算出被测件的腔长值以及其平均厚度。

其中Δλ为单步移相的波长移相干涉仪中激光器的调谐量，人为进行设置。λ₀为起始调谐波长。n₁是被测件的折射率，此值与材料有关，是已知量。在计算腔长和被测件的平均厚度d时，因为后表面的光程h+n₁×d是与腔长h和被测件光学厚度n₁×d的相加值，因此若其中一个值估算不准确则会影响总体估算值的精度，因此本实施例所采取的方法是，先将前表面信号的干涉频率带入到上式计算腔长值，而后将后表面信号计算的腔长与被测件光学厚度的和的值减去通过前表面干涉频率计算的腔长值，剩下的部分与厚度变化信号干涉频率所估算的光学厚度值进行平均值的计算，从而进一步减少估算误差，得到更为精确的被测件的光学厚度，从而与折射率做差求得被测件的平均厚度d。

至此，透明被测件的多表面测量就完成了。图7所示即为求解得到的面形分布结果，纵向高度方向的单位是微米。若取仿真时的真值作为对比，则各表面求解的面形的误差如图4所示，其高度方向的单位是微米。所求解得到的对比度的分布和误差分布如图6和图5所示。由以上的求解结果和误差结果可看出，本实施例所提出的方法不仅可以完成多表面被测件的一次性干涉测量，并且可在被测件平均厚度未知的情况下准确得到各表面的面形分布。

上面对本发明实施例结合附图进行了说明，但本发明不限于上述实施例，还可以根据本发明的发明创造的目的做出多种变化，凡依据本发明技术方案的精神实质和原理下做的改变、修饰、替代、组合或简化，均应为等效的置换方式，只要符合本发明的发明目的，只要不背离本发明的技术原理和发明构思，都属于本发明的保护范围。

Claims (8)

1.一种基于频率盲估计的多表面面形测量方法，其特征在于：通过对采集得到的干涉图中的各点干涉频率的盲估计和幅值、相位的求解，完成具有多表面的被测件的一次性无接触测量；使用波长移相干涉仪对被测件进行干涉图的采集，所述基于频率盲估计的多表面面形测量方法包括如下步骤：

a.对叠加的余弦信号的光强变化条纹图进行基于各信号干涉频率独立分布的盲频率估计，从而确定各个干涉信号的精确频率分布；

b.通过所构建的新型最小二乘初始相位和幅值求解方法进行计算，从而恢复各个干涉信号中所包含的面形信息，对透明的平行平板进行多表面分离和重建，并对被测件的腔长值和平均厚度进行求取，从而完成被测件的多表面干涉测量。

2.根据权利要求1所述基于频率盲估计的多表面面形测量方法，其特征在于：在进各个信号的干涉频率的求解时，考虑由前后表面和厚度变化的叠加的干涉图，根据激光传播特性，采集的到的条纹图表示为三组余弦信号相叠加的形式，将

表示为由K个频率不同的余弦信号及噪声e(n)叠加组成的信号，写为以下表达式：

式中，实采数据的序数为n＝1，...，N，{A_k，ω_k：k＝1，...K}是未知实量，θ_k：k＝1，...K是在-π～π之间的独立随机变量，具有均匀分布特性，并且对应各个信号的初始相位；ω_k是各个信号的频率；A_k是各信号的对比度；所述盲频率估计是通过余弦检测从给定的N个实采数据

中精确计算各未知数；假设e(n)的均值为零、方差为σ²的噪声分量；采用消除噪声的方法为：对每一个像素点计算各帧干涉图下的平均值，将采集到的干涉中各点减去此值，消除背景光强和噪声；移相值设置为π/4，为了保证测量精度，总采集帧数选择为200帧，即N＝200；

对因此将无噪声谐波y(n)及其组成表示为：

改写为下式：

式中，

其中e是自然底数，j是虚数单位；因此当频率分量ω_k求出后，以上N个方程组是b_k的线性方程，因此求出各个b_k是具有可行性的，并由实采序列

的各组成信号的频率进行计算和分析；进而，任意频率的K个无噪声谐波的实数余弦信号可用线性预测模型中的预测系数a_k进行表达，并且其表达式具有对称性；由实采序列可求出预测系数a_k的估计值，因此对预测模型表达式做Z变换得到下式：

式中，Ψ(z)是无噪声谐波表达式通过线性预测模型进行表示的形式，进而使用Z变换后得到的原式的表达式，其2K个根均为z平面上单位圆的共轭复根；因此可由Z域变量z的表达式得到余弦信号频率ωk可由下式求得：

3.根据权利要求1所述基于频率盲估计的多表面面形测量方法，其特征在于：在进行各个信号的幅值和相位求解时，考虑由3组余弦信号组成的混叠的干涉信号，此时引入一个N×N维的正交矩阵Q_H，为6个初等反射变换阵的乘积，是总的初等反射变换矩阵，使得：

因为||Q_H||＝I，所以有如下公式：

||Ax-b||₂＝||Q_HAx-Q_Hb||₂＝||Tx-Q_Hb||₂

令

且取

的前6个分量，有如下公式：

式中，m＝6；选择采集200帧干涉图，即N＝200；上式的右边的第2项不含有待估计参量，为使上式达到最小值，取

由此可得待估参量的最小二乘估计式，用符号x表示为

然后进一步得到各个信号的幅值和初始相位的估计值，其中角标为1，2，3分别对应后表面、前表面和厚度的信息，参见如下公式：

求解得到的初始相位；在得到各表面的初始相位以后，因为所得到的是包裹结果，因此需要进行一次解包裹操作；解包裹以后，得到前表面、后表面、厚度变化三个信号所对应的初始相位为θ_f、θ_r、θ_f-r，通过如下操作得到三个信号的对应的面形分布公式：

其中，n1代表被测件折射率，h代表干涉仪参考镜的后表面到达被测件外表面的距离，W_f、W_r、W_f-r分别代表被测件的前表面、后表面和厚度变化信号的面形分布。

4.根据权利要求1所述基于频率盲估计的多表面面形测量方法，其特征在于：在进行被测件的腔长值和平均厚度的求解时，通过求取的各信号的频率分布，通过光程差的计算，反算出被测件的腔长值以及其平均厚度，参见如下公式：

其中Δλ为单步移相的波长移相干涉仪中激光器的调谐量，人为进行设置；λ₀为起始调谐波长；n₁是被测件的折射率，此值与材料有关，是已知量；后表面的光程h+n₁×d是与腔长h和被测件光学厚度n₁×d的相加值，在计算腔长和被测件的平均厚度d时，将前表面信号的干涉频率带入到上式计算腔长值，而后将后表面信号计算的腔长与被测件光学厚度的和的值减去通过前表面干涉频率计算的腔长值，剩下的部分与厚度变化信号干涉频率所估算的光学厚度值进行平均值的计算，从而进一步减少估算误差，得到更为精确的被测件的光学厚度，从而与折射率做差求得被测件的平均厚度d。

5.根据权利要求2所述基于频率盲估计的多表面面形测量方法，其特征在于：根据在多表面干涉过程中的光程差对于干涉频率的影响，当求解出频率的数据组以后，在进行频率匹配时，应当把该数据组进行降序排列，此时各表面信号所对应的频率分量为后表面、前表面、和厚度变化信号。

6.根据权利要求3所述基于频率盲估计的多表面面形测量方法，其特征在于：进行初始相位和幅值的求解，将所构造的最小二乘系数矩阵经过N次初等反射变换后得到新元素，N为系数矩阵A的列数；经过若干次初等反射变换以后，所得到的新的矩阵的范数与进行变换之前是不改变的，但是形式已经化为部分三角矩阵，在新组成的矩阵中，非零元素的列数不变，但是得到的经过初等反射变换以后的矩阵的非零行数M＜m。

7.根据权利要求1所述基于频率盲估计的多表面面形测量方法，其特征在于：所述新型最小二乘初始相位和幅值求解方法的目标在于寻求a_k的最佳值使得下式中的扰动矩阵E和扰动项r的总误差

最小，这里

为矩阵的F-范数：

(A+E)a＝(b+r)

上式中扰动矩阵E和扰动项r的总体误差

最小，A，b为数据矩阵；将上式表示为增广形式：

并且：

C＝[A，b]∈C^{(N-2K)×(K+1)}

L＝[E，r]

因在实际测量过程中不可避免地存在噪声，因此C是满秩的；记C₁＝C+L为扰动矩阵的增广形式；d＝[a，-1]^T为未知向量；上述增广矩阵又可表示为：

(C+L)d＝C₁d＝0

根据上式，使用所构造的新型最小二乘方法求解未知量d的最佳值的问题又表述为：计算或搜寻最小F-范数矩阵L使C₁＝C+L为非满秩形式；与此同时，应当使C₁的零空间向量d的最后一个元素非0，从而通过d的最后一个元素按照

的求解方法归一化地求解a的值；

在矩阵

中，假设(N-2K)≥(K+1)，因此矩阵C₁的秩最大为(K+1)，将矩阵C₁最后一列全部化为0，这样矩阵C₁便为非满秩矩阵，从而解出增广形式矩阵的最佳值；

上述最佳值的求解问题由矩阵C的奇异值分解进行计算，将C写成奇异值分解表示：

上述表示形式有最小F-范数，可使C降秩的矩阵L便是上式中的最后一项的相反数

C的零空间向量d＝v_K+1与

都是正交的，k＝1，...，K；因此，将v_K+1的最后一个元素对v_k进行归一化处理，即可求出上述方程组的最小二乘的最佳值a。

8.根据权利要求7所述基于频率盲估计的多表面面形测量方法，其特征在于，采用所构建的新型最小二乘方法进行叠加的余弦信号中的盲频率估计的过程包括以下步骤：

第一步：设置移相值为π/4，使用波长移相干涉仪采集干涉图200帧；

第二步：计算每个像素点的平均值，并且每个像素点减去该值以消除噪声和背景光强，得到纯谐波形式的混叠的干涉光强；

第三步：将纯谐波形式的数据矩阵A，b构造其增广表达式C＝[A，b]∈C^{(N-2K)×(K+1)}，其元素C(i，j)为：

第四步：对上述C矩阵进行奇异值分解C＝U∑V^T，并对第(K+1)个向量的最后一个元素

进行归一化操作得到a_k；

第五步：按式所构造的ψ(z)求其根z_k，最后由式

求出各信号的频率，带入到通过对采集得到的干涉图中的各点干涉频率的盲估计和幅值、相位的求解过程中，根据在多表面干涉过程中的光程差对于干涉频率的影响，当求解出频率的数据组以后，在进行频率匹配时，将该数据组进行降序排列，此时各表面信号所对应的频率分量为：后表面、前表面、厚度变化信号。

Images