Optical surveying and mapping satellite constellation three-dimensional imaging simulation method

Abstract

The invention relates to an optical surveying and mapping satellite constellation three-dimensional imaging simulation method. The method comprises the following steps: constructing a live-action three-dimensional model based on an aviation oblique photogrammetry technology and a three-dimensional illumination rendering mode; the method comprises the following steps: acquiring preset satellite imaging data, constructing an image space vector and an object space vector, converting the image space vector and the object space vector to a body coordinate system, and constructing an imaging geometric model considering a satellite platform error parameter, a satellite camera error parameter and an atmospheric error parameter; traversing all pixels of the imaging plane in the live-action three-dimensional model, substituting image space coordinates corresponding to the pixels into the imaging geometric model, and generating a simulation image; a simulation image is previewed and output in real time through an interface, and satellite auxiliary data is output according to the imaging geometric model during output. Compared with the prior art, the method has the advantages that simulation imaging of ground objects with vertical facades such as buildings and steep high mountains can be achieved, and application of three-dimensional modeling is met.

Description

Technical Field

[0001] This invention relates to the field of satellite imaging simulation, and in particular to a three-dimensional imaging simulation method for optical mapping satellite constellations. Background Technology

[0002] Traditional mapping satellites employ a multi-line array system, maintaining a stable Earth-attitude orbit. To avoid impacting positioning accuracy, the satellite's maneuvering angle during on-orbit mapping and imaging must not exceed 10°, significantly constraining its operational flexibility and mapping efficiency. The mission of next-generation optical stereo mapping satellites has shifted from large-scale global mapping to precise mapping of key areas, and mapping products are transitioning from traditional 2D and 2.5D products to 3D reality products. This necessitates faster and more efficient mapping capabilities, including the ability to create multi-angle 3D images of urban buildings. Optical satellites urgently need to address the issue of ensuring mapping accuracy under large-angle side-swing conditions, enabling agile maneuvering mapping. During large-angle side-swing, atmospheric refraction, ground object obstruction, and resolution differences all affect the accuracy of mapping processing. However, limited satellite resources and a lack of on-orbit image samples from different angles mean that the relationship between side-swing angle and positioning accuracy is only qualitative, lacking quantitative analysis. Therefore, imaging simulation is needed for agile mapping and 3D modeling tasks to support the design and analysis of integrated satellite-ground performance indicators.

[0003] For satellite imaging simulation methods and software systems, published patents mainly focus on satellites used for traditional imaging, lacking simulations for 3D imaging of urban buildings. Urban buildings are characterized by rapid height changes and numerous vertical facades. Chinese patent application CN117593943A discloses a dynamic scene optical simulation system and imaging method for agile satellites, proposing a semi-physical satellite imaging mission simulation system that projects array images onto a dome screen and uses a camera simulator to image the simulated scene. However, this invention has limited ability to simulate scenes and satellite states, and cannot construct 3D realistic scenes with vertical facade features, such as urban buildings. Chinese patent CN105138756B, a satellite agile imaging simulation and positioning accuracy evaluation method, and Chinese patent application CN106126839A, a three-line array stereo mapping satellite imaging simulation method and system, use terrain products such as DOM data, DEM data, and laser point cloud data as scenes. These can only reflect horizontal or terrain undulation information, failing to reflect the vertical facade details of urban buildings and thus unable to support 3D realistic modeling and simulation.

[0004] Meanwhile, existing technologies often focus on simulating or evaluating single elements, neglecting the coupled effects of multiple factors. For example, Chinese patent CN105528500B, a decimeter-level satellite-borne TDICCD stereo mapping camera imaging simulation method and system, uses Worldview-3 stereo images as the data source, and Chinese patent application CN106126839A, a three-line array stereo mapping satellite imaging simulation method and system, is geared towards traditional steady-state pushbroom satellites. Neither of these technologies considers the imaging characteristics under large-angle side-swing conditions, cannot simulate the special effects of atmospheric refraction and ground object obstruction under large-angle conditions, and cannot provide support for the accuracy analysis of satellites with large-angle side-swing. Summary of the Invention

[0005] The purpose of this invention is to overcome the shortcomings of the existing technology and provide a three-dimensional imaging simulation method for optical mapping satellite constellations, which can generate three-dimensional images of urban buildings containing information on each facade, and meet the requirements of three-dimensional imaging simulation under different imaging angles and other conditions.

[0006] The objective of this invention can be achieved through the following technical solutions:

[0007] A method for 3D imaging simulation of an optical mapping satellite constellation, the method comprising:

[0008] Constructing a realistic 3D model based on aerial oblique photogrammetry and 3D lighting rendering;

[0009] Acquire preset satellite imaging data, construct image-space vector and object-space vector, transform the image-space vector and object-space vector to the body coordinate system, and construct an imaging geometric model that considers satellite platform error parameters, satellite camera error parameters and atmospheric error parameters;

[0010] In the real-world 3D model, all pixels of the imaging plane are traversed, and the image-side coordinates corresponding to the pixels are substituted into the imaging geometry model to generate a simulation image.

[0011] The simulated images are previewed in real time through the interface, and satellite auxiliary data is also output based on the imaging geometry model during the output process.

[0012] Furthermore, the three-dimensional lighting rendering method is as follows: when constructing the real-world three-dimensional model, lighting rendering is performed on the three-dimensional scene;

[0013] The process of rendering the lighting of the three-dimensional scene includes: calculating the altitude angle and azimuth angle of the sun based on the preset satellite imaging time and the imaging target position; calculating the spectral radiance received at various locations in the three-dimensional scene based on the altitude angle and azimuth angle of the sun; and simulating the lighting of the three-dimensional scene based on the spectral radiance.

[0014] Furthermore, the expression for calculating the spectral radiance is as follows:

[0015]

[0016] Among them, L λ,sensor F represents the spectral radiance. S This is the projection shadow coefficient. E is the solar zenith angle. λ,Direct V represents the direct solar spectral irradiance. S ρ is the sky observation factor. λ τ is the surface reflectance. λ,Atm For upward atmospheric transmittance, L λ,path Radiance;

[0017] In the calculation of the spectral radiance, the projection shadow coefficient is reduced to simulate shadow effects, taking into account the degree of occlusion of the face away from the sun and ground objects in the three-dimensional scene.

[0018] Furthermore, the process of generating the simulated image includes:

[0019] Based on the three-dimensional ground points of the object in the real scene 3D model, the elevation values ​​of each geographical location are obtained by bilinear interpolation, and then combined with the planar coordinates of each geographical location to obtain the three-dimensional coordinates of each geographical location.

[0020] Based on the three-dimensional coordinates of each geographical location, the pixel grayscale value corresponding to each geographical location is obtained by interpolation;

[0021] Define the imaging plane, traverse all pixels of the imaging plane in the real scene 3D model, and substitute the image-side coordinates of the pixels into the imaging geometry model. Establish a one-to-one correspondence between the image-side pixels and the object-side 3D ground points through the imaging geometry model.

[0022] The pixel grayscale values ​​of each geographical location are assigned to the corresponding pixels on the imaging plane to generate a simulated image.

[0023] Furthermore, the preset satellite imaging data includes preset attitude information and orbit information of the satellite during its on-orbit operation, as well as the imaging time of each row during satellite pushbroom imaging.

[0024] Furthermore, the process of constructing the image-space vector and the object-space vector includes:

[0025] The imaging time of each row during satellite pushbroom imaging is obtained, and interpolation is performed based on the row number to obtain the accurate imaging time of the image points in each row. Based on the accurate imaging time of each image point, the preset orbit information and attitude data are interpolated to obtain the accurate attitude and orbit information of the satellite at the imaging time. Based on the attitude and orbit information, image-space and object-space vectors are constructed.

[0026] Furthermore, the satellite auxiliary data output according to the imaging geometry model includes: attitude information including measurement errors, orbit information including measurement errors, and imaging time including measurement errors.

[0027] Furthermore, the satellite platform error parameters include the installation parameters of the platform sensors, the platform thermal deformation parameters, the satellite orbit parameters, the satellite attitude parameters, and the platform micro-vibration parameters;

[0028] The installation parameters of the platform sensors include the preset installation angle of the satellite camera, the installation angle of the star sensor, and the installation position of the GPS phase center;

[0029] The platform's thermal deformation parameters include the thermal deformation error of the star sensor, and its calculation expression is as follows:

[0030]

[0031] Among them, V LFE It is the thermal error of the star sensor, ω is the fundamental frequency of the Fourier series, and a j b represents the magnitude of the thermal error in the cosine component. j Let be the amplitude of the thermal error in the sinusoidal component, and k be the independent variable;

[0032] The satellite orbital parameters are the preset number of satellite orbital elements;

[0033] The satellite attitude parameters include the number of imaging angles used by multiple satellites in a single mission (preset);

[0034] The micro-vibration parameters of the platform are obtained based on the typical vibration frequency bands of the satellite's moving parts, and their calculation expression is as follows:

[0035]

[0036] Among them, J t Let A be the satellite flutter at time i. i Let ω be the amplitude of satellite flutter at time i. i Let i be the frequency of satellite flutter at time i. Let be the phase angle of satellite flutter at time i, and t be the moment of satellite flutter.

[0037] Furthermore, the satellite camera error parameters include camera principal distance distortion and camera lens geometric distortion;

[0038] The calculation expression for the image point change caused by the camera principal distance distortion in the satellite image is as follows:

[0039]

[0040] Δf=f′-f

[0041] Where f is the laboratory calibration value, f′ is the camera principal distance after focusing, and Δx f Let Δy be the x-axis translation of the image point. f Let x be the y-axis translation of the image point, Δf be the variation parameter, x be the ideal x-coordinate of the image point, and y be the ideal y-coordinate of the image point.

[0042] The geometric distortion of the camera lens includes radial distortion and tangential distortion.

[0043] The calculation expression for the image point change caused by the radial distortion in the satellite image is as follows:

[0044]

[0045] Where, Δx r Let Δy be the radial distortion along the x-axis of the image point. r Let x be the radial distortion along the y-axis of the image point. fp Let x be the x-axis coordinate of the image point, y be the y-axis coordinate of the image point fp Let k1, k2, and k3 be the y-axis coordinates of the image point, and k1, k2, and k3 be the radial distortion coefficients of the camera lens.

[0046] The calculation expression for the image point change caused by the tangential distortion in the satellite image is as follows:

[0047]

[0048] Where p1 and p2 are tangential distortion coefficients.

[0049] Furthermore, the atmospheric error parameters include image point offsets in satellite images caused by atmospheric refraction;

[0050] The expression for calculating the image point offset is:

[0051] Δx = (n-1)·tanθ·(x-x0)

[0052] Δy=(n-1)·tanθ·(y-y0)

[0053] Where Δx and Δy represent the changes in the x and y coordinates of the image point, respectively, n is the refractive index of atmospheric refraction, θ is the incident angle of light entering the atmosphere, x represents the original x coordinate of the image point, y represents the original y coordinate of the image point, and x0 and y0 represent the x and y coordinates of the atmospheric refraction reference point, respectively.

[0054] Compared with the prior art, the beneficial effects of the present invention include:

[0055] 1. This invention targets 3D imaging satellite constellations, possessing simulation capabilities for large-angle, multi-view imaging conditions. It establishes a rigorous geometric imaging model considering all factors across the entire process, using realistic high-resolution, high-precision 3D scenes as a base. This allows for flexible simulation of complete 3D images of urban buildings, including information on each facade, accurately reflecting the parallax of ground feature heights on the image plane. Unrestricted by imaging conditions, it can obtain simulated images under any geometric state, providing ample samples for evaluating the 3D modeling capabilities of ultra-high-resolution optical mapping satellite constellations. Furthermore, this invention constructs a realistic 3D model, in which 3D lighting rendering is performed, enabling the acquisition of images under any geometric state, regardless of the actual imaging conditions of the satellite in orbit. The simulated images provided can be directly used to support the platform performance demonstration and application capability evaluation of large-angle agile mapping satellites, such as the quantitative analysis of the side-sway angle and positioning accuracy. This solves the pain point of existing technologies that lack on-orbit image samples from different angles and cannot quantitatively analyze the impact of side-sway on accuracy. Based on real high-resolution, high-precision 3D scenes, this invention takes into account satellite platform errors, camera payload errors, and atmospheric transmission errors to establish an imaging geometric model for simulation imaging. The simulation results can truly reflect the geometric information of the 3D imaging of the ultra-high resolution optical mapping satellite constellation, especially the parallax of ground object height on the image plane, providing high-precision simulation samples from multiple angles for subsequent geometric positioning processing and 3D modeling capability evaluation.

[0056] 2. This invention uses aerial oblique photogrammetry to construct a real-world 3D model, extending the simulation scene to a full 3D real-world scene, and realizing a complete simulation of the three-dimensional features of ground objects.

[0057] 3. This invention considers the coupled effects of multiple factors such as atmospheric refraction, platform thermal deformation, and camera distortion, and incorporates the quantification of errors in each stage into the simulation, making the results closer to the actual imaging process of the satellite in orbit, and providing more comprehensive support for the design and analysis of integrated space-ground performance indicators.

[0058] 4. This invention generates a high-precision simulation image and a complete dataset containing error auxiliary data. The output results can be directly used for quantitative analysis of lateral tilt angle and positioning accuracy, and evaluation of 3D modeling accuracy. Attached Figure Description

[0059] Figure 1 This is a flowchart of the method of the present invention;

[0060] Figure 2 This is a flowchart of the process of generating simulated images in this invention. Detailed Implementation

[0061] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort should fall within the scope of protection of the present invention.

[0062] Example 1

[0063] This embodiment discloses a three-dimensional imaging simulation method for an optical mapping satellite constellation, the specific method being as follows: Figure 1 As shown, the method includes steps S1-S5. A constellation specifically refers to a group of satellites consisting of multiple ultra-high resolution optical mapping satellites that work collaboratively in a specific orbital configuration. The specific steps of this method are as follows:

[0064] Step S1: Construct a real-world 3D model based on aerial oblique photogrammetry and 3D lighting rendering.

[0065] The satellite constellation consists of multiple ultra-high resolution optical mapping satellites. Traditional two-dimensional terrain scenes only require one satellite to support single-view imaging. However, when the satellite constellation collects images of cities, due to the many details of the three-dimensional vertical facades of urban images (the height differences of urban buildings), multiple satellites need to be deployed on different orbital planes to form a side-by-side configuration. Multiple satellites synchronously cover the same imaging area at similar times and collect images of the same imaging area. Multiple satellites combine to collect images from different angles such as frontal, side, and oblique views to ensure that the simulated images can capture different facade features of the ground objects (such as the front, side, and top of buildings). However, during simulation, imaging simulation cannot be carried out simply from the perspective of a single satellite. Therefore, in this embodiment, a real-scene three-dimensional model is established to simulate the three-dimensional scene of the ground objects.

[0066] The purpose of constructing a 3D scene of ground features is to provide the most realistic signal source possible for imaging simulation. The geometric and texture information of the 3D scene must be as consistent as possible with the real scene. The resolution of the 3D scene should be more than twice the minimum resolution of the imaging simulation task. This simulation method is aimed at urban 3D modeling satellite constellations. The resolution of satellite imaging reaches below 0.3m. Therefore, the resolution of the 3D scene should be less than 0.15m.

[0067] The 3D scene is a high-resolution, realistic 3D model built using aerial oblique photogrammetry. This model boasts extremely high resolution and geometric accuracy. While creating the model using aerial oblique photogrammetry, open-source realistic 3D models from the internet were also downloaded as a supplement to the simulation system. These open-source realistic 3D models are generally also based on aerial oblique photogrammetry, featuring wide coverage and meeting the required resolution and geometric accuracy.

[0068] After establishing the real-scene 3D model, it is also necessary to manage the real-scene 3D model and generate corresponding Shapefile files for the acquired 3D scene. In the simulation system, the distribution and coverage of the existing real-scene 3D model from a global perspective can be viewed intuitively, which facilitates simulation applications.

[0069] In step S1, 3D lighting rendering specifically means rendering the 3D scene with lighting according to simulation requirements when constructing a real-world 3D model.

[0070] The process of rendering lighting in a 3D scene includes:

[0071] Calculate the solar elevation angle and azimuth angle based on the preset satellite imaging time and the imaging target location;

[0072] Calculate the spectral radiance received at various locations in the three-dimensional scene based on the sun's elevation angle and azimuth angle;

[0073] Simulated lighting for 3D scenes based on spectral radiance.

[0074] Assuming the Earth's surface is equivalent to a Lambertian reflector and the surface reflectivity is ρ λ The expression for calculating spectral radiance is:

[0075]

[0076] Among them, L λ,sensor F represents the spectral radiance. S This is the projection shadow coefficient. E is the solar zenith angle. λ,Direct V represents the direct solar spectral irradiance. S ρ is the sky observation factor. λ τ is the surface reflectance. λ,Atm For upward atmospheric transmittance, L λ,path Radiance;

[0077] In the calculation of spectral radiance, the projection shadow coefficient is reduced to simulate shadow effects, taking into account the degree of occlusion of the face away from the sun and ground objects in the 3D scene.

[0078] In this embodiment, if the scene is a shadow area, then F S =0, otherwise F S =1.

[0079] The sky observation factor is defined as the ratio of the diffuse reflection of the sky received at a point to the diffuse reflection received by an unobstructed horizontal plane, and it is between 0 and 1.

[0080] Solar direct spectral irradiance E λ,Direct Diffuse spectral irradiance E λ,DiffuseAscending atmospheric transmittance τ λ,Atm and process radiance L λ,path The four atmospheric physical quantities can be calculated directly using existing software based on the solar zenith angle, time, and geographical location.

[0081] Step S2: Obtain preset satellite imaging data, construct image-space vector and object-space vector, transform the image-space vector and object-space vector to the body coordinate system, and construct an imaging geometric model that considers satellite platform error parameters, satellite camera error parameters and atmospheric error parameters.

[0082] The preset satellite imaging data includes the preset attitude information and orbit information of the satellite during its on-orbit operation, as well as the imaging time of each row during satellite pushbroom imaging.

[0083] Before constructing the simulation imaging model, it is necessary to unify the time system of each preset satellite imaging data, that is, to construct the image-space vector and the object-space vector.

[0084] The process of constructing the image-space vector and the object-space vector includes:

[0085] The imaging time of each row during satellite pushbroom imaging is obtained, and interpolation is performed based on the row number to obtain the accurate imaging time of the image points in each row.

[0086] Based on the accurate imaging time of each image point, the satellite's accurate attitude and orbit information at the imaging time is obtained by interpolation of the preset orbit information and attitude data. Based on the attitude and orbit information, image-space and object-space vectors are constructed.

[0087] Since linear array satellites perform pushbroom imaging along the rows, the accurate imaging time of the image points can be obtained by interpolating the row numbers. This step usually uses the Lagrange interpolation method or the polynomial interpolation method.

[0088] The image-space vector can be converted to the image space coordinate system based on the row and column numbers in the image coordinate system. Therefore, the image-space vector can be represented as:

[0089] The object vector in the WGS84 coordinate system can be represented as: Transform the object vector to the body coordinate system:

[0090]

[0091] This allows for the construction of an ideal imaging geometry model:

[0092]

[0093] This ideal imaging geometry model is an imaging geometry model that does not contain errors. Using this imaging geometry model for simulation can obtain simulated images under ideal conditions.

[0094] Since there are various errors in the actual imaging process, we will introduce various errors into the imaging geometry model, including satellite platform error parameters, satellite camera error parameters and atmospheric error parameters.

[0095] The satellite platform error parameters include the installation parameters of the platform sensors, the platform thermal deformation parameters, the satellite orbit parameters, the satellite attitude parameters, and the platform micro-vibration parameters. These parameters are given in the form of an installation matrix.

[0096] The platform sensor installation parameters include the preset installation angle of the satellite camera, the installation angle of the star sensor, and the installation position of the GPS phase center.

[0097] The thermal deformation parameters of the platform include the thermal error of the star sensor. The main cause of platform thermal deformation is the drastic temperature fluctuations caused by the satellite entering and exiting the sun's shadow. The onboard thermal control system can only control the temperature within a certain range, not maintain a constant temperature, thus causing thermal deformation of the star sensor mounting structure and resulting in thermal error. Since the angle of solar illumination varies with the orbital period, the thermal error can be considered a periodic signal, with its period approximately equal to the satellite's orbital period. Traditional EKF (Extended Kalman Filter) methods can only filter out high-frequency noise components and cannot eliminate low-frequency components. Therefore, new methods must be used to reduce the impact of thermal error on satellite attitude determination.

[0098] Since the frequency of the thermal error is unknown, we assume its period is the orbital period and half of the orbital period, i.e., ω and 2ω. Let the thermal error of the star sensor be V. LFE Because of V LFE It changes periodically, therefore it can be written in the form of a Fourier series, that is:

[0099]

[0100] Where ω is the fundamental frequency of the Fourier series, and V LFE ,a j ,b j Both are 3×1 vectors, corresponding to the three directions of the star sensor coordinate system, a j b represents the magnitude of the thermal error in the cosine component. j Let ω be the amplitude of the thermal error in the sinusoidal component, and k be the independent variable. If ω and a can be identified... j and b jThis allows for compensation of thermal error. In practical applications, when using the EKF method to determine the attitude of a gyroscope / star sensor system, the thermal error of the star sensor will be reflected in the estimation result of the constant drift. That is, the estimation of the constant drift will also show the characteristic of periodic change. Therefore, it is possible to identify the parameters of thermal error from the estimation of the constant drift and compensate for them.

[0101] The satellite orbital parameters are preset orbital elements. During the design phase, the satellite's orbital elements are designed according to its type and mission objectives. Multiple satellites in the constellation operate on different orbital planes, forming a side-by-side configuration, and simultaneously pass over the imaging area, acquiring images from different angles.

[0102] The satellite attitude parameters represent the number of imaging angles used by multiple satellites in a single mission. During a single imaging session, the attitude angles of each satellite are fixed. After imaging at the current angle is completed, the satellite quickly maneuvers to the next fixed attitude angle to complete the pushbroom imaging of the next viewpoint. The imaging area is the same for each imaging attitude, and the attitude angles are given in the form of three-axis Euler angles or quaternions. During simulation, these are converted into attitude matrices and substituted into the simulation imaging model.

[0103] Satellite attitude parameters also include attitude stability, which affects the modulation transfer function (MTF) index of imaging. A model is established to illustrate the impact of attitude stability on the MTF. For a TDI (Time-Delay Integration) CCD optical system, the MTF error caused by image shift matching error occurs at the Nyquist frequency f. N The MTF at this location is calculated as follows:

[0104] When the attitude angle occurs at a rate of When the time changes linearly, the image displacement on the image plane caused by the attitude angular rate during the integration time is:

[0105]

[0106] The impact on MTF is calculated as follows:

[0107]

[0108] Where, N TDI For the TDICCD integral series, f N For the Nyquist frequency, d is the cell size, t e is the integration time, and f is the camera focal length.

[0109] The platform's micro-vibration parameters are obtained based on the typical vibration frequency bands of the satellite's moving parts. According to the Fourier analysis principle, the platform flutter waveform is decomposed into a superposition of multiple simple harmonics, and its calculation expression is as follows:

[0110]

[0111] Among them, J t Let A be the satellite flutter at time i. i Let ω be the amplitude of satellite flutter at time i. i Let i be the frequency of satellite flutter at time i. Let be the phase angle of satellite flutter at time i, and t be the moment of satellite flutter.

[0112] Satellite camera error parameters include camera principal distance distortion and camera lens geometric distortion;

[0113] Camera principal distance distortion can cause image point translation. Let the principal distance of the focused camera be f', the laboratory calibration value be f, and the focal length error be:

[0114] Δf=f′-f

[0115] For a point p on the camera's focal plane, its coordinates on the focal plane are (x, y). The change in image point p caused by the camera's principal distance distortion (Δx) is... f ,Δy f )for:

[0116]

[0117] Where f is the laboratory calibration value, f′ is the camera principal distance after focusing, Δf is the focal length error, x is the ideal x-coordinate of the image point, and y is the ideal y-coordinate of the image point.

[0118] Camera lens geometric distortion is a non-linear error that includes radial and tangential distortion of the camera lens.

[0119] The formula for calculating the image point change caused by radial distortion in satellite images is as follows:

[0120]

[0121] Where, Δx r Let Δy be the radial distortion along the x-axis of the image point. r Let x be the radial distortion along the y-axis of the image point. fp Let x be the x-axis coordinate of the image point, y be the y-axis coordinate of the image point fp Let k1, k2, and k3 be the y-axis coordinates of the image point, and k1, k2, and k3 be the radial distortion coefficients of the camera lens.

[0122] r is the distance from the image point to the principal point, i.e.

[0123] The expression for calculating the image point change caused by tangential distortion in satellite images is as follows:

[0124]

[0125] Where p1 and p2 are tangential distortion coefficients.

[0126] Atmospheric error parameters include image point shifts in satellite images caused by atmospheric refraction;

[0127] The expression for calculating the image point offset is:

[0128] Δx = (n-1)·tanθ·(x-x0)

[0129] Δy=(n-1)·tanθ·(y-y0)

[0130] Where Δx and Δy represent the changes in the x and y coordinates of the image point, respectively, n is the refractive index of atmospheric refraction, θ is the incident angle of light entering the atmosphere, x represents the original x coordinate of the image point, y represents the original y coordinate of the image point, and x0 and y0 represent the x and y coordinates of the atmospheric refraction reference point, respectively.

[0131] In this embodiment, the refractive index of atmospheric refractive difference is 1.0001787 for the troposphere and 1.0000167 for the stratosphere.

[0132] Step S3: Traverse all pixels of the imaging plane in the real-world 3D model, substitute the image coordinates corresponding to the pixels into the imaging geometry model, and generate a simulation image.

[0133] The specific process of generating the simulated image includes:

[0134] Step S301, calculate the elevation values ​​of the three-dimensional ground points on the object:

[0135] The three-dimensional ground points in the object area refer to the ground features within the simulation area, such as building vertices and ground feature points. The basic elevation framework of each geographical location in the scene is stored in the real-scene three-dimensional model. For any ground point in the real-scene three-dimensional model that does not directly provide an elevation value, the elevation is calculated using bilinear interpolation. That is, by using the coordinates and elevation data of the four known elevation points around the point, the accurate elevation value of the point is fitted to ensure that all ground points in the object area have complete three-dimensional coordinates. Through this step, the three-dimensional coordinates of each geographical location (i.e., the three-dimensional coordinates of the three-dimensional ground points in the object area) can be obtained.

[0136] Step S302, interpolate to obtain pixel grayscale values ​​in the 3D model:

[0137] In step S1, the 3D scene lighting rendering has given the 3D model realistic texture and radiation characteristics, such as building wall texture, ground material reflectivity and radiation intensity attenuation in shadow areas. Based on the 3D coordinates of each geographical location, coordinate matching and interpolation are performed in the texture data of the real-world 3D model. If the geographical location point is exactly located at the center of the texture pixel of the 3D model, the gray value of the corresponding pixel is directly obtained. If the point is located in the pixel gap, the gray value of the point is generated by interpolating the gray values ​​of the surrounding pixels (such as bilinear interpolation) to ensure the continuity and realism of radiation information. The pixel gray value corresponding to each geographical location can be obtained through interpolation.

[0138] Step S303: Traverse the pixels of the imaging plane and substitute them into the imaging geometry model to generate a simulation image:

[0139] Define the imaging plane and specify parameters such as the number of rows and columns of pixels and the pixel size. The imaging plane is the focal plane of the camera on the satellite, corresponding to the pixel array of the simulated image. Traverse all pixels of the imaging plane in the real-world 3D model and substitute the image-side coordinates of the pixels into the imaging geometry model. Establish a one-to-one correspondence between image-side pixels and object-side 3D ground points through the imaging geometry model. The difference between image-side coordinates and commonly used 3D coordinates is that the z-coordinate is replaced by the camera focal length f. Assign the pixel grayscale values ​​of each geographical location to the corresponding pixels on the imaging plane to generate the simulated image under ideal conditions. Then, introduce the error term in the imaging geometry model to obtain an accurate simulated image.

[0140] In the imaging simulation of step S3, since the satellite constellation contains multiple satellites, steps S301-S303 will be repeated repeatedly during the imaging process, requiring the input of orbital elements, satellite attitude angles, and camera parameters for different satellites. For example, for the same 3D scene, satellite 1 (orbital plane 1, attitude angle α), satellite 2 (orbital plane 2, attitude angle β), and satellite 3 (orbital plane 3, attitude angle γ) are used as imaging simulation objects. Their respective attitude and orbit information and attitude matrices are used for imaging simulation, generating multiple sets of simulated images from different perspectives in parallel. The simulated image of satellite 2 is a frontal view showing the top of ground features, while the simulated images of satellites 1 and 3 are side views showing the sides of buildings. These multi-view images collectively contain complete 3D information of the ground features, realistically simulating the process and results of multi-satellite image acquisition.

[0141] Step S4: Simulate satellite auxiliary data based on the imaging geometry model.

[0142] The satellite auxiliary data output by the imaging geometry model includes: attitude information with measurement errors, orbit information with measurement errors, and imaging time with measurement errors.

[0143] The attitude measurement data is based on the accurate attitude and orbit information in step S2. The theoretical attitude data is output at an update frequency of 10Hz. The measurement error of the attitude measurement sensor is then added to this theoretical attitude to obtain attitude information including the measurement error, simulating the measurements obtained by the satellite in orbit through the attitude measurement sensor. The measurement error is simulated according to three different variation periods: low frequency, mid frequency, and high frequency. The low-frequency error term represents the thermal low-frequency term of the star sensor, consistent with the orbital period, and is modeled using a sine function. The mid-frequency error term represents the spatial low-frequency error term of the star sensor, consistent with the time it takes for a star to pass through the star sensor's field of view, and is modeled using a sine function. The high-frequency error term is the measurement noise of the star sensor, exhibiting random characteristics, and is simulated using random numbers generated according to a Gaussian distribution.

[0144] The simulation of imaging time including measurement errors, i.e., time synchronization errors, involves navigation receiver measurement data, attitude measurement data, and image data. All three types of measurement data have time stamp bits, with a time resolution of 1 microsecond. Using the theoretical time corresponding to the orbital position as a reference, the orbital position timescale is obtained by adding a normally distributed navigation receiver clock error; the attitude measurement system time synchronization error, which conforms to a normally distributed timescale, is used as the attitude measurement parameter timescale; and the camera time synchronization error, which conforms to a normally distributed timescale, is used as the image timescale.

[0145] The orbital information, including measurement errors, is calculated using the number of orbits provided in step S2 to obtain the satellite's theoretical coordinates in an inertial coordinate system or a fixed Earth coordinate system, and output at a frequency of 1Hz. Random noise conforming to a normal distribution is added to the satellite's theoretical coordinates to simulate the measurement data from the navigation receiver.

[0146] Step S5: Preview and output simulation images and satellite-aided data in real time through the interface.

[0147] The output simulated imagery consists of 3D images of urban buildings containing information on each facade. By overlaying simulated multi-view images from multiple satellites, the parallax of ground feature heights across different image planes can be clearly presented. This result provides a multi-angle sample library for subsequent 3D modeling capability assessment. For example, by analyzing the parallax of simulated images from different satellites at different side-swing angles, the 3D image acquisition accuracy of the satellite constellation can be quantitatively evaluated.

[0148] The real-time preview interface is based on a preset parameter setting interface. This interface is designed to provide an intuitive, easy-to-use, and fully functional environment for adjusting parameters, meeting the needs of users performing satellite imaging simulations. Simulated images are previewed in real-time through the interface. Simultaneously, auxiliary data and image data are framed and stored on disk according to the agreed-upon image data protocol.

[0149] This simulation method is based on a real-world 3D model built using pre-set satellite data and aerial oblique photogrammetry technology. It does not require actual satellite data acquisition for simulation, and the simulation results are highly comprehensive, allowing for the evaluation of satellite acquisition accuracy both before and after launch. During the satellite design phase, the simulation system can be used to repeatedly adjust input parameters, quantifying the impact of different platform errors on the final image quality. This provides precise data support for the design of satellite platform and payload specifications, avoiding over- or under-design. After launch, before executing important mapping tasks, mission simulations can be conducted in the system to simulate satellite transit, set different imaging angles, and predict the image coverage, overlap, and stereoscopic imaging effects under different planning schemes. This allows for the selection of the optimal imaging scheme, improving satellite efficiency and data acquisition success rate. Furthermore, simulation can accurately predict the final image quality, providing a reliable guarantee for market promotion and application implementation.

[0150] Example 2

[0151] Based on Embodiment 1, this embodiment provides an electronic device, including: one or more processors and a memory, wherein the memory stores one or more programs, the one or more programs including instructions for executing the aforementioned optical mapping satellite constellation three-dimensional imaging simulation method.

[0152] At the hardware level, the electronic device includes a processor, internal bus, network interface, memory, and non-volatile memory, and may also include other hardware required for business operations. The processor reads the corresponding computer program from the non-volatile memory into memory and then runs it to implement the aforementioned three-dimensional imaging simulation method for optical mapping satellite constellations. Of course, in addition to software implementation, this invention does not exclude other implementation methods, such as logic devices or a combination of hardware and software, etc. That is to say, the execution subject of the following processing flow is not limited to individual logic units, but can also be hardware or logic devices.

[0153] Memory may include non-persistent storage in computer-readable media, such as random access memory (RAM) and / or non-volatile memory, such as read-only memory (ROM) or flash RAM. Memory is an example of computer-readable media.

[0154] Computer-readable media includes both permanent and non-permanent, removable and non-removable media that can store information using any method or technology. Information can be computer-readable instructions, data structures, modules of programs, or other data. Examples of computer storage media include, but are not limited to, phase-change memory (PRAM), static random access memory (SRAM), dynamic random access memory (DRAM), other types of random access memory (RAM), read-only memory (ROM), electrically erasable programmable read-only memory (EEPROM), flash memory or other memory technologies, CD-ROM, digital versatile optical disc (DVD) or other optical storage, magnetic tape, disk storage or other magnetic storage devices, or any other non-transferable medium that can be used to store information accessible by a computing device. As defined herein, computer-readable media does not include transient computer-readable media, such as modulated data signals and carrier waves.

[0155] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any person skilled in the art can easily conceive of various equivalent modifications or substitutions within the technical scope disclosed in the present invention, and these modifications or substitutions should all be covered within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of the claims.

Claims

1. A three-dimensional imaging simulation method for an optical mapping satellite constellation, characterized in that, The method includes: Constructing a realistic 3D model based on aerial oblique photogrammetry and 3D lighting rendering; Acquire preset satellite imaging data, construct image-space vector and object-space vector, transform the image-space vector and object-space vector to the body coordinate system, and construct an imaging geometric model that considers satellite platform error parameters, satellite camera error parameters and atmospheric error parameters; In the real-world 3D model, all pixels of the imaging plane are traversed, and the image-side coordinates corresponding to the pixels are substituted into the imaging geometry model to generate a simulation image. The simulated images are previewed in real time through the interface, and satellite auxiliary data is also output based on the imaging geometry model during the output process.

2. The method for 3D imaging simulation of an optical mapping satellite constellation according to claim 1, characterized in that, The three-dimensional lighting rendering method is as follows: when constructing the real-world three-dimensional model, the three-dimensional scene is rendered with lighting. The process of rendering the lighting of the three-dimensional scene includes: calculating the altitude angle and azimuth angle of the sun based on the preset satellite imaging time and the imaging target position; calculating the spectral radiance received at various locations in the three-dimensional scene based on the altitude angle and azimuth angle of the sun; and simulating the lighting of the three-dimensional scene based on the spectral radiance.

3. The method for 3D imaging simulation of an optical mapping satellite constellation according to claim 2, characterized in that, The formula for calculating the spectral radiance is as follows: Among them, L λ,sensor F represents the spectral radiance. S θ is the projection shadow coefficient. solar E is the solar zenith angle. λ,Direct V represents the direct solar spectral irradiance. S ρ is the sky observation factor. λ τ is the surface reflectance. λ,Atm For upward atmospheric transmittance, L λ,path Radiance; In the calculation of the spectral radiance, the projection shadow coefficient is reduced to simulate shadow effects, taking into account the degree of occlusion of the face away from the sun and ground objects in the three-dimensional scene.

4. The method for 3D imaging simulation of an optical mapping satellite constellation according to claim 1, characterized in that, The process of generating the simulated image includes: Based on the three-dimensional ground points of the object in the real scene 3D model, the elevation values ​​of each geographical location are obtained by bilinear interpolation, and then combined with the planar coordinates of each geographical location to obtain the three-dimensional coordinates of each geographical location. Based on the three-dimensional coordinates of each geographical location, the pixel grayscale value corresponding to each geographical location is obtained by interpolation; Define the imaging plane, traverse all pixels of the imaging plane in the real scene 3D model, and substitute the image-side coordinates of the pixels into the imaging geometry model. Establish a one-to-one correspondence between the image-side pixels and the object-side 3D ground points through the imaging geometry model. The pixel grayscale values ​​of each geographical location are assigned to the corresponding pixels on the imaging plane to generate a simulated image.

5. The method for 3D imaging simulation of an optical mapping satellite constellation according to claim 1, characterized in that, The preset satellite imaging data includes preset attitude information and orbit information of the satellite during its on-orbit operation, as well as the imaging time of each row during satellite pushbroom imaging.

6. The three-dimensional imaging simulation method for an optical mapping satellite constellation according to claim 5, characterized in that, The process of constructing the image-space vector and the object-space vector includes: The imaging time of each row during satellite pushbroom imaging is obtained, and interpolation is performed based on the row number to obtain the accurate imaging time of the image points in each row. Based on the accurate imaging time of each image point, the preset orbit information and attitude data are interpolated to obtain the accurate attitude and orbit information of the satellite at the imaging time. Based on the attitude and orbit information, image-space and object-space vectors are constructed.

7. The three-dimensional imaging simulation method for an optical mapping satellite constellation according to claim 5, characterized in that, The satellite auxiliary data output from the imaging geometry model includes: attitude information with measurement errors, orbit information with measurement errors, and imaging time with measurement errors.

8. The method for 3D imaging simulation of an optical mapping satellite constellation according to claim 1, characterized in that, The satellite platform error parameters include the installation parameters of the platform sensors, the platform thermal deformation parameters, the satellite orbit parameters, the satellite attitude parameters, and the platform micro-vibration parameters. The installation parameters of the platform sensors include the preset installation angle of the satellite camera, the installation angle of the star sensor, and the installation position of the GPS phase center; The platform's thermal deformation parameters include the thermal deformation error of the star sensor, and its calculation expression is as follows: Among them, V LFE It is the thermal error of the star sensor, ω is the fundamental frequency of the Fourier series, and a j b represents the magnitude of the thermal error in the cosine component. j Let be the amplitude of the thermal error in the sinusoidal component, and k be the independent variable; The satellite orbital parameters are the preset number of satellite orbital elements; The satellite attitude parameters include the number of imaging angles used by multiple satellites in a single mission (preset); The micro-vibration parameters of the platform are obtained based on the typical vibration frequency bands of the satellite's moving parts, and their calculation expression is as follows: Among them, J t Let A be the satellite flutter at time i. i Let ω be the amplitude of satellite flutter at time i. i Let i be the frequency of satellite flutter at time i. Let be the phase angle of satellite flutter at time i, and t be the moment of satellite flutter.

9. The three-dimensional imaging simulation method for an optical mapping satellite constellation according to claim 1, characterized in that, The satellite camera error parameters include camera principal distance distortion and camera lens geometric distortion; The calculation expression for the image point change caused by the camera principal distance distortion in the satellite image is as follows: Δf=f ′ -f Where f is the laboratory calibration value, f ′ Δx represents the camera's principal distance after focusing. f Let Δy be the x-axis translation of the image point. f Let x be the y-axis translation of the image point, Δf be the variation parameter, x be the ideal x-coordinate of the image point, and y be the ideal y-coordinate of the image point. The geometric distortion of the camera lens includes radial distortion and tangential distortion. The calculation expression for the image point change caused by the radial distortion in the satellite image is as follows: Where, Δx r Let Δy be the radial distortion along the x-axis of the image point. r Let x be the radial distortion along the y-axis of the image point. fp Let x be the x-axis coordinate of the image point, y be the y-axis coordinate of the image point fp Let k1, k2, and k3 be the y-axis coordinates of the image point, and k1, k2, and k3 be the radial distortion coefficients of the camera lens. The calculation expression for the image point change caused by the tangential distortion in the satellite image is as follows: Where p1 and p2 are tangential distortion coefficients.

10. The method for 3D imaging simulation of an optical mapping satellite constellation according to claim 1, characterized in that, The atmospheric error parameters include the image point offset in the satellite image caused by atmospheric refraction; The expression for calculating the image point offset is: Δx = (n-1)·tanθ·(x-x0) Δy=(n-1)·tanθ·(y-y0) Where Δx and Δy represent the changes in the x and y coordinates of the image point, respectively, n is the refractive index of the atmospheric refraction difference, θ is the incident angle of the light entering the atmosphere, x represents the original x coordinate of the image point, y represents the original y coordinate of the image point, and x0 and y0 represent the x and y coordinates of the atmospheric refraction difference reference point, respectively.

Images