High-orbit SAR permanent shadow region interpretation model based on orbit information and DEM product

Abstract

The invention relates to the technical field of radar remote sensing, and discloses a high-orbit SAR permanent shadow region interpretation model based on orbit information and a DEM product, which comprises an orbit position calculation module, an observation geometry construction module, an instantaneous shadow interpretation module and a permanent shadow determination module. According to the model, orbit description information of different types is uniformly converted into satellite space position coordinates under an earth-centered earth-fixed coordinate system; secondly, calculating a three-dimensional coordinate and a local normal vector of an earth surface pixel by using a digital elevation model, constructing a radar sight line vector and calculating an incident angle; then, judging a self-shadow state according to the change of an incident angle within synthetic aperture time, and detecting terrain shielding on a path by utilizing a line-of-sight propagation model so as to determine an instantaneous shadow; and finally, carrying out time sequence statistics on the instantaneous shadow state covering the complete orbital period, and when a threshold condition is met, determining that the shadow area is a permanent shadow area. According to the method, accurate calculation of the permanent shadow region of the high-orbit SAR under multiple orbit types is realized.

Description

Technical Field

[0001] This invention relates to the field of radar remote sensing technology, specifically to a high-orbit SAR permanent shadow area interpretation model based on orbit information and DEM products. Background Technology

[0002] Synthetic Aperture Radar (SAR) systems, as an all-weather, all-day active microwave remote sensing method, play a crucial role in Earth observation. Traditional SAR systems mostly operate in low Earth orbit, possessing high imaging resolution and stable global coverage. However, limited by their low orbital altitude, their temporal coverage exhibits periodic discontinuities, making it difficult to meet the needs of continuous, high-time-efficiency observation of specific areas. With the development of remote sensing technology, high-orbit platforms such as inclined geosynchronous orbits, quasi-geosynchronous orbits, Earth-Moon Lagrange points, and lunar bases have gradually become important carriers for wide-area monitoring. High-orbit SAR, with its orbital altitude advantage, can form stable observation geometry and remain near target areas for extended periods, showing application potential in wide-area disaster monitoring and strategic reconnaissance. However, the imaging geometry of high-orbit SAR differs from that of traditional low-orbit SAR. High-orbit observation involves signal propagation over extremely long distances, and the changes in observation angle and incident angle are more complex, exhibiting non-stationary, non-moving, and asymmetrical geometric characteristics. Existing shadow area interpretation methods and coverage analysis tools are mostly based on low-orbit imaging geometry, typically assuming that the distance between the observation platform and the surface target changes little and the incident angle range is limited. When these methods are applied to high-orbit scenarios, they are difficult to accurately handle the geometric distortions caused by large oblique angles and long-distance ground occlusion chains, resulting in the inability to accurately define the boundary between the effective imaging area and the terrain occlusion area.

[0003] Current coverage analysis tools primarily focus on the geometric line-of-sight coverage of the satellite over the ground, with less emphasis on incorporating terrain data for effectiveness assessment. In practical mission planning, calculating only the geometric coverage area while ignoring distortion effects such as shading, perspective shortening, and overlay caused by terrain can lead to assessment results that deviate from actual imaging capabilities. This is especially true under high-orbit observation conditions, where terrain occlusion caused by surface undulations is more pronounced. Furthermore, because high-orbit satellites have a slower angular velocity relative to the ground, some areas may remain in shadow for extended periods. Existing technologies lack the ability to perform temporal comprehensive analysis of multi-moment observation conditions, making it difficult to distinguish between occasional transient shadows and permanent shadows with engineering constraints, thus failing to provide accurate data for constellation design and energy strategies.

[0004] Furthermore, high-orbit SAR platforms employ diverse orbit description formats, including orbital elements and precise ephemeris data. Existing technologies are often designed for single data formats, lacking a universal model capable of accommodating multiple orbital inputs and uniformly performing spatial calculations. Simultaneously, when dealing with high-precision digital elevation models on a global scale, existing algorithms suffer from high computational complexity, hindering efficient and stable shadow region interpretation in large-scale high-orbit observation scenarios, thus limiting their engineering application in system validation and mission planning. Summary of the Invention

[0005] To address the shortcomings of existing technologies, this invention provides a permanent shadow area interpretation model for high-orbit SAR based on orbit information and DEM products. This model solves the problem that existing shadow interpretation methods are difficult to adapt to the complex imaging geometry and multi-source orbit data of high-orbit SAR, and cannot accurately define the permanent shadow area of ​​terrain under long-period observation.

[0006] To achieve the above objectives, the present invention is implemented through the following technical solution: a high-orbit SAR permanent shadow zone interpretation model based on orbit information and DEM products.

[0007] The model uses modules for orbital position calculation, observation geometry construction, instantaneous shadow interpretation, and permanent shadow determination to extract the permanent shadow region of high-orbit synthetic aperture radar under specific terrain conditions.

[0008] For orbital description information from different sources, this invention employs a unified coordinate calculation method to convert the position of a high-orbit SAR platform into spatial coordinates in a geocentric-geofixed coordinate system. For orbital element data, the geographic latitude tangent is obtained through corrections including Earth's oblateness, and longitude changes are analyzed using the westward retreat rate of the ascending node and the angular velocity of motion, thus completing the conversion from spherical coordinates to rectangular coordinates. For ephemeris data, a conversion between the time system and coordinate system is performed. By accumulating the jump seconds and time deviations published by the International Earth Rotation Service, and combining relativistic correction terms, the mechanical time of the solar system's center of mass is calculated. The inertial frame state vector is mapped to the geocentric-geofixed coordinate system using a rotation matrix including polar motion, precession, nutation, and Earth's rotation. For lunar-based or ground-based fixed stations, the position is further calculated using a station-centric coordinate transformation algorithm. For the Earth-Moon Lagrange point platform, an effective potential energy function incorporating both gravitational and centrifugal potential energy is established. A system of nonlinear equations is numerically solved to determine the dimensionless position, and the geocentric coordinate vector is obtained by scaling up the lunar position vector.

[0009] After determining the spatial position of the satellite and the Earth's surface, this invention constructs the observation geometry based on vector operations. Using global digital elevation model data and Earth ellipsoid parameters, the three-dimensional coordinates of the target surface pixels are calculated, constructing a radar line-of-sight vector pointing from the satellite to the target surface pixels, and calculating the surface unit normal vector based on the terrain elevation gradient. The cosine of the radar beam incident angle is obtained by performing a dot product operation between the line-of-sight vector and the surface unit normal vector, which serves as the geometric basis for subsequent shadow determination.

[0010] This invention employs a combination of geometric analysis and line-of-sight tracking for instantaneous shadow interpretation. First, self-shadowing is determined based on the cosine of the incident angle. If the cosine value is less than or equal to zero, the target surface pixel is determined to be in the self-shadowing zone of the back slope. Second, terrain-projected shadows are detected based on a radar line-of-sight propagation model. A parameterized equation for the radar line-of-sight path is established and discretely sampled. The geometrical height of the line-of-sight at the sampling points is calculated using linear interpolation and compared with the corresponding actual terrain elevation. If the actual terrain elevation is greater than the geometrical height of the line-of-sight, the target surface pixel is determined to be in the terrain-projected shadow zone.

[0011] To determine permanent shadow areas, this invention performs temporal statistics on instantaneous shadow states within an observation time window. Discrete observation times are selected within a range covering the entire orbital period, and the shadow state sequence of target surface pixels across the entire scene is calculated. The arithmetic mean of these sequences is then used as the permanent shadow criterion. When this criterion meets a set threshold condition, the target surface pixel is designated as a permanent shadow area.

[0012] Furthermore, this invention establishes a mapping relationship between matrix index and geographic space through a six-parameter affine transformation, outputting a raster image with a coordinate system; and utilizes boundary tracing and contour extraction algorithms to generate closed polygon vector features containing area and center coordinate attributes for data application.

[0013] This invention provides a model for interpreting permanent shadow areas in high-orbit SAR based on orbital information and DEM products. It offers the following advantages: 1. This invention constructs a unified orbital position calculation architecture that is compatible with various orbital description formats, including orbital elements, ephemeris data, and Earth-Moon Lagrange points. By performing time system transformations and coordinate system rotations with relativistic corrections on the ephemeris data, and establishing a potential energy function solution model for the Lagrange points, it achieves unified spatial position calculation for various high-orbit SAR platforms, including geostationary orbits, highly elliptical orbits, and deep-space orbits. This solves the compatibility problem of shadow interpretation for heterogeneous orbital data under the same model, improving the model's versatility.

[0014] 2. This invention employs a geometric analysis method based on three-dimensional spatial vectors. It determines the incident angle by calculating the dot product of the local surface normal vector and the radar line-of-sight vector, and combines this with parametric line-of-sight path tracking technology for terrain occlusion detection. This approach not only considers the slope and aspect variations of the micro-topography but also solves the problem of judging long-distance terrain occlusion through a line-of-sight propagation model. It effectively addresses the complex geometric projection distortion caused by the long observation distance and large viewing angle variations of high-orbit SAR, ensuring the geometric accuracy of shadow area interpretation under undulating terrain conditions.

[0015] 3. This invention introduces a time-series statistical mechanism covering the entire orbital cycle, transforming discrete observations of instantaneous shadow states into quantitative determinations of permanent shadow areas. By performing arithmetic averaging of shadow states throughout the entire cycle and setting a determination threshold, it can accurately identify surface areas that are long-term radar beam blind zones throughout the entire orbital cycle, and generate vectorized data products with geographic coordinates and geometric attributes, providing direct data support for mission planning, energy strategy formulation, and thermal control design of high-orbit SAR satellites. Attached Figure Description

[0016] Figure 1 This is a structural block diagram of the high-orbit SAR permanent shadow region interpretation model of the present invention; Figure 2 This is a flowchart of the high-orbit SAR permanent shadow region interpretation method of the present invention.

[0017] Among them, 10 is the orbital position calculation module; 20 is the observation geometry construction module; 30 is the instantaneous shadow interpretation module; and 40 is the permanent shadow determination module. Detailed Implementation

[0018] The technical solutions in 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 embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0019] To better understand the present invention, the above content will be described in detail below with reference to specific embodiments.

[0020] Please see the appendix Figure 1 This invention provides a high-orbit SAR permanent shadow area interpretation model based on orbit information and DEM products. The model may include: an orbit position calculation module 10, an observation geometry construction module 20, an instantaneous shadow interpretation module 30, and a permanent shadow determination module 40.

[0021] The orbit position calculation module 10 is configured to receive orbit description information from the high-orbit SAR platform. The orbit description information includes orbital elements, ephemeris data, or equivalent orbital parameters. Based on the input orbit description information, the orbit position calculation module 10 establishes a unified time system and coordinate transformation model to calculate the satellite spatial position coordinates of the high-orbit SAR platform at a specified observation time. These coordinates are typically represented as a three-dimensional vector in a geocentric coordinate system.

[0022] The observation geometry construction module 20 is connected to the orbital position calculation module 10. The observation geometry construction module 20 is configured to read global digital elevation model data and perform spatial vectorization processing on the surface raster pixels. Based on the satellite spatial position coordinates and the three-dimensional coordinates of the surface pixels, the observation geometry construction module 20 constructs the radar line-of-sight direction vector. The observation geometry construction module 20 is also used to calculate the local normal vectors of the surface pixels and the incident angle of the radar beam to determine the initial geometric visibility.

[0023] The transient shadow interpretation module 30 is connected to the observation geometry construction module 20. The transient shadow interpretation module 30 is configured to perform terrain occlusion detection on surface pixels based on a radar line-of-sight propagation model. The transient shadow interpretation module 30 traverses terrain height data along the radar line-of-sight direction to determine whether there are terrain points on the line-of-sight path that meet the occlusion conditions. If occlusion exists, the transient shadow interpretation module 30 marks the pixel as shadowed at the current moment; otherwise, it marks it as observable.

[0024] The permanent shadow determination module 40 is connected to the instantaneous shadow interpretation module 30. The permanent shadow determination module 40 is configured to perform temporal statistics on the instantaneous shadow states at multiple observation times. Based on a preset time window or orbital period, the permanent shadow determination module 40 accumulates and statistically analyzes the shadow state sequence of surface pixels. When the shadow criterion of a pixel within the statistical period meets a preset threshold condition, the permanent shadow determination module 40 determines that pixel as a permanent shadow area and generates the corresponding raster mask or vector boundary data product.

[0025] Please see the appendix Figure 2 The method mainly includes the following steps: S1. Acquire orbital information of the high-orbit SAR platform and global digital elevation model data. Orbital information includes parameters describing the motion of inclined geosynchronous orbit SAR, Earth-Moon Lagrange point platform SAR, or lunar-based SAR. Digital elevation model data is rasterized elevation data covering the target area.

[0026] S2. Establish the spatial coordinates of the high-orbit SAR platform at the observation time. Based on the type of input orbital information, select the appropriate orbital dynamics model or ephemeris interpolation algorithm to uniformly calculate the SAR platform's position into the Earth-centered Earth-fixed coordinate system. For a Lunar-Earth space platform, this step includes unified correction of the time system and coordinate transformation from the inertial frame to the Earth-fixed frame.

[0027] S3. Construct the three-dimensional observation geometry of surface pixels. For each raster pixel in the digital elevation model, convert it into a position vector in the geocentric geofixed coordinate system by combining the Earth ellipsoid parameters. Using the platform position coordinates obtained in step S2, calculate the line-of-sight vector from the SAR platform to the surface pixel, and calculate the radar incidence angle based on the local surface slope.

[0028] S4. Perform terrain occlusion detection based on the line-of-sight propagation model. For each surface pixel, establish a parameterized path equation along its corresponding radar line-of-sight direction. Retrieve the terrain elevation values ​​at each point along the line-of-sight path and determine whether the line of sight is truncated by the terrain ahead using geometric inequalities. If the line of sight is determined to be truncated, then the pixel is determined to be in a transient shadow area at the current observation time.

[0029] S5. Statistically analyze multi-temporal shadow states and determine permanent shadow areas. Select multiple discrete observation times covering the entire orbital cycle or a specific mission cycle, and repeat steps S2 to S4. Generate time-series shadow state data for surface pixels. Calculate the shadow accumulation criterion for each pixel, identify pixels that meet the long-term invisibility condition as permanent shadow areas, and output the interpretation results.

[0030] During the calculation of the location of a high-orbit SAR platform, for observation platforms described in the form of orbital elements, the orbital location calculation module 10 needs to establish a mapping relationship between orbital dynamic parameters and the geocentric coordinate system. This embodiment employs an analytical algorithm based on spherical trigonometry to directly calculate the satellite's nadir trajectory and spatial position to meet the requirements of high-precision coverage analysis.

[0031] The orbit position calculation module 10 first obtains the orbital elements information of the SAR platform, including at least the orbital inclination angle. Earth's oblateness Longitude of the sub-satellite point when the satellite passes through the ascending node angular velocity of translational motion and the rate of westward retreat of the ascending node The system calculates the time increment relative to the rising node time. We will construct parametric equations to describe the motion state of the satellite.

[0032] In calculating the geocentric latitude of satellites The module is determined according to the following formula: ; in, Indicates geocentric latitude, Represents the angular velocity of translational motion. This represents the time increment of the current moment relative to the time of the rising node. Indicates the orbital inclination angle.

[0033] Considering the influence of the Earth's ellipsoidal shape on the definition of latitude, the orbital position calculation module 10 will calculate the geocentric latitude. Convert to geographic latitude This conversion process is based on the Earth's oblateness. The following tangent relationship is satisfied after correction: ; This step ensures that the calculated latitude conforms to the definition of the geodetic coordinate system, thus maintaining spatial consistency with subsequent digital elevation model data.

[0034] In determining the change in longitude of the satellite At this time, the module utilizes spherical trigonometric transformation relationships to simultaneously calculate the sine and cosine components of the longitude variation, thereby eliminating quadrant ambiguity. Longitude variation The following system of equations must be satisfied: ; in, This represents the westward retreat rate of the ascending node relative to the Earth-fixed system (including the Earth's rotation effect). Solving the system of equations simultaneously yields... The precise value is then used. Based on this, the satellite's current geographic longitude is determined by overlaying the initial longitude. : ; Obtain the geographical longitude of the satellite's nadir point. with geographical latitude Then, combining the orbital altitude information, the orbital position calculation module 10 uniformly transforms the satellite's spatial position to the geocentric coordinate system. For an altitude of... The geocentric coordinate vector of a high-orbit SAR satellite. The result is calculated through coordinate transformation. For conventional elliptical or circular orbits, the orbital radius must be considered. (In the geocentric coordinate system) and the aforementioned latitude and longitude, perform a spherical coordinate to rectangular coordinate conversion. If a simplified homology conversion is used, the geocentric and Earth-fixed coordinates... It can be represented as: ; It should be noted that the above geographical longitudes If the Earth's rotation correction (i.e., the Earth-fixed system longitude) is already included, there is no need to superimpose the Greenwich apparent sidereal time angle again. Through the above steps, the system transforms the abstract orbital elements into three-dimensional Cartesian coordinates that can be used for geometric calculations, providing accurate radiation source location input for subsequent line-of-sight vector construction and shadow interpretation. For conventional celestial mechanics processing steps involved in orbital element calculation, such as solving Kepler's equations and converting mean anomalies to true anomalies, those skilled in the art can implement them using standard celestial dynamics algorithms, and will not be elaborated upon here.

[0035] For calculating the position of the Earth-Moon space platform based on ephemeris data, in this embodiment, the orbital position calculation module 10 uses high-precision numerical ephemeris as the basis for position acquisition to adapt to the complex motion characteristics of deep space observation nodes such as the Earth-Moon L1 / L2 point platform and the lunar base platform under a non-two-body gravitational field. Since numerical ephemeris typically uses the solar system's center-of-mass mechanical time as an independent variable, while the mission planning and data acquisition of SAR systems are usually based on Coordinated Universal Time (UTC), the orbital position calculation module 10 performs a strict time system conversion before querying the ephemeris.

[0036] The orbital position calculation module 10 first obtains the Coordinated Universal Time (UTC) at the observation time. And calculate the corresponding solar system center of mass mechanical time according to the following timescale conversion model. : ; in, This indicates the jump second time difference published by the International Earth Rotation Service, used to eliminate errors accumulated due to uneven Earth rotation. This represents a fixed deviation between Earth time and atomic time, with a value of 32.184 seconds. This is a relativistic correction term used to compensate for the periodic changes between coordinate time and original time caused by the gravitational field of the solar system.

[0037] Regarding the relativistic correction term The module uses the following sine series approximation formula for calculation: ; in, and These are empirical constant coefficients. The value is approximately 0.001657. The value is approximately 0.000022. The angle of average anomaly representing the Earth's orbit along an elliptical path. Represents the Pinghuang Classic, This represents Jupiter's mean right ascension. Through the above calculations, the system achieves time synchronization with nanosecond-level precision, ensuring the accuracy of the time reference for ephemeris lookup.

[0038] After unifying the time base, the orbital position calculation module 10 uses the calculated data... As an index, the module retrieves the state vectors of celestial bodies (including the Earth, Moon, and Sun) in the inertial frame from the DE430 ephemeris database. For lunar-based platforms, if they are located at a specific site on the lunar surface, the module also needs to combine lunar libration parameters to convert the lunar-centric inertial coordinates to lunar-fixed coordinates, and then obtain the inertial spatial position of the site through coordinate translation.

[0039] After obtaining the position vector in the inertial frame, the orbital position calculation module 10 performs a spatial coordinate transformation from the geocentric celestial reference frame to the geocentric Earth-fixed coordinate system in order to perform geometric calculations with the digital elevation model of the Earth's surface. This transformation process comprehensively considers the effects of Earth's rotation, polar motion, precession, and nutation, and its mathematical model is described as follows: ; in, This represents the transformed geocentric and geofixed coordinate vector. The inertial coordinate vector is obtained from the ephemeris query. This is a rotation matrix combining polar motion, precession, and nutation. The parameters of this matrix are calculated in real time from the Earth orientation parameter file provided by IERS. To characterize the Greenwich right ascension rotation matrix representing the Earth's rotation, let the Greenwich true sidereal hour angle at the observation time be . Then the Greenwich right ascension rotation matrix The specific form is: ; Through the aforementioned matrix operations, the orbital position calculation module 10 accurately maps the instantaneous position of the deep-space high-orbit SAR platform to a fixed Earth coordinate space consistent with the surface DEM data, providing a unified spatial reference for subsequent line-of-sight vector construction and visibility analysis. For the Earth-Moon Lagrange point platform, after obtaining the ephemeris positions of the Earth and the Moon, the module further calculates its spatial coordinates based on the dynamic equilibrium conditions of the restricted three-body problem. This part involves the introduction of the Earth-Moon mass ratio parameter and the solution of the effective potential energy gradient.

[0040] When processing spatial position calculations involving multiple celestial systems, especially for lunar-based SAR or deep space exploration platforms, the orbital position calculation module 10 needs to perform coordinate translations between different celestial body center inertial frames, as well as transformations from the station-centered coordinate system to the geocentric / lunar-centered geo-fixed coordinate system. This embodiment uses a multi-level coordinate transformation strategy to uniformly calculate orbital or station data under different reference bases into a spatial rectangular coordinate system with the Earth's center of mass as the origin.

[0041] For deep-space high-orbit SAR platforms, the raw ephemeris data is provided based on either the solar system's center-of-mass reference frame or the lunar-centered celestial reference frame. To establish the geometric line-of-sight relationship between the platform and targets on the Earth's surface, the orbital position calculation module 10 uses the vector translation principle to transform the origin of the inertial frame. The module obtains the position vector of the Earth's center of mass in the solar system's center-of-mass reference frame. and the position vector of the Moon's center of mass in the same reference frame. For a given high-orbit SAR platform, if its position is known in the lunar inertial frame, then... Or in a geocentric inertial frame of reference Then the vector transformation relationship between the coordinate systems satisfies: ; in, This represents the platform's absolute position in the solar system's center-of-mass reference frame. Through the aforementioned vector addition and subtraction operations, the system can convert the motion state centered on the moon (such as lunar orbit or a fixed point on the lunar surface) into an inertial state vector centered on the earth, and then use the earth's rotation matrix in the aforementioned embodiment to convert it to the geocentric coordinate system.

[0042] For observation nodes with fixed or quasi-fixed locations, such as lunar-based SAR stations deployed on the lunar surface or ground-based tracking stations, their positions are typically described in the local station-centered coordinate system (East-North-Sky coordinate system, ENU). To participate in global coverage analysis, the orbital position calculation module 10 performs a transformation from the station-centered coordinate system to the geocentric (or lunar-centered) geo-fixed coordinate system. Let the reference ellipsoidal longitude of the station be... Latitude is Its coordinates at the origin in the reference frame are If the coordinates of a certain observation point or sensor phase center in the station-centered coordinate system are... Then its corresponding Earth-fixed coordinates It is obtained through two rotation transformations and one translation transformation: ; in, and These are rotation matrices around the Y-axis and Z-axis, respectively, used to eliminate directional differences caused by the latitude and longitude of the station. The specific definitions of the rotation matrices are as follows: ; ; in, Based on the station's latitude, longitude, and elevation The coordinates are calculated using the standard ellipsoidal transformation formula. This transformation step applies not only to the position calculation of lunar-based radar stations but also to the coordinate unification of ground calibration stations or ground receiving stations, ensuring the geometric consistency of all relevant entities within the Earth-Moon space under the same mathematical framework. Through the aforementioned layer-by-layer transformation and unification of coordinate systems, this invention achieves a seamless transition from a local observation perspective to a global geocentric perspective, supporting the interpretation and calculation of shadow regions across celestial scales.

[0043] For SAR platforms with special orbital positions such as the Earth-Moon Lagrange point, the orbital position calculation module 10 calculates the position based on the dynamic equilibrium conditions of the restricted three-body problem. The position of such platforms no longer follows the simple Keplerian orbital law, but is determined by the combined gravitational field of the Earth and Moon and the centrifugal force generated by the system's rotation. This embodiment constructs the effective potential energy function of the Earth-Moon system, solves for the dynamic equilibrium point, and then determines the platform's real-time position in the Earth-centered Earth-fixed coordinate system.

[0044] The orbital position calculation module 10 first defines the dimensionless mass parameters of the Earth-Moon system. This parameter represents the proportion of the Moon's mass in the total mass of the Earth and Moon system, and its calculation formula is as follows: ; in, Represents Earth's mass. Represents the mass of the Moon. A dimensionless mass parameter based on standard astronomical constants. The value is approximately 0.01215. Setting the average Earth-Moon distance as a normalized unit of length and placing the origin of the coordinate system at the system's centroid, the positions of the Earth and the Moon in the rotating coordinate system are respectively located at... and Place.

[0045] In this rotating coordinate system, the module establishes an effective potential energy function describing the motion of the test particle. This function combines the gravitational potential energy of the two-body system with the centrifugal potential energy of the rotating reference frame. The specific expression is as follows: ; in, This represents the distance from a point in space to the Earth's center of mass. This represents the distance from a point in space to the Moon's center of mass. Let be dimensionless spatial coordinates in a rotating coordinate system. The Lagrange point, as the dynamic translational point in the system, satisfies the condition that the partial derivatives of the effective potential energy function in all three spatial directions are zero, i.e.: ; By numerically solving the aforementioned nonlinear equations, the orbital position calculation module 10 determines the dimensionless position coordinates of each Lagrange point. Taking the Earth-Moon L1 point as an example, this point lies on the line connecting the Earth and the Moon, and its distance from the Earth is proportional to the ratio of the Earth's distance to the Moon's distance. The calculated value is approximately 0.8369. To convert this relative position to geocentric-fixed coordinates for Earth observation analysis, the module utilizes the Moon's current position vector in the geocentric-fixed coordinate system. As a spatial reference.

[0046] Under the approximation of neglecting small perturbations in the lunar orbital eccentricity, the geocentric and Earth-fixed coordinate vector of the L1 point platform is... Calculated using the following formula: ; This calculation process leverages the high precision of lunar ephemeris data to directly map the abstract dynamic solution of the Lagrange point into physical space coordinates. For the L2 point or other translational points, the module only needs to adjust the scaling factor. Alternatively, the position can be calculated using the corresponding vector synthesis formula. This method can quickly determine the geometric position of an ultra-long-range high-orbit SAR platform relative to the Earth's surface, providing the necessary radiation source coordinate input for subsequent line-of-sight occlusion analysis.

[0047] To analyze the geometric visibility of surface targets by a high-orbit SAR platform in a unified three-dimensional space, the observation geometry construction module 20 first reads global digital elevation model data covering the target area. This data is typically stored in raster form and defines surface function relationships. Each raster cell contains geographic longitude. Geographical latitude and the corresponding geodetic elevation Since the SAR platform position output by the aforementioned orbit position calculation module 10 is described in a geocentric coordinate system, this module needs to convert the surface pixels described in geodetic coordinates into three-dimensional position vectors in a geocentric coordinate system one by one to achieve a unified coordinate reference.

[0048] In this spatial vectorization process, the observation geometry construction module 20 introduces the geometric parameters of the Earth reference ellipsoid, mainly including the semi-major axis of the ellipsoid. and first eccentricity For each DEM cell, the module first calculates the radius of curvature of the zonal circle at its latitude. This parameter reflects the curvature characteristics of the reference ellipsoid surface and is a key intermediate variable for converting geodetic coordinates to spatial rectangular coordinates. Its calculation is based on the following formula: ; in, Represents the semi-major axis of the Earth's ellipsoid. Indicates the eccentricity of the ellipsoid. Indicates the geographic latitude of the target pixel.

[0049] Obtaining the radius of curvature of the zonal circle Subsequently, the observation geometry construction module 20 combines the geodetic elevation of the pixels. Geographical longitude and geographical latitude Construct the three-dimensional position vector of surface pixels in the geocentric-geostatic coordinate system. This vector uniquely determines the precise spatial coordinates of a point on the Earth's surface relative to the Earth's center of mass; its matrix expression is as follows: ; in, It is a three-dimensional position vector. The radius of curvature of the ramusoidal circle calculated in the previous step is denoted as . Through this step, the original two-dimensional elevation raster data is transformed into a three-dimensional spatial point matrix with a rigorous analytical form, allowing the surface terrain to directly participate in subsequent line-of-sight construction and occlusion calculations based on vector algebra. For batch processing of massive DEM data, this embodiment can employ a parallel computing architecture to accelerate the above coordinate transformation process, thereby improving the efficiency of large-area interpretation.

[0050] After completing the spatial vectorization processing of surface pixels, the observation geometry construction module 20 establishes the relative geometric relationship between the high-orbit SAR platform and the surface target based on a unified geocentric coordinate system. This step calculates the propagation path direction of the radar electromagnetic waves through vector algebra operations. The observation geometry construction module 20 retrieves the position vector of the SAR platform at a specific observation time, calculated by the orbit position calculation module 10. and the three-dimensional position vectors of the surface pixels currently being processed. .

[0051] To accurately describe the geometric line of sight from the radar antenna phase center to a ground target, the module calculates the difference vector between the two points, i.e., the line-of-sight direction vector. Mathematically, this vector is defined as the target position vector minus the platform position vector. Subsequently, to ensure consistent scaling in subsequent angle calculations and step detection, the module normalizes the line-of-sight vector to obtain a unit line-of-sight vector. The above calculation process follows the following vector equation: ; in, The magnitude of the unnormalized line-of-sight vector is... This represents the straight-line slant distance between the SAR platform and the surface pixel. This represents the three-dimensional position vector of the Earth's core and Earth-fixed system for a surface pixel. This represents the geocentric and ground-fixed coordinate vector of the SAR platform. This is the normalized unit line-of-sight vector. This unit vector contains only directional information and is used to subsequently determine the angle between the radar beam and the ground surface normal vector, as well as to perform terrain occlusion retrieval along the line-of-sight direction in three-dimensional space. Through this step, the system simplifies the complex orbit-surface geometry into a standard three-dimensional vector description, providing the basic geometric input for shadow interpretation based on the principle of ray projection.

[0052] After establishing the line-of-sight direction vector, the observation geometry construction module 20 further calculates the local terrain parameters and the incident angle of the radar beam to evaluate the geometric scattering characteristics of surface pixels and potential self-shadowing. To accurately describe the impact of terrain undulations on the electromagnetic wave incident angle, this module first calculates the surface unit normal vector at the target pixel based on the elevation gradient distribution of the target pixel and its neighboring pixels in the digital elevation model. This normal vector is perpendicular to the local tangent plane and points away from the Earth's interior. Its component form in the geocentric coordinate system reflects the local slope and aspect information.

[0053] In obtaining the unit line-of-sight vector With the surface unit normal vector Based on this, the observation geometry construction module 20 uses vector dot product operations to determine the radar incident angle. The angle of incidence is defined as the angle between the direction of the radar incident beam (or the opposite of the line-of-sight direction) and the local normal to the Earth's surface. The cosine of this angle is calculated using the following formula: ; in, Indicates the local angle of incidence. This represents the line-of-sight vector per unit distance from the SAR platform to the surface pixel. This is the unit normal vector at a pixel on the Earth's surface. Because... The direction points towards the ground, while Pointing outwards from the ground, the dot product of the two is usually negative. The cosine of the incident angle can be obtained by correcting the sign by taking the negative sign.

[0054] This calculation step is crucial for determining the initial visibility of a pixel. If the calculated... Less than or equal to zero (i.e.) This indicates that the radar beam is tangential to the ground surface or incident from the back side of the ground, and the pixel is in a self-shadowed area created by the back slope. In this case, the system can directly mark the pixel as unobservable without subsequent long-distance line-of-sight occlusion detection, thus reducing the computational load of the entire scene calculation. For The system retains the pixels for the next level of processing module to perform projection shadow interpretation based on the line-of-sight propagation model.

[0055] To detect the occlusion effect (i.e., projected shadow) of distant terrain undulations on the radar beam, the instantaneous shadow interpretation module 30 constructs a line-of-sight path parameter equation connecting the phase center of the high-orbit SAR platform and the center of the ground target pixel. This process is based on the principle of geometric optics approximation, modeling the propagation trajectory of radar electromagnetic waves in non-dispersive media as a straight line segment in three-dimensional Euclidean space, thereby transforming the complex electromagnetic scattering problem into a spatial geometric intersection detection problem.

[0056] Instantaneous Shadow Interpretation Module 30 Defines Parameter Variables As the propagation distance along the line of sight, establish the spatial position vector of any point on the line of sight path. The parameterized equation describes the ray trajectory originating from the SAR platform and pointing towards a specific surface pixel. Its mathematical expression is: ; in, The position vector of the SAR platform in the geocentric coordinate system output by the orbit position calculation module 10 is used as the starting point of the line-of-sight path; The path parameter represents the straight-line distance from a point on the line-of-sight path to the SAR platform.

[0057] To ensure that the detection range covers the entire beam propagation path and does not extend beyond the surface target, parameters The range of values ​​is limited to the open interval. Internally, through this parametric modeling, the instantaneous shadow interpretation module 30 can change the path parameters. The numerical values ​​are used to traverse spatial points along the radar beam direction at a preset sampling step size to obtain a three-dimensional coordinate sequence along the line-of-sight path. These coordinate sequences are then used to index digital elevation model data to extract the corresponding terrain elevation below the line of sight, thus providing geometric input for subsequent inequality criteria. This parameterized description method is applicable to SAR observation geometry at any orbital altitude and any viewing angle, ensuring the universality of the shadow interpretation algorithm under different orbital configurations.

[0058] After establishing a parameterized model of the line-of-sight path, the transient shadow interpretation module 30 executes a point-by-point terrain occlusion detection algorithm, which can combine the synthetic aperture time information designed by the sensor. To determine whether the radar beam can reach the target pixel on the ground, the transient shadow interpretation module 30 performs discretization sampling or traversal retrieval of the digital elevation model along the line-of-sight direction. This process aims to detect whether there is a spatial intersection between the line-of-sight path and the surface terrain profile.

[0059] During the traversal, the instantaneous shadow interpretation module 30 calculates the sampling points on the line of sight path. The module establishes a quantitative relationship between the theoretical elevation of a point and the corresponding actual terrain elevation. It first indexes the Global Digital Elevation Model (BDMA) based on the spatial coordinates of the sampling point to obtain the surface terrain elevation at that location. At the same time, the module calculates the geometric height of the line of sight at the same location. If the terrain elevation at any point on the line of sight path is higher than the geometric height of the radar line of sight at that point, it is determined that the radar beam is truncated and the target pixel is in the terrain shadow area.

[0060] Based on the principle of geometric linear interpolation, the instantaneous shadow interpretation module 30 uses the following inequality as the occlusion criterion: ; in, Indicates the midpoint on the line of sight. The actual terrain elevation at the location; Indicates the elevation of the target surface pixel Indicates the elevation of the SAR platform (or the geometric height of the line-of-sight starting point); Indicates the midpoint The distance is relative to the horizontal projection distance of the target pixel or the path distance along the line of sight. It should be noted that this formula describes the geometric characteristic that the line of sight height increases linearly with distance, that is, the line of sight height should always be higher than the terrain height along the way as the distance extends from the target point to the satellite point.

[0061] If at any position during the traversal process If the above inequality condition is met, the instantaneous shadow interpretation module 30 immediately stops the subsequent search of the path and changes the current observation time. The target pixel's three-dimensional position vector Instantaneous shadow state The pixel is marked as 1 (representing a shadow state). Conversely, if no point satisfying the above conditions is found along the entire line of sight, it indicates that the radar beam is not obstructed by terrain, and the module marks the instantaneous shadow state of that pixel as 0 (representing an observable state). By detecting each pixel in the entire scene one by one, the system generates a binary instantaneous shadow distribution map for the current moment.

[0062] The permanent shadow determination module 40, based on the instantaneous shadow interpretation results obtained at a single moment, further constructs a multi-temporal shadow state sequence. Because the observation geometry of a high-orbit SAR platform (especially an inclined geosynchronous orbit or an Earth-Moon Lagrange point orbit) relative to the Earth's surface dynamically changes over time, the shadow distribution at a single moment cannot characterize the long-term observability of surface targets. Therefore, this module sets an observation time window covering the complete orbital regression cycle or a specific mission planning interval, and discretizes a series of observation moments within this window based on the Nyquist sampling theorem or the rate of change of orbital position. .

[0063] For each discrete time section The permanent shadow determination module 40 calls upon the aforementioned orbital position calculation, observation geometry construction, and instantaneous shadow interpretation functions to independently calculate the shadow distribution of all surface pixels in the entire scene. The module then calculates the three-dimensional position vector of each pixel. The interpretation results at each time point are stacked and organized in chronological order to generate a time-series state dataset for that pixel. A three-dimensional position vector for the pixel is defined. At any moment Instantaneous shadow state variable The numerical quantization rules are as follows: ; Among them, when When the value is 1, it indicates that at time 1 The surface pixel is in the radar beam blind zone due to terrain shading (projected shadow) or back slope effect (self-shadow); when the value is 0, it indicates that the pixel meets the geometric visibility conditions at the current moment. By traversing all At each sampling time, the system establishes a binary state vector for each grid cell within the target area, with the same length as the number of sampling times. This vector comprehensively records the visibility change history of the surface target throughout the entire orbital cycle, thus transforming instantaneous geometric relationships into temporal statistical characteristics. In terms of data organization, this state sequence is typically stored in a three-dimensional matrix structure to facilitate rapid temporal indexing and cumulative statistical operations in subsequent modules.

[0064] After obtaining the time-series shading status of surface pixels, the permanent shading determination module 40 performs statistical interpretation to identify permanently shaded areas. This module focuses not only on visibility at a single moment but also on evaluating the comprehensive accessibility of the target area over a specific observation period. The permanent shading determination module 40 performs temporal integration or averaging on the previously constructed binarized state sequence to calculate the three-dimensional position vector of each surface pixel. Permanent shadow criterion function This function characterizes the percentage of time a given Earth surface location is unobservable under a specific orbital configuration in Earth-Moon space, and its mathematical definition is: ; in, This represents the total number of sampling times within the observation period, corresponding to the length of the time series. Pixel 3D position vector In the Each sampling time The instantaneous shadow state value (1 indicates shadow, 0 indicates visibility). The calculated... It is a range of values The closer the value is to 1, the higher the frequency of the location being blocked during the orbital period.

[0065] After obtaining the criterion function values ​​of all scene pixels, the permanent shadow determination module 40 determines the permanent shadow based on a preset determination threshold. Pixels are classified. For strictly defined permanent shadow areas (i.e., areas never illuminated by the SAR beam throughout the entire orbital period), a threshold is set. It is usually set to 1. At this point, if the condition is met... If the pixel is in a state of permanent shadow, it is considered a permanent shadow area. In practical applications, considering track control errors, DEM data noise, or the allowance for short-term incidental visibility, the system can set a high confidence threshold close to 1 (e.g., ...). or ).when At that time, the module classifies the area as a permanent shadow area. Based on this classification result, the system generates the final permanent shadow area distribution mask or vector boundary file. This data product clearly identifies the geographical blind spots that cannot be effectively imaged and detected by a specific high-orbit SAR system, providing direct spatial data support for subsequent mission planning, orbit optimization design, or multi-baseline collaborative blind spot filling.

[0066] After completing the threshold determination, the permanent shadow determination module 40 transforms the logical interpretation result into a spatial data product with georeferenced information. The module first performs geocoding of the raster data. Using the geotransformation parameters carried by the input digital elevation model data, the system establishes a mapping relationship between the row and column indices in the result matrix and the actual geospatial coordinates. This mapping relationship is typically defined by a six-parameter affine transformation matrix, which includes the starting coordinates of the top-left corner of the image, the spatial resolution of the pixel in the X and Y directions, and rotation coefficients. By combining the binarized determination result or probability distribution matrix of the permanent shadow with this georeferenced information, the module generates a single-band raster image file with a defined coordinate system. In this output image, the pixel value of a pixel directly represents the attribute status of whether the corresponding geographical location belongs to the permanent shadow area of ​​the high-orbit SAR.

[0067] To meet the analytical needs of vector-based geographic information systems, the permanent shadow determination module 40 further performs boundary extraction processing for raster-to-vector conversion. The module employs boundary tracing algorithms or contour extraction algorithms to automatically identify the geometric boundaries between shadowed and unshaded pixels in the raster image. By connecting, smoothing, and constructing topological relationships among the discrete pixels on the boundaries, the system generates a series of closed polygonal vector features. These vector features accurately delineate the spatial outline of the permanent shadow area and are encapsulated and stored in a standard vector data format. Each vector polygon object is associated with an attribute table that records the area, center coordinates, and corresponding interpretation confidence information of the region.

[0068] Furthermore, the output data product includes a complete metadata description, clearly defining the horizontal and vertical datum, as well as the map projection method. The spatial resolution of the output data is strictly consistent with the grid size of the input digital elevation model, ensuring the geometric accuracy of the interpretation results matches the terrain data source. The generated permanent shadow distribution raster map and vector boundary files constitute the final technical output of this invention. This data can be directly loaded into space mission planning software or coverage analysis tools to assist in the selection of landing sites for deep space probes, the analysis of solar panel illumination, and the formulation of blind spot filling strategies for multi-base radar collaborative observation.

Claims

1. A high-orbit SAR permanent shadow area interpretation model based on track information and DEM products, characterized in that, The method comprises the following steps: An orbit position resolving module (10) is used to receive orbit description information, unify space-time reference, and calculate the satellite spatial position coordinates of a high-orbit SAR platform in a geocentric coordinate system at an input observation time; An observation geometry construction module (20) is connected with the orbit position resolving module (10) and is used to calculate the three-dimensional coordinates and local normal vector of a target ground surface element according to global digital elevation model data, construct a radar line-of-sight direction vector by using the satellite spatial position coordinates and the three-dimensional coordinates of the target ground surface element, and calculate the incidence angle of a radar beam based on the local normal vector of the target ground surface element and the radar line-of-sight direction vector; An instantaneous shadow interpretation module (30) is connected with the observation geometry construction module (20) and is used to determine the self-shadow state of the target ground surface element according to the incidence angle, detect the terrain obstruction on the path of the radar line-of-sight direction vector based on a radar line-of-sight propagation model, and determine the instantaneous shadow state of the target ground surface element; A permanent shadow determination module (40) is connected with the instantaneous shadow interpretation module (30) and is used to perform time series statistics on the instantaneous shadow states at multiple observation times, and determine the target ground surface element as a permanent shadow area when the statistical result meets the shadow determination threshold condition. 2.The high orbit SAR permanent shadow area interpretation model based on track information and DEM product according to claim 1, wherein, When processing orbit information in the form of orbit elements, the orbit position resolving module (10) performs the following processing: The tangent value corresponding to the geographic latitude is calculated by using a correction formula containing the Earth's flattening; The sine component and the cosine component of the longitude change amount are calculated by using a trigonometric equation set containing the ascending node regression rate and the angular velocity of the plane motion, and then the current geographic longitude of the satellite is determined; The three-dimensional rectangular coordinates of the satellite in the geocentric coordinate system are calculated by using a spherical coordinate to rectangular coordinate conversion formula in combination with the orbit height information, the geographic latitude, and the geographic longitude. 3.The high-orbit SAR permanent shadow area interpretation model based on track information and DEM product according to claim 1, wherein, When processing orbit information based on ephemeris data, the orbit position resolving module (10) performs time system conversion: The coordinated universal time at the observation time is obtained, the heliocentric system time is calculated by using a summation formula of the jump second time difference published by the International Earth Rotation and Reference Systems Service, the fixed deviation between the Earth time and the atomic time, and the relativistic correction term; The relativistic correction term is calculated by using a sine series approximation formula containing the Earth's mean anomaly, the mean ecliptic longitude, and the mean longitude of Jupiter. 4.The high-orbit SAR permanent shadow area interpretation model based on track information and DEM product according to claim 3, wherein, The orbit position resolving module (10) performs conversion from the inertial system to the geocentric coordinate system: The state vector in the inertial system is obtained by using the heliocentric system time to retrieve the ephemeris data; A combined rotation matrix of the polar motion, the precession, and the nutation, and a Greenwich right ascension rotation matrix representing the Earth's rotation are constructed; The state vector in the geocentric coordinate system is calculated by using a matrix multiplication operation formula; For the lunar base and the ground fixed station, the position in the geocentric coordinate system is further calculated by using a station-centered coordinate conversion algorithm containing the longitude rotation matrix and the latitude rotation matrix. 5.The high-orbit SAR permanent shadow area interpretation model based on track information and DEM product according to claim 1, wherein, For the Lagrange point platform, the orbit position resolving module (10) performs the following processing: A dimensionless mass parameter of the Earth-Moon system is defined, and an effective potential function is established, which contains the potential energy of the binary gravitational system and the centrifugal potential energy of the rotating reference frame; The dimensionless position coordinates of the Lagrange points are determined by solving the nonlinear equations whose partial derivatives of the effective potential function in three spatial directions are zero through a numerical algorithm; The position vector of the Moon in the Earth-Centered Earth-Fixed coordinate system is used as a spatial reference, and the dimensionless position coordinates of the Lagrange points are calculated by vector scaling operation formula combined with the scaling factor of the Lagrange points relative to the Earth-Moon distance. 6.The high-orbit SAR permanent shadow area interpretation model based on track information and DEM product according to claim 1, wherein, When calculating the three-dimensional coordinates of the target ground surface element and constructing the radar line-of-sight direction vector, the observation geometry construction module (20) performs the following steps: The colure radius at the latitude of the target ground surface element is calculated by the curvature radius formula containing the long semi-axis and the eccentricity of the Earth ellipsoid; The position vector of the target ground surface element in the Earth-Centered Earth-Fixed coordinate system is calculated by the conversion formula from geodetic coordinates to Earth-Centered Earth-Fixed coordinates combined with the colure radius, the elevation value, the longitude and the latitude of the target ground surface element; The unit line-of-sight vector from the satellite to the target ground surface element is calculated by vector subtraction and normalization calculation formula.

7. The orbit information and DEM product based high-orbit SAR permanent shadow region interpretation model according to claim 6, characterized in that, The observation geometry construction module (20) calculates the local normal vector and the incidence angle of the radar beam in the following way: The ground surface unit normal vector of the target ground surface element is calculated by the elevation gradient calculation algorithm according to the elevation gradient distribution of the target ground surface element and its neighboring elements in the digital elevation model; The cosine value of the radar incidence angle is calculated by the dot product operation formula of the unit line-of-sight vector and the ground surface unit normal vector of the target ground surface element; If the cosine value is less than or equal to zero, it is directly determined that the target ground surface element and its neighboring elements are in the self-shadow area. 8.The high orbit SAR permanent shadow area interpretation model based on orbit information and DEM product according to claim 1, wherein, The specific way of detecting terrain occlusion by the instantaneous shadow interpretation module (30) is as follows: A parameterized equation of the radar line-of-sight path is established, which describes the spatial point positions along the line-of-sight direction from the satellite spatial position coordinates; The geometric height of the sampling points is calculated by the linear interpolation formula, and the actual terrain elevation corresponding to the sampling points is obtained by traversing the sampling points along the line-of-sight path; Whether the target ground surface element is in the terrain shadow area is calculated by the inequality criterion of comparing the actual terrain elevation with the line-of-sight geometric height; The line-of-sight geometric height is calculated by the linear interpolation formula containing the elevation of the target ground surface element, the height difference between the line-of-sight start and end points, and the path projection component. 9.The high orbit SAR permanent shadow area interpretation model based on orbit information and DEM product according to claim 1, wherein, The permanent shadow determination module (40) performs timing statistics in the following way: A series of observation times are discretely selected within the observation time window covering the complete orbital period, and for each time, the instantaneous shadow state variable of each target ground surface element in the full scene is generated, where the shadow state is assigned a first value and the non-shadow state is assigned a second value; The permanent shadow criterion function is calculated by the arithmetic mean value calculation formula, which is the arithmetic mean value of the instantaneous shadow state variables of the target ground surface element at all observation times; When the value of the permanent shadow criterion function is greater than or equal to the shadow determination threshold, the target ground surface pixel is classified as a permanent shadow area. 10.The high-orbit SAR permanent shadow area interpretation model based on track information and DEM product according to claim 1, wherein, The permanent shadow determination module (40) is further configured to generate a data product: A mapping relationship between a result matrix row and column index and a geospatial coordinate is calculated by a six-parameter affine transformation algorithm, and a raster image file with a coordinate system is generated; A geometric boundary line between the permanent shadow area and the non-permanent shadow area is calculated by a boundary tracking algorithm and an isopleth extraction algorithm, a closed polygon vector element is generated, and an area and a center coordinate attribute are associated.

Images