Simulation method for multi-source coupling imaging of laser radar based on inverse monte carlo path tracing

Abstract

The application discloses a kind of based on laser radar multi-source coupling imaging simulation method of reverse Monte Carlo path tracking.It is unified to build target and background triangle facet three-dimensional scene and mark medium type;Support three kinds of atmospheric modes of vacuum, Malkmus narrow band and complete atmospheric radiation transmission and construct atmospheric radiation lookup table.Shoot ray from each pixel of detector reversely, after intersection, differential processing is carried out according to medium type, combined with bidirectional reflection distribution function importance sampling and Russian roulette control multiple ejection, coupling target, background, atmosphere and solar radiation;According to medium type, light source type and ejection order, the radiation contribution is decomposed into six signal components of laser direct reflection, laser background reflection, sky path radiation, solar reflection, atmospheric backscattering and multiple ejection;According to time-of-flight mapping and Gaussian pulse convolution synthesis time-domain full waveform, output intensity image, distance image, time-of-flight image, signal decomposition file and full waveform data.

Description

Technical Field

[0001] This invention relates to the field of lidar imaging simulation technology. Background Technology

[0002] LiDAR imaging simulation has significant application value in aerospace, remote sensing and mapping, and reconnaissance fields. Simulation allows for the pre-evaluation of imaging effects under different atmospheric conditions, target materials, and detection parameters before actual hardware development, significantly reducing testing costs and timelines.

[0003] Currently, the mainstream lidar imaging simulation methods mainly include the following three categories: (1) Simulation methods based on forward ray tracing: such as forward ray tracing software like HELIOS++, which emits rays from a laser source into the scene and traces the propagation path of each ray. This method is highly efficient when dealing with single scattering, but in complex scenes, it requires explicit enumeration of all possible scattering paths. The computational complexity increases exponentially with the number of scatterings, making it difficult to efficiently simulate the global illumination effect caused by multiple scatterings.

[0004] (2) Simplified simulation method based on analytical radar equations: The received power of each pixel is directly calculated using the LiDAR (Light Detection and Range) radar equations. This method is extremely fast, but it cannot handle multiple scattering, shadow effects, and complex bidirectional reflection distribution function materials, resulting in low simulation fidelity.

[0005] (3) Simulation method based on pre-calculated lookup table: Radiation transfer lookup tables under different conditions are pre-calculated, and interpolation is performed during simulation. This method has poor flexibility and is difficult to adapt to arbitrary geometric scenarios and non-Lambertian surfaces.

[0006] More importantly, existing simulation methods generally suffer from the following problems: (1) Separation of target simulation and background simulation: Most simulation tools process the aircraft target and the terrain / sea background separately, without simultaneously coupling the radiation transmission processes of the target, background, atmosphere and sun within a unified physical framework. This results in the inaccurate capture of multi-source coupling effects (such as multiple reflections between the target and background, the contribution of sun-background scattering to the target, etc.).

[0007] (2) Single dimension of output data: Most simulation tools only output intensity images and distance images, lacking the ability to quantitatively decompose the radiation contribution of each source during the simulation process, which limits subsequent physical analysis and algorithm development.

[0008] Therefore, there is an urgent need to develop a lidar imaging simulation method that can uniformly handle multi-source coupling of target, background, atmosphere, and sun, has high-fidelity multiple ejection simulation capability, and outputs rich dimensions. Summary of the Invention

[0009] The purpose of this invention is to overcome the shortcomings of existing technologies and provide a multi-source coupled lidar imaging simulation method based on inverse Monte Carlo path tracing. This method involves emitting rays backward from the detector pixel, simulating the radiation transmission process of multiple ejections through bidirectional reflection distribution function importance sampling and Russian roulette strategy, realizing the four-element coupling of target, terrain background, atmosphere and sun within a unified physical framework, and automatically performing quantitative decomposition of the total received power of each pixel into six source components.

[0010] To achieve the above objectives, the technical solution adopted by this invention is as follows: a simulation method for multi-source coupled imaging of lidar based on inverse Monte Carlo path tracing, comprising the following steps: Step 1: Construct a triangular mesh model of the imaging target and its corresponding spatial acceleration structure, and label the medium type attributes of the triangular mesh to form a 3D scene that is consistent with the background and can be ray-traced; at the same time, according to the simulation requirements of the imaging target, set the atmospheric mode to vacuum mode, Malkmus statistical narrowband model mode, or complete atmospheric radiative transfer model mode. The medium types include target surface, background surface, atmospheric path, and sky direction. When the atmospheric model is the complete atmospheric radiative transfer model, the complete atmospheric radiative transfer model automatically calls the external atmospheric radiative transfer process, extracts atmospheric transmittance, path radiation, and sky background radiation data, and constructs a multidimensional atmospheric radiation lookup table. Step 2: Based on the inverse Monte Carlo method, sample each pixel of the imaging target in the detector imaging plane. The sampling of any pixel starts from any pixel position and randomly reverses the ray direction within the detector's field of view to generate an initial ray. Furthermore, by utilizing a space acceleration structure, the intersection calculation between the reverse ray and the triangular surface element is performed in a unified scene containing the imaging target and the background. The location of the hit point, the geometric normal direction of the surface of the hit triangular surface element, and the medium type attribute of the triangular surface element are obtained. Based on the medium type attribute, different reverse ray radiative transmission processing is performed to obtain the transmittance and radiation of the atmospheric path through which the reverse ray passes. Next, at each hit point location of the reverse ray tracing, the importance of the reflection contribution distribution determined by the bidirectional reflection distribution function of the hit triangle surface element is sampled for the next emission direction to ensure the unbiasedness of the radiation contribution estimation. Finally, the ejection and termination of the reverse ray are controlled by Russian roulette, and the radiation contribution of the corresponding pixel is obtained after the reverse ray terminates. Step 3: During the reverse ray path tracing process in Step 2, the radiation contribution is classified and accumulated to the corresponding signal component according to the medium type, light source type and ejection order of the hit point location, so as to achieve multi-source signal decomposition. Step 4: Based on the time-of-flight distribution of the inverse ray of each pixel, synthesize the full temporal waveform data of the corresponding pixel through time bin mapping and Gaussian pulse convolution, and output the simulation results including intensity image, distance image, time-of-flight image, signal decomposition file and full waveform data file.

[0011] Furthermore, in step one, the three-dimensional scene includes a background scene constructed from a terrain digital elevation model and land cover classification data, wherein the triangular surface element is mapped using a bidirectional reflection distribution function based on the land cover classification data.

[0012] Furthermore, in step two, the execution of different reverse ray radiative transfer processes based on the medium type attributes specifically includes: If the reverse ray hits the imaging target surface, the bidirectional reflectance distribution function of the imaging target surface is sampled and the ray throughput is updated; if the reverse ray hits the background surface, the bidirectional reflectance distribution function of the background surface is sampled and the background signal is accumulated; if the reverse ray passes through the atmospheric path, the multidimensional atmospheric radiation lookup table is queried to obtain the transmittance and radiation of the atmospheric path traversed by the reverse ray; if the reverse ray escapes to the sky direction, the sky radiance is calculated based on the current reverse ray direction and the sky radiation signal is accumulated.

[0013] Furthermore, the calculation of the sky radiance includes the calculation of two parts: the basic radiance of the sky and the angular distribution. The baseline radiance of the sky is calculated from the solar irradiance at the top of the atmosphere, atmospheric transmittance, and surface albedo; the angular distribution is determined based on the CIE standard cloudy sky model, the model expression of which is: L(θ) = L_zenith×(1+2cosθ) / 3, In the formula, L(θ) is the sky radiance in the direction of the zenith angle θ, and L_zenith is the sky radiance in the zenith direction; a nearest neighbor enhancement factor in the direction of the sun is superimposed on the angular distribution to simulate the scattering enhancement effect of the solar corona.

[0014] Furthermore, in step two, the importance sampling specifically involves sampling the emission direction of the reverse ray based on the bidirectional reflection distribution function of the hit triangle surface element, so that the direction of the reverse ray that contributes more is sampled with a higher probability. The throughput of the reverse ray is updated based on the bidirectional reflection distribution function value, the cosine of the emission direction, and the sampling probability density function. Other sampling strategies are used to calculate the importance sampling weight to achieve weighted correction of the Monte Carlo estimate, so as to ensure the unbiasedness or low variance of the estimate. The Russian roulette game controls the ejection and termination of the reverse ray. Specifically, when the reverse ray continues to eject, a new ray is emitted from the point of impact along the emission direction obtained from the sampling, and the intersection calculation and different radiation transmission processing are repeated until the ray terminates.

[0015] Furthermore, the bidirectional reflection distribution function of the hit triangular surface element includes three modes to enable different components in the same scene to use different bidirectional reflection distribution functions, specifically including: Mode 1, the Lambertian diffuse reflection model, with the corresponding bidirectional reflection distribution function f r for: f r = ρ / π, Where ρ is the surface hemispherical reflectivity, which is the hemispherical uniform sampling with cosθ weighting, and the probability density function pdf = cosθ / π; Mode 2, Improved Von Model, with corresponding bidirectional reflection distribution function f r for: f r = ρ d / π+ρ s ×(n+2) / (2π)×cos n α, Where, ρ d ρ is the diffuse reflectance coefficient. s denoted as the specular reflection coefficient, n as the Phong exponent, and α as the angle between the outgoing direction and the ideal specular reflection direction. During sampling, hemispherical uniform sampling or Phong specular lobe sampling is selected with a preset probability, and balance heuristic is used to allocate multiple importance sampling weights. Mode 3 is the Stanford-Robertson physical two-way reflection distribution function model, with the corresponding two-way reflection distribution function f. r for: f r = (ε+ρ d )×(b / ((cosθ i +cosθ o ) e )), Where ε is the specular reflection coefficient, ρ d θ is the diffuse reflection coefficient, b is the scale parameter, e is the roughness exponent, and θ is the surface roughness index. i Let θ be the angle of incidence of the ray. o The ray exit angle is used; during sampling, the sampling distribution width is controlled according to the roughness parameter e, and powerheuristic is used to allocate the sampling weights for multiple importance samples; The three bidirectional reflection distribution functions f r The unified normalization convention model DRF is: DRF = π×f r .

[0016] Furthermore, the execution method of the Russian roulette is as follows: setting a maximum ejection depth, a minimum ejection depth, and a Russian roulette probability; when the number of ejections of the reverse ray is less than the minimum ejection depth, the reverse ray continues to eject; when the number of ejections of the reverse ray is greater than or equal to the minimum ejection depth but less than the maximum ejection depth, the Russian roulette probability is used to determine whether the reverse ray continues to eject; when the number of ejections reaches the maximum ejection depth, the reverse ray terminates.

[0017] Furthermore, the signal components include six components, which are defined as follows: Component 1 is the direct reflection signal from the laser target: the radiation contribution from the laser beam that is reflected from the target surface after atmospheric attenuation and enters the detector. Component two is the laser background reflection signal: the radiation contribution of the laser after atmospheric attenuation to the background surface, which is reflected into the detector. Component three is the sky path radiation signal: when the ray escapes to the sky direction, the sky background radiance directly enters the radiation contribution of the detector; Component four is the solar surface reflection signal: the radiation contribution of sunlight after direct sunlight shines on the target surface or background surface and enters the detector after reflection; Component 5 is the atmospheric backscattering signal: the radiation contribution of atmospheric aerosol backscattering directly into the detector, which is used as a noise monitoring quantity; Component six is ​​the sum of radiation contributions from multiple ejection signals: the total radiation contribution of reverse ray paths with an ejection number greater than or equal to 1; each reverse ray path is added to component six while also being added to component one or component two based on the medium type attribute of the surface hit. The radiation contributions of components one, two, three, four, and six constitute the total signal output in the simulation, ensuring the physical completeness of the signal decomposition. Each of the six components is output as a separate binary data file.

[0018] Furthermore, in step four, the specific process of time-domain full waveform synthesis is as follows: preset the time bin width, Gaussian pulse full width at half maximum and maximum detection distance, and calculate the number of time bins; map all inverse Monte Carlo samples of each pixel to the corresponding time bin according to the flight time, and synthesize the time-domain full waveform data of each pixel through Gaussian pulse convolution.

[0019] Furthermore, the triangular surface grid of the imaging target is a CFD grid for an aircraft, and the corresponding spatial acceleration structure is a uniform grid acceleration structure or a hierarchical enclosing volume acceleration structure.

[0020] The beneficial effects of this invention are: 1. High-fidelity multi-source coupling simulation: Through the inverse Monte Carlo path tracing framework, the radiation transmission processes of the target, background, atmosphere and sun are coupled simultaneously in a unified physical model, which can accurately simulate complex physical effects such as multiple reflections between the target and background and the influence of atmospheric scattering on the target signal.

[0021] 2. Flexible support for bidirectional reflection distribution function materials: Provides multi-level material description capabilities from simple Lambertian models to physically bidirectional reflection distribution function models, combined with bidirectional reflection distribution function importance sampling and MIS, achieving a balance between computational efficiency and sampling quality.

[0022] 3. Rich output data dimensions: In addition to the conventional intensity and range images, it provides quantitative decomposition of six-source signals and full waveform synthesis output, providing comprehensive data support for subsequent physical analysis, feature extraction and algorithm development.

[0023] 4. High-efficiency computation driven by media type: By using a differentiated processing strategy driven by media type tags, different computational models with different complexities are adopted for different material paths, thereby optimizing computational efficiency while ensuring physical accuracy.

[0024] 5. High-fidelity coupling of atmosphere and terrain: By integrating the MODTRAN atmospheric radiative transfer model and CopernicusDEM / ESA WorldCover terrain data, physical consistency coupling of atmospheric transmittance, aerosol scattering and terrain reflection characteristics is achieved, supporting simulations under arbitrary geographical locations and atmospheric conditions. Attached Figure Description

[0025] Figure 1 This is the overall flowchart of the lidar imaging simulation method of the present invention; Figure 2 This is a quantitative decomposition diagram of the six source signal components in a specific embodiment of the present invention; Figure 3 This is the intensity image output from the simulation in Specific Example 1 of the present invention. Detailed Implementation

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

[0027] To achieve the above objectives, the present invention provides the following specific embodiments: For ease of understanding and implementation, several conventions in the embodiments are explained as follows. RMCM stands for Reverse Monte Carlo Method (RMCM) path tracing; "hit" refers to the intersection of the reverse ray with a triangular element in the scene; "ejection" refers to the behavior of the reverse ray continuing to track at a hit point according to the sampled emission direction; "ejection count" refers to the number of times the reverse ray hits the scene surface again after its first hit from the detector pixel: if the ray terminates its path after the first hit (reaching the laser source direction, escaping to the sky direction, or meeting the termination condition), its ejection count is recorded as 0; each time the above-mentioned second hit occurs, the ejection count is incremented by 1; "ejection order" refers to the sequential numbering of each hit point traversed by the reverse ray along the propagation direction.

[0028] In multi-source signal decomposition, the direct reflection signal of the laser target (component 1) and the laser background reflection signal (component 2) correspond to the direct path (ejection count 0) after the first hit of the reverse ray on the target surface or background surface. All multiple ejection paths with an ejection count greater than or equal to 1 are uniformly accumulated into component 6, and are assigned according to the medium type attribute (target surface or background surface) corresponding to the last surface reflection of the path (i.e., the first hit of the reverse ray on the surface), forming a target class / background class multiple ejection sub-statistic for physical path analysis. Multiple ejection paths are not repeatedly counted in component 1 or component 2, thus ensuring that the sum of component 1, component 2, component 3, component 4 and component 6 equals the total signal output by the simulation, avoiding double counting.

[0029] The ray throughput W is used to record the accumulated radiation weight of a reverse ray path in each ejection, with an initial value of W0=1. When the next ejection direction is sampled for importance according to the bidirectional reflection distribution function fr at the hit point and the tracking continues, the throughput is updated according to W=W×fr×|cosθ| / p(ω), where θ is the angle between the ejection direction and the surface normal of the hit point, and p(ω) is the sampling probability density of the direction. After the ray terminates, the source radiation corresponding to the termination event is multiplied by the throughput and accumulated to the corresponding signal component. The source radiation is taken as the irradiance of the laser after atmospheric attenuation, the direct solar irradiance, the sky background irradiance, or the atmospheric path radiation according to the light source type. The pixel received radiation value is the arithmetic mean of the contributions of all sampling paths of the pixel, and is converted by coefficients such as pixel solid angle, receiving aperture area, and optical efficiency.

[0030] Example 1: A simulation method for multi-source coupled imaging of lidar based on inverse Monte Carlo path tracing, comprising the following steps: S01. Construct a triangular mesh model of the imaging target and the corresponding spatial acceleration structure, and label the medium type attributes of the triangular mesh to form a three-dimensional scene that is consistent with the imaging target and the background and can be ray-traced. The three-dimensional scene includes a background scene constructed from a terrain digital elevation model and land cover classification data, wherein the triangular facets are mapped using a bidirectional reflection distribution function based on the land cover classification data.

[0031] S02. Based on the simulation requirements of the imaging target, set the atmospheric mode to vacuum mode, Malkmus statistical narrowband model mode, or complete atmospheric radiative transfer model mode. The medium types include target surface, background surface, atmospheric path, and sky direction. When the atmospheric model is the complete atmospheric radiative transfer model, the complete atmospheric radiative transfer model automatically calls the external atmospheric radiative transfer process, extracts atmospheric transmittance, path radiation, and sky background radiation data, and constructs a multidimensional atmospheric radiation lookup table. S03. Based on the inverse Monte Carlo method, each pixel of the imaging target in the detector imaging plane is sampled. The sampling of any pixel is based on the position of any pixel, and the initial ray is generated by randomly reversing the ray direction within the detector's field of view. Furthermore, by utilizing a space acceleration structure, the intersection calculation between the reverse ray and the triangular surface element is performed in a unified scene containing the imaging target and the background. The location of the hit point, the geometric normal direction of the surface of the hit triangular surface element, and the medium type attribute of the triangular surface element are obtained. Based on the medium type attribute, different reverse ray radiative transmission processing is performed to obtain the transmittance and radiation of the atmospheric path through which the reverse ray passes. If the reverse ray hits the imaging target surface, the bidirectional reflectance distribution function of the imaging target surface is sampled and the ray throughput is updated; if the reverse ray hits the background surface, the bidirectional reflectance distribution function of the background surface is sampled and the background signal is accumulated; if the reverse ray passes through the atmospheric path, the multidimensional atmospheric radiation lookup table is queried to obtain the transmittance and radiation of the atmospheric path traversed by the reverse ray; if the reverse ray escapes to the sky direction, the sky radiance is calculated based on the current reverse ray direction and the sky radiation signal is accumulated. The calculation of sky radiance includes the calculation of two parts: the basic sky radiance and the angular distribution. The baseline radiance of the sky is calculated from the solar irradiance at the top of the atmosphere, atmospheric transmittance, and surface albedo; the angular distribution is determined based on the CIE standard cloudy sky model, the model expression of which is: L(θ) = L_zenith×(1+2cosθ) / 3, In the formula, L(θ) is the sky radiance in the direction of the zenith angle θ, and L_zenith is the sky radiance in the zenith direction; a nearest neighbor enhancement factor in the direction of the sun is superimposed on the angular distribution to simulate the scattering enhancement effect of the solar corona.

[0032] S04. At each hit point location in the reverse ray tracing, the importance of the reflection contribution distribution determined by the two-way reflection distribution function of the hit triangle element surface in the next emission direction is sampled to ensure the unbiasedness of the radiation contribution estimation. Importance sampling specifically involves sampling the emission direction of the reverse ray based on the bidirectional reflection distribution function of the hit triangle surface element, so that the direction of the reverse ray that contributes more is sampled with a higher probability. The throughput of the reverse ray is updated based on the bidirectional reflection distribution function value, the cosine of the emission direction, and the sampling probability density function. Other sampling strategies are used to calculate the importance sampling weight, thereby achieving a weighted correction of the Monte Carlo estimate to ensure the unbiasedness or low variance of the estimate. The bidirectional reflection distribution function for hitting the surface of a triangular element includes three modes, used to allow different parts in the same scene to use different bidirectional reflection distribution functions, specifically including: Mode 1, the Lambertian diffuse reflection model, with the corresponding bidirectional reflection distribution function f r for: f r = ρ / π, Where ρ is the surface hemispherical reflectivity, which is the hemispherical uniform sampling with cosθ weighting, and the probability density function pdf = cosθ / π; Mode 2, Improved Von Model, with corresponding bidirectional reflection distribution function f r for: f r = ρ d / π+ρ s ×(n+2) / (2π)×cos n α, Where, ρ d ρ is the diffuse reflectance coefficient. s denoted as the specular reflection coefficient, n as the Phong exponent, and α as the angle between the outgoing direction and the ideal specular reflection direction. During sampling, hemispherical uniform sampling or Phong specular lobe sampling is selected with a preset probability, and balance heuristic is used to allocate multiple importance sampling weights. Mode 3 is the Stanford-Robertson physical two-way reflection distribution function model, with the corresponding two-way reflection distribution function f. r for: f r = (ε+ρ d)×(b / ((cosθ i +cosθ o ) e )), Where ε is the specular reflection coefficient, ρ d θ is the diffuse reflection coefficient, b is the scale parameter, e is the roughness exponent, and θ is the surface roughness index. i Let θ be the angle of incidence of the ray. o The ray exit angle is used; during sampling, the sampling distribution width is controlled according to the roughness parameter e, and powerheuristic is used to allocate the sampling weights for multiple importance samples; The three bidirectional reflection distribution functions f r The unified normalization convention model DRF is: DRF = π×f r .

[0033] S05. The Russian roulette wheel is used to control the ejection and termination of the reverse ray. When the reverse ray continues to eject, a new ray is emitted from the hit point along the emission direction obtained by sampling, and the intersection calculation and different radiation transmission processing are repeated until the ray terminates. After the reverse ray terminates, the radiation contribution of the corresponding pixel is obtained. The execution method of Russian Roulette is as follows: set the maximum ejection depth, the minimum ejection depth, and the Russian Roulette probability; when the number of ejections of the reverse ray is less than the minimum ejection depth, the reverse ray continues to eject; when the number of ejections of the reverse ray is greater than or equal to the minimum ejection depth but less than the maximum ejection depth, the Russian Roulette probability is used to determine whether the reverse ray continues to eject; when the number of ejections reaches the maximum ejection depth, the reverse ray terminates.

[0034] S06. During the reverse ray path tracing process from S03 to S05, the radiation contribution is classified and accumulated to the corresponding signal component according to the medium type, light source type and ejection order of the hit point location, so as to achieve multi-source signal decomposition. The signal components consist of six components, which are defined as follows: Component 1 is the direct reflection signal from the laser target: the radiation contribution from the laser beam that is reflected from the target surface after atmospheric attenuation and enters the detector. Component two is the laser background reflection signal: the radiation contribution of the laser after atmospheric attenuation to the background surface, which is reflected into the detector. Component three is the sky path radiation signal: when the ray escapes to the sky direction, the sky background radiance directly enters the radiation contribution of the detector; Component four is the solar surface reflection signal: the radiation contribution of sunlight after direct sunlight shines on the target surface or background surface and enters the detector after reflection; Component 5 is the atmospheric backscattering signal: the radiation contribution of atmospheric aerosol backscattering directly into the detector, which is used as a noise monitoring quantity; Component six is ​​the sum of radiation contributions from multiple ejection signals: the total radiation contribution of reverse ray paths with an ejection number greater than or equal to 1; each reverse ray path is added to component six while also being added to component one or component two based on the medium type attribute of the surface hit. The radiation contributions of components one, two, three, four, and six constitute the total signal output in the simulation, ensuring the physical completeness of the signal decomposition. Each of the six components is output as a separate binary data file.

[0035] S07. Based on the time-of-flight distribution of the inverse ray of each pixel, synthesize the full temporal waveform data of the corresponding pixel through time bin mapping and Gaussian pulse convolution, and output the simulation results including intensity image, distance image, time-of-flight image, signal decomposition file and full waveform data file.

[0036] The specific process of time-domain full waveform synthesis is as follows: preset the time bin width, Gaussian pulse full width at half maximum and maximum detection distance, and calculate the number of time bins; map all inverse Monte Carlo samples of each pixel to the corresponding time bin according to the flight time, and synthesize the time-domain full waveform data of each pixel through Gaussian pulse convolution.

[0037] Example 2: Same as Example 1, except that the triangular surface grid of the imaging target is a spacecraft CFD grid, and the corresponding spatial acceleration structure is a uniform grid acceleration structure or a hierarchical enclosing volume acceleration structure.

[0038] like Figures 1-3 As shown, to further illustrate the technical solution and technical effects of the present invention, the following specific examples are provided: Specific example 1: Reference Figure 1 This specific example provides a simulation method for multi-source coupled lidar imaging based on inverse Monte Carlo (RMCM) path tracing. The example uses a flight model (FLUENT CFD mesh, 152,862 triangular elements) under US Standard 1976 atmospheric conditions (visibility 23 km, rural aerosols, CO2 415 ppmv) for 1064 nm lidar imaging simulation.

[0039] S1. Scene Construction: Import the FLUENT CFD mesh file of the aircraft, containing 152,862 triangular elements, including nine components: fuselage, air intake, wings, control surfaces, control surface wingtips, tail, tail section, air intake wingtips, and wingtips. Construct a uniform mesh to accelerate the structure. Import the Copernicus 30m digital elevation model (30m resolution) and ESA WorldCover 10m (2021) land cover classification data covering the target area, with the region center located at 34.22°N, 108.8°E. Automatically generate a terrain background triangular mesh and automatically map the bidirectional reflectance distribution function parameters of each triangular element according to the ESA WorldCover classification system.

[0040] Label the medium type for each triangular element: the target surface element is labeled TARGET and associated with the reflectivity of its components (aircraft body ρ = 0.30, air intake ρ = 0.55, wing ρ = 0.30, control surface ρ = 0.35, control surface wingtip ρ = 0.25, tail fin ρ = 0.25, tail section ρ = 0.40, air intake wingtip ρ = 0.28, wingtip ρ = 0.25), and the terrain surface element is labeled BACKGROUND.

[0041] S2. Atmospheric model initialization: Set atm mode = 2 to enable MODTRAN mode. Select MODEL = 5 (US Standard Atmosphere 1976), VIS = 23km, IHAZE = 1 (rural aerosols), CO2 = 415ppmv, wavelength = 1.064μm.

[0042] The system automatically invokes the MODTRAN process to generate the tape5 input file and parse the tape6 output file, extracting the atmospheric transmittance sequence and path radiation sequence, and constructing an atmospheric radiation lookup table. At a sensor height of 0.02 km and a slant distance of 0.02 km, the path transmittance is approximately 0.9985.

[0043] S3, RMCM Main Simulation Loop: Detector parameters are set as follows: image resolution 480×240 pixels, field of view 60° horizontal × 45° vertical, sensor located at (0, 0, 20) m, pointing towards (4, 0, 0) m. Laser parameters: wavelength 1064nm, transmit power 1000W, receiver aperture 0.5m, optical efficiency 0.8. RMCM settings: samples per pixel spp = 64, maximum ejection depth max depth = 8, Russian roulette probability rr prob = 0.33, minimum ejection depth rr min depth = 3. In this specific example, the target bidirectional reflectance distribution function uses mode = 0 (Lambertian diffuse reflection), and each component is allocated according to the actual reflectivity.

[0044] The sun / sky was set to enabled on June 15th, at 4:00 AM UTC (12:00 PM Beijing time), with direct solar radiation (DNI) of 800 W / m². Parallel computation was performed using 36 OpenMP threads, with a total simulation time of approximately 314 seconds, processing 115,200 pixels.

[0045] S4. Six-Source Signal Decomposition: During RMCM tracing, the radiation contribution of each ray path is automatically classified and accumulated into six signal source components based on its ejection trajectory and source type. (See reference...) Figure 2 The quantitative decomposition results of the six source signal components show the radiation contribution distribution of each component.

[0046] In this specific example, under the condition of short-range detection at noon, the solar-background reflection component accounts for 99.5% of the total intensity, the laser-target direct reflection component accounts for approximately 0.48%, and the multibounce component accounts for approximately 0.0065%, reflecting the multiple reflection effects between the target surface and the terrain or between target components. The sum of all radiation components is equivalent to the total radiation intensity, ensuring the completeness of signal decomposition.

[0047] S5. Full Waveform Synthesis: In this specific example, waveform synthesis is not enabled (waveform enabled = 0), but the method framework supports waveform synthesis—preset bin width = 0.5ns, FWHM = 5ns, maximum detection distance 3000m, and number of bins 4096. When enabled, each sample is mapped to the corresponding time bin according to the time of flight (ToF), and the temporal full waveform data of each pixel is synthesized through Gaussian pulse convolution.

[0048] S6. Output Results: The simulation output includes 480×240 pixel intensity image, range image, time-of-flight image, six-source signal decomposition file, and full waveform data file. The intensity image output from the simulation is shown below. Figure 3As shown, this reflects the distribution of radiated power received by each pixel on the detector plane.

[0049] Specific example 2: This specific example, building upon the method in Specific Example 1, focuses on illustrating the specific implementation of different bidirectional reflection distribution function models. The method of this invention supports three bidirectional reflection distribution function models to meet the needs of various application scenarios.

[0050] Bidirectional reflectance distribution function mode = 0 (Lambertian diffuse reflectance): Bidirectional reflectance distribution function f r = ρ / π, where ρ is the hemispherical reflectivity of the surface. In this specific example, the reflectivity of the aircraft fuselage is ρ = 0.30, the reflectivity of the air intake is ρ = 0.55, the reflectivity of the wing is ρ = 0.30, and the reflectivity of the control surface is ρ = 0.35. The sampling strategy is cosθ-weighted uniform random sampling in the hemispherical space, with the probability density function pdf = cosθ / π. This algorithm is simple, computationally efficient, and suitable for materials with predominantly diffuse reflection; it is the mode currently used in experimental operations.

[0051] Two-way reflection distribution function mode mode = 1 (Modified Phong model): Two-way reflection distribution function f r = ρ d / π+ρ s ×(n+2) / (2π)×cos n α, where ρ d ρ is the diffuse reflectance coefficient. s Let be the specular reflection coefficient, n be the Phong exponent, and α be the angle between the outgoing direction and the ideal specular reflection direction. During sampling, hemispherical uniform sampling (diffuse component) is selected with a 50% probability, and sampling within the Phong specular lobe (spectral component) is selected with a 50% probability. MIS multi-importance sampling employs a balance heuristic weighting strategy. This mode is suitable for metal or painted surfaces with medium gloss.

[0052] Two-way reflection distribution function mode mode = 2 (Stanford-Robertson physical two-way reflection distribution function model): Two-way reflection distribution function f r = (ε+ρ d )×(b / ((cosθ i +cosθ o ) e It includes four physical parameters—ε is the specular reflection coefficient, ρ is the reflectance coefficient, and ρ is the reflectance coefficient. d Let be the diffuse reflection coefficient, b be the scale parameter, and e be the roughness exponent. Directional reflection factor DRF = π×(ε+ρ) d)×(b / ((cosθ i +cosθ o ) e During sampling, the sampling distribution width is controlled according to the roughness parameter e: the smaller the e value, the smoother the surface, and the more concentrated the sampling direction is near the specular reflection direction.

[0053] Taking this aircraft model as an example, referring to the Stanford-Robertson material parameters—aircraft fuselage configuration (ε = 0.12, ρ...) d = 0.025, b = 0.060, e = 0.0250), wing configuration (ε = 0.10, ρ d = 0.020, b = 0.080, e = 0.0150), intake configuration (ε = 0.08, ρ d = 0.015, b = 0.040, e = 0.0100). The three bidirectional reflection distribution function models can be flexibly selected according to the material characteristics of the target surface, supporting the use of different bidirectional reflection distribution function models for different parts in the same scene.

[0054] Specific example 3: This specific example details the six-source signal decomposition mechanism and its implementation principle.

[0055] Six-source signal decomposition is one of the core features of this invention. In traditional lidar imaging simulation, only the total intensity image is usually output, which cannot distinguish the radiation contributions from different sources. This invention automatically classifies and accumulates the radiation contributions generated by each ray path at each ejection point during RMCM path tracing. The classification criteria include combinations of light source type (laser source, solar source, sky source, atmospheric path radiation) and radiation transmission path (direct target hit, background hit, atmospheric scattering, multiple ejections between targets, etc.).

[0056] The specific definitions of the six signal components are as follows: Component 1 (target, laser target direct reflection signal) is the radiation contribution of the laser beam after being reflected from the target surface and entering the detector. Its value depends on whether the medium type property of the hit surface is the target surface. Component 2 (background, laser background reflection signal) is the radiation contribution of the laser beam after being reflected by the background surface and entering the detector. Its value depends on whether the medium type property of the hit surface is the background surface. Component 3 (sky, sky path radiation signal) is the radiation contribution of the sky radiance that directly enters the detector when the ray escapes to the sky direction; Component 4 (sun, solar surface reflection signal) is the radiation contribution of sunlight reflected into the detector after direct sunlight shines on the target surface or background surface; Component 5 (backscatter, atmospheric backscatter signal) is the radiation contribution of atmospheric aerosol backscatter directly into the detector (in RMCM path tracing, it is monitored as a noise term rather than a signal component and is added together). Component 6 (multibounce signal) is the radiation contribution of ray paths with a number of ejections greater than or equal to 1. Each ray path is simultaneously classified as either Component 1 or Component 2 depending on the medium type of the surface it hits.

[0057] The radiation contributions of components one, two, three, four, and six constitute the total signal output by the simulation, and their sum ensures the physical completeness of the signal decomposition.

[0058] The correspondence between the above six components and the physical radiation transmission path is as follows: the laser light is attenuated by the atmosphere and then reflects off the target surface into the detector (physical path, corresponding to component one); the laser light is attenuated by the atmosphere and then reflects off the background surface into the detector (physical path, corresponding to component two); the background radiation from the sky directly enters the detector (physical path, corresponding to component three); sunlight is reflected off the target or background surface into the detector (physical path, corresponding to component four); atmospheric backscattering enters the detector (physical path, corresponding to component five); multiple reflections between the target and the background, or between different components of the target (physical effect, recorded in component six, and depending on whether the final hit surface belongs to component one or component two).

[0059] During RMCM path tracing, the system automatically adds the radiation contribution to the corresponding signal source component each time it records the radiation contribution, based on the medium type, light source type, and ejection order of the current ejection point. This online decomposition method requires no additional post-processing steps; the calculation results are available immediately after tracing is complete.

[0060] Reference Figure 2 This specific example demonstrates the quantitative decomposition results of a six-source signal based on actual simulation data. The total simulated intensity comprises six components: target component, background component, sky component, solar component, atmospheric backscattering component, and multiple ejection component. In the noon detection scenario of this specific example, the solar-background reflection component constitutes the main part of the total signal, which is consistent with the physical expectation of high direct solar irradiance. The multiple ejection component reflects the multiple reflection effects between the target surface and the terrain or between target components, and its existence confirms that the method of this invention can effectively capture the contribution of multiple ejection paths to the overall radiative transfer.

[0061] The decomposition of the six-source signal provides important information for subsequent physical analysis. By analyzing the proportion of each component, the relative importance of each radiation source under different atmospheric conditions, detection angles, and solar geometry can be assessed. In the development of target detection and recognition algorithms, the information of each component can be used for signal separation and feature extraction. In multi-sensor fusion applications, each component can serve as prior knowledge of the responses of different sensors.

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

Claims

1. A simulation method for multi-source coupled imaging of lidar based on inverse Monte Carlo path tracing, characterized in that, Includes the following steps: Step 1: Construct a triangular mesh model of the imaging target and its corresponding spatial acceleration structure, and label the medium type attributes of the triangular mesh to form a 3D scene that is consistent with the background and can be ray-traced; at the same time, according to the simulation requirements of the imaging target, set the atmospheric mode to vacuum mode, Malkmus statistical narrowband model mode, or complete atmospheric radiative transfer model mode. The medium types include target surface, background surface, atmospheric path, and sky direction. When the atmospheric model is the complete atmospheric radiative transfer model, the complete atmospheric radiative transfer model automatically calls the external atmospheric radiative transfer process, extracts atmospheric transmittance, path radiation, and sky background radiation data, and constructs a multidimensional atmospheric radiation lookup table. Step 2: Based on the inverse Monte Carlo method, sample each pixel of the imaging target in the detector imaging plane. The sampling of any pixel starts from any pixel position and randomly reverses the ray direction within the detector's field of view to generate an initial ray. Furthermore, by utilizing a space acceleration structure, the intersection calculation between the reverse ray and the triangular surface element is performed in a unified scene containing the imaging target and the background. The location of the hit point, the geometric normal direction of the surface of the hit triangular surface element, and the medium type attribute of the triangular surface element are obtained. Based on the medium type attribute, different reverse ray radiative transmission processing is performed to obtain the transmittance and radiation of the atmospheric path through which the reverse ray passes. Next, at each hit point location of the reverse ray tracing, the importance of the reflection contribution distribution determined by the bidirectional reflection distribution function of the hit triangle surface element is sampled for the next emission direction to ensure the unbiasedness of the radiation contribution estimation. Finally, the ejection and termination of the reverse ray are controlled by Russian roulette, and the radiation contribution of the corresponding pixel is obtained after the reverse ray terminates. Step 3: During the reverse ray path tracing process in Step 2, the radiation contribution is classified and accumulated to the corresponding signal component according to the medium type, light source type and ejection order of the hit point location, so as to achieve multi-source signal decomposition. Step 4: Based on the time-of-flight distribution of the inverse ray of each pixel, synthesize the full temporal waveform data of the corresponding pixel through time bin mapping and Gaussian pulse convolution, and output the simulation results including intensity image, distance image, time-of-flight image, signal decomposition file and full waveform data file.

2. The LiDAR multi-source coupled imaging simulation method for inverse Monte Carlo path tracing as described in claim 1, characterized in that, In step one, the three-dimensional scene includes a background scene constructed from a terrain digital elevation model and land cover classification data, wherein the triangular surface element is mapped to a bidirectional reflection distribution function based on the land cover classification data.

3. The LiDAR multi-source coupled imaging simulation method for inverse Monte Carlo path tracing as described in claim 1, characterized in that, In step two, the specific steps of performing different reverse ray radiative transfer processes based on the medium type properties are as follows: If the reverse ray hits the imaging target surface, the bidirectional reflectance distribution function of the imaging target surface is sampled and the ray throughput is updated; if the reverse ray hits the background surface, the bidirectional reflectance distribution function of the background surface is sampled and the background signal is accumulated. If the reverse ray passes through the atmospheric path, then the transmittance and radiation of the atmospheric path through which the reverse ray passes can be obtained by consulting the multidimensional atmospheric radiation lookup table. If the reverse ray escapes to the sky, the sky radiance is calculated based on the current reverse ray direction, and the sky radiation signal is accumulated.

4. The LiDAR multi-source coupled imaging simulation method for inverse Monte Carlo path tracing as described in claim 3, characterized in that, The calculation of sky radiance includes the calculation of two parts: the basic sky radiance and the angular distribution. The baseline radiance of the sky is calculated from the solar irradiance at the top of the atmosphere, atmospheric transmittance, and surface albedo; the angular distribution is determined based on the CIE standard cloudy sky model, the model expression of which is: L(θ) = L_zenith×(1+2cosθ) / 3, In the formula, L(θ) is the sky radiance in the direction of the zenith angle θ, and L_zenith is the sky radiance in the zenith direction; Nearest neighbor enhancement factors in the solar direction are superimposed on the angular distribution to simulate the scattering enhancement effect of the solar corona.

5. The LiDAR multi-source coupled imaging simulation method for inverse Monte Carlo path tracing as described in claim 1, characterized in that, In step two, the importance sampling specifically involves sampling the emission direction of the reverse ray based on the bidirectional reflection distribution function of the surface of the hit triangle element, so that the direction of the reverse ray that contributes more is sampled with a higher probability. The throughput of the reverse ray is updated based on the bidirectional reflection distribution function value, the cosine of the emission direction, and the sampling probability density function. Other sampling strategies are used to calculate the importance sampling weight, thereby achieving a weighted correction of the Monte Carlo estimate to ensure the unbiasedness or low variance of the estimate. The Russian roulette game controls the ejection and termination of the reverse ray. Specifically, when the reverse ray continues to eject, a new ray is emitted from the point of impact along the emission direction obtained from the sampling, and the intersection calculation and different radiation transmission processing are repeated until the ray terminates.

6. The LiDAR multi-source coupled imaging simulation method for inverse Monte Carlo path tracing as described in claim 5, characterized in that, The bidirectional reflection distribution function of the hit triangular surface element includes three modes to enable different components in the same scene to use different bidirectional reflection distribution functions, specifically including: Model one, Lambertian diffuse reflection model, the corresponding bidirectional reflection distribution function f r is: f r = p / π, Where ρ is the surface hemispherical reflectivity, which is the hemispherical uniform sampling with cosθ weighting, and the probability density function pdf = cosθ / π; Mode 2, Improved Von Model, with corresponding bidirectional reflection distribution function f r for: f r = p d / p+r s ×(n+2) / (2π)×cos n a, Where, ρ d ρ is the diffuse reflectance coefficient. s denoted as the specular reflection coefficient, n as the Phong exponent, and α as the angle between the outgoing direction and the ideal specular reflection direction. During sampling, hemispherical uniform sampling or Phong specular lobe sampling is selected with a preset probability, and balance heuristic is used to allocate multiple importance sampling weights. Mode 3 is the Stanford-Robertson physical two-way reflection distribution function model, with the corresponding two-way reflection distribution function f. r for: f r = (e+r) d )×(b / ((cosθ i +cosθ o ) e )), Where ε is the specular reflection coefficient, ρ d θ is the diffuse reflection coefficient, b is the scale parameter, e is the roughness exponent, and θ is the surface roughness index. i Let θ be the angle of incidence of the ray. o The ray exit angle is used; during sampling, the sampling distribution width is controlled according to the roughness parameter e, and power heuristic is used to allocate the sampling weights for multiple importance samples; The three bidirectional reflection distribution functions f r The unified normalization convention model DRF is: DRF = π×f r 。 7. The LiDAR multi-source coupled imaging simulation method for inverse Monte Carlo path tracing as described in claim 5, characterized in that, The execution method of the Russian roulette is as follows: set the maximum ejection depth, the minimum ejection depth, and the Russian roulette probability; when the number of ejections of the reverse ray is less than the minimum ejection depth, the reverse ray continues to eject; when the number of ejections of the reverse ray is greater than or equal to the minimum ejection depth but less than the maximum ejection depth, the Russian roulette probability is used to determine whether the reverse ray continues to eject; when the number of ejections reaches the maximum ejection depth, the reverse ray terminates.

8. The LiDAR multi-source coupled imaging simulation method for inverse Monte Carlo path tracing as described in claim 1, characterized in that, The signal components include six components, which are defined as follows: Component 1 is the direct reflection signal from the laser target: the radiation contribution from the laser beam that is reflected from the target surface after atmospheric attenuation and enters the detector. Component two is the laser background reflection signal: the radiation contribution of the laser after atmospheric attenuation to the background surface, which is reflected into the detector. Component three is the sky path radiation signal: when the ray escapes to the sky direction, the sky background radiance directly enters the radiation contribution of the detector; Component four is the solar surface reflection signal: the radiation contribution of sunlight after direct sunlight shines on the target surface or background surface and enters the detector after reflection; Component 5 is the atmospheric backscattering signal: the radiation contribution of atmospheric aerosol backscattering directly into the detector, which is used as a noise monitoring quantity; Component six is ​​the sum of radiation contributions from multiple ejection signals: the total radiation contribution of reverse ray paths with an ejection number greater than or equal to 1; each reverse ray path is added to component six while also being added to component one or component two based on the medium type attribute of the surface hit. The radiation contributions of components one, two, three, four, and six constitute the total signal output in the simulation, ensuring the physical completeness of the signal decomposition. Each of the six components is output as a separate binary data file.

9. The LiDAR multi-source coupled imaging simulation method for inverse Monte Carlo path tracing as described in claim 1, characterized in that, In step four, the specific process of time-domain full waveform synthesis is as follows: preset the time bin width, Gaussian pulse full width at half maximum and maximum detection distance, and calculate the number of time bins; map all inverse Monte Carlo samples of each pixel to the corresponding time bin according to the flight time, and synthesize the time-domain full waveform data of each pixel through Gaussian pulse convolution.

10. The lidar multi-source coupled imaging simulation method for inverse Monte Carlo path tracing as described in any one of claims 1-9, characterized in that, The triangular surface grid of the imaging target is the CFD grid of the aircraft, and the corresponding spatial acceleration structure is a uniform grid acceleration structure or a hierarchical enclosing volume acceleration structure.

Images