A method for measuring the three-dimensional topography of an object surface based on optical interferometry

Abstract

This invention relates to the field of optical detection technology and discloses a method for measuring the three-dimensional topography of an object surface based on optical interferometry. The method includes: acquiring a mixed temporal-domain interferometric signal sequence during axial scanning; extracting seed pixels and constructing a spatial topological connectivity feature matrix; introducing a fractional Fourier transform to construct a joint decoupling objective function containing a frequency-domain nonlinear phase distortion evaluation index and a temporal-domain physical invariant evaluation index; performing multi-level optimization for the seed pixels to lock the globally optimal rotation angle; performing array adaptive decoupling in conjunction with the spatial topological connectivity feature matrix to generate the optimal rotation angle matrix; and performing signal decoupling and reconstruction on the mixed signal to extract the absolute surface height of each pixel to generate three-dimensional topography data. This invention solves the phase decoupling problem caused by thin-film interferometric aliasing through multi-dimensional joint decoupling and spatial prior constraints, improving measurement accuracy while reducing the computational overhead of large pixel arrays.

Description

Technical Field

[0001] This invention relates to the field of optical detection technology, specifically to a method for measuring the three-dimensional morphology of an object surface based on optical interference. Background Technology

[0002] In precision manufacturing fields such as semiconductor manufacturing, microelectromechanical systems (MEMS), and optical processing, three-dimensional topography measurement technology based on the principles of white light interferometry or phase-shifting interferometry is a key means to achieve micro- and nano-level precision detection. However, in actual detection scenarios, the surface of the object being measured is often covered with single or multiple transparent thin films, such as silicon oxide passivation layers or photoresist coatings. When the physical thickness of the thin film is close to or less than the coherence length of the light source, the optical signals generated by reflections from the upper and lower surfaces of the thin film will undergo severe coherent overlap in the time domain, forming mutually overlapping interference envelopes.

[0003] Due to the presence of this physical aliasing, the frequency domain phase of the interference signal exhibits a nonlinear distribution, making it impossible for traditional frequency domain analysis methods to accurately extract the phase slope reflecting the surface height through linear fitting, thus leading to significant deviations in measurement results. To eliminate the phase distortion introduced by the thin film, some existing technologies attempt to introduce time-frequency analysis techniques such as fractional Fourier transform to parametrically decouple the aliased signal. However, such methods have obvious limitations in application. Because large-scale iterative optimization is required for parameters such as the rotation angle of the interference signal, when the pixel scale of the area array detector reaches millions, pixel-by-pixel global blind search will cause huge computational overhead, resulting in the processing time of a single frame of data far exceeding the allowable range of industrial online monitoring.

[0004] Existing decoupling schemes typically treat each pixel in the array as an independent computational unit, failing to fully utilize the continuous characteristics of the object's surface morphology in physical space. This isolated computational mode is prone to getting stuck in local extrema due to fluctuations in the signal-to-noise ratio of interference signals when faced with complex conditions such as surface steps, scratches, or abrupt changes in film thickness. This leads to high jumps between sampling points and data distortion, making it difficult to simultaneously achieve both the accuracy of morphology reconstruction and the processing efficiency of large-scale arrays. Summary of the Invention

[0005] To address the shortcomings of existing technologies, this invention provides a method for measuring the three-dimensional topography of object surfaces based on optical interference, which solves the problems of limited single-surface phase decoupling accuracy caused by thin-film interference signal aliasing and low computational efficiency when performing parameter optimization for large-scale pixel arrays.

[0006] To achieve the above objectives, the present invention provides the following technical solution: a method for measuring the three-dimensional topography of an object surface based on optical interferometry, comprising:

[0007] S10, acquire the mixed temporal interference signal sequence covering the coordinates of each pixel in the two-dimensional array during the axial scanning process of the object under test;

[0008] S20, extract seed pixels and construct a spatial topological connectivity feature matrix that characterizes the spatial physical continuity of the surface of the object under test;

[0009] S30 introduces fractional Fourier transform to construct a joint decoupling objective function that includes frequency domain nonlinear phase distortion evaluation index and time domain physical invariant evaluation index;

[0010] S40, for the seed pixel, perform multi-level optimization with the minimization of the joint decoupling objective function as the convergence condition, and lock the global optimal rotation angle of the seed pixel;

[0011] S50, starting with the seed pixel, perform array adaptive decoupling in combination with the spatial topological connectivity feature matrix to generate the optimal rotation angle matrix corresponding to each pixel;

[0012] S60, based on the optimal rotation angle matrix, perform signal decoupling and reconstruction on the mixed temporal interference signal sequence of each pixel, and extract the absolute surface height of each pixel to generate the three-dimensional shape data of the object under test.

[0013] In one embodiment, the process of acquiring the hybrid time-domain interference signal sequence is as follows: the actuator is controlled to perform a physical scan with a preset step length along the axis, the detector is triggered to collect spatial interference patterns at each scanning step node, and the interference patterns at each scanning position are stacked according to the axial coordinates. The light intensity values ​​of the same pixel at different axial positions are extracted to form a one-dimensional data vector.

[0014] In one implementation, the process of constructing the spatial topological connectivity feature matrix includes: calculating the contrast of the interference fringes of the entire array of pixels, selecting the pixel with the largest contrast as the seed pixel; performing a fast Fourier transform on the mixed time-domain interference signal sequence to extract the initial frequency domain phase, and calculating the spatial phase gradient between adjacent pixels to obtain a first-order spatial phase gradient modulus array characterizing the degree of drastic change in the morphology of the measured surface.

[0015] In one embodiment, the extraction process of the frequency domain nonlinear phase distortion evaluation index includes: projecting the mixed time-domain interference signal sequence onto the fractional Fourier domain corresponding to a specific rotation angle, restoring it to the conventional wavenumber domain after energy filtering to obtain the frequency domain phase sequence, calculating the second derivative of the phase using a differential filter to obtain the phase curvature, and calculating the weighted phase curvature variance with power spectral density as the weight, which serves as an index for quantifying the degree of nonlinear distortion introduced by thin-film crosstalk.

[0016] In one implementation, the extraction process of the evaluation index of time-domain physical invariants includes: performing an inverse transform on the filtered signal in the fractional Fourier domain to restore it to the reconstructed time-domain signal; extracting the peak coordinates of the interference envelope and the coordinates of the carrier zero-phase point of the reconstructed time-domain signal respectively; and calculating the axial physical displacement deviation between the two as an index for quantifying the residual asymmetric distortion in the time domain.

[0017] In one implementation, the process of constructing the joint decoupling objective function includes: assigning normalized weight coefficients to the weighted phase curvature variance and the axial physical displacement deviation respectively and summing them; the normalized weight coefficients are dynamically determined based on the reciprocal of the statistical mean of the two indicators in the parameter space.

[0018] In one implementation, the multi-level optimization includes: performing global coarse discrete step size optimization to locate coarse optimization candidate angles; introducing linear frequency modulated Z-transform in the neighborhood defined by the coarse optimization candidate angles; achieving local spectrum amplification by mapping the complex plane calculation trajectory to a narrow band spiral in the target frequency band; and performing high-resolution refinement optimization to lock the globally optimal rotation angle.

[0019] In one implementation, the process of performing array adaptive decoupling includes: planning a layer-by-layer traversal path that radiates outwards; mapping the element values ​​in the spatial topological connectivity feature matrix to the dynamic search radius of the current pixel to be calculated using a nonlinear exponential saturation growth model; and calling a local refinement optimization algorithm to determine the optimal rotation angle of the current pixel to be calculated within the boundary defined by the optimal rotation angle of the reference pixel and the dynamic search radius.

[0020] In one implementation, the reconstruction process includes: performing a fractional Fourier transform on the signal using the optimal rotation angle corresponding to each pixel; using a narrow-band filter window with a cosine roll-off coefficient in the fractional domain to suppress crosstalk components; and then performing an inverse fractional Fourier transform with opposite rotation angles to recover the single-surface pure time-domain interference signal.

[0021] In one implementation, the process of extracting the absolute surface height of each pixel includes: extracting a continuous linear phase sequence with respect to the wavenumber from the pure time-domain interference signal, calculating the first-order partial derivative of the continuous linear phase sequence with respect to the wavenumber using the least squares method to obtain the phase-wavenumber slope, and calculating the axial physical height of each pixel relative to the zero optical path difference position based on the phase-wavenumber slope.

[0022] This invention provides a method for measuring the three-dimensional topography of an object surface based on optical interferometry. It has the following advantages:

[0023] 1. This invention solves the envelope separation problem caused by the high overlap of thin-film interference signals in the time domain by constructing a joint decoupled objective function that includes a frequency-domain nonlinear phase distortion evaluation index and a time-domain physical invariant evaluation index. The frequency-domain phase curvature variance characterizes the nonlinear distortion caused by dispersion, while the time-domain axial physical displacement deviation provides positional constraints. The synergistic mechanism of the two can lock the optimal rotation angle under a highly degenerate group delay, achieving accurate reconstruction of the pure signal from a single surface and improving the absolute accuracy of three-dimensional topography measurement.

[0024] 2. This invention introduces a linear frequency-modulated Z-transform within the candidate neighborhood determined by coarse optimization, mapping the computational trajectory of the complex plane to a narrow-band spiral within the target frequency band, thus achieving local physical amplification in the frequency domain. This mechanism avoids the large-scale zero-padding operations performed by traditional methods to improve resolution, improving the accuracy and computational efficiency of parameter search without increasing the length of the original data.

[0025] 3. This invention establishes a linkage logic between two-dimensional spatial features and one-dimensional time-frequency signal optimization by introducing prior constraints on spatial topological connectivity, significantly improving the processing speed of large-scale pixel arrays. By utilizing a nonlinear exponential saturation growth model to dynamically map surface topography gradients to the search radius of the pixels to be calculated, adaptive adjustment of the optimization boundary is achieved. In regions with gentle topography, the search range is narrowed to reduce redundant computation, while in regions with abrupt topography changes, the boundary is dynamically widened to ensure that physical features are not lost, thus meeting the real-time requirements of online detection of complex surfaces. Attached Figure Description

[0026] Figure 1 This is a flowchart of the method of the present invention;

[0027] Figure 2 The graph shows a comparison of measurement error and calculation time for different film thicknesses according to the present invention, where (a) is a comparison of measurement error and (b) is a comparison of calculation time. Detailed Implementation

[0028] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0029] Please see the appendix Figure 1 This invention provides a method for measuring the three-dimensional topography of an object surface based on optical interferometry, comprising the following steps:

[0030] S10: Acquire a time-domain hybrid interference sequence array; control the piezoelectric ceramic actuator to perform axial scanning along the Z-axis with a preset step size, trigger the two-dimensional area array detector to acquire spatial interference patterns at each scanning step node, stack and combine the interference patterns of all scanning positions according to the Z-axis coordinates, and obtain a hybrid time-domain interference signal sequence covering the coordinates of each pixel in the two-dimensional area array.

[0031] S20: Extract seed pixels and initialize spatial topology; calculate the zero-order interference fringe contrast of all array pixels, select the pixel with the largest contrast as the starting point for optimization calculation and define it as seed pixel, perform fast Fourier transform on the array mixed signal to extract the initial frequency domain phase, calculate the first-order spatial phase gradient matrix between adjacent pixels, and construct the spatial topology connectivity feature matrix.

[0032] S30, construct a joint decoupling objective function; introduce a fractional Fourier transform under a specific rotation angle for the mixed signal, use a power spectrum weighted Savitzky-Gore differential filter in the frequency domain branch to process the separated signal to extract the noise-resistant phase curvature variance, and extract the deviation between the peak coordinates of the single-surface interference envelope and the zero-phase point coordinates through inverse transform in the time domain branch. The frequency domain parameters and time domain parameters are weighted and fused to construct a joint objective function.

[0033] S40, perform multi-level optimization of seed pixels; solve for seed pixels with the joint objective function minimization as the convergence condition, perform global coarse resolution optimization to obtain coarse optimization candidate angles, and introduce a local refinement optimization algorithm based on linear frequency modulated Z-transform based on the coarse optimization candidate angles to perform high-resolution refinement optimization, and lock the global optimal rotation angle of the seed pixels.

[0034] S50, perform array adaptive decoupling; traverse and calculate from the seed pixel to the surrounding adjacent pixels to be calculated, map the spatial phase gradient matrix to the optimization dynamic search radius of the pixels to be calculated, and call the local refinement optimization algorithm to lock the optimal rotation angle of each pixel in the full field of view array within the adaptive boundary defined by the optimal rotation angle of the adjacent known reference pixels and the dynamic search radius, and generate the optimal rotation angle matrix.

[0035] S60 reconstructs the pure interference signal and extracts the three-dimensional morphology; substitutes the optimal rotation angle matrix into the array operation channel, performs a refined fractional Fourier transform, narrowband filtering and inverse reconstruction at the corresponding optimal angle for each pixel, removes the crosstalk signal introduced by thin film interference, extracts the first-order partial derivative of the pure single-surface frequency domain phase with respect to the wavenumber after reconstruction to calculate the absolute surface height of each pixel, and generates the three-dimensional morphology data of the measured object by combining the absolute surface heights of each pixel.

[0036] The specific implementation principles of each step of the present invention will be explained in detail below with reference to the accompanying drawings.

[0037] In this embodiment, the specific execution process of step S10 can be implemented through the following sub-steps to complete the synchronous control of the underlying hardware of the white light interferometry system and the structured generation of basic raw data. Step S10 specifically includes the following sub-steps:

[0038] S101 controls the piezoelectric ceramic actuator to perform axial scanning along the Z-axis with preset step sizes. Before performing the physical scan, the single-step displacement amount satisfying the Nyquist sampling theorem must be determined based on the center wavelength of the broadband light source to prevent spatial frequency aliasing of subsequent interference signals. The white light interferometry system uses a low-coherence broadband light source for illumination. After passing through the interferometer objective, the beam is split into a measurement beam and a reference beam. The piezoelectric ceramic actuator is fixed to the external mounting base or the stage below the interferometer objective. The system control unit sends a stepped voltage drive signal to the piezoelectric ceramic actuator. Based on this drive signal, the piezoelectric ceramic actuator drives the interferometer objective to move in equal step sizes along the Z-axis direction perpendicular to the surface of the object being measured. The starting position of the Z-axis scan is set as... The end position of the scan is The physical displacement of a single step is .

[0039] The specific values ​​of the start and end positions are determined by the estimated thickness and surface morphology range of the object being measured, so as to ensure that the scanning path can completely cover the interference envelope of all physical interfaces to be measured.

[0040] For the closed-loop feedback control of piezoelectric ceramic actuators and the linear calibration of driving voltage and physical displacement, those skilled in the art can use existing capacitive displacement sensors or strain gauges for closed-loop calibration. The displacement closed-loop control is a well-known technology in the field and will not be elaborated here.

[0041] S102, at each scanning step node, triggers the two-dimensional area array detector to acquire spatial interferometry patterns. To ensure strict alignment of spatial coordinates and time series during multi-source data acquisition, a single pass is completed by the piezoelectric ceramic actuator. After the stepping motion reaches a mechanically stable state, the system control unit outputs a hard trigger level signal to the two-dimensional area array detector to avoid phase measurement errors caused by motion blur due to mechanical vibrations during the scanning process. The underlying hardware implementation of the two-dimensional area array detector is a charge-coupled device or a complementary metal-oxide-semiconductor image sensor. Within a preset exposure time after receiving the trigger level signal, the two-dimensional area array detector records the spatial interference pattern formed by the superposition of the measurement beam and the reference beam at the current Z-axis position. This spatial interference pattern is a two-dimensional grayscale image containing phase and intensity information across the entire field of view. (Corresponding to the...) Optical path difference of each step node The expression is as follows:

[0042] ;

[0043] In the formula, It is a non-negative integer representing the index number of the scan step. Its maximum value is determined by rounding down the ratio of the total scan travel to the single step displacement, and is used to define the sequence length of the discrete signal.

[0044] S103, the interferograms of all scanned positions are stacked and combined according to the Z-axis coordinate to obtain a mixed temporal interferometric signal sequence covering the coordinates of each pixel in the two-dimensional array. After obtaining the single-frame image sequence, a three-dimensional data space suitable for the underlying algorithm to extract one-dimensional vectors needs to be further constructed. The system control unit will move the two-dimensional array detector to each The 2D grayscale image output by the node is transmitted to the computer memory via the data bus. The computer orthogonally stacks the 2D images in 3D space according to the time sequence of the Z-axis scan, forming a 3D data cube.

[0045] In this three-dimensional data cube, the coordinates of any physical pixel in the spatial plane are extracted along the Z-axis. From the one-dimensional data vector, the temporal interference signal sequence corresponding to that pixel can be obtained. Since the surface of the object being measured has a transparent or semi-transparent film, the measurement beam will be reflected from both the upper and lower surfaces of the film. The temporal interference signal sequence appears on the time axis as a linear superposition of the interference envelopes of the upper and lower surfaces; this superimposed signal is defined as a hybrid temporal interference signal sequence. This sequence serves as the core data carrier characterizing the optical properties of the object under test, and it fully maps the physical quantities of depth and reflectivity of each reflective interface. The physical model expression of the hybrid time-domain interferometric signal sequence is as follows:

[0046] ;

[0047] In the formula, Representing coordinates Optical path difference is Interference light intensity at time; This indicates the DC background light intensity of that pixel; These refer to the upper and lower surfaces of the film on the object being tested, respectively. This represents the contrast coefficient of the interference fringes on the corresponding surface. Its value ranges from 0 to 1, and its specific value is positively correlated with the reflectivity of the corresponding physical interface. The low-coherence envelope function representing a broadband light source is often approximated as a Gaussian distribution function in practical engineering models; This indicates the physical position of the Z-axis when the corresponding surface produces zero optical path difference; The center wave number of a broadband light source satisfies ,in The center wavelength of a known broadband light source; This represents the initial phase introduced by the corresponding surface reflection. This hybrid time-domain interferometric signal sequence serves as the basic raw data, which is then input into the subsequent algorithm processing channel to perform noise reduction and signal decoupling.

[0048] In this embodiment, step S20 specifically includes the following sub-steps:

[0049] S201, calculate the zero-order interference fringe contrast of the entire array of pixels, and extract the pixel with the highest contrast as the seed pixel. This is applied to the hybrid temporal interference signal sequence obtained in the previous step. The system iterates through the coordinates of each pixel in the two-dimensional array. The modulation index of the corresponding interference signal is calculated. In this data processing stage, to eliminate invalid data interference caused by strong reflections, pixels that cause photoelectric sensor saturation need to be pre-removed. Specifically, the saturation light intensity threshold of the detector is set to... This threshold is determined based on the camera's analog-to-digital conversion bit depth; for example, it is 255 for an 8-bit detector. When the maximum light intensity value in the pixel signal reaches or exceeds... When the contrast is forced to zero, for each unsaturated effective pixel, the maximum and minimum light intensities in the interference signal envelope are extracted, and the zero-order interference fringe contrast is calculated. To prevent detector defects or extremely dark areas from causing the denominator to approach zero and triggering a division-by-zero anomaly, a very small positive bias constant is introduced into the denominator term. The formula for calculating contrast is as follows:

[0050] ;

[0051] In the formula, Representing coordinates Contrast of interference fringes at the location; This represents the peak light intensity in the temporal interference sequence of that pixel; This indicates the corresponding nearest valley light intensity or the global minimum light intensity; The bias constant, used to prevent calculation overflow, is typically set to a value range of 10. -6 Up to 10 -4 The decimal places are determined by the level of background noise in the system.

[0052] Global extremum optimization is performed across the entire effective pixel range to find the contrast ratio. To obtain the physical location of the global maximum contrast, and to avoid misjudgment of a single extreme value caused by isolated noise points, as a preferred method, after determining the candidate pixel for the global maximum contrast, it is necessary to further calculate the contrast variance within a preset spatial neighborhood (e.g., a 3×3 pixel window) centered on that pixel. Only when this local variance is less than a set threshold can it be confirmed as a valid seed pixel. If the variance exceeds the limit, the next candidate point will be evaluated in descending order of contrast value.

[0053] S202 performs a Fast Fourier Transform on the array's mixed signals to extract the initial frequency domain phase. Using a two-dimensional spatial grid as a reference, it performs a mixed time-domain interference signal sequence on each pixel of the entire array. A Discrete Fast Fourier Transform (DFT) is performed along the Z-axis optical path difference direction. This general transformation process converts the one-dimensional time-domain interference signal to the wavenumber domain, generating a frequency-domain complex signal containing multiple frequency components. To extract fundamental phase information that characterizes surface morphology, the system uses a broadband light source with a central wavenumber of... At this point, the real and imaginary parts of the complex signal are separated, and the initial frequency domain phase is then calculated. In specific programming operations, to avoid computational overflow when the real part approaches zero and to determine the accurate quadrant of the phase, the four-quadrant arctangent function is used for solving the problem. Its mathematical expression is:

[0054] ;

[0055] In the formula, This represents the operation of extracting the imaginary part of a complex signal in the frequency domain; This indicates the operation of extracting the real part; The arctangent function is defined in the fourth quadrant, and its output range is... ; The specific values ​​are obtained based on the calibration parameters of the system's light source. For discretization operations based on the Fast Fourier Transform, those skilled in the art can use existing radix-2 time-domain decimation algorithms, the underlying mathematical implementation of which is a well-known technology in the field and will not be elaborated here.

[0056] S203, calculate the first-order spatial phase gradient matrix between adjacent pixels to construct the spatial topological connectivity feature matrix. Based on the obtained initial frequency domain phase of the entire field of view. In a two-dimensional spatial domain, a difference operation is performed to quantize the spatial phase change rate between adjacent physical pixels. It is particularly important to note that the initial frequency domain phase obtained through direct solution... exist Phase truncation necessitates performing phase unwrapping along the two-dimensional spatial path before spatial differencing, or performing modulo operations after the differencing operation to eliminate spurious gradients caused by truncation jumps. After eliminating the truncation effect, the first-order spatial partial derivatives of the initial frequency domain phase along the horizontal and vertical directions of the two-dimensional image are calculated, and the square root of the sum of their squares is taken to obtain the magnitude of the first-order spatial phase gradient. Specifically, this partial derivative can be obtained through the central finite difference method between pixels or by using the Sobel operator and convolution with the image matrix. This magnitude array is defined as a feature matrix characterizing the spatial topological connectivity of the measured surface. The corresponding formulas for continuous domain theory are as follows:

[0057] ;

[0058] In the formula, and These represent the partial derivative operators along the X-axis and Y-axis of the image, respectively. Feature matrix This matrix is ​​used to numerically characterize the spatial continuity of the physical state of the surface of the measured object. When the value of an element at a certain coordinate is small, it indicates that the surface morphology of the corresponding physical region is gentle and the film thickness is uniform; when the value is large, it indicates that there is abrupt change in morphology or a jump in film thickness in that region. This feature matrix establishes a correlation between one-dimensional signal optimization and two-dimensional spatial physical structure, and will serve as a dynamic constraint input parameter for the adaptive optimization boundary in subsequent calculation steps.

[0059] In this embodiment, to achieve effective blind separation of aliased interference signals, a mathematical metric capable of accurately evaluating the purity of the separation must be established. Based on this, step S30 specifically includes the following sub-steps:

[0060] S301 projects the time-domain interference signal into the fractional Fourier domain. For the extracted mixed time-domain interference signal sequence, to separate the overlapping interference envelope caused by the thin film of the object under test, a fractional Fourier transform is introduced to rotate its time-frequency characteristics. At the underlying physical mechanism, the interference envelopes generated by the upper and lower surfaces of the transparent film exhibit slight differences in group delay rates due to material dispersion differences. In the conventional time or frequency domain, this difference is masked by aliasing energy; however, the fractional Fourier transform, by introducing a rotation angle parameter, is equivalent to rotating the time-frequency coordinate system, causing the signal component with a specific group delay rate to appear as a highly concentrated pulse in the corresponding fractional domain. Let the discrete time-domain interference signal of a single pixel be... Introducing a rotation angle of The fractional Fourier transform maps it to fractional domain coordinates. In the generated fractional spectrum The expression is as follows:

[0061] ;

[0062] In the formula, The physical coordinates of the axial optical path difference sequence of the input signal; These are the coordinates of the fractional domain after projection; Let be the rotation angle, and its value range is within . Takes values ​​consecutively between; Let be the integral kernel function of the fractional Fourier transform. This represents an infinite integral operation along the entire theoretical optical path difference domain; Let be an integral infinitesimal element along the axial optical path difference direction. To cover all projection transformation cases and avoid singularities, the integral kernel function... Based on rotation angle Segmentation definition:

[0063] ;

[0064] In the formula, The imaginary unit (satisfying) ); Pi is a constant. The cotangent function of the rotation angle; The cosecant function of the rotation angle; For the natural constant Exponential functions with base 0; It is an integer; is the Dirac impulse function, used to characterize the ideal focusing state of a signal at a specific angle.

[0065] Specifically, system computing Obtain the energy spectrum, pinpoint the position of the highest global peak as the center frequency, and decay the energy to the value of that peak. The horizontal coordinate range at a certain point is used as the dynamic filtering bandwidth to retain the principal component range with the highest energy and suppress crosstalk components at other thin film interfaces.

[0066] S302 employs a power spectral weighted Savitzky-Gore differential filter to extract the noise-resistant phase curvature variance. In the frequency domain branch, to calculate the physical parameters required for evaluation, the system first performs an inverse fractional Fourier transform (i.e., applies) to the filtered single-peak signal in the fractional domain. The opposite rotation angle is then restored to the normal time domain, and a fast Fourier transform is performed to map it to the normal wavenumber domain, thus obtaining the corresponding rotation angle. frequency domain phase sequence With power spectral density sequence ,in This represents the discrete wavenumber index within the effective light source frequency band. To evaluate the nonlinear phase distortion caused by dispersion and residual crosstalk, the second derivative of the phase, i.e., the phase curvature, needs to be calculated. Directly performing a quadratic center difference operation on the discrete phase will drastically amplify high-frequency shot noise, masking the true physical morphology. Therefore, a Savitzky-Gore differential filter is introduced to address this issue. The local window data is fitted using polynomial least squares, and the smoothed, noise-resistant phase curvature is directly output. Set the half-width of the sliding window to [value]. The order of the polynomial is To preserve low-frequency phase topography information and filter out high-frequency shot noise, the parameters are typically set to... and .

[0067] Furthermore, to quantify the degree of global frequency domain nonlinear distortion, power spectral density is used. As point-by-point confidence weights, the weighted phase curvature variance is calculated. The core physical purpose of this weighting mechanism is to reduce the interference of low signal-to-noise ratio regions at the edges of the light source spectrum on the overall evaluation. Its calculation model is as follows:

[0068] ;

[0069] In the formula, It is the weighted average of the power spectrum of the anti-noise phase curvature corresponding to all wavenumber points at the current rotation angle; This means iterating through and summing all qualified wavenumber bands within the set effective bandwidth threshold for the current pixel; To prevent extremely small normal values ​​from causing division by zero anomalies when the denominator approaches zero in the background noise range, their value is fixed at 10. -8 Weighted phase curvature variance The smaller the value, the higher the frequency domain phase linearity separated at that rotation angle, and the lower the residual thin film crosstalk.

[0070] S303 extracts the deviation between the peak coordinates and the zero-phase point coordinates of the single-surface interference envelope in the time domain. In the time domain branch, the system simultaneously performs an inverse fractional Fourier transform on the filtered signal in the fractional domain to restore the reconstructed time-domain signal. For a pure, ideal single-surface white light interference signal, the energy peak position of its low-coherence envelope must strictly coincide with the zero-phase point position of the interference carrier on the Z-axis physical coordinate system. When there is crosstalk from thin film layers that has not been completely separated in the signal, group delay dispersion will cause the envelope peak to shift, thus violating this physical invariant.

[0071] The system Perform a Hilbert transform to construct an analytic signal, extract the maximum index of the absolute value sequence of the analytic signal's amplitude, and map it to the axial absolute coordinate of the single-surface envelope energy peak. Simultaneously, within a threshold neighborhood set by the envelope peak value, a zero-crossing point is searched where the reconstructed interferometric carrier amplitude is zero and the slope is positive. To avoid erroneous locking to adjacent carrier periods, this search neighborhood is strictly limited to... Within the closed interval, where To determine the center wavelength of the system's broadband light source, the system calculates the axial physical distances between all positive slope zero-crossing points and the envelope peak coordinates within this interval, uniquely identifying the nearest zero-crossing point, and then precisely locating the axial absolute coordinates of the zero-phase point through local linear interpolation. Calculate the absolute deviation between the two. , used to quantify temporal residual asymmetric distortion.

[0072] S304. Construct a joint objective function based on time-frequency dual-domain invariants. This involves system-level fusion of frequency-domain distortion evaluation metrics and time-domain physical invariant evaluation metrics to construct a joint objective function that guides subsequent optimization. This function requires that the separation operation simultaneously satisfy the conditions for the smoothest phase curvature in the frequency domain and the highest overlap between the envelope and phase in the time domain. The mathematical closed-loop model expression for the joint objective function is:

[0073] ;

[0074] In the formula, and The corresponding normalized weight coefficients are determined by dynamically normalizing the variance scale of both indicators, since the magnitude of the phase curvature variance is usually much smaller than that of the axial coordinate deviation. This is to ensure the balance of the multi-dimensional evaluation logic along the optimization gradient direction and to avoid the algorithm falling into biased judgments due to the numerical influx of a single dimension. Specifically, the system pre-samples several points randomly in the parameter space and calculates the mean of each indicator. To ensure computational robustness, when the statistical mean is greater than a preset minimum value... (e.g. 10) -5 When ), assign the reciprocal of its mean to and If the mean is less than or equal to If so, then directly set the corresponding weight coefficient to 1, and the objective function... As the core convergence criterion, when When the global minimum value is obtained, the corresponding rotation angle is determined to be the optimal transformation angle that enables the current physical pixel to achieve decoupling and separation.

[0075] In this embodiment, considering the physical limit of group delay degeneracy in the extreme thin-film interference signal, performing a high-resolution blind search directly in the global domain would lead to an exponential increase in computational load. Based on this consideration, in order to obtain high-precision decoupling parameters with controllable computational overhead, step S40 specifically includes the following sub-steps:

[0076] S401 optimizes the coarse discrete step size of a single pixel in the global range by minimizing the joint objective function as the convergence condition. For the seed pixels extracted in the previous step, a rotation angle is set. The global search physical domain. To prevent the cotangent term in the fractional-order transformation kernel function... When the value approaches infinity at the boundary, causing a floating-point overflow, the system forcibly limits the search domain to [value missing]. Closed interval, where This is the minimum safety margin constant, typically taken as 10. -6To achieve an initial balance between computational efficiency and search accuracy, the global coarse discrete step size is configured as follows: The specific value of this coarse step length is determined based on the ratio of the coherence length of the broadband light source to the estimated thickness of the film being measured, and is usually set to... to Between radians, to ensure the physical step size does not exceed half the width of the main lobe of the signal. The system control unit drives the internal algorithm module to make the rotation angle between radians. The incremental iteration is performed within the search domain using discretization. For each discrete angle state, the previously constructed joint objective function model is invoked to calculate the corresponding evaluation value. After the traversal cycle ends, the discrete angles that enable the joint objective function to reach its global minimum are extracted. To prevent spurious extremum interference caused by sudden high-frequency noise in local areas, the system synchronously calculates the first-order difference of the function at adjacent step sizes before and after the minimum point. Only when the local region satisfies the physical condition of a stable valley bottom where the first-order difference changes sign and the second-order difference is positive can it be defined as a coarse-optimization candidate angle. At the physical level, this coarse optimization operation can quickly traverse the parameter space containing a large number of non-target aliasing signals and roughly locate the approximate time-frequency coordinate domain where the group delay dispersion introduced by the reflected light from the upper and lower surfaces of the thin film begins to decouple.

[0077] S402, candidate angles and candidate interval ranges determined based on the coarse optimization results. Stable coarse optimization candidate angles locked in the above stages. A local search domain for high-resolution analysis is constructed centered on this angle. Specifically, the candidate interval range for local refinement optimization is set as follows: Closed interval. As a preferred approach, considering the periodicity of the fractional Fourier transform in the angular domain, if the calculated interval boundary exceeds the aforementioned... The effective safety domain is defined by adding or subtracting parts that exceed the boundary. The period is folded and mapped back to the valid interval. This interval delineation logic is based on the local unimodal physics assumption, ensuring that the real physical extreme points are completely enveloped within the search range. This reduces the computational scale by hundreds of times while mathematically avoiding the omission of the global optimal solution.

[0078] S403 introduces a linear frequency modulated Z-transform within the candidate interval to adjust the sampling scaling factor for high-resolution local refinement optimization. Once within the local candidate interval, a very small refinement step size must be used. The process continues iterating through angles, and the resulting fractional-order spectra must possess extremely high wavenumber resolution to distinguish minute dispersion differences. Traditional Fast Fourier Transform (FFT) is constrained by uniform sampling and full-band computation, making it impossible to locally improve frequency domain resolution without increasing the length of the original time-domain sequence. Therefore, a linear frequency-modulated Z-transform is introduced into the underlying computational mechanism. This general signal processing technique maps the computational trajectory on the complex plane from a uniform unit circle to an arbitrarily specified narrowband spiral, allowing for the configuration of high-density independent sampling points within a defined target frequency band, thus achieving localized physical amplification in the frequency domain. Let the given discrete rotation angle within the local candidate interval be... The mathematical model for the linear frequency modulated Z-transform is as follows:

[0079] ;

[0080] In the formula, Indicates the rotation angle The high-resolution fractional-order spectral sequence after local magnification; Represents the original discrete-time domain interference sequence; The total number of data points in the original time-domain sequence is directly determined by the number of hardware scanning layers. The non-negative integer index of the discrete sampling point in the time domain; This is the optical path difference step; This is the discrete point index of the target high-resolution spectral domain, with a value range of [value range missing]. to ,in The preset number of local spectrum analysis points; It is the imaginary unit.

[0081] parameter and The core complex variable controlling the frequency sampling start point and sampling resolution is defined as follows:

[0082] ;

[0083] ;

[0084] In the formula, The starting physical coordinates of the principal energy of the interference spectrum of the thin film under test; This is the spatial sampling frequency of the original system scanning along the Z-axis, and its value is strictly equal to the physical step size of a single scan. The reciprocal of (i.e.) ); This is the sampling scaling factor, used to characterize the magnification of local resolution. To ensure the validity of matrix operations, the system enforces verification and constraints. And it must be a positive integer, while ensuring the sampling frequency. Sampling scaling factor The determination is based on the target observation bandwidth. and the number of analysis points The correlation calculation shows that it satisfies By dynamically adjusting this factor, it is equivalent to constructing a local high-resolution observation window in the frequency domain. Without increasing the amount of original data, it performs high-density analysis only on the weak features of generation delay separation within the candidate frequency band, significantly improving the sensitivity of the objective function in this interval.

[0085] S404 optimizes the output through local refinement and locks the globally optimal rotation angle of the seed pixel. Under the aforementioned finely constructed high-resolution data channel, the system uses a refinement step size... Gradually update the rotation angle within the candidate interval. The refinement step size is typically forced to be on the order of one-thousandth of the coarse step size, with its specific lower limit determined by the floating-point arithmetic precision of the microcontroller, thus penetrating the physical resolution limits of the instrument's hardware. At each refined angle node, the high-resolution spectrum is recalculated, and the joint objective function is extracted. The system extracts candidate refinement angles that minimize the joint objective function after traversing all subdivision nodes, and simultaneously calculates the average gradient of the objective function within a set window (e.g., containing 5 sampling points) on both sides of the extreme point. The validity of the extreme point is confirmed only when the objective function value corresponding to the extreme point is lower than a certain percentage (e.g., a set empirical threshold of 50%) of the global average value of the candidate interval, and the absolute value of the average gradient on both sides is greater than a preset smoothing threshold. The angle after the above multi-dimensional physical verification is then confirmed by the system as the optimal decoupling parameter representing the physical position corresponding to the seed pixel, and is locked as the globally optimal rotation angle. These optimal parameters not only perform noise reduction and desensitization themselves, but also serve as the core transmission benchmark for subsequent adaptive decoupling of the two-dimensional spatial array, and are written into the system's high-speed cache for adjacent pixels to read and access.

[0086] In this embodiment, performing a global blind search on each of the millions of pixels across the entire field of view would cause the system to encounter a computational bottleneck. To control the overall computational complexity, the system constructs a cooperative computing mechanism that introduces prior constraints of two-dimensional optical space into one-dimensional time-frequency signal optimization. Step S50 specifically includes the following sub-steps:

[0087] S501, a propagation path is planned layer by layer from the seed pixel as the origin to the adjacent pixels to be calculated on the periphery of the two-dimensional array. In optical interferometry, since the surface of the sample being measured usually has physically continuous growth properties, the film thickness and spatial morphology of adjacent pixels are strongly correlated in most areas. This constitutes the spatial prior basis for array decoupling. The seed pixel coordinates extracted in step S20 are used as the basis for this process. As the starting reference node, a region growing algorithm or a breadth-first search algorithm is used to plan the outward radiating computation path in a two-dimensional discrete pixel grid. As a preferred approach, the system employs the eight-neighbor connectivity criterion, extracting adjacent pixels to be computed layer by layer from the inside out. In actual processing, to avoid including invalid background areas in the calculation, which could cause singularity errors or unnecessary computational waste, the system pre-reads the zero-order interference fringe contrast matrix of all array pixels. Only when the contrast of adjacent candidate pixels is greater than a set effective signal threshold (e.g., an empirical value of 0.1 calibrated based on the system's background noise) can a pixel be pushed into the processing queue. To avoid the algorithm generating infinite loops in closed-loop paths or repeatedly calculating the same physical coordinates, the system synchronously constructs a binary marker matrix in memory, the same size as the image resolution. Whenever a pixel node completes decoupling calculation, its corresponding matrix element is marked as processed, thus guiding the optimization pointer to continuously advance along the uncalculated effective physical region until the effective data boundary of the entire field of view is covered.

[0088] S502 transforms the spatial phase gradient matrix into a dynamic search radius for the pixel to be calculated using a nonlinear exponential mapping model. This is done after determining the current pixel to be calculated. Subsequently, the system does not directly assign a fixed search range, but instead uses pre-extracted spatial topology priors for adaptive constraints. Specifically, it extracts the spatial topological connectivity feature matrix (i.e., the magnitude matrix of the first-order spatial phase gradient) constructed in step S20. The element value corresponds to the coordinates. This value intuitively represents the degree of drastic change in the physical topography of the current local region. The system introduces a nonlinear exponential saturation growth model to map this physical gradient as the dynamic search radius in the one-dimensional time-frequency optimization process. The calculation formula is as follows:

[0089] ;

[0090] In the formula, This represents the angle search half-width for the currently calculated pixel to be dynamically allocated; This is the preset maximum allowable search half-width increment, and its specific value is set based on the maximum possible thickness jump range of the film under test between adjacent pixels. It is typically set to a value of [value missing]. to radian; The minimum basic search half-width is used to compensate for the quantization error of the system hardware and avoid the collapse of the search interval when the local gradient approaches zero. Its value is usually set to 10 to 20 times the local refinement optimization step size. The curvature adjustment coefficient controls the sensitivity of the search radius to changes in spatial gradient. Its value is usually preset to the reciprocal of the average gradient of the surface of similar historical samples, and the range is usually between 0.1 and 1.0.

[0091] The above computational model establishes a non-proportional mapping between the gradient space and the angle optimization interval at the mathematical level. This ensures that an extremely narrow search window is used in areas with gentle topography to accelerate computation, while the search window is automatically widened in areas with steps or abrupt changes in topography to prevent the loss of physical features.

[0092] S503 invokes a local refinement optimization algorithm within the adaptive boundary defined by the optimal rotation angle of the reference pixel and the dynamic search radius. It extracts the optimal rotation angle locked by the reference pixel (i.e., the parent node of the current pixel) that is adjacent to the currently calculated pixel in the traversal propagation path and has already been solved. As a preferred approach, when multiple neighboring pixels that have already been solved exist around the pixel to be calculated, the system does not select randomly. Instead, it extracts the spatial phase gradient values ​​corresponding to these candidate reference pixels and selects the neighboring pixel with the smallest gradient (i.e., the smoothest physical topography and the strongest correlation) as the sole reference node. This approach, at the physical level, minimizes the propagation of solving errors across step edges. Combined with the dynamic search radius obtained from the above calculations... The system dynamically constructs a one-dimensional adaptive search boundary for the current pixel. In terms of boundary handling logic, if the calculated upper or lower limit exceeds the effective safety definition domain... If this happens, a hard truncation is forced to prevent overflow in subsequent cotangent function operations. Within this dynamically established local boundary, the system directly skips the time-consuming coarse discrete global traversal and directly calls the high-resolution local refinement optimization algorithm with linear frequency modulated Z-transform described in step S40 for calculation.

[0093] To ensure the robustness of the decoupling logic, the output of this local algorithm is monitored in multiple dimensions. If, within the currently constructed adaptive boundary, the joint decoupling objective function does not exhibit an effective minimum value that satisfies the conditions of first-order difference sign change and second-order difference positive equal stable valley bottom verification, it usually indicates that there is an unexpected topographic fault or foreign object interference at the current physical location, causing prior failure. In this case, the system triggers an anomaly self-recovery mechanism, automatically discarding the current local boundary and reverting to perform coarse-to-fine two-level optimization of the global range for the isolated pixel, thereby eliminating the risk of error accumulation and region misjudgment from an engineering mechanism perspective.

[0094] S504 traverses the entire field of view and generates the optimal rotation angle matrix for the entire array. As the aforementioned region growing and adaptive calculation logic spreads from the seed pixel to the edge of the 2D array, it calculates the optimal rotation angle parameters for each effective physical pixel within the entire field of view under the condition of achieving the extreme value of separation purity. The result data corresponding to these discrete coordinate points are then aligned row-by-row and assembled in memory. For invalid pixel holes that were not pushed into the calculation queue due to extremely low reflectivity, surface contamination, or exceeding the measurement range, numerical filling is performed using inverse distance weighted interpolation based on the optimal rotation angles of the surrounding effective pixels to ensure the continuity of the output physical field. When the marker matrix shows that all pixels determined to be effective have been traversed and the holes have been repaired, the system officially outputs the complete optimal rotation angle matrix for the entire array. At the overall data structure level, this matrix fully characterizes the optimal dispersion separation trajectory distribution of the multi-layer optical interface of the tested object in the fractional Fourier transform time-frequency space, providing high-resolution source data for subsequent reconstruction of the absolute three-dimensional morphology and thin film physical thickness of each separation interface.

[0095] In this embodiment, to achieve inverse physical metric recovery after signal dealiasing and ultimately complete the calculation of absolute topographic features, step S60 specifically includes the following sub-steps:

[0096] S601, substitute the full-array optimal rotation angle matrix into the data processing channel and execute the reconstruction instruction of thinning fractional Fourier transform pixel by pixel. For the full-array optimal rotation angle matrix generated in step S50... The system's main control module loads the data into the high-speed cache of the parallel data processing channel. Based on the discrete stepping scanning characteristics of the optical interferometer system in the depth direction, the acquired 3D interferometric data blocks and 2D spatial feature matrices belong to different data structure levels in memory. To ensure data alignment, the system uses the original spatial physical grid coordinates... Perform strict memory addressing mapping to extract the corresponding discrete-time domain hybrid interference signals one by one. Then, the optimal rotation angle parameters locked to the corresponding coordinates are called. The process initiates a high-resolution fractional Fourier transform instruction that incorporates linear frequency modulated Z-transform. This calculation projects the time-domain interferometric data corresponding to different pixels onto their respective fractional-order coordinate systems, each possessing optimal group delay decoupling characteristics, generating a specific fractional-order spectral sequence for each pixel. Through this pixel-by-pixel customized orthogonal transformation mapping, the multi-interface interference components that were originally mixed in the conventional time-frequency domain are transformed into discrete pulses with highly concentrated energy in the parameter space.

[0097] S602, a narrowband filter is applied in the fractional-order domain to remove the thin film interference crosstalk signal and perform inverse reconstruction. After the above projection operation is completed, due to the decoupling effect of dispersion at the optimal rotation angle, the interference energy peak generated by the target single surface and the crosstalk energy peak introduced by other interfaces of the thin film are aligned on the fractional-order coordinate axis. The above points appear as mutually separated local extrema. To effectively separate non-target crosstalk signals, the system constructs an adaptive narrowband filtering window within the current fractional-order domain. As a preferred method, system locking... The center coordinates of the main energy peak with the largest intermediate amplitude Furthermore, by combining the system hardware calibration constant with the reciprocal proportional relationship between the effective coherence length of the broadband light source, the dynamic filtering half-width is set. Based on this, to avoid the Gibbs effect and pseudo-envelope ringing distortion caused by hard truncation during time-domain reconstruction when directly using a rectangular window, the system generation center is located at... Half width is Furthermore, a graphical window with a cosine roll-off factor fixed between 0.1 and 0.25 is used. This filter window is multiplied by the original spectrum to suppress crosstalk components outside the main peak bandwidth.

[0098] Based on this, the filtered target single extreme value spectrum sequence Perform an inverse fractional Fourier transform. Its physical purpose is to remap the clean frequency domain features, filtered out of crosstalk noise, back to the normal coordinate domain. The rotation angle of the inverse transform is fixed at the inverse of the aforementioned optimal rotation angle. Its mathematical integral expression is as follows:

[0099] ;

[0100] In the formula, This is to obtain a pure time-domain interference signal after reconstruction, eliminating the thin-film coupling effect; The optical path difference axis is the physical coordinate. Fractional domain coordinates; This is the inverse transform kernel function carrying the opposite rotation angle. After this inverse reconstruction, the system recovers the target surface interference signal with high symmetry and a single envelope characteristic in the time domain channel.

[0101] S603, the first-order partial derivative of the reconstructed pure single-surface frequency-domain phase with respect to wavenumber is calculated to solve for the absolute surface height. Based on the physical principles of white-light interferometry, the interference phase generated by a single surface under broadband illumination exhibits a strictly linear mapping relationship with wavenumber, and the slope of this linear relationship is directly equivalent to the physical height of the target surface relative to the zero optical path difference position. The system uses the above-acquired pure time-domain interferometric signal... Perform a standard Fast Fourier Transform to extract the continuous phase within the effective light source frequency band, and then perform phase unwrapping operations to obtain the true physical phase sequence. ,in For the first A valid discrete wavenumber index, where the wavenumber is defined as the reciprocal of the wavelength.

[0102] To combat quantization noise introduced by discretization sampling and random phase disturbances caused by system mechanical vibration, the system abandons the method of directly dividing the single-point phase by the wavenumber. Instead, it calculates the altitude by taking the first-order partial derivative of the phase with respect to the wavenumber over a local continuous effective frequency band. The system uses the least squares method to linearly fit the phase-wavenumber curve to obtain the slope. Considering that sudden noise in actual operating conditions can easily lead to local phase unwrapping, the system takes into account this. Jump error, the system calculates and extracts the coefficient of determination while performing fitting ( (Test value). The system only adopts the fitted output when the coefficient of determination is greater than a preset linearity safety threshold (e.g., 0.95). The absolute physical height is constructed by combining a dual-path reflection model. The mathematical calculation model is as follows:

[0103] ;

[0104] In the formula, Represents pixel coordinates in two-dimensional space The absolute physical height of the target surface obtained from the solution; and These represent the horizontal and vertical discrete pixel index coordinates of the two-dimensional detector target surface, respectively. Pi is a constant. These are mathematical operators for partial derivatives; For continuous wavenumber variables, its physical meaning is defined as the reciprocal of the spatial wavelength; Wave number at corresponding pixel coordinates Continuously varying reconstructed pure frequency domain phase function; The positive integer index for summing the discrete sampling points participating in the linear fitting, with values ​​ranging from 1 to... ; This represents the total number of qualified wavenumber sampling points within the effective power spectral density threshold range of the light source; For the first Discrete wavenumber points; For all within the effective wavenumber interval The arithmetic mean of wavenumber points; For the corresponding spatial pixel coordinates in the th Pure unwrapped phase values ​​at each wavenumber point; Within the effective wavenumber interval, the above The arithmetic mean phase of each effective phase.

[0105] To ensure the validity and robustness of the calculation, the system performs a check on the denominator term before performing the division. Forced verification will be performed. If the effective wavenumber bandwidth is detected to be too narrow, causing the denominator value to be less than the preset safety lower limit threshold (e.g., 10), the verification will be performed. -6 If the error occurs (indicating local loss of the light source spectrum or signal-to-noise ratio not meeting requirements), or if the aforementioned determination coefficients fail to meet the standards, the system will directly intercept the division command and mark the coordinate point as invalid, thereby blocking the propagation of division-to-zero anomalies and divergent noise. For the specific unfolding algorithm of phase unwrapping, those skilled in the art can use the well-known branching method or minimum norm method to implement it; their basic derivations are well-known techniques in this field and will not be elaborated upon here.

[0106] S604 is the output format for generating 3D topographic data of the measured object by integrating the absolute surface height data of the entire array. After the above point-by-point independent calculation process is completed, the system will obtain the discrete one-dimensional height scalar. Spatial stitching is performed in system memory. Strictly adhering to the physical pixel arrangement rules of the front-end hardware, the height values ​​of each coordinate point are assembled into a two-dimensional floating-point matrix with consistent resolution. .

[0107] For discrete coordinate points marked as invalid in the aforementioned verification mechanism, median filtering or inverse distance weighted interpolation algorithms based on the effective neighborhood of the perimeter are used for smooth filling. Finally, the system generates a standard point cloud data file or renders it as a pseudo-color depth map based on this two-dimensional matrix. At this point, the three-dimensional topography data containing the absolute height distribution of the top single surface of the measured object is reconstructed, and the output reflects the true physical shape after removing parasitic interference from the bottom thin film.

[0108] To enable those skilled in the art to better understand the present invention, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments. It should be noted that, unless otherwise specified, the embodiments and features in the embodiments of this application can be combined with each other. This embodiment uses a surface covered with transparent silicon oxide (… The application scenario is the measurement of the three-dimensional morphology of microsteps on semiconductor silicon substrates in thin films.

[0109] During the data acquisition phase, the underlying piezoelectric ceramic actuator performs a vertical scan along the Z-axis with a constant physical step size of 20 nm. A two-dimensional area array detector is synchronously triggered at each step node, acquiring a spatial interferometric pattern with a resolution of 1024×1024. Because the silicon substrate surface is covered with approximately... Thick transparent The thin film and the broadband beam are reflected at the upper and lower surfaces of the thin film respectively. The time-domain interference sequence acquired by the system is a mixed state of double envelope superposition.

[0110] During parameter initialization and seed extraction, the system traverses the global contrast matrix, identifies the coordinates (512, 512) in the flat region at the center of the field of view as the pixel with the maximum contrast, and sets it as the starting point for seed optimization. Based on the extracted initial frequency domain phase, the system calculates the spatial phase gradient matrix of adjacent physical pixels. In the flat region, this gradient value remains low, approaching 0; while at the physical edges of the microsteps, the gradient value exhibits spatial jumps.

[0111] To address the aliasing signals of seed pixels, the system introduces fractional Fourier transform for time-frequency projection. After locking the candidate angle interval in the coarse optimization stage, the system calls a high-resolution local thinning algorithm based on linear frequency modulated Z-transform, effectively amplifying the local spectral resolution of the fractional domain by a factor of 100 without increasing the original mechanical scanning layer count. Combining the weighted phase curvature variance and the absolute deviation of the zero-crossing point in the time domain, the system constructs a joint decoupling objective function, ultimately achieving the desired effect at the rotation angle. The global minimum value is obtained at the current coordinate, which precisely locks the optimal decoupling parameters for the current coordinate.

[0112] Under normal, smooth operating conditions: the system uses the seed pixel as the origin and utilizes the eight-neighbor connectivity to diffuse the calculation layer by layer outwards from the periphery of the two-dimensional array. Since most regions have smooth spatial topography and uniform thin films, the spatial phase gradient is extremely small. Therefore, the system assigns an extremely narrow dynamic search radius to adjacent pixels based on a nonlinear exponential mapping model (e.g., ...). Based on this adaptive constraint boundary, the system directly skips a large-scale blind search step, significantly reducing computing power consumption.

[0113] Simulation of abrupt changes in morphology: When the diffusion calculation path extends to the edge of a microstep with a 500nm height difference in the substrate, the sudden change in the underlying physical morphology leads to a change in the spatial phase gradient matrix. The corresponding elements in the dataset experience a sudden, precipitous increase. Based on a dynamic feedback mechanism, the system automatically widens the search radius tolerance of the pixel to be calculated at this edge to the maximum limit. The local refinement interval is reset based on the angle of adjacent known effective pixels. Through this boundary adjustment action, the system successfully overcomes the discrete phase abrupt change region caused by topographic steps, and mathematically avoids the loss of the optimal decoupling angle and the cross-regional propagation of divergence error.

[0114] After the optimal rotation angle matrix of the entire array is constructed, the system substitutes the data of each coordinate node into the calculation channel and generates an adaptive graph basis narrowband filter window in the specified fractional domain to filter out the data generated by the matrix. The parasitic interference energy peak introduced by the reflection from the thin film surface. Subsequently, the system performs an inverse projection transformation with opposite rotation angles to restore the pure features to the normal time domain, and reconstructs the true three-dimensional height distribution matrix of the silicon substrate microsteps after completely eliminating thin film coupling interference by solving the absolute slope of the phase-wavenumber curve.

[0115] Comparison group settings:

[0116] Control group A: The traditional frequency domain analysis method (FDA) is used to directly perform fast Fourier transform and linear fitting to extract the height of the aliased interference signal, but it lacks decoupling and separation logic for transparent film signals.

[0117] Control group B: The global high-resolution fractional Fourier transform blind search method is adopted. The calculation process does not introduce the local amplification mechanism of linear frequency modulated Z-transform, and there is no spatial gradient adaptive boundary constraint when the two-dimensional array is extended.

[0118] Experimental group (this invention): Decoupling is achieved by using a joint objective function of time-frequency dual-domain invariants. The underlying calculation introduces linear frequency-modulated Z-transform to achieve resolution scaling. The array optimization coupling is based on the dynamic boundary propagation path of spatial phase gradient.

[0119] The experimental data are shown in Table 1:

[0120] Table 1 Comparison of Measurement Accuracy and Solving Efficiency of Microsteps in Complex Thin Films Table 1 Comparison of Measurement Accuracy and Solving Efficiency of Microsteps in Complex Thin Films

[0121]

[0122] in conclusion:

[0123] According to Table 1 and Figure 2 In complex microscopic surface measurement scenarios with transparent thin films, the traditional FDA solution (control group A) suffers from severe distortion in its high-resolution calculation due to its inability to isolate the phase superposition interference of parasitic interference components. The measurement error deteriorates exponentially as the thin film thins (leading to severe coherent superposition). While the global blind search solution (control group B) can theoretically achieve signal decoupling to some extent, its computational cost is extremely high when dealing with megapixel arrays, with a single frame calculation time exceeding 5 minutes, failing to meet the time-sensitive requirements of industrial-grade online inspection. This invention constructs an adaptive optimization boundary based on physical invariant features, combined with underlying linear frequency modulation scaling operations, to strictly limit the matrix computational cost of high-resolution decoupling to a local minimum neighborhood. While ensuring absolute topography reconstruction accuracy at the nanometer level (with errors controlled within 2.5nm even under extreme aliasing), the computation time of the overall two-dimensional array was reduced by about 96%. This successfully broke through the technical bottleneck of high precision and high computational load in thin-film interferometric aliasing measurement, providing a highly confident and feasible technical means for three-dimensional topography measurement of semiconductor wafers and microelectromechanical systems (MEMS).

[0124] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A method for measuring the three-dimensional topography of an object surface based on optical interferometry, characterized in that, include: S10, acquire the mixed temporal interference signal sequence covering the coordinates of each pixel in the two-dimensional array during the axial scanning process of the object under test; S20, extract seed pixels and construct a spatial topological connectivity feature matrix that characterizes the spatial physical continuity of the surface of the object under test; S30 introduces fractional Fourier transform to construct a joint decoupling objective function that includes frequency domain nonlinear phase distortion evaluation index and time domain physical invariant evaluation index; S40, for the seed pixel, perform multi-level optimization with the minimization of the joint decoupling objective function as the convergence condition, and lock the global optimal rotation angle of the seed pixel; S50, starting with the seed pixel, perform array adaptive decoupling in combination with the spatial topological connectivity feature matrix to generate the optimal rotation angle matrix corresponding to each pixel; S60, based on the optimal rotation angle matrix, perform signal decoupling and reconstruction on the mixed temporal interference signal sequence of each pixel, and extract the absolute surface height of each pixel to generate the three-dimensional shape data of the object under test.

2. The method for measuring the three-dimensional shape of an object surface based on optical interferometry according to claim 1, characterized in that, In step S10, the process of obtaining the hybrid time-domain interference signal sequence is as follows: The control actuator performs a physical scan with a preset step length along the axis. At each scanning step node, the detector is triggered to collect spatial interference patterns. The interference patterns at each scanning position are stacked according to the axial coordinates, and the light intensity values ​​of the same pixel at different axial positions are extracted to form a one-dimensional data vector.

3. The method for measuring the three-dimensional topography of an object surface based on optical interferometry according to claim 1, characterized in that, Step S20, the process of constructing the spatial topological connectivity feature matrix includes: Calculate the contrast of the interference fringes of the entire array of pixels, and select the pixel with the largest contrast as the seed pixel; perform a fast Fourier transform on the hybrid time-domain interference signal sequence to extract the initial frequency domain phase, and calculate the spatial phase gradient between adjacent pixels to obtain a first-order spatial phase gradient modulus array that characterizes the degree of drastic change in the morphology of the measured surface.

4. The method for measuring the three-dimensional topography of an object surface based on optical interferometry according to claim 1, characterized in that, In step S30, the extraction process of the frequency domain nonlinear phase distortion evaluation index includes: The mixed time-domain interference signal sequence is projected onto the fractional Fourier domain corresponding to a specific rotation angle. After energy filtering, it is restored to the conventional wavenumber domain to obtain the frequency domain phase sequence. The second derivative of the phase is calculated using a differential filter to obtain the phase curvature. The weighted phase curvature variance is calculated with power spectral density as the weight, which serves as an indicator to quantify the degree of nonlinear distortion introduced by thin-film crosstalk.

5. The method for measuring the three-dimensional shape of an object surface based on optical interferometry according to claim 1, characterized in that, In step S30, the extraction process of the time-domain physical invariant evaluation index includes: The filtered signal in the fractional Fourier domain is subjected to an inverse transform to restore it to a reconstructed time-domain signal. The peak coordinates of the interference envelope and the zero-phase point coordinates of the carrier are extracted from the reconstructed time-domain signal, and the axial physical displacement deviation between the two is calculated as an index for quantifying the residual asymmetric distortion in the time domain.

6. The method for measuring the three-dimensional morphology of an object surface based on optical interferometry according to claim 5, characterized in that, In step S30, the construction process of the joint decoupling objective function includes: The weighted phase curvature variance and the axial physical displacement deviation are assigned normalized weighting coefficients and summed. The normalized weighting coefficients are dynamically determined based on the reciprocal of the statistical mean of the two indicators in the parameter space.

7. The method for measuring the three-dimensional shape of an object surface based on optical interferometry according to claim 1, characterized in that, In step S40, the multi-level optimization includes: Perform a global coarse discrete step search to locate coarse search candidate angles; A linear frequency modulated Z-transform is introduced within the neighborhood defined by the coarse optimization candidate angle. By mapping the complex plane calculation trajectory to a narrow band spiral in the target frequency band, the spectrum is locally amplified, and high-resolution refinement optimization is performed to lock the globally optimal rotation angle.

8. The method for measuring the three-dimensional shape of an object surface based on optical interferometry according to claim 1, characterized in that, In step S50, the process of performing array adaptive decoupling includes: Plan a layer-by-layer traversal path that radiates outwards; use a nonlinear exponential saturation growth model to map the element values ​​in the spatial topological connectivity feature matrix to the dynamic search radius of the pixel to be calculated. Within the boundary defined by the optimal rotation angle of the reference pixel and the dynamic search radius, the local refinement optimization algorithm is invoked to determine the optimal rotation angle of the current pixel to be calculated.

9. The method for measuring the three-dimensional topography of an object surface based on optical interferometry according to claim 1, characterized in that, In step S60, the reconstruction process includes: The optimal rotation angle corresponding to each pixel is used to perform a fractional Fourier transform on the signal. In the fractional domain, a narrow-band filter window with a cosine roll-off coefficient is used to suppress crosstalk components. Then, an inverse fractional Fourier transform with opposite rotation angles is performed to recover the pure time-domain interference signal of a single surface.

10. The method for measuring the three-dimensional shape of an object surface based on optical interferometry according to claim 1, characterized in that, In step S60, the process of extracting the absolute surface height of each pixel includes: A continuous linear phase sequence with respect to the wavenumber is extracted from the pure time-domain interferometric signal. The first-order partial derivative of the continuous linear phase sequence with respect to the wavenumber is calculated using the least squares method to obtain the phase-wavenumber slope. Based on the phase-wavenumber slope, the axial physical height of each pixel relative to the zero optical path difference position is calculated.

Images