阅读说明：本技术 一种基于相移特征多项式高精度拟合的多表面测量方法 (Multi-surface measurement method based on phase shift characteristic polynomial high-precision fitting ) 是由 常林 王陈 闫恪涛 郑维伟 徐瞿磊 孙涛 于瀛洁 于 2020-06-28 设计创作，主要内容包括：本发明涉及一种基于相移特征多项式高精度拟合的多表面测量方法,通过多加权多步相移算法中的特征多项式进行计算,得到其自变量的系数分布,为了测量时应用的便捷性,对该系数分布进行不同参数下的拟合,并且给出在不同的腔长系数下的可用相移参数的算法适用范围,实现测量方案的制定与算法设计的统一。此外还能够根据测量人员的实际需求和主观意愿进行测量方案中窗函数的选择,提供更多的测量方案制定的可行性,节约测量成本,能够实现一次性的多表面非接触式测量。本发明方法最大程度减少必须采集帧数,简化测量过程,并且尽可能避免纳入误差,并减少应用和计算成本。(The invention relates to a multi-surface measurement method based on phase shift characteristic polynomial high-precision fitting, which is characterized in that the coefficient distribution of independent variables of the multi-weighted multi-step phase shift algorithm is obtained by calculating the characteristic polynomial, the coefficient distribution is fitted under different parameters for the convenience of application during measurement, the algorithm application range of available phase shift parameters under different cavity length coefficients is given, and the formulation of a measurement scheme and the unification of algorithm design are realized. In addition, the selection of a window function in a measurement scheme can be carried out according to the actual requirements and subjective wishes of measurement personnel, more feasibility for making the measurement scheme is provided, the measurement cost is saved, and one-time multi-surface non-contact measurement can be realized. The method of the invention reduces the number of frames to be acquired to the maximum extent, simplifies the measurement process, avoids the inclusion of errors as far as possible and reduces the application and calculation cost.)

1. A multi-surface measurement method based on phase shift characteristic polynomial high-precision fitting is characterized by comprising the following steps:

(1) calculating through a characteristic polynomial in a multi-weighting multistep phase shift algorithm to obtain coefficient distribution of independent variables of the characteristic polynomial, namely ideal window function distribution, and expanding coefficients from the polynomial, so that the coefficient distribution cannot be directly utilized in engineering;

(2) for convenience of application in measurement, fitting coefficient distribution of independent variables of the characteristic polynomial under different parameters, and analyzing distribution of residual errors solved by an algorithm of available phase shift parameters N under different cavity length coefficients M;

(3) setting algorithm parameters by a numerical analysis method, obtaining the range of the position of the detected piece corresponding to different algorithm parameters, digitizing the range of the position of the detected piece based on the algorithm characteristics and the numerical analysis result, and calculating the general formula of the cavity length coefficient corresponding to different phase shift parameters, so as to be directly applied during measurement and realize the unification of the measurement scheme and the algorithm design;

(4) according to the actual requirements and subjective wishes of measuring personnel, the window function in the measuring scheme is selected, more feasibility for making the measuring scheme is provided, and the measuring cost is saved, so that the one-time multi-surface non-contact measurement is customized and completed.

2. The method for multi-surface measurement based on high-precision fitting of phase shift feature polynomials as claimed in claim 1, characterized in that: in the step (2), gaussian fitting, fourier fitting and polynomial fitting are performed on the characteristic polynomial through the distribution of the independent variable coefficients of the characteristic polynomial, and three fitting numerical results are provided for selection.

3. The method for multi-surface measurement based on high-precision fitting of phase shift feature polynomials as claimed in claim 2, characterized in that: in the step (2), the three fitted window functions are calculated through fitting, and the general formula is expressed as:

where k is the number of phase shift ordinal,for sampling window functions using a polynomial fitting method,

4. The method for multi-surface measurement based on high-precision fitting of phase shift feature polynomials as claimed in claim 2, characterized in that: in the step (2), when polynomial fitting is performed, a phase shift ordinal number, that is, an independent variable coefficient distribution ordinal number, is used as a variable of a fitting value, and a corresponding current polynomial distribution coefficient is used as a fitting target value; sampling window function based on polynomial fitting

sampling window function based on polynomial fittingIs expressed as:

wherein a is₀～a₅Is a parameter to be solved; the minimum Q of the squared error is found and expressed as:

wherein Z is the total acquisition frame number of the interferogram; by using the principle of least square, solving the square value of the above formula and making the square value tend to 0, thus obtaining the distribution of each coefficient; specifically, the partial derivative of each coefficient value is calculated to be 0, so as to obtain the optimal solution of the extreme value:

by solving the above formula, the equation expression of least square solution of each coefficient can be obtained respectively, and the least square solution can be obtained by solving the equation to complete parameter fitting.

5. The polynomial of claim 2 based on phase shift characteristicsThe high-precision fitting multi-surface measurement method is characterized by comprising the following steps of: in the step (2), when fourier fitting is performed, the distribution of the independent variable coefficients of the characteristic polynomial is also used as a fitting ideal value, and a sampling window function based on the fourier fitting method is used

wherein b is₀,b₁,b₂,

wherein

6. The method for multi-surface measurement based on high-precision fitting of phase shift feature polynomials as claimed in claim 2, characterized in that: in the step (2), the Gaussian fitting form is constructed by analyzing the distribution of the target window function during the Gaussian fittingThe expression is as follows:

wherein c is_m,

On the circle near each step iteration point, the current point is used to approximate the objective functionSecond order approximation function of numbers

WhereinAndare the first and second order gradients of the objective function; h () represents the Hessian matrix of the acquisition object;

J^T(k) is the transposition of the Jacobian matrix of the objective function; first, consider thatFor the minimum solution, p_kThe fixed direction of (a) is the gradient descent direction, along which the minimum in the trust domain is found; from the initial point, the iteration vector can be represented as:

wherein Δ_kThe radius of the confidence domain can be adjusted by the ratio of the descending value of the quadratic approximation function to the descending value of the objective function; if the ratio is close to 1, the radius is increased; if the ratio is equal to 1, the radius remains unchanged; if the ratio is less than 1, the radius is reduced; tau is_kFrom the following eachItem determination:

the initial point may be based on experience or set to a random number between 0-1; in order to ensure the iteration precision and the result reliability under the maximum iteration step length, the iteration parameters are set as follows: the maximum and minimum changes of the finite difference gradient coefficients are set to 0.1 and 10, respectively^-8(ii) a The maximum evaluation times and the allowed iteration times of the iteration model are 600 and 400 respectively; minimum iteration error of 10^-6(ii) a And if the maximum iteration number or the minimum iteration error is met, stopping iteration and outputting a current value.

7. The method for multi-surface measurement based on high-precision fitting of phase shift feature polynomials as claimed in claim 1, characterized in that: in the step (3), a quantitative analysis result under the condition of available cavity length coefficients under different phase shift parameters N is summarized according to a comprehensive analysis of error distribution results by using a numerical analysis method, and is represented by using a general formula:

M：

comprehensively considering the cavity length position and the wavelength minimum tuning range of the measured piece, taking N as an integer of 8-14, and taking r as a non-negative integer of 0,1,2 and 3 … …; when the available N values are addressed at the cavity length factor M, the r values increase from 0 and N increases from 8, addressing is done stepwise.

Technical Field

The invention relates to a multi-surface measuring method based on phase shift characteristic polynomial high-precision fitting, in particular to a calculation method for realizing convenient representation of a window function through high-precision fitting of an independent variable coefficient of a characteristic polynomial of a phase shift algorithm.

Background

The high-precision parallel flat plate with high surface quality has important significance in the design and construction process of an optical system, so that the accurate measurement of the surface morphology of the high-precision parallel flat plate has high application value. How to realize the non-contact simultaneous measurement of a multi-surface measured piece by a modern measurement method and the unification of algorithm design and measurement schemes are important problems faced by researchers in the field.

The traditional measuring method of the multi-surface measured piece comprises the following steps: the method comprises the steps of collecting an interferogram of a measured piece through a hardware phase-shifting interferometer, coating Vaseline or delustering paint on the surface of the interferogram of each surface of the measured piece in sequence in the process of collecting the interferogram of each surface of the measured piece to inhibit reflection signals of multiple surfaces, processing the interferogram of each surface obtained in sequence, and obtaining an initial phase of the interferogram through an algorithm so as to obtain the surface appearance. The drawbacks of this technique are however very evident:

1. non-contact measurement cannot be realized, which means that the high-precision surface of a measured piece is easily damaged in the process of smearing and cleaning the extinction material;

2. the measurement cost is high, because the extinction material needs to be cleaned and smeared for multiple times, the measurement time needs to be long, multiple times of algorithm processing are needed, the appearance results of multiple surfaces of the measured piece are obtained in sequence, and the calculation cost is also high;

3. errors are easy to be brought in, and the traditional hardware phase-shift interferometer is used, so that the phase-shift precision in the measurement process is limited, and the hardware errors and the hysteresis errors are easy to be brought in; in addition, because the extinction material needs to be coated and cleaned for many times, and the tested piece needs to be clamped and taken down for many times, the clamping position and the clamping inclination of each time cannot be strictly guaranteed to be the same, and a larger error is brought in;

4. the thickness cannot be measured, and because the traditional interference pattern acquisition mode is that each surface of the front surface and the back surface is acquired once, the thickness information of the measured piece cannot be obtained, and the specific numerical value of the thickness change signal cannot be obtained once through the algorithm.

In recent years, measurement methods based on wavelength phase shift interferometry have been developed. The main characteristics are as follows: the wavelength is tuned through the wavelength tunable laser, so that phase shift is realized, and a related dephasing algorithm is designed based on the difference of phase shift frequencies of all surfaces of the tested piece, so that dephasing and measurement of the multi-surface tested piece are realized.

There are two main types of mainstream algorithms used in the current field:

one method is a multi-surface measurement method based on a least square algorithm, and the method can accurately demodulate the phase of the multi-surface of the measured piece under the condition of small error. The disadvantages of this algorithm are however also very significant, namely: measurement under the condition of large error cannot be realized; the phase shift frequency needs to be accurately evaluated, otherwise, the measurement result cannot be demodulated; the initial conditions need to be accurately given, otherwise the measurement results are very inaccurate.

Another method is a multi-step weighted phase shift algorithm, which can obtain the initial phase values of the surfaces of the measured object simultaneously in a simpler weighting calculation mode. Meanwhile, the disadvantage of this algorithm is that: the weighted multi-step phase shift algorithm is developed based on discrete fourier transform, so that the distribution of the phase shift frequency of the target signal has higher requirements, namely that multi-surface phase demodulation of the measured piece at any measuring position cannot be realized. The current weighting multi-step algorithm does not give out a specific measurement range of the algorithm, in other words, no numerical measurement position is available as a guide for formulating a measurement scheme, so that the distribution of numerical measured piece positions is not generated in the actual measurement process, the algorithm is easy to fail, and demodulation cannot be performed. In addition, the existing weighted multi-step algorithm does not numerically give the applicable range of the algorithm and the general formula thereof, and is very inconvenient in actual measurement. Another disadvantage is that the algorithm determines the specific value of the sampling window function based on the coefficient expansion of each order variable of the characteristic polynomial, but the current algorithm processes this step by using the coefficient expansion of the characteristic polynomial and then performing a piecewise polynomial fitting on the characteristic polynomial to obtain the general formula of the window function. But the segmented polynomial fit has the disadvantages that: the calculation cost required by the piecewise fitting is high, and more operation operations are required; secondly, since the fitting accuracy of polynomial fitting is not high, measurement errors occur in nature, and this is a technical problem to be solved urgently.

Disclosure of Invention

In order to solve the technical problems, the invention provides a multi-surface measuring method based on phase shift characteristic polynomial high-precision fitting, which utilizes a multi-step weighting algorithm to carry out measurement and optimization of multi-surface separation, aims to solve the defects in the existing measuring method, particularly improves the fitting precision of a window function in the weighting multi-step phase shift algorithm, and can be represented by a general formula (general formula) of the position of a measured piece. The measured piece aimed at by the invention is a multi-surface transparent flat plate measured piece, and the acquired interference pattern is an aliasing light intensity distribution diagram formed by mutual interference among all surfaces, so that the phase resolution and measurement cannot be directly carried out.

The principle of the invention is as follows: the coefficient distribution of the independent variable is obtained by calculating the characteristic polynomial in the multi-weighting multistep phase shift algorithm, the coefficient distribution is fitted under different parameters for the convenience of application in measurement, the application range of the algorithm under different measurement positions is given, and the unification of the formulation of the measurement scheme and the algorithm design is realized. The measuring method provided by the invention can save the measuring cost and realize one-time multi-surface non-contact measurement. For convenience of expression and calculation, the cavity length coefficient M is defined as: the actual distance/(the thickness of the measured piece x the refractive index of the measured piece material) between the front surface of the measured piece and the reference mirror; defining the phase shift parameter N as: the parameter can determine the phase shift value to be 2 pi/N, and simultaneously, the parameter is related to the highest power of independent variables of a synchronous detection polynomial in a phase shift algorithm to determine the acquisition frame number of the interferogram.

In order to achieve the purpose, the invention adopts the following inventive concept:

a multi-surface measurement method based on phase shift characteristic polynomial high-precision fitting comprises the following main steps:

1. solving the direct independent variable coefficient expression of the characteristic polynomial based on the characteristic polynomial, wherein the direct independent variable coefficient expression is an ideal window function numerical value; fitting coefficient expansion of an independent variable through Gaussian fitting, polynomial fitting and Fourier fitting methods with superior performance to obtain a fitting form of a window function of an optimal multistep weighted phase shift algorithm under different algorithm parameters; meanwhile, the result is calculated, and the selected fitting method has higher fitting precision; and simultaneously, giving a general formula by using the three fitting forms so as to select the three fitting forms according to different measurement conditions and error tolerance requirements during measurement.

2. The method comprises the steps of designing key parameters of an algorithm by a numerical analysis method, obtaining the range of the position of an available measured piece corresponding to different algorithm parameters, digitizing the available range based on the algorithm characteristics and the numerical analysis result, and calculating the general expression of the corresponding available cavity length coefficient under different phase shift parameters, so as to be directly applied during measurement and realize the unification of an experimental scheme and the algorithm design.

Therefore, under the condition of knowing the thickness and the refractive index of the measured piece, according to the current cavity length value, a proper phase shift parameter can be selected, and the minimum acquisition frame number and the single-step phase shift value necessary for the interferogram can be determined through the coefficient.

According to the inventive concept, the technical scheme adopted by the invention is as follows:

a multi-surface measurement method based on phase shift characteristic polynomial high-precision fitting comprises the following steps:

(1) calculating through a characteristic polynomial in a multi-weighting multistep phase shift algorithm to obtain coefficient distribution of independent variables of the characteristic polynomial, namely ideal window function distribution, and expanding coefficients from the polynomial, so that the coefficient distribution cannot be directly utilized in engineering;

(2) for convenience of application in measurement, fitting coefficient distribution of independent variables of the characteristic polynomial under different parameters, and analyzing distribution of residual errors solved by an algorithm of available phase shift parameters N under different cavity length coefficients M;

(3) setting algorithm parameters by a numerical analysis method, obtaining the range of the position of the detected piece corresponding to different algorithm parameters, digitizing the range of the position of the detected piece based on the algorithm characteristics and the numerical analysis result, and calculating the general formula of the cavity length coefficient corresponding to different phase shift parameters, so as to be directly applied during measurement and realize the unification of the measurement scheme and the algorithm design;

(4) according to the actual requirements and subjective wishes of measuring personnel, the window function in the measuring scheme is selected, more feasibility for making the measuring scheme is provided, and the measuring cost is saved, so that the one-time multi-surface non-contact measurement is customized and completed.

Preferably, in the step (2), gaussian fitting, fourier fitting and polynomial fitting are performed on the characteristic polynomial through the distribution of the independent variable coefficients of the polynomial, and numerical results of three fitting types are provided for selection.

Preferably, in the step (2), the three fitted window functions are calculated by fitting them, and the general formula is expressed as:

where k is the phase shift ordinal number, i.e. the frame order number of the interferogram,

for sampling window functions using a polynomial fitting method,for the sampling window function using the fourier fitting method,is a sampling window function using a gaussian fitting method.

Preferably, in the step (2), when performing polynomial fitting, a phase shift ordinal number k, that is, an independent variable coefficient distribution ordinal number, is used as a variable of a fitting value, and a corresponding current polynomial distribution coefficient is used as a fitting target value; sampling window function based on polynomial fittingIs expressed as:

wherein a is₀～a₅Is a parameter to be solved; the minimum Q of the squared error is found and expressed as:

wherein Z is the total acquisition frame number of the interferogram; by using the principle of least square, solving the square value of the above formula and making the square value tend to 0, thus obtaining the distribution of each coefficient; specifically, the partial derivative of each coefficient value is calculated to be 0, so as to obtain the optimal solution of the extreme value:

by solving the above formula, the equation expression of least square solution of each coefficient can be obtained respectively, and the least square solution can be obtained by solving the equation to complete parameter fitting.

Preferably, in the step (2), when fourier fitting is performed, the distribution of the independent variable coefficients of the characteristic polynomial is also used as a fitting ideal value, and the sampling window function of the fourier fitting method is based onCan be expressed as:

wherein b is₀,b₁,b₂,Omega is a parameter to be solved; the solving of the coefficient can be completed by a least square method, but in order to ensure the solving precision, the Newton iteration method is adopted for solving; constructing an iterative function

When the point of the next iteration is x_k(is a vector including the above-mentioned parameters to be solved), and the next iteration point is x_k+1Wherein x is_k+1In relation to the derivative of the iterative function, it is written in detail as:

whereinFor the current iteration point; therefore, the direction of the next iteration point is related to the function value and the derivative value of the current iteration point; selecting between initial iteration points 0-1And selecting the iteration error of a vector consisting of random numbers to be 0.0001, selecting the maximum iteration number to be 1000 steps, stopping iteration when the iteration error or the maximum iteration number is met, and outputting a current value as a coefficient of solving the fitting.

Preferably, in the step (2), the gaussian fitting is performed by analyzing the distribution of the target window function to construct a gaussian fitting form thereof

The expression is as follows:

wherein c is_m,(m ═ 1,2) is the fitting parameter to be solved; wherein each parameter in the gaussian fitting needs to be solved accurately to ensure the fitting accuracy. Therefore, the invention constructs an error function by utilizing a successive iteration algorithm of a nonlinear least square method, takes the minimum error function as an iteration target, adopts a confidence domain algorithm, determines the radius of the region by the descent rate of a quadratic approximation function, and then obtains the second-order approximate minimum value of the target function in the radius range. If the minimum value is sufficient to drop the objective function, then the next iteration is entered and the radius of the confidence domain is expanded. If the minimum value does not allow the objective function to drop sufficiently, it indicates that the second order approximation of the current confidence domain is unreliable, and therefore the radius of the confidence domain needs to be reduced and the minimum value recalculated. This iteration will continue until the conditions required for convergence are met. First, the squared error is defined as the objective function f (k) (in the form of a complete representation function whose argument k is shown in parentheses in the body of the function, whose meaning is unchanged, e.g., w)_k＝w_k(k))：

On the circle around the point of iteration of each step,second order approximation function for approximating objective function by current point

WhereinAndare the first and second order gradients of the objective function. H () represents the Hessian matrix of the acquisition object.

J^T(k) Is the transposition of the jacobian matrix of the objective function. First, consider that

For the minimum solution, p_kIs the gradient descent direction along which the minimum in the trust domain is found. From the initial point, the iteration vector can be represented as:

wherein Δ_kThe confidence domain radius is adjusted by the ratio of the fall of the quadratic approximation function to the fall of the objective function. If the ratio is close to 1, the radius is increased; if the ratio is equal to 1, the radius remains unchanged; if the ratio is less than 1, the radius is decreased. Tau is_kIs determined by:

the initial point may be based on experience or set to a random number between 0-1. In order to ensure the iteration precision and the result reliability under the maximum iteration step length, the iteration parameters are set as follows: the maximum and minimum changes of the finite difference gradient coefficients are set to 0.1 and 10, respectively^-8(ii) a The maximum evaluation times and the allowed iteration times of the iteration model are 600 and 400 respectively; minimum iteration error of 10^-6. And if the maximum iteration number or the minimum iteration error is met, stopping iteration and outputting a current value.

Thus, the three fitted window functions can be expressed as:

preferably, in the step (3), a quantitative analysis result under the condition of available cavity length coefficients under different phase shift parameters N is summarized according to a comprehensive analysis of error distribution results by using a numerical analysis method, and is represented by using a general formula:

comprehensively considering the cavity length position and the wavelength minimum tuning range of the measured piece, taking N as an integer of 8-14, and taking r as a non-negative integer of 0,1,2 and 3 … …; when the available N values are addressed at the cavity length factor M, the r values increase from 0 and N increases from 8, addressing is done stepwise.

Preferably, when different N overlaps with a single cavity length coefficient M, selecting smaller value N, and performing optimized selection of the value N. When the method is used for addressing, in order to reduce the cost of algorithm and measurement implementation to the maximum extent, when different N overlaps with a single cavity length coefficient M, N with a smaller value should be selected, and the processing can reduce the number of acquisition frames necessary for an interferogram, reduce the possibility of error inclusion and realize the optimal selection of the N value.

Preferably, the acquisition of the interferogram is performed by using an INF600-LP interferometer.

Preferably, during the evaluation of the expression of the coefficients of the characteristic polynomial, it detects synchronously the highest power n of the polynomial_maxCan inhibit the maximum n_max-phase shift error of order 2, this parameter being customizable.

Preferably, the clamp of the tested piece is an adjusting frame clamp of a large self-centering adjusting frame FPSTA-4SCML-12 model.

Preferably, after solving for a specific value of the polynomial coefficient, fitting the value to facilitate use in measurement, the fitting method selects:

a. fitting an integral polynomial;

b. fourier fitting;

c. performing Gaussian fitting;

and the calculated amount and the fitting precision are considered, and the optimal fitting expression is solved through a corresponding solving method respectively for selection.

Preferably, after solving each fitting expression, expressing the fitting expressions through digitization based on error distribution to give values of available measurement positions corresponding to each algorithm parameter; it is particularly noted that when different algorithm parameters correspond to overlapping available measurement positions, the minimum algorithm parameter corresponding to the current measurement position should be selected, that is, the algorithm parameter corresponding to the minimum necessary interferogram acquisition frame number is selected as the optimal parameter.

Compared with the prior art, the invention has the following obvious and prominent substantive characteristics and remarkable technical progress:

1. according to the invention, three fitted sampling window functions are provided for selection through fitting of the independent variable coefficient distribution of the characteristic polynomial, and higher fitting precision is ensured;

2. the invention carries out visual error analysis on the algorithm key parameter-phase shift parameter under different cavity length coefficients, and obtains the general expression of the numerical distribution range of the available cavity length coefficients under different phase shift parameters, thereby reducing the application and calculation cost;

3. according to the method for judging the optimal phase shift parameter, the selection result of the optimal phase shift parameter under the current cavity length coefficient can be obtained, the number of frames to be acquired is reduced to the greatest extent, the measurement process is simplified, and errors are avoided as much as possible;

4. the invention can select the window function in the measuring scheme according to the actual demand and subjective will of the measuring personnel, provides more feasibility for making the measuring scheme, saves the measuring cost and can realize one-time multi-surface non-contact measurement.

Drawings

FIG. 1 is a schematic diagram of RMSE for three fitting approaches.

FIG. 2 is a schematic representation of SSE for three fitting approaches.

FIG. 3 is a graph of front surface error distribution for different cavity length coefficients and phase shift parameters.

FIG. 4 is a plot of the back surface error distribution for different cavity length coefficients and phase shift parameters.

FIG. 5 is a graph of thickness variation error distribution for different cavity length coefficients and phase shift parameters.

Fig. 6 is a diagram of the results of the dephasing.

Detailed Description

The invention will be further described with reference to the drawings and preferred embodiments. The following examples are only for illustrating the technical solutions of the present invention more clearly, and the protection scope of the present invention is not limited thereby.