Method and device for computing light field transmission in ultra-high precision modeling of inter-satellite laser link

Abstract

The application discloses an optical field transmission calculation method and device in inter-satellite laser link ultra-high precision modeling, mainly solves the problems that in the inter-satellite laser link ultra-long distance optical field propagation modeling of space gravitational wave detection, the calculation result precision is difficult to quantitatively confirm, the integral precision parameter depends on artificial experience setting and the multiple working condition high precision repeated calculation overhead is large. The method takes the target equivalent optical path error threshold as a unified precision criterion, realizes the high precision transmission calculation of the exit pupil optical field of the transmitting end telescope to the entrance pupil area of the receiving end telescope through numerical integral node and weight coefficient calling, exit pupil sampling node mapping, complex amplitude contribution task mapping between exit pupil sampling nodes and grid points, complex amplitude contribution calculation, complex accumulation and complex calculation verification under different integral precision levels. The application can improve the reliability, repeated calculation efficiency and engineering deployment capability of the calculation result in the inter-satellite laser link optical field transmission modeling.

Description

Technical Field

[0001] This invention relates to a method for calculating optical field transmission, specifically to a method and apparatus for calculating optical field transmission in ultra-high precision modeling of inter-satellite laser links. Background Technology

[0002] Space-based gravitational wave detection missions typically employ a constellation of multiple spacecraft spanning hundreds of thousands to millions of kilometers, establishing inter-satellite laser links between adjacent spacecraft. The light field in these inter-satellite laser links, after propagating over extremely long distances, reaches the entrance pupil of the receiving telescope. Its complex amplitude, phase, and energy distributions affect the phase extraction accuracy of subsequent interferometric measurements, thus impacting the reliability of the gravitational wave signal calculation results. Therefore, during the design, error allocation, assembly verification, and performance evaluation of space-based gravitational wave telescopes, it is essential to perform high-precision modeling of the light field distribution at the entrance pupil of the receiving telescope after its extremely long propagation.

[0003] In existing technologies, far-field optical propagation can typically be calculated using Fresnel-Fraunhofer diffraction integrals, fast Fourier transform (FFT) methods, or analytical approximation methods. The FFT method offers high computational efficiency and is suitable for wave optics simulations at conventional scales and with conventional accuracy requirements. However, in applications requiring ultra-long transmission distances and extremely high phase measurement accuracy, such as space gravitational wave detection, this method is easily limited by uniform sampling, periodic boundary assumptions, spectral aliasing, and increased grid size, making it difficult to achieve a stable balance between computational accuracy and computational resource requirements. Analytical approximation methods usually rely on regular apertures, ideal beams, or simple wavefront structures. Their applicability is limited for actual exit pupil optical fields that contain non-ideal amplitude distributions, complex wavefront distortions, central obstruction, irregular pupils, assembly residuals, or platform pointing errors.

[0004] Numerical integration methods can directly calculate the optical field distribution in the entrance pupil region of the receiving telescope based on the complex amplitude distribution at the exit pupil of the transmitting telescope, exhibiting good physical applicability. Existing research has been able to model the far-field propagation process of the complex amplitude at the exit pupil of the transmitting telescope to the entrance pupil region of the receiving telescope based on Fresnel-Fronhoeffer diffraction integrals or other numerical integration methods, and can be used to analyze the influence of different wavefront errors, beam distributions, and platform pointing errors on the optical field distribution at the receiving end.

[0005] However, such methods typically focus on physical propagation mechanism analysis and offline simulation calculations. Their computational flow generally involves pre-determining the integration order, the number of sampling nodes at the exit pupil, the grid density at the receiver, and the data precision level, then directly outputting the receiver's amplitude or phase distribution. Existing procedures often lack unified accuracy criteria and self-verification mechanisms to determine whether the current computational parameters can meet the picometer-level or sub-picometer-level equivalent optical path error requirements. Researchers often need to repeatedly adjust the integration order, the number of sampling nodes, the grid density, or the data precision level based on experience, and then compare the results through multiple trial calculations to determine whether the calculation results meet the target accuracy requirements.

[0006] Especially when the inter-satellite transmission distance, receiver observation area, platform pointing state, or beam parameters change, the original calculation parameters may no longer meet the same equivalent optical path error requirements, necessitating parameter reselection and accuracy verification. This leads to problems in existing numerical integration simulation processes, such as numerous manual trial calculations, low efficiency of parameter reconfiguration, and indirect confirmation of calculation result accuracy in scenarios involving multiple distances, multiple operating conditions, and high-precision repetitive calculations. Consequently, these issues make it difficult to meet the requirements of high-precision engineering modeling for inter-satellite laser links in space gravitational wave detection, which demands high computational reliability and efficiency for repetitive calculations.

[0007] Current numerical integration simulation workflows typically treat integration nodes, integration weights, exit pupil sampling nodes, and receiver grid points as internal computational variables. They lack a readily accessible coefficient repository for integration nodes and weights corresponding to different integration orders, accuracy levels, and exit pupil integration intervals. Furthermore, they lack a mechanism for retrieving integration accuracy levels corresponding to the target's equivalent optical path error threshold. Therefore, when multiple inter-satellite optical field transmission calculations are required under varying computational accuracies, optical system parameters, or exit pupil conditions, problems such as repetitive generation of integration nodes and weights, cumbersome parameter configuration, insufficient consistency in computational workflows, and unstable computational efficiency can easily arise.

[0008] Furthermore, existing inter-satellite light field propagation simulation workflows typically treat complex amplitude calculation as a single mathematical integration process, lacking a computational resource organization structure geared towards engineering deployment. For example, existing workflows generally do not decompose the complex amplitude contribution relationship between the exit pupil sampling nodes and the receiving grid points into schedulable computational tasks, and also lack a device architecture for modularizing and parallelizing the complex amplitude contribution calculation, complex accumulation, and precision verification processes. Since inter-satellite light field propagation calculations involve a large number of complex amplitude contribution calculations between the exit pupil sampling nodes and the receiving grid points, the lack of a clear task decomposition, task scheduling, and accumulation structure is detrimental to engineering deployment on multi-core CPUs, GPUs, or dedicated computing hardware.

[0009] Therefore, there is an urgent need for a high-precision optical field transmission calculation method and device for inter-satellite laser links for space gravitational wave detection. This device should not only be able to complete the optical field calculation at the receiver based on the complex amplitude distribution of the transmitter's exit pupil, but also automatically select or adjust the integration accuracy level according to the target equivalent optical path error threshold, quickly call the corresponding integration nodes and weight coefficients, establish the complex amplitude contribution task mapping relationship between the exit pupil sampling nodes and the receiver grid points, and achieve self-verification by comparing the phase difference of the complex amplitude calculation results under different integration accuracy levels. This would improve the computational reliability, computational efficiency, and engineering feasibility of optical field transmission modeling for inter-satellite laser links. Summary of the Invention

[0010] The purpose of this invention is to solve the problems of difficulty in quantifying and confirming the accuracy of calculation results, dependence on manual experience in setting integral accuracy parameters, and high overhead of repeated high-precision calculations under multiple operating conditions in the ultra-long distance optical field propagation modeling of inter-satellite laser links for space gravitational wave detection. The invention provides a method and apparatus for optical field transmission calculation in ultra-high precision modeling of inter-satellite laser links.

[0011] To achieve the above objectives, the technical solution provided by this invention is as follows:

[0012] A method for calculating the optical field transmission in ultra-high precision modeling of inter-satellite laser links, characterized by the following steps:

[0013] Step 1: Obtain the calculation input parameters and integration accuracy level of the inter-satellite laser link;

[0014] Step 2: Determine the radius of the observation area at the receiving end based on the calculated input parameters. And discretize the observation area into grid points ; and All are positive integers. , ;

[0015] Step 3: Based on the calculation input parameters, retrieve the coefficients from the pre-established coefficient repository that match the current integration accuracy level. The corresponding numerical integration nodes and their weight coefficients;

[0016] Step 4: Establish a polar coordinate system on the exit pupil of the transmitter. , , ; within the standard range Radial numerical integration nodes are selected above. and its weighting coefficients Angular numerical integration nodes and its weighting coefficients And radial numerical integration nodes and angular numerical integration nodes Mapping to polar coordinates Determine the discrete sampling nodes on the exit pupil Generate a set of discrete sampling nodes for the exit pupil; , , and These represent the number of radial and angular integration nodes, respectively;

[0017] Step 5: Based on the exit pupil function of the transmitter. The validity of the discrete sampling node set of the exit pupil is determined to obtain the valid exit pupil sampling node set;

[0018] Step 6: Establish effective exit pupil sampling nodes and grid points The mapping relationship between nodes generates a node-grid point complex amplitude contribution calculation task to calculate the contribution of each effective exit pupil sampling node. For all grid points Complex amplitude contribution ;

[0019] Step 7: For each grid point All corresponding complex amplitude contributions The complex amplitude is calculated by performing weighted summation on each number separately. ;

[0020] Step 8: Adopt the integral accuracy level The preset key grid points are recalculated. This yields complex amplitude recalculation results with a higher integration accuracy level. and according to and Calculate the equivalent optical path error Meanwhile, the largest one is denoted as ;

[0021] Step 9, when Preset target equivalent optical path error threshold If so, increase the integration order, increase the integration accuracy level, or increase the data accuracy in the calculation of complex amplitude contribution and complex weighted accumulation process, and return to step 3; when Preset target equivalent optical path error threshold When the time is right, the current complex amplitude calculation result will be output. As the final result of the optical field transmission calculation, the optical field transmission calculation is completed.

[0022] Furthermore, in step 5, when the discrete sampling node When a node is located in a central occlusion region, a non-transparent region, or an invalid pupil region, the node should be removed or its corresponding two-dimensional integral weight should be removed. Set to zero; when discrete sampling nodes If the node is located within the effective exit pupil area of ​​the transmitting telescope, then the node is retained. For effective exit pupil sampling nodes... Its two-dimensional integral weight satisfy:

[0023] ;

[0024] In the formula: The exit pupil radius of the transmitting end; These are the radial coordinates of the exit pupil after mapping;

[0025] In step 6, the mapping relationship satisfies:

[0026] ;

[0027] ;

[0028] In the formula: These are the mapped exit pupil angular coordinates;

[0029] The contribution of complex amplitude satisfy:

[0030] ;

[0031] In the formula: The exit pupil of the transmitter is at the sampling node. Complex amplitude at the location; For sampling nodes from the exit pupil propagation to grid points The diffraction propagation function;

[0032] In step 7, the complex amplitude calculation result The expression is as follows:

[0033] ;

[0034] In step 8, first according to and Calculate the phase difference between the two Then, based on this phase difference The equivalent optical path error was obtained through conversion. :

[0035] ;

[0036] ;

[0037] In the formula: To map the phase difference to Phase wrapping within the interval; This refers to the equivalent optical path error of the same grid point under different integration accuracy levels.

[0038] Furthermore, in step 7, the complex weighted accumulation includes calculating the contributions of the complex amplitudes respectively. real part and the virtual part The weighted summation satisfies:

[0039] ;

[0040] .

[0041] Furthermore, when Preset target equivalent optical path error threshold When this is the case, prioritize increasing the integration order or integration accuracy level; if increasing the integration order or integration accuracy level still does not meet the target equivalent optical path error threshold requirement, then increase the data accuracy in the complex amplitude contribution calculation process or increase the data accuracy in the complex weighted accumulation process.

[0042] Furthermore, in step 1, the calculation input parameters include at least the operating wavelength. Inter-satellite transmission distance Exit pupil radius of the transmitting telescope Receiver telescope entrance pupil radius Maximum pointing error of the platform Target equivalent optical path error threshold Complex amplitude distribution at the exit pupil of the transmitting telescope ;

[0043] in, ; ;

[0044] In the formula: For the exit pupil amplitude distribution, The exit pupil pupil function, This is the exit pupil wavefront error. For wave number, The operating wavelength;

[0045] In step 2, based on the entrance pupil radius of the receiving telescope... Inter-satellite transmission distance and the platform's maximum pointing error Determine the radius of the observation area at the receiving end. :

[0046] ;

[0047] In the formula: This is for engineering sampling margin;

[0048] Step 3 specifically involves, based on the target equivalent optical path error threshold... Retrieves coefficients from a pre-established coefficient repository that match the current calculation accuracy level. The corresponding numerical integration nodes and their weight coefficients.

[0049] Furthermore, in step 3, the coefficient repository stores numerical integration nodes and their weight coefficients corresponding to different integration orders, different computational accuracy levels, and different exit pupil integration intervals in advance.

[0050] The coefficient repository adopts a hierarchical index structure, which includes at least an integration method type index, an integration order index, a target equivalent optical path error level index, an exit pupil radius adaptation index, and a data precision level index.

[0051] Furthermore, in step 3, the numerical integration nodes and their weight coefficients are pre-generated based on a high-precision numerical integration method, and after variable transformation and adaptation for the integration interval of the telescope exit pupil, they are stored in the coefficient storage library in the form of a floating-point array.

[0052] The high-precision numerical integration method includes one or more of the Gauss-Legend integration method, the Gauss-Chebyshev integration method, or the adaptive Simpson integration method.

[0053] The key receiver grid points include one or more of the following: the center grid point of the receiver observation area, the edge grid point of the entrance pupil, the grid point in the neighborhood of the maximum amplitude of the light spot, the grid point in the neighborhood of the centroid of the light spot, and the grid point whose phase gradient is greater than a preset gradient threshold.

[0054] In step 8, the key grid points include one or more of the following: the center grid point of the observation area, the edge grid point of the entrance pupil, the neighborhood grid point of the maximum amplitude of the light spot, the neighborhood grid point of the centroid of the light spot, and the grid point whose phase gradient is greater than a preset gradient threshold.

[0055] In step 6, the node-grid point complex amplitude contribution calculation task includes the effective exit pupil sampling node number, grid point number, and two-dimensional integral weight. Transmission distance Operating wavelength Exit pupil polar coordinate parameters and and grid point coordinate parameters .

[0056] Meanwhile, this invention also provides an optical field transmission calculation device for ultra-high precision modeling of inter-satellite laser links, used to implement the aforementioned optical field transmission calculation method, which is characterized by:

[0057] It includes a coefficient storage library, a computing resource configuration module, a node scheduler, a parallel computing coordinator, a diffraction calculation unit, an accumulator array, and a precision verification module;

[0058] The coefficient repository stores the calculation input parameters and integration accuracy levels;

[0059] The computing resource configuration module is used to determine the observation area radius of the receiving end based on the computing input parameters. And discretize the observation area into grid points Retrieves the coefficient from the coefficient repository with the current integration accuracy level. The corresponding numerical integration nodes and their weight coefficients;

[0060] The node scheduler is used in standard intervals. Radial numerical integration nodes are selected above. and its weighting coefficients Angular numerical integration nodes and its weighting coefficients And radial numerical integration nodes and angular numerical integration nodes Mapping to polar coordinates Determine the discrete sampling nodes on the exit pupil Generate a set of discrete sampling nodes for the exit pupil;

[0061] The parallel computing coordinator is used to establish effective exit pupil sampling nodes and grid points. The mapping relationship between nodes and grid points is used to generate complex amplitude contribution calculation tasks for nodes and grid points, and at the same time, complex amplitude contribution calculation tasks are assigned to diffraction calculation units.

[0062] The diffraction calculation unit is used to perform the complex amplitude contribution calculation task to calculate each effective exit pupil sampling node. For all grid points Complex amplitude contribution ;

[0063] The accumulator array is used for each grid point All corresponding complex amplitude contributions The complex amplitude is calculated by performing weighted summation on each number separately. ;

[0064] The accuracy verification module is used to employ integral accuracy levels. Recalculate the preset key grid points. This yields complex amplitude recalculation results with a higher integration accuracy level. and according to and Calculate the equivalent optical path error Meanwhile, the largest one is denoted as ;when Preset target equivalent optical path error threshold If necessary, increase the integration order, improve the integration accuracy level, or improve the data accuracy in the calculation of complex amplitude contribution and complex weighted accumulation process, and recalculate the optical field transmission; when Preset target equivalent optical path error threshold When the time is right, the current complex amplitude calculation result will be output. As the final result of the optical field transmission calculation, the optical field transmission calculation is completed.

[0065] Furthermore, the accumulator array includes real accumulators and imaginary accumulators;

[0066] The real part accumulator is used to contribute to the complex amplitude. The real parts are weighted and accumulated, and the imaginary part accumulator is used to calculate the contribution of the complex amplitude. The imaginary parts are weighted and accumulated. Compared with the prior art, the present invention has the following beneficial technical effects:

[0067] First, the present invention uses the target equivalent optical path error threshold. As a unified accuracy criterion, the selection of integration accuracy level, self-verification of integration order, and judgment of calculation result output are all unified under the same error standard. By comparing the phase difference between the current integration accuracy level and the complex amplitude calculation results at the receiver at a higher integration accuracy level, and converting this phase difference into an equivalent optical path error, it is possible to quantitatively determine whether the current calculation result meets the target accuracy requirements. This avoids relying solely on manual experience to judge whether the integration order or the number of sampling nodes is sufficient, thereby improving the reliability of picometer-level or sub-picometer-level phase calculation results.

[0068] Secondly, this invention performs recalculation verification on key receiver grid points at a higher integration accuracy level, ensuring the verifiability of the calculation results while avoiding blind high-order recalculation on all receiver grid points. These key receiver grid points may include the center point of the receiver entrance pupil, entrance pupil edge points, spot peaks or centroid neighborhood points, and grid points with large phase gradients. Therefore, it can improve the accuracy verification capability of calculation results in phase error-sensitive areas while controlling computational overhead.

[0069] Furthermore, this invention provides a data organization foundation for the aforementioned accuracy self-verification and feedback adjustment process through an integral node and weight coefficient library, node mapping, complex amplitude contribution task mapping, and complex number accumulation structure. This foundation can be quickly invoked, scheduled, and implemented in parallel, which is beneficial for improving the computational efficiency, stability, and engineering feasibility when performing repeated calculations over multiple transmission distances, multiple optical parameters, and multiple operating conditions. Attached Figure Description

[0070] Figure 1 This is a schematic diagram of the overall process of the optical field transmission calculation method in the ultra-high precision modeling of inter-satellite laser links described in this invention.

[0071] Figure 2 This is a schematic diagram of the self-verification and feedback adjustment process of equivalent optical path error based on different integration accuracy levels as described in this invention.

[0072] Figure 3 This is a schematic diagram of the far-field receiver optical field distribution calculated by the method described in this invention, wherein (a) is the receiver optical field amplitude distribution diagram, (b) is the receiver optical field phase distribution diagram, and (c) is the feedback convergence result of the equivalent optical path error of key grid points under different integration accuracy levels. Detailed Implementation

[0073] To make the objectives, advantages, and features of the present invention clearer, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments. Those skilled in the art should understand that these embodiments are merely used to explain the technical principles of the present invention and are not intended to limit the scope of protection of the present invention.

[0074] This embodiment presents a method for calculating the optical field transmission in ultra-high precision modeling of inter-satellite laser links. It is geared towards far-field optical field propagation modeling of inter-satellite laser links for space gravitational wave detection, and is particularly suitable for high-precision calculation scenarios that require obtaining complex amplitude distribution, phase distribution, intensity distribution, or energy distribution at the entrance pupil of the receiver.

[0075] The core of this embodiment lies in: using the target equivalent optical path error threshold. As a unified accuracy criterion, it combines coefficient storage, node scheduling, node-grid point task mapping, complex amplitude contribution calculation, complex number accumulation, and integral accuracy self-verification to realize inter-satellite optical field transmission calculation for engineering deployment.

[0076] like Figure 1 As shown, the optical field transmission calculation method in ultra-high precision modeling of inter-satellite laser links in this embodiment specifically includes the following steps:

[0077] Step 1: Calculate and obtain input parameters

[0078] Obtain the computational input parameters and integration accuracy level for the inter-satellite laser link. The computational input parameters must include at least the operating wavelength. Inter-satellite transmission distance Exit pupil radius of the transmitting telescope Receiver telescope entrance pupil radius Maximum pointing error of the platform Target equivalent optical path error threshold and the complex amplitude distribution at the exit pupil of the transmitting telescope .

[0079] In one implementation, when the target error is at a target phase error threshold When representing, the target phase error threshold can be used. Converted to target equivalent optical path error threshold .

[0080] Complex amplitude distribution at the exit pupil of the transmitting telescope It can be determined jointly by the exit pupil amplitude distribution, the exit pupil pupil function, the transmittance distribution, and the wavefront error. In one embodiment, the complex amplitude distribution satisfies:

[0081] ;

[0082] In the formula, For the exit pupil amplitude distribution, The exit pupil pupil function, This is the exit pupil wavefront error. This refers to the wave number. It should be noted that... It can also come from measured pupil complex amplitude data, complex field data exported from optical design software, interferometric reconstruction results, or complex amplitude data constructed from Zernike aberration coefficients.

[0083] Step 2: Determining the observation area and discretizing the grid at the receiving end.

[0084] Based on the entrance pupil radius of the receiving telescope Inter-satellite transmission distance and the platform's maximum pointing error Determine the radius of the observation area at the receiving end. .

[0085] In one implementation, to avoid the impact of boundary truncation effect or insufficient edge sampling on subsequent energy calculations, spot centroid calculations, phase distribution analysis, and light field calculations within the receiver's entrance pupil region, a method can be used... Set a sampling margin outside the defined basic observation area radius. The radius of the observation area at the receiving end satisfies: . Not used as the target equivalent optical path error threshold The direct conversion result, similarly, the target equivalent optical path error threshold. It is not used directly to determine the radius of the observation area at the receiving end, but is used for subsequent selection of integration accuracy level and self-verification of calculation results.

[0086] After determining the observation area of ​​the receiving end, the observation area is discretized into... There are 10 receiver grid points. Each receiver grid point is denoted as . ,in and All are positive integers. , The grid points at the receiving end can be represented in rectangular coordinates, or they can be converted to polar coordinates as needed.

[0087] In engineering implementation, the observation area at the receiving end can cover the entrance pupil area of ​​the receiving telescope and the spot offset range caused by the platform's maximum pointing error. For scenarios requiring further calculation of the bounding energy, spot centroid, or edge interpolation results, an additional sampling margin can be set in the specific implementation; however, this sampling margin is not used as the target equivalent optical path error threshold. A direct function of .

[0088] Step 3: Calling the coefficient repository and integral node

[0089] like Figure 2 As shown, the optical field transmission computing device in this embodiment includes a coefficient storage library, a computing resource configuration module, a node scheduler, a parallel computing coordinator, a diffraction computing unit, an accumulator array, and a precision verification module.

[0090] The coefficient repository is used to pre-store numerical integration nodes and their weight coefficients corresponding to different integration orders, accuracy levels, and exit pupil integration intervals. The numerical integration nodes and their weight coefficients can be pre-generated based on the Gauss-Legendé integration method, the Gauss-Chebyshev integration method, the adaptive Simpson integration method, or other high-precision numerical integration methods. The Gauss-Legendé integration method is preferred for generating standard intervals. The integral nodes and weight coefficients on the surface.

[0091] To avoid repeatedly generating integration nodes and weighting coefficients during each optical field transmission calculation, this embodiment pre-generates integration nodes and weighting coefficients corresponding to different integration orders and different precision levels, and stores them in a coefficient repository as floating-point arrays. The coefficient repository adopts a hierarchical index structure, which may include an integration method type index, an integration order index, a target equivalent optical path error level index, an exit pupil radius adaptation index, and a data precision level index.

[0092] The computational resource allocation module calculates the target equivalent optical path error threshold. Retrieve the current integration accuracy level from the coefficient repository. The corresponding radial numerical integration nodes, angular numerical integration nodes, and their weighting coefficients. Among them, Indicates the current integration accuracy level.

[0093] Step 4: Mapping from standard interval nodes to the polar coordinate region of the transmitter's exit pupil

[0094] Establish a polar coordinate system at the exit pupil of the telescope at the transmitting end. ,in , Within the standard range Radial numerical integration nodes are selected above. and its weighting coefficients and angular numerical integration nodes and its weighting coefficients .in, and These represent the number of radial and angular integration nodes, respectively. , .

[0095] The node scheduler performs node mapping, mapping the radial and angular numerical integration nodes under the standard interval to the polar coordinate region of the transmitter telescope's exit pupil, thus obtaining the discrete sampling nodes of the exit pupil. The mapping relationship satisfies: , In the formula, These are the radial coordinates of the exit pupil after mapping. The coordinates of the exit pupil angle after mapping.

[0096] Consider the area element in polar coordinate integrals Mapped two-dimensional integral weights satisfy: In the formula, The two-dimensional integral weights are mapped to the polar coordinate region of the telescope's exit pupil. Radial numerical integration nodes The corresponding weighting coefficients, Angular numerical integration nodes The corresponding weighting coefficients, The exit pupil radius of the transmitting telescope. These are the radial coordinates of the exit pupil after mapping.

[0097] Step 5: Determination of valid exit pupil sampling nodes

[0098] After completing node mapping, the node scheduler determines the output pupil function of the transmitting telescope. The validity of the discrete sampling node set at the exit pupil is determined to obtain the valid exit pupil sampling node set. When the discrete sampling node... When the node is located within the effective exit pupil region of the transmitting telescope, retain the node and generate the corresponding complex amplitude contribution calculation task; when the node is discretely sampled... If a node is located in a central occlusion region, a non-transparent region, a partially occluded region, or another invalid pupil region, the node should be removed, or the corresponding two-dimensional integral weight should be removed. Set to zero.

[0099] Based on the above validity determination, this embodiment can accommodate central obstruction, irregular pupils, locally invalid pupil regions, non-uniform transmittance distribution, and actual exit pupil boundary variations introduced by assembly and adjustment errors. In engineering implementation, the actual exit pupil may deviate from the ideal circular aperture due to factors such as optical assembly and adjustment errors, central obstruction, support structure blockage, or measured aperture boundary deviation. Therefore, aperture boundary data, transmittance distribution data, or pupil function data obtained from ground testing or optical design software can be used as input for determining the effective exit pupil sampling node.

[0100] Step 6: Node-Grid Point Complex Amplitude Contribution Task Mapping

[0101] The parallel computing coordinator establishes a task mapping relationship between valid exit pupil sampling nodes and receiving grid points. For each valid exit pupil sampling node... and each receiving grid point Generate a complex amplitude contribution calculation task .

[0102] Complex amplitude contribution calculation task At least the following should be included: valid exit pupil sampling node number, receiver grid point number, and two-dimensional integral weight. Transmission distance Operating wavelength Exit pupil polar coordinate parameters and and the grid point coordinate parameters at the receiving end .

[0103] The diffraction calculation unit calculates the complex amplitude contribution corresponding to each complex amplitude contribution calculation task based on the diffraction propagation function. In one implementation, the first... The effective exit pupil sampling node for the first Complex amplitude contribution of each receiving grid point satisfy: In the formula, Indicates the first The effective exit pupil sampling node for the first The complex amplitude contribution of each grid point at the receiving end; The exit pupil of the transmitting telescope at the sampling node Complex amplitude at the location; For sampling nodes from the exit pupil propagation to receiving grid points The diffraction propagation function; This represents the two-dimensional integral weight corresponding to the exit pupil sampling node.

[0104] Among them, the diffraction propagation function It can be determined based on the Fresnel diffraction form, Fraunhofer diffraction form, or other equivalent far-field propagation model used, and may include propagation phase factor, propagation distance normalization factor, coordinate transformation factor, or other factors corresponding to the selected propagation model.

[0105] Step 7: Complex summation of complex amplitude contributions

[0106] Accumulator array for the same receiver grid point Corresponding multiple complex amplitude contributions By performing complex summation, the current integration accuracy level can be obtained. The complex amplitude calculation results at the grid point of the receiving end are as follows The summation of complex numbers satisfies: .

[0107] If the two-dimensional integral weight corresponding to the invalid exit pupil sampling node If the value is already zero, the above summation can be performed on all discrete sampling nodes of the exit pupil; if invalid exit pupil sampling nodes have been removed, the above summation can be performed only on valid exit pupil sampling nodes.

[0108] To ensure numerical consistency in the coherent superposition of complex amplitudes, the accumulator array can include real-part accumulators and imaginary-part accumulators. The real-part accumulator is used to accumulate the real parts of each complex amplitude contribution, and the imaginary-part accumulator is used to accumulate the imaginary parts of each complex amplitude contribution, satisfying the following: , .

[0109] After the accumulation is complete, it can be based on and This yields the complex amplitude, phase distribution, intensity distribution, or energy distribution at the grid points of the receiver. For example, it can be obtained as follows: , , In the formula, This represents the complex vibration amplitude at the current integration accuracy level. The phase at the current integration accuracy level, This represents the intensity at the current integration accuracy level.

[0110] Step 8: Self-verification of integration accuracy

[0111] After obtaining the complex amplitude calculation results for each grid point within the observation area of ​​the receiving end, the accuracy verification module executes the integration accuracy self-verification process. For example... Figure 3 As shown, the accuracy verification module calls a higher-level integral accuracy function than the current one. Second integral accuracy level By recalculating the grid points at the key receiver end, we obtain complex amplitude recalculation results with a higher integration accuracy level. .

[0112] The key receiver grid points constitute a set In one implementation, the set It may include one or more of the following: the center grid point of the observation area at the receiving end, the edge grid point of the entrance pupil, the grid point in the neighborhood of the maximum amplitude of the light spot, the grid point in the neighborhood of the centroid of the light spot, and the grid point whose phase gradient is greater than a preset gradient threshold.

[0113] For sets At the same key receiver grid point, the accuracy verification module compares the current integration accuracy level. Calculation results of complex amplitude under the following conditions Compared with the second integral accuracy level Complex amplitude recalculation results phase difference between The phase difference, after phase wrapping, satisfies:

[0114] ;

[0115] in, This indicates that the phase difference is mapped to... Phase wrapping within the interval.

[0116] The accuracy verification module will check the phase difference. Converted to equivalent optical path error It satisfies: .

[0117] Step 9: Feedback and Adjustment

[0118] When the key receiver grid point set The equivalent optical path error in the equation satisfies: Then determine the current integration accuracy level. The target equivalent optical path error threshold requirement is met, and the complex amplitude calculation result at the current integration accuracy level is output. This serves as the result of the optical field transmission calculation at the entrance pupil of the receiving end.

[0119] When the key receiver grid point set The equivalent optical path error in the equation satisfies: Then determine the current integration accuracy level. The target equivalent optical path error threshold requirement was not met. At this point, the accuracy verification module sends feedback to the computational resource configuration module indicating insufficient computational accuracy. Based on this feedback, the computational resource configuration module increases the integration order, improves the integration accuracy level, or increases the data accuracy in the complex amplitude contribution calculation and complex number accumulation process. It then re-calls the corresponding integration nodes and weight coefficients, repeating the node mapping, validity determination, task mapping, complex amplitude contribution calculation, complex number accumulation, and accuracy verification processes until the target equivalent optical path error threshold requirement is met.

[0120] It should be noted that the target equivalent optical path error threshold It is used for selecting the integration accuracy level, self-checking the integration order, and determining the output of calculation results; it is not directly used to determine the radius of the observation area at the receiving end.

[0121] Example of specific parameters

[0122] In one specific embodiment, the operating wavelength for The laser power is 2W, the beam waist radius at the exit point is 150mm, and the exit pupil diameter of the transmitting telescope is [missing information]. Inter-satellite transmission distance The maximum pointing error of the platform is 2.5 million kilometers. The value is on the order of nanoarcs, set to 20 nard, and further adjusted to the order of 0.5 pm based on the required computational precision. The computational resource configuration module determines the radius of the observation area at the receiver based on the above parameters. and receiving end grid points and retrieve the current integration accuracy level from the coefficient repository. The corresponding integration nodes and weight coefficients.

[0123] In this embodiment, the complex amplitude distribution of the exit pupil of the transmitting telescope can be determined by the Gaussian beam, the telescope wavefront error, and the pupil function. The wavefront error can be given by the Zernike aberration coefficient, measured wavefront data, or wavefront data exported from optical design software. The pupil function can represent an ideal circular pupil, a centrally occluded pupil, an irregular pupil, or the actual pupil boundary reconstructed from test data.

[0124] As a specific example, the wavefront error distribution may include a Y-direction tilt term, a field curvature term, an astigmatism term, and a coma term, with corresponding wavefront error coefficients of 2.35 nm, 2 nm, 0.6 nm, and 1.5 nm, respectively. These values ​​are provided only as a specific embodiment to illustrate that the present invention can receive wavefront error data from a real optical system and construct the complex amplitude distribution of the transmitter's exit pupil.

[0125] Get as Figure 3 The diagram shows the light field distribution at the entrance pupil of the far-field receiver. Figure 3 (a) is a diagram showing the amplitude distribution of the light field at the entrance pupil of the receiver; Figure 3 (b) is a phase distribution diagram of the light field at the entrance pupil of the receiver. Figure 3 (c) shows the feedback convergence results of the equivalent optical path error of key grid points under different integration accuracy levels.

[0126] In this embodiment, several representative grid points are selected from the receiver's observation grid to form a key grid point set. The key grid points include: one entrance pupil center point, used to characterize the baseline state of the optical field energy and phase distribution at the receiver; four edge points located above, below, left, and right of the entrance pupil, used to check the integration error and phase change in the edge region of the entrance pupil; one point with the maximum optical field amplitude, used to characterize the complex amplitude calculation result at the location of the strongest optical field energy; one point near the centroid of the light spot, used to characterize the optical field state at the actual center of the light spot; twelve grid points with the largest phase gradients, used to focus on checking regions where phase changes are rapid and numerical integration is prone to non-convergence; and four neighboring grid points located above, below, left, and right of the amplitude peak point, used to check whether the local complex amplitude changes near the optical field peak are accurately calculated.

[0127] According to the above selection rules, theoretically 23 candidate key grid points can be obtained. Since candidate points of different categories may be located at the same grid position, for example, the center point of the entrance pupil may be both the point of maximum amplitude and the centroid of the spot, the program performs deduplication on the candidate grid points, and the final number of key grid points can be less than 23.

[0128] The results in the figure, after multiple rounds of integration accuracy feedback adjustments, show that the number of key grid points is 22 (with 2 overlapping), while the maximum equivalent optical path error decreases from the initial 1359.7 pm to 4.301 pm, and finally to 0.0016833 pm. The results demonstrate that this invention can adjust the integration accuracy level based on the computational differences between different integration accuracy levels, enabling the receiver's optical field calculation results to gradually and automatically converge towards the required high accuracy. Especially in the final round of calculation, the maximum equivalent optical path error between the current integration configuration and a higher integration configuration has been reduced to 1.6833 × 10⁻⁶. −15 m, thus verifying the effectiveness of the self-verification and feedback adjustment method for integral accuracy based on the target equivalent optical path error.

[0129] Parallel computing and device implementation

[0130] In one implementation, the parallel computing coordinator can allocate complex amplitude contribution calculation tasks to multi-core CPU threads, GPU computing units, FPGAs, or dedicated computing hardware based on the hardware platform configuration. For the same receiver grid point, the complex amplitude contributions corresponding to multiple exit pupil sampling nodes can be calculated in parallel and the complex numbers are accumulated in the accumulator array; for different receiver grid points, parallel scheduling can also be performed according to grid point number, calculation task number, or data block number.

[0131] Under the condition of meeting the target equivalent optical path error threshold requirement, the diffraction calculation unit can adopt a double-precision floating-point number, extended-precision floating-point number, or a mixed-precision calculation method that meets the accuracy requirements after recalculation and verification. The accumulator array can adopt double-precision accumulation, compensated summation, block accumulation, or tree-structured reduction accumulation to improve the numerical stability in the large-scale complex amplitude contribution accumulation process.

[0132] Through the above method, this invention organizes the calculation of complex amplitude contributions from numerous exit pupil sampling nodes and receiver grid points in inter-satellite optical field transmission calculations into a schedulable, accumulative, and recalculated engineering calculation process. This process aligns with the target equivalent optical path error threshold. The combination of precision verification and feedback adjustment mechanisms, under the target equivalent optical path error threshold requirement, is beneficial to improving the reliability of calculation results, efficiency of repeated calculations, and engineering feasibility of inter-satellite laser link optical field transmission modeling.

[0133] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and are not intended to limit them. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features therein. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the present invention.

Claims

1. A method for calculating optical field transmission in ultra-high precision modeling of inter-satellite laser links, characterized in that, Includes the following steps: Step 1: Obtain the calculation input parameters and integration accuracy level of the inter-satellite laser link; Step 2: Determine the radius of the observation area at the receiving end based on the calculated input parameters. And discretize the observation area into grid points ; and All are positive integers. , ; Step 3: Based on the calculation input parameters, retrieve the coefficients from the pre-established coefficient repository that match the current integration accuracy level. The corresponding numerical integration nodes and their weight coefficients; Step 4: Establish a polar coordinate system on the exit pupil of the transmitter. , , ; within the standard range Radial numerical integration nodes are selected above. and its weighting coefficients Angular numerical integration nodes and its weighting coefficients And radial numerical integration nodes and angular numerical integration nodes Mapping to polar coordinates Determine the discrete sampling nodes on the exit pupil Generate a set of discrete sampling nodes for the exit pupil; , , and These represent the number of radial and angular integration nodes, respectively; Step 5: Based on the exit pupil function of the transmitter. The validity of the discrete sampling node set of the exit pupil is determined to obtain the valid exit pupil sampling node set; Step 6: Establish effective exit pupil sampling nodes and grid points The mapping relationship between nodes generates a node-grid point complex amplitude contribution calculation task to calculate the contribution of each effective exit pupil sampling node. For all grid points Complex amplitude contribution ; Step 7: For each grid point All corresponding complex amplitude contributions The complex amplitude is calculated by performing weighted summation on each number separately. ; Step 8: Adopt the integral accuracy level Recalculate the preset key grid points. This yields complex amplitude recalculation results with a higher integration accuracy level. and according to and Calculate the equivalent optical path error Meanwhile, the largest one is denoted as ; Step 9, when Preset target equivalent optical path error threshold If so, increase the integration order, increase the integration accuracy level, or increase the data accuracy in the calculation of complex amplitude contribution and complex weighted accumulation process, and return to step 3; when Preset target equivalent optical path error threshold When the time is right, the current complex amplitude calculation result will be output. As the final result of the optical field transmission calculation, the optical field transmission calculation is completed.

2. The optical field transmission calculation method in ultra-high precision modeling of inter-satellite laser links according to claim 1, characterized in that, In step 5, when the discrete sampling node When a node is located in a central occlusion region, a non-transparent region, or an invalid pupil region, the node should be removed or its corresponding two-dimensional integral weight should be removed. Set to zero; when discrete sampling nodes If the node is located within the effective exit pupil area of ​​the transmitting telescope, then the node is retained. For effective exit pupil sampling nodes... Its two-dimensional integral weight satisfy: ； In the formula: The exit pupil radius of the transmitting end; These are the radial coordinates of the exit pupil after mapping; In step 6, the mapping relationship satisfies: ； ； In the formula: These are the mapped exit pupil angular coordinates; The contribution of complex amplitude satisfy: ； In the formula: The exit pupil of the transmitter is at the sampling node. Complex amplitude at the location; For sampling nodes from the exit pupil propagation to grid points The diffraction propagation function; In step 7, the complex amplitude calculation result The expression is as follows: ； In step 8, first according to and Calculate the phase difference between the two Then, based on this phase difference The equivalent optical path error was obtained through conversion. : ； ； In the formula: To map the phase difference to Phase wrapping within the interval; This refers to the equivalent optical path error of the same grid point under different integration accuracy levels.

3. The optical field transmission calculation method in ultra-high precision modeling of inter-satellite laser links according to claim 1 or 2, characterized in that: In step 7, the complex weighted accumulation includes the contributions of each complex amplitude. real part and the virtual part The weighted summation satisfies: ； 。 4. The optical field transmission calculation method in ultra-high precision modeling of inter-satellite laser links according to claim 3, characterized in that: when Preset target equivalent optical path error threshold When this is the case, prioritize increasing the integration order or integration accuracy level; if increasing the integration order or integration accuracy level still does not meet the target equivalent optical path error threshold requirement, then increase the data accuracy in the complex amplitude contribution calculation process or increase the data accuracy in the complex weighted accumulation process.

5. The optical field transmission calculation method in ultra-high precision modeling of inter-satellite laser links according to claim 1, characterized in that: In step 1, the calculation input parameters include at least the operating wavelength. Inter-satellite transmission distance Exit pupil radius of the transmitting telescope Receiver telescope entrance pupil radius Maximum pointing error of the platform Target equivalent optical path error threshold Complex amplitude distribution at the exit pupil of the transmitting telescope ; in, ; ; In the formula: For the exit pupil amplitude distribution, The exit pupil pupil function, This is the exit pupil wavefront error. For wave number, The operating wavelength; In step 2, based on the entrance pupil radius of the receiving telescope... Inter-satellite transmission distance and the platform's maximum pointing error Determine the radius of the observation area at the receiving end. : ； In the formula: This is for engineering sampling margin; Step 3 specifically involves, based on the target equivalent optical path error threshold... Retrieves coefficients from a pre-established coefficient repository that match the current calculation accuracy level. The corresponding numerical integration nodes and their weight coefficients.

6. The optical field transmission calculation method in ultra-high precision modeling of inter-satellite laser links according to claim 5, characterized in that: In step 3, the coefficient repository stores numerical integration nodes and their weight coefficients corresponding to different integration orders, different computational accuracy levels, and different exit pupil integration intervals in advance. The coefficient repository adopts a hierarchical index structure, which includes at least an integration method type index, an integration order index, a target equivalent optical path error level index, an exit pupil radius adaptation index, and a data precision level index.

7. The optical field transmission calculation method in ultra-high precision modeling of inter-satellite laser links according to claim 6, characterized in that: In step 3, the numerical integration nodes and their weight coefficients are pre-generated based on a high-precision numerical integration method, and after variable transformation and adaptation for the integration interval of the telescope exit pupil, they are stored in the coefficient storage library in the form of floating-point arrays. The high-precision numerical integration method includes one or more of the Gauss-Legend integration method, the Gauss-Chebyshev integration method, or the adaptive Simpson integration method. The key receiver grid points include one or more of the following: the center grid point of the receiver observation area, the edge grid point of the entrance pupil, the grid point in the neighborhood of the maximum amplitude of the light spot, the grid point in the neighborhood of the centroid of the light spot, and the grid point whose phase gradient is greater than a preset gradient threshold. In step 8, the key grid points include one or more of the following: the center grid point of the observation area, the edge grid point of the entrance pupil, the neighborhood grid point of the maximum amplitude of the light spot, the neighborhood grid point of the centroid of the light spot, and the grid point whose phase gradient is greater than a preset gradient threshold. In step 6, the node-grid point complex amplitude contribution calculation task includes the effective exit pupil sampling node number, grid point number, and two-dimensional integral weight. Transmission distance Operating wavelength Exit pupil polar coordinate parameters and and grid point coordinate parameters .

8. A light field transmission calculation device for ultra-high precision modeling of inter-satellite laser links, used to implement the light field transmission calculation method according to any one of claims 1 to 7, characterized in that: It includes a coefficient storage library, a computing resource configuration module, a node scheduler, a parallel computing coordinator, a diffraction calculation unit, an accumulator array, and a precision verification module; The coefficient repository stores the calculation input parameters and integration accuracy levels; The computing resource configuration module is used to determine the observation area radius of the receiving end based on the computing input parameters. And discretize the observation area into grid points Retrieves the coefficient from the coefficient repository with the current integration accuracy level. The corresponding numerical integration nodes and their weight coefficients; The node scheduler is used in standard intervals. Radial numerical integration nodes are selected above. and its weighting coefficients Angular numerical integration nodes and its weighting coefficients And radial numerical integration nodes and angular numerical integration nodes Mapping to polar coordinates Determine the discrete sampling nodes on the exit pupil Generate a set of discrete sampling nodes for the exit pupil; The parallel computing coordinator is used to establish effective exit pupil sampling nodes and grid points. The mapping relationship between nodes and grid points is used to generate complex amplitude contribution calculation tasks for nodes and grid points, and at the same time, complex amplitude contribution calculation tasks are assigned to diffraction calculation units. The diffraction calculation unit is used to perform the complex amplitude contribution calculation task to calculate each effective exit pupil sampling node. For all grid points Complex amplitude contribution ; The accumulator array is used for each grid point All corresponding complex amplitude contributions The complex amplitude is calculated by performing weighted summation on each number separately. ; The accuracy verification module is used to employ integral accuracy levels. Recalculate the preset key grid points. This yields complex amplitude recalculation results with a higher integration accuracy level. and according to and Calculate the equivalent optical path error Meanwhile, the largest one is denoted as ;when Preset target equivalent optical path error threshold If necessary, increase the integration order, improve the integration accuracy level, or improve the data accuracy in the calculation of complex amplitude contribution and complex weighted accumulation process, and recalculate the optical field transmission; when Preset target equivalent optical path error threshold When the time is right, the current complex amplitude calculation result will be output. As the final result of the optical field transmission calculation, the optical field transmission calculation is completed.

9. The optical field transmission calculation device for ultra-high precision modeling of inter-satellite laser links according to claim 8, characterized in that: The accumulator array includes real accumulators and imaginary accumulators; The real part accumulator is used to contribute to the complex amplitude. The real parts are weighted and accumulated, and the imaginary part accumulator is used to calculate the contribution of the complex amplitude. The imaginary part is weighted and accumulated.

Images