Computer vision based remote sensing image processing system

Abstract

This invention discloses a computer vision-based remote sensing image processing system, belonging to the field of remote sensing image processing technology. It includes: a multi-source data spatiotemporal registration and parameter analysis module; a land surface illumination prior construction module; a variational energy functional construction module; a partial differential equation iterative solution module; and a geoscientific threshold classification module. In this invention, a rational polynomial model is used to complete orthorectification and spatial registration of multi-source remote sensing images and elevation data. Combining the land surface normal vector matrix and the solar illumination direction vector at the imaging time, a physical illumination prior that fits the real terrain imaging environment is constructed, distinguishing between direct illumination, diffuse reflection components, and terrain-occluded shadow areas. Then, a multi-constraint optimization objective is constructed through a logarithmic domain variational energy functional, and the alternating direction multiplier method is used for iterative solution to accurately separate the illumination field from the real land surface reflectivity.

Description

Technical Field

[0001] This invention relates to the field of remote sensing image processing technology, and more particularly to a remote sensing image processing system based on computer vision. Background Technology

[0002] With the rapid development of remote sensing satellite technology, high-resolution optical remote sensing imagery has been widely applied in core areas of national economy and people's livelihood, such as land surveys, geological disaster monitoring, ecological environment assessment, and agricultural and rural management. However, under complex terrain conditions such as cloudy skies and mountainous areas, remote sensing images generally suffer from spectral reflectance distortion caused by topographic relief, uneven illumination, and atmospheric scattering. Specifically, the reflectance of objects on the sunny side is overestimated, while the reflectance of objects on the shady side and in valley shadow areas is severely underestimated, resulting in the phenomenon of "different spectra for the same object and the same spectrum for different objects," which greatly limits the interpretation accuracy and application value of remote sensing images.

[0003] In existing technologies, atmospheric correction is performed using a radiative transfer model, and topographic correction is performed using DEM data. However, this approach has shortcomings: First, most of them use affine transformation for image and DEM registration, which has serious projection difference in undulating mountainous areas. Insufficient registration accuracy leads to the failure of subsequent correction. Secondly, it can only achieve simple slope correction and cannot accurately simulate complex lighting effects such as terrain-projected shadows and terrain mutual reflection, nor can it decouple lighting from the actual surface reflectivity.

[0004] For example, patent application CN119624781B discloses a spatiotemporal super-resolution method for remote sensing images based on surface process modeling. This method involves interpolating low-resolution surface remote sensing images using time- and space-based interpolation to obtain interpolated surface remote sensing images. It then constructs control equations to describe spatiotemporal changes in the surface and solves these equations using LASSO regression to obtain differential equations for calculating high-resolution predicted images. Finally, it calculates high-resolution predicted images for a specified time and spatial location based on these differential equations. However, the initial images obtained in this method cannot guarantee that they effectively reflect illumination and the true reflectivity of the surface. Summary of the Invention

[0005] To overcome the problems of the prior art, the present invention aims to provide a computer vision-based remote sensing image processing system. This system constructs a physical illumination prior that conforms to the real terrain imaging environment by using the surface normal vector matrix and the sunlight direction vector at the imaging time, distinguishing between direct illumination components, diffuse reflection components, and terrain illumination interference. Then, it constructs a multi-constraint optimization objective by using variational energy functionals and iteratively solves the problem using the alternating direction multiplier method, thereby improving the registration accuracy of remote sensing images and effectively decoupling illumination from the real surface reflectance.

[0006] To achieve the above objectives, the present invention adopts the following technical solution: Computer vision-based remote sensing image processing systems include: The multi-source data spatiotemporal registration and parameter analysis module is configured to: use rational polynomials to orthorectify and resample the original remote sensing images and digital elevation model data to generate orthorectified images; calculate the surface normal vector matrix and analyze the solar illumination direction vector at the imaging time; The surface illumination prior construction module is configured to: calculate the direct illumination component and the total diffuse reflection component based on the surface normal vector matrix and the solar illumination direction vector, and superimpose them to generate a geophysical radiation prior penalty term. The variational energy functional construction module is configured to: treat the orthophoto as the product of reflectivity and illumination field and map it to the logarithmic domain, and construct a variational energy functional that includes a data fidelity term, a reflectivity smoothing term based on the first norm, an illumination field continuity term based on the second norm, and the geophysical radiation prior penalty term. The partial differential equation iterative solution module is configured to: use the alternating direction multiplier method to iteratively solve the variational energy functional, update the illumination field and reflectivity in the frequency domain, update the boundary using the soft shrinkage threshold operator, and after convergence, perform inverse mapping to output the surface physical reflectivity matrix that has removed terrain illumination interference. The geoscience threshold classification module is configured to: calculate a multispectral feature index based on the surface physical reflectance matrix, and after global threshold binarization, output a vector polygon representing the spatial topological relationship through morphological optimization.

[0007] Preferably, the step of generating orthorectified images by orthorectifying and resampling the original remote sensing images and digital elevation model data using rational polynomials includes: Input data includes single-band or multi-band raw remote sensing image data. and digital elevation models covering the same geographical area ,in, For remote sensing image column numbers, For remote sensing image row numbers, Longitude Latitude; Construct three-dimensional coordinates of normalized ground longitude, latitude, and elevation. Two-dimensional pixel coordinates of the normalized values ​​of the row number and column number of the original remote sensing image. The mapping relationship; With digital elevation model The target image grid is generated based on the geographic grid, and any geographic coordinate point in the target image grid is extracted. Corresponding elevation value ; Will After scale normalization, the geographic points in the original remote sensing image data are calculated. The corresponding two-dimensional pixel coordinates A bicubic convolution interpolation algorithm is used, utilizing two-dimensional cell coordinates. The grayscale value of the target point is calculated by resampling the surrounding integer pixel values, and an orthophoto with spatial consistency with the target image grid is generated based on the grayscale value. .

[0008] Preferably, the process of calculating the surface normal vector matrix and resolving the sunlight direction vector at the imaging time is as follows: Digital Elevation Model The Gaussian smoothing operator is applied, and then the elevation gradient at each spatial location is solved. , ; Calculate and After obtaining the partial derivatives of elevation in the direction, a digital elevation model is constructed based on the definition of the normal vector of a space surface in calculus. The unnormalized 3D normal vector at each grid point in the array ; By calculating the three-dimensional normal vector The modulus length is used to obtain the unit surface normal vector matrix at each pixel. ; From raw remote sensing image data The latitude and longitude of the image center point were extracted. and imaging Coordinated Universal Time ; Based on imaging Coordinated Universal Time The Julian Day is calculated, and then the declination angle of the sun at the moment of image formation is determined. Sum of time angles ; Combined with the latitude of the target area Solve for the solar altitude angle. and solar azimuth ; and convert it into a solar illumination direction vector based on the geographic northeast-sky coordinate system. .

[0009] Preferably, based on the surface normal vector matrix and the solar illumination direction vector, the direct sunlight component and the total diffuse reflection component are calculated, as follows: By analyzing the normal vector matrix and the direction vector of sunlight Perform pixel-by-pixel dot product operations to solve for the initial direct cosine coefficients. ;when At that time, its direct sunlight illuminance was set to 0; For digital elevation models Each cell in The system emits a ray in three-dimensional space along the direction opposite to the sun's rays. , These are the step parameters; In each step, calculate the current ray. Projected coordinates on a two-dimensional plane And query the location in the digital elevation model. The corresponding actual surface elevation ; If there is any step point within the set maximum tracking distance satisfy Then the pixel Being physically occluded by terrain and placed in cast shadow; generate a Boolean occlusion mask matrix. 1 represents an unobstructed value, and 0 represents an obstructed value. Based on the direct cosine coefficient Boolean occlusion mask matrix Obtain the direct light component ,in, This is the solar constant at the top of the atmosphere. The transmittance parameter of the entire atmosphere at the time of imaging; In this pixel At that location, multiple probe rays are emitted into the three-dimensional hemispherical space at the azimuth angle; Define azimuth sampling interval ; along the first Each azimuth direction The indicated horizontal direction in the digital elevation model The search proceeds from top to bottom, calculating the elevation angles of all surface undulations relative to the center point, and recording the maximum elevation angle. Solve for the sky view factor. ; This is the total number of azimuth angle samples; In obtaining the sky view factor Then, calculate the ambient diffuse light component. , where constant This represents the energy of diffused light over flat ground. By introducing an approximate cross-reflection term gain, the corrected total diffuse reflection component is obtained. ,in, This is an estimate of the average albedo of the background area. The average direct irradiance of the surrounding terrain.

[0010] Preferably, the direct light component With the total diffuse component By performing linear superposition, we obtain the geophysical radiation prior penalty term. .

[0011] Preferably, the variational energy functional is a global energy functional, expressed as: ; in, For data fidelity items, For logarithmic domain reflectance, For logarithmic domain illumination field; This is a reflectance smoothing term based on the first norm; For the illumination field continuity term based on the 2-norm; For illumination field Anchoring Geophysical Radiation Prior Penalty ; , , This is the preset Lagrange hyperparameter.

[0012] ; for Logarithmic domain reflectance The spatial gradient vector; Represents the logarithmic domain reflectance along Partial derivatives in direction; Represents the logarithmic domain reflectance along Partial derivatives in direction; ; for Logarithmic domain illumination field at the location The spatial gradient vector; Represents the logarithmic domain illumination field along Partial derivatives in direction; Represents the logarithmic domain illumination field along Partial derivatives in direction; The spatial domain representing the entire remote sensing image. , This represents a two-dimensional spatial area element.

[0013] These are prior constraints on the illumination field, used to constrain the logarithmic domain illumination field to be solved. Logarithmic field form approaching the geophysical radiation prior field , Geophysical radiation a priori penalty term The logarithmic field form.

[0014] Preferably, the partial differential equation iterative solution module further includes constructing an augmented Lagrangian function, used to transform the variational energy functional. The minimization problem is transformed into an unconstrained subproblem that can be solved iteratively using the alternating direction multiplier method; The variational energy functional There exists a non-differentiable norm term in it. ; Introducing auxiliary variables By substitution, the non-differentiable gradient term is removed from the reflectivity term. Separation from the middle; ; Auxiliary variable The two components; Therefore, we construct the augmented Lagrangian function. : in, For the scaled dual variable, The penalty parameters are set. For gradient operators, For orthophotos Logarithmic domain observation image obtained after logarithmic mapping.

[0015] Preferably, the iterative solution process for the augmented Lagrangian function is as follows: Extracting the augmented Lagrangian function that is only related to the logarithmic domain illumination field Related items, constructing the updated values ​​of iterative variables Logarithmic domain illumination field Find the variational derivative and set it equal to zero to obtain the light field in the logarithmic domain. The Euler-Lagrange differential equation is solved by two-dimensional discrete Fourier transform to obtain the updated values ​​of the iterative variables. ; This is the current iteration number; Extracting the logarithmic domain reflectance to be solved Related terms, construct sub-problems Solve the Euler-Lagrange equations for it; rearrange to obtain the equations containing the divergence operator. and Laplace operator The equations are given; subproblems are also given. The analytical update equation; Determine the form of the subproblem Based on the proximal gradient algorithm, a classic soft contraction threshold operator from one-dimensional signal processing is introduced, and a closed-form solution is given at the matrix level. , ; Scaled dual variable Perform gradient ascent updates to accumulate violations of equality constraints. historical error , ; At the end of each iteration, the principal residual is calculated. and dual residual When both are less than the set minimum physical tolerance value The partial differential equation system is considered convergent when the preset maximum number of analytical iterations is reached. Finally, regarding the logarithmic domain reflectance Perform an exponential inverse mapping to output the surface physical reflectance matrix. .

[0016] Preferably, the geoscience threshold classification module receives the surface physical reflectance matrix to calculate the multispectral feature index, and the specific steps of global threshold binarization are as follows: Based on the surface physical reflectance matrix after removing topographic lighting interference, multispectral feature indices are reconstructed pixel by pixel. The multispectral feature indices include at least the normalized vegetation index and the normalized water index. The system calls a pre-defined set of global rigid thresholds, performs Boolean judgments on the multispectral feature indices, and generates binarized target masks for various land features.

[0017] Preferably, after generating the binarized target mask, the specific steps for morphologically optimizing the output vector polygon are as follows: An isotropic cross-shaped kernel with a set size of 3×3 is used as the morphological structural element to perform mathematical morphological opening and closing operations on the binarized target mask. The topology tracing algorithm is invoked to convert the rasterized binary target mask, which has been processed by mathematical morphology opening and closing operations, into a vector polygon containing spatial topological relationships, and then output it with the corresponding spatial reference frame.

[0018] In summary, due to the adoption of the above technical solution, the beneficial effects of the present invention are: 1. This invention completes orthorectification and spatial registration of multi-source remote sensing images and elevation data through a rational polynomial model. It combines the surface normal vector matrix and the sunlight direction vector at the imaging time to construct a physical illumination prior that fits the real terrain imaging environment, distinguishing between direct illumination, diffuse reflection components and terrain occlusion shadow areas. Then, it constructs a multi-constraint optimization objective through a logarithmic domain variational energy functional and uses the alternating direction multiplier method to solve iteratively, accurately separating the illumination field from the real surface reflectivity.

[0019] 2. This invention, based on a high-precision physical reflectance matrix that has removed terrain and lighting interference, calculates multispectral feature indices and combines global threshold binarization and morphological optimization to avoid the misclassification and omission of ground features caused by lighting distortion and shadow interference in traditional methods. This significantly improves the accuracy of ground feature classification and the precision of ground feature boundaries. Simultaneously, it can directly output vector polygon results with spatial topological relationships and corresponding spatial references, making it adaptable to various geographic information system application scenarios and reducing the manual cost of remote sensing image interpretation.

[0020] In summary, this invention can accurately separate the illumination field from the actual surface reflectance, significantly improve the efficiency of remote sensing image processing, greatly enhance the accuracy of land cover classification and the accuracy of land cover boundaries, and is adaptable to various geographic information system application scenarios, reducing the manual cost of remote sensing image interpretation. Attached Figure Description

[0021] Figure 1 This is a system structure diagram of the present invention. Detailed Implementation

[0022] Several embodiments of this application will now be described in more detail with reference to the accompanying drawings to enable those skilled in the art to implement this application. This application may be embodied in many different forms and for various purposes and should not be limited to the embodiments set forth herein. These embodiments are provided to make this application thorough and complete, and to fully convey the scope of this application to those skilled in the art. The embodiments described do not limit this application.

[0023] Unless otherwise defined, all terms used herein (including technical and scientific terms) shall have the same meaning as commonly understood by one of ordinary skill in the art to which this application pertains. It will be further understood that terms such as those defined in commonly used dictionaries shall be interpreted as having a meaning consistent with their meaning in the relevant field and / or the context of this specification, and shall not be interpreted in an idealized or overly formal sense unless expressly defined herein.

[0024] Example 1 Its specific implementation method is combined with the appendix Figure 1 Please provide a detailed explanation.

[0025] Appendix Figure 1 The diagram below shows the structure of a computer vision-based remote sensing image processing system according to an embodiment of the present invention. It illustrates the connection between the multi-source data spatiotemporal registration and parameter analysis module and the geoscientific threshold classification module, and marks the main functional interaction flow of each module.

[0026] In this embodiment, the computer vision-based remote sensing image processing system includes: Module 1, Multi-Source Data Spatiotemporal Registration and Parameter Analysis Module, is configured as follows: using rational polynomials to orthorectify and resample the original remote sensing images and digital elevation model data to generate orthorectified images; calculating the surface normal vector matrix and analyzing the solar illumination direction vector at the imaging time; If there is a subpixel-level deviation between the image's pixels and its corresponding digital elevation model (DEM) elevation points, or if the solar angle at the time of imaging is calculated incorrectly, it will directly cause the radiative transfer equation to diverge in subsequent solutions.

[0027] In existing technologies, multi-source data registration often employs affine transformation, which can produce severe projection differences in undulating mountainous areas.

[0028] Input data includes uncorrected Level-1 single-band or multi-band raw remote sensing image data. and external high-precision digital elevation models covering the same geographical area. ,in, For remote sensing image column numbers, For remote sensing image row numbers, Longitude Latitude; At the same time, the general rational polynomial coefficient set included in the remote sensing image header file is read.

[0029] First, the system constructs normalized three-dimensional coordinates of ground longitude, latitude, and elevation using the RPC forward equation. Two-dimensional pixel coordinates derived from the normalized values ​​of the row and column numbers of the original remote sensing image. Mapping relationship: ; ; In the formula, , , , These represent third-order polynomials composed of 20 RPC polynomial coefficients.

[0030] To achieve imagery and digital elevation model For registration, the system executes the back projection resampling algorithm.

[0031] Specifically, the system uses a digital elevation model The target image grid is generated based on the geographic grid. For any geographic coordinate point in the target image grid... The system is based on the digital elevation model. Extract its corresponding elevation value ; Then After scaling, input the above RPC formula ( ; ;), which inversely calculates the geographic points in the original remote sensing image data. The corresponding precise floating-point level two-dimensional cell coordinates ; Since the calculated coordinates are usually non-integer, the system uses a bicubic convolution interpolation algorithm, utilizing two-dimensional cell coordinates. The grayscale value of the target point is calculated by resampling the surrounding 16 integer pixel values.

[0032] At this point, an orthophoto with spatial consistency with the target image grid has been generated. In this image, coordinates It represents both the image pixel index and directly corresponds to the digital elevation model. A grid cell; Solution of the microscopic three-dimensional normal vector field of the Earth's surface based on the central difference method: To calculate the cross-sectional area of ​​each pixel receiving direct sunlight, the normal vector of the land surface corresponding to each pixel is calculated, which falls under the category of geographic differential geometry. The system is based on the registered actual elevation of the land surface. Calculate the partial derivatives of the elevation of the Earth's surface in the horizontal and vertical directions; To avoid the impact of sudden elevation changes caused by local noise on the stability of partial derivative calculations, the system first optimizes the digital elevation model. The Gaussian smoothing operator is applied, followed by the central finite difference method to solve for the elevation gradient at each spatial location. , : ; ; In the formula, and These represent the actual physical resolution of the DEM grid in the longitudinal and latitudinal directions (e.g., 10 meters or 30 meters).

[0033] Calculate and After calculating the partial derivatives of elevation in the direction, the system constructs a digital elevation model based on the definition of the normal vector of a spatial surface in calculus. The unnormalized 3D normal vector at each grid point in the array : .

[0034] Subsequently, the three-dimensional normal vector was calculated. The modulus length is used to obtain the normalized unit surface normal vector matrix at each pixel. : ; The normal vector matrix It is stored as a floating-point three-channel tensor with the same dimensions as the original image, and its components represent the precise orientation of each surface micro-element pointing towards the zenith in three-dimensional Euclidean space.

[0035] The system uses raw remote sensing image data Two core temporal and spatial parameters were extracted from the image: the latitude and longitude of the image center point. And imaging Coordinated Universal Time (UTC time) .

[0036] Based on imaging Coordinated Universal Time Using astronomical calendar and physical formulas, the Julian Day is first calculated, and then the declination angle of the sun at the moment of image formation is determined. Sum of time angles Combined with the latitude of the target area The system solves for the solar altitude angle. and solar azimuth .

[0037] In practical implementation, in order to perform dot product operations with subsequent three-dimensional spatial vectors, the system needs to convert the solar altitude angle and azimuth angle in the spherical coordinate system into the solar illumination direction vector based on the geographic northeast-sky coordinate system (North-East-Down system converted to the corresponding Euclidean space). .

[0038] The conversion formula is: , , , ; At this point, Module 1 has completed all the preliminary preparations and is ready to output to the next module: orthophotos. Surface normal vector matrix and the direction vector of sunlight at the moment of imaging. .

[0039] Module 2, Surface Illumination Prior Construction Module, is configured as follows: Based on the surface normal vector matrix and the solar illumination direction vector, calculate the direct illumination component and the total diffuse reflection component, and superimpose them to generate a geophysical radiation prior penalty term; This module is used to distinguish between real low-reflectivity materials on the ground surface (such as black asphalt roads or coal mine surfaces) and shadowed areas formed by the occlusion of high-reflectivity materials. Existing techniques typically use image grayscale histograms to threshold and extract shadows, but this is a data-empirical approach without a physical basis.

[0040] This embodiment abandons the use of orthophotos. Instead of relying on its own observations, the system utilizes the spatial geometric parameters output by Module 1 to deduce a theoretical prior illumination distribution field through radiometric laws. .

[0041] The total solar radiation received by the Earth's surface can be decomposed into two parts: the direct irradiance component from direct sunlight, and the diffuse irradiance component from atmospheric scattering and reflection by surrounding terrain. This module solves for these two parts independently in steps.

[0042] Calculation of direct illumination field based on Lambert's cosine law and ray tracing: First, the attenuation of effective light-receiving area due to terrain undulations is calculated. According to Lambert's cosine law in radiometry, the irradiance received by an ideal diffuse surface is proportional to the cosine of the angle between the incident light direction and the normal to the ground surface. The system calculates this by adjusting the normal vector field matrix. and the direction vector of sunlight Perform pixel-by-pixel dot product operations to solve for the initial direct cosine coefficients. : ; when When the light source is at a certain time, it indicates that the light is incident from the back of the micro-element on the ground, meaning that the micro-element is in a "self-shaded" area, and its direct light intensity is set to 0.

[0043] However, the above calculations cannot identify projected shadows, that is, although the normal of a certain point is pointing towards the sun ( However, the light beam is blocked by a higher mountain in front of it. To address this, the system introduces a rigid physical ray tracing algorithm based on DEM.

[0044] For digital elevation models Each cell in The system is along the direction opposite to the sun's rays (i.e., vector). A ray is emitted in three-dimensional space from the direction of ( ), and the equation of the ray is: ; ; ; in For step parameters, The set tracking step size (usually half the DEM resolution to ensure that it does not penetrate the mountain grid).

[0045] In each step, calculate the current ray. Projected coordinates on a two-dimensional plane And by using bilinear interpolation, the location in the digital elevation model can be found. The corresponding actual surface elevation .

[0046] The system executes the collision detection logic as follows: if any step point exists within the set maximum tracking distance (e.g., the tracking distance reaches the atmospheric boundary), satisfy This proves that the ray penetrated underground in three-dimensional space, meaning that the pixel... The object is physically occluded by the terrain and is in a cast shadow. The system generates a Boolean occlusion mask matrix. 1 represents an unobstructed value, and 0 represents an obstructed value.

[0047] Then, the direct light component The analytical expression is: ; In the formula, This is the solar constant at the top of the atmosphere. This is the transmittance parameter of the entire atmospheric layer at the time of imaging (which can be read from the default values ​​of the general MODTRAN atmospheric model).

[0048] Continuous integral estimation of diffuse scattering field based on sky-view factor (SVF): Features within the shadow area are not completely black because they still receive scattered light from the sky dome (such as Rayleigh and Mie scattering). However, pixels at the bottom of deep valleys receive far less diffuse light than pixels at the top of mountains because the surrounding mountains reduce the proportion of visible sky.

[0049] The system uses the Sky View Factor (SVF) to quantify this physical phenomenon. SVF is defined as the ratio of the solid angle of the sky visible from a point (not obscured by terrain) to the solid angle of the hemisphere on open flat ground, with a value range between [0,1].

[0050] To calculate SVF, the system needs to be in this cell At this location, multiple probe rays are emitted into the three-dimensional hemispherical space of the azimuth angle; the system defines the azimuth angle sampling interval. (For example, a probe line is launched every 10 degrees, for a total of 36 directions). The system along the... Each azimuth direction The indicated horizontal direction in the digital elevation model The search proceeds from top to bottom, calculating the elevation angles of all surface undulations relative to the center point, and recording the maximum elevation angle. The maximum elevation angle represents the limit of how far the mountain can obscure the skyline in that direction.

[0051] The system solves for the sky view factor by discrete numerical integration of the azimuth angle. ; among them, of This is the total number of azimuth angle samples, and its value is determined by the azimuth angle sampling interval. Decide; In obtaining the sky view factor Then, the system calculates the ambient diffuse light component. Assume that the energy of diffused light over a flat surface is constant under ideal conditions. The diffuse field model adjusted by terrain is then: ; Considering that the surrounding mountains, after receiving sunlight, will reflect the detection rays back to this pixel, (Terrain mutual reflection), the system introduces an approximate mutual reflection term gain. The mutual reflection factor of adjacent terrain is... Proportional; therefore, the total diffuse component is corrected to: ; in, This is an estimate of the average albedo of the background area. The average direct irradiance of the surrounding terrain.

[0052] Finally, the system derives the direct light component from purely physical derivations. With the total diffuse component By performing linear superposition, the final geophysical radiation prior penalty term is obtained. : ; To ensure that the data volume is consistent with the original image grayscale in the subsequent partial differential energy functional, the system performs... Perform global maximum and minimum value normalization to map the values ​​to the same range as the image pixels (e.g., a floating-point number of [0, 255] for an 8-bit image, or directly map it to a relative illumination coefficient tensor of [0.0, 1.0].

[0053] Module 3, Variational Energy Functional Construction Module, is configured as follows: the orthophoto is regarded as the product of reflectivity and illumination field and mapped to the logarithmic domain, and a variational energy functional is constructed including a data fidelity term, a reflectivity smoothing term based on the first norm, an illumination field continuity term based on the second norm, and the geophysical radiation prior penalty term. After completing the one-to-one mapping of physical quantities and the forward evolution of the a priori field in Modules 1 and 2, the task of Module 3 is to transform the typical ill-conditioned inverse problem of extracting the true reflectance from a single observation image into a mathematical variational extremum problem with a unique global optimum or a stable local optimum.

[0054] The system establishes a fundamental physical observation model based on Retinex theory and the Lambert reflection hypothesis in computer vision: The orbital radiance received by a remote sensing sensor in a specific band (i.e., the observed value after radiometric calibration of the original image). After stripping away atmospheric path radiation, it can be rigorously equivalent to the true spectral reflectance of the Earth's surface. Combined with the total illumination energy received at that pixel The pixel-by-pixel product. That is: ; Because multiplicative noise easily leads to gradient explosion in mathematical solutions and is difficult to solve linearly, the system first performs a logarithmic domain mapping on the above observation model. Let , , Then the multiplicative physical model is transformed into an additive linear model: ; In order to start from a known single variable Two unknown independent variables were solved. and The system introduces strong prior constraints based on geophysical laws.

[0055] The system constructs a global energy functional consisting of four energy sub-terms. : ; The purpose of this functional is to find a set of The combination of these factors makes the global energy functional... It reaches the preset minimum value. The following describes the global energy functional. Explanation of each item: For data fidelity items, For logarithmic domain reflectance, For the illumination field: This ensures that the decomposed reflectivity and illumination field, after recombination, do not deviate from the actual satellite observations. The system uses the squared L2 norm to penalize the residuals: ; In the formula, The spatial domain of the entire remote sensing image represents the effective pixel region that participates in the variational energy functional calculation. , This represents a two-dimensional spatial area element. The reason for using the L2 norm instead of the L1 norm is that it is assumed that the quantization noise of the satellite sensor follows a Gaussian distribution, and the L2 norm can provide the best Gaussian noise smoothing capability.

[0056] For the reflectivity smoothing term based on the 1-norm: Key constraints regarding surface physical properties. In real-world surface environments (excluding water bodies), the spectral reflectance should be uniform within the same geological structure or crop cover area. Reflectance values ​​should only exhibit abrupt changes at the boundaries of different land features (such as the boundary between grassland and road).

[0057] To mathematically describe this step function characteristic of being internally flat and having steep edges, the system abandons the L2 norm, which causes edge blurring, and instead adopts the L1 norm of the first-order spatial gradient, i.e., the total variational regularization term: ; for Logarithmic domain reflectance The spatial gradient vector; Represents the logarithmic domain reflectance along Partial derivatives in direction; Represents the logarithmic domain reflectance along Partial derivatives in direction; The introduction of the first norm is crucial, as it allows reflectance gradients to have extremely non-zero values ​​at a few boundaries (i.e., preserving sharp physical boundaries) while forcing gradients to zero in other regions (i.e., eliminating small spectral fluctuations within similar land features).

[0058] For the continuity term of the illumination field based on the 2-norm: In stark contrast to reflectivity, the distribution of light in nature (especially diffuse scattering and the transition from direct sunlight on gentle slopes) is spatially continuous and slowly varied. Except for the edges of cast shadows created by steep cliffs, the light field should not contain high-frequency abrupt changes.

[0059] Therefore, for the illumination field The spatial gradient is penalized using the L2 norm, i.e., Tikhonov regularization: ; for Logarithmic domain illumination field at the location The spatial gradient vector; Represents the logarithmic domain illumination field along Partial derivatives in direction; Represents the logarithmic domain illumination field along Partial derivatives in direction; The L2 norm severely punishes any sharp abrupt changes in the illumination field, forcing the calculated... The formation of a gently curved surface conforms to the physical propagation laws of photon scattering and reflection on macroscopic terrain.

[0060] For geophysical radiation a priori penalty terms: Traditional intrinsic image decomposition often suffers from severe semantic confusion due to the lack of references, resulting in the allocation of reflectance and illumination (e.g., identifying black coal mines as deep shadows).

[0061] The system will use the geophysical radiation prior penalty term derived from module two based on elevation and solar angle in a purely positive direction. Introduce a functional. Then take its logarithm. The system mandates that the illumination field be calculated mathematically. Anchored to the logarithm of the geophysical radiation prior penalty term in the global trend. : ; at last, , , This refers to the Lagrange hyperparameter. In this embodiment, it is automatically assigned based on the physical reliability of the input data: when the input DEM resolution is extremely high and the imaging weather is clear, it is increased by a preset amount. The weighting of the physical priors allows the system to highly trust the physical priors; conversely, if the DEM is coarser, the weighting is reduced. improve This makes the system more reliant on the TV smoothness of the image itself.

[0062] Module 4, Partial Differential Equation Iterative Solution Module, is configured as follows: Iteratively solve the variational energy functional using the alternating direction multiplier method, update the illumination field and reflectivity in the frequency domain, update the boundary using the soft contraction threshold operator, and after convergence, perform inverse mapping to output the surface physical reflectivity matrix that has removed terrain illumination interference. After constructing the global energy functional Then, the problem is transformed into a pure optimization solution. However, due to the global energy functional... There exists a non-differentiable norm term in it. This means that extreme value problems cannot be solved directly using the conventional Gauss-Newton method or gradient descent method.

[0063] This system designs a deterministic iterative solution engine based on the alternating direction multiplier method, which obtains a unique and exact solution through matrix operations on the central processing unit (CPU).

[0064] The system first introduces auxiliary variables. By substitution, the non-differentiable gradient term is removed from the reflectivity term. Peeling from the middle. Let. This transforms the original unconstrained extremum problem into an optimization problem with equality constraints. Auxiliary variable The two components.

[0065] Among them, auxiliary variables When solving ill-conditioned inverse problems with L1 total variational regularization using the Alternating Direction Multiplier Method (ADMM), gradient splitting variables, also called boundary feature auxiliary variables, are introduced to decouple variables and split non-differentiable terms.

[0066] Subsequently, the augmented Lagrangian function is constructed. : ; The augmented Lagrangian function is used to transform functions with equality constraints. The problem of minimizing the variational energy functional E(r,s) is transformed into an unconstrained subproblem that can be solved iteratively using the alternating direction multiplier method. Specifically, the augmented Lagrangian function introduces scaling dual variables into the original energy functional. and penalty parameters Constraining auxiliary variables through penalty terms With reflectivity gradient Maintaining consistency, thereby ensuring the logarithmic domain illumination field Logarithmic domain reflectance Auxiliary variables and scaled dual variables It can be updated separately in each iteration.

[0067] in, For the scaled dual variable, The penalty parameters are set. For gradient operators, For the logarithmic domain illumination field. For orthophotos Logarithmic domain observation image obtained after logarithmic mapping.

[0068] The core idea of ​​ADMM is to decompose complex multivariate joint optimization into multiple simple single-variable subproblems. The system iteratively solves the problem using the following four steps, assuming the current iteration number is... : fixed , , Update the illumination field : Extracting the augmented Lagrangian function that is only related to the logarithmic domain illumination field Related items, constructing the updated values ​​of iterative variables : ; Logarithmic domain illumination field Find the variational derivative and set it equal to zero to obtain the light field in the logarithmic domain. Euler-Lagrange differential equations: ; Summarized as follows: ; In the formula, This is a two-dimensional Laplace operator. Since this is a standard partial differential equation involving the Laplace operator, solving it using Jacobi iteration in the spatial domain results in an astonishingly high computational cost and extremely slow convergence.

[0069] The system calls the Fast Fourier Transform (FFT) within the loop to map the equation to the frequency domain for solution. Based on the derivative property theorem of the Fourier transform, the Laplace operator becomes a scalar multiplication in the frequency domain.

[0070] make This represents the two-dimensional discrete Fourier transform, whose analytical closed-form solution directly yields the updated illumination field of the iterative variables. : ; In the formula, For the discretized forward difference convolution kernel, 1 represents an all-one matrix, and this division is an element-wise division of the matrix elements. Through this purely mathematical technique, each... The update no longer requires thousands of iterations, but instead obtains an accurate analytical solution instantly through two FFT operations, meeting the needs of real-time processing of massive remote sensing data.

[0071] fixed , , Update log reflectance : Similarly, extract the logarithmic domain reflectance to be solved. Related terms, construct sub-problems : ; Find the Euler-Lagrange equation for it: ; The result is obtained by sorting and including the divergence operator. and Laplace operator The equation: ; The system also uses the frequency domain acceleration principle to give sub-problems. Analytical update equation: ; This step removes the illumination residual and injects the physical boundary sharpness information retained from the previous iteration into the reflectivity of this round through the divergence term.

[0072] fixed Equal variables, using the soft threshold operator to update auxiliary boundary variables. : Determine the form of the subproblem. : ; Based on the proximal gradient algorithm, the classic soft shrinkage threshold operator from one-dimensional signal processing is directly introduced, and a pixel-level closed solution is given at the matrix level. , : ; ; in .

[0073] If the reflectivity gradient at the current location Less than the set threshold (This means that this is a small undulation caused by sensor noise within the same type of ground feature), the operator directly clears it to zero, forcing it to become smooth; if the gradient is greater than the threshold (meaning that this is the real boundary of the ground feature), the operator preserves its gradient direction and shrinks the value, protecting the geometric topological boundary of the ground surface.

[0074] Update dual Lagrange multipliers : As the conclusion of the standard ADMM process, the scaled dual variable Perform gradient ascent updates to accumulate violations of equality constraints. historical error , : ; ; Convergence criteria and inverse logarithmic mapping output: At the end of each iteration, the system calculates the principal residual. and dual residual When both are less than the system's set minimum physical tolerance value. (For example When the maximum number of analyses is reached, the partial differential equation system is considered to have converged.

[0075] Finally, the obtained logarithmic domain result is subjected to an exponential inverse mapping to output the surface physical reflectivity matrix: ; ; At this point, the output of module four is... It is a pure surface physical reflectivity matrix that does not contain any mountain shadows, cloud cover, or slope light attenuation.

[0076] Module 5, Geoscience Threshold Classification Module, is configured as follows: calculate multispectral feature index based on the surface physical reflectance matrix, and output vector polygons representing spatial topological relationships after global threshold binarization and morphological optimization. After rigorous physical spatial alignment, forward evolution of illumination priors, and inverse solution of partial differential equations in the first four modules, the system successfully eliminated all non-geological / non-biological physical radiation interference caused by topographic relief, changes in solar altitude angle, and atmospheric scattering. This module receives the pure surface physical reflectivity matrix output by module four. Then, perform the final business-level target extraction.

[0077] Pixel-by-pixel reconstruction of multispectral feature indices: Since the partial differential systems of Modules 3 and 4 act independently on each spectral band of the remote sensing image, the surface physical reflectance matrix input to this module... It actually includes multiple bands (such as the red light band). Near-infrared band Shortwave infrared band Green light spectrum band The three-dimensional physical tensor of ).

[0078] Based on these calculated spectral reflectance matrices, the system recalculates classical remote sensing physical indices. Taking the extraction of forest cover and hidden water bodies in mountainous areas as an example, the system performs the following deterministic algebraic operations: Calculate the normalized vegetation index to characterize the true physical abundance of chlorophyll in vegetation: ; Calculate the normalized water index to accurately extract the water surface after eliminating hill shadows: ; In the raw images not processed by this system, deep valley shadow areas have extremely low radiant energy across all bands (close to the sensor noise floor), resulting in meaningless random values ​​for the calculated traditional NDVI, leading to large areas of voids in vegetation extraction. Simultaneously, the strong absorption characteristics of shadows in the shortwave infrared band cause abnormally high traditional MNDWI values, leading to a high probability of misidentification as water bodies by existing algorithms. However, in the output of this system… and In the matrix, since both the denominator and numerator restore the true reflectance of the matter, even in a valley that is physically completely backlit, pine trees will still exhibit a strong red-edge effect, and water will still exhibit standard infrared absorption characteristics.

[0079] Binarization segmentation based on a global rigid threshold: Because the physical properties of the same material are uniform in space, the system does not require the introduction of any complex adaptive thresholding algorithms or superpixel segmentation. The system directly calls a pre-defined set of rigid global thresholds with universal geochemical significance. ; in, A rigid global threshold is set for classifying healthy vegetation; A rigid global threshold is set for classifying pure water bodies; A rigid global threshold is used to classify bare soil / bedrock.

[0080] For example, performing a Boolean judgment: for the spatial domain any coordinate in , if (For example, setting) If the value is 0, then the pixel category is determined to be "healthy vegetation"; if and If so, the pixel category is determined to be "pure water". As a result, the system generates binary target masks for various ground features. .

[0081] Topology optimization and vector output based on mathematical morphology: Despite the surface physical reflectance matrix The image is already clean, but to eliminate salt-and-pepper noise from individual dead pixels in the sensor, the system applies mathematical morphological opening and closing operations to the binary mask: First, a morphological opening operation (erosion followed by dilation) is performed to eliminate isolated fragment noise points smaller than a set area threshold. Then, a closing operation (dilation followed by erosion) is performed to fill in the tiny holes inside the extracted features caused by physical mixing pixels. The morphological structuring element is forcibly set to a 3×3 isotropic cross-shaped core.

[0082] Finally, the system calls a topology tracing algorithm (such as Moore's neighbor tracing) in the underlying interface of the Geographic Information System (GIS) to rasterize the data. It is directly converted into a vector polygon (Shapefile or GeoJSON format) containing spatial topological relationships, and assigned a corresponding spatial reference system (such as WGS84) before being output to the terminal business system.

[0083] The foregoing has only described certain exemplary embodiments of the present invention by way of illustration. Undoubtedly, those skilled in the art can modify the described embodiments in various ways without departing from the spirit and scope of the present invention. Therefore, the foregoing drawings and descriptions are illustrative in nature and should not be construed as limiting the scope of protection of the claims of the present invention.

[0084] It should be noted that, in this document, the use of relational terms such as "first" and "second" is merely for distinguishing one entity or operation from another, and does not necessarily require or imply any such actual relationship or order between these entities or operations. Furthermore, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or apparatus. Without further limitations, an element defined by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, method, article, or apparatus that includes the element.

[0085] It should be understood that in the various embodiments of this application, the order of the above-mentioned processes does not imply the order of execution. The execution order of each process should be determined by its function and internal logic, and should not constitute any limitation on the implementation process of the embodiments of this application.

[0086] Those skilled in the art will recognize that the units and algorithm steps of the various examples described in conjunction with the embodiments disclosed herein can be implemented in electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are implemented in hardware or software depends on the specific application and design constraints of the technical solution. Those skilled in the art can use different methods to implement the described functions for each specific application, but such implementation should not be considered beyond the scope of this application.

[0087] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

Claims

1. A computer vision-based remote sensing image processing system, characterized in that, include: The multi-source data spatiotemporal registration and parameter analysis module is configured to: use rational polynomials to orthorectify and resample the original remote sensing images and digital elevation model data to generate orthorectified images; calculate the surface normal vector matrix and analyze the solar illumination direction vector at the imaging time; The surface illumination prior construction module is configured to: calculate the direct illumination component and the total diffuse reflection component based on the surface normal vector matrix and the solar illumination direction vector, and superimpose them to generate a geophysical radiation prior penalty term. The variational energy functional construction module is configured to: treat the orthophoto as the product of reflectivity and illumination field and map it to the logarithmic domain, and construct a variational energy functional that includes a data fidelity term, a reflectivity smoothing term based on the first norm, an illumination field continuity term based on the second norm, and the geophysical radiation prior penalty term. The partial differential equation iterative solution module is configured to: use the alternating direction multiplier method to iteratively solve the variational energy functional, update the illumination field and reflectivity in the frequency domain, update the boundary using the soft shrinkage threshold operator, and after convergence, perform inverse mapping to output the surface physical reflectivity matrix that has removed terrain illumination interference. The geoscience threshold classification module is configured to: calculate a multispectral feature index based on the surface physical reflectance matrix, and after global threshold binarization, output a vector polygon representing the spatial topological relationship through morphological optimization.

2. The computer vision-based remote sensing image processing system of claim 1, wherein, The process of generating orthorectified images by orthorectifying and resampling the original remote sensing images and digital elevation model data using rational polynomials includes: The input data comprises single or multi-band raw remote sensing image data and a digital elevation model covering the same geographical area wherein, is a remote sensing image column number, is a remote sensing image row number, is a longitude, is a latitude; Construct three-dimensional coordinates of normalized ground longitude, latitude, and elevation. Two-dimensional pixel coordinates of the normalized values ​​of the row number and column number of the original remote sensing image. The mapping relationship; With digital elevation model The target image grid is generated based on the geographic grid, and any geographic coordinate point in the target image grid is extracted. Corresponding elevation value ; Will After scale normalization, the geographic points in the original remote sensing image data are calculated. The corresponding two-dimensional pixel coordinates A bicubic convolution interpolation algorithm is used, utilizing two-dimensional cell coordinates. The grayscale value of the target point is calculated by resampling the surrounding integer pixel values, and an orthophoto with spatial consistency with the target image grid is generated based on the grayscale value. .

3. The computer vision-based remote sensing image processing system according to claim 2, characterized in that, The process of calculating the surface normal vector matrix and resolving the sunlight direction vector at the imaging time is as follows: Digital Elevation Model The Gaussian smoothing operator is applied, and then the elevation gradient at each spatial location is solved. , ; Calculate and After obtaining the partial derivatives of elevation in the direction, a digital elevation model is constructed based on the definition of the normal vector of a space surface in calculus. The unnormalized 3D normal vector at each grid point in the array ; By calculating the three-dimensional normal vector The modulus length is used to obtain the unit surface normal vector matrix at each pixel. ; From raw remote sensing image data The latitude and longitude of the image center point were extracted. and imaging Coordinated Universal Time ; Based on imaging Coordinated Universal Time The Julian Day is calculated, and then the declination angle of the sun at the moment of image formation is determined. Sum of time angles ; Combined with the latitude of the target area Solve for the solar altitude angle. and solar azimuth ; It is then converted into a solar illumination direction vector based on the geographic northeast-sky coordinate system. .

4. The computer vision-based remote sensing image processing system according to claim 3, characterized in that, Based on the aforementioned surface normal vector matrix and solar illumination direction vector, the direct sunlight component and the total diffuse reflection component are calculated as follows: By analyzing the normal vector matrix and the direction vector of sunlight Perform pixel-by-pixel dot product operations to solve for the initial direct cosine coefficients. ;when At that time, its direct sunlight illuminance was set to 0; For digital elevation models Each cell in The system emits a ray in three-dimensional space along the direction opposite to the sun's rays. , These are the step parameters; In each step, calculate the current ray. Projected coordinates on a two-dimensional plane And query the location in the digital elevation model. The corresponding actual surface elevation ; If there is any step point within the set maximum tracking distance satisfy Then the pixel It is physically obscured by the terrain and is in cast shadow; Generate a Boolean occlusion mask matrix 1 represents an unobstructed value, and 0 represents an obstructed value. Based on the direct cosine coefficient Boolean occlusion mask matrix Obtain the direct light component ,in, This is the solar constant at the top of the atmosphere. The transmittance parameter of the entire atmosphere at the time of imaging; In this pixel At that location, multiple probe rays are emitted into the three-dimensional hemispherical space at the azimuth angle; Define azimuth sampling interval ; along the first Each azimuth direction The indicated horizontal direction in the digital elevation model The search proceeds from top to bottom, calculating the elevation angles of all surface undulations relative to the center point, and recording the maximum elevation angle. Solve for the sky view factor. ; This is the total number of azimuth angle samples; In obtaining the sky view factor Then, calculate the ambient diffuse light component. , where constant This represents the diffuse light energy over flat ground; and by introducing an approximate cross-reflection term gain, we obtain the corrected total diffuse reflection component. ,in, This is an estimate of the average albedo of the background area. The average direct irradiance of the surrounding terrain.

5. The computer vision-based remote sensing image processing system according to claim 4, characterized in that, Direct light component With the total diffuse component By performing linear superposition, we obtain the geophysical radiation prior penalty term. .

6. The computer vision-based remote sensing image processing system according to claim 1, characterized in that, The variational energy functional is a global energy functional. , is represented as: ； in, For data fidelity items, For logarithmic domain reflectance, For logarithmic domain illumination field; This is a reflectance smoothing term based on the first norm; For the illumination field continuity term based on the 2-norm; These are prior constraints on the illumination field, used to constrain the logarithmic domain illumination field to be solved. Logarithmic field form approaching the geophysical radiation prior field , Geophysical radiation a priori penalty term The logarithmic field form; , , Preset the Lagrange hyperparameters; ; for Logarithmic domain reflectance The spatial gradient vector; Represents the logarithmic domain reflectance along Partial derivatives in direction; Represents the logarithmic domain reflectance along Partial derivatives in direction; ; for Logarithmic domain illumination field at the location The spatial gradient vector; Represents the logarithmic domain illumination field along Partial derivatives in direction; Represents the logarithmic domain illumination field along Partial derivatives in direction; The spatial domain representing the entire remote sensing image; , This represents a two-dimensional spatial area element.

7. The computer vision-based remote sensing image processing system according to claim 6, characterized in that, The partial differential equation iterative solution module also includes constructing an augmented Lagrangian function, used to transform the variational energy functional. The minimization problem is transformed into an unconstrained subproblem that can be solved iteratively using the alternating direction multiplier method; The variational energy functional There exists a non-differentiable norm term in it. ; Introducing auxiliary variables By substitution, the non-differentiable gradient term is removed from the reflectivity term. Separation from the middle; ; Auxiliary variable The two components; Therefore, we construct the augmented Lagrangian function. : in, For the scaled dual variable, The penalty parameters are set. For gradient operators, For orthophotos Logarithmic domain observation image obtained after logarithmic mapping.

8. The computer vision-based remote sensing image processing system according to claim 7, characterized in that, The iterative solution process for the augmented Lagrange function is as follows: Extracting the augmented Lagrangian function that is only related to the logarithmic domain illumination field Related items, constructing the updated values ​​of iterative variables Logarithmic domain illumination field Find the variational derivative and set it equal to zero to obtain the light field in the logarithmic domain. The Euler-Lagrange differential equation is solved by two-dimensional discrete Fourier transform to obtain the updated values ​​of the iterative variables. ; This is the current iteration number; Extracting the logarithmic domain reflectance to be solved Related terms, construct sub-problems ; Solve for its Euler-Lagrange equations; rearrange to obtain the equation containing the divergence operator. and Laplace operator The equations are given; subproblems are also given. The analytical update equation; Determine the form of the subproblem ; Based on the proximal gradient algorithm, a classic soft contraction threshold operator from one-dimensional signal processing is introduced, and a closed-form solution is given at the matrix level. , ; Scaled dual variable Perform gradient ascent updates to accumulate violations of equality constraints. historical error , ; At the end of each iteration, the principal residual is calculated. and dual residual When both are less than the set minimum physical tolerance value The partial differential equation system is considered convergent when the preset maximum number of analytical iterations is reached. Finally, regarding the logarithmic domain reflectance Perform an exponential inverse mapping to output the surface physical reflectance matrix. .

9. The computer vision-based remote sensing image processing system according to claim 1, characterized in that, The geoscience threshold classification module receives the surface physical reflectance matrix and calculates the multispectral feature index. The specific steps of global threshold binarization are as follows: Based on the surface physical reflectance matrix after removing topographic lighting interference, multispectral feature indices are reconstructed pixel by pixel. The multispectral feature indices include at least the normalized vegetation index and the normalized water index. The system calls a pre-defined set of global rigid thresholds, performs Boolean judgments on the multispectral feature indices, and generates binarized target masks for various land features.

10. The computer vision-based remote sensing image processing system according to claim 9, characterized in that, After generating the binarized target mask, the specific steps for morphological optimization to output the vector polygon are as follows: An isotropic cross-shaped kernel with a set size of 3×3 is used as the morphological structural element to perform mathematical morphological opening and closing operations on the binarized target mask. The topology tracing algorithm is invoked to convert the rasterized binary target mask, which has been processed by mathematical morphology opening and closing operations, into a vector polygon containing spatial topological relationships, and then output it with the corresponding spatial reference frame.

Images