Non-sequence optical system stray light propagation analysis method

Abstract

The present application belongs to the technical field of optical simulation, and particularly relates to a stray light propagation analysis method for a non-sequential optical system. The method comprises the following steps: S1: establishing a three-dimensional geometric model of the optical system to be analyzed; S2: establishing a directed propagation edge of an effective propagation module pair, and obtaining a non-sequential module propagation graph; S3: based on the directed propagation edge of the non-sequential module propagation graph, extracting paths to be expanded, abstracting each path to be expanded into a path state to construct an initial flood fill queue, and performing priority sorting on the initial flood fill queue to obtain a flood fill queue; S4: performing dynamic pruning according to the contribution upper bound of each path to be expanded to obtain an effective path set; and S5: based on the stray light energy of each effective path contained in the effective path set, calculating the total stray light of the detector contained in the optical system to be analyzed. The present application aims to realize accurate solution of stray light full path, intelligent pruning, and effective identification of key paths and sensitive modules for a complex system.

Description

Technical Field

[0001] This invention belongs to the field of optical simulation technology, and in particular relates to a method for analyzing stray light propagation in non-sequential optical systems. Background Technology

[0002] Stray light is undesirable radiated energy in an optical system that deviates from the main imaging and transmission channels and ultimately reaches the detection surface, image plane, or other evaluation surfaces. It mainly originates from residual reflections from optical surfaces, micro-scattering, scattering from the lens barrel and mechanical structure, and off-axis orders of diffraction elements. The presence of stray light can cause problems such as background lifting, ghost images, decreased imaging contrast and signal-to-noise ratio, and increased measurement deviation, severely affecting the imaging quality and detection accuracy of optical equipment.

[0003] Currently, stray light analysis in engineering mainly includes two mainstream technical approaches: The first is the Monte Carlo random sampling method based on non-sequential ray tracing. This method is highly versatile and adaptable to optical systems with complex geometries. It detects the surface energy distribution by emitting a massive number of rays and is currently the mainstream analysis method. However, this method is highly dependent on the number of samples. For low-probability critical stray light paths formed by multiple weak reflections, scattering, and inefficient diffraction, the sampling hit rate is extremely low, easily leading to missed detection of critical interference paths and difficulty in convergence. The second approach is a module-level solution method based on radiative transfer and surface element coupling. This method has a clear physical mechanism and can accurately analyze the directional energy coupling relationship between components. However, this method cannot effectively handle complex propagation paths with combined properties of transmission, reflection, scattering, and diffraction. Furthermore, the number of paths in complex systems tends to expand exponentially, lacking effective pruning strategies and resulting in extremely high computational costs.

[0004] The two existing technologies still have many common defects: the simulation results can only output the energy distribution of the detection surface and cannot trace the stray light propagation path, source module and contribution ratio of each segment; there is a lack of dynamic pruning mechanism based on physical contribution, and invalid branches consume a lot of computing power; at the same time, there is a lack of a unified framework to handle the propagation path coupled with multiple physical mechanisms, and most of them are modeled independently of a single mechanism. Summary of the Invention

[0005] In view of this, the present invention aims to provide a stray light propagation analysis method for non-sequential optical systems, addressing the technical problems of existing stray light analysis techniques, such as easy omission of low-probability critical paths, difficulty in efficiently solving complex non-sequential paths, inability to perform path source contribution analysis, lack of dynamic pruning mechanisms, and inability to uniformly adapt to multi-physics coupling propagation mechanisms. The present invention proposes a stray light propagation analysis method for non-sequential optical systems, aiming to achieve accurate full-path solving of stray light in complex systems, intelligent pruning, and effective identification of critical paths and sensitive modules.

[0006] To achieve the above objectives, the technical solution created by this invention is implemented as follows: A method for analyzing stray light propagation in a non-sequential optical system, specifically including the following steps: S1: Establish a three-dimensional geometric model of the optical system to be analyzed, and discretize the three-dimensional geometric model into no less than three propagation modules, and configure the attributes of each propagation module; S2: Perform pairwise matching verification on all propagation modules and establish directed propagation edges for valid propagation module pairs to obtain the non-sequence module propagation graph; S3: Based on the directed propagation edges of the non-sequential module propagation graph, extract the paths to be expanded, abstract each path to be expanded into a path state to construct an initial flood filling queue, and sort the initial flood filling queue by priority to obtain the flood filling queue. S4: Calculate the total contribution of each path to be expanded contained in the flood filling queue, and perform dynamic pruning based on the upper bound of the contribution of each path to be expanded to obtain a set of effective paths. S5: Calculate the total stray light of the detectors in the optical system to be analyzed based on the stray light energy of each effective path contained in the effective path set.

[0007] Furthermore, in step S1, the types of propagation modules include at least an optical working surface module, a detector module, and a light source module.

[0008] Furthermore, in step S1, the optical properties of the optical working surface module include at least one or more of the following parameters: reflectivity, transmittance, absorptivity, total integrated scattering, two-way scattering distribution function, two-way transmission distribution function, diffraction order efficiency, and medium transmission factor. The boundary properties of the light source module include at least one or more of the following parameters: initial radiant flux, angular radiation distribution characteristics, wavelength or band, polarization state, and effective radiation aperture size of the light source; The boundary properties of the detector module include at least one or more of the following parameters: detector surface position, effective receiving size, pixel region discretization method, receiving direction constraint, and statistical rules.

[0009] Furthermore, step S2 specifically includes the following steps: S21: Calculate the single-segment propagation parameters of the current propagation module pair: ; ； Where i and j are the index numbers of the two propagation modules contained in the current propagation module. For the single-segment propagation parameters of the current propagation module pair, The visibility factor for the current propagation module pair. The medium transmittance factor for the current propagation module pair. Let i be the solid angle corresponding to propagation module i at propagation module j. The projection geometry factor of the current propagation module pair. These are optional modifications for the current propagation module pair. Let j be the equivalent area of ​​the propagation module. Let be the angle between the direction of the line connecting the center of propagation module i to the center of propagation module j and the normal direction of propagation module j. The center distance of the current propagation module pair; S22: Determine whether the single-segment propagation parameter of the current propagation module pair is zero. If not, then take the current propagation module pair as a valid propagation module pair and establish a directed propagation edge for the current valid propagation module pair. Otherwise, execute step S23. S23: Replace the current propagation module pair with the next propagation module pair, repeat steps S21-S22, obtain the directed propagation edges of all valid propagation module pairs, and construct a non-sequence module propagation graph based on all valid propagation module pairs and their directed propagation edges.

[0010] Furthermore, step S3 specifically includes the following steps: S31: Traverse all directed propagation edges, take the propagation module corresponding to the light source module as the starting point of the path, and extend the initial path along the directed propagation edge as the initial path to be extended, and obtain the initial path to be extended set. The initial path to be extended includes at least: the initial path to be extended containing only the propagation module corresponding to the light source module, and the initial path to be extended consisting of the propagation module corresponding to the light source module and the downstream propagation module connected to the propagation module corresponding to the light source module. S32: Each path to be expanded is abstracted into a path state using the following formula, where the path state corresponding to the Kth path to be expanded is... : ; in, Let K be the set of propagation modules contained in the Kth path to be expanded; This refers to the end propagation module that has been extended to the Kth path to be extended; This represents the direction domain or direction bucket for the Kth path to be expanded. Let K be the channel label of the Kth path to be expanded; Let K be the path layer number of the Kth path to be expanded; Let be the minimum number of remaining propagation layers from the end propagation module to the detection module in the Kth path to be extended; The cumulative contribution to the Kth path to be expanded; This is the maximum upper bound estimate of the state propagation for the Kth path to be expanded; Let $\overall_bound$ be the state bound corresponding to the $K$-th path to be expanded, and $\overall_bound$ be the state bound. ; This is the historical encoding of the Kth path to be expanded; and These are the wavelength label and polarization label for the Kth path to be extended, respectively. S33: Starting from the propagation module corresponding to the detector module, perform a reverse breadth-first traversal of the propagation graph of the non-sequence module and record the minimum number of remaining propagation layers from each propagation module to the detector module. S34: Add the paths to be expanded to the initial flood fill queue in ascending order of propagation layer number; S35: Based on the path status corresponding to each path to be expanded, calculate the priority score of each path to be expanded included in the initial flood filling queue using the following formula, where the priority score of the Kth path to be expanded is... for: ; S36: Based on the calculation results of step S35, sort and adjust the paths to be expanded in the initial flood filling queue according to the priority score from large to small, and arrange the paths to be expanded with the same priority score or the priority score difference is less than the preset tolerance in the order of the number of propagation layers from small to large to obtain the flood filling queue.

[0011] Furthermore, step S4 specifically includes the following steps: S41: Select a path to be expanded from the flood fill queue and calculate the individual contribution of each directed propagation edge contained in the current path to be expanded. S42: Calculate the total contribution of the current path to be expanded based on the individual contributions of each directed propagation edge: ; in, The total contribution of the Pth path to be expanded. Let n be the light source flux, and n be the total number of propagation modules. For the k-th propagation module, For the surface interaction response of the k-th propagation module, For the single-segment propagation coupling term of the k-th propagation module; S43: Replace the current path to be expanded with the next path to be expanded, and repeat steps S41-S42 until the total contribution of all paths to be expanded is obtained.

[0012] Furthermore, in step S41, when the directed propagation edge is from the propagation module... To the propagation module Extending this, the formula for calculating the individual contribution of the current directed propagation edge is as follows: The current contribution of the directed propagation edge to the mirror reflection channel is: ; in, This represents the current interactive response of the directed propagation edge under the mirror reflection channel. For propagation module i at wavelength λ and incident angle The reflectivity function under the condition of polarization state pol; The current contribution of the directed propagation edge to the transmission channel is: ; in, This represents the interactive response of the current directed propagation edge in the transmission channel. For propagation module i at wavelength λ and incident angle Transmittance function under polarization state pol; Set up a propagation module The incident flux is The total scattering flux of the current directed propagation edge is: ; in, The total scattered flux is obtained by converting the incident flux into the propagation module i. The total integral scattering coefficient of propagation module i represents the proportion of the incident flux allocated to the scattering channel; Propagation Module The corner domain The scattering response is normalized as follows: ; in, For the scattered light from propagation module i to enter the angular region where propagation module j is located. The normalized scattering response, The total integrated scattering coefficient of propagation module i, For the scattered light rays from propagation module i, the full-space output solid angle is given. For propagation module i in the incident direction and the two-way scattering distribution function at the emission direction ω, The incident direction, The direction of launch; The current contribution of the directed propagation edge in the discrete diffraction order sub-channels is: ; in, The interactive response corresponds to the discrete diffraction order m. For propagation module i at wavelength λ and incident angle The m-th order diffraction efficiency function under the condition of polarization state pol.

[0013] Furthermore, step S4 specifically includes: S41: Estimate the upper bound of the contribution of the Kth path to be expanded using the following formula. : ; ; in, Let K be the current cumulative contribution value of the Kth path to be expanded. To estimate the upper bound of the maximum product for the Kth path to be expanded, Let t be the minimum number of remaining propagation layers from the end propagation module of the Kth path to be extended to the corresponding propagation module of the detector, where t is the index of the remaining propagation layer. For the first Maximum single-layer multiplier upper bound estimate within each remaining propagation layer; S42: If the upper bound of the contribution of the Kth path to be expanded is less than the first preset threshold. If the current path to be expanded is not found, the expansion of the current path to be expanded will be terminated, and the Kth path to be expanded will be considered an invalid path; otherwise, the Kth path to be expanded will be considered a valid path. S43: Replace the k-th path to be expanded with the (k+1)-th path to be expanded, and repeat steps S41-S42 until all valid paths are obtained; S44: Determine if there are duplicate valid paths among all valid paths. If so, merge the duplicate valid paths and update the valid paths. Otherwise, proceed directly to step S45. S45: Pair up each propagation module in the current valid path, count the number of back propagations between the paired propagation modules, and delete valid paths with a number of back propagations greater than the second preset threshold. S46: Replace the current valid path with the next valid path, repeat step S45, complete the further update of valid paths, and obtain the final set of valid paths.

[0014] Furthermore, in step S42, the formula for calculating the first preset threshold is: ; in, The upper bound of the contribution to the flood filling queue is the maximum value. α is a preset proportionality constant, 0 < α < 1.

[0015] Furthermore, in step S5, the formula used to calculate the total stray light of the detectors contained in the optical system to be analyzed is as follows: ; in, This represents the total stray light accumulated by the detector at pixel (x,y). Let P be the pixel index, and P be the valid path to the pixel (x, y). This is the set of valid paths that hit the pixel (x, y).

[0016] Compared with the prior art, the present invention can achieve the following beneficial effects: (1) The stray light propagation analysis method for non-sequential optical systems described in this invention can systematically discover key stray light paths with low probability and high impact. The flood-filled path expansion is not a random sample hit, but rather starts from the light source module and advances layer by layer on the non-sequential module propagation map. The path expansion priority is controlled by a combined detector-guided scoring mechanism. Thus, it can find direct and short paths first, and can also continue to track key paths formed by multiple weak reflections, multiple scatterings and inefficient diffraction when needed.

[0017] (2) The stray light propagation analysis method for non-sequential optical systems described in this invention can provide a quantitative interpretation of each path to be extended. The single-segment propagation and module interaction in the path to be extended are both calculated explicitly with radiative transfer, thus the path-level contribution ranking and module-level contribution ranking can be given directly.

[0018] (3) The stray light propagation analysis method for non-sequential optical systems described in this invention can uniformly handle surfaces with mixed properties of transmission, specular reflection, scattering, and diffraction. Composite property surfaces can generate multi-channel sub-states under the same parent state, enabling different physical mechanisms to be consistently expressed within the same framework.

[0019] (4) The non-sequential optical system stray light propagation analysis method described in this invention can perform dynamic pruning and priority control based on the current contribution and the remaining upper bound, thereby improving the efficiency of engineering calculation and avoiding wasting resources on a large number of low contribution branches.

[0020] (5) The non-sequential optical system stray light propagation analysis method described in this invention can directly serve engineering rectification, that is, by outputting the main contributing modules and critical paths, it can directly guide the extinction of the lens barrel, the adjustment of the position of the light-blocking ring, the blackening of the edge, the improvement of the coating and the design of the structure to avoid light. Attached Figure Description

[0021] The accompanying drawings, which form part of this invention, are used to provide a further understanding of the invention. The illustrative embodiments and descriptions of the invention are used to explain the invention and do not constitute an undue limitation of the invention. In the drawings: Figure 1 A schematic flowchart illustrating the stray light propagation analysis method for a non-sequential optical system as described in an embodiment of the present invention; Figure 2 The non-sequential propagation graph between modules described in the embodiments of the present invention; Figure 3 This is a schematic diagram of the Cooke three-lens imaging system with a lens barrel, as described in an embodiment of the present invention.

[0022] Explanation of reference numerals in the attached figures: 1. Off-axis stray light source; 2. Lens tube structure; 3. First positive lens; 4. Second negative lens; 5. Third positive lens; 6. Image plane detector. Detailed Implementation

[0023] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and specific embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and do not constitute a limitation thereof.

[0024] It should be noted that, unless otherwise specified, the embodiments and features described in the present invention can be combined with each other.

[0025] In the description of this invention, it should be understood that the terms "center," "longitudinal," "lateral," "upper," "lower," "front," "rear," "left," "right," "vertical," "horizontal," "top," "bottom," "inner," and "outer," etc., indicating orientations or positional relationships based on the orientations or positional relationships shown in the accompanying drawings, are only for the convenience of describing this invention and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation, and therefore should not be construed as a limitation on this invention. Furthermore, the terms "first," "second," etc., are used for descriptive purposes only and should not be construed as indicating or implying relative importance or implicitly specifying the number of indicated technical features. Thus, features defined with "first," "second," etc., may explicitly or implicitly include one or more of that feature. In the description of this invention, unless otherwise stated, "a plurality of" means two or more.

[0026] In the description of this invention, it should be noted that, unless otherwise explicitly specified and limited, the terms "installation," "connection," and "linking" should be interpreted broadly. For example, they can refer to a fixed connection, a detachable connection, or an integral connection; they can refer to a mechanical connection or an electrical connection; they can refer to a direct connection or an indirect connection through an intermediate medium; and they can refer to the internal connection of two components. Those skilled in the art will understand the specific meaning of the above terms in this invention based on the specific circumstances.

[0027] The invention will now be described in detail with reference to the accompanying drawings and embodiments.

[0028] like Figures 1-2As shown, this invention proposes a method for analyzing stray light propagation in a non-sequential optical system, specifically including the following steps: S1: Establish a three-dimensional geometric model of the optical system to be analyzed, and discretize the three-dimensional geometric model into no less than three propagation modules, and configure the attributes of each propagation module; S2: Perform pairwise matching verification on all propagation modules and establish directed propagation edges for valid propagation module pairs to obtain the non-sequence module propagation graph; S3: Based on the directed propagation edges of the non-sequential module propagation graph, extract the paths to be expanded, abstract each path to be expanded into a path state to construct an initial flood filling queue, and sort the initial flood filling queue by priority to obtain the flood filling queue. S4: Calculate the total contribution of each path to be expanded contained in the flood filling queue, and perform dynamic pruning based on the upper bound of the contribution of each path to be expanded to obtain a set of effective paths. S5: Calculate the total stray light of the detectors in the optical system to be analyzed based on the stray light energy of each effective path contained in the effective path set.

[0029] It should be noted that this invention relates to stray light analysis technology for optical systems, and is applicable to imaging lenses, infrared detection systems, space remote sensing optical systems, projection systems, lithography illumination and projection systems, laser transmission systems, AR / VR optical systems, and other optomechanical systems that require analysis of complex non-sequential stray light propagation.

[0030] Furthermore, the present invention aims to propose a stray light propagation analysis method for non-sequential optical systems that combines the flood-fill algorithm with radiative transfer theory. This method discretizes the complex optical system into propagation modules at the geometric modeling level, expands reachable paths layer by layer using the flood-fill algorithm at the path management level, and recursively calculates the single-segment propagation and total contribution of each path using the radiative transfer formula at the physical calculation level. It also utilizes the upper bound of the contribution to achieve priority expansion, pruning, and state merging, thereby obtaining the stray light distribution on the detection surface, the critical path ranking, and the sensitive modules (referring to...) or The identification results of propagation modules that are particularly high and significantly contribute to the stray light flux reaching the image plane are presented. This invention does not rely solely on random ray hit paths, nor does it only estimate local radiative exchange between static modules. Instead, it establishes a unified solution framework of "module discretization - propagation map construction - path state definition - flood filling expansion - radiative transfer recursion - pruning and merging - detector surface accumulation".

[0031] In some embodiments, in step S1, the types of propagation modules include at least an optical working surface module, a detector module, and a light source module.

[0032] It should be noted that, firstly, a three-dimensional geometric model of the optical system to be analyzed is established, including at least the spatial relationships of the optical surface, mechanical structural components, inner wall of the lens barrel, light-blocking ring, aperture, detector, reference surface, evaluation surface, and light source. Then, according to surface type, normal variation, roughness variation, functional zoning, and engineering requirements, the system is discretized into multiple propagation modules. A propagation module can be a whole surface or a local element.

[0033] To ensure the representativeness of parameters within a propagation module, it is preferable to maintain approximately consistent normal variations, surface shape characteristics, and roughness levels within the same module. Each propagation module should record at least the following: propagation module number, geometric center, normal, area or equivalent area, surface category, applicable wavelength range, and occlusion information with other propagation modules.

[0034] In some embodiments, in step S1, the optical properties of the optical working surface module include at least one or more of the following parameters: reflectivity, transmittance, absorptivity, total integrated scattering, bidirectional scattering distribution function, bidirectional transmission distribution function, diffraction order efficiency, and medium transmission factor. The boundary properties of the light source module include at least one or more of the following parameters: initial radiant flux, angular radiation distribution characteristics, wavelength or band, polarization state, and effective radiation aperture size of the light source; The boundary properties of the detector module include at least one or more of the following parameters: detector surface position, effective receiving size, pixel region discretization method, receiving direction constraint, and statistical rules.

[0035] It should be noted that each propagation module (specifically the optical surface module) must be configured with at least one or more of the following parameters: reflectivity transmittance Absorption rate Total integral scattering Two-way scattering distribution function Two-way transmission distribution function Diffraction order efficiency Medium permeability factor Wavelength tags and polarization response tags. For surfaces with composite properties, the same propagation module can simultaneously have transmission, specular reflection, scattering, and diffraction channels.

[0036] For the light source module, its initial radiant flux, emission angle distribution, wavelength or band, polarization state, and source size also need to be defined. For the detector module, the detector surface position, size, pixel division method, receiving direction constraints, and statistical rules need to be defined.

[0037] In some embodiments, step S2 specifically includes the following steps: S21: Calculate the single-segment propagation parameters of the current propagation module pair: ; ； Where i and j are the index numbers of the two propagation modules contained in the current propagation module. For the single-segment propagation parameters of the current propagation module pair, The visibility factor for the current propagation module pair. The medium transmittance factor for the current propagation module pair. Let i be the solid angle corresponding to propagation module i at propagation module j. The projection geometry factor of the current propagation module pair. These are optional modifications for the current propagation module pair. Let j be the equivalent area of ​​the propagation module. Let be the angle between the direction of the line connecting the center of propagation module i to the center of propagation module j and the normal direction of propagation module j. The center distance of the current propagation module pair; S22: Determine whether the single-segment propagation parameter of the current propagation module pair is zero. If not, then take the current propagation module pair as a valid propagation module pair and establish a directed propagation edge for the current valid propagation module pair. Otherwise, execute step S23. S23: Replace the current propagation module pair with the next propagation module pair, repeat steps S21-S22, obtain the directed propagation edges of all valid propagation module pairs, and construct a non-sequence module propagation graph based on all valid propagation module pairs and their directed propagation edges.

[0038] It should be noted that for any two propagation modules and If the following conditions are met: geometric visibility, path unobstructed or quantifiable partial occlusion, directional coverage, and upper bound of propagation contribution greater than a threshold, then a propagation graph from the non-sequence module is established. arrive The directed propagation edge. Preferably, the propagation graph adopts a "two-layer representation" method: the first layer edge only represents the geometric reachability relationship between propagation modules and saves basic parameters such as visibility, effective solid angle, projection factor, and medium transmission factor; the second layer dynamically expands the same geometric edge into sub-channels such as transmission, specular reflection, scattering, and diffraction according to the incident direction, wavelength, and polarization of the current state when the path is expanded.

[0039] The single-segment propagation parameters between storage modules on a finite propagation edge are as follows. The single-segment propagation coupling term is defined as: ; in, The visibility factor for propagation pairs; The medium permeability factor for propagation pairs; For the propagation module In the propagation module The effective solid angle corresponding to the location; The projection geometry factor for the propagation pair; Optional corrections for propagation pairs are used to represent aperture trimming, occlusion correction, polarization correction, or engineering empirical corrections. In engineering implementation, these can be... Adopted The normalized effective solid angle is used as a reference, or incorporated into Unified processing, thereby It is used as a dimensionless coupled multiplier in the path recursion process.

[0040] Under the surface element approximation condition, if the propagation module The equivalent area is The angle between the normal of the propagation module and the direction of the connecting line is The center distance between the two propagation modules is Then the effective solid angle can be approximately expressed as: ; The above nodes and propagation edges constitute the non-sequential module propagation graph. The set of nodes Represents the set of propagation modules and the set of edges. This represents a directed propagation edge that satisfies the propagation conditions.

[0041] In some embodiments, step S3 specifically includes the following steps: S31: Traverse all directed propagation edges, take the propagation module corresponding to the light source module as the starting point of the path, and extend the initial path along the directed propagation edge as the initial path to be extended, and obtain the initial path to be extended set. The initial path to be extended includes at least: the initial path to be extended containing only the propagation module corresponding to the light source module, and the initial path to be extended consisting of the propagation module corresponding to the light source module and the downstream propagation module connected to the propagation module corresponding to the light source module. S32: Each path to be expanded is abstracted into a path state using the following formula, where the path state corresponding to the Kth path to be expanded is... : ; in, Let K be the set of propagation modules contained in the Kth path to be expanded; This refers to the end propagation module that has been extended to the Kth path to be extended; This represents the direction domain or direction bucket for the Kth path to be expanded. Let K be the channel label of the Kth path to be expanded; Let K be the path layer number of the Kth path to be expanded; Let be the minimum number of remaining propagation layers from the end propagation module to the detection module in the Kth path to be extended; The cumulative contribution to the Kth path to be expanded; This is the maximum upper bound estimate of the state propagation for the Kth path to be expanded; Let $\overall_bound$ be the state bound corresponding to the $K$-th path to be expanded, and $\overall_bound$ be the state bound. ; This is the historical encoding of the Kth path to be expanded; and These are the wavelength label and polarization label for the Kth path to be extended, respectively. S33: Starting from the propagation module corresponding to the detector module, perform a reverse breadth-first traversal of the propagation graph of the non-sequence module and record the minimum number of remaining propagation layers from each propagation module to the detector module. S34: Add the paths to be expanded to the initial flood fill queue in ascending order of propagation layer number; S35: Based on the path status corresponding to each path to be expanded, calculate the priority score of each path to be expanded included in the initial flood filling queue using the following formula, where the priority score of the Kth path to be expanded is... for: ; S36: Based on the calculation results of step S35, sort and adjust the paths to be expanded in the initial flood filling queue according to the priority score from large to small, and arrange the paths to be expanded with the same priority score or the priority score difference is less than the preset tolerance in the order of the number of propagation layers from small to large to obtain the flood filling queue.

[0042] It should be noted that, in order to implement path search, each path to be expanded is abstracted as a path state. The path status should at least include the current path to be extended, the current end-propagation module, the current directional domain or directional bucket, the current channel label, the current path layer number, the minimum remaining layer number to the detection module, the current cumulative contribution, the upper bound of the remaining maximum potential contribution, the wavelength label, the polarization label, and path history information. It can take the following form:

[0043] in, Let K be the set of propagation modules contained in the Kth path to be expanded; This is the end-propagation module for the Kth path to be extended; This represents the direction domain or direction bucket for the Kth path to be expanded. Let K be the channel label of the Kth path to be expanded; Let K be the path layer number of the Kth path to be expanded; Let be the minimum number of remaining propagation layers from the end propagation module to the detection module in the Kth path to be extended; The cumulative contribution to the Kth path to be expanded; This is the maximum upper bound estimate of the state propagation for the Kth path to be expanded; Let $\overall_bound$ be the state bound corresponding to the $K$-th path to be expanded, and $\overall_bound$ be the state bound. ; This is the historical encoding of the Kth path to be expanded; and These are the wavelength label and polarization label for the Kth path to be extended, respectively.

[0044] During initialization, in addition to generating an initial state containing only the light source module, it is preferable to first perform a reverse traversal on the propagation graph starting from the detector module to obtain the minimum remaining propagation layers from any module to the detector. This layer information serves as target guidance information for subsequent flood filling, prioritizing the discovery of paths from the light source module to the detector module more quickly. If a direct propagation edge exists between the light source module and the detector module, a "light source module — detector module" path state is directly generated and prioritized for queuing; otherwise, paths with lower next-hop propagation module values ​​are prioritized. The path status. Through this reverse layer number marker, flood filling is no longer a aimless spread, but forms a directional path search mechanism that "starts from the light source module and prioritizes advancing towards the front edge that is closer to the detector module".

[0045] From the initial flood filling queue Extract the current state from Read the end module of the current path The set of all downstream adjacent modules in the non-sequential module propagation graph. For each adjacent module First, the basic coupling is calculated along the directed propagation edge of the geometry. Then, based on the composite properties of the current propagation module and the current state label, it is dynamically expanded into sub-channels such as transmission, specular reflection, scattering, and diffraction, thereby forming sub-states.

[0046] This invention employs a physically weighted, target-guided floodfill algorithm, rather than ordinary unweighted floodfill. Specifically, it first advances layer by layer according to the number of path layers, and then sorts the paths by their upper bound scores within the same layer. Since each additional propagation module in the path to be expanded typically introduces an additional interaction response term or propagation coupling term less than 1, paths directly from the light source to the detector, and paths that pass through only a few propagation modules and have a shorter remaining number of layers, usually have a larger theoretical upper bound and a higher practical contribution, and can be prioritized for searching and evaluation. Thus, the floodfill algorithm in this invention does not play a simple traversal role, but rather serves as the search backbone for "prioritizing the discovery of paths that are effective and contribute significantly to the detector."

[0047] Initial flood fill queue The scheduling rule adopted is "first expand layer by layer according to the number of paths, and then sort by detector guidance score within the same layer." Let... For state The number of propagation modules already experienced, This represents the minimum number of remaining propagation layers from the current end module to the detector. Given an upper bound on the contribution of this state, the priority score can be defined as follows: ; First press Layered advancement, prioritizing pop-ups within the same layer. Larger states; when scores are similar, pop up first. Smaller states. Since each additional propagation module typically introduces an interaction response term and a propagation coupling term less than 1, paths directly from the light source to the detector, as well as paths that pass through only a few propagation modules and have a short remaining number of layers, are usually searched earlier and generally have a larger energy contribution. Therefore, short paths are not only more likely to become dominant paths in engineering experience, but their theoretical contribution upper bound is also usually larger, making them suitable as priority expansion targets for flood filling. Through this mechanism, the algorithm first searches for direct paths, single-interaction paths, and other high-contribution short paths; if these short paths cannot explain the main energy of the detector, it continues to expand to deeper, weakly contributing paths. At the same time, the high-contribution short paths found can also serve as a reference benchmark for subsequent pruning, thereby further compressing the invalid search space.

[0048] In some embodiments, step S4 specifically includes the following steps: S41: Select a path to be expanded from the flood fill queue and calculate the individual contribution of each directed propagation edge contained in the current path to be expanded. S42: Calculate the total contribution of the current path to be expanded based on the individual contributions of each directed propagation edge: ; in, The total contribution of the Pth path to be expanded. Let n be the light source flux, and n be the total number of propagation modules. For the k-th propagation module, For the surface interaction response of the k-th propagation module, For the single-segment propagation coupling term of the k-th propagation module; S43: Replace the current path to be expanded with the next path to be expanded, and repeat steps S41-S42 until the total contribution of all paths to be expanded is obtained.

[0049] In some embodiments, in step S41, when the directed propagation edge is from the propagation module To the propagation module Extending this, the formula for calculating the individual contribution of the current directed propagation edge is as follows: The current contribution of the directed propagation edge to the mirror reflection channel is: ; in, This represents the current interactive response of the directed propagation edge under the mirror reflection channel. For propagation module i at wavelength λ and incident angle The reflectivity function under the condition of polarization state pol; The current contribution of the directed propagation edge to the transmission channel is: ; in, This represents the interactive response of the current directed propagation edge in the transmission channel. For propagation module i at wavelength λ and incident angle Transmittance function under polarization state pol; Set up a propagation module The incident flux is The total scattering flux of the current directed propagation edge is: ; in, The total scattered flux is obtained by converting the incident flux into the propagation module i. The total integral scattering coefficient of propagation module i represents the proportion of the incident flux allocated to the scattering channel; Propagation Module The corner domain The scattering response is normalized as follows: ; in, For the scattered light from propagation module i to enter the angular region where propagation module j is located. The normalized scattering response, The total integrated scattering coefficient of propagation module i, For the scattered light rays from propagation module i, the full-space output solid angle is given. For propagation module i in the incident direction and the two-way scattering distribution function at the emission direction ω, The incident direction, The direction of launch; The current contribution of the directed propagation edge in the discrete diffraction order sub-channels is: ; in, The interactive response corresponds to the discrete diffraction order m. For propagation module i at wavelength λ and incident angle The m-th order diffraction efficiency function under the condition of polarization state pol.

[0050] In some embodiments, step S4 specifically includes: S41: Estimate the upper bound of the contribution of the Kth path to be expanded using the following formula. : ; ; in, Let K be the current cumulative contribution value of the Kth path to be expanded. To estimate the upper bound of the maximum product for the Kth path to be expanded, Let t be the minimum number of remaining propagation layers from the end propagation module of the Kth path to be extended to the corresponding propagation module of the detector, where t is the index of the remaining propagation layer. For the first Maximum single-layer multiplier upper bound estimate within each remaining propagation layer; S42: If the upper bound of the contribution of the Kth path to be expanded is less than the first preset threshold. If the current path to be expanded is not found, the expansion of the current path to be expanded will be terminated, and the Kth path to be expanded will be considered an invalid path; otherwise, the Kth path to be expanded will be considered a valid path. S43: Replace the k-th path to be expanded with the (k+1)-th path to be expanded, and repeat steps S41-S42 until all valid paths are obtained; S44: Determine if there are duplicate valid paths among all valid paths. If so, merge the duplicate valid paths and update the valid paths. Otherwise, proceed directly to step S45. The conditions for merging states include at least the following: the current end propagation modules of the two states are the same, the current channel labels are the same, and the difference in the direction domain is not greater than a preset threshold. The wavelength difference is not greater than the preset threshold. Furthermore, the polarization labels must be consistent. If the above conditions are met, the two states are merged into an equivalent state, the cumulative contribution is accumulated, and the upper bound of the contribution is updated.

[0051] S45: Pair up each propagation module in the current valid path, count the number of back propagations between the paired propagation modules, and delete valid paths with a back propagation count greater than the second preset threshold (a value set manually according to requirements). S46: Replace the current valid path with the next valid path, repeat step S45, complete the further update of valid paths, and obtain the final set of valid paths.

[0052] In some embodiments, in step S42, the formula for calculating the first preset threshold is: ; in, The upper bound of the contribution to the flood filling queue is the maximum value. α is a preset proportionality constant, 0 < α < 1.

[0053] It should be noted that if the current path to be extended originates from the propagation module... Extended to the propagation module Therefore, it is necessary to first calculate the interaction response between the two propagation modules, and then multiply it by the single-segment propagation coupling term to obtain the contribution of that segment. The interaction response of different surface mechanisms is defined as follows.

[0054] For the specular reflection channel: ; For the transmission channel: ; For the scattering channel, use and Joint constraints on total quantity and angular distribution. Assume a propagation module. The incident flux is The total scattering flux is: ; like Indicates a given incident direction Downwards points to the direction of emission The angular distribution points towards the downstream propagation module. The corner domain The scattering response can be normalized as: ; in, This represents the effective half-space or allowed scattering space corresponding to the propagation module. Using this form ensures that the total scattered energy is controlled... Constraints, and can also be passed It reflects the differences in angular distribution.

[0055] For discrete diffraction order channels: ; If a path P to be extended consists of a module sequence Composition, in which For light source module, If the target evaluation module is the detector module or the target evaluation module (a propagation module of interest to the user other than the detector), then the total contribution of the path P to be extended can be expressed as: ; Whenever the path P to be expanded expands by a module, the current cumulative contribution is multiplied by the newly added contribution. and This yields new path contribution values. and After normalization, they are all used as dimensionless multipliers, therefore With source flux Keep the same dimensions.

[0056] To avoid wasting computational resources on a large number of low-contribution branches, this invention continuously evaluates the total upper bound of the state during the path expansion process. .set up For the current state in the remaining The maximum possible upper bound of the product within each propagation layer is: ; ; in, Indicates the first The maximum upper bound of the single-layer multipliers that may occur within each remaining propagation layer can be given by a combination of the maximum reflectivity, maximum transmittance, maximum scattering response, maximum diffraction efficiency, and maximum propagation coupling term. If the current state satisfies:

[0057] Therefore, it is assumed that even if this branch continues to expand, its final contribution to the detection module will be below the threshold. The branch can be directly pruned to terminate the process. This threshold can be either an absolute threshold or a relative threshold; in a preferred embodiment, after obtaining the first direct path or the first batch of short paths, the maximum contribution of the two is recorded as... And order: ; In this way, once a high-contribution short path is identified early, it can be used to immediately update the pruning threshold, and weak branch pruning can be performed on a large number of theoretically achievable but practically low-contribution branches. If two states reach the same module and their direction domain, band, polarization, and channel labels are approximately consistent within a preset tolerance (set manually as needed), then state merging is performed; for repeated module sequences that have already appeared or backpropagation exceeding the allowed number of times, loop closure suppression is performed.

[0058] In some embodiments, when the terminal module of the effective path is a detector module or a target evaluation module, the contribution of the effective path is accumulated to the corresponding pixel, region, or receiving channel. If the detector coordinates are used... Indicating the receiving position, the total stray light from the detector can be expressed as: ; in, For all final hit detector locations The effective path set. The final output includes: stray light distribution on the detection surface, ranking of major contributing paths, ranking of major contributing modules, comparison results of stray light under different bands or polarization conditions, and location results of sensitive modules for engineering rectification.

[0059] The above steps constitute the complete technical solution of this invention. Its core lies in: the flood filling algorithm is used to organize and control the system expansion of complex non-sequential paths; the radiative transfer formula is used to quantitatively calculate the energy contribution of each effective path; and contribution upper bound pruning and state merging are used to ensure computational efficiency and engineering interpretability.

[0060] Example 1 The following describes in detail the specific implementation of the present invention using the Cooke three-lens imaging system with a lens barrel as an example.

[0061] 1) System Composition like Figure 3 As shown, a Cooke three-lens imaging system with a telescope barrel is used as the analysis object. The Cooke three-lens imaging system with a telescope barrel, from the object side to the image side, includes: an off-axis stray light source 1, a first positive lens 3, a second negative lens 4, a third positive lens 5, a telescope barrel structure 2 surrounding the three lenses, and an image plane detector. Besides the cylindrical section, the telescope barrel also includes lens mounting steps and a local constriction structure; these structures can serve as important scattering or blocking modules in stray light propagation. For ease of explanation, this embodiment uses an axisymmetric structure for illustration, and the axisymmetric surface is discretized into several propagation modules in the three-dimensional model.

[0062] The propagation module is constructed using the following approach: "Optical working surface module + mechanical scattering module + detector module + light source module". Specifically, the six refractive surfaces serve as optical working modules; the outer edges of the three lenses serve as lens edge scattering modules; the inner wall of the lens barrel is divided into several inner wall scattering modules along the axial direction; the lens barrel steps and local constriction points serve as independent mechanical scattering modules; the image plane serves as the detector module; and off-axis stray light source 1 serves as the light source module. This design effectively reflects both residual reflections and micro-scattering from the lens surfaces, as well as typical stray light paths caused by the lens barrel. The main parameters of the six refractive surfaces and the image plane are described in Table 1.

[0063] Table 1

[0064] Furthermore, the glass material of the first positive lens 3 is SK16_SCHOTT, the glass material of the second negative lens 4 is F4_HOYA, and the glass material of the third positive lens 5 is SK16_SCHOTT.

[0065] 2) Optical property configuration In this embodiment, each refractive surface module simultaneously possesses three types of properties: primary transmission, residual specular reflection, and surface micro-scattering. For each refractive surface module, at least transmittance, residual reflectance, and total integrated scattering are configured. Applicable wavelength bands and polarization labels. For surface micro-scattering, empirical methods can be used. Model or actual measurement The data is described in a table.

[0066] The inner wall module of the lens barrel, the stepped module of the lens barrel, and the edge module of the lens are mainly for scattering and are primarily configured with... and Absorption rate and surface blackening correction coefficients can be added if necessary. The image plane detector module only receives path contributions and does not continue to emit. Off-axis stray light source 1 defines its position, initial radiant flux, emission angle distribution, wavelength label, and polarization label.

[0067] In practical engineering, It can be directly provided from measured data tables, supplier data, or empirical models; It can be estimated from surface roughness and process grade, or obtained directly from scattering test data. If complete measurement data is temporarily unavailable, a normalized empirical model can be used as an approximation to ensure that the total scattering and emission angle distribution can be included in the propagation calculation.

[0068] In this embodiment, for ease of explanation, the following simplified empirical scattering model can be used: ; in, The angle at which the scattering deviates from the principal direction; and These are empirical parameters; according to Constraints are determined through normalization. This empirical scattering model is mainly used for scattering modules such as the inner wall of the lens barrel, lens barrel steps, and lens edges; for the main transmission and residual specular reflection from the refractive surface, transmittance and reflectance functions are still used respectively. This expression is only for illustrative purposes and does not limit the invention to using higher-precision measured data. Table or Model.

[0069] Then, a directed propagation edge determination is performed on any two modules. If the two modules satisfy the conditions of geometric visibility, directional coverage, non-zero effective solid angle, and contribution upper bound higher than the threshold, then a geometric propagation edge is established from the former to the latter, and saved. , , , These are the basic parameters. For the same geometric propagation edge, it is dynamically expanded into sub-channels such as transmission, specular reflection, scattering, and diffraction in the subsequent path expansion stage based on the current state.

[0070] 3) Initial state settings and flood filling expansion An initial path state is generated based on the off-axis light source module (stray light is often caused by off-axis light sources outside the imaging field of view, hence the specific term "off-axis light source"). The initial path state only contains the light source module, and the initial contribution is taken as the radiant flux or normalized unit flux of the light source module. The initial upper bound is estimated based on the maximum possible transmittance, reflectance, scattering rate, and solid angle envelope of the system. The initial state is then added to the queue to be expanded. middle.

[0071] Simultaneously, starting from the detector module, a reverse traversal is performed on the non-sequential module propagation graph to calculate the minimum remaining propagation layers from each propagation module to the detector module. The minimum remaining propagation layer number is used as the target guidance label for the corresponding propagation module. The information of the i-th propagation layer is written into the state or attached to the module as target guidance information. Then, forward flood filling is performed: first, it searches whether the light source can reach the detector directly through the lens tube gap or edge leakage; if not, it prioritizes searching for paths that only pass through one propagation module, such as "off-axis light source → local surface element of the inner wall of the lens tube → detector" and "off-axis light source → rear lens surface → detector"; then it continues to expand the two-stage and three-stage interaction paths.

[0072] Then from the queue to be expanded The current propagation module's state is extracted, and all downstream adjacent modules of the current end module of the path to be extended are read. These are then combined with the current propagation module's interaction channels to generate sub-states. For example, for a lens surface, transmission sub-states, residual reflection sub-states, and surface scattering sub-states can be generated simultaneously; for the inner wall module of the lens barrel, only scattering sub-states are generated; for lens barrel steps and local constriction locations, scattering sub-states or termination states can be generated based on edge illumination conditions. Each newly generated sub-state simultaneously updates three quantities: the new cumulative path contribution; the upper bound of the contribution for continued propagation from this sub-state to the detector; and the minimum residual propagation layers from the current end module to the detector.

[0073] In this embodiment, flood filling does not evenly expand all branches, but instead employs a weighted expansion strategy of "path layer priority + target-oriented scoring sorting". Specifically, queue scheduling prioritizes path layer priority. Progress through tiers, then within each tier, based on scores. The sorting method, among which Since short paths tend to have less decay and a larger upper bound, the algorithm can quickly identify dominant direct paths and short paths with one interaction, and use their contribution values ​​to update the pruning threshold of subsequent weak branches.

[0074] In the Cooke three-lens imaging system with a lens barrel, the algorithm first checks whether an off-axis source can directly reach the detector through lens barrel gaps or edge leakage. Secondly, it prioritizes searching for short paths with only one interaction, such as "off-axis source—local surface element on the inner wall of the lens barrel—detector" and "off-axis source—rear lens surface—detector." Once these paths are found, their contribution values ​​are immediately used as references for subsequent pruning. Only when the short paths cannot adequately explain stray light on the detector does the flood filling front continue to advance towards paths with two or three interactions, such as "off-axis source—inner wall of the lens barrel—rear lens—detector" and "off-axis source—front lens—inner wall of the lens barrel—rear lens—detector." Therefore, in this embodiment, the flood filling algorithm does not simply traverse paths indiscriminately, but actually undertakes the hierarchical, directional, and priority search task from the light source module to the detector module.

[0075] 4) Single-segment propagation calculation and path contribution recursion When the path is from the propagation module Extended to the propagation module First, calculate the interactive response according to the current channel type. Then coupled with single-segment propagation term Multiply. If the current channel is specular reflection, use the reflectivity function; if the current channel is transmissive, use the transmissivity function; if the current channel is scattering, use the multiplication function. and The normalized integral calculation points to The scattering response in the angular domain.

[0076] For example, when the inner wall module of the lens barrel When illuminated by an off-axis source, it points to a certain lens surface module. scattering contribution for: ; in, Off-axis light source flux, For the inner wall module of the lens barrel Pointing to a certain lens surface module Surface scattering interaction response, When the inner wall module of the lens barrel When illuminated by an off-axis source, it points to a certain lens surface module. The single-segment propagation coupling term, i.e., the inner wall module of the cylinder. With lens surface module Propagation coupling terms between segments.

[0077] When the lens surface module Then, it is coupled to the detector module through residual reflection or transmission. When this happens, the contribution of the corresponding path continues to be calculated recursively. If the path is... = (light source module, inner wall module of lens barrel, lens surface module, detector module), then its contribution is obtained by the following formula: ; in, The stray light flux that ultimately contributes to detector module Md by the effective path P. This represents the initial radiant flux or normalized source flux of module 1 of the off-axis stray light source. The outgoing response or source direction weight of the off-axis stray light source 1 pointing to the inner wall module Mb of the lens barrel. This is the single-segment propagation coupling term from the off-axis stray light source module 1 to the inner wall module Mb of the lens barrel. The interaction response of the lens surface module Ms to the detector module Md after residual reflection or transmission. This is the single-segment propagation coupling term from the lens surface module Ms to the detector module Md.

[0078] The above expression illustrates that this invention does not separately count the local quantities of each hop, but rather unifies and concatenates all interactive responses and propagation coupling terms of the entire path to obtain the truly effective path contribution to the detector.

[0079] 5) Pruning, state merging, and loop control As the expansion depth increases, the number of paths may increase rapidly. To improve the efficiency of the engineering solution, this embodiment performs three types of control immediately after each sub-state is generated. The first type is upper bound pruning: when Less than the threshold When this happens, the expansion of the state is stopped; where It can be dynamically updated after obtaining the first high-contribution short path. The second category is state merging: if two states are in the same current module and have similar direction domains, bands, and polarization categories, they are merged into an equivalent state, and their contributions are accumulated. The third category is loop closure control: suppression is performed on repeating module sequences or backpropagation exceeding the allowed number of times.

[0080] These three types of control significantly reduce theoretically existing but practically insignificant propagation branches, allowing algorithmic resources to be focused on more meaningful paths. In particular, the first discovered direct paths and single-interaction paths not only provide major contributions but also directly participate in threshold updates and weak branch pruning. Therefore, the floodfill algorithm and radiative transfer recursion form a closed-loop synergy in this embodiment.

[0081] 6) Detector contribution accumulation and engineering output When an effective path reaches the detector module, its contribution is projected onto the corresponding pixel, region, or receiving channel. For array detectors, contributions can be recorded by pixel grid; for integral evaluation surfaces, contributions can be recorded by the integral of the receiving region. After the path contributions are accumulated, the following engineering results can be output: (1) Stray light distribution map of the detection surface; (2) Ranking table of major contribution paths; (3) Ranking table of major contributing modules; (4) Comparison results of stray light under different wavelengths, different polarizations, or different operating conditions; (5) Location results of sensitive modules used for structural rectification.

[0082] The above information can be used to more quickly identify high-contribution dominant paths. By employing layered flood filling, detector inverse layer labeling, and contribution upper bound sorting, paths that run directly from the light source to the detector or pass through only a few propagation modules can be prioritized. These paths typically have less attenuation and greater contribution, thus facilitating the rapid identification of the main stray light source.

[0083] This invention has completed the engineering simulation process design and program prototype construction, and can complete the complete calculation process from module discretization, propagation map establishment, path expansion, single-segment radiative transfer recursion to detector contribution accumulation on the Cooke three-piece lens imaging system with lens barrel.

[0084] Theoretically, this invention has a clear physical basis and is feasible. On the one hand, the calculation of the single-segment propagation contribution is directly based on mature optical quantities such as BSDF, TIS, reflectivity, transmittance, diffraction efficiency, effective solid angle, and occlusion relationship; on the other hand, flood filling, priority queues, state hash merging, and graph search are all mature and feasible data structures and algorithms. Therefore, this invention is feasible at the software implementation level.

[0085] The suggested simulation verification method is as follows: (1) Select a typical imaging system that includes a lens surface, inner wall of the lens barrel, light-blocking ring, cover plate and detector, and establish its three-dimensional geometric model; (2) Assign reflectivity, transmittance, BSDF, TIS and diffraction efficiency parameters to each surface; (3) Simulate the same external stray light condition using the method of this invention and the conventional non-sequential Monte Carlo method respectively; (4) Compare the differences between the two methods in terms of stray light distribution, consistency, critical path discovery capability and computational efficiency on the detector surface; (5) Statistically analyze the top N high contribution paths output by the method of this invention, and check whether these paths are consistent with the path analysis results of engineering experience and structural hotspot areas.

[0086] Based on the aforementioned path state definition, detector reverse layer number marking, contribution upper bound ranking, and dynamic pruning mechanism, this invention can prioritize expanding propagation paths with higher theoretical contributions and fewer remaining layers reaching the detector in the non-sequential module propagation map, and can output the effective path set, path contribution ranking, and module contribution ranking. Compared with the non-sequential Monte Carlo method that relies solely on random sample hits, this invention does not use random ray hits as the sole condition for discovering low-probability critical paths, thus reducing the risk of missing multiple weak reflections, scattering, or inefficient diffraction paths. By accumulating the effective path contributions, this invention can not only obtain the stray light distribution on the detector surface but also provide path attribution and module attribution information, providing a basis for subsequent tube extinction, baffle ring adjustment, edge blackening, and coating improvement.

[0087] This invention is not limited to the above-described embodiments. Various modified designs and alternative solutions exist, all of which should be considered to fall within the protection scope of this invention.

[0088] (1) At the discrete level of the module, the propagation module can be defined as a whole surface, a local surface element, a volume scattering unit or a functional substructure of interest to the user, and is not limited to this embodiment.

[0089] (2) At the search strategy level, flood filling can adopt first-in-first-out queue, priority queue, hierarchical queue or multi-queue scheduling method; priority can be determined by current contribution, remaining upper bound, target coupling strength, path depth or a combination thereof.

[0090] (3) At the propagation direction modeling level, either the inter-module discrete edge method or the directional bucket, solid angle grid or angular spectrum discrete method can be used to represent the directional expansion from the current module to the downstream module.

[0091] (4) In terms of physical modeling, in addition to BSDF / TIS, BRDF, BTDF, polarization Mueller matrix, wavelength discrete model, temperature condition model or surface contamination degradation model can also be introduced.

[0092] (5) At the path organization level, it can be filled from the light source module in the forward direction or from the detection module in the reverse direction, and then meet and match with the forward path of the source module at the intermediate module.

[0093] (6) At the algorithm fusion level, the method of this invention can be used in combination with the Monte Carlo method. For example, the method of this invention can be used to search and lock high contribution paths first, and then random refinement sampling can be used on local paths.

[0094] (7) In addition to being applicable to stray light analysis of imaging lenses, the present invention can also be used in infrared systems, laser systems, projection systems, illumination systems, space remote sensing systems, micro-nano diffraction element systems, photolithography illumination systems, AR / VR optical systems, etc.

[0095] It should be understood that the various forms of processes shown above can be used to reorder, add, or delete steps. For example, the steps described in this invention disclosure can be executed in parallel, sequentially, or in different orders, as long as the desired result of the technical solution disclosed in this invention can be achieved, and this is not limited herein.

[0096] The specific embodiments described above do not constitute a limitation on the scope of protection of this invention. Those skilled in the art should understand that various modifications, combinations, sub-combinations, and substitutions can be made according to design requirements and other factors. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of this invention should be included within the scope of protection of this invention.

Claims

1. A method for analyzing stray light propagation in a non-sequential optical system, characterized in that: Specifically, the steps include the following: S1: Establish a three-dimensional geometric model of the optical system to be analyzed, and discretize the three-dimensional geometric model into no less than three propagation modules, and configure the attributes of each propagation module; S2: Perform pairwise matching verification on all propagation modules and establish directed propagation edges for valid propagation module pairs to obtain the non-sequence module propagation graph; S3: Based on the directed propagation edges of the non-sequential module propagation graph, extract the paths to be expanded, abstract each path to be expanded into a path state to construct an initial flood filling queue, and sort the initial flood filling queue by priority to obtain the flood filling queue. S4: Calculate the total contribution of each path to be expanded contained in the flood filling queue, and perform dynamic pruning based on the upper bound of the contribution of each path to be expanded to obtain a set of effective paths. S5: Calculate the total stray light of the detectors in the optical system to be analyzed based on the stray light energy of each effective path contained in the effective path set.

2. The method for analyzing stray light propagation in a non-sequential optical system according to claim 1, characterized in that: In step S1, the types of propagation modules include at least an optical working surface module, a detector module, and a light source module.

3. The method for analyzing stray light propagation in a non-sequential optical system according to claim 2, characterized in that: In step S1, the optical properties of the optical working surface module include at least one or more of the following parameters: reflectivity, transmittance, absorptivity, total integrated scattering, two-way scattering distribution function, two-way transmission distribution function, diffraction order efficiency, and medium transmission factor. The boundary properties of the light source module include at least one or more of the following parameters: initial radiant flux, angular radiation distribution characteristics, wavelength or band, polarization state, and effective radiation aperture size of the light source; The boundary properties of the detector module include at least one or more of the following parameters: detector surface position, effective receiving size, pixel region discretization method, receiving direction constraint, and statistical rules.

4. The method for analyzing stray light propagation in a non-sequential optical system according to claim 1, characterized in that: Step S2 specifically includes the following steps: S21: Calculate the single-segment propagation parameters of the current propagation module pair: ； ； Where i and j are the index numbers of the two propagation modules contained in the current propagation module. For the single-segment propagation parameters of the current propagation module pair, The visibility factor for the current propagation module pair. The medium transmittance factor for the current propagation module pair. Let i be the solid angle corresponding to propagation module i at propagation module j. The projection geometry factor of the current propagation module pair. These are optional modifications for the current propagation module pair. Let j be the equivalent area of ​​the propagation module. Let be the angle between the direction of the line connecting the center of propagation module i to the center of propagation module j and the normal direction of propagation module j. The center distance of the current propagation module pair; S22: Determine whether the single-segment propagation parameter of the current propagation module pair is zero. If not, then take the current propagation module pair as a valid propagation module pair and establish a directed propagation edge for the current valid propagation module pair. Otherwise, execute step S23. S23: Replace the current propagation module pair with the next propagation module pair, repeat steps S21-S22, obtain the directed propagation edges of all valid propagation module pairs, and construct a non-sequence module propagation graph based on all valid propagation module pairs and their directed propagation edges.

5. The method for analyzing stray light propagation in a non-sequential optical system according to claim 1, characterized in that: Step S3 specifically includes the following steps: S31: Traverse all directed propagation edges, take the propagation module corresponding to the light source module as the starting point of the path, and extend the initial path along the directed propagation edge as the initial path to be extended, and obtain the initial path to be extended set. The initial path to be extended includes at least: the initial path to be extended containing only the propagation module corresponding to the light source module, and the initial path to be extended consisting of the propagation module corresponding to the light source module and the downstream propagation module connected to the propagation module corresponding to the light source module. S32: Each path to be expanded is abstracted into a path state using the following formula, where the path state corresponding to the Kth path to be expanded is... : ; in, Let K be the set of propagation modules contained in the Kth path to be expanded; This refers to the end propagation module that has been extended to the Kth path to be extended; This represents the direction domain or direction bucket for the Kth path to be expanded. Let K be the channel label of the Kth path to be expanded; Let K be the path layer number of the Kth path to be expanded; Let be the minimum number of remaining propagation layers from the end propagation module to the detection module in the Kth path to be extended; The cumulative contribution to the Kth path to be expanded; This is the maximum upper bound estimate of the state propagation for the Kth path to be expanded; Let $\overall_bound$ be the state bound corresponding to the $K$-th path to be expanded, and $\overall_bound$ be the state bound. ; This is the historical encoding of the Kth path to be expanded; and These are the wavelength label and polarization label for the Kth path to be extended, respectively. S33: Starting from the propagation module corresponding to the detector module, perform a reverse breadth-first traversal of the non-sequence module propagation graph and record the minimum number of remaining propagation layers from each propagation module to the detector module. S34: Add the paths to be expanded to the initial flood fill queue in ascending order of propagation layer number; S35: Based on the path status corresponding to each path to be expanded, calculate the priority score of each path to be expanded included in the initial flood filling queue using the following formula, where the priority score of the Kth path to be expanded is... for: ； S36: Based on the calculation results of step S35, sort and adjust the paths to be expanded in the initial flood filling queue according to the priority score from large to small, and arrange the paths to be expanded with the same priority score or the priority score difference is less than the preset tolerance in the order of the number of propagation layers from small to large to obtain the flood filling queue.

6. The method for analyzing stray light propagation in a non-sequential optical system according to claim 1, characterized in that: Step S4 specifically includes the following steps: S41: Select a path to be expanded from the flood filling queue and calculate the individual contribution of each directed propagation edge contained in the current path to be expanded. S42: Calculate the total contribution of the current path to be expanded based on the individual contributions of each directed propagation edge: ； in, The total contribution of the Pth path to be expanded. Let n be the light source flux, and n be the total number of propagation modules. For the k-th propagation module, For the surface interaction response of the k-th propagation module, For the single-segment propagation coupling term of the k-th propagation module; S43: Replace the current path to be expanded with the next path to be expanded, and repeat steps S41-S42 until the total contribution of all paths to be expanded is obtained.

7. The method for analyzing stray light propagation in a non-sequential optical system according to claim 1, characterized in that: In step S41, when the directed propagation edge is from the propagation module To the propagation module Extending this, the formula for calculating the individual contribution of the current directed propagation edge is as follows: The current contribution of the directed propagation edge to the mirror reflection channel is: ； in, This represents the current interactive response of the directed propagation edge under the mirror reflection channel. For propagation module i at wavelength λ and incident angle The reflectivity function under the condition of polarization state pol; The current contribution of the directed propagation edge to the transmission channel is: ； in, This represents the interactive response of the current directed propagation edge in the transmission channel. For propagation module i at wavelength λ and incident angle Transmittance function under polarization state pol; Set up a propagation module The incident flux is The total scattering flux of the current directed propagation edge is: ； in, The total scattered flux is obtained by converting the incident flux into the propagation module i. The total integral scattering coefficient of propagation module i represents the proportion of the incident flux allocated to the scattering channel; Propagation Module The corner domain The scattering response is normalized as follows: ； in, For the scattered light from propagation module i to enter the angular region where propagation module j is located. The normalized scattering response, The total integrated scattering coefficient of propagation module i, For the scattered light rays from propagation module i, the full-space output solid angle is given. For propagation module i in the incident direction and the two-way scattering distribution function at the emission direction ω, The incident direction, The direction of launch; The current contribution of the directed propagation edge in the discrete diffraction order sub-channels is: ； in, The interactive response corresponds to the discrete diffraction order m. For propagation module i at wavelength λ and incident angle The m-th order diffraction efficiency function under the condition of polarization state pol.

8. The method for analyzing stray light propagation in a non-sequential optical system according to claim 1, characterized in that: Step S4 specifically includes: S41: Estimate the upper bound of the contribution of the Kth path to be extended using the following formula. : ； ； in, Let K be the current cumulative contribution value of the Kth path to be expanded. To estimate the upper bound of the maximum product for the Kth path to be expanded, Let t be the minimum number of remaining propagation layers from the end propagation module of the Kth path to be extended to the corresponding propagation module of the detector, where t is the index of the remaining propagation layer. For the first Maximum single-layer multiplier upper bound estimate within each remaining propagation layer; S42: If the upper bound of the contribution of the Kth path to be expanded is less than the first preset threshold. If the current path to be expanded is not found, the expansion of the current path to be expanded will be terminated, and the Kth path to be expanded will be considered an invalid path; otherwise, the Kth path to be expanded will be considered a valid path. S43: Replace the k-th path to be expanded with the (k+1)-th path to be expanded, and repeat steps S41-S42 until all valid paths are obtained; S44: Determine if there are duplicate valid paths among all valid paths. If so, merge the duplicate valid paths and update the valid paths. Otherwise, proceed directly to step S45. S45: Pair up each propagation module in the current valid path, count the number of back propagations between the paired propagation modules, and delete valid paths with a number of back propagations greater than the second preset threshold. S46: Replace the current valid path with the next valid path, repeat step S45, complete the further update of valid paths, and obtain the final set of valid paths.

9. The method for analyzing stray light propagation in a non-sequential optical system according to claim 8, characterized in that: In step S42, the formula for calculating the first preset threshold is: ； in, The upper bound of the contribution to the flood filling queue is the maximum value. α is a preset proportionality constant, 0 < α < 1.

10. The method for analyzing stray light propagation in a non-sequential optical system according to claim 6, characterized in that: In step S5, the formula used to calculate the total stray light of the detectors in the optical system to be analyzed is as follows: ； in, This represents the total stray light accumulated by the detector at pixel (x,y). Let P be the pixel index, and P be the valid path to the pixel (x, y). This is the set of valid paths that hit the pixel (x, y).

Images