A dynamic beam shaping method and system based on a spatial light modulator

CN121679915BActive Publication Date: 2026-10-23WUHAN JINDUN LASER TECH CO LTD +1
View PDF 2 Cites 0 Cited by

Patent Information

Application Number
CN202610054740.6
Authority / Receiving Office
CN · China
Patent Type
Patents(China)
Current Assignee / Owner
Filing Date
2026-01-15
Publication Date
2026-10-23
Estimated Expiration
2046-01-15

AI Technical Summary

Technical Problem

[0004]基于 SPGD 算法的无波前传感光束整形方案虽然结构简单,但仍存在以下缺点:第一,SPGD 算法作为并行优化算法,高维度的优化自变量会显著拖慢其收敛速度,而在光束整形中,相位图的控制变量往往需要采用几十甚至上百个泽尼克多项式的权重作为优化自变量

Benefits of technology

[0079]本发明的基于空间光调制器的动态光束整形系统及方法,先根据目标平顶光斑的参数和输出面相机像素尺寸生成目标强度分布图像,根据输入面和输出面坐标系计算所述目标平顶光斑的理想光束整形相位,并构建光束调制器的优化相位的表达式,然后通过输出面探测与反馈模块采集的光斑强度分布图像与目标强度分布图像构建复合性能指标表达式;在多阶段SPGD算法的每个阶段中,生成两组方向随机的扰动向量并分别对所述理想光束整形相位进行更新,将得到两组新的光束整形相位加载至光束调制器上,并生成两束新的激光束;最后采集两束新的激光束对应的光斑强度分布图像,并根据所述复合性能指标表达式分别计算其复合性能指标,兼顾了能量分布与形状特征,对所述优化相位表达式中的目标参数进行优化,并由所述光束调制器根据优化后的目标参数对激光束进行动态整形,从聚焦光学组件输出均匀的平顶光束。本发明通过 SPGD 算法结合泽尼克多项式作为控制变量,能够校正输出平顶光束的均匀性,因此可更好地适配光源,降低平顶光束整形的难度,提高系统的鲁棒性和工程应用性,同时在光束整形应用中具有更高的自由度,可实现多种平顶光束图案的生成与整形,对光源光束质量的要求较低,采用了多阶段逐级提升泽尼克多项式阶数的方式加速了平顶光束的校正过程,并以复合性能指标的收敛趋势作为阶段切换的判据;在超参数选择方面,设计了抗测量噪声的扰动策略与学习率自动选择方案,使算法在实际应用中能够快速稳定地收敛;优化器应用方面,采用Adam优化器加速整个优化过程,同时针对多阶段的算法应用对优化器的行为进行了调整。

✦ Generated by Eureka AI based on patent content.

Smart Images

  • Figure CN121679915B_ABST
    Figure CN121679915B_ABST
Patent Text Reader

Abstract

The application relates to a dynamic beam shaping system and method based on a spatial light modulator, which comprises a laser light source module, a beam expansion and collimation optical path, a beam modulator, a focusing optical assembly, a detection and feedback module and a control and optimization module. The application can correct the uniformity of the output flat-top beam by means of a stochastic parallel gradient descent (SPGD) algorithm and in combination with Zernike polynomials as control variables. Meanwhile, a composite performance index considering the energy distribution and shape characteristics is constructed, the correction process of the flat-top beam is accelerated by adopting a multi-stage gradual promotion of the Zernike polynomial order, the convergence trend of the performance index is taken as the judgment basis for stage switching, the dimension expansion of the independent variable is implemented in stages, and the convergence speed is significantly improved. In addition, the adaptive selection of the disturbance amplitude and the learning rate parameter in the SPGD algorithm is realized by using noise estimation and local linearity index, and the robustness of the system to noise and environmental disturbance is enhanced.
Need to check novelty before this filing date? Find Prior Art

Description

Technical Field

[0001] This invention relates to the field of laser beam shaping technology, and in particular to a dynamic beam shaping method and system based on a spatial light modulator. Background Technology

[0002] Flat-top beam shaping has wide applications in ultrafast laser processing, microscopic illumination, and optical tweezers. Commonly used methods to obtain flat-top beams include diffractive optical elements, microlens arrays, freeform surface optical elements, and metasurfaces. However, traditional flat-top beam shaping devices are designed for specific light sources and cannot dynamically adjust the output. If the light source is changed or the beam quality degrades, they may malfunction. Compared to fixed-structure beam shaping devices, spatial light modulators, as dynamic beam modulators, have the advantage of programmability. They can generate various flat-top beams by loading phase maps and can compensate for incident beam distortion to a certain extent.

[0003] In existing technologies, methods that directly measure the incident light wavefront distortion using a wavefront sensor and apply phase compensation require adding a beam splitter and a wavefront sensor to the optical path. This not only reduces system energy utilization but also increases system complexity and cost, hindering its widespread application in industrial laser processing and other similar scenarios. Therefore, wavefront-free sensing technology, which measures the light intensity distribution at the output surface using a camera and then updates the phase map using the Stochastic Parallel Gradient Descent (SPGD) algorithm, is more widely used.

[0004] While the wavefront-less beamforming scheme based on the SPGD algorithm is structurally simple, it still suffers from the following drawbacks: First, as a parallel optimization algorithm, the high-dimensional optimization variables of SPGD significantly slow down its convergence speed. Furthermore, in beam shaping, the control variables of the phase diagram often require weights of tens or even hundreds of Zernike polynomials as optimization variables. Second, in practical optical systems, the measurement of algorithm performance indicators is affected by various noises such as camera readout noise, laser power fluctuations, and mechanical vibrations. Therefore, the selection of hyperparameters for the SPGD algorithm requires an interference-resistant automatic calibration method to ensure the convergence speed and stability of the algorithm during operation. Summary of the Invention

[0005] The technical problem to be solved by the present invention is to provide a dynamic beam shaping system and method based on a spatial light modulator, which addresses the shortcomings of the prior art.

[0006] The technical solution of the present invention to solve the above-mentioned technical problems is as follows: A dynamic beam shaping system based on a spatial light modulator, comprising a laser source module, a beam expanding and collimating optical path, a beam modulator, a focusing optical component, a detection and feedback module, and a control and optimization module connected in sequence;

[0007] The laser source module is located at the input surface of the system's optical path and is used to output a laser beam;

[0008] The beam expanding and collimating optical path is used to expand and collimate the laser beam.

[0009] The beam modulator is used to modulate the laser beam after it has been shaped and collimated.

[0010] The focusing optical component is used to form a flat-top beam of the target shape on the focal plane or the target surface after the modulation of the laser beam.

[0011] The detection and feedback module is located on the output surface of the system's optical path and is used to acquire the light spot intensity distribution image of the flat-top beam and transmit it to the control and optimization module.

[0012] The control and optimization module is used to update the phase distribution of the beam modulator in real time according to the light spot intensity distribution image, so as to correct the output flat-top beam.

[0013] This invention also provides a dynamic beam shaping method based on a spatial light modulator, employing the aforementioned dynamic beam shaping system based on a spatial light modulator, comprising the following steps:

[0014] Generate a target intensity distribution image based on the parameters of the target flat-top light spot and the output camera pixel size. The ideal beam-shaping phase of the target flat-top spot is calculated based on the input and output plane coordinate systems. And construct the optimized phase of the beam modulator. The expression;

[0015] The light spot intensity distribution image acquired by the output surface detection and feedback module Image of target intensity distribution Construct a composite performance index expression;

[0016] In each stage of the multi-stage SPGD algorithm, two sets of random perturbation vectors are generated and the ideal beam shaping phase is updated respectively. The resulting two sets of new beam shaping phases are loaded onto the beam modulator and two new laser beams are generated.

[0017] The intensity distribution images of the two new laser beams are acquired, and their composite performance indices are calculated according to the composite performance index expression. The target parameters in the optimized phase expression are optimized, and the beam modulator dynamically shapes the laser beams according to the optimized target parameters, so as to output a uniform flat-top beam from the focusing optics.

[0018] Based on the above technical solution, the present invention can be further improved as follows:

[0019] Further: The ideal beam shaping phase of the target flat-top light spot is calculated based on the input and output plane coordinate systems. And construct the optimized phase of the beam modulator. The expression specifically includes the following steps:

[0020] The shape, size, and energy distribution range of the target flat-top light spot are defined, and a corresponding target intensity image is generated by combining this with the pixel size of the output camera. ;

[0021] Ideal beam shaping phase for achieving the target flat-top light spot is pre-calculated based on Fourier optics and the standing phase method. ,in The physical coordinate system of the input plane;

[0022] Calculate Zernike polynomials And construct the optimized phase of the beam modulator. expression:

[0023]

[0024]

[0025]

[0026]

[0027]

[0028]

[0029] Where p is the sorting index of the Zernike polynomial excluding constant terms (n=0, m=0). The polar coordinate system is the input plane of the system's optical path. Here, represents the orthogonal normalization coefficients, n is the radial order of the Zernike polynomial, and m is the azimuth frequency of the Zernike polynomial. Let be a radial polynomial of Zernike polynomials, with order (n, m). Let be the objective parameters of the Zernike polynomial to be optimized. Let be a vector composed of P Zernike polynomials.

[0030] Further: the light spot intensity distribution image acquired through the output surface detection and feedback module Image of target intensity distribution Constructing a composite performance index expression specifically includes the following steps:

[0031] Calculate the light spot intensity distribution image Image of target intensity distribution Normalized mean square error terms between The calculation formula is:

[0032]

[0033] in, Let F be the first balancing parameter used to balance the order-of-magnitude difference between the two images caused by the measurement, and let F be the Frobenius norm. The inner product of two image matrices;

[0034] Calculate the structural similarity between the light spot intensity distribution image and the target intensity distribution image. The calculation formula is:

[0035]

[0036] in, Image of light spot intensity distribution Image of target intensity distribution covariance, Image of light spot intensity distribution standard deviation For target intensity distribution image standard deviation This is a tiny constant used to avoid calculation errors caused by dividing the function by zero;

[0037] According to the normalized mean square error term With structural similarity Construct a composite performance index expression:

[0038]

[0039] in, This is a second balancing parameter used to balance the order-of-magnitude difference between the normalized mean square error term and the structural similarity.

[0040] Further: The step of generating two sets of randomly oriented perturbation vectors and updating the phase of the ideal beam shaping accordingly specifically includes the following steps:

[0041] Generate a set of perturbation vectors with random positive and negative values ​​but uniform size. :

[0042]

[0043] in, It is a length of Vectors that satisfy the Rademacher distribution Let be the amplitude of the disturbance in stage s;

[0044] According to a set of disturbance vectors and the target parameter vector for the current stage The ideal beam shaping phase is updated to obtain new beam shaping phases with two opposite directions under perturbation. :

[0045]

[0046] in, Let be the number of Zernike polynomials in the s-th stage. Let p be the target parameter with sorting index p and iteration number t in stage s. Let p be the perturbation with sorting index p and iteration number t in stage s. Let p be the p-th Zernikal polynomial.

[0047] Further: The optimization of the target parameter in the optimized phase expression specifically includes the following steps:

[0048] Numerical gradient estimated using composite performance index of the intensity distribution images of the two new laser beams. The calculation formula is:

[0049]

[0050] in, It is a composite performance index of the spot intensity distribution image corresponding to the new beam shaping phase under positive perturbation vector perturbation. It is a composite performance index of the spot intensity distribution image corresponding to the new beam shaping phase under negative perturbation vector perturbation. Let be the disturbance amplitude in stage s. Direction vector The p-th value;

[0051] Based on the numerical gradient Composition of numerical gradient vector The Adam optimizer is used to optimize the target parameters in the optimized phase expression.

[0052] Further: the numerical gradient vector The optimization of the target parameters in the optimized phase expression using the Adam optimizer specifically includes the following steps:

[0053] Set the number of preheating iterations When the current iteration number is less than or equal to the preheating iteration number Then, the objective parameters in the optimized phase expression are updated using the gradient descent method, and the calculation formula is as follows:

[0054]

[0055] in, Let be the target parameter vector at the (t+1)th iteration in the S-th stage. Let be the target parameter vector at the t-th iteration. The learning rate for stage S. Let be the numerical gradient vector at the t-th iteration;

[0056] When the current iteration number is greater than the preheating iteration number At that time, the first-order momentum at the t-th iteration of the Adam optimizer is used. With second momentum The target parameters in the optimized phase expression are updated using the following formula:

[0057]

[0058] in, These are the decay coefficients of the Adam optimizer. It is a constant. , These are the first-order and second-order momentum of the normalized Adam optimizer, respectively;

[0059] Read the most recent The composite performance index is calculated based on the number of iterations, and the change in the number of iterations is fitted using least squares. The slope of the fitted curve is then calculated. and in The ratio of the initial slope to the initial slope is less than a preset proportional threshold. Then proceed to the next stage, repeating the above steps until all stages are completed and the final target parameters are obtained.

[0060] Furthermore, before each stage begins, the following steps are also included:

[0061] Regarding the disturbance amplitude mentioned in the previous stage and learning rate Adaptive calibration is performed to obtain the optimal perturbation amplitude for each stage. and learning rate .

[0062] Further: the disturbance amplitude of the previous stage and learning rate Adaptive calibration includes the following steps:

[0063] Repeated measurements under the same phase diagram Subcomposite performance index measurement value , where ε is zero-mean noise, and the variance of the noise is estimated. ;

[0064] Based on the measured values ​​of the composite performance index The Local Linearity Index (LLI) is calculated using the following formula:

[0065]

[0066]

[0067]

[0068]

[0069] in, These are the m-th random direction vectors. The measured values ​​of the composite performance index corresponding to the positive and negative candidate perturbation amplitudes. These are the initial target parameters used in the current optimization phase. The mean of the composite performance index measurements obtained from the corresponding M0 measurements, for The standard deviation of measurement For the median operation, This refers to the local linearity calculated using the m-th composite performance index measurement value under the current disturbance amplitude;

[0070] The optimal disturbance amplitude is determined based on the local linearity index LLI. and based on Adaptively set the learning rate.

[0071] Further: the local linearity index Determine the optimal disturbance amplitude, and based on Adaptive learning rate setting specifically includes the following steps:

[0072] In the first stage of the multi-stage SPGD algorithm, the initial value of the perturbation amplitude is... And in subsequent stages, the initial value of the disturbance amplitude will be... Set to the disturbance amplitude of the previous stage ;

[0073] Search for the local linearity index (LLI) under the current perturbation amplitude and check if it exceeds the threshold. If the threshold is exceeded, the Local Linearity Index (LLI) under the current perturbation amplitude is reduced to half of its original value; conversely, if the threshold is exceeded, the LLI is increased to twice its original value, until a set of values ​​at the threshold values ​​is obtained. Up and down and closest to the threshold Disturbance amplitude ;

[0074] Based on the disturbance amplitude The optimal disturbance amplitude was calculated. The calculation formula is:

[0075]

[0076] based on and optimal perturbation amplitude Adaptive learning rate setting The calculation formula is:

[0077]

[0078] in, For use in controlling the learning rate The corresponding target parameter update magnitude is a normal number within the expected range.

[0079] The present invention discloses a dynamic beam shaping system and method based on a spatial light modulator. First, a target intensity distribution image is generated based on the parameters of the target flat-top spot and the pixel size of the output surface camera. Then, the ideal beam shaping phase of the target flat-top spot is calculated according to the input and output surface coordinate systems, and an expression for the optimized phase of the beam modulator is constructed. Next, a composite performance index expression is constructed using the spot intensity distribution image acquired by the output surface detection and feedback module and the target intensity distribution image. In each stage of the multi-stage SPGD algorithm, two sets of randomly oriented perturbation vectors are generated and updated respectively on the ideal beam shaping phase. The resulting two new beam shaping phases are loaded onto the beam modulator, generating two new laser beams. Finally, spot intensity distribution images corresponding to the two new laser beams are acquired, and their composite performance indices are calculated according to the composite performance index expression, taking into account both energy distribution and shape characteristics. The target parameters in the optimized phase expression are optimized, and the beam modulator dynamically shapes the laser beams according to the optimized target parameters, outputting a uniform flat-top beam from the focusing optics. This invention uses the SPGD algorithm combined with Zernike polynomials as control variables to correct the uniformity of the output flat-top beam. This allows for better adaptation to the light source, reduces the difficulty of flat-top beam shaping, and improves the system's robustness and engineering applicability. It also offers greater freedom in beam shaping applications, enabling the generation and shaping of various flat-top beam patterns with lower requirements for light source beam quality. A multi-stage, progressively increasing Zernike polynomial order accelerates the flat-top beam correction process, using the convergence trend of composite performance indicators as the criterion for stage switching. Regarding hyperparameter selection, a disturbance strategy to resist measurement noise and an automatic learning rate selection scheme are designed, enabling the algorithm to converge quickly and stably in practical applications. For the optimizer application, the Adam optimizer is used to accelerate the entire optimization process, and the optimizer's behavior is adjusted for multi-stage algorithm applications. Attached Figure Description

[0080] Figure 1 This is a schematic diagram of a dynamic beam shaping system based on a spatial light modulator according to an embodiment of the present invention;

[0081] Figure 2 This is a schematic diagram of the structure of a simulated optical system according to an embodiment of the present invention;

[0082] Figure 3A This is a schematic diagram comparing the simulation results of single-stage and multi-stage SPGD algorithms according to an embodiment of the present invention.

[0083] Figure 3B This is a schematic diagram comparing the simulation convergence curves of single-stage and multi-stage SPGD algorithms according to an embodiment of the present invention.

[0084] Figure 4AThis is a schematic diagram of the light spot formed by the laser beam emitted from the laser source module in a dynamic beam shaping system based on a spatial light modulator according to an embodiment of the present invention.

[0085] Figure 4B This is a schematic diagram comparing the uniformity of the flat-top beam before and after dynamic beam shaping according to an embodiment of the present invention. Detailed Implementation

[0086] The principles and features of the present invention are described below with reference to the accompanying drawings. The examples given are only for explaining the present invention and are not intended to limit the scope of the present invention.

[0087] like Figure 1 As shown, a dynamic beam shaping system based on a spatial light modulator includes a laser source module, a beam expanding and collimating optical path, a beam modulator, a focusing optical component, a detection and feedback module, and a control and optimization module connected in sequence.

[0088] The laser source module is located at the input surface of the system's optical path and is used to output a laser beam;

[0089] Here, the output beam of the laser source module can be a near-Gaussian distribution, an elliptical Gaussian distribution, or other complex amplitude distribution, and errors such as beam center offset and beam waist diameter size error are allowed.

[0090] The beam expanding and collimating optical path is used to expand and collimate the laser beam.

[0091] In embodiments of the present invention, the beam expanding and collimating optical path includes components such as a beam expander, a polarizer, and a reflector, used to expand and collimate the laser beam.

[0092] The beam modulator is used to modulate the laser beam after it has been shaped and collimated.

[0093] Optionally, in one or more embodiments of the present invention, the beam modulator employs a phase-type spatial light modulator (SLM).

[0094] The focusing optical component is used to form a flat-top beam of the target shape on the focal plane or the target surface after the modulation of the laser beam.

[0095] In embodiments of the present invention, the focusing optical assembly comprises a 4f system and one or more objective lenses and lenses. The 4f system serves to filter the 0th-order light of the SLM and to reduce the beam size. The objective lens then shapes the beam. Simultaneously, the objective lens and the lens in front of the camera in the detection and feedback module together form a magnifying optical path for observing the shaped light spot. The focusing optical assembly can form a flat-top beam of the target shape on the focal plane or target surface after the beam modulated by the SLM.

[0096] The detection and feedback module is located on the output surface of the system's optical path and is used to acquire the light spot intensity distribution image of the flat-top beam and transmit it to the control and optimization module.

[0097] In embodiments of the present invention, the detection and feedback module includes a CCD, a CMOS camera, etc., which transmits the acquired light spot intensity distribution image to the control and optimization module.

[0098] The control and optimization module is used to update the phase distribution of the beam modulator in real time according to the light spot intensity distribution image, so as to correct the output flat-top beam.

[0099] In embodiments of the present invention, the control and optimization module is implemented by a computer or embedded processor, and has a built-in multi-stage SPGD algorithm module to update the phase distribution of the SLM in real time according to the acquired pattern.

[0100] This invention also provides a dynamic beam shaping method based on a spatial light modulator, employing the aforementioned dynamic beam shaping system based on a spatial light modulator, comprising the following steps:

[0101] S1: Generate a target intensity distribution image based on the parameters of the target flat-top light spot and the output camera pixel size. The ideal beam-shaping phase of the target flat-top spot is calculated based on the input and output plane coordinate systems. And construct the optimized phase of the beam modulator. The expression;

[0102] S2: Spot intensity distribution image acquired by the output surface detection and feedback module Image of target intensity distribution Construct a composite performance index expression;

[0103] S3: In each stage of the multi-stage SPGD algorithm, two sets of random perturbation vectors are generated and the ideal beam shaping phase is updated respectively. The resulting two sets of new beam shaping phases are loaded onto the beam modulator and two new laser beams are generated.

[0104] S4: Acquire the spot intensity distribution images corresponding to the two new laser beams, calculate their composite performance index according to the composite performance index expression, optimize the target parameters in the optimized phase expression, and have the beam modulator dynamically shape the laser beams according to the optimized target parameters to output a uniform flat-top beam from the focusing optical component.

[0105] In one or more embodiments of the present invention, the ideal beam shaping phase of the target flat-top spot is calculated based on the input and output surface coordinate systems. And construct the optimized phase of the beam modulator. The expression specifically includes the following steps:

[0106] S11: Set the shape, size, and energy distribution range of the target flat-top light spot, and generate a corresponding target intensity image based on the output camera pixel size. ;

[0107] For example, if the light intensity within the edge region of a square with a side length of 2R is set to be the same, then the target intensity image... for:

[0108]

[0109] S12: Ideal beam shaping phase for achieving the target flat-top light spot based on Fourier optics and the standing phase method through pre-calculation. ,in The physical coordinate system of the input plane;

[0110] Here, for the one-dimensional case, the ideal flat-top spot generation phase is:

[0111] ,

[0112]

[0113] in Let λ be the beam waist radius of the incident Gaussian beam, and R be the radius of the uniformly distributed target. λ is the incident light wavelength, and f is the focal length of the lens. If the target pattern is a one-dimensional separable pattern (such as a square), the generation phase of the two corresponding one-dimensional beams can be obtained by superimposing their shaping phases in the vertical direction. .

[0114] S13: Calculate the Zernik polynomial And construct the optimized phase of the beam modulator. expression:

[0115]

[0116]

[0117]

[0118]

[0119]

[0120]

[0121] Where p is the sorting index of the Zernike polynomial excluding constant terms (n=0, m=0). The polar coordinate system is the input plane of the system's optical path. Here, represents the orthogonal normalization coefficients, n is the radial order of the Zernike polynomial, and m is the azimuth frequency of the Zernike polynomial. Let be a radial polynomial of Zernike polynomials, with order (n, m). Let be the objective parameters of the Zernike polynomial to be optimized. Let be a vector composed of P Zernike polynomials.

[0122] Here, the order of the Zernike polynomial is determined by the radial order n and the azimuth frequency m. p is the noll sorting index of the Zernike polynomial excluding the constant term (n=0, m=0).

[0123] It should be noted that in step S1, the optimization variable allocation method of the multi-stage SPGD algorithm needs to be preset. The multi-stage SPGD algorithm has a total of K stages, and the number of Zernike polynomials used in each stage satisfies the following condition: ,For example This means that the first stage only uses the first 14 Zernike polynomials for optimization, the second stage expands to 65, the third stage expands to 135, and so on.

[0124] In one or more embodiments of the present invention, the light spot intensity distribution image acquired by the output surface detection and feedback module Image of target intensity distribution Constructing a composite performance index expression specifically includes the following steps:

[0125] S21: Calculate the light spot intensity distribution image Image of target intensity distribution Normalized mean square error terms between The calculation formula is:

[0126]

[0127] in, Let F be the first balancing parameter used to balance the order-of-magnitude difference between two images caused by the measurement, and let F be the Frobenius norm. For a matrix, the F-norm is defined as the square root of the sum of the squares of each element in the matrix. The inner product of two image matrices;

[0128] Here, the normalized mean square error term Used to compare the light spot intensity distribution images acquired by the output surface detection and feedback module. and target image The difference between each pixel is balanced by α to account for the order-of-magnitude difference between the two images caused by the measurement.

[0129] S22: Calculate the structural similarity between the spot intensity distribution image and the target intensity distribution image. The calculation formula is:

[0130]

[0131] in, Image of light spot intensity distribution Image of target intensity distribution covariance, Image of light spot intensity distribution standard deviation For target intensity distribution image standard deviation This is a tiny constant used to avoid calculation errors caused by dividing the function by zero;

[0132] In this embodiment, structural similarity It is used to compare structural differences between two images and is not sensitive to brightness, contrast, or differences between individual pixels.

[0133] S23: Based on the normalized mean square error term With structural similarity Construct a composite performance index expression:

[0134]

[0135] in, This is a second balancing parameter used to balance the order-of-magnitude difference between the normalized mean square error term and the structural similarity, and can be set through simulation or pre-experimentation.

[0136] Below, we will take the s-th stage (using the first P) of the multi-stage SPGD algorithm as an example. s Taking a Zernike polynomial as an example, we will introduce the process of the multi-stage SPGD algorithm in detail.

[0137] In one or more embodiments of the present invention, the step of generating two sets of randomly oriented perturbation vectors and updating the ideal beam shaping phase respectively includes the following steps:

[0138] S31: Generate a set of perturbation vectors with random positive and negative values ​​but uniform size. :

[0139]

[0140] in, It is a length of Vectors that satisfy the Rademacher distribution This is the amplitude of the disturbance in stage s, which is also a parameter that needs to be automatically calibrated before the start of each subsequent stage.

[0141] S32: Based on a set of the aforementioned perturbation vectors and the target parameter vector for the current stage The ideal beam shaping phase is updated to obtain new beam shaping phases with two opposite directions under perturbation. :

[0142]

[0143] in, Let be the number of Zernike polynomials in the s-th stage. Let p be the target parameter with sorting index p and iteration number t in stage s. Let p be the perturbation with sorting index p and iteration number t in stage s. Let p be the p-th Zernikal polynomial.

[0144] Because of the added perturbation, the ideal beam-shaping phase will change; therefore, this is addressed here by using the perturbation vector. and the target parameter vector for the current stage The ideal beam shaping phase is updated, allowing the output of the beam modulator to be adjusted based on the updated beam shaping phase, and the corresponding output beam spot to be acquired. and And calculate the corresponding composite performance index according to step S23 above. and .

[0145] In one or more embodiments of the present invention, optimizing the target parameter in the optimized phase expression specifically includes the following steps:

[0146] S41: Estimate the numerical gradient using the composite performance index of the spot intensity distribution images corresponding to the two new laser beams. The calculation formula is:

[0147]

[0148] in, It is a composite performance index of the spot intensity distribution image corresponding to the new beam shaping phase under positive perturbation vector perturbation. It is a composite performance index of the spot intensity distribution image corresponding to the new beam shaping phase under negative perturbation vector perturbation. Let be the disturbance amplitude in stage s. Direction vector The p-th value;

[0149] S42: Based on the numerical gradient Composition of numerical gradient vector The Adam optimizer is used to optimize the target parameters in the optimized phase expression.

[0150] Specifically, in one or more embodiments of the present invention, the step based on the numerical gradient vector The optimization of the target parameters in the optimized phase expression using the Adam optimizer specifically includes the following steps:

[0151] S421: Set the number of preheating iterations When the current iteration number is less than or equal to the preheating iteration number Then, the objective parameters in the optimized phase expression are updated using the gradient descent method, and the calculation formula is as follows:

[0152]

[0153] in, Let be the target parameter vector at the (t+1)th iteration in the S-th stage. Let be the target parameter vector at the t-th iteration. The learning rate for stage S. Let be the numerical gradient vector at the t-th iteration;

[0154] Here, to avoid slow convergence due to Adam's suppression of large gradient calculations in the first stage (s=1), this invention sets the number of preheating iterations. During warm-up (i.e., the number of iterations) At this point, the parameters are updated using simple gradient descent.

[0155] S422: When the current iteration number is greater than the preheating iteration number At that time, the first-order momentum at the t-th iteration of the Adam optimizer is used. With second momentum The target parameters in the optimized phase expression are updated using the following formula:

[0156]

[0157] in, These are the decay coefficients for the Adam optimizer, initially set to 0.9 and 0.99, respectively. It is a constant, taking the value 10. -8 To prevent the occurrence of extremely small values ​​when divided by 0, , These are the first and second momentum of the normalized Adam optimizer, respectively, to prevent excessive errors in gradient estimation during initialization; this ensures convergence stability in the initial stage while improving convergence speed in later stages.

[0158] S423: Read the most recent The composite performance index is calculated based on the number of iterations, and the change in the number of iterations is fitted using least squares. The slope of the fitted curve is then calculated. As iterations continue, the curve gradually converges, and... The ratio of the initial slope to the initial slope is less than a preset proportional threshold. (In this embodiment, the value is 0.3) When the value is reached, proceed to the next stage, repeat the above steps, until all K stages are completed, and the final target parameter is obtained. Here, the most recent value is read. When calculating the composite performance index for the number of iterations, the composite performance index from the previous stage is not included. Furthermore, when switching stages, the dimension of the variable parameters is first expanded to... The newly added parameters can be initialized to 0, and the second-order momentum of the Adam optimizer will be reset to 0 to enhance its adaptability to the new dimensions.

[0159] Preferably, in one or more embodiments of the present invention, the following steps are further included before the start of each stage:

[0160] S420: Regarding the disturbance amplitude mentioned in the previous stage and learning rate Adaptive calibration is performed to obtain the optimal perturbation amplitude for each stage. and learning rate This approach balances noise robustness with the local linear approximation of numerical gradient estimation.

[0161] Specifically, in one or more embodiments of the present invention, the disturbance amplitude of the previous stage and learning rate Adaptive calibration includes the following steps:

[0162] S4201: Repeat measurements under the same phase diagram Subcomposite performance index measurement value , where ε is zero-mean noise, and the variance of the noise is estimated. ;

[0163] Here, we assume the performance index measurement value is , where ε is zero-mean noise. The variance of the noise. This can be achieved by repeatedly acquiring data under the same phase map. Sub-performance index estimation. The larger the disturbance amplitude, the greater the change in the evaluation index, and the higher the signal-to-noise ratio.

[0164] S4202: Based on the measured values ​​of the composite performance indicators Calculate the local linearity index (Local Linearity Index, LLI), calculated using the following formula:

[0165]

[0166]

[0167]

[0168]

[0169] in, These are the m-th random direction vectors. The measured values ​​of the composite performance index corresponding to the positive and negative candidate perturbation amplitudes. These are the initial target parameters used in the current optimization phase. The mean value of the composite performance index obtained from the corresponding M0 measurements is used for local linearity estimation. for The standard deviation of measurement To calculate the standard deviation of the system's detection noise, For the median operation, This refers to the local linearity calculated using the m-th composite performance index measurement value under the current disturbance amplitude;

[0170] First, set M random direction vectors. For each candidate perturbation amplitude Measured M times and obtained a set of evaluation indicators , for The measurement standard deviation can be obtained from the results of the previous measurement using the error propagation theorem. LLI represents the final assessment of local linearity under the current perturbation size. A larger LLI indicates poorer local linearity and less reliable numerical gradient estimation under the current perturbation, requiring a reduction in the perturbation amplitude. However, due to measurement noise, LLI evaluation may be inaccurate, leading to very large LLI measurements even under small perturbations. Therefore, an additional factor is added to the denominator of the LLI measurement formula. It can achieve a 95% confidence level under the noise estimation conditions of the Gaussian model to avoid the denominator of LLI approaching 0.

[0171] The local linearity index (LLI) measures the nonlinearity of performance index J, and this index changes with the amplitude of the disturbance. The increase of [variable name] leads to an increase in [variable name]. During gradient estimation, the smaller the perturbation amplitude, the more accurate the gradient estimate in that region. In the formula for calculating the Local Linearity Index (LLI), the numerator corresponds to [variable name] in the current [variable name]. The second-order remainder of the Taylor expansion nearby, the denominator It is a first-order difference term.

[0172] S4203: Determine the optimal disturbance amplitude based on the local linearity index LLI. and based on Adaptively set the learning rate.

[0173] The selection of perturbation amplitude in the SPGD algorithm needs to consider both signal-to-noise ratio and local linearity. Therefore, a bisection method is designed to automatically search for the optimal perturbation amplitude.

[0174] In one or more embodiments of the present invention, the optimal perturbation amplitude is determined based on the Local Linearity Index (LLI), and based on... Adaptive learning rate setting specifically includes the following steps:

[0175] S42031: Initial value of the perturbation amplitude in the first stage of the multi-stage SPGD algorithm. And in subsequent stages, the initial value of the disturbance amplitude will be... Set to the disturbance amplitude of the previous stage ;

[0176] S42032: Search for the local linearity index (LLI) under the current perturbation amplitude and whether it exceeds the threshold. (In this embodiment, the value is 0.1). If the threshold is exceeded, the local linearity index (LLI) under the current disturbance amplitude is reduced to half of its original value; conversely, the local linearity index (LLI) under the current disturbance amplitude is increased to twice its original value, until a set of values ​​at the threshold is obtained. Up and down and closest to the threshold Disturbance amplitude ;

[0177] S42033: Based on the disturbance amplitude The optimal disturbance amplitude was calculated. The calculation formula is:

[0178]

[0179] based on and optimal perturbation amplitude Adaptive learning rate setting The calculation formula is:

[0180]

[0181] in, For use in controlling the learning rate The corresponding target parameter update magnitude is a normal number within the expected range.

[0182] The dynamic beam shaping system and method based on a spatial light modulator of the present invention has the following advantages:

[0183] a) This system does not rely on wavefront sensors; it only measures the intensity distribution of the output surface through a camera and constructs scalar performance indicators to achieve wavefront sensor-independent closed-loop beam shaping.

[0184] b) Multi-stage Zernikal polynomial expansion strategy. By controlling the order of the Zernikal polynomials involved in the optimization in stages, a small number of low-order optimizations are used to guide the rapid convergence of high-dimensional independent variables.

[0185] c) A multi-stage SPGD algorithm combining stochastic parallel gradients and Adam momentum. This algorithm uses stochastic perturbations to estimate gradients and is simple to implement in hardware. By introducing optimizer warm-up in the first stage and second-order momentum reset during stage switching, the convergence efficiency and algorithm robustness are improved.

[0186] d) Adaptive calibration mechanism for perturbation amplitude and learning rate. The LLI index is systematically introduced to quantitatively characterize the degree of nonlinearity at a given perturbation amplitude. A bisection search strategy is employed to automatically measure the perturbation amplitude and learning rate, avoiding empirical parameter tuning and enabling the algorithm to exhibit good adaptability under different target patterns and noise levels.

[0187] e) A composite performance index was designed, which combines the normalized mean square error with structural similarity. This constrains the overall energy distribution while emphasizing the overall shape similarity, making it more global than a single energy error constraint.

[0188] f) Automatic stage switching strategy driven by stage convergence trend. The slope of the fitted line of the performance index is used as an evaluation of the convergence trend to determine whether to switch to the next stage. Compared with optimization with a fixed number of iterations, this avoids wasting convergence times or premature termination.

[0189] g) The versatility and scalability of the device structure. The wavefront controller can be an SLM or a digital micromirror device (DMD), etc.

[0190] The following examples illustrate the dynamic beam shaping system and method based on spatial light modulators of this invention:

[0191] Example 1:

[0192] 1. Simulation Model Description:

[0193] Simulated optical systems are as follows Figure 2The optical path is shown. The pixel width of the input surface is set to 12.5 μm, and the lens parameters are set to focal length f = 10 mm and numerical aperture NA = 0.45. It is assumed that the laser beam at the input surface is a non-ideal elliptic Gaussian beam, such as... Figure 3A As shown: The elliptic Gaussian beam has an ellipticity of 0.6, an effective Gaussian radius of r0 = 2.1 mm, and a wavefront phase with aberrations randomly composed of the first 12 Zernike modes. The desired output is a square flat-top beam with a half-width of R = 0.1 mm, and the initial phase can be simulated. Light intensity distribution on the output surface, such as Figure 3A The initial values ​​are shown.

[0194] 2. Algorithm parameter settings:

[0195] This example compares a single-stage optimization strategy (S-ASPGD) using the Adam optimizer with the multi-stage optimization strategy (M-ASPGD) proposed in this invention. To ensure a fair comparison of the optimization processes, the Adam preheating startup method is also adopted in the S-ASPGD scheme. And the same hyperparameter selection strategy as M-ASPGD, M=8.

[0196] Other settings for the M-ASPGD algorithm: Evaluation metric: downward trend; Evaluation range: W = 40; Zernike polynomial order used in each stage. .

[0197] The S-ASPGD algorithm uses a P = 65 Zernike polynomial in a single operation.

[0198] 3. Performance Indicators and Comparison Results:

[0199] like Figure 3B As shown, the performance index of the M-ASPGD algorithm is close to convergence when the number of Zernike polynomials is increased to 65 in the third stage. However, the S-ASPGD algorithm is still being optimized at this time, and its performance index is worse than that of M-ASPGD.

[0200] The uniformity of the output flat-top beam is evaluated by the following two metrics:

[0201] Peak-valley values:

[0202]

[0203] And the flat-top uniformity as defined in ISO standard 13694:

[0204]

[0205] in and These are the maximum and minimum light intensity values ​​within the flat-top profile. The energy density histogram N(I) represents the energy density histogram. x,y ) middle The full width at half maximum (FWHM) of the nearest dominant peak. Output energy is represented by the measured intensity I, while N(I) x,y ) is an energy density of I x,y The number of image pixels. A perfectly uniform flat-top beam of U... p =0.

[0206] like Figure 3A As shown, S-ASPGD and M-ASPGD performed PV=0.62 and U for the same correction task. p =0.07 and PV=0.33, U p =0.05.

[0207] Example 2:

[0208] 1. Optical system structure:

[0209] The optical setup for this experiment is shown in the diagram below. Figure 1 As shown, it consists of an SLM with a pixel size of 12.5 μm × 12.5 μm and a continuous laser with a wavelength of λ = 532 nm. Figure 4A The image shows the laser spot, which is an elliptic Gaussian beam with an ellipticity of 0.8 and a maximum diameter d along its major axis. max = 1.48 mm, minimum diameter in the minor axis direction is d min = 1.19 mm. The beam is first filtered by a linear polarizer, then expanded by a 3x beam expander, and subsequently modulated by an SLM. To suppress the 0th order of the SLM, a 4th order grating phase is superimposed on each beam-shaping phase. The SLM-modulated beam is imaged by a 4f system, and the 0th order component is blocked at the Fourier plane by an aperture stop; simultaneously, the 4f system reduces the beam to a factor of 0.5 to fit the entrance pupil of the objective lens. After passing through the 4f system, the beam is split at a beam splitter and fed into the objective lens (numerical aperture NA = 0.25, focal length f1 = 16.5 mm). A mirror is placed at the focal plane of the objective lens, and the reflected light is magnified by the objective lens and a lens (focal length f2 = 100 mm) before being detected by a camera with a pixel size of 2.4 μm × 2.4 μm.

[0210] 2. Algorithm parameter settings:

[0211] Noise estimation related parameter M0=20, other hyperparameters are selected the same as in simulation M=8, W warm =W=40.

[0212] The phased strategy of the M-ASPGD algorithm is as follows: .

[0213] 3. Experimental phenomena and results:

[0214] like Figure 4B As shown, the target pattern is set as a square flat-top beam with a half-width of R=0.05 mm. Before dynamic shaping, PV=0.94, U p =0.38, the uniformity after dynamic shaping is PV=0.56, U p =0.29, average time taken: 20 minutes.

[0215] The above experimental results demonstrate that the dynamic beam shaping method based on a spatial light modulator of this invention...

[0216] a) The convergence speed is faster than traditional single-stage SPGD or single-stage multi-stage SPGD algorithms using the Adam optimizer. This invention employs staged expansion of the independent variable dimension, pre-startup of the optimizer, and staged second-order momentum reset, resulting in a significant improvement in convergence speed with the same number of iterations.

[0217] b) By using noise estimation and local linearity index, the perturbation amplitude and learning rate parameters in the multi-stage SPGD algorithm can be adaptively selected, which improves the robustness to noise and environmental perturbations.

[0218] c) No wavefront sensor required, system is simple and energy-efficient. This invention relies solely on a camera to measure light intensity distribution. The absence of an additional wavefront sensor and beam splitter allows for a more compact system layout and easier debugging.

[0219] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A dynamic beam shaping method based on a spatial light modulator, comprising a dynamic beam shaping system based on a spatial light modulator, consisting of a laser source module, a beam expanding and collimating optical path, a beam modulator, a focusing optical component, a detection and feedback module, and a control and optimization module connected in sequence, characterized in that: The laser source module is located at the input surface of the system's optical path and is used to output a laser beam; The beam expanding and collimating optical path is used to expand and collimate the laser beam. The beam modulator is used to modulate the laser beam after it has been shaped and collimated. The focusing optical component is used to form a flat-top beam of the target shape on the focal plane or the target surface after the modulation of the laser beam. The detection and feedback module is located on the output surface of the system's optical path and is used to acquire the light spot intensity distribution image of the flat-top beam and transmit it to the control and optimization module. The control and optimization module is used to update the phase distribution of the beam modulator in real time according to the beam intensity distribution image, so as to correct the output flat-top beam. The method includes the following steps: Generate a target intensity distribution image based on the parameters of the target flat-top light spot and the output camera pixel size. The ideal beam-shaping phase of the target flat-top spot is calculated based on the input and output plane coordinate systems. And construct the optimized phase of the beam modulator. The expression; The light spot intensity distribution image acquired by the output surface detection and feedback module Image of target intensity distribution Construct a composite performance index expression; In each stage of the multi-stage SPGD algorithm, two sets of random perturbation vectors are generated and the ideal beam shaping phase is updated respectively. The resulting two sets of new beam shaping phases are loaded onto the beam modulator and two new laser beams are generated. Acquire the spot intensity distribution images corresponding to the two new laser beams, and calculate their composite performance index according to the composite performance index expression. Optimize the target parameters in the optimized phase expression, and the beam modulator dynamically shapes the laser beams according to the optimized target parameters to output a uniform flat-top beam from the focusing optical component. The ideal beam-shaping phase of the target flat-top light spot is calculated based on the input and output plane coordinate systems. And construct the optimized phase of the beam modulator. The expression specifically includes the following steps: The shape, size, and energy distribution range of the target flat-top light spot are defined, and a corresponding target intensity image is generated by combining this with the pixel size of the output camera. ; Ideal beam shaping phase for achieving the target flat-top light spot is pre-calculated based on Fourier optics and the standing phase method. ,in The physical coordinate system of the input plane; Calculate Zernike polynomials And construct the optimized phase of the beam modulator. expression: ; ; ; ; ; ; Where p is the sorting index of the Zernike polynomial excluding constant terms (n=0, m=0). The polar coordinate system is the input plane of the system's optical path. Here, represents the orthogonal normalization coefficients, n is the radial order of the Zernike polynomial, and m is the azimuth frequency of the Zernike polynomial. Let be a radial polynomial of Zernike polynomials, with order (n, m). Let be the objective parameters of the Zernike polynomial to be optimized. Let P be a vector composed of Zernike polynomials; The light spot intensity distribution image acquired by the output surface detection and feedback module Image of target intensity distribution Constructing a composite performance index expression specifically includes the following steps: Calculate the light spot intensity distribution image Image of target intensity distribution Normalized mean square error terms between The calculation formula is: ; in, Let F be the first balancing parameter used to balance the order-of-magnitude difference between the two images caused by the measurement, and let F be the Frobenius norm. The inner product of two image matrices; Calculate the structural similarity between the light spot intensity distribution image and the target intensity distribution image. The calculation formula is: ; in, Image of light spot intensity distribution Image of target intensity distribution covariance, Image of light spot intensity distribution standard deviation For target intensity distribution image standard deviation This is a tiny constant used to avoid calculation errors caused by dividing the function by zero; According to the normalized mean square error term With structural similarity Construct a composite performance index expression: ; in, This is a second balancing parameter used to balance the order-of-magnitude difference between the normalized mean square error term and the structural similarity; The process of generating two sets of randomly oriented perturbation vectors and updating the phase of the ideal beam shaping vectors accordingly includes the following steps: Generate a set of perturbation vectors with random positive and negative values ​​but uniform size. : ; in, It is a length of Vectors that satisfy the Rademacher distribution Let be the amplitude of the disturbance in stage s; According to a set of disturbance vectors and the target parameter vector of the current stage The ideal beam shaping phase is updated to obtain new beam shaping phases with two opposite directions under perturbation. : ; in, Let be the number of Zernike polynomials in the s-th stage. Let p be the target parameter with sorting index p and iteration number t in stage s. Let p be the perturbation with sorting index p and iteration number t in stage s. Let p be the p-th Zernikal polynomial.

2. The dynamic beam shaping method based on a spatial light modulator according to claim 1, characterized in that: The optimization of the target parameter in the optimized phase expression specifically includes the following steps: Numerical gradient estimated using composite performance index of the intensity distribution images of the two new laser beams. The calculation formula is: ; in, It is a composite performance index of the spot intensity distribution image corresponding to the new beam shaping phase under positive perturbation vector perturbation. It is a composite performance index of the spot intensity distribution image corresponding to the new beam shaping phase under negative perturbation vector perturbation. Let be the disturbance amplitude in stage s. Direction vector The p-th value; Based on the numerical gradient Composition of numerical gradient vector The Adam optimizer is used to optimize the target parameters in the optimized phase expression.

3. The dynamic beam shaping method based on a spatial light modulator according to claim 2, characterized in that: The numerical gradient vector The optimization of the target parameters in the optimized phase expression using the Adam optimizer specifically includes the following steps: Set the number of preheating iterations When the current iteration number is less than or equal to the preheating iteration number Then, the objective parameters in the optimized phase expression are updated using the gradient descent method, and the calculation formula is as follows: ; in, Let be the target parameter vector at the (t+1)th iteration in the S-th stage. Let be the target parameter vector at the t-th iteration. The learning rate for stage S. Let be the numerical gradient vector at the t-th iteration; When the current iteration number is greater than the preheating iteration number At that time, the first-order momentum at the t-th iteration of the Adam optimizer is used. With second momentum The target parameters in the optimized phase expression are updated using the following formula: ; in, These are the decay coefficients of the Adam optimizer. It is a constant. , These are the first-order and second-order momentum of the normalized Adam optimizer, respectively; Read the most recent The composite performance index is calculated based on the number of iterations, and the change in the number of iterations is fitted using least squares. The slope of the fitted curve is then calculated. and in The ratio of the initial slope to the initial slope is less than a preset proportional threshold. Then proceed to the next stage, repeating the above steps until all stages are completed and the final target parameters are obtained.

4. The dynamic beam shaping method based on a spatial light modulator according to claim 3, characterized in that: Before each stage begins, the following steps are also included: Regarding the disturbance amplitude mentioned in the previous stage and learning rate Adaptive calibration is performed to obtain the optimal perturbation amplitude for each stage. and learning rate .

5. The dynamic beam shaping method based on a spatial light modulator according to claim 4, characterized in that: The disturbance amplitude of the previous stage and learning rate Adaptive calibration includes the following steps: Repeated measurements under the same phase diagram Subcomposite performance index measurement value , where ε is zero-mean noise, and the variance of the noise is estimated. ; Based on the measured values ​​of the composite performance index The Local Linearity Index (LLI) is calculated using the following formula: ; ; ; ; in, These are the m-th random direction vectors. The measured values ​​of the composite performance index corresponding to the positive and negative candidate perturbation amplitudes. These are the initial target parameters used in the current optimization phase. The mean of the composite performance index measurements obtained from the corresponding M0 measurements, for The standard deviation of measurement For the median operation, This refers to the local linearity calculated using the m-th composite performance index measurement value under the current disturbance amplitude; Based on the local linearity index Determine the optimal disturbance amplitude and based on Adaptively set the learning rate.

6. The dynamic beam shaping method based on a spatial light modulator according to claim 5, characterized in that: The local linearity index Determine the optimal perturbation amplitude, and based on Adaptive learning rate setting specifically includes the following steps: In the first stage of the multi-stage SPGD algorithm, the initial value of the perturbation amplitude is... And in subsequent stages, the initial value of the disturbance amplitude will be... Set to the disturbance amplitude of the previous stage ; Search for the local linearity index (LLI) under the current perturbation amplitude and check if it exceeds the threshold. If the threshold is exceeded, the Local Linearity Index (LLI) under the current perturbation amplitude is reduced to half of its original value; conversely, if the threshold is exceeded, the LLI is increased to twice its original value, until a set of values ​​at the threshold values ​​is obtained. Up and down and closest to the threshold Disturbance amplitude ; Based on the disturbance amplitude The optimal disturbance amplitude was calculated. The calculation formula is: ; based on and optimal perturbation amplitude Adaptive learning rate setting The calculation formula is: ; in, For use in controlling the learning rate The corresponding target parameter update magnitude is a normal number within the expected range.

Citation Information

Patent Citations

  • Digital holographic microscopy system phase distortion correction method based on adaptive optics

    CN120163745A

  • SAR photoelectric processor phase error compensation method based on Zernike polynomial

    CN121114938A