Wavefront restoration method based on multi-path integration

Abstract

The invention discloses a wavefront restoration method based on multi-path integration, and relates to the technical field of wavefront measurement and restoration, and the method comprises the following specific steps: data collection and initialization: carrying out the preliminary detection of wavefront through a Hartmann sensor microlens array, dividing the space where the wavefront is located into a plurality of rectangular small regions, according to the method, a multi-path calculation phase evaluation value set is introduced, all shortest paths from an initial point to each to-be-estimated point are found by using a shortest path search algorithm, it is ensured that multiple possibilities of wavefront propagation can be considered in the calculation process, and by introducing a path weight distribution mechanism, the estimation accuracy of the wavefront propagation is improved. According to the method, geometric characteristics, optical characteristics and correlation with surrounding paths of paths are comprehensively considered, each path is endowed with different weights, forward integration is carried out on discrete data results by adopting an innovative mixed integration algorithm, local detail loss or error amplification is effectively avoided, and thus a more accurate wavefront phase value is obtained.

Description

Technical Field

[0001] The present invention relates to the technical field of wavefront measurement and restoration, and particularly to a wavefront restoration method based on multi-path integration. Background Art

[0002] Wavefront measurement and restoration technology is a key technology in the field of optical engineering, and is widely used in fields such as astronomical observation, adaptive optical systems, optical imaging, and laser engineering. Its core lies in measuring parameters such as the slope and phase of the wavefront, and then reconstructing the complete form of the wavefront, providing key data support for the performance optimization of optical systems. The wavefront restoration technology not only requires high precision, but also needs to have good real-time performance and robustness to cope with complex and changing optical environments and system requirements. With the rapid development of optical technology, the requirements for wavefront restoration technology are also getting higher and higher, especially in application scenarios with high precision and large data volume. Traditional wavefront restoration methods face many challenges.

[0003] Traditional regional methods, such as algorithms related to the Southwell model, although have achieved certain results in wavefront reconstruction, their limitations are becoming increasingly prominent. When calculating the phase evaluation value of sub-apertures, such methods usually calculate the phase value of the point to be estimated only through a single specific path, which results in limited calculation accuracy. In the case of large data volume, the traditional least squares fitting method has high requirements for the computing performance of the computer, not only increasing the calculation cost, but also restricting the wide application of the algorithm in complex optical systems. In addition, when traditional regional methods process wavefront slope data, they do not fully utilize the correlation between sub-apertures, resulting in the reconstruction accuracy being easily affected in areas with large wavefront curvature changes. These problems seriously restrict the application effect of wavefront restoration technology in high-precision and large-data-volume scenarios.

[0004] In view of the above problems, it is necessary to optimize an existing wavefront restoration method based on multi-path integration. By designing a multi-path mean algorithm, the accuracy of wavefront reconstruction by a Hartmann sensor is improved, and the dependence on computer performance is reduced to a certain extent. Therefore, it is of great significance to develop a wavefront restoration method based on multi-path integration that can comprehensively achieve the above characteristics. Summary of the Invention

[0005] The objective of the present invention is to make up for the deficiencies of the prior art and provide a wavefront restoration method based on multi-path integration. It can achieve precise acquisition and efficient processing of wavefront slope data through an adaptive sub-aperture screening strategy, effectively reducing data redundancy and computational complexity. At the same time, by constructing a sub-aperture network topology data structure and applying the shortest path search algorithm and path weight assignment mechanism, precise calculation of the wavefront phase evaluation value is realized. During the process of obtaining the precise phase value through forward integration, a hybrid integration algorithm combining piecewise integration based on local data characteristics and overall weighted averaging is adopted to further improve the calculation accuracy of the phase value. In addition, the present invention also integrates the wave aberration calculation methods based on the geometric optical model and the wave optical model to achieve comprehensive and precise calculation of wave aberration information.

[0006] To solve the above technical problems, the present invention provides the following technical solution: A wavefront restoration method based on multi-path integration, the method comprising the following specific steps:

[0007] Data acquisition and initialization: Use the micro-lens array of the Hartmann sensor to preliminarily detect the wavefront, divide the space where the wavefront is located into multiple small rectangular regions, for each divided small region, select multiple sampling points, calculate the curvature of the region by using a fitting algorithm by detecting the difference in wavefront height between the sampling points and set a curvature threshold. For regions with curvature values greater than the threshold, increase the sampling quantity, and for regions with curvature values less than or equal to the threshold, reduce the sampling quantity. Integrate the wavefront average slope data corresponding to each sub-aperture collected to construct a vector matrix S, and it is stipulated that the slope data in the X direction is arranged first and then the slope data in the Y direction. At the same time, based on the detection of the wavefront energy distribution, determine the sub-aperture in the energy concentration region or adjacent region as the initial point with coordinates (0, 0) to provide a data basis and starting reference for subsequent calculations;

[0008] Multi - path calculation of phase evaluation value set: Construct a data structure of the connection relationship between sub - apertures in the Hartmann sensor microlens array in the form of an adjacency matrix or adjacency list. Use the shortest - path search algorithm based on the constructed data structure to search for all the shortest paths from the initial point (0, 0) to each point to be estimated (i, j). During the search process, record the node sequence information of each shortest path. Based on the comprehensive evaluation of the geometric characteristics of the path, the optical characteristics of the sub - apertures passed through, and the correlation with the surrounding paths, assign different weights to each path. For each found shortest path, calculate the phase evaluation values of adjacent sub - apertures in sequence along the weighted path based on the Southwell model. During the calculation process, according to the path order, calculate the phase evaluation value of the current sub - aperture based on the calculation result of the previous sub - aperture and the slope data corresponding to the current sub - aperture until the calculation of all sub - apertures on the entire path is completed, obtaining a phase evaluation value result for the point to be estimated under this path, and repeat this calculation process for all found shortest paths, thereby obtaining the set of phase evaluation value results for the point to be estimated;

[0009] Forward integration to obtain the phase value: According to the local change trend of the wavefront, divide the result set into multiple sub - regions according to the change rate of the phase difference between adjacent data points. Use the numerical integration method to calculate the integration result within the sub - region, determine its importance weight by combining the area ratio and energy contribution of the sub - region, and obtain the wavefront phase value of the point to be estimated through weighted average;

[0010] Calculation of wave aberration information: Use the geometric - optics model to estimate the approximate range and main characteristics of the wave aberration based on the phase value and the slope matrix through the ray - tracing principle, and perform precise correction calculations on the preliminary estimation results considering diffraction and interference phenomena at the wavefront edge and local details based on the wave - optics model. Integrate the two results to obtain complete and accurate wave aberration information to achieve precise wavefront reconstruction.

[0011] Furthermore, in the data acquisition and initialization step, by detecting the difference in wavefront height between these sampling points, use the fitting algorithm to calculate the curvature of this region. The algorithm formula is: where, C r represents the curvature value of a certain region divided, which is used to measure the degree of wavefront curvature in this region and is a key index determining the sampling density of sub - apertures in this region. W(x i , y j ) represents the wavefront height value at the sampling point with coordinates (x i , y j ) in this region. i and j respectively represent the sampling - point serial numbers in the x - direction and y - direction, and n and m are the number of sampling points in the x - direction and y - direction respectively, and respectively represent the wavefront height values W(x i , yj )Second-order partial derivatives with respect to x and y.

[0012] Furthermore, in the data acquisition and initialization step, the sub-aperture in the energy concentration region or adjacent region is determined based on the wavefront energy distribution detection as the initial point with coordinates (0, 0), and its calculation formula is: where P init is the determined initial point coordinate, m is the number of regions participating in the calculation, E i is the total energy of the i-th region, P i is the phase consistency vector of the i-th region, and the calculation method is: where k is the number of sampling points in the region, is the phase value of the j-th sampling point, V j is the vector pointing from this sampling point to the center of the region, P i reflects the overall trend and consistency degree of the phase within the region, θ i is the angle between the energy vector and the phase consistency vector of the i-th region.

[0013] Furthermore, in the step of calculating the phase evaluation value set for multipath, the shortest path search algorithm is used to search for all shortest paths from the initial point (0, 0) to each point to be estimated (i, j) based on the constructed data structure. Specifically, an optical coherence function is defined to measure the degree of coherence maintained during light propagation between two sub-apertures, and its calculation formula is: where C mn represents the optical coherence degree between sub-aperture m and sub-aperture n, ω is the coherence reference coefficient, β is the coherence attenuation coefficient, represents the phase difference between sub-aperture m and sub-aperture n, ρ nn is the optical topological distance deviation between sub-aperture m and sub-aperture n, d mn is the physical space distance between sub-aperture m and sub-aperture n, α is the distance influence index, and the coherence weight of the path is defined based on the optical coherence function. For any path P = {s, v 1 , v 2 , …, v k-1 , t} from the initial point s(0, 0) to the point to be estimated t(i, j), where v i represents the intermediate sub-aperture passed on the path, and its path coherence weight is defined as the product of the optical coherence degrees between all adjacent sub-apertures on the path, that is: where represents sub-aperture v i and sub-aperture v i+1The optical coherence degree between them. At the same time, considering the connectivity and complexity of the sub-aperture network topology itself, a topology weight function is introduced. For path P, its topology weight is defined as: Where is a measure of the complexity of the connection between sub-aperture v i and v i+1 based on the topology. Based on the coherence weight and topology weight of the path, the comprehensive weight of the path is defined as: For all the shortest paths to be searched from the initial point (0, 0) to the point to be estimated (i, j), let the set of all possible paths from the initial point s to the point to be estimated t be Π st , then the shortest path set P shortest (s, t) is determined by the formula , that is, among all the paths from the initial point s to the point to be estimated t, the path comprehensive weight is greater than or equal to the comprehensive weight of any other path P' in the set . The paths form the set of all the shortest paths from the initial point (0, 0) to the point to be estimated (i, j).

[0014] Furthermore, in the step of calculating the phase evaluation value set of multiple paths, based on the comprehensive evaluation of the geometric characteristics of the path, the optical characteristics of the sub-apertures passed through, and the correlation with the surrounding paths, different weights are assigned to each path. Its comprehensive weight formula is: W(P) = W G (P)·W O (P)·W R (P), where W G (P) is the weight factor based on geometric characteristics, W O (P) is the weight factor based on the optical characteristics of the sub-apertures passed through, and W R (P) is the weight factor based on the correlation with the surrounding paths. For the weight factor based on geometric characteristics, its calculation formula is: Where L P is the sum of the distances between all adjacent sub-apertures passed through on the path, α is the length influence coefficient, β is the bending degree influence coefficient, is the average angle change value. For the weight factor based on the optical characteristics of the sub-apertures passed through, let path P pass through m sub-apertures. For the k-th sub-aperture, its optical characteristic index is obtained through a special optical detection method and is represented by O k , and the weight factor based on the optical characteristics of the sub-apertures passed through is defined as: For the weight factor based on the correlation with the surrounding paths, by calculating the total length L overlap of the overlapping line segments of path P and other paths and the total length L P of path P, that is Define the weight factor based on the correlation with the surrounding paths as follows: where γ is the correlation influence coefficient.

[0015] Furthermore, in the step of calculating the set of phase evaluation values for multiple paths, for each found shortest path, based on the Southwell model, the phase evaluation values of adjacent sub-apertures are sequentially calculated along the weighted path. During the calculation process, according to the path order, based on the calculation result of the previous sub-aperture and the slope data corresponding to the current sub-aperture, the phase evaluation value of the current sub-aperture is calculated. Its calculation formula is: where i and j represent the sub-aperture numbers, S represents the slope, P represents the phase, and h is a constant. Due to discrete sampling, the formula can be written in matrix form: S = AP, where S is the vector matrix containing all slope measurement data, with the x-direction slope first and the y-direction slope second, P is the vector matrix containing all phase evaluation points, and A is the coefficient matrix.

[0016] Furthermore, in the step of obtaining the phase value by forward integration, the numerical integration method is used to calculate the integration result within the sub-region, and the importance weight is determined by combining the area ratio of the sub-region and the energy contribution. The wavefront phase value of the point to be estimated is obtained through weighted averaging. Specifically, let the set of phase evaluation value results be N is the number of phase evaluation values, and the corresponding spatial position coordinates are (x 1 , y 1 ), (x 2 , y 2 ), …, (x N , y N ). For two adjacent phase evaluation values and their phase change rate is Set the wavefront local change threshold ε. When , this place is determined as the boundary of the sub-region, thus realizing the division of the sub-region. Let the number of divided sub-regions be M, and the set of phase evaluation values included in the j-th sub-region be For each divided sub-region, the numerical integration method is used to calculate its integral value. For the j-th sub-region, its integral value I j The calculation formula is: where Δs j，k represents the spatial interval corresponding to two adjacent phase evaluation values within the j-th sub-region.

[0017] Furthermore, in the step of obtaining the phase value by forward integration, the importance weight is determined by combining the area ratio of the sub-region and the energy contribution. The wavefront phase value of the point to be estimated is obtained through weighted averaging. Specifically, let the total area corresponding to the entire wavefront be S total , and the area corresponding to the j-th sub-region be Sj , then the weight factor \(w\) of the sub-region area ratio S，j is: For the \(j\)-th sub-region, let the relationship function between the optical field energy and the phase be Then the energy \(E\) of the \(j\)-th sub-region j can be expressed as: Then the weight factor \(w\) of the sub-region energy contribution E，j is: Combining the two weight factors of the sub-region area ratio and the energy contribution, the importance weight \(w\) of the sub-region \(j\) is obtained j , and the calculation formula is Based on the obtained weights, the wavefront phase value of the point to be estimated is obtained by weighted averaging the integration results of each sub-region The calculation formula is:

[0018] Furthermore, in the wave aberration information calculation step, the geometric optical model is used to estimate the approximate range and main characteristics of the wave aberration according to the phase value and the slope matrix through the ray tracing principle. Specifically, for the preliminary calculation of the wave aberration of the geometric optical model, let the sub-aperture coordinates be \((x\) i , \(y\) i ), and the corresponding phase value is The slope value in the \(x\) direction is In the \(y\) direction is According to the slope value, calculate the direction vector of the ray at the sub-aperture. Then the unit direction vector component of the ray in the \(x\) direction and the unit direction vector component in the \(y\) direction can be expressed as: For the ray propagation path deviation, let the ray start from the reference plane and the path deviation generated after passing through the sub-aperture in the \(x\) direction be \(\Delta x\) i , and in the \(y\) direction be \(\Delta y\) i , and calculate through the formula where \(\Delta z\) j represents the distance that the ray propagates between adjacent sub-apertures. Combine the path deviations at all sub-apertures to obtain the wave aberration vector based on the geometric optical model

[0019] Furthermore, in the wave aberration information calculation step, based on the wave optical model, the preliminary estimation results are accurately corrected by considering diffraction and interference phenomena at the wavefront edge and local details. According to the calculated wave aberration vector and the initial light source conditions, construct the wavefront complex amplitude distribution \(U\) geo (x, y), and its calculation formula is: where \(A(x, y)\) is the amplitude distribution, For the phase distribution estimated based on geometric optics, substituting it into the correction formula gives the corrected complex amplitude distribution U wave (x′, y′), that is

[0020] where the wave number λ is the wavelength of light, i is the imaginary unit, z represents the axial distance between the wavefront plane and the observation plane. Extract the phase information from the corrected complex amplitude distribution U wane (x′, y′) Calculate the exact wave aberration vector under the wave optics model according to the phase information Thus, the geometric optics estimation results are accurately corrected at the wavefront edge and local details.

[0021] Compared with the prior art, the wavefront restoration method based on multi-path integration has the following beneficial effects:

[0022] First, by introducing a multi-path calculation of the phase evaluation value set and using the shortest path search algorithm to find all the shortest paths from the initial point to each point to be estimated, it ensures that various possibilities of wavefront propagation can be considered during the calculation process. And by introducing a path weight assignment mechanism, comprehensively considering the geometric characteristics, optical properties of the path and its correlation with the surrounding paths, different weights are assigned to each path. An innovative hybrid integration algorithm is used to perform forward integration on the discrete data results, effectively avoiding the loss of local details or error amplification, thereby obtaining a more accurate wavefront phase value.

[0023] Second, through the combination of the adaptive sub-aperture screening strategy and the multi-path integration algorithm, the calculation complexity is effectively reduced. The adaptive sub-aperture screening strategy dynamically adjusts the selection range and density of the sub-apertures according to the local characteristics of the wavefront, optimizes the data acquisition volume, reduces unnecessary data redundancy. At the same time, the multi-path integration algorithm only calculates the shortest path and combines the weight assignment mechanism, avoiding the complex matrix operations and a large amount of data processing in the traditional algorithm. It not only reduces the calculation complexity, but also improves the practicality and scalability of the algorithm.

[0024] Other advantages, objectives and features of the present invention will be described to some extent in the subsequent specification, and to some extent, will be obvious to those skilled in the art based on the study of the following text, or can be taught from the practice of the present invention. BRIEF DESCRIPTION OF THE DRAWINGS

[0025] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the following will briefly introduce the drawings required for the description of the embodiments or the prior art. Obviously, the drawings in the following description are only some embodiments of the present invention. For those of ordinary skill in the art, without creative efforts, other drawings can also be obtained based on these drawings.

[0026] Figure 1 It is a flow operation diagram of a wavefront restoration method based on multi-path integration;

[0027] Figure 2 It is a flow chart of a wavefront restoration method based on multi-path integration. Specific embodiments

[0028] To further elaborate on the technical means and effects adopted by the present invention to achieve the predetermined invention purpose, the following will, in conjunction with the drawings and preferred embodiments, detail the specific embodiments, structures, features and their effects of the present invention as follows.

[0029] Embodiment 1

[0030] Apply the wavefront restoration method based on multi-path integration of the present invention in the adaptive optical system of an astronomical telescope to improve the clarity and resolution of telescope imaging.

[0031] Start the Hartmann sensor to detect the wavefront corresponding to the starlight brightness received by the telescope. Divide the space where the wavefront is located into many small square regions with equal areas (for example, square regions with a side length of 1 mm). For each region, select multiple sampling points. By analyzing the wavefront height differences between the sampling points, use a specific fitting algorithm to calculate the curvature of each region. The algorithm formula is: Among them, C r represents the curvature value of a certain region divided, which is used to measure the bending degree of the wavefront in this region and is a key index determining the sampling density of sub-apertures in this region. W(x i , y j ) represents the wavefront height value at the sampling point with coordinates (x i , y j ) in this region. i and j respectively represent the sampling point numbers in the x direction and the y direction, and n and m are the numbers of sampling points in the x direction and the y direction respectively. and respectively represent the wavefront height values W(x i , y j)For the second-order partial derivatives with respect to x and y, set a curvature threshold. For regions where the curvature is greater than this threshold, it means that the wavefront changes drastically here. For example, in some regions at the edge of the telescope's field of view, due to factors such as oblique incidence of light, the wavefront bends and changes greatly, so the sampling density is increased. For relatively flat regions like the center of the field of view (where the curvature is less than the threshold), the sampling is appropriately reduced. In this way, the wavefront average slope data corresponding to each sub-aperture is collected and organized into a vector matrix S, with the slope data in the X direction in the front and the slope data in the Y direction in the back. At the same time, using a wavefront energy distribution detection device, the entire wavefront of the telescope is divided into multiple energy detection regions in the form of concentric rings, and the average energy value in each region is calculated separately. It is found that the energy is most concentrated in the central ring region. In this central region or in an adjacent sub-aperture, a sub-aperture is selected as the initial point (0, 0) to prepare for subsequent path calculation.

[0032] Construct the connection relationship data structure between sub-apertures in the Hartmann sensor microlens array in the form of an adjacency matrix. The rows and columns of the matrix correspond to each sub-aperture. If two sub-apertures are adjacent, the corresponding position is assigned a value of 1, indicating connectivity. If they are not adjacent, the value is assigned 0. In this way, the adjacent relationship between sub-apertures is clearly presented, facilitating subsequent path search. Based on the above constructed data structure, start searching for all the shortest paths from the initial point (0, 0) to each point to be estimated (such as sub-apertures at different positions like the point to be estimated (5, 5), etc.) through an algorithm. Specifically, define an optical coherence function to measure the degree of coherence maintained during light propagation between two sub-apertures. Its calculation formula is: where C mn represents the optical coherence between sub-aperture m and sub-aperture n, ω is the coherence reference coefficient, β is the coherence attenuation coefficient, represents the phase difference between sub-aperture m and sub-aperture n, ρ nn is the optical topological distance deviation between sub-aperture m and sub-aperture n, d mn is the physical space distance between sub-aperture m and sub-aperture n, and α is the distance influence index. Define the coherence weight of the path based on the optical coherence function. For any path P = {s, v 1 , v 2 , …, v k-1 , t} from the initial point s(0, 0) to the point to be estimated t(i, j), where v i represents the intermediate sub-aperture passed on the path, and its path coherence weight is defined as the product of the optical coherence between all adjacent sub-apertures on the path, that is: where represents the optical coherence between sub-aperture v i and sub-aperture v i+1The optical coherence degree between them. At the same time, considering the connectivity and complexity of the sub-aperture network topology itself, a topological structure weight function is introduced. For path P, its topological structure weight is defined as: where is a measure of the complexity of the connection between sub-aperture v i and v i+1 based on the topological structure. Based on the coherence weight and topological structure weight of the path, the comprehensive weight of the path is defined as: For all the shortest paths to be searched from the initial point (0, 0) to the point to be estimated (i, j), let the set of all possible paths from the initial point s to the point to be estimated t be Π st , then the shortest path set P shortest (s, t) is determined by the formula , that is, among all the paths from the initial point s to the point to be estimated t, the path with the comprehensive weight greater than or equal to the comprehensive weight of any other path P' in the set constitutes all the shortest path sets from the initial point (0, 0) to the point to be estimated (i, j). During the search process, the sequence information of the sub-aperture nodes passed by each shortest path is recorded in detail, such as the sub-aperture numbers passed in sequence starting from the initial point.

[0033] For each found shortest path, its weight is comprehensively evaluated. From the geometric feature aspect, the path length and curvature are calculated. For example, a more tortuous and relatively long path may accumulate more errors during light propagation. According to the preset geometric feature weight calculation rule, it will be given a relatively low weight; from the optical characteristics of the sub-apertures passed, by detecting indicators such as the optical medium uniformity at the sub-apertures, if a path passes through sub-apertures with high uniformity and low light absorption, its weight in this regard will be higher. Then considering the correlation with surrounding paths, the overlap situation between the path and other paths is statistically analyzed. The weight of a path with more overlap and greater possibility of being interfered is correspondingly reduced. Considering these three aspects of factors, according to the established weight distribution algorithm W(P) = W G (P)·W O (P)·W R (P), a suitable weight value is determined for each shortest path.

[0034] Calculating the phase evaluation value based on the Southwell model: Based on the relationship between the sub-aperture slope and the phase described by the Southwell model, along each weighted shortest path, starting from the initial point, according to the order of the sub-apertures on the path, using the phase evaluation value of the previous sub-aperture and the slope data corresponding to the current sub-aperture (obtained from the vector matrix S), the phase evaluation values of adjacent sub-apertures are calculated in sequence until the phase evaluation value of the point to be estimated is calculated. The calculation formula is: Where i and j represent the sub-aperture numbers, S represents the slope, P represents the phase, and h is a constant. Due to discrete sampling, the formula can be written in matrix form: S = AP, where S is a vector matrix containing all slope measurement data, with the x-direction slope first and the y-direction slope second, P is a vector matrix containing all phase evaluation points, and A is the coefficient matrix. For each point to be estimated, the above operations are repeated, and the phase evaluation values obtained from all the shortest paths are aggregated to form a set of phase evaluation values for the point to be estimated.

[0035] Analyze the obtained set of discrete phase evaluation values to check the change in phase values between adjacent data points. When the change rate of the phase difference between adjacent data points exceeds the set threshold, this location is determined as the boundary of the sub-region, and the entire result set is divided into multiple sub-regions. For example, some phase evaluation values near the edge of the telescope optical correction element may be divided into relatively more and smaller sub-regions due to their relatively frequent changes, while relatively larger sub-regions are divided in relatively stable regions. Specifically, let the set of phase evaluation values be N is the number of phase evaluation values, and the corresponding spatial position coordinates are (x 1 , y 1 ), (x 2 , x 2 ), …, (x N , y N ). For two adjacent phase evaluation values and their phase change rate is Set the local wavefront change threshold ε. When , this location is determined as the boundary of the sub-region, thus realizing the division of the sub-region. Let the number of divided sub-regions be M, and the set of phase evaluation values contained in the j-th sub-region is For each divided sub-region, the numerical integration method is used to calculate its integral value. For the j-th sub-region, its integral value I j The calculation formula is: Where Δs j，k represents the spatial interval corresponding to two adjacent phase evaluation values within the j-th sub-region. Then, calculate the importance weight of each sub-region. From the perspective of area ratio, calculate the ratio of the area of each sub-region to the total area of the entire wavefront detection region. From the perspective of energy contribution, according to the pre-set energy calculation method, combine the phase evaluation values to calculate the energy situation of the sub-region, and then calculate its ratio to the total wavefront energy. Combine these two aspects and determine the importance weight of the sub-region according to a certain weighting rule. Finally, sum the integral results of each sub-region multiplied by their respective weights and divide by the sum of the weights of all sub-regions. Through such a weighted average operation, a more accurate wavefront phase value of the point to be estimated is obtained.

[0036] Based on the phase value vector matrix and wavefront slope vector matrix obtained from the previous calculations, using the principle of ray tracing, for each sub-aperture, determine the change in the outgoing direction of the ray at that sub-aperture according to its phase value and slope, simulate the change in the propagation path of the ray on the wavefront, and thus estimate the approximate range and main characteristics of the wave aberration. Specifically, for the preliminary calculation of the wave aberration of the geometric optical model, assume the sub-aperture coordinates are (x i , y i ), and the corresponding phase value is The slope value in the x direction is In the y direction is Calculate the direction vector of the ray at the sub-aperture according to the slope value. Then, the unit direction vector component of the ray in the x direction And the unit direction vector component in the y direction Can be expressed as: For the ray propagation path deviation, assume the ray starts from the reference plane, and the path deviation generated after passing through the sub-aperture in the x direction is Δx i , and in the y direction is Δy i , and calculate through the formula , where Δz j Represents the distance that the ray propagates between adjacent sub-apertures. Combine the path deviations at all sub-apertures to obtain the wave aberration vector based on the geometric optical model For example, if it is initially determined that the aberrations that may occur in telescope imaging mainly concentrate on the edge of the field of view and certain specific angular directions, using wave optical methods such as Fresnel diffraction integral, take the wavefront complex amplitude situation corresponding to the wave aberration estimation result obtained from the geometric optical model as the input, fully consider the diffraction, interference and other phenomena brought by the wave nature of light, and conduct fine analysis and calculation on the wavefront edge and local details. According to the calculated wave aberration vector And the initial light source conditions to construct the wavefront complex amplitude distribution U geo (x, y), and its calculation formula is: Where A(x, y) is the amplitude distribution, Is the phase distribution estimated based on geometric optics. Substitute it into the correction formula to obtain the corrected complex amplitude distribution U wave (x′, y′), that is

[0037] Among them, the wave number λ is the wavelength of light, i is the imaginary unit, and z represents the axial distance between the wavefront plane and the observation plane. Extract the phase information wane From the corrected complex amplitude distribution U Calculate the accurate wave aberration vector under the wave optical model according to the phase information Thus, precise corrections are made to the results of geometric-optics estimation at the wavefront edge and local details, such as correcting the wave aberration information of the blurred part caused by diffraction at the field edge to obtain more accurate wave aberration information.

[0038] According to the weights determined in advance according to the characteristics of the telescope optical system, the wave aberration information preliminarily estimated based on the geometric-optics model and the wave aberration information precisely corrected based on the wave-optics model are weighted and fused to form the final complete and accurate wave aberration information.

[0039] After calculating the accurate wave aberration information through the above wavefront restoration method based on multi-path integration, the control module in the adaptive optics system drives optical correction elements such as deformable mirrors to make corresponding adjustments based on this information. After a period of actual observation and comparison, it is found that the resolution of the telescope imaging is improved, and the clarity and contrast of the image are also significantly improved, effectively verifying the practicability and effectiveness of the present invention in the actual astronomical telescope adaptive optics system, and realizing the precise restoration of the wavefront and the improvement of the imaging quality.

[0040] Embodiment 2

[0041] The wavefront restoration method based on multi-path integration of the present invention is applied in an adaptive optics system for optical microscopy imaging, aiming to improve the quality of microscopic images and make the microscopic structure presented more clearly.

[0042] The Hartmann sensor is turned on to detect the object wavefront focused by the microscope objective. The microscopic field region where the wavefront is located is divided into several small rectangular regions. For each small region, multiple sampling points are selected to analyze the wavefront situation. By calculating the wavefront height difference between the sampling points, the curvature situation of each region is obtained using a suitable fitting algorithm, and a curvature threshold is set. For example, when observing a cell sample, the curvature of the wavefront region corresponding to the edge of structures such as cell nuclei in the sample changes greatly and exceeds the set curvature threshold, then the sampling density is increased to more finely capture the wavefront changes. In the relatively flat background region around the sample where the wavefront curvature changes little, unnecessary sampling is reduced. According to such a strategy, the wavefront average slope data corresponding to each sub-aperture are collected and constructed into a vector matrix S, where the slope data in the X direction are arranged in the front and the slope data in the Y direction are arranged in the back. At the same time, using a specially designed device for detecting the energy of the microscopic wavefront, the entire microscopic field is divided into 9 energy detection regions in the form of a nine-square grid, and the average energy value in each region is calculated. It is found that the energy distribution in the middle region is relatively concentrated due to factors such as light convergence. Therefore, a sub-aperture is selected as the initial point (0, 0) in this middle region or in the sub-aperture adjacent to it, which is convenient for subsequent operations such as path calculation based on this.

[0043] The connection relationship data structure between sub-apertures in the Hartmann sensor microlens array is constructed in the form of an adjacency list. Each structure corresponds to a sub-aperture. The structure contains the sub-aperture's own number information and a list of pointers pointing to the structures of its adjacent sub-apertures. Through this structure, the adjacent sub-aperture situation of each sub-aperture can be retrieved and traversed conveniently, providing efficient data support for subsequent search paths. Based on the constructed adjacency list data structure, the algorithm starts to search for all the shortest paths from the initial point (0, 0) to each point to be estimated (such as sub-apertures at different positions like the point to be estimated (3, 4), etc.). During the search, the sequence of sub-aperture nodes passed by each shortest path is carefully recorded, just like recording key information such as the order of sub-aperture numbers passed through successively starting from the initial point, for preparation for subsequent calculations.

[0044] For each found shortest path, its weight situation is comprehensively considered. In terms of geometric features, the length and bending degree of the path are calculated. For example, a path that winds a lot and is overall long is more likely to generate errors during light propagation. According to the preset geometric feature weight calculation rules, it will be given a lower weight. From the perspective of the optical properties of the sub-apertures passed through, indicators such as the optical medium uniformity and light scattering degree at the sub-aperture are obtained through optical detection means. For paths passing through sub-apertures with good optical properties (uniform medium, less light scattering), their weights in this regard will be correspondingly increased. For the correlation with surrounding paths, the length of overlapping line segments and the number of intersection points between this path and other paths are statistically analyzed. If there are many overlaps and intersections with many other paths, it means a high possibility of interference, then the weight of this path will be reduced. Considering the factors in these three dimensions, according to the established weight distribution mechanism, an appropriate weight value is assigned to each shortest path.

[0045] Based on the inherent relationship between the sub-aperture slope and the phase described by the Southwell model, along each shortest path that has been assigned a weight, starting from the initial point, in accordance with the order of sub-apertures on the path, the phase evaluation value of the previous sub-aperture and the slope data corresponding to the current sub-aperture (obtaining the corresponding data from the vector matrix S) are used to calculate the phase evaluation value of the adjacent sub-aperture one by one until the phase evaluation value of the point to be estimated is calculated. For each point to be estimated, the above calculation process is repeated, and then the phase evaluation values obtained through all the shortest paths are summarized to form the result set of the phase evaluation value of the point to be estimated.

[0046] The obtained set of discrete phase evaluation values is analyzed, with a focus on the amplitude of the phase value change between adjacent data points. When the change rate of the phase difference between adjacent data points exceeds a pre-set threshold, this location is taken as the boundary of the sub-region, and the entire result set is divided into multiple sub-regions. For example, in the wavefront region corresponding to some fine structures inside the cell, the phase change is relatively complex, and more relatively small sub-regions may be divided; while in the relatively uniform background region outside the cell, the divided sub-regions are relatively large and fewer in number.

[0047] Integral operations are carried out for each divided sub-region. Combining the phase evaluation values within the sub-region and information such as their corresponding spatial positions, the integral results of the sub-regions are calculated. Subsequently, the importance weights of each sub-region are calculated. From the perspective of area proportion, it is determined by analyzing the area ratio of the sub-region in the entire microscopic field of view region. From the perspective of energy contribution, based on a pre-set energy calculation method adapted to the characteristics of the microscopic imaging optical system, the energy situation of the sub-region is calculated in combination with the phase evaluation values, and then the ratio of its energy to the total energy of the entire wavefront is calculated. The importance weights of the sub-regions are determined according to the established weighting rules by comprehensively considering these two aspects. Through weighted average operations, a more accurate wavefront phase value of the point to be estimated is obtained.

[0048] According to the previously calculated phase value vector matrix and wavefront slope vector matrix, using the principle of ray tracing, for each sub-aperture in the microscope field of view, the change in the propagation direction of the light ray after passing through the sub-aperture is determined based on its phase value and slope, simulating the change in the propagation path of the light ray on the wavefront, and then a preliminary estimate of the approximate range and main characteristics of the wave aberration is made. For example, when observing a cell sample, it is initially judged that due to the non-uniform refractive index of the sample itself and the small deviations of the microscope optical elements, the wave aberration is mainly concentrated near the cell edge and the wavefront regions corresponding to some organelles. Using the calculation method based on wave optics, taking the wavefront complex amplitude distribution corresponding to the wave aberration estimation result obtained from the geometric optical model as the basic input information, fully considering physical phenomena such as diffraction and interference brought by the wave nature of light, precise analysis and calculation are carried out for the edge part of the wavefront and the local detail regions corresponding to the fine structures inside the cell, and the wave aberration information is carefully corrected so that the calculation result of the wave aberration can more accurately reflect the actual microscopic optical situation.

[0049] According to the weights determined in advance according to the characteristics of the microscope adaptive optical system, the wave aberration information preliminarily estimated based on the geometric optical model and the wave aberration information precisely corrected based on the wave optical model are weighted and fused to form the final complete and accurate wave aberration information.

[0050] With the above wavefront restoration method based on multi-path integration, after calculating the accurate wave aberration information, the controller in the adaptive optical system will adjust the parameters of optical correction elements such as liquid crystal spatial light modulators according to this information. Through actual microscopic imaging experiment comparison, it is found that the fine internal structure of cells after imaging is more clearly distinguishable, and the resolution of the image has been significantly improved compared with that without using this method, fully verifying the practicability and effectiveness of the present invention in the adaptive optical system of optical microscopy, and achieving the accurate restoration of the microscopic wavefront and the effective improvement of the imaging quality.

[0051] The above are only the preferred embodiments of the present invention, and do not impose any form of limitation on the present invention. Although the present invention has been disclosed above with the preferred embodiments, it is not intended to limit the present invention. Any person skilled in the art can make some changes or modifications to equivalent embodiments by using the above-disclosed technical content within the scope of the technical solution of the present invention. However, as long as it does not depart from the content of the technical solution of the present invention, any brief modifications, equivalent changes and modifications made to the above embodiments based on the technical essence of the present invention still fall within the scope of the technical solution of the present invention.

Claims

1. A wavefront restoration method based on multipath integration, characterized in that: The method comprises the following specific steps: Data collection and initialization: Use the Hartmann sensor microlens array to perform preliminary detection of the wavefront, divide the space where the wavefront is located into multiple rectangular small areas, select multiple sampling points from each divided small area, detect the difference in wavefront height between sampling points, use the fitting algorithm to calculate the curvature of the area and set the curvature threshold, increase the number of samples for areas where the curvature value is greater than the threshold, and reduce the number of samples for areas where the curvature value is less than or equal to the threshold, integrate the collected wavefront average slope data corresponding to each sub-aperture, construct a vector matrix S, and stipulate that the slope data in the X direction is arranged first, and then the slope data in the Y direction is arranged. At the same time, based on the wavefront energy distribution detection, the sub-aperture of the energy concentration area or the adjacent area is determined as the initial point with coordinates (0, 0), providing a data basis and starting reference for subsequent calculations; Multi-path calculation of phase evaluation value set: constructing a data structure of the connection relationship between subapertures in the Hartmann sensor microlens array in the form of an adjacency matrix or an adjacency list, and using the shortest path search algorithm based on the constructed data structure to search for all shortest paths from the initial point (0, 0) to each point to be estimated (i, j). During the search process, the node sequence information of each shortest path is recorded, and different weights are assigned to each path based on the geometric characteristics of the path, the optical characteristics of the subaperture, and the correlation with the surrounding paths. For each shortest path found, the phase evaluation values ​​of adjacent subapertures are calculated in sequence along the weighted path based on the Southwell model. During the calculation process, according to the path sequence, the phase evaluation value of the current subaperture is calculated based on the calculation result of the previous subaperture and the slope data corresponding to the current subaperture, until the calculation of all subapertures on the entire path is completed, and a phase evaluation value result of the point to be estimated under the path is obtained, and this calculation process is repeated for all the shortest paths found, thereby obtaining a phase evaluation value result set of the point to be estimated; Forward integration to obtain phase value: Based on the local change trend of the wavefront, the result set is divided into multiple sub-regions according to the change rate of the phase difference between adjacent data points. The integral result is calculated by numerical integration method in the sub-region, and its importance weight is determined by combining the area proportion and energy contribution of the sub-region. The wavefront phase value of the point to be estimated is obtained by weighted average; Calculation of wave aberration information: The geometric optics model is used to estimate the approximate range and main characteristics of the wave aberration based on the phase value and slope matrix through the principle of ray tracing. The preliminary estimation results are accurately corrected and calculated based on the wave optics model at the edge of the wavefront and local details taking into account diffraction and interference phenomena. The results of the two are integrated to obtain complete and accurate wave aberration information to achieve accurate wavefront reconstruction.

2. The wavefront restoration method based on multipath integration according to claim 1, characterized in that: In the data collection and initialization step, the difference in wavefront height between these sampling points is detected, and a fitting algorithm is used to calculate the curvature of the area. The algorithm formula is: Among them, C r Represents the curvature value of a certain area, which is used to measure the curvature of the wavefront in the area and is the key indicator for determining the sub-aperture sampling density in this area. i ,y j ) means that the coordinates in this region are (x i ,y j ), i and j represent the sampling point numbers in the x and y directions respectively, n and m represent the number of sampling points in the x and y directions respectively, and Respectively represent the wavefront height value W(x i ,y j ) with respect to x and y.

3. The wavefront restoration method based on multipath integration according to claim 1, characterized in that: In the data acquisition and initialization step, the sub-aperture of the energy concentration area or the adjacent area is determined as the initial point with coordinates (0, 0) based on the wavefront energy distribution detection, and the calculation formula is: Among them, P init is the coordinate of the initial point, m is the number of regions involved in the calculation, E i is the total energy of the ith region, P i is the phase consistency vector of the ith region and is calculated as: Where k is the number of sampling points in the region, is the phase value of the jth sampling point, V j is the vector pointing from the sampling point to the center of the region, P i It reflects the overall trend and consistency of the phase in the region, θ i is the angle between the energy vector and the phase consistency vector of the ith region.

4. The wavefront restoration method based on multipath integration according to claim 1, characterized in that: In the multipath calculation phase evaluation value set step, the shortest path search algorithm is used based on the constructed data structure to search for all shortest paths from the initial point (0, 0) to each point to be estimated (i, j). Specifically, an optical coherence function is defined to measure the degree of coherence maintained during light propagation between two sub-apertures, and its calculation formula is: Among them, C mn represents the optical coherence between subaperture m and subaperture n, ω is the coherence reference coefficient, β is the coherence attenuation coefficient, represents the phase difference between subaperture m and subaperture n, ρ nn is the optical topological distance deviation between subaperture m and subaperture n, d mn is the physical space distance between subaperture m and subaperture n, α is the distance influence index, and the coherence weight of the path is defined based on the optical coherence function. For any path P = {s, v1, v2, ..., v k-1 , t}, where v i represents the intermediate sub-aperture passed on the path, and its path coherence weight It is defined as the product of the optical coherence between all adjacent sub-apertures on the path, that is: in represents the subaperture v i With subaperture v i+1 At the same time, considering the connectivity and complexity of the sub-aperture network topology itself, a topology weight function is introduced. For path P, its topology weight is defined as: in, The subaperture v is based on the topological structure i With v i+1 The complexity measurement index of the connection between nodes is based on the coherence weight and topological structure weight of the path, and the comprehensive weight of the path is defined as: To search for all the shortest paths from the initial point (0, 0) to the point to be estimated (i, j), suppose all possible paths from the initial point s to the point to be estimated t constitute a set Π st , then the shortest path set P shortest (s, t) is expressed by the formula Determine, that is, in the set of all paths from the initial point s to the estimated point t, the path comprehensive weight Greater than or equal to the comprehensive weight of any other path P′ in the set The paths constitute the set of all shortest paths from the initial point (0, 0) to the point to be estimated (i, j).

5. The wavefront restoration method based on multipath integration according to claim 1, characterized in that: In the multipath calculation phase evaluation value set step, different weights are assigned to each path based on the geometric characteristics of the path, the optical characteristics of the sub-aperture, and the correlation with the surrounding paths. The comprehensive weight formula is: W(P) = W G (P)·W O (P)·W R (P), where W G (P) is the weight factor based on geometric features, W O (P) is a weighting factor based on the optical properties of the subaperture, W R (P) is a weight factor based on the correlation with the surrounding paths. For the weight factor based on geometric features, its calculation formula is: Among them, L P is the sum of the distances between all adjacent sub-apertures on the path, α is the length influence coefficient, β is the curvature influence coefficient, is the average angle change value. For the weight factor based on the optical characteristics of the sub-aperture, assume that the path P passes through m sub-apertures. For the kth sub-aperture, its optical characteristic index is obtained by a special optical detection method, and O k It is expressed as follows, and the weight factor based on the optical properties of the sub-aperture is defined as: For the weight factor based on the correlation with the surrounding paths, the total length L of the overlapping line segments of the statistical path P and other paths is calculated. overlap and the total length L of path P P ,Right now The weight factor based on the correlation with surrounding paths is defined as: Among them, γ is the correlation influence coefficient.

6. The wavefront restoration method based on multipath integration according to claim 1, characterized in that: In the multipath calculation phase evaluation value set step, for each shortest path found, the phase evaluation values ​​of adjacent subapertures are calculated in sequence along the weighted path based on the Southwell model. During the calculation process, the phase evaluation value of the current subaperture is calculated based on the calculation result of the previous subaperture and the slope data corresponding to the current subaperture in the order of the paths, and the calculation formula is: Where i and j represent the subaperture number, S represents the slope, P represents the phase, and h is a constant. Due to discrete sampling, the formula can be written in matrix form: S = AP, where S is a vector matrix containing all slope measurement data, with the slope in the x direction in front and the slope in the y direction in the back, P is a vector matrix containing all phase evaluation points, and A is a coefficient matrix.

7. The wavefront restoration method based on multipath integration according to claim 1, characterized in that: In the forward integration step of obtaining the phase value, the numerical integration method is used to calculate the integral result in the sub-region, and the importance weight is determined by combining the area ratio and energy contribution of the sub-region. The wavefront phase value of the point to be estimated is obtained by weighted average. Specifically, the phase evaluation value result set is N is the number of phase evaluation values, and the corresponding spatial position coordinates are (x1, y1), (x2, y2), ..., (x N ,y N ), for two adjacent phase evaluation values and Its phase change rate is Set the wavefront local change threshold ε, when When , this is determined as the boundary of the sub-region, so as to realize the division of the sub-region. Assume that the number of divided sub-regions is M, and the phase evaluation value set contained in the j-th sub-region is For each divided sub-region, the numerical integration method is used to calculate its integral value. For the j-th sub-region, its integral value I j The calculation formula is: where Δs j，k Represents the spatial interval corresponding to two adjacent phase evaluation values ​​in the jth sub-region.

8. The wavefront restoration method based on multipath integration according to claim 1, characterized in that: In the forward integration step of obtaining the phase value, the importance weight of the sub-region is determined by combining the area ratio and energy contribution of the sub-region, and the wavefront phase value of the point to be estimated is obtained by weighted average. Specifically, the total area corresponding to the entire wavefront is S total , the area corresponding to the jth sub-region is S j , then the weight factor of the sub-region area ratio is w S，j for: For the jth sub-region, the relationship function between the light field energy and the phase is: Then the energy E of the jth sub-region is j It can be expressed as: Then the sub-region energy contribution weight factor w E，j for: Combining the two weight factors of sub-region area ratio and energy contribution, the importance weight w of sub-region j is obtained j , the calculation formula is Based on the obtained weights, the wavefront phase value of the point to be estimated is obtained by weighted averaging the integral results of each sub-region. The calculation formula is:

9. The wavefront restoration method based on multipath integration according to claim 1, characterized in that: In the wave aberration information calculation step, the geometric optical model is used to estimate the approximate range and main characteristics of the wave aberration based on the phase value and the slope matrix through the ray tracing principle. Specifically, for the preliminary calculation of the wave aberration of the geometric optical model, the sub-aperture coordinate is set to (x i ,y i ), the corresponding phase value is The slope value in the x direction is In the y direction Calculate the direction vector of the light at the sub-aperture according to the slope value, then the unit direction vector component of the light in the x direction is and the unit direction vector component in the y direction It can be expressed as: For the light propagation path deviation, assume that the light starts from the reference plane and the path deviation after passing through the sub-aperture in the x direction is Δx i , in the y direction is Δy i , through the formula Calculate, where Δz j Represents the distance that light travels between adjacent sub-apertures. By combining the path deviations at all sub-apertures, we can obtain the wave aberration vector based on the geometric optics model:

10. The wavefront restoration method based on multipath integration according to claim 1, characterized in that: In the wave aberration information calculation step, the preliminary estimation result is accurately corrected and calculated based on the wave optics model at the edge of the wavefront and local details taking into account diffraction and interference phenomena. And the initial light source conditions to construct the wavefront complex amplitude distribution U geo (x, y), the calculation formula is: Where A(x, y) is the amplitude distribution, is the phase distribution estimated based on geometric optics. Substituting it into the correction formula, we get the corrected complex amplitude distribution U wave (x′, y′), that is Among them, the wave number λ is the wavelength of light, i is the imaginary unit, z represents the axial distance between the wavefront plane and the observation plane, and the corrected complex amplitude distribution U wane Extract phase information from (x′, y′) Calculate the exact wave aberration vector in the wave optics model based on the phase information This allows accurate corrections to be made to the geometrical optics estimation results at the wavefront edge and local details.

Images