Disturbance quantification model and rooting method combined multi-parameter estimation algorithm in interference measurement

Abstract

The invention relates to a multi-parameter estimation algorithm combining a disturbance quantification model and a rooting method in interference measurement. In order to solve the problems that in the prior art, the disturbance degree of a sequence is difficult to quantify, and harmonic key parameters are influenced by spectrum leakage, the invention provides the following innovative scheme: firstly, a sequence multi-classification measurement entropy quantification model is constructed, and quantitative evaluation of the signal noise degree is realized; secondly, through a Hankel matrix construction and rooting parameter estimation method, realizing accurate estimation of the harmonic frequency of the real-value signal; meanwhile, an amplitude correction algorithm based on window function characteristics is designed, and the signal sampling condition is effectively judged; and finally, through matrix weighting operation of a phase demodulation function, efficient reconstruction of an optical surface shape is realized, the problems of spectrum leakage and disturbance quantification in a traditional method are solved, and the method has relatively high measurement precision.

Description

Technical Field

[0001] The present invention relates to the field of optical measurement technology, and specifically to a multi-parameter estimation algorithm that combines a perturbation quantization model with a root-finding method in interferometric measurement. In response to the problems in the prior art where the degree of perturbation of a sequence is difficult to quantify and where key harmonic parameters are affected by spectrum leakage, the present invention proposes the following innovative solutions: first, a multi-classification metric entropy quantization model for the sequence is constructed to achieve quantitative evaluation of the degree of signal noise; second, by constructing a Hankel matrix and a root-finding parameter estimation method, accurate estimation of the harmonic frequency of a real-valued signal is achieved; at the same time, an amplitude correction algorithm based on window function characteristics is designed to effectively determine the signal sampling situation; finally, through the matrix weighting operation of a phase demodulation function, efficient reconstruction of the optical surface shape is achieved, solving the spectrum leakage and perturbation quantization problems in traditional methods and achieving high measurement accuracy. Background Art

[0002] Surface inspection of key equipment components is a crucial step in assessing machining accuracy. In non-contact measurement and 3D imaging technologies, demodulating the parameters of the photoelectric signal carrying the measured object is crucial for reconstruction performance. Optical interferometry, in particular, can theoretically capture the topographic distribution of the measured surface at the nanometer level across the entire measuring aperture, and is rapidly gaining traction in the field of high-precision machining.

[0003] Optical interferometry is based on the principle of interference. It uses an interferometer to collect images of the interference signal between the reference mirror and the north mirror inside the device. The frequency, phase and harmonic amplitude of the interference signal are obtained through a certain parameter demodulation method. This technology has extremely high surface reconstruction accuracy.

[0004] Single-surface measurement of optical components has been extensively studied. However, when the DUT has two parallel surfaces, front and back, the interference of the reflected light beams from each surface can cause the acquired fringe pattern to exhibit harmonic superposition characteristics. In this case, the harmonic signals corresponding to each surface must be analyzed and parameter demodulated to obtain the initial phase information containing the measured surface shape.

[0005] Interferometry can be divided into single-frame measurement technology and multi-frame measurement technology, and the latter can be roughly divided into the hardware phase shifting method achieved by piezoelectric ceramics driving the reference mirror and the wavelength tuning method achieved by changing the wavelength. In the interference pattern obtained by the hardware phase shifting method, since the microscopic displacement is changed, the individual harmonics cannot show independent frequency differences. The wavelength tuning method allows the light source to be changed to achieve phase shifting. The harmonics of each surface show their own time domain frequency due to the difference in optical path difference, avoiding the complex structure introduced by the hardware interferometry method. It can also process the interference pattern of multiple surfaces superimposed, which is suitable for the simultaneous detection of each surface of parallel plate optical elements.

[0006] In optical measurement algorithms, Kim et al. developed a series of weighted sampling algorithms (doi:10.1016 / j.optlastec.2023.110082, doi:10.1007 / s12541-024-01134-0, doi:10.1016 / j.measurement.2020.107870) using characteristic polynomials for window function design. These algorithms are robust to harmonic errors and explore the application of deep learning methods in multi-surface optical testing. Servin et al. developed a phase-shifting interferometry algorithm through the design and analysis of filters and frequency-domain transfer functions, and studied the problem of non-uniform phase shifting (doi:10.1364 / AO.58.001134, doi:10.1364 / ao.58.001134). Least-squares techniques can construct an iterative structure using ideal and actual interferometric signal models to solve for phase and phase shift values. Sun et al. developed this into a multi-surface iterative algorithm, but the initial estimate of the iteration significantly affects the solution (doi:10.1016 / j.optlaseng.2018.11.008). In recent years, this type of iterative algorithm has optimized the iteration structure, selecting only a limited region to calculate the phase shift value, thereby improving computational efficiency. However, when the interferometric signal is adversely affected by the sampling process, how to assess the degree of interference to achieve better measurement results remains a research topic. Eigenvalue decomposition was introduced by Bai et al. to estimate optical parameters for wavenumber tuning techniques, enabling detection of the number of signal sources and constructing a least-squares solution model. This algorithm also incorporates the Estimating Signal Parameter via Rotational Invariance technique (ESPRIT). This type of spectral estimation algorithm has been shown to have high parameter accuracy, but its relationship with the real-valued optical intensity signal needs further clarification and validation, and its efficiency should be considered. Traditional FFT (Fast Fourier Transform) algorithms can usually obtain relatively accurate results when the signal is subject to weak disturbances. However, once the noise and errors in the signal are strong, the resulting spectrum leakage and fence effect problems will greatly reduce the accuracy of solving harmonic parameters such as frequency and amplitude. At this time, the parameter solution process needs to be corrected to obtain more accurate results. However, it is currently difficult to quantify the extent to which the signal is affected by noise and errors.

[0007] Based on the above analysis, the current problems in multi-surface measurement can be summarized as follows: (1) A quantitative evaluation model should be constructed when the sequence is disturbed in order to determine the necessity of parameter estimation correction; (2) When there is measurement disturbance, amplitude information cannot be obtained when solving the frequency, and it is difficult to take both into account (such as ESPRIT and least squares iteration); (3) Traditional spectral correction methods (such as ESPRIT) only use complex exponential signals containing real and imaginary parts as the solution object, which makes it difficult to achieve accurate frequency amplitude solution of real-valued light intensity signals and efficient phase demodulation. Summary of the Invention

[0008] The purpose of the present invention is to overcome the shortcomings of the existing technology and propose a multi-parameter estimation algorithm that combines a perturbation quantization model with a root-finding method in interferometry. The purpose is to quantitatively evaluate the degree of disturbance of the interference signal during the measurement process, accurately solve the frequency amplitude and phase, and thus achieve simultaneous high-precision measurement of the surface morphology of multi-surface optical components. The designed multi-parameter estimation algorithm that combines a perturbation quantization model with a root-finding method in interferometry mainly includes: a multi-classification metric entropy model of the optical interference signal, a harmonic frequency correction method based on root-finding parameters, an amplitude correction method, and a phase demodulation algorithm. The designed combined algorithm first constructs a multi-classification metric entropy model to quantitatively evaluate the degree of disturbance of the interference signal during measurement, and then introduces a spectral estimation method based on the root-finding matrix to obtain the frequency of each harmonic, and uses the window function characteristics to correct the harmonic amplitude. Finally, the initial phase of each surface is efficiently demodulated by weighted calculation through the phase demodulation function.

[0009] The principle of the present invention is that the acquisition of interference patterns in the process of optical interferometry is a key step in the measurement, but it is easily affected by factors such as environmental disturbances and instrument errors during the acquisition process. When the interference signal obtained by the acquisition is directly subjected to FFT (Fast Fourier Transform), the accuracy of the frequency and amplitude is low due to the presence of noise and error, and there is currently little research on the quantitative evaluation model for the degree of interference signal disturbance. To this end, this patent designs a multi-classification metric entropy model that associates the entropy distribution with the degree of sequence disturbance to evaluate the size of noise and error in the interference pattern acquisition process, and accordingly determines whether further harmonic parameter correction is needed. If the judgment result is that the entropy value is less than the set entropy threshold, the current sequence disturbance degree is small, and there is no need to further correct the harmonic parameters. The FFT spectrum of the interference signal is directly used to solve the harmonic frequency and amplitude; if the judgment result is that the entropy value is greater than the set entropy threshold, the current sequence disturbance degree is large, and further parameter correction of the harmonic amplitude and frequency is needed. The designed parameter correction algorithm includes frequency correction and amplitude correction. In the frequency correction part of the parameter correction algorithm, the light intensity distribution of the five pixels at the center of the interference pattern at different sampling frames is first de-averaged. The de-averaged light intensity sequence is then used to construct a Hankel matrix. The covariance and eigenvalues ​​of the Hankel matrix are then calculated, and a root parameter estimation polynomial is constructed to accurately solve the harmonic frequencies of the three signals. In the amplitude correction part of the parameter correction algorithm, the spectrum modulus function is derived based on the window function. The amplitude correction model is then constructed using the frequency and amplitude near the peak in the FFT spectrum to accurately solve the harmonic amplitudes of the three signals. The amplitude solution results are analyzed. If the amplitude of the second harmonic peak from left to right in the spectrum is found to be greater than the amplitude of the third harmonic peak, it is determined that frequency folding has occurred. The true frequency distribution of the second and third harmonics can be obtained by using half the sampling frequency as the axis of symmetry. Finally, a phase demodulation function is constructed based on the obtained harmonic frequency, and a weighted operation is performed on the interference pattern to obtain the wrapped form of the harmonic initial phase corresponding to each surface. The wrapped phase is then unwrapped and de-tilted, and the surface shape is reconstructed using the direct linear correspondence between the wavefront height and the phase, thereby achieving accurate measurement of each surface.

[0010] The wavelength-tunable phase-shifting interferometry measurement technology utilized in this patent achieves phase shifting by changing the wavelength of the light source. When measuring an optical test piece with a high degree of parallelism between its front and rear surfaces, interference occurs between the reflected signals from each test surface in the measurement system, so the collected interference pattern is the result of the superposition of multiple surfaces. Specifically, after the light beam emitted from the light source is collimated and expanded by the optical system in the interferometer, a portion of the beam is reflected by the reference mirror surface (recorded as signal PA0, reflected back into the interference system), while another portion passes through the reference mirror in the interferometer and is reflected by the front surface of the test piece. At this time, a portion of the signal is still reflected (recorded as signal PA1, reflected back into the interference system), and another portion of the signal passes through the front surface of the test piece and reaches the back surface of the test piece and is reflected (recorded as signal PA2, reflected back into the interference system). At this time, the three sets of interference signals interfere with each other. The interference signal between PA0 and PA1 can be recorded as the front surface harmonic signal, the interference signal between PA0 and PA2 can be recorded as the back surface harmonic signal, and the interference signal between PA1 and PA2 can be recorded as the thickness change harmonic signal. Since the optical path differences of the interference signals in the optical interferometer cavity are different, when the wavelength of the light source is changed, each harmonic has its own independent phase shift value, and also exhibits different frequencies in the time domain (or frame domain). The frequency difference is the core difference between the harmonics. Therefore, for a certain point in the interference pattern as the observation object, the interference light intensity at this point presents the form of a superposition of three cosine signals of different frequencies at different frame numbers, and the initial phase parameter is linearly corresponding to the height of the point on the test piece corresponding to this point. At this point, the core problem of optical interferometry is transformed into the problem of solving the frequency and phase parameters in the superimposed cosine signal.

[0011] The designed multi-parameter estimation algorithm combining the disturbance quantization model and the root-finding method in interferometry mainly includes: preliminary step: interference pattern acquisition, step 1: multi-classification metric entropy model of optical interference signal, step 2: harmonic frequency correction method based on root-finding parameters, step 3: amplitude correction method, step 4: phase demodulation algorithm. The specific contents of each step are as follows:

[0012] The implementation process of the interferogram acquisition in the preliminary step is as follows: the measured object is placed in front of the interferometer with K = 4N PS and Δλ=λ0 2 / (2L min N PS ) to collect interference patterns, N PS The value is an integer not less than 20, L min is the optical thickness T0 of the measured element;

[0013] The implementation process of the multi-classification metric entropy model of the optical interference signal in step 1 is as follows: select 5 pixel points near the center of the interference pattern to average the interference light intensity as the analysis object and record it as ISPD , use the entropy calculation process formula to calculate I SPD Multi-classification metric entropy F E , if F E Greater than the entropy threshold F ET Then perform frequency and amplitude correction, the set F ET At least 0.7, this patent sets F ET is 0.75;

[0014] The implementation process of the harmonic frequency correction method based on root-finding parameters in step 2 is as follows: SPD Construct the Hankel matrix of the light intensity signal, calculate the unit root matrix of the Hankel matrix of the light intensity signal using the frequency solution process formula, and obtain the frequency of the surface harmonic signal of each device under test;

[0015] The implementation process of the amplitude correction method in step 3 is as follows: use the amplitude correction formula to calculate the amplitude of each harmonic and determine whether a frequency folding problem has occurred. If the amplitude height of the second harmonic peak from left to right in the spectrum is greater than the amplitude height of the third harmonic peak, it is determined that frequency folding has occurred. At this time, if frequency folding occurs, the corresponding harmonic is N PS / 2 is folded in the first spectrum period, which can be obtained by h p +N PS / 2 to obtain the correct harmonic frequency. If spectrum folding does not occur, the harmonic frequency obtained by the harmonic frequency correction method based on the root parameter is directly output for subsequent calculation;

[0016] The implementation process of the phase demodulation algorithm in step four is as follows: substitute the frequency values ​​obtained in steps two and three into the phase demodulation function design formula, use the phase demodulation formula to obtain the wrapped phase, unwrap and de-tilt the wrapped phase, and then input the wavefront height and phase distribution linear formula to obtain the wavefront of each surface, and finally obtain the surface shape distribution of each surface.

[0017] Compared with the prior art, the combined algorithm designed by the present invention has the following beneficial effects:

[0018] (1) The designed multi-classification metric entropy model can establish a correlation between the entropy value distribution and the degree of disturbance of the light intensity sequence to evaluate the size of the noise and error in the interference pattern acquisition process, and based on this, the necessity of parameter estimation correction can be determined; (2) When there is measurement disturbance, the designed root-finding algorithm can accurately obtain the frequency solution results of each harmonic, and the developed model can process the real-valued light intensity signal, overcoming the defect of the traditional method that relies on the complex signal model; (3) The designed amplitude correction algorithm can accurately obtain the amplitude distribution of each harmonic under disturbance conditions, and can obtain the true harmonic frequency distribution results; (4) The designed combined algorithm can simultaneously realize accurate frequency amplitude solution and efficient phase demodulation for the real-valued light intensity signal, and the surface shape distribution of each surface can be obtained through a set of interference patterns. BRIEF DESCRIPTION OF THE DRAWINGS

[0019] Figure 1 Schematic diagram of wavelength-tuned interferometry: (a) Internal components of the interferometer; (b) Multi-surface interference superposition process; (c) Spatial harmonic superposition; (d) Time domain harmonic superposition;

[0020] Figure 2 This is the implementation flow chart of the multi-parameter estimation algorithm that combines the disturbance quantization model with the root-finding method in the designed interferometry;

[0021] Figure 3 It is the simulated optical measurement model, the simulated surfaces, the interference pattern of each surface simulation, and the harmonic signal;

[0022] Figure 4 The superimposed noisy signal: (a) noisy signal in the time domain; (b) superimposed signal in the time domain; (c) superimposed interference pattern in the spatial domain;

[0023] Figure 5 is the distribution of entropy under different signal-to-noise ratios;

[0024] Figure 6 This is the result of detecting harmonic peaks using this algorithm and the FFT method;

[0025] Figure 7 is the amplitude and frequency solution error obtained by this algorithm and FFT method;

[0026] Figure 8 is the surface reconstruction error and wrapping phase obtained using this algorithm;

[0027] Figure 9 The frequency and amplitude solution errors of the FFT method under different phase shift error conditions are: (a) frequency; (b) amplitude;

[0028] Figure 10 The frequency and amplitude solution errors of this algorithm under different phase shift error conditions are: (a) frequency; (b) amplitude;

[0029] Figure 11 The solution errors under different signal-to-noise ratios are: (a) amplitude error; (b) surface error; (c) frequency error;

[0030] Figure 12 The solution errors under different signal-to-noise ratios are: (a) surface shape error; (b) frequency error;

[0031] Figure 13 is the interferometer used in the experiment;

[0032] Figure 14 Experimental data: (a) interference pattern; (b) interference light intensity used to solve the frequency and amplitude; (c) time domain curve of the window function used;

[0033] Figure 15 is the spectrum analysis result in the experiment;

[0034] Figure 16 The experimental results are: (a) the surface shape of each surface; (b) the wrapping phase;

[0035] Figure 17 These are the results of repeated measurements: (a) PV curve of surface shape; (b) RMS curve of surface shape; (c) PVQ curve of surface shape; (d) curve of harmonic amplitude; and (e) curve of harmonic frequency. DETAILED DESCRIPTION

[0036] The present invention will be further described below in conjunction with the accompanying drawings and preferred embodiments. The following embodiments are only used to more clearly illustrate the technical solutions of the present invention and are not intended to limit the scope of protection of the present invention. This section will further illustrate the technical solutions in the above invention content with reference to specific embodiments. The preferred embodiments of the present invention are described in detail as follows:

[0037] Example 1

[0038] In this embodiment, see Figures 1 and 2 The designed multi-parameter estimation algorithm combining the disturbance quantization model and the root-finding method in interferometry mainly includes: preliminary step: interference pattern acquisition, step one: multi-classification metric entropy model of optical interference signal, step two: harmonic frequency correction method based on root-finding parameters, step three: amplitude correction method, step four: phase demodulation algorithm, the algorithm implementation flow chart is as follows Figure 2 As shown. The implementation process of the designed algorithm can be organized as follows:

[0039] The implementation process of the interferogram acquisition in the preliminary step is as follows: the measured object is placed in front of the interferometer with K = 4N PS and Δλ=λ0 2 / (2L min N PS) to collect interference patterns, N PS The value is an integer not less than 20, L min is the optical thickness T0 of the measured element;

[0040] The implementation process of the multi-classification metric entropy model of the optical interference signal in step 1 is as follows: select 5 pixel points near the center of the interference pattern to average the interference light intensity as the analysis object and record it as I SPD , use the entropy calculation process formula (corresponding to formula (6) to formula (11)) to calculate I SPD Multi-classification metric entropy F E , if F E Greater than the entropy threshold F ET Then perform frequency and amplitude correction, the set F ET At least 0.7, this patent sets F ET is 0.75;

[0041] The implementation process of the harmonic frequency correction method based on root-finding parameters in step 2 is as follows: SPD Construct the Hankel matrix of the light intensity signal (corresponding to formula (12)), use the frequency solution process formula (corresponding to formulas (13) to (19), where formula (19) is the final step of solving the frequency) to calculate the unit root matrix of the Hankel matrix of the light intensity signal and obtain the frequency of the surface harmonic signal of each device under test;

[0042] The implementation process of the amplitude correction method in step 3 is as follows: use the amplitude correction formula (corresponding to formulas (24) to (28)) to calculate the amplitude of each harmonic and determine whether a frequency folding problem has occurred. If the amplitude height of the second harmonic peak from left to right in the spectrum is greater than the amplitude height of the third harmonic peak, it is determined that frequency folding has occurred. At this time, if frequency folding occurs, the corresponding harmonic is N PS / 2 is folded in the first spectrum period, which can be obtained by h p +N PS / 2 to obtain the correct harmonic frequency. If spectrum folding does not occur, the harmonic frequency obtained by the harmonic frequency correction method based on the root parameter is directly output for subsequent calculation;

[0043] The implementation process of the phase demodulation algorithm in step 4 is as follows: substitute the frequency values ​​obtained in steps 2 and 3 into the phase demodulation function design formula (corresponding to formula (5)), use the phase demodulation formula (corresponding to formula (4)) to obtain the wrapped phase, unwrap and de-tilt the wrapped phase, and then substitute it into the linear formula of wavefront height and phase distribution (corresponding to formula (3)) to obtain the wavefront of each surface, and finally obtain the surface shape distribution of each surface.

[0044] The technical and instrumental diagram of the wavelength tuning phase-shifting interferometry method is shown in the attached figure. Figure 1As shown. The typical process of this technology is to drive the light source to tune the wavelength while collecting a series of interference patterns, and then locate the frequency characteristics of each surface harmonic sub-signal from the interference pattern, and then demodulate the initial phase and wavefront height accordingly. Figure 1 Where L is the distance from the front surface (SF) of the DUT to the reference surface (SRF), which can also be called the measurement distance; T represents the average length (average thickness) between the front surface (SF) and the back surface (SR) of the DUT. In this case, its optical thickness T0 can be expressed as n z T(n z is the refractive index).

[0045] From a single pixel (x w ,y w ) The interference light intensity at different frame numbers is the observation object, and the light intensity of the time domain harmonic superposition I SP As shown in formula (1), k is the phase shift sequence (k = 0, 1, 2, ..., K-1, K is the total number of phase shift frames), p is the harmonic sequence (p = 1, 2, 3 corresponds to the thickness variation interference signal, the front surface interference signal and the back surface interference signal respectively), G p is the harmonic contrast related to the number of reflections on each surface, L p is the optical path difference of each harmonic (also called harmonic cavity length, as shown in the figure Figure 1 shown), V p is the wavefront height of the corresponding surface, λ0 is the starting tuning wavelength (λ0 is 632.67nm for the selected interferometer), Δλ is the single-step wavelength tuning amount, (x w ,y w ) are pixel coordinates.

[0046]

[0047] According to L p Adjusting Δλ is the key measure to determine the harmonic phase shift interval and harmonic frequency. The calculation model of single-step wavelength tuning is Δλ=λ0 2 / (2L min N PS ), where L min is the minimum optical path difference, N PS It is possible to determine that the phase shift interval of the harmonic with the shortest optical path difference is 2π / N PS The phase shift division coefficient.

[0048] Therefore, collecting the interference pattern is the first step of the measurement. After the piece to be measured is placed in front of the interferometer, K = 4N PS and Δλ=λ0 2 / (2L min N PS ) to collect interference patterns, N PSThe value is an integer not less than 20, L min is the optical thickness T0 of the measured element.

[0049] Furthermore, based on formula (1) and the calculation model of the single-step wavelength tuning amount, by Taylor expansion and ignoring the high-order terms, the expression of each harmonic frequency can be obtained as shown in formula (2). Q in formula (2) A and Q B They are the number of times the laser beam is reflected between the front surface of the DUT and the reference mirror, and the number of times it is reflected inside the DUT.

[0050]

[0051] Furthermore, the wavefront height of each surface can be determined by the initial phase The linear relationship between the two is shown in formula (3).

[0052]

[0053] One of the effective methods to obtain the initial phase is to use Fourier transform technology to demodulate the light intensity signal based on the precise harmonic frequency obtained, thereby obtaining the initial phase as shown in formula (4): After unwrapping and de-tilting the initial phase of the wrapped form, the accurate surface shape distribution of each surface can be obtained using formula (3).

[0054]

[0055] The harmonic phase demodulation functions Ap and Bp in formula (4) can be designed using formula (5). W in formula (5) R (k) is the selected window function, which helps to improve the phase demodulation accuracy.

[0056]

[0057] The above calculation process is the calculation step of the phase demodulation part.

[0058] For signals with low disturbance levels and near-pure targets, accurate frequency estimation can be achieved using the classic FFT method. However, under strong disturbance conditions, the inherent spectral leakage and picket fence effects of the FFT method can significantly reduce parameter estimation accuracy. Therefore, in interferometry, especially when faced with the problem of multi-surface interference superposition in measuring DUTs with parallel front and back surfaces, it is necessary to construct a quantitative model for the disturbance of the measurement and acquisition process signal.

[0059] The complexity of a signal sequence can be characterized by a statistical model called entropy, which increases significantly under the influence of measurement disturbances or sampling errors. Given the high complexity of measurement processes in real-world environments, which often involve multiple noise sources, the boundary between measurement error and sampling error can be ambiguous. Therefore, the introduction of a multi-class metric function is considered a preferred option, aiming to enable the constructed quantitative model to be applicable to the assessment of the degree of disturbance of signal sequences under complex measurement conditions.

[0060] When the signal is affected by sampling errors, the complexity of the sequence will increase. Therefore, using the entropy model to evaluate the degree of interference of the sequence is an effective method, which helps to decide whether to modify the key harmonic parameters (because parameter modification requires additional calculations). ae )(1≤i ae ≤K), given the pattern dimension m ae (usually 2 or 3), the mode vector U(i ae ) can be expressed as:

[0061] U(i ae )=[u(i ae ),u(i ae +1),…,u(i ae +m ae -1)]-u0(i ae ), (6)

[0062] In formula (6) And i ae ≠j ae .

[0063] definition is the mode vector and U(i ae ) and U(j ae ), expressed as formula (7):

[0064]

[0065] Traditional entropy theory relies on the Heaviside function to construct statistical models, resulting in the generated entropy function having obvious step characteristics and poor continuity. To overcome this limitation, the multi-class metric entropy theory is designed to use an exponential multi-class metric function to evaluate the similarity between patterns, thereby achieving continuous multi-level quantization of similarity measurement, as shown in Formula (8).

[0066]

[0067] In formula (8), r = g × S td[U], Std is the standard deviation of U, g is the similarity tolerance coefficient, r is between 0.1 and 0.25Std(U); n is the gradient of the multi-classification metric function, and the function shown in formula (9) is defined:

[0068]

[0069] Therefore, the multi-classification metric entropy function can be defined as formula (10):

[0070]

[0071] For a finite time series, the multi-classification metric entropy can be expressed as formula (11):

[0072]

[0073] The interference light intensity after averaging the five pixel points near the center of the interference pattern is selected as the analysis object and recorded as I SPD , substitute into formula (6) to formula (11) to calculate I SPD Multi-classification metric entropy F E , if F E Greater than the entropy threshold F ET Then perform frequency and amplitude correction, the set F ET At least 0.7, this patent sets F ET is 0.75.

[0074] The above are the steps for constructing a multi-classification metric entropy model for optical interference signals.

[0075] The accuracy of the harmonic signal frequency has a decisive influence on phase demodulation. The accurate harmonic frequency h is calculated. p Only then can the accurate initial phase distribution be obtained. SPD As the analysis object, the observation vector matrix X of the signal can be constructed from the K-point sampling data in the form of formula (12):

[0076]

[0077] When constructing the model shown in Equation (12), m > 2P, where the number of signal sources is P. For multi-surface interference signals, P = 3 when undersampling is not present. The light intensity data in the harmonic signal analysis model shown in Equation (12) is constructed as a Hankel matrix of real-valued data, requiring less data than traditional algorithms based on complex models.

[0078] According to the root parameter estimation principle, the covariance matrix Rx of the observation vector matrix X is obtained. Rx can be expressed as the outer product expectation function of X, that is, the form shown in formula (13):

[0079] R x =E{XX T} . (13)

[0080] Performing eigenvalue decomposition on the matrix Rx, we can obtain formula (14):

[0081]

[0082] The matrices in formula (14) can be expressed as:

[0083] U N =(S1,…,S m-2p ), (15a)

[0084] U S =(S m-2p+1 ,…,S m ), (15b)

[0085] Σ N =diag(λ1,…,λ m-2p ), (15c)

[0086] Σ S =diag(λ m-2p+1 ,…,λ m ), (15d)

[0087] Among them, λ l and S l R x The lth eigenvalue and its corresponding eigenvector in , and λ1≈λ2≈…≈λ m-2P ≤λ m-2P+1 ≤…≤λ m , the number of signal sources can be determined according to the size of the eigenvalues. If the first three eigenvalues ​​are significantly larger than the following eigenvalues, it can be determined that the current number of sources is 3, that is, there are 3 main signals in the analysis object; U N By R x The subspace spanned by the eigenvectors corresponding to the first m-2P small eigenvalues ​​of , that is, the noise subspace; U S By R x The subspace spanned by the eigenvectors corresponding to the last 2P largest eigenvalues ​​of is also the signal subspace. As can be seen, the larger m is, the higher the dimension of the noise space is, the higher the parameter estimation accuracy is, but the computational complexity is also higher.

[0088] In order to extract spatial parameter information using noise feature vectors, a root-finding parameter estimation polynomial is constructed:

[0089]

[0090] In formula (16), q(z)=[1, z,…, z m-1 ] T ,The superscript H indicates conjugate transpose.

[0091] The frequency of the signal can be estimated based on the zero point of formula (16). However, considering that formula (16) also contains the power term of z*, this makes the zero-finding process complicated. Therefore, formula (16) can be corrected to the form of formula (17):

[0092]

[0093] Formula (17) is the root parameter polynomial. It can be seen that f(z) is a 2(m-1) degree polynomial. It has (m-1) pairs of mutually conjugate roots that are mirror-symmetric about the unit circle, and the 4P roots with the largest amplitude are exactly distributed on the unit circle (P pairs of conjugate double roots). The number of roots is linearly corresponding to the P frequencies of the collected signal. Therefore, the root matrix R on the unit circle is f It can be expressed as formula (18):

[0094] R f =[z1,z2,…,z 4P ] T . (18)

[0095] Then the estimated harmonic frequency of the sampling signal is given by formula (19):

[0096] h p =arg(z 4i-3 ) / 2π . (19)

[0097] The arg function in formula (19) refers to the phase of the complex number. The above is a harmonic frequency correction method based on root parameters. The harmonic frequency can be solved using formula (19).

[0098] Harmonic amplitude can assist in determining the correspondence between harmonic types and spectrum peaks in interferometry. When the signal sequence is disturbed, spectrum leakage will cause the frequency energy that should be concentrated at the peak position to spread to adjacent frequency points, so the amplitude needs to be corrected. This patent takes the Hanning window as an example and introduces an amplitude solution algorithm based on ratio correction. The normalized Hanning window is defined in the form of formula (20):

[0099]

[0100] The spectral function of formula (20) can be written as formula (21):

[0101]

[0102] In formula (21), a=0.5, and the normalized frequency is set to The simplified form of the spectrum modulus function formula (22) can be obtained:

[0103]

[0104] In formula (22), when f 1 →0, W1(f 1 )→a; when f 1 →±1, W1(f 1 )→(1-a) / 2. Let Then formula (22) can be written in the form of formula (23).

[0105]

[0106] Therefore, the correction function is constructed as shown in formula (24):

[0107]

[0108] y in formula (24) k and y k+1 is the harmonic amplitude of the corresponding peak in the FFT discrete spectrum, that is, the peak height before correction. Since the Hanning window a = 0.5, then c → ∞, then formula (24) can be written as formula (25):

[0109]

[0110] Solve The inverse function of can be obtained later as formula (26):

[0111]

[0112] Then we can get formula (27):

[0113]

[0114] From Equation (26), we can see that the main lobe function of the Hanning window has the following properties: Take any two points on the curve and The center of gravity of these two points is at the origin of the coordinate system. When the origin corresponds to the amplitude spectrum, the frequency at the center of gravity is the true frequency of the harmonic signal window normalized by the length of the harmonic signal. This can be called the center of gravity theorem for discrete spectrum correction by adding Hanning window truncation to the harmonic signal. The Hanning window amplitude correction model can then be obtained as shown in formula (28):

[0115]

[0116] Perform the above operations on the frequency points near the spectrum peaks corresponding to the three harmonics to obtain the harmonic amplitude G p (G in the above formula is the amplitude of a single harmonic. After calculating the amplitudes of the three harmonics respectively, they can be marked as G p , p=1,2,3).

[0117] If the amplitude of the second harmonic peak from left to right is found to be greater than the amplitude of the third harmonic peak in the spectrum, it is determined that frequency folding has occurred, and the true frequency distribution of the second and third harmonics can be obtained by using half of the sampling frequency as the symmetry axis.

[0118] It should be noted that in the combined algorithm designed in this patent, the window functions used in the amplitude solution and the phase solution can be different. For this patent, the amplitude solution uses a Hanning window, and the phase solution uses a Blackman window.

[0119] Example 2:

[0120] In this embodiment, see Figures 3 to 17 This section will quantitatively analyze the performance of the multi-parameter estimation algorithm combined with the perturbation quantization model and the root-finding method in the interference measurement designed in this patent based on the multi-surface interference simulation model shown in formula (29), and conduct experimental verification to analyze the effectiveness of this algorithm in the actual measurement environment. First, construct the Zernike polynomial to simulate the reflected wavefront V' of each surface. s , in the multi-surface interference simulation model, V' s The subscripts s=A, B represent the measured front and back surfaces respectively, and the reference mirror is set to an absolute plane.

[0121]

[0122] In formula (29), k0~k4, (x0, y0) and R W The contour coefficient and distribution coefficient of the wavefront are shown in Table 1. The interference pattern size is set to 500 pixels × 500 pixels.

[0123] Table 1. Simulation parameters

[0124]

[0125] N PS =20, K=4N, the signal-to-noise ratio (SNR) is set to 50dB, and a 10% wavelength tuning error is introduced to make the simulation closer to the actual conditions, n z =1.5, T=10mm, L=46.875mm (simulating non-integer multiple measurement positions).

[0126] exist Figure 3 In the figure, the light intensity of the interference pattern is normalized between -1 and 1, and its unit is no longer displayed (the interference signal below is also processed in the same way). Figure 2 The simulation models of each surface can obtain the interference results of time domain and space domain superposition as shown in the attached figure. Figure 4 shown.

[0127] In order to verify the effectiveness of the entropy model in quantifying the degree of interference of the sequence, based on the simulation conditions given in Example 1, the intensity sequence of the five pixels in the center of the interference pattern after removing the light intensity is selected as the analysis object (the entropy, frequency and amplitude analysis in this patent are all processed in this way), and the entropy under different signal-to-noise ratios is calculated. The results are shown in the accompanying figure. Figure 5 As shown. Figure 5 From the results shown, we can see that the entropy is larger when the signal-to-noise ratio is smaller, which also verifies the possibility of the entropy model to quantify the errors in the sequence.

[0128] Then, the frequency and amplitude performance are comprehensively analyzed. First, based on the simulated interference model, using FFT as the comparison algorithm, the spectrum results can be obtained as shown in the attached figure. Figure 6 As shown in the figure, harmonics 1 to 3 refer to the harmonic signals with increasing frequency distribution in the light intensity signal corresponding to each surface. Therefore, harmonic 1 corresponds to the thickness change signal, and harmonics 2 and 3 correspond to the front surface signal and the back surface signal, respectively.

[0129] Furthermore, the absolute and relative errors of harmonic frequencies and amplitudes can be extracted to obtain Figure 7 The results shown.

[0130] from Figure 7 The results shown in the figure show that the accuracy of this algorithm is 1 to 2 orders of magnitude higher than that of the FFT method in terms of frequency and amplitude estimation errors. Based on the obtained harmonic frequencies, the phase of each surface can be demodulated and the wave surface can be reconstructed. The results are shown in the figure. Figure 8 shown. Figure 8 The surface reconstruction error shown demonstrates that this algorithm can achieve sub-nanometer theoretical accuracy and accurately separate and demodulate surface shape information. Excluding the drawing function, the average runtime for 10 calculations was 0.73 seconds (CPU-AMD4800H, 16GB RAM).

[0131] Furthermore, in order to verify the performance of the multi-parameter estimation algorithm combining the disturbance quantization model and the root-finding method in the designed interferometry and its improvement over the existing methods, L E As the linear phase shift error, we can get Figure 9 and Figure 10 The result of Figures 9 and 10It can be seen from the results that this algorithm can still obtain relatively accurate solution results under the condition of phase shift error.

[0132] At the same time, in order to verify the robustness of the multi-parameter estimation algorithm combining the disturbance quantization model and the root-finding method in the designed interferometry under Gaussian noise conditions, the error distribution of frequency, amplitude and surface shape (the root mean square of the error surface, that is, the RMS value, is extracted as an indicator) is analyzed under different signal-to-noise ratios. The results are shown in the attached figure. Figure 11 As shown. As a comparison object, the double iterative least squares algorithm (doi:10.1016 / j.optlaseng.2018.11.008) is selected as the comparison object, and the solution result can be obtained as shown in the figure Figure 12 However, it should be noted that this method cannot obtain the amplitude parameters. Figures 11 and 12 It can be seen from the results that the multi-parameter estimation algorithm combining the disturbance quantization model and the root-finding method in the designed interferometry measurement has good robustness to additive Gaussian noise, and the effectiveness of this method is further proved.

[0133] In order to verify the effectiveness of the multi-parameter estimation algorithm combining the perturbation quantization model and the root-finding method in the designed interferometry measurement under experimental conditions, a Fizeau interferometer equipped with a wavelength-tunable laser (TLB-6804, NEW FOCUS, with a total tuning linewidth of 200 GHz, i.e., 0.26 nm) was used to collect the interferogram of a transparent parallel flat mirror with an average thickness of 30 mm, as shown in the attached figure. Figure 13 Random disturbances were minimized during the measurement, and the ambient temperature during interferogram acquisition was controlled at 25 ± 0.3°C.

[0134] In the acquisition plan, N PS =14, K=56, Δλ=3.2728e-04nm (equivalent tuning linewidth is 0.2346GHz). The test aperture is 100mm×100mm. To avoid the influence of invalid interference area on the measurement, the interference pattern is cropped from 1200 pixels×1200 pixels to 400 pixels×400 pixels. The data involved in the experiment are shown in the attached figure. Figure 14 Taking the selected light intensity as the object, the algorithm can be used to obtain the corrected frequency and amplitude as shown in the attached figure. Figure 14 shown.

[0135] Then, as shown in the accompanying drawings Figure 16 As shown in the figure, the phase demodulation function is constructed using the obtained frequency, which can realize the phase demodulation and surface reconstruction of each surface. Figure 15The results show that the light intensity signal in the experiment has spectrum leakage problem, that is, the peak energy is scattered in the adjacent frequency points, and the multi-parameter estimation algorithm combining the perturbation quantization model and the root-finding method in the designed interferometry can accurately solve the frequency and amplitude. Figure 16 The results show a clear distribution of phase and surface shape, demonstrating that the algorithm can separate and extract information from each surface without residual harmonics. The computation time is 0.32 seconds, completing the frequency and amplitude calculations and surface reconstruction for each surface.

[0136] In order to verify the repeatability of this algorithm, repeated measurements were performed under the same conditions. The peak-to-valley value (PV), root mean square, and PVQ (the peak-to-valley value obtained after removing the 5% maximum and minimum peaks) of the reconstructed surface shape are shown in the attached figure. Figure 17 As shown in the figure. At the same time, the frequency and amplitude of the solution are also analyzed. Figure 17 As can be seen in the figure, the designed algorithm demonstrates good performance in the experiment due to low repeated measurement errors. The graph shows that the maximum PV, RMS, and PVQ errors in repeated measurements are 4.45nm, 0.49nm, and 2.41nm, respectively, while the errors in the amplitude and frequency solutions are both on the order of 10e-3, demonstrating accurate and stable measurement performance.

[0137] The above describes the embodiments of the present invention in conjunction with the accompanying drawings, but the present invention is not limited to the above examples. Various changes can be made according to the purpose of the invention. Any changes, modifications, substitutions, combinations or simplifications made according to the spirit and principles of the technical solution of the present invention should be equivalent replacement methods. As long as they comply with the purpose of the invention and do not deviate from the technical principles and inventive concepts of the present invention, they belong to the scope of protection of the present invention.

Claims

1. The multi-parameter estimation algorithm designed by the present invention, which combines a perturbation quantization model with a root-finding method in interferometry, mainly includes: Preliminary steps: interference pattern acquisition, step one: multi-classification metric entropy model of optical interference signal, step two: harmonic frequency correction method based on root parameters, step three: amplitude correction method, step four: phase demodulation algorithm. The implementation process of the designed algorithm can be summarized as follows: The interferogram acquisition process of the preliminary step is as follows: the test piece is placed in front of the interferometer and K = 4N PS and Δλ=λ0 2 / (2L min N PS ) to collect interference patterns, N PS The value is an integer not less than 20, L min is the optical thickness T0 of the measured element; The implementation process of the multi-classification metric entropy model of the optical interference signal in step 1 is as follows: select 5 pixel points near the center of the interference pattern to average the interference light intensity as the analysis object and record it as I SPD , use the entropy calculation process formula to calculate I SPD Multi-classification metric entropy F E , if F E Greater than the entropy threshold F ET Then perform frequency and amplitude correction, the set F ET At least 0.7, this patent sets F ET is 0.75; The specific algorithm calculation process involved in this step is: For the K-point time series u(i ae )(1≤i ae ≤K), given the pattern dimension m ae Later, the mode vector U(i ae ) can be expressed as: U(i ae )=[u(i ae ),u(i ae +1),…,u(i ae +m ae -1)]-u0(i ae ) (6) In formula (6) i ae ,j ae =1,2,…,Km ae +1, and i ae ≠j ae ;definition is the mode vector and U(i ae ) and U(j ae ), expressed as formula (7): The exponential multi-classification metric function is used to quantify the similarity metric continuously at multiple levels, as shown in formula (8): In formula (8), r = g × Std[U], Std is the standard deviation of U, g is the similarity tolerance coefficient, and r is between 0.1 and 0.25Std(U); n is the gradient of the multi-classification metric function, which further defines formula (9): Therefore, the multi-classification metric entropy function F of a finite time series is E It can be defined as formula (11): Using formula (11), we can get I SPD Multi-classification metric entropy F E ; The implementation process of the harmonic frequency correction method based on root-finding parameters in step 2 is as follows: SPD Construct the Hankel matrix of the light intensity signal and use the frequency solution process formula to calculate the light intensity signal I SPD The unit root matrix of the Hankel matrix and the frequency of the surface harmonic signal of each DUT are obtained; The specific algorithm calculation process involved in this step is: Will I SPD The observation vector matrix X of the constructed signal is in the form of formula (12): According to the root parameter estimation principle, find the covariance matrix R of the observation vector matrix X x , R x It can be expressed as the outer product expectation function of X, R x That is, the form shown in formula (13): R x =E{XX T } (13) Performing eigenvalue decomposition on the matrix Rx, we can obtain formula (14): The matrices in formula (14) can be expressed as: U N =(S1,…,S m-2p ) (15a) U S =(S m-2p+1 ,…,S m ) (15b) S N =diag(λ1,…,λ m-2p ) (15c) S S =diag(λ m-2p+1 ,…,l m ) (15d) The root parameter estimation polynomial can be further constructed as shown in formula (16): In formula (16), q(z)=[1, z,…, z m-1 ] T ,Superscript H represents conjugate transpose; Considering that there is a power term of z* in formula (16), formula (16) can be corrected to the form of formula (17): The root matrix R on the unit circle f It can be expressed as formula (18): R f =[z1,z2,…,z 4P ] T (18) Then the harmonic frequency is formula (19): h p =arg(z 4i-3 ) / 2π (19) The arg function in formula (19) refers to the phase of the complex number, and the harmonic frequency can be solved using formula (19); The implementation process of the amplitude correction method in step 3 is: use the amplitude correction formula to calculate the amplitude of each harmonic and determine whether the frequency folding problem occurs. If the frequency folding problem occurs, the h p +N PS / 2 to obtain the correct harmonic frequency. If spectrum folding does not occur, the harmonic frequency obtained by the harmonic frequency correction method based on the root parameter is directly output for subsequent calculation; The specific algorithm calculation process involved in this step is: The correction function is constructed as shown in formula (24): y in formula (24) k and y k+1 is the harmonic amplitude of the corresponding peak in the FFT discrete spectrum, that is, the frequency peak height before correction. Then, formula (24) can be written as formula (25): Solve The inverse function of can be obtained later as formula (26): Then we can get formula (27): Therefore, the Hanning window amplitude correction model can be obtained as shown in formula (28): Perform the above operations on the frequency points near the spectrum peaks corresponding to the three harmonics to obtain the harmonic amplitude G p ; The implementation process of the phase demodulation algorithm in step 4 is as follows: the frequency values ​​obtained in steps 2 and 3 are substituted into the phase demodulation function design formula, the wrapped phase is obtained using the phase demodulation formula, the wrapped phase is unwrapped and de-tilted, and then the wavefront height and phase distribution linear formula is input to obtain the wavefront of each surface, and finally the surface shape distribution of each surface is obtained; The specific algorithm calculation process involved in this step is: The wavefront height of each surface can be determined by the initial phase As shown in formula (3): Using formula (4) to get the initial phase After unwrapping and de-tilting the initial phase of the wrapped form, the accurate surface shape distribution V of each surface can be obtained using formula (3) p (x w ,y w ): The harmonic phase demodulation function A in formula (4) p and B p The design can be done using formula (5), where W R (k) is the selected window function: After step 4 is completed, the surface shape distribution of each surface can be obtained.

2. The feature of step 2 of the multi-parameter estimation algorithm combining the disturbance quantization model and the root-finding method in the interference measurement according to claim 1 is that: when constructing the model shown in formula (12), m>2P is required, where the number of signal sources is P, and when no undersampling occurs, P=3 for the multi-surface interference signal; the light intensity data in the harmonic signal analysis model shown in formula (12) are all constructed as Hankel matrices of real-valued data, which requires less data than the traditional algorithm based on the complex model.

3. The characteristic of step 2 of the multi-parameter estimation algorithm combining the perturbation quantization model and the root-finding method in interferometry according to claim 1 is: l and S l R x The lth eigenvalue and its corresponding eigenvector in , and λ1≈λ2≈…≈λ m-2P ≤λ m-2P+1 ≤…≤λm, the number of signal sources can be determined based on the size of the eigenvalues. If the first three eigenvalues ​​are significantly larger than the following eigenvalues, it can be determined that the current number of sources is 3, that is, there are 3 main signals in the analysis object; U N By R x The subspace spanned by the eigenvectors corresponding to the first m-2P small eigenvalues ​​of , that is, the noise subspace; U S By R x The subspace spanned by the eigenvectors corresponding to the last 2P large eigenvalues ​​is the signal subspace.

4. The feature of step three of the multi-parameter estimation algorithm combining the perturbation quantization model and the root-finding method in the interference measurement according to claim 1 is: if the amplitude height of the second harmonic peak from left to right in the spectrum is found to be greater than the amplitude height of the third harmonic peak, it is determined that frequency folding has occurred, and the true frequency distribution of the second harmonic and the third harmonic can be obtained by taking half of the sampling frequency as the symmetry axis. In the algorithm designed in this patent, the Hanning window is used in the amplitude solution part, and the Blackman window is used in the phase solution part.

Images