本发明涉及一种基于频偏修正的多干涉谐波叠加信号频率及相位参数重建方法，结合近似熵量化计算了移相干涉测量中的非同步采样问题，利用序列重排技术实现了具有更低频谱泄漏的谐波频率分析，并且据此实现了多表面干涉谐波的幅值、频率修正和欠采样识别，应用在精密激光移相干涉测量领域。
背景技术：
1、干涉测量是一种高精度的测量技术，通过将被测波形与参考波形进行干涉，进而通过一定的方法提取被测波面与参考波面的初始相位，最终实现对被测面表面的精确测量。其过程如下：首先，将被测件放置在干涉仪前方，使用干涉仪将被测波面与已知波面进行干涉，形成干涉图样，并通过ccd相机等硬件进行干涉图的存储。然后，通过测量干涉图样的变化(通常需要进行移相)，可以确定被测波的波长、相位和振幅等参数。
2、具体而言，当被测波面与已知波面(也可称为参考波面)在干涉仪中相遇时，它们会发生干涉，形成稳定的干涉条纹图样，人为引入主动的相位变化值(也即移相)获得一系列具有相同初始相位、同时也具有不同的主动相位的干涉条纹图(也可以称为干涉强度图，通常简称为干涉图)作为数据基础。由于干涉图样的光强变化值与被测波面的波长、相位和振幅等参数密切相关，因此基于采样得到的具有不同相位值的干涉光强图，通过相应的参数解调方法，可以实现上述参数的精确求解。干涉测量术属于非接触式光学测量法，具有高精度、高灵敏度等优点，广泛应用于微纳结构识别与声、光检测等领域。
3、在主动相位引入方式方面，传统的硬件移相法通过改变干涉仪前方的参考镜的微观距离从而改变被测件和参考镜之间的光程差(距离值与材料折射率的乘积)，可以实现主动相位的改变。但是这种移相方式存在较多问题，具体如下：
4、(1)相移器稳定性问题：硬件移相法需要使用压电结构推动参考镜(称为移相器)，而相移器的稳定性直接影响干涉测量的精度。但是，在实际测量时，相移器容易受到温度、机械振动和老化等因素的影响，从而导致测量误差；
5、(2)硬件移相的实现较为复杂：硬件移相法需要进行对压电结构进行主动电压输入和输出修正等操作，复杂度较高，需要使用专门的电路和设备，增加了测量成本和测量难度；
6、(3)测量范围有限：硬件移相法的测量范围受到相移器移相范围的限制，对于某些需要大范围或者大口径相位调制的干涉测量，硬件移相法可能无法满足要求；
7、(4)难以实现高精度相位测量：受到上述不良因素的影响，硬件移相法的相位测量精度也受到限制，难以实现高精度相位测量；
8、(5)多表面叠加信号无法求解的问题，当被测件是透明且前后表面具有较高平行度的元件时，采集得到的干涉图是各表面干涉的叠加结果，硬件移相法无法处理该种干涉条纹图。而这种多表面元件在光学系统中应用广泛，正是本专利所针对的测量目标。
9、而近些年获得较快发展的波长移相法可以较好地克服上述问题。波长移相干涉测量技术是通过改变光源的波长进行主动相位的改变(波长与主动相位呈线性关系)，也即移相，并在改变波长的同时进行干涉图采集。
10、这种移相方法的优势在于：激光器作为移相器和光源，不需要引入具有复杂结构的硬件移相装置，因此干涉仪结构较为简单，可以克服硬件移相所具有的机械迟滞效应。并且由于避免了参考镜的推动，因此该种移相方法可以实现大口径被测件的测量。最重要的是，在对具有多个表面的透明平板元件进行检测时，因为平板元件前后表面平行度较高，各被测面与参考镜中的参考件均会产生干涉，采集得到的干涉图是各个被测面表面干涉信息叠加的结果，传统的硬件移相法得到的干涉图无法进行参数解调。但通过波长移相法得到的干涉图中，各个被测面对应干涉信号(称为谐波)因为具有不同的光程差，因此可以在时域中展现出不同的谐波频率，这为多个被测表面的同时重建提供了可能性。
11、在干涉参数的求解方面，通过传统的离散傅里叶变换对谐波信号的频率和相位进行求解是一种常用的经典方法。但是这种方法的缺陷在于：由于波长移相过程中(也即采样过程中)，波长的改变量并非完全线性的，这就会造成非同步采样问题。直观而言，这是由于谐波信号的采样长度与信号周期不成整数倍关系所导致的，完全的同步采样在实际测量中受到各种误差的影响，也是不可能存在的。而非同步采样问题的存在，会导致传统离散傅里叶变换后的频率幅值谱的频谱泄漏，致使谐波频率、相位和幅值的求解精度受到损失，非常不利于高精度检测的实现。同时，在实际测量过程中，非同步采样问题是难以量化识别的。
12、此外，根据多频谐波采样定理，当谐波频率大于的采样频率一半时，则在单边频率谱中，该种高频信号的频率曲线会以二分之采样频率为轴线折叠到探测到的频谱中，这种现象可以称之为欠采样频率模糊现象或频谱折叠问题，此时探测到的频率是假频，基于这种频谱得到的频率和相位等参数信息，会致使测量的失效，无法反映被测件的真实面形特征。
13、在目前已存在的技术和装置中，专利发明人常林等人在发明专利《一种用于光学多表面透明板夹持的可旋转夹具-申请号：202010595035.x》、《一种用于光学透镜的一体式夹持及夹紧机构-申请号：202110177033.3》以及论文《lin chang，yingjie yu，et al.ageneral auto-shift minimal-step phase-shifting algorithm for arbitrary cavitylength[j].optics and lasers in engineering，2022，149:106791.(sci二区，第一作者)》、《lin chang，yingjie yu，et al.multi-surface phase-shifting algorithm usingthe window function fitted by the nonlinear least squares method[j]journal ofmodern optics，2022，69(3):160-171.(sci四区，第一作者)》中所设计的新型系统装置和测量方法实现了被测件的便捷夹持、灵活位置下的各表面光学检测等目标。
14、但是，上述专利和论文中的创新结构和方法无法解决上面所指出的问题，具体包括：
15、(1)存在非同步采样问题的量化识别；(2)具有更低频谱泄漏的频率分析的实现；(3)精确的干涉光强对应频率、幅值修正；(4)精确的各测量表面的同时重建；(5)欠采样频峰模糊问题发生时的高精度识别。
技术实现思路
1、为解决上述问题，本发明设计一种基于频偏修正的多干涉谐波叠加信号频率及相位参数重建的方法，可以以更低的频谱泄漏实现高精度的谐波频率求解和各表面同时面形重建，并可以通过近似熵对非同步采样程度进行量化。
2、对材质透明且前后表面具有高平行度的光学元件进行检测的大致过程可以归纳为：将被测件夹持在夹具上并且放置在干涉仪前方使待测元件与干涉仪中的参考镜平行以产生强度叠加干涉图、打开干涉仪并且利用外部控制器驱动光源改变波长实现移相同时采集干涉图、对干涉图进行处理获得各表面对应的包裹相位、对包裹相位进行解包裹和消倾斜等操作最终实现各表面对应波面高度的重建。
3、本发明所设计的测量方法的主要流程如下：
4、步骤1：通过波长调谐移相干涉仪采集2n-1帧干涉图，其中n应为不小于20的整数，选取干涉图中像素点(xp、yp)的光强记为ip进行分析，设定近似熵阈值且对ip进行近似熵计算，步骤1主要包括(1.1)向量序列分组、(1.2)向量距离量化计算、(1.3)量化向量相似度、(1.4)近似熵量化求解、(1.5)主频率搜索5个主要过程，具体计算步骤如下：
5、优选地，近似熵(approximate entropy，apen)主要用于评估信号序列的复杂性和不规则性的参数。时间序列信号复杂度与近似熵的结果成正比，也即，越规则的序列则近似熵越小。这种计算方法的一个显著的优势在于它对数据量的需求较小，且计算结果稳健、可靠，这意味着大多数实际测量的时间序列都满足其要求，从而可以获得稳定可靠的结果，是对干涉测量中的光强序列规则性进行量化的一种合适的方法。
6、(1.1)向量序列分组：将光强ip带入x(na)中，对于由na个光强数据组成的时间序列x(na)＝x(1),x(2),…,x(na)，近似熵的计算方法如下：
7、解释上述过程如下，各光强数据也即：某像素点在不同帧序数(na)下的光强变化即为一时间序列，na即为此时对应的总帧数。
8、按序号组成一组维数为m的向量序列g，具体表示为其中向量序列g代表从第ia点开始的m个连续的x(na)中的元素值，例如当m＝2时，
9、(1.2)向量距离量化计算：优选地，定义向量与的距离为：
10、
11、上式中的max为求解最大值，||为计算绝对值。
12、(1.3)量化向量相似度：优选地，对于一组与计算小于等于误差半径ra的数量并将结果记为ci，进而计算εm,r：
13、
14、其中，m和ra可以根据需求进行设定，总体误差变动范围用ga表示，数值上ga＝2ra，但设定ra时应满足ra处于0.1fsd(x)和0.25fsd(x)之间的要求，其中fsd(x)为序列的x(na)的标准差；m一般设定为2。
15、(1.4)近似熵量化求解：优选地，增加维数到m+1，重复步骤(1.1)～(1.3)，得到εm+1,r，并通过下式计算近似熵fapen(m,r)：
16、fapen(m，r)＝εm，r-εm+1，r (式-3)
17、(1.5)主频率搜索：对光强信号进行ip补零操作，补零长度为5n，补零后序列的总长度为7n-1且记为ip0，对ip0去均值后进行傅里叶变换并取单边幅值谱进行3个主要谐波峰值对应横坐标搜索，且将各个谐波峰值的横坐标记为主频率值。
18、至此，步骤1计算过程结束，如果近似熵大于预设熵值，则转至步骤2。
19、步骤1所设计的干涉光强序列的近似熵计算方法，其有益效果是能够有效识别非同步采样问题，可以为后续的进一步频谱分析和处理提供可靠的量化基础。
20、步骤2：优选地，将光强ip带到x(n)中以后进行下一步处理，步骤2主要包括(2.1)构造重排序列、(2.2)序列重排fft函数求解、(2.3)加窗序列频谱函数求解、(2.4)构造相位差、(2.5)频率和幅值修正五个具体步骤，流程如下：
21、以x(n)(此时，只截取n≥0的部分进行分析，n＝0,1,…,n-1为时间序列的序号)为单频复指数信号作为示例展示序列重排傅里叶变换的计算过程：
22、
23、上式中，信号的数字频率ω0表示为β倍频率间隔2π/n的形式，其中β可以是小数、θ为初始相位、n为总序列点数，j为虚数单位，则经过归一化以后，x(n)的原始fft频谱函数(fft即为快速傅里叶变换，是离散傅里叶变换的一种快速计算形式)可表示为：
24、
25、上式中，k＝0,1,…,n-1。进而可将上式展开为：
26、
27、(2.1)构造重排序列：序列重排fft是考虑包含序列中某点在所有循环移位后的分段傅里叶变换而产生的一种分析方法。
28、优选地，所述的序列重排是对包含着序列中的某个点的向量进行重新排序，进而进行循环抽取和排列，以x(0)为例，存在且只存在n个包含该点的n维向量：
29、
30、优选地，上述公式中，x0到xn-1右下角的角标表示序列重排后的向量序号，将每个向量进行循环排序且把序列值x(0)移到序列中的第一位，则可得到另外的n个n维向量：
31、
32、优选地，对准x(0)相加并取其平均值，则可得到重排后的数据向量，以xap表示如下：
33、
34、(2.2)序列重排fft函数求解：根据傅里叶变换的移位性质，以角标ka表示向量序号，则x’ka(n)的离散傅里叶变换x’ka(k)和x(n)的原始fft频谱函数xka(k)之间关系如下：
35、
36、优选地，对上式的x′k(k)进行求和平均即为序列重排fft函数：
37、
38、
39、(2.3)加窗序列频谱函数求解：对于长度为n的序列，若用窗函数wf对单频复指数信号进行加窗fft(注意，这里仍然是以n≥0为例截取长度为n的序列值进行分析和展示)，以函数fg为窗函数的频域形式，则得到的离散传统fft的谱表达式为：
40、
41、上式中，δω＝2π/n为频率分辨率单位，群延迟量为τ＝(n-1)/2，根据fft的lti特性，对幅值为a、初相为θ、时间序列延迟为n0、角频率为ω*的复指数信号再进行加窗fft谱分析，则其离散谱表示为：
42、
43、若对序列进行n个值和后n个值都执行加窗操作，则前、后窗函数均为wf的双窗序列重排fft的频谱函数表示为：
44、
45、则根据序列重排fft的线性特性，对信号进行同样的双窗序列重排fft谱分析，其离散谱表示为：
46、
47、(2.4)构造相位差：优选地，以k*代表频率幅值谱中主频率谐波峰值的主谱线序号值，进而将传统相位谱值减去序列重排后的相位谱值则可得到如下的相位差表达式如下：
48、
49、(2.5)频率和幅值修正：优选地，测出传统fft和序列重排的两种谱分析在主谱线上的相位大小并作差，便可计算出频谱分析得到的频率与实际频率之间的偏移值，进而通过下式对信号真实频率的估计：
50、
51、此处注意，若不只取主谱线的相位，也可以绘制相位差值的曲线。
52、分别取(式-13)和(式-15)的主谱线上的模值，有
53、
54、将(式-15)进行进一步变换，即可得到频率幅值的修正值：
55、
56、由(式-17)至(式-19)就可完成光强信号ip中3个谐波的频率、幅值估计，以频率值增序频排列，则经修正后的3个谐波的频率和幅值分别记为fp1、fp2、fp3和ap1、ap2、ap3。
57、在上面的计算中，窗序列wf可以根据用户需求进行选取，常用的窗函数有汉宁窗、海明窗、矩形窗、三角窗、布莱克曼窗等，默认选择汉宁窗。
58、步骤2所设计的干涉光强序列的序列重排fft以及频率、幅值修正方法，其有益效果是能够显著降低频谱泄漏，并且可以得到精确的谐波频率和幅值信息，可以为后续的相位解调、波面重建以及欠采样的判定打下坚实基础。
59、至此，步骤2计算过程结束，转至步骤3。
60、步骤3：频谱折叠判定及相位求解：当(ap2-ap3)/ap1大于折叠阈值ejs时，则判定不存在频谱折叠，若(ap3-ap2)/ap1大于折叠阈值ejs则存在频谱折叠的问题，若满足|(ap2-ap3)/ap1|≤ejs则无法准确判断是否发生频谱折叠，将像素点移动到(xp+δp、yp+δp)后再次执行步骤1和步骤2；
61、折叠阈值可根据用户需求进行设定，一般应不小于0.2且不大于0.5。
62、优选地，当存在频谱折叠时，fqs-1＝fp1为厚度变化谐波的相位计算频率值、fqs-2＝fp2为前表面谐波的相位计算频率值、fqs-3＝fp3为后表面谐波的相位计算频率值；
63、优选地，当不存在频谱折叠时，fqs-1＝fp1仍然为厚度变化谐波的相位计算频率值、fqs-2＝fp3为前表面谐波的相位计算频率值、fqs-3＝fp2为后表面谐波的相位计算频率值；
64、以为各谐波的初始相位、bp和ap分别是各谐波对应的相位求解函数，进而干涉图中的初始相位的求解公式可以写为：
65、
66、上式中u0表示干涉图中的光强，u0中包含着5n帧元素全为0的干涉图。各谐波对应的bp和ap分别通过下面式子确定：
67、
68、
69、上式中re和im分别表示求实部和求虚部，fqs-p、bp和ap中的右下角角标p＝1、2、3时分别代表厚度变化信号、前表面信号和后表面信号。应注意，此时(式-20)和(式-21)中的n＝0，1，2，...7n-1。
70、求解出各谐波的包裹相位以后，对包裹相位进行解包裹、消倾斜以后便可根据波面高度和相位之间的关系进行各谐波面形的重建。
71、步骤3所设计的相位解调和欠采样识别方法，其有益效果是能够高效实现干涉图中的初始相位提取，并且可以根据前面计算得到的精确的谐波频率和幅值信息，对欠采样问题进行有效识别和处理，对于精确的多表面波面重建是非常便利的。
72、至此，步骤3计算过程结束，测量完成。
73、本发明与现有技术相比较，具有如下显而易见的突出实质性特点和显著优点：
74、(1)结合所采用的近似熵计算法和频率幅值修正法，能够实现非同步采样问题的量化评估；
75、(2)所设计的欠采样时频谱折叠的判定方法使所设计的技术的适用范围更加广泛，有效克服欠采样频谱折叠问题发生时的相位解调和频率对应；
76、(3)所采用的频率和幅值修正方法能够显著提高谐波参数重建精度；
77、(4)所采用的序列重排数据组合方法能够有效降低频谱泄漏；
78、(5)通过本专利所采用的近似熵求解、频率幅值修正和相位解调的组合方法，可以通过一组干涉图同时对被测件的多个波面进行高精度重建。
1.一种基于频偏修正的多干涉谐波叠加信号频率及相位参数重建方法，其特征在于，包括以下步骤：
2.根据权利要求1中所述的一种基于频偏修正的多干涉谐波叠加信号频率及相位参数重建方法，其特征在于：步骤1中的近似熵判定的预设熵值可根据用户需求选择，默认选择0.3；步骤2中的窗序列wf可以根据用户需求进行选取，默认选择汉宁窗。
3.根据权利要求1中所述的一种基于频偏修正的多干涉谐波叠加信号频率及相位参数重建方法，其特征在于：在步骤3中，若满足|(ap2-ap3)/ap1|≤ejs则无法准确判断是否发生频谱折叠，将像素点移动到(xp+δp、yp+δp)后再次执行步骤1和步骤2；同时，折叠阈值可根据用户需求进行设定，一般应不小于0.2且不大于0.5。
