Acquisition method of laser interference light spot simulation model and application thereof

Abstract

The invention relates to the technical field of optical simulation imaging, in particular to an acquisition method and application of a laser interference spot simulation model, and the method comprises the following steps: carrying out image preprocessing, edge detection, Hough circle detection and radial intensity analysis on a ray tracing result to obtain a characteristic circular ring and the position and size thereof; a related mathematical model is established according to different laser parameters (such as a laser incident angle, a laser transmission distance and laser emission power) and characteristic circular ring parameters (such as different positions, internal power density and peak power density), and finally a laser interference spot simulation model is obtained. The calculation complexity of the simulation model is remarkably reduced, the simulation process carried out by utilizing the simulation model has the characteristics of low calculation complexity, high calculation speed, strong universality and the like, and compared with a ray tracing structure, the fitting precision of a light spot image obtained by simulation is high.

Description

Technical Field

[0001] The present invention relates to the technical field of optical simulation imaging, and particularly relates to a method for obtaining a laser interference spot simulation model and its application. Background Art

[0002] In modern optoelectronic detection systems, laser interference is an important optoelectronic countermeasure method, which is widely used in fields such as military, aerospace, and remote sensing measurement. Laser interference not only reduces the imaging quality of the detection system, but may also form complex stray light through reflection, scattering, and diffraction within the optical system, thereby further weakening the detection ability of the system.

[0003] After the laser enters the optical system, its propagation path is affected by the reflection, refraction, and surface scattering of optical elements, and finally forms a non-uniform light energy distribution on the imaging surface. The shape and energy distribution of the spot are not only affected by the laser incident angle, but are also closely related to the focal length, lens curvature, material properties, and internal structure of the optical system. Therefore, by accurately simulating the laser interference effect at different incident angles and analyzing the spatial distribution characteristics of the spot, it can provide a key basis for evaluating the anti-interference ability of optoelectronic detection systems.

[0004] Currently, there are mainly two simulation methods for laser interference spots on the detector target surface: one is the simulation method based on ray tracing, which generates a spot energy distribution map by tracking the propagation paths (reflection, refraction, scattering, etc.) of a large number of discrete rays in the optical system and statistically analyzing the ray distribution on the target surface; the other method is to perform calculations based on a theoretical model, which models and simulates the laser energy distribution caused by diffraction, interference, and scattering effects respectively according to the laser propagation theory, and obtains the total laser energy distribution acting on the optoelectronic detector after superposition. However, this type of technology has the following disadvantages:

[0005] (1) Low calculation efficiency. Traditional ray tracing methods require a large number of ray simulation calculations, which are difficult to meet the real-time simulation requirements. Especially when high-precision simulation is required, the number of rays needs to reach the million level or even higher, resulting in too long calculation time;

[0006] (2) Limited fitting accuracy. The simulation calculation method based on the theoretical model has idealized condition assumptions and does not consider actual factors such as optical aberration and diffraction at the non-ideal aperture edge, and cannot reflect the diffuse distribution characteristics of the real spot;

[0007] (3) Insufficient versatility. When analyzing different simulation conditions, the simulation method based on ray tracing needs to adjust the model parameters. Ray tracing needs to adjust the model parameters for different systems or scenarios, and the reusability is limited.

[0008] In 2023, Li Haobo selected the Schlick BRDF model in the article "Research on Laser Interference Imaging and Stray Light Simulation Based on Ray Tracing" (Research on Laser Interference Imaging and Stray Light Simulation Based on Ray Tracing [D]. Xidian University, 2023). Based on the rough surface tracing algorithm process, the beam landing point and energy distribution on the image plane were simulated. Although the rough surface tracing algorithm can effectively simulate the interaction between light and the rough surface, this method still has the following disadvantages: large computational amount, low efficiency when dealing with complex surfaces, and inability to meet real-time requirements; relying on simplified mathematical models, with low accuracy.

[0009] In 2019, Wang Y, Chen Q, Lei H and others proposed in the article "The formation simulation of laser disturbing effect image" (14th National Conference on Laser Technology and Optoelectronics (LTO 2019)) that the complex interference effect of laser in the optical system is decomposed into three main physical processes, namely diffraction, interference and scattering. Theoretical models are established for diffraction and interference respectively for independent calculation, and the scattering effect is modeled and simulated through TracePro. Finally, the total energy distribution on the detector is obtained through energy superposition; however, this method has the following disadvantages: a large number of parameters need to be calibrated for the model. In practical applications, if there is a lack of prior experimental data, it may lead to a decrease in simulation accuracy. The scattering effect is not modeled, and the scattering part still relies on TracePro for ray tracing simulation, which is time-consuming in calculation. Summary of the Invention

[0010] In order to overcome the above-mentioned disadvantages of the prior art, the purpose of the present invention is to propose a method for obtaining a laser interference spot simulation model and its application. This method preprocesses the image, performs edge detection, Hough circle detection and radial intensity analysis on the ray tracing results to obtain the characteristic ring and its position and size. A mathematical model is established between different laser parameters (such as laser incident angle, laser transmission distance, laser emission power) and characteristic ring parameters (such as different positions, internal power density, peak power density), and finally a laser interference spot simulation model is obtained; the computational complexity of this simulation model is significantly reduced, and the simulation process using the simulation model has the characteristics of low computational complexity, fast calculation speed and strong versatility. Compared with the ray tracing structure, the fitting accuracy of the spot image obtained by simulation is high.

[0011] To achieve the above object, the technical solution adopted by the present invention is as follows:

[0012] In the first aspect, a method for obtaining a laser interference spot simulation model includes the following steps:

[0013] S1: Simulate the ray tracing process of the laser interference source in the optical system, obtain the numerical values of the laser interference spot power density distribution matrix for different laser incident angles, different laser transmission distances, and different laser emission powers, and convert the numerical values of the laser interference spot power density distribution matrix into an EXR image;

[0014] S2: Perform image preprocessing, edge detection, Hough circle detection, and radial intensity analysis on the EXR image obtained in step S1 in sequence to obtain characteristic rings, and extract the positions and sizes of the characteristic rings. Among them, the threshold calculation of the full width at half maximum method (FWHM) in the radial intensity analysis is as follows:

[0015] T k = 0.5 × I(r peak , θ k )

[0016] where I(r peak , θ k ) is the power density value at the point with a radial distance r k on the ray with an angle θ peak , that is, the peak power density value on the ray with an angle θ k ;

[0017] S3: Using the positions and sizes of the characteristic rings described in step S2, respectively: extract the different positions corresponding to each characteristic ring at different laser incident angles, and establish a mathematical model of the change of the different positions of the characteristic rings with the laser incident angle; extract the average power density of each radial distance inside each characteristic ring and perform Gaussian fitting, and establish a mathematical model of the change of the power density inside the characteristic rings with the radial distance; and extract the internal peak power density of each characteristic ring under different working conditions of laser incident angles, different laser transmission distances, and different laser emission powers, and establish mathematical models of the change of the peak power density of the characteristic rings with the laser incident angle, laser transmission distance, and laser emission power, respectively.

[0018] Further, the image preprocessing in step S2 is Gaussian blur processing and normalization processing in sequence. The Gaussian kernel weight in the Gaussian blur processing is a two-dimensional Gaussian distribution, and the specific formula is as follows:

[0019]

[0020] where: G(x, y) is the Gaussian weight at the coordinate point (x, y); σ is the standard deviation of the Gaussian distribution, which determines the degree of diffusion of the distribution;

[0021] The normalization processing is to map the pixel values of the image after Gaussian blur processing to the standard range (0 - 255) through the following formula:

[0022]

[0023] Wherein: I N (x, y) is the pixel value of the normalized image at the coordinate point (x, y), and I G (x, y) is the pixel value of the blurred image at the coordinate point (x, y); min(I G ) is the minimum pixel value of the blurred image, and max(I G ) is the maximum pixel value of the blurred image.

[0024] Furthermore, the Hough circle detection described in step S2 uses a parameter space voting mechanism to transform the geometric shape detection problem in the image space into an extreme value search problem in the parameter space.

[0025] Furthermore, the radial intensity analysis described in step S2 includes the following steps:

[0026] S21: Generation of radial sampling lines:

[0027] Taking the center (x c , y c ) of the Hough circle detection as the center, uniformly generate 4 or more radiation lines within the range of 0 ≤ θ ≤ 2π, and the maximum sampling length of each radiation line along the radius direction is 2 times the radius of the Hough circle detection;

[0028] S22: Extraction of power density distribution curve:

[0029] For each radiation line k described in step S21, the position of each sampling point along the direction θ k is:

[0030] x i = x c + i·Δr·cosθ k , y i = y c + i·Δr·sinθ k

[0031] Wherein, i.e., the sampling point spacing along each radiation line, i ≥ 10;

[0032] Use bilinear interpolation to obtain the power density value I(r i , θ k ):

[0033]

[0034] Wherein, I(r i , θk) is the power density value at the polar coordinates (r i , θk) obtained by interpolation; (xi, yi) are the pixel coordinates of the interpolation target point; is the integer coordinate adjacent to the upper left corner of the target point; ω mn is the interpolation weight, calculated based on the relative position of the target point among 4 adjacent pixels; m and n range from 0 to 1, representing the index offsets in the 2×2 neighborhood; are the original power density values of 4 adjacent pixels in the 2×2 grid;

[0035] S23: Peak detection and full width at half maximum analysis:

[0036] For the power density value distribution I(r,θ k ) of each radiation line described in step S22, first find the global peak position, i.e., the radial distance r peak :

[0037]

[0038] Secondly, calculate the threshold T of the full width at half maximum method (FWHM): k :

[0039] T k = 0.5×I(r peak ,θ k )

[0040] where I(r peak ,θ k ) is the power density value of the point with radial distance r on the radiation line at angle θ, i.e., the peak power density value on the radiation line at angle θ k ; peak k k ;

[0041] S24: Inner radius detection, search from the center of the global peak position described in step S23 towards the peak direction to find the first position r k that satisfies I(r,θ k )≥T in,k :

[0042] r in,k = max{r|I(r,θ k )≥T k and r<r peak}

[0043] Outer radius detection, search from the peak outwards to find the first position r k that satisfies I(r,θ k )≤T out,k :

[0044] r out,k = min{r|I(r,θ k )≤Tk and r>r peak}

[0045] S25: Based on the inner radius and outer radius calculated in step S24, calculate the median of the inner radius and outer radius for all angles to obtain the stable inner radius and outer radius.

[0046] Furthermore, the characteristic circular ring described in step S2 is obtained through the following process:

[0047] Generate an annular mask image based on the stable inner radius and outer radius described in step S25. The annular mask image is a binary image, with the annular region being 1 and the rest being 0:

[0048]

[0049] where Mask(x, y) is the pixel value of the annular mask image, indicating whether the point (x, y) belongs to the annular region (1 means inside the ring, 0 means outside the ring); (x c , y c ) is the center coordinate of the characteristic circular ring; rin is the inner radius of the characteristic circular ring; rout is the outer radius of the characteristic circular ring.

[0050] Furthermore, the mathematical model of the change of different positions of the characteristic circular ring with the laser incident angle in step S3 is:

[0051] x = aθ + b

[0052] where x is the center position of the characteristic circular ring; θ is the incident angle; a is the proportionality factor, indicating the change rate of the center position x when θ changes; b is the offset, indicating the initial value of the center position x when the incident angle θ = 0;

[0053] Use the least squares method to solve the sum of the squares of the minimum errors S, thereby determining a and b of the fitting curve equation for the position transformation of each characteristic circular ring:

[0054]

[0055] where S is the sum of the squares of the minimum errors; n represents the number of data points, that is, the number of different incident angles simulated in step S2; x i represents the center position of the current characteristic circular ring at the i-th incident angle; θ i is the incident angle value of the i-th incident angle.

[0056] Furthermore, the average power density described in step S3 satisfies the following formula:

[0057]

[0058] where Iavg (r) is the average value of the power density corresponding to the pixels in all sector areas on the radius r; N is the number of pixel points in all sector areas on the radius r; I char_ring (x i , y i ) is the power density value at the characteristic circular ring independent image coordinate point (x i , y i ), and x i , y i satisfies:

[0059] In step S3, the Gaussian fitting mentioned uses a Gaussian function to fit the radial mean change curve I avg_fit (r) of each selected characteristic circular ring to obtain a mathematical model of the power density change with the radial distance inside the characteristic circular ring:

[0060]

[0061] where I avg_fit (r) is the power density value of the pixel points on the fitted radius r, a is the maximum amplitude, r is the radius, b is the center position of the Gaussian distribution, c is the standard deviation, and the least squares method is used for nonlinear fitting to solve the optimal parameters a, b, c.

[0062] Furthermore, the calculation process of the peak power density in step S3 is as follows:

[0063] S31. Calculate the radial mean: For the simulation results at different angles of the same characteristic circular ring, the starting and ending angles of the statistical sector area are set respectively, and the parts overlapping with other circular rings are removed. For each radius r, r in ≤ r ≤ r out , and the mean value of the pixel values in all sector areas on this radius is statistically calculated:

[0064]

[0065] To reduce the influence of noise, the extreme value removal mean method (such as removing the upper and lower 10% of the pixel values) is adopted:

[0066] I trimmed (r) = trimmean(I avg (r), 10%)

[0067] where I trimmed (r) is the average value of the remaining data calculated after removing the highest and lowest 10% of the data in the power density data set corresponding to the pixels in all sector areas on the radius r;

[0068] S32. Extract the peak power density: For each incident angle θ i, take the peak power density from the radial mean:

[0069] I max (θ i ) = max{I trimmed_i (r)}

[0070] where I max (θ i ) is the peak power density of the characteristic ring when the incident angle is θ i .

[0071] Furthermore, the mathematical model of the peak power density of the characteristic ring described in step S3 changing with the laser incident angle is:

[0072]

[0073] where I max (θ) represents the peak power density of the characteristic ring when the incident angle is θ; A is the power density when incident at the reference incident angle θ0; B is the attenuation coefficient; θ0 is the reference incident angle;

[0074] The mathematical model of the peak power density of the characteristic ring described in step S3 changing with the laser transmission distance is:

[0075] I(d) = I0·d -α

[0076] where I(d) is the power density at the transmission distance d, I0 is the power density at the reference distance of 1 km, d is the transmission distance, and α is the attenuation exponent;

[0077] The mathematical model of the peak power density of the characteristic ring described in step S3 changing with the laser emission power is:

[0078] I(P) = k·P

[0079] where I(P) represents the power density when the laser emission power is P; k is the proportionality coefficient, and P is the laser emission power.

[0080] The mathematical models obtained in step S3 for the different positions changing with the laser incident angle, the power density inside the characteristic ring changing with the radial distance, and the peak power density of the characteristic ring changing with the laser incident angle, laser transmission distance, and laser emission power are collectively called the laser interference spot simulation model.

[0081] In the second aspect, an application of a laser interference spot simulation model, based on the laser interference spot simulation model obtained by the acquisition method.

[0082] Compared with the prior art, the beneficial effects of the present invention are:

[0083] 1. In step S3 of this method, by establishing a non-linear model of the positions of each ring varying with the incident angle and a model of the power density distribution in the radial distance within the ring, the aberration effect is fully retained, and the fitting accuracy of the real optical system effect is high. Comparing the ray tracing results of TracePro software with the model simulation results, the average error is less than 10%.

[0084] 2. This method extracts the characteristic rings of the laser interference spot and performs hierarchical parametric modeling, converting the ray tracing of millions of rays into the mathematical description of a small number of characteristic rings. Compared with the traditional ray-by-ray calculation, only the position, radial power density, and parameter relationship model of the finite characteristic rings need to be solved, and the computational complexity is significantly reduced. Using TracePro for ray tracing, a single simulation takes more than half an hour, while the simulation calculation of the laser interference spot simulation model obtained by this method only takes a few seconds, and the simulation calculation efficiency is significantly improved.

[0085] 3. This method decouples the influence of variables such as the incident angle, transmission distance, and emission power of the laser interference light source on the spot through a hierarchical modeling framework. Through physically driven parameter separation modeling, adjusting model parameters is avoided, making the model universal.

[0086] In summary, the present invention determines the stable inner and outer radii of the characteristic rings through the power density in the radial distance, and further establishes a mathematical model by layering according to the incident angle, transmission distance, and emission power of the laser interference light source. Finally, the computational complexity of the obtained laser interference spot simulation model is significantly reduced. The simulation process using the simulation model is fast, highly versatile, and the fitting accuracy of the simulated spot image is high. Description of the Drawings

[0087] Figure 1 is the derivation flow chart of the laser interference spot simulation model in this method.

[0088] Figure 2 is the simulation flow chart using the simulation model obtained by this method.

[0089] Figure 3 is the Cassegrain optical system with a sunshade in the embodiment of this method.

[0090] Figure 4(a) is the original data image of the characteristic ring extraction result in step S2 of this method.

[0091] Figure 4(b) is the central bright spot of the characteristic ring extraction result in step S2 of this method.

[0092] Figure 4(c) is the leftmost ring of the characteristic ring extraction result in step S2 of this method.

[0093] Figure 4(d) is the second left ring of the characteristic ring extraction result in step S2 of this method.

[0094] Figure 4(e) is the rightmost ring in the extraction result of the feature rings in step S2 of this method.

[0095] Figure 4(f) is the second right ring in the extraction result of the feature rings in step S2 of this method.

[0096] Figure 5(a) is the original spot image when selecting a fan-shaped area for data statistics in step S3 of this method.

[0097] Figure 5(b) is the acquisition of the fan-shaped area when selecting a fan-shaped area for data statistics in step S3 of this method.

[0098] Figure 6(a) is the simulated laser ray tracing spot image in the first embodiment of this method.

[0099] Figure 6(b) is the simulated spot image of the laser interference spot simulation model in the first embodiment of this method.

[0100] Figure 7(a) is the simulated laser ray tracing spot image in the second embodiment of this method.

[0101] Figure 7(b) is the simulated spot image of the laser interference spot simulation model in the second embodiment of this method.

[0102] Figure 8 It is a comparison chart of the fitting results between the ray tracing spot image and the simulated spot image in the second embodiment of this method.

[0103] Figure 9(a) is the simulated laser ray tracing spot image in the third embodiment of this method.

[0104] Figure 9(b) is the simulated spot image of the laser interference spot simulation model in the third embodiment of this method.

[0105] Figure 10 It is a comparison chart of the fitting results between the ray tracing spot image and the simulated spot image in the third embodiment of this method. Specific implementation manners

[0106] The following Figures 1 to 10 is a further detailed description of the present invention:

[0107] Figure 1 Show the derivation flow chart of the laser interference spot simulation model in this method.

[0108] Selection of the laser interference source and the optical system in step S1:

[0109] (1) Modeling of the laser interference source: In the present invention, a circular lattice light source with a certain divergence angle and Gaussian distribution is used to simulate the Gaussian laser interference source. The intensity distribution of the Gaussian beam can be expressed by the following formula:

[0110]

[0111] Wherein: I laser (r) is the light intensity at a radial distance r from the optical axis, I0 is the maximum light intensity at the center, r is the radial distance from the optical axis, and ω is the beam radius (i.e., the 1 / e 2 width).

[0112] By setting the position and orientation of the laser interference source, the laser incident angle and transmission distance can be adjusted; by setting the total laser luminous flux, the emission power can be adjusted; the number of lattice graphic rings is set according to the simulation accuracy requirements. The higher the number of rings, the more the number of traced rays and the higher the simulation accuracy.

[0113] (2) Optical system modeling: The structure of the optical system is as Figure 3 shown. In the present invention, the simulation optical system adopts the Cassegrain system structure with a correction lens group, and the primary and secondary mirrors and the outer light shield are designed to suppress the main stray light.

[0114] Obtaining the ray tracing result and the EXR image in step S1:

[0115] Use TracePro software to simulate emitting a large number of rays from the laser interference source described in step S1 into the conical region. By tracing the position where each ray hits the detector target surface after passing through the optical system transmission, and finally by counting the number of rays hitting each pixel point on the detector target surface, the power density value of each pixel point is calculated, and then the laser interference spot image on the entire detector target surface is obtained. In the present invention, the control variable method is used to set different laser incident angles, laser transmission distances, and laser emission powers for ray tracing, obtain the laser interference spot images on the detector target surface under different laser parameters, use the built-in function of TracePro software to save the laser interference spot power density distribution matrix values in TXT format, and finally by reading the saved TXT results, write the laser interference spot power density distribution matrix values into an EXR single-channel image, namely the EXR image, for convenient subsequent processing.

[0116] The process of extracting the characteristic rings in step S2 is as follows:

[0117] (1) Image preprocessing. Since the dynamic range of the EXR image is large, directly performing edge detection may lead to unstable results. Therefore, the following preprocessing is required:

[0118] a. Gaussian blur processing: Suppress high-frequency noise through convolution with a Gaussian kernel to improve the overall smoothness, which is a pre-step of the Canny edge detection algorithm; the Gaussian kernel weights are a two-dimensional Gaussian distribution as follows:

[0119]

[0120] Where: G(x, y) is the Gaussian weight at the coordinate point (x, y); σ is the standard deviation of the Gaussian distribution, which determines the degree of diffusion of the distribution.

[0121] Perform a discrete convolution operation on the original EXR image and the Gaussian kernel to obtain the blurred image:

[0122]

[0123] Where: I G (x, y) is the pixel value of the blurred image at the coordinate point (x, y); k is the Gaussian kernel radius; G(i, j) is the Gaussian weight, representing the contribution of the pixel point at the offset (i, j) to the central pixel; I R (x + i, y + j) is the pixel value of the original image at the coordinate point (x + i, y + j).

[0124] b. Normalization processing: Map the pixel values of the image after Gaussian blurring to the standard range (0 - 255) through the following formula for subsequent processing:

[0125]

[0126] Where: I N (x, y) is the pixel value of the normalized image at the coordinate point (x, y), I G (x, y) is the pixel value of the blurred image at the coordinate point (x, y); min(I G ) is the minimum pixel value of the blurred image, max(I G ) is the maximum pixel value of the blurred image.

[0127] (2) Edge detection. The Canny algorithm extracts continuous and refined edges through multiple steps, solving the problems of wide edges and false edges existing in other edge detection algorithms, and is a prerequisite for the Hough transform.

[0128] The intensity values after the previous Gaussian smoothing need to be processed as follows: a. Gradient calculation: The commonly used Sobel operator is an edge detection method based on convolution. It uses two 3×3 convolution kernels to calculate the gradients in the x - direction and y - direction respectively:

[0129] Horizontal gradient kernel G x :

[0130]

[0131] Vertical gradient kernel The gradient magnitude is calculated as:

[0132]

[0133] The gradient direction is:

[0134]

[0135] b. Non-maximum suppression, refine the edge to a single-pixel width, and eliminate the thick-edge response. For each pixel, compare the adjacent pixels along the gradient direction. If the gradient amplitude of the current pixel is not the local maximum in this direction, set it to zero.

[0136] c. Double-threshold detection and edge connection, using hysteresis thresholding to suppress noise while preserving the continuity of real edges. First, use the high threshold T H to select strong-edge pixel points. These pixel points have a relatively high gradient amplitude and can almost be determined to belong to real edges. Then, use the low threshold T L to select weak-edge pixel points. These pixel points may be part of the edge or may be noise. To remove the noise and retain the real edges, through edge connection, only those weak-edge pixels that can be connected to strong-edge pixels are retained, while isolated weak-edge pixels are suppressed.

[0137] (3) Hough circle detection, transform the geometric shape detection problem in the image space into an extreme value search problem in the parameter space through the parameter space voting mechanism. The specific process is as follows:

[0138] a. Candidate center voting for the circle. For each edge point (x i , y i ), vote accumulatively in the parameter space (xc, yc) along the reverse direction of its gradient direction (i.e., the center direction of the circle). The possible center coordinates satisfy:

[0139] x c = x i - r·cos(θ)

[0140] y c = y i - r·sin(θ)

[0141] where r is the unknown radius.

[0142] All edge points vote for the potential center positions along the reverse direction of the gradient, forming an accumulator matrix. The peak value (exceeding the threshold) of the accumulator is the candidate center (x c , y c ).

[0143] b. Radius determination. For each candidate center (x c , y c ), collect the corresponding edge points and calculate the distance ri from each point to the center:

[0144]

[0145] Count all r i histogram, and the radius corresponding to the peak value is the final result. Through the above steps, the Hough circle detection algorithm can efficiently locate the center and radius in the parameter space, that is, determine the center and radius of the laser interference spot image.

[0146] (4) Radial intensity analysis: Determine the center and preliminary radius through Hough circle detection, then sample the power density distribution along 36 radiation lines, and use bilinear interpolation to obtain data. Through peak detection and full width at half maximum method, search from the center outwards to determine the inner radius respectively, and search from the peak outwards to determine the outer radius. Finally, take the median of all angles calculated to ensure the radius stability, so as to locate the physical boundary of the annular structure.

[0147] The accuracy and robustness of spot feature extraction are improved through the ring analysis method: 1. More accurate energy distribution evaluation: The ring can accurately depict the inner and outer boundaries of the spot, while a single circle usually only provides an approximate center position and cannot reflect the actual expansion range of the spot; 2. More stable feature extraction: Using peak detection and full width at half maximum method to calculate the ring boundary can reduce the influence of noise, while the traditional circle fitting method is vulnerable to local outliers.

[0148] The radial intensity analysis includes the following steps:

[0149] a. Generation of radial sampling lines:

[0150] With the center (xc, yc) detected by the Hough circle as the center, uniformly generate 4 or more radiation lines in the range of 0 ≤ θ ≤ 2π, and the maximum sampling length of each radiation line along the radius direction is 2 times the radius detected by the Hough circle; in this embodiment, 36 radiation lines are uniformly generated in the range of 0 ≤ θ ≤ 2π;

[0151] b. Extraction of power density distribution curve:

[0152] For each radiation line k, along the direction θ k the position of each sampling point is:

[0153] x i = x c + i·Δr·cosθ k , y i = y c + i·Δr·sinθ k

[0154] where i.e., the sampling point spacing along each radiation line, i ≥ 10; in this embodiment, i takes 100.

[0155] Use bilinear interpolation to obtain the intensity value I(ri , θ k ):

[0156]

[0157] wherein, I(r i , θk) is the power density value at the polar coordinates (r i , θk) obtained by interpolation; (xi, yi) are the pixel coordinates of the interpolation target point; are the integer coordinates adjacent to the upper left corner of the target point; ω mn is the interpolation weight, calculated based on the relative position of the target point among 4 adjacent pixels; m, n traverse 0 and 1, representing the index offsets in the 2×2 neighborhood; are the original power density values of 4 adjacent pixels in the 2×2 grid. The construction of the polar coordinate power density distribution I(r, θ) of the annular region is completed, thereby reflecting the variation of the power density with the radius and the angle.

[0158] c. Peak detection and full width at half maximum analysis:

[0159] For the power density distribution I(r, θ k ) of each radiation line, first find the global peak position:

[0160]

[0161] Calculate the threshold of the full width at half maximum (FWHM) method:

[0162] T k = 0.5 × I(r peak , θ k )

[0163] wherein, I(r peak , θ k ) is the power density value at the point with the radial distance r k on the radiation line at the angle θ peak , that is, the peak power density value on the radiation line at the angle θ k ;

[0164] Inner radius detection, search from the center towards the peak direction to find the first position r k that satisfies I(r, θ k ) ≥ T in,k :

[0165] r in,k = max{r | I(r, θ k ) ≥ T k and r < r peak}

[0166] Outer radius detection: Search outward from the peak to find the first position \(r\) that satisfies \(I(r,\theta)\leq T\). k )\leq T k : out,k :

[0167] r out,k = \min\{r|I(r,\theta)\leq T k )\leq T k and \(r > r peak \}\)

[0168] Finally, by calculating the median of the inner and outer radii for all angles, more stable inner and outer radii are obtained.

[0169] (5) Feature ring extraction.

[0170] First, generate a circular mask image, which is a binary image with the circular area being 1 and the rest being 0:

[0171]

[0172] where \(Mask(x,y)\) is the pixel value of the mask image, indicating whether the point \((x,y)\) belongs to the circular area (1 means inside the ring, 0 means outside the ring); \((x c ,y c ) is the center coordinate of the feature ring; \(r in \) is the inner radius of the feature ring; \(r out \) is the outer radius of the feature ring.

[0173] By multiplying the spot image and each mask image pixel by pixel, each feature ring is extracted from the original spot image under each simulation condition, that is, different laser emission powers, laser incident angles, and laser transmission distances, to obtain independent images of each feature ring for subsequent processing.

[0174] As Figures 4(a) to 4(f) shown, the extraction results of some feature rings on the detector target surface are given when the laser transmission distance is 1 km, the emission power is 5 W, and the incident angle is 0.06°. From the extraction results, through the implementation method in step 2, the complete areas of each feature ring can be extracted from the original image, facilitating further analysis later.

[0175] In step S3, establish a mathematical model for the variation of different positions of each feature ring with the laser incident angle.

[0176] After extracting each characteristic ring (i.e., the independent image of each characteristic ring) from the interference light spots with different incident angles simulated by TracePro software, the same characteristic ring is analyzed. In the simulation, for the convenience of fitting, all incident directions and the optical axis are in the same plane, and since the fitting angle is small, the center position of the ring is approximately linearly moved horizontally along the center of the image.

[0177] The center position x of each characteristic ring and the incident angle θ satisfy a linear relationship:

[0178] x = aθ + b

[0179] where x is the center position of the characteristic ring; θ is the incident angle; a is the proportionality factor, representing the change rate of the center position x when θ changes; b is the offset, representing the initial value of the center position x when the incident angle θ = 0.

[0180] The least squares method is used to solve the sum of the squares of the minimum errors S, so as to determine a and b of the fitting curve equation for the position transformation of each characteristic ring:

[0181]

[0182] where S is the sum of the squares of the minimum errors; n represents the number of data points, that is, the number of different incident angles simulated in step 2; x i represents the center position of the current characteristic ring at the i-th incident angle; θ i is the numerical value of the incident angle of the i-th incident angle.

[0183] In step S3, a mathematical model for the change of the internal power density of each characteristic ring with the radial distance is established.

[0184] (1) Calculate the radial mean value of each characteristic ring:

[0185] For each characteristic ring, select the independent image of the characteristic ring obtained in step S3 with less overlap with other rings, set the start and end angles of the statistical sector area, and further eliminate the influence of the overlapping part.

[0186] As shown in Figure 5(a), there is an overlapping area with other characteristic rings inside the right area of the characteristic ring image extracted by step S3, resulting in this area being significantly brighter. In order not to be affected by this area, by setting the start angle to 90° and the end angle to 270°, the area wrapped by the green border shown in Figure 5(b) is specified as the statistical sector area. For each radius r, r in ≤ r ≤ r out , count the pixel values in all sector areas on this radius, and calculate the mean value of the power density:

[0187]

[0188] Among them, I avg (r) is the average power density corresponding to the pixels in all sector areas on the radius r; N is the number of pixel points in all sector areas on the radius r; I char_ring (x i , y i ) is the power density value at the independent image coordinate point (x i , y i ) of the characteristic ring, and x i , y i satisfies:

[0189] (2) Gaussian fitting:

[0190] Use the Gaussian function to fit the change curve I avg_fit (r) of the radial mean value of each selected ring:

[0191]

[0192] Among them, I avg_fit (r) is the power density value of the pixel point on the fitted radius r, a is the maximum amplitude, r is the radius, b is the center position of the Gaussian distribution, and c is the standard deviation. Use the least squares method for nonlinear fitting to solve the optimal parameters a, b, c.

[0193] In step S3, establish a mathematical model of the power density of each characteristic ring changing with different parameters.

[0194] (1) The calculation process of the peak power density of each characteristic ring is as follows:

[0195] a. Calculate the radial mean value: For the simulation results at different angles of the same characteristic ring, respectively set the starting and ending angles of the statistical sector area, and remove the part overlapping with other rings (the same as step S3). For each radius r, r in ≤ r ≤ r out , and statistically calculate the mean value of the pixel values in all sector areas on this radius:

[0196]

[0197] To reduce the influence of noise, use the method of removing extreme values and taking the mean (such as removing the upper and lower 10% of the pixel values):

[0198] I trimmed (r) = trimmean(I avg (r), 10%)

[0199] Among them, I trimmed(r) is the average value of the remaining data calculated after removing the highest and lowest 10% of the data in the power density data set corresponding to the pixels in all fan-shaped regions with radius r.

[0200] b. Extract the peak power density: For each incident angle θ i , take the peak power density from the radial mean value:

[0201] I max (θ i ) = max{I trimmed_i (r)}

[0202] where I max (θ i ) is the peak power density of the characteristic ring when the incident angle is θ i .

[0203] (2) Fit the power density of each characteristic ring with the change of the laser incident angle.

[0204] Based on the original data determined by steps a and b in step S3(1), perform a fit of the peak power density with the change of the incident angle. The peak power density changes exponentially with the incident angle. The mathematical model of the peak power density of the characteristic ring with the change of the laser incident angle is:

[0205]

[0206] where I max (θ) represents the peak power density of the characteristic ring at the incident angle θ; A is the power density when incident at the reference incident angle θ0; B is the attenuation coefficient; θ0 is the reference incident angle.

[0207] (3) Fit the power density of each characteristic ring with the change of the laser transmission distance.

[0208] For the independent images of each characteristic ring at different transmission distances extracted in step S2, based on the original data determined by steps a and b in step S3(1), the power density changes according to a power function distribution. The mathematical model of the peak power density of the characteristic ring with the change of the laser transmission distance is:

[0209] I(d) = I0·d -α

[0210] where I(d) is the power density at the transmission distance d, I0 is the power density at the reference distance of 1 km, d is the transmission distance, and α is the attenuation exponent.

[0211] (4) Fit the power density of each characteristic ring with the change of the laser emission power.

[0212] For the independent images of each characteristic ring under different laser emission powers extracted in step S2, based on the original data determined by steps a and b in step S3(1), the power density has a linear relationship with the change in laser emission power. The mathematical model of the peak power density of the characteristic ring varying with the laser emission power is as follows:

[0213] I(P) = k·P

[0214] where I(P) represents the power density when the laser emission power is P; k is the proportionality coefficient, and P is the laser emission power.

[0215] The mathematical models obtained in step S3 for the variation of different positions with the laser incident angle, the variation of the power density inside the characteristic ring with the radial distance, and the variation of the peak power density of the characteristic ring with the laser incident angle, laser transmission distance, and laser emission power are collectively referred to as the laser interference spot simulation model.

[0216] In steps 2 to 3 of this method, physical reference data is generated based on TracePro high-precision ray tracing. The characteristic rings of the spot are extracted through the inversion extrapolation technique, and the laser interference spot simulation model, that is, the inversion extrapolation model, is constructed. This not only retains the real optical effects but also avoids the empirical deviation of pure numerical methods, forming a closed-loop link of "physical simulation - feature inversion - model extrapolation".

[0217] In step 3 of this method, the spot is decomposed into multiple characteristic rings, and the positions, radial power density distributions, and parameter responses of each ring are independently modeled.

[0218] Figure 2 The fitting process of the simulation using the simulation model obtained by this method is shown: The laser interference spot simulation model is established by this method. When the laser emission power, transmission distance, and incident angle are given, substituting them into the corresponding mathematical model to determine the positions of each characteristic ring and the power density values at each radial distance; by superimposing the power density distributions of each characteristic ring, the overall fitted spot can be obtained. The specific simulation process is as follows:

[0219] (1) Fitting of the reference radial power density of each characteristic ring. Through the internal radial power density attenuation model of each characteristic ring, the reference radial power distribution data of each characteristic ring can be directly generated.

[0220] (2) Calculation of the peak power density of each characteristic ring. Substitute the incident angle, transmission distance, and emission power into the calculation model of the peak power density change of the characteristic ring respectively, calculate the product of each influencing factor, determine the gain of the peak power density value of each characteristic ring, and multiply the reference radial power distribution data by the gain to obtain the radial power density distribution of the characteristic ring.

[0221] (3) Determine the positions of each characteristic ring. Given the incident angle, the center positions of the characteristic rings can be calculated respectively through the model x = aθ + b of the change of each characteristic ring position with the incident angle, and the whole of each characteristic ring is moved.

[0222] (4) Generate the simulated spot image. Linearly superimpose the spot-fitted images of each characteristic ring to generate the final spot.

[0223] Compare the simulated spot image with the ray tracing result to verify the accuracy.

[0224] Establish a general laser interference spot simulation model through the inversion of finite TracePro samples, realize the rapid prediction of the spot distribution across systems, reduce the calculation time from several hours to several seconds, and ensure that the accuracy error is less than 10%.

[0225] Example 1

[0226] Figure 6(a) shows the simulated laser ray tracing spot image obtained by using TracePro software when the simulation parameters are laser emission power of 2W, transmission distance of 5km, and incident angle of 0.07°; Figure 6(b) shows the spot image obtained by simulation calculation through the laser interference spot simulation model in this method when the simulation parameters are laser emission power of 2W, transmission distance of 5km, and incident angle of 0.07°.

[0227] Example 2

[0228] Figure 7(a) and Figure 7(b) respectively give the simulated laser ray tracing spot image obtained by using TracePro software and the simulated calculation spot image obtained by simulation calculation through the laser interference spot simulation model in this method when the laser emission power is 5W, the transmission distance is 1km, and the incident angle is 0°;

[0229] Example 3

[0230] Figure 9(a) and Figure 9(b) respectively give the simulated laser ray tracing spot image obtained by using TracePro software and the simulated calculation spot image obtained by simulation calculation through the laser interference spot simulation model in this method when the laser emission power is 5W, the transmission distance is 1km, and the incident angle is 0.03°.

[0231] From Figure 7(a), Figure 7(b), Figure 9(a), and Figure 9(b), it can be seen respectively that the distributions of the corresponding two in each characteristic ring and the gray scale display of the rings are approaching consistency; Figure 8 、 Figure 10Furthermore, the corresponding numerical comparison curve of the light spot image power density between the light ray tracing result of TracePro software and the simulation calculation along the centrosymmetric direction is given. According to the drawn curves, it can be seen that except for a small fluctuation caused by the influence of the number of light ray tracings in TracePro software at individual positions, the two curves basically coincide and have the same trend. Through the above judgment, it shows that the laser interference light spot simulation model (i.e., the inversion extrapolation model) constructed by this method has reliability and authenticity.

[0232] The working principle of the present invention is as follows:

[0233] Based on the light ray tracing result of the simulated laser interference source in the optical system, this method obtains the characteristic ring and its position and size by performing image preprocessing, edge detection, Hough circle detection, and radial intensity analysis on the light ray tracing result, establishes a mathematical model between different laser parameters (such as laser incident angle, laser transmission distance, laser emission power) and characteristic ring parameters (such as different positions, internal power density, peak power density), and finally obtains the laser interference light spot simulation model; through simulation experiments, the simulation calculation results of this simulation model are consistent with the light ray tracing results in the distribution of each characteristic ring and the grayscale display of the ring, indicating that the physical authenticity and reliability of this simulation model are relatively strong.

Claims

1. A method for obtaining a simulation model of laser interference spot, characterized in that It includes the following steps: S1: Simulate the ray tracing process of the laser interference source in the optical system, obtain the numerical values of the laser interference spot power density distribution matrix for different laser incident angles, different laser transmission distances, and different laser emission powers, and convert the numerical values of the laser interference spot power density distribution matrix into an EXR image; S2: Perform image preprocessing, edge detection, Hough circle detection, and radial intensity analysis on the EXR image obtained in step S1 in sequence to obtain characteristic rings, and extract the positions and sizes of the characteristic rings. Among them, the threshold calculation of the full width at half maximum method (FWHM) in the radial intensity analysis is as follows: T k = 0.5 × I(r peak , θ k ) Among them, I(r peak ,θ k ) is the angle θ k Radial distance r peak The power density value of the point, that is, the angle θ k The peak power density value of the radiation line; S3: Using the positions and sizes of the characteristic rings in step S2, respectively: extract the different positions corresponding to each characteristic ring under different laser incident angles, and establish a mathematical model of the change of the different positions of the characteristic rings with the laser incident angle; extract the average power density of each radial distance inside each characteristic ring and perform Gaussian fitting to establish a mathematical model of the change of the power density inside the characteristic ring with the radial distance; and extract the internal peak power density of each characteristic ring under different laser incident angles, different laser transmission distances, and different laser emission power conditions, and establish mathematical models of the change of the peak power density of the characteristic rings with the laser incident angle, laser transmission distance, and laser emission power respectively.

2. The acquisition method according to claim 1, wherein The image preprocessing in step S2 is Gaussian blur processing and normalization processing in sequence. The Gaussian kernel weight in the Gaussian blur processing is a two-dimensional Gaussian distribution, and the specific formula is as follows: Where: G(x, y) is the Gaussian weight at the coordinate point (x, y); σ is the standard deviation of the Gaussian distribution, which determines the degree of diffusion of the distribution; The normalization processing is to map the pixel values of the image after Gaussian blur processing to the standard range (0 - 255) through the following formula: Where: I N I(x,y) is the pixel value of the normalized image at the coordinate point (x,y), and I G I'(x,y) is the pixel value of the blurred image at the coordinate point (x,y); min(I G ) is the minimum pixel value of the blurred image, and max(I G ) is the maximum pixel value of the blurred image.

3. The acquisition method according to claim 1, wherein The Hough circle detection in step S2 uses a parameter space voting mechanism to transform the geometric shape detection problem in the image space into an extreme value search problem in the parameter space.

4. The acquisition method according to claim 1, wherein The radial intensity analysis in step S2 includes the following steps: S21: Generation of radial sampling lines: With the center of the circle detected by the Hough circle detection (x c , y c ) as the center, generate 4 or more radiation lines evenly within the range of 0 ≤ θ ≤ 2π, and the maximum sampling length of each radiation line along the radial direction is 2 times the radius detected by the Hough circle detection; S22: Extraction of power density distribution curves: For each ray k described in step S21, along the direction θ k the position of each sampling point is: x i = x c + i·Δr·cosθ k , y i = y c + i·Δr·sinθ k Among them, that is, the sampling point spacing along each radiation line, i≥10; Obtain the power density value I(r i , θ k ) at non-integer coordinates using bilinear interpolation: where I(r i , θ k ) is the power density value at the polar coordinates (r i , θ k ) obtained by interpolation; (xi, yi) are the pixel coordinates of the interpolation target point; are the integer coordinates adjacent to the upper left corner of the target point; ω mn is the interpolation weight, calculated based on the relative position of the target point among 4 neighboring pixels; m, n traverse 0 and 1, representing the index offsets in the 2×2 neighborhood; are the original power density values of 4 adjacent pixels in the 2×2 grid; S23: Peak detection and full width at half maximum analysis: For the power density numerical distribution I(r, θ k ) of each radiation ray described in step S22, first find the global peak position, that is, the radial distance r peak : Next, calculate the threshold T of the Full Width at Half Maximum (FWHM) k : T k = 0.5 × I(r peak , θ k ) where I(r peak , θ k ) is the power density value at the point with a radial distance r k from the ray in the direction of the angle θ peak , i.e., the peak power density value on the ray at the angle θ k ; S24: Inner radius detection. Search from the center of the global peak position described in step S23 in the peak direction to find the first position r k ) ≥ T k that satisfies the condition in,k : r in,k = max{r | I(r, θ k ) ≥ T k and r < r peak} Outer radius detection. Search outward from the peak to find the first position r that satisfies I(r,θ k ) ≤ T k : out,k ​ r out,k = min{r | I(r, θ k ) ≤ T k and r > r peak} S25: Based on the inner radius and outer radius calculated in step S24, calculate the median of the inner radius and outer radius for all angles to obtain stable inner radius and outer radius.

5. The acquisition method according to claim 4, wherein The characteristic rings in step S2 are obtained through the following process: Based on the stable inner radius and outer radius in step S25, generate an annular mask image. The annular mask image is a binary image, with the annular area being 1 and the rest being 0: Among them, Mask(x,y) is the pixel value of the annular mask image, indicating whether the point (x,y) belongs to the annular region (1 means inside the ring, 0 means outside the ring); (x c ,y c ) is the center coordinate of the feature circular ring; rin is the inner radius of the feature circular ring; r out is the outer radius of the feature circular ring.

6. The acquisition method according to claim 1, wherein The mathematical model of the change of the different positions of the characteristic rings with the laser incident angle in step S3 is: x = aθ + b Where, x is the center position of the characteristic ring; θ is the incident angle; a is the proportionality factor, indicating the change rate of the center position x when θ changes; b is the offset, indicating the initial value of the center position x when the incident angle θ = 0; Use the least squares method to solve the sum of the squares of the minimum errors S, so as to determine a and b of the fitting curve equation of the position transformation of each characteristic ring; Where S is the sum of the squares of the minimum errors; n represents the number of data points, that is, the number of different incident angles simulated in step S2; x i represents the central position of the current feature ring at the i-th incident angle; θ i is the incident angle value at the i-th incident angle.

7. The acquisition method according to claim 1, wherein The average power density in step S3 satisfies the following formula: Among them, I avg (r) is the average power density corresponding to the pixels in all sector areas on the radius r; N is the number of pixel points in all sector areas on the radius r; I char_ring (x i , y i ) is the power density value at the independent image coordinate point (x i , y i ) of the characteristic circular ring, and x i , y i satisfies: In step S3, the Gaussian fitting uses a Gaussian function to fit the radial mean change curve I avg_fit (r) of each selected characteristic ring, and obtains a mathematical model of the internal power density of the characteristic ring varying with the radial distance: Among them, I avg_fit (r) is the power density value of the pixel points on the fitted radius r, a is the maximum amplitude, r is the radius, b is the center position of the Gaussian distribution, c is the standard deviation, and the least squares method is used for nonlinear fitting to solve the optimal parameters a, b, and c.

8. The acquisition method according to claim 1, wherein The calculation process of the peak power density described in step S3 is as follows: S31: Calculate the radial mean: For the simulation results at different angles of the same characteristic ring, set the start and end angles of the statistical sector area respectively, and remove the parts overlapping with other rings. For each radius r, r in ≤ r ≤ r out , and calculate the mean pixel value in all sector areas on this radius: To reduce the influence of noise, the method of removing extreme values and taking the mean is adopted: I trimmed (r) = trimmean(I avg (r), 10%) Among them, I trimmed (r) is the average value of the remaining data calculated after removing the highest and lowest 10% of the data in the pixel corresponding power density data set in all sector areas with radius r; S32: Extract the peak power density: For each incident angle θ i , take the peak power density from the radial mean value: I max (θ i ) = max{I trimmed_i (r)} where I max (θ i ) is the peak power density of the characteristic ring when the incident angle is θ i .

9. The acquisition method according to claim 1, wherein The mathematical model of the peak power density of the characteristic ring described in step S3 varying with the laser incident angle is: Among them, I max (θ) represents the peak power density of the characteristic ring at the incident angle of θ; A is the power density at the reference incident angle θ0; B is the attenuation coefficient; θ0 is the reference incident angle; The mathematical model of the peak power density of the characteristic ring described in step S3 varying with the laser transmission distance is: I(d) = I0·d -α Where, I(d) is the power density at the transmission distance d, I0 is the power density at the reference distance of 1 km, d is the transmission distance, and α is the attenuation index; The mathematical model of the peak power density of the characteristic ring described in step S3 varying with the laser emission power is: I(P) = k·P Among them, I(P) represents the power density when the laser emission power is P; k is a proportionality coefficient, and P is the laser emission power.

10. An application of a laser interference spot simulation model prepared by the acquisition method according to claim 1.

Images