Method for monitoring resident-free island by adopting satellite remote sensing

Abstract

The invention relates to the technical field of satellite remote sensing monitoring, in particular to a method for monitoring a resident-free island by adopting satellite remote sensing, which comprises the following steps of: acquiring optical images and cloud layer data of a target island through optical and meteorological satellites, predicting a cloud shielding probability, and scheduling a radar satellite when the cloud shielding probability exceeds a threshold value; comparing the effective area of the optical image with a reference demand, triggering radar imaging when a difference value reaches a preset proportion, extracting a coastline feature point set from a historical image, calculating a spatial distribution entropy screening subset, constructing an optimized affine transformation matrix to compensate a radar image, and obtaining a new optical image when the cloud cover is low. Coastline feature point errors of the calibration radar image and the new optical image are calculated, secondary compensation is carried out, global time sequence monitoring data are generated, monitoring continuity and data consistency in a cloud shielding scene are improved, and the method is suitable for resident-free island dynamic monitoring.

Description

Technical Field

[0001] The present invention relates to the technical field of satellite remote sensing monitoring, and in particular to a method for monitoring uninhabited islands using satellite remote sensing. Background Art

[0002] In the field of satellite remote sensing monitoring of uninhabited islands, cloud cover is the core problem that leads to the loss of optical image data. Existing technologies mainly use infrared and visible light band data from meteorological satellites to build cloud cover probability models, combine historical cloud motion trajectories to predict short-term cloud dynamics, and use this information to dispatch radar satellites as an alternative to optical imagery. However, this method lacks dynamic adaptability to complex meteorological systems such as monsoon vortices. Simple correlation models of cloud top temperature, cloud albedo, and sea surface humidity make it difficult to accurately quantify cloud fragmentation index and vertical development gradient, leading to large deviations in cloud obscuration probability predictions, which in turn can cause ineffective scheduling of radar satellites or missing data during critical periods. Although radar images have the advantage of penetrating clouds, they face technical bottlenecks in monitoring uninhabited islands. Tidal changes cause coastline feature points, such as reef apex and turning points at the base of sea cliffs, to drift dynamically. Traditional feature extraction algorithms do not distinguish between tidal-sensitive and tidal-insensitive features, resulting in inconsistent spatial distribution of feature point sets between historical optical images and radar images. This in turn leads to the polarization scattering characteristics of radar images, such as the scattering entropy and scattering angle of the VV / VH channel, which are essentially different from the texture characteristics of optical images. The existing calculation of spatial distribution entropy values ​​does not integrate the coastline curvature segmentation constraint mechanism, which makes the feature point screening results deviate from the actual space. distribution, which makes the affine transformation optimization process rely on the inverse distance weight matrix, but does not introduce the terrain shielding angle correction of the digital elevation model, resulting in the error of displacement field reconstruction being significantly amplified in areas with sudden terrain changes such as steep slopes and valleys, further deriving the problem of continuity break in time series data. The unoptimized radar image coordinate offset compensation will accumulate residuals, causing the error of the average Euclidean distance between the new optical image and the radar image coastline feature points in the secondary calibration stage to exceed the limit, ultimately affecting the accuracy of the global change monitoring of uninhabited islands. To solve this technical problem, we provide a method for monitoring uninhabited islands using satellite remote sensing. Summary of the Invention

[0003] The object of the present invention is to provide a method for monitoring uninhabited islands using satellite remote sensing to solve the problems raised in the above background technology.

[0004] Because existing technologies have large deviations in cloud obscuration probability predictions and lack adaptability to complex weather conditions, resulting in ineffective radar scheduling or data loss, this case constructs a three-dimensional cloud field model using multi-source meteorological data, combines it with monsoon vortex characteristic vector matching to calculate probability, and triggers radar scheduling, which can improve the accuracy of cloud obscuration predictions and reduce resource waste.

[0005] Because existing feature point extraction techniques do not distinguish between tidal sensitivity and do not incorporate terrain masking corrections, resulting in large image calibration errors, this case uses an anti-tidal interference algorithm to extract stable feature points, integrates terrain masking coefficients to construct a weight matrix to optimize the affine transformation, which can improve image matching accuracy and ensure the consistency of time series data.

[0006] To achieve the above objectives, one of the objectives of the present invention is to provide a method for monitoring uninhabited islands using satellite remote sensing, comprising the following steps: Build a three-dimensional cloud field dynamic model and output the cloud obscuration probability and the difference between the effective area of ​​the optical image and the reference required area. When the cloud obscuration probability and the difference reach the corresponding preset thresholds, the radar satellite is driven to acquire radar images. When the cloud cover probability and difference both exceed the corresponding preset thresholds: The historical optical image coastline feature point set and radar image feature point set are extracted from the historical optical image and radar image, respectively. The spatial distribution entropy values ​​of the historical optical image coastline feature point set and the radar image feature point set are calculated. A feature point subset is screened based on the spatial distribution entropy values. An inverse distance weight matrix is ​​constructed based on the feature point subset. The L1 norm is introduced to optimize the inverse distance weight matrix to obtain an affine transformation matrix. The affine transformation matrix is ​​applied to perform coordinate offset compensation on the radar satellite to obtain a calibrated radar image. When the cloud cover probability and difference are both lower than the corresponding preset thresholds: Reacquire the new optical image, calculate the average Euclidean distance error of the coastline feature points of the calibration radar image and the new optical image, perform secondary compensation calibration on the calibration radar image based on the average Euclidean distance error, and generate time series monitoring data for the entire uninhabited island area.

[0007] Compared with the prior art, the present invention has the following beneficial effects: An anti-tidal interference algorithm is used to extract tide-insensitive feature points such as reef apex. Stable point sets are screened through spatial density clustering and gradient direction test to avoid the drift of feature points caused by tides. For radar images, matching points are extracted through dual-polarization scattering feature enhancement. The spatial distribution entropy is calculated by combining the coastline curvature segmentation to ensure that the feature point screening fits the actual coastal morphology. At the same time, the terrain shielding coefficient of the digital elevation model is integrated to construct an inverse distance weight matrix. The affine transformation matrix is ​​solved by L1 norm optimization to effectively correct errors in areas with sudden terrain changes such as steep slopes and valleys, thereby improving the coastline alignment accuracy between radar images and optical images.

[0008] Gradient domain fusion compensation is used to achieve fine correction of radar image coordinate offset, and then secondary compensation calibration is performed in combination with new optical images. The average Euclidean distance error of coastline feature points is calculated through multi-level confidence control, and residual deviations are dynamically corrected. This solves the residual accumulation problem of traditional methods, ensures the consistency of long-term time series monitoring data, and can accurately reflect subtle changes in island topography and coastline, providing reliable data support for the ecological protection and rights protection of uninhabited islands. BRIEF DESCRIPTION OF THE DRAWINGS

[0009] Figure 1 A diagram showing the steps of the method of the present invention; Figure 2 It is the overall workflow diagram of the present invention. DETAILED DESCRIPTION

[0010] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of the present invention.

[0011] See also Figure 1-Figure 2 As shown, this embodiment provides a method for monitoring uninhabited islands using satellite remote sensing, including the following steps: Build a three-dimensional cloud field dynamic model and output the cloud obscuration probability and the difference between the effective area of ​​the optical image and the reference required area. When the cloud obscuration probability and the difference reach the corresponding preset thresholds, the radar satellite is driven to acquire radar images. Predict the probability of cloud obscuration, specifically: In the process of using satellite remote sensing to monitor uninhabited islands, accurate prediction of cloud cover probability is the key prerequisite for deciding whether to dispatch radar satellites. This link is achieved through a multi-source meteorological data fusion engine, which not only integrates multi-dimensional meteorological observation data, but also improves prediction accuracy by combining historical cloud movement laws. It forms a coherent technical chain with the subsequent image acquisition and calibration process, and simultaneously accesses the infrared band cloud top temperature of the meteorological satellite, the visible light band cloud albedo and the sea surface humidity of the ground meteorological station. Among them, the infrared band cloud top temperature is inverted by measuring the infrared radiation intensity at the top of the cloud layer. The lower the value, the warmer the temperature. The higher the clear cloud layer, the thicker it is. The cloud albedo in the visible light band refers to the proportion of visible light reflected by the cloud layer. The higher the albedo, the thicker the cloud layer. The sea surface humidity measured by ground meteorological stations reflects the water vapor content of the sea surface air. The higher the humidity, the easier it is to form or maintain clouds. Based on this data, a three-dimensional cloud field dynamic model is constructed. The cloud top temperature and albedo data are gridded according to spatial coordinates. The cloud formation potential of each grid is calculated in combination with the sea surface humidity. Then, numerical simulation is used to predict the spatial distribution and movement trend of clouds in the next 1 to 3 hours. The final output is real-time cloud field data including cloud position, thickness, and movement direction. The monsoon vortex characteristic vector of the target island sea area is extracted based on the historical cloud motion trajectory database. That is, from the cloud motion trajectory data of the past five years, the monsoon vortex events affecting the sea area are screened out, and their typical features are extracted to form a vector. Among them, the angular velocity of the cloud movement direction reflects the speed of cloud rotation or translation (unit: degree / hour, positive value indicates clockwise, negative value indicates counterclockwise), the cloud fragmentation index is the degree to which the cloud is divided into small pieces (with a value of 0-1, 0 indicates a complete cloud, and 1 indicates complete fragmentation), and the vertical development gradient reflects the vertical development of the cloud. The thickening rate in the direction of the cloud is calculated (in kilometers per hour, a positive value indicates upward development). The real-time cloud field data is dynamically time-warped and matched with the monsoon vortex characteristic vector. The probability value of cloud coverage of the target island and sea area is calculated. The real-time cloud field data is decomposed into a data point every 10 minutes according to the time series. Each point contains three features: cloud center coordinates, moving speed, and albedo mean. At the same time, the monsoon vortex characteristic vector is also interpolated at the same time interval to form a comparable feature sequence. The sea area weight factor is introduced to give higher weight to the cloud characteristics close to the target island. For example, the weight within 50 kilometers from the island is 1.2, and 50-100 kilometers is 0.8. The dynamic time warping algorithm is used to find the matching path that minimizes the difference between the two sets of sequence features. For example, the cloud movement speed in a certain period of time in the real-time cloud field is 5 kilometers per hour, which is matched with the period of 4-6 kilometers per hour in the historical vortex characteristics. The Euclidean distance between the two is calculated and accumulated, and the accumulated distance is converted into similarity. Then, the probability of cloud existence in the grid where the target island is located in the real-time cloud field is calculated by albedo and humidity, and the final result is obtained by fusion through the Bayesian formula. For example, if the real-time similarity is 0.8 and the probability of cloud existence in the grid is 0.7, the final obscuration probability is 0.8×0.7+(1-0.8)×0.3=0.62. When the probability value exceeds the preset probability threshold, a radar satellite dispatch instruction is generated. This process avoids the deviation of single data prediction through multi-source data fusion and matching with historical rules, and strengthens the targeting of target islands through sea area weight adjustment, providing an accurate basis for the subsequent timely dispatch of radar satellites in the event of cloud obscuration, ensuring the continuity of monitoring data.

[0012] When the cloud cover probability and difference both exceed the corresponding preset thresholds: The historical optical image coastline feature point set and radar image feature point set are extracted from the historical optical image and radar image, respectively. The spatial distribution entropy values ​​of the historical optical image coastline feature point set and the radar image feature point set are calculated. A feature point subset is screened based on the spatial distribution entropy values. An inverse distance weight matrix is ​​constructed based on the feature point subset. The L1 norm is introduced to optimize the inverse distance weight matrix to obtain an affine transformation matrix. The affine transformation matrix is ​​applied to perform coordinate offset compensation on the radar satellite to obtain a calibrated radar image. The stability feature extraction algorithm against tidal interference is adopted to extract the feature point set of the historical optical image coastline, specifically: After completing the cloud shielding probability prediction and determining whether to dispatch the radar satellite, to ensure the accuracy of subsequent image calibration, the stable coastline feature point set needs to be extracted from the historical optical image. This process adopts the stability feature extraction algorithm against tidal interference, filters the feature points that are not affected by tidal fluctuations, and provides a reliable reference for radar image and optical image matching. It forms a progressive link in data processing with the cloud field monitoring described above. When identifying the reef vertex, the base of the sea cliff turning point, and the permanent artificial feature corner on the historical optical image, the system first performs multi-scale edge detection on the image: For the reef vertex, the gray value mutation region of high-resolution optical image is identified. The reef is above the sea surface and has hard texture, so the junction line between its top and seawater will form a clear gray gradient. The algorithm determines the vertex coordinates by finding the maximum gradient point. For example, a certain reef appears as a triangular outline in the image, with a gray value of 50 units higher than the surrounding seawater at the top pixel, and the gradients around it all point to this point, which is determined as the reef vertex. For the base of the sea cliff turning point, focus on the junction between the cliff and the beach. The shadow formed by the vertical height of the cliff and the bright area of the beach form a strong contrast, and the turning point is the intersection of the shadow edge and the beach texture. For example, the base of a certain cliff has a pixel with a gray value of less than 50 on the cliff side and greater than 150 on the beach side, and this point has the maximum edge curvature. The permanent artificial feature corner is identified through shape matching, such as the right angle of a lighthouse base and the polygon corner of an observation platform. The algorithm locates the corner features of the pre-stored artificial features in the image through template matching and verifies their positional stability in multiple historical images. The coordinate drift in three consecutive images is less than 2 pixels. These identified points form the tidal-insensitive feature set. Since the corresponding entities are all above the high tide level, they will not be submerged or change shape due to tidal changes, ensuring the stability of the feature points. When performing spatial density clustering analysis on each type of feature point in the tidal-insensitive feature set, the coastal zone partition clustering mechanism is introduced: First, the coastline is divided into multiple sub-areas with a length of 500 meters. In each sub-area, density clustering is performed on similar feature points, such as reef vertices. A neighborhood with a radius of 30 meters is set with each feature point as the center. It is dynamically adjusted according to the complexity of the island terrain. The number of feature points in the neighborhood is counted. If the number is greater than or equal to 5, it is marked as a high-density clustering core. The feature points within 30 meters from the core point are then included in the same cluster to form multiple clustering units. For example, there are 12 reef vertices in a sub-area, and two clusters are formed after clustering, containing 7 and 5 points respectively. Feature points with spatial density lower than the preset threshold are then eliminated. The threshold is set according to the average point density in the cluster, such as the whole The average number of points in each cluster in the region is 6, so the threshold is set to 3. Clusters with less than 3 points are determined to be isolated points or noise points and are removed. For example, a cluster contains only 2 reef vertices. Because its density is insufficient to reflect the coastal characteristics, it is deleted. The main direction of the remaining feature points is screened by edge gradient direction consistency test. Edge gradient direction consistency test means calculating the gradient direction of edge pixels around the feature point to verify whether these directions are concentrated around a dominant direction, so as to determine whether the feature points are located on the same structural edge. The main direction of the feature point refers to the main extension direction of the edge where the feature point is located, which is used to ensure the direction consistency of the feature points during image matching. The specific screening process is as follows: For each retained feature point, the edge gradient direction within the 3×3 pixel window around it is extracted, and the distribution histogram of these directions is statistically analyzed. If the pixel proportion of a certain direction in the histogram exceeds 60%, that is, most of the gradient directions are consistent, then this direction is determined to be the main direction of the feature point. If the proportion is less than 60%, it indicates that the edge around the point is messy and it may be a noise point, so it is removed. For example, the gradient direction around a reef apex mostly points to the northwest-southeast direction, accounting for 70%, then its main direction is the northwest-southeast direction, and the point is retained. If the gradient direction around a point is scattered, with a maximum proportion of only 40%, it is removed, which lays a high-quality foundation for subsequent matching with the radar image feature point set, makes the calibration of cross-sensor images more accurate, and thus ensures the temporal continuity of the uninhabited island monitoring data.

[0013] The multi-polarization scattering feature enhancement method is used to extract the radar image feature point set, specifically: When extracting radar image feature point sets, a multi-polarization scattering feature enhancement method is used. This process complements the extraction of historical optical image coastline feature point sets in the previous article. By enhancing the difference in scattering characteristics between radar waves and island objects, it provides stable features for cross-image matching and is closely connected with the subsequent spatial distribution entropy calculation link. The VV polarization channel of the radar image refers to the polarization mode of vertical transmission and vertical reception, which has a strong response to the scattered signals of sea surface roughness and metal artificial objects. The VH polarization channel is a polarization mode of vertical transmission and horizontal reception, which is easier to capture the scattering characteristics of complex structures such as vegetation and reefs. The combination of the two can fully reflect the scattering characteristics of the objects. When performing sea clutter suppression filtering on the VV and VH polarization channels respectively, a polarization adaptive threshold filtering mechanism is adopted: First, the statistical characteristics of sea clutter are calculated for the VV channel, and a dynamic threshold is set. Signals with a variance greater than three times the mean are identified as targets, retaining strong scattering signals. For the VH channel, the polarization ratio is introduced, and the power ratio of VV to VH is used to assist in filtering. When the polarization ratio is greater than 2, it indicates that the ground object is more likely to be a rigid structure, and the filtering threshold is lowered to retain more weak scattering features. For example, in a certain area, the mean of sea clutter in the VV channel is 20, the variance is 5, and the dynamic threshold is set to 35. Pixels with signal strengths exceeding 35 are retained. At the same time, the polarization ratio in this area is calculated to be 3, and the threshold of the VH channel is lowered by 20% to ensure that weak scattering targets such as reefs are not filtered out. After filtering, the effective signals of the two channels are fused according to pixel position to form a dual-polarization image. This not only suppresses more than 90% of sea clutter but also retains the polarization scattering differences of the ground object. When locating strong scattering target points on the dual-polarization image, multi-window sliding detection is used to achieve the following: Two 3×3 and 5×5 windows are slid across the image, and the average scattering intensity within each window is calculated. When the average intensity of both windows is twice the average of the entire image, the center pixel of the window is marked as a strong scattering target point. For example, if the average intensity of a pixel in the 3×3 window is 80 (the average of the entire image is 30) and in the 5×5 window is 75, it is located as a strong scattering target point. When calculating the polarization scattering entropy value, it is necessary to first obtain the polarization scattering matrix of the point to describe the scattering power relationship of polarization channels such as VV and VH. The three eigenvalues ​​are obtained by matrix eigenvalue decomposition, and then the entropy value is calculated according to the formula. The value ranges from 0 to 1, where 0 represents a single scattering mechanism and 1 represents a mixture of multiple scattering mechanisms. For example, the apex of a reef is dominated by surface scattering, with an entropy value of about 0.3. The vegetation area has a high proportion of volume scattering, and the entropy value can reach 0.8. The scattering angle is the angle between the incident direction of the radar wave and the surface normal of the object, ranging from 0 to 90 degrees. It is calculated based on the orbital parameters of the radar satellite, such as the incident angle, and the altitude data of the target point. For example, the scattering angle of near-shore reefs is mostly between 30 and 45 degrees. When screening candidate feature points, the preset conditions are: The polarization scattering entropy value is ≤0.5, ensuring that surface scattering is dominant. This corresponds to rigid structures such as reefs and artificial features, and the scattering angle is between 20 and 60 degrees. Signal distortion caused by too steep or too slow an angle of incidence is eliminated. For example, a strong scattering target point with an entropy value of 0.4 and a scattering angle of 35 degrees meets the requirements and is selected as a candidate point. Points with an entropy value of 0.7 are eliminated due to their complex scattering mechanism. After projecting the candidate feature points to the geographic coordinate system, they are spatially matched with the historical optical image coastline feature point set, using a buffer level matching strategy: First, using the minimum bounding rectangle of the historical optical feature point set as a benchmark, a first-level buffer zone is formed by expanding outward 500 meters to cover islands and nearshore areas. Candidate points falling within this zone are initially retained. The feature types within the first-level buffer zone are then subdivided, and a second-level buffer zone with a radius of 200 meters is established for each optical feature point. The spatial overlap between the candidate point and optical feature points of the same type is calculated. For example, if a candidate point falls within the second-level buffer zone of three reef-type optical feature points, the overlap is 3, and candidate points with an overlap of 1 or more are retained. For example, reef vertices in historical optical imagery are concentrated near 121°30′E and 30°20′N, with a second-level buffer zone of 200 meters. Radar candidate points within this zone whose scattering characteristics match reef properties are retained, ultimately forming the radar image feature point set. This provides a high-quality matching foundation for the subsequent calculation of spatial distribution entropy and the construction of the affine transformation matrix, making the calibration of radar and optical images more reliable.

[0014] The coastline morphological constraint factor is introduced to calculate the spatial distribution entropy value, specifically: After obtaining the historical optical image coastline feature point set and the radar image feature point set, in order to accurately screen the feature point subset that can reflect the true shape of the coastline, it is necessary to calculate the spatial distribution entropy of the two. The introduction of the coastline morphological constraint factor can make the entropy calculation more suitable for the morphological characteristics of different coastal segments. This process inherits the results of the previous feature point extraction and lays a quantitative foundation for subsequent feature point matching. When segmenting the coastline into straight segments and curved segments according to curvature changes, the sliding window curvature detection method is used: First, connect the coastline feature points into a broken line in spatial order, and use every 5 consecutive points as a sliding window to calculate the average curvature of the broken line in the window. The larger the curvature value, the more curved the line segment. Set a curvature threshold (such as 0.01 radians / meter). When the average curvature in the window is less than the threshold, the segment is judged to be a straight line segment. If the average curvature of 3 consecutive windows is greater than the threshold, it is marked as a curved segment, and the starting point of the first window and the end point of the last window are used as the boundaries of the curved segment. For example, the average curvature of the 5-point window of a certain coastline is 0.008, 0.007, and 0.009, respectively, which are all less than the threshold. , is divided into a straight line segment, while the window curvature of the other segment is 0.012, 0.015, and 0.013, all exceeding the threshold, and is divided into a curved segment. Through this dynamic segmentation method, the transition of the coastline from flat to curved can be accurately captured. The Euclidean space grid division method is used to calculate the Shannon entropy value in the straight line segment. The Euclidean space grid division method refers to dividing the two-dimensional space where the straight line segment is located into equidistant grids, so that each grid corresponds to a fixed area in the actual geographic space. The Shannon entropy value is used to measure the degree of discreteness of the distribution of feature points in the grid. The larger the value, the more dispersed the distribution, and the smaller the value, the more concentrated the distribution. The specific calculation process is as follows: Several 10m x 10m grids are drawn along the straight line segment. The number of feature points in each grid is counted. According to the Shannon entropy formula, the ratio of the number of points in each grid to the total number of points is used as the probability value. The negative value of the sum of the product of the probability values ​​and the logarithm of the probability is calculated as the Shannon entropy value. For example, a straight line segment has 50 feature points distributed in 10 grids, of which 3 grids contain 8 points each, 5 grids contain 6 points each, and 2 grids contain 2 points each. The calculated Shannon entropy value is 1.8, reflecting that the distribution of feature points in this segment is relatively uniform. The morphological deviation entropy value is calculated for the curve segment using the Fréchet distance matching method. The Fréchet distance matching method is a method for measuring the similarity of two curves. By calculating the minimum and maximum distances between the points on the two curves, the degree of similarity of the curve morphology is judged. The morphological deviation entropy value quantifies the morphological difference between the feature point sequence and the standard curve template based on the Fréchet distance. The larger the value, the more significant the deviation, and the smaller the value, the closer the morphology. The specific calculation process is as follows: First, a standard curve template is selected based on the curvature characteristics of the curve segment. The sequence of feature points is connected to form the actual curve. The Fréchet distance between the actual curve and the standard template is calculated, and the ratio of this distance to the curve segment length is used as the deviation probability. Then, combined with the calculation logic of Shannon entropy, the deviation probabilities of multiple segments are converted into entropy values ​​to obtain the morphological deviation entropy value. For example, a curve segment is divided into three segments, and the ratios of the Fréchet distance to the length are 0.1, 0.15, and 0.08, respectively. The calculated morphological deviation entropy value is 0.6, indicating that the morphology of this curve segment deviates significantly from the standard template. Finally, the weighted sum of the Shannon entropy value and the morphological deviation entropy value is used as the spatial distribution entropy value. The Shannon entropy value of the straight segment has a weight of 0.6, and the morphological deviation entropy value of the curved segment has a weight of 0.4, because the morphological characteristics of the curved segment are more critical to coastline matching. For example, the Shannon entropy value of a straight segment is 1.8, and the morphological deviation entropy value of a curved segment is 0.6. The comprehensive spatial distribution entropy value is 1.8×0.6+0.6×0.4=1.32. This entropy calculation method combined with coastline morphology can more accurately reflect the spatial distribution characteristics of the feature point set, provide a scientific basis for the subsequent screening of feature point subsets, and make the feature matching of optical images and radar images more consistent with the actual coastal morphology.

[0015] The feature point subset is screened through bidirectional stability verification, specifically: After completing the calculation of the spatial distribution entropy of the feature point sets of the historical optical image and the radar image, it is necessary to screen the feature point subset through bidirectional stability verification. This process is the core link to ensure the accuracy of the subsequent affine transformation matrix. It not only inherits the entropy quantification results in the previous article, but also strengthens the stability of the feature points through dual verification in the geometric and temporal dimensions, providing reliable control points for image calibration. Point pairs refer to a combination of feature points with similar spatial positions selected from the historical optical image coastline feature point set and the radar image feature point set. Each point pair contains one optical feature point and one radar feature point. The two correspond to the same feature in geographic space, such as the same reef vertex, and are the basic units for cross-sensor image matching. After calculating the spatial distribution entropy difference between the feature point sets of the historical optical image and the radar image, point pairs with a difference within the threshold range (such as a difference ≤ 0.3) are retained. The spatial distribution characteristics of these point pairs are closer and have a matching basis. The point pairs are then verified for geometric invariance using the local grid residual analysis method: A 50m x 50m local grid is defined with each point pair as the center. The coordinates of the optical feature points are used as the reference. A preliminary affine transformation, including rotation and translation parameters, is performed to calculate the predicted coordinates of the radar feature points in the optical coordinate system. The pixel difference between the two is the local affine transformation residual. 1 pixel corresponds to 1 meter of actual distance. The preset residual threshold is 3 pixels, or 3 meters. If the residual of a point pair is 4 meters, it exceeds the threshold and is eliminated. If the residual is 2 meters, the point pair is retained. This step ensures the consistency of the point pairs in spatial transformation by quantifying the geometric deviation. When performing temporal invariance verification on the remaining point pairs, a three-phase image drift tracking mechanism is introduced: The coordinate data for the point pair was extracted from three historical imagery sessions over the past year, separated by four months. The coordinate drift between each two-period period was calculated. For example, the drift between the first and second periods was 0.5 meters, and between the second and third periods was 0.3 meters. The standard deviation of these drifts was then calculated, with a threshold of 0.4 meters. If a point pair's drift standard deviation was 0.2 meters, indicating stable position, it passed the validation. If the standard deviation was 0.5 meters, indicating significant position fluctuation, the point pair was deemed unstable and discarded. For example, a point pair at the top of a reef exhibited minimal coordinate change across the three imagery sessions, with a standard deviation of 0.15 meters, significantly below the threshold, and was retained. However, a point pair in a nearshore sandy area exhibited a drift standard deviation of 0.6 meters due to scouring and was discarded. Ultimately, the point pairs that passed both geometric and temporal invariance verification constituted the feature point subset. For example, 100 pairs of points were initially screened, 70 pairs were retained after entropy difference screening, 50 pairs remained after geometric residual verification, and 35 stable pairs were ultimately retained after temporal drift verification, forming a subset of feature points. This layer-by-layer screening process ensures both the spatial matching accuracy of the point pairs and their stability in the temporal dimension, providing high-quality control points for the subsequent construction of the affine transformation matrix. This ensures that the calibration of radar and optical images better matches the actual terrain changes on the island, thereby improving the consistency of global time-series monitoring data.

[0016] The terrain shielding coefficient is integrated to construct the distance inverse weight matrix, specifically: After screening out the feature point subset, in order to construct a distance inverse weight matrix that better fits the island terrain, it is necessary to integrate the terrain shielding coefficient for weighted correction. This process not only inherits the results of the previous feature point stability verification, but also improves the rationality of weight distribution by introducing terrain factors, providing an accurate weight basis for the subsequent optimization solution of the affine transformation matrix. The digital elevation model data is the target island terrain elevation data obtained through satellite remote sensing or aerial measurement. The altitude of each point is recorded in a grid form, such as a grid point every 10 meters, and the vertical distance of the point above sea level is recorded. It can intuitively reflect the terrain characteristics of the island, such as the undulating terrain and the direction of the mountain. When calculating the terrain shielding angle of the connecting line direction of each feature point subset, first extract the two end points of the point in the digital elevation model. The altitude of the optical feature points and radar feature points and the altitude of the midpoint of the connecting line. For example, if point A is 50 meters above sea level, point B is 60 meters above sea level, and the midpoint of the connecting line is 45 meters above sea level, then the midpoint of the connecting line is used as the vertex to calculate the angle formed by the connecting line between the two end points and the midpoint. Specifically, the slope angles from point A to the midpoint and from point B to the midpoint are calculated using the inverse tangent function based on the altitude difference and horizontal distance between the two points. The difference between the two slope angles is the terrain shielding angle. For example, if the slope angle from point A to the midpoint is 10° and the slope angle from point B to the midpoint is 8°, then the terrain shielding angle is 2°. This angle reflects the degree of shielding of the connecting line by the terrain. The larger the angle, the more significant the shielding. The terrain shielding coefficient is obtained by comprehensively considering the altitude difference, horizontal distance, and terrain shielding angle in the connecting line direction: First calculate the altitude difference between the two points. For example, if the altitude difference between point A and point B is 10 meters and the horizontal distance, calculate it through geographic coordinates, such as 100 meters, to obtain the relative altitude difference ratio (10 / 100=0.1). Then, combined with the terrain shielding angle (such as 2°), calculate it through the formula terrain shielding coefficient = 0.5×relative altitude difference ratio + 0.5×(terrain shielding angle / 90°). For example, in the above example, the coefficient is 0.5×0.1+0.5×(2° / 90°)≈0.05+0.0 11=0.061. The larger the coefficient value, the stronger the shielding effect of the terrain on the line connecting the point pair. When constructing the inverse distance weight matrix, the inverse of the spatial distance between the point pairs in the feature point subset is used as the initial weight. For example, if the horizontal distance between the two points is 100 meters, the initial weight is 1 / 100=0.01. The initial weight is then multiplied by the inverse of the terrain shielding coefficient (1 / 0.061≈16.39) for weighted correction. The corrected weight is 0.01×16.39≈0.164. The same calculation is performed on all point pairs in the feature point subsets, and the corrected weights are arranged in order of the point pairs to form an inverse distance weight matrix. For example, a feature point subset contains 3 pairs of point pairs, and the corrected weights are 0.164, 0.21, and 0.18, respectively. The matrix is ​​constructed with these three values ​​as elements to ensure that point pairs with large terrain masking effects receive higher weights in subsequent matrix solutions, making the affine transformation matrix more in line with the image matching requirements under actual terrain conditions, providing a more accurate weight basis for the subsequent optimization solution of the affine transformation matrix, and making the compensation calibration of radar images more in line with the complex terrain characteristics of the island, thereby improving the accuracy of cross-sensor image matching.

[0017] When introducing the L1 norm to optimize the distance inverse weight matrix, an adaptive smoothing constraint mechanism is adopted, specifically: After constructing the distance-inverse weight matrix that incorporates the terrain masking coefficient, the L1 norm needs to be introduced to optimize and solve it to obtain the optimal affine transformation matrix. This process balances the geometric consistency and spatial distribution characteristics of the feature points through an adaptive smoothing constraint mechanism. It not only inherits the quantization results of the weight matrix in the previous article, but also provides a mathematical basis for the precise compensation of radar images, forming a close connection with the subsequent image calibration process. The L1 norm refers to the sum of the absolute values ​​of each element in the quantity. In this method, it is used to measure the overall deviation of the weighted displacement residual. Compared with the L2 norm, it is less sensitive to outliers and can effectively suppress the interference of noise points on the solution. The affine transformation matrix describes the geometric transformation relationship between radar and optical images. It is derived from the coordinate correspondence between points in a feature point subset and uses mathematical modeling to achieve spatial alignment between the two types of images. The matrix contains three parameters: rotation, translation, and scaling. The rotation parameter reflects the image's rotation angle around the origin (in degrees), the translation parameter reflects the image's offset in the x and y directions (in meters), and the scaling parameter represents the image's horizontal and vertical scaling ratio (unitless, 1 indicates uniform scaling). When initializing these parameters, a coarse alignment + empirical value correction strategy is used: First, the initial translation parameter is calculated by the average coordinate difference of the point pairs in the feature point subset. For example, if the average x-coordinate of the optical point is 50 meters larger than that of the radar point, the x-direction translation parameter is initialized to -50. Based on the overall direction of the island coastline, such as the east-west direction, the rotation parameter is initialized to 0 degrees. The scaling parameter is based on the resolution ratio of the optical image and the radar image, such as 1 meter / pixel for optical images and 2 meters / pixel for radar images. The scaling parameter is initialized to 0.5. The inverse distance weight matrix is ​​used as the basis for residual weighting because the matrix has integrated terrain masking and spatial distance factors and can reflect the matching reliability of different point pairs. Point pairs with higher weights, such as those with less terrain masking and closer distances, have higher matching accuracy and should be given greater weight in the residual calculation. The residual weighting basis refers to using the elements in the weight matrix to weight the displacement residual of each point pair so that the residual of the reliable point pair accounts for a higher proportion in the objective function. When constructing the objective function, a dynamic residual threshold is introduced: For each point pair, the difference between the predicted coordinates of the radar feature point after affine transformation and the actual coordinates of the optical feature point is calculated, that is, the displacement residual. This residual is multiplied by the weight of the corresponding point pair in the inverse distance weight matrix, and the sum of the absolute values ​​is taken as the objective function value. By minimizing this value, the overall residual is minimized under the constraint of reliable point pairs. For example, a point pair with a weight of 0.2 has a residual of 3 meters and a weighted residual of 0.6. Another point pair with a weight of 0.1 has a residual of 5 meters and a weighted residual of 0.5. The objective function will prioritize reducing the residual of high-weighted key pairs. Since the spatial distribution density of the feature point subset may be uneven, the smoothing constraint strength needs to be dynamically adjusted according to the density. This step is achieved through density-constraint mapping: The smoothing constraint strength refers to the coefficient for controlling the smoothness of the affine transformation parameter change. The greater the value, the smoother the parameter change, avoiding over-fitting local point pairs. First, calculate the point pair density in each 100m x 100m grid. For example, if a grid contains 8 point pairs, the density is 8. Calculate the average density of the entire area (e.g., 5). When the local density is higher than 30% of the average value (e.g., ≥6.5), reduce the smoothing constraint strength (e.g., from 1.0 to 0.7) to allow greater parameter changes in that area to fit the details of dense point pairs. When the local density is lower than 30% of the average value (e.g., ≤3.5), increase the smoothing constraint strength (e.g., from 1.0 to 1.3) to avoid parameter fluctuations caused by sparse point pairs. This adjustment ensures the fitting accuracy of dense areas and prevents over-distortion in sparse areas, naturally linking the construction and solution of the objective function. Finally, solve the optimal affine transformation matrix using the iterative reweighted least squares method: First, calculate the weighted residual of each point pair with the initial parameters, giving smaller iteration weights to high-weight residual points. For example, if the residual is greater than twice the average, the weight of the point pair is reduced to 0.5. Second, recalculate the affine transformation parameters based on the adjusted weights and update the objective function value. Third, repeat the first two steps, checking the change in the objective function value after each iteration. Stop iterating when the change is less than 0.001. Fourth, output the optimal affine transformation matrix composed of the final parameters. For example, after 5 iterations, the objective function value decreases from 120 to 15, with a change of 0.0008. The matrix parameters at this time are the optimal solution, laying the foundation for the consistency of global time series monitoring data.

[0018] The gradient domain fusion compensation is used to perform coordinate offset compensation for radar satellites, including the following steps: After obtaining the optimal affine transformation matrix, the radar image needs to be calibrated through gradient domain fusion compensation. This process is a key link between mathematical transformation models and actual image coordinate correction, connecting the solution of the affine transformation matrix and improving the matching accuracy of radar images and optical images through global and local compensation, laying the foundation for subsequent secondary compensation and calibration. When decomposing the affine transformation matrix into global affine components and local elastic components, the matrix hierarchical analysis method is used: The global affine component corresponds to the linear transformation part of the matrix, including rotation, translation, and scaling parameters. These parameters describe the overall geometric deviation of the radar image relative to the optical image. The local elastic component is the non-linear deviation part remaining after removing the global affine component from the affine transformation matrix, reflecting the local subtle deviations caused by terrain undulations and imaging distortion. During decomposition, the rotation angle, translation amount, and scaling factor in the matrix are extracted as global affine components. Then, the local elastic component is obtained through matrix subtraction to achieve the separation of the two. The global affine component is resampled using bilinear interpolation to obtain the global displacement field. The process is as follows: First, the radar image is divided into a 10×10 pixel grid. The target coordinates of each grid vertex in the optical image coordinate system are calculated based on the parameters of the global affine component. For example, a vertex should move from (x1, y1) to (x2, y2) after rotation and translation. Then, within each grid, the displacement of all pixels in the grid is calculated using a bilinear interpolation algorithm. That is, based on the target coordinates of the four vertices of the grid, the corresponding displacement values ​​are weighted according to the relative positions of the pixels in the grid. Finally, a global displacement field covering the entire image is formed. This displacement field can correct the overall geometric deviation of the radar image and make the image roughly aligned with the optical image. The gradient field of the radar image is extracted for the local elastic component and Poisson reconstruction is performed to generate the local displacement field. The steps are as follows: First, the gradient values ​​of the radar image in the x and y directions are calculated using the Sobel operator to reflect the rate and direction of pixel grayscale changes, forming a gradient field. Then, the gradient field is subjected to Poisson reconstruction using the deviation matrix of the local elastic component as a constraint. That is, by solving the Poisson equation, the gradient field is converted into a continuous displacement distribution, so that the reconstructed displacement field can meet the gradient constraint while conforming to the nonlinear deviation of the local elastic component. For example, in the reef edge area, the gradient field shows drastic changes. Poisson reconstruction will generate corresponding local stretching displacement to correct the edge distortion of the radar image caused by the imaging angle. After superimposing the global displacement field with the local displacement field, coordinate offset compensation is performed on the radar image. This process lies in the dynamic weight superposition + multi-scale compensation mechanism: First, dynamic weights are assigned to the global and local displacement fields based on the terrain complexity of the imaged area. In flat areas (such as beaches), the global displacement field weight is set to 0.8, and the local displacement field weight is 0.2 to ensure overall alignment accuracy. In areas with complex terrain (such as reefs), the local displacement field weight is increased to 0.7, and the global displacement field weight is 0.3 to enhance local detail correction. After superposition, a comprehensive displacement field is generated. Coordinate offset compensation is then performed on each pixel of the radar image. The target coordinates of the pixel are calculated based on the comprehensive displacement field. Then, a sub-pixel offset and edge sharpening strategy is applied. Sub-pixel interpolation is used to preserve detail in the offset pixels. Key edge areas, such as coastlines, are also sharpened to avoid blurry edges after compensation. For example, a reef pixel needs to be offset by 2.3 pixels after comprehensive displacement field calculation. It is accurately located at the target position through sub-pixel interpolation, and its edge contrast with the sea water is enhanced. The final calibrated radar image can be highly matched with the optical image both overall and locally. The alignment error of the coastline feature points can be controlled within 2 pixels, providing high-quality basic data for subsequent secondary compensation calibration with the new optical image, effectively improving the data consistency of cross-sensor images in uninhabited island monitoring.

[0019] When the cloud cover probability and difference are both lower than the corresponding preset thresholds: Reacquire new optical image, and calculate the average Euclidean distance error of the coastline feature points of the calibrated radar image and the new optical image. According to the average Euclidean distance error, the calibrated radar image is secondarily compensated and calibrated to generate the global time-series monitoring data of the uninhabited islands.

[0020] The average Euclidean distance error of the coastline feature points of the calculated calibrated radar image and the new optical image is realized through multi-level confidence control, specifically: After completing the gradient domain fusion compensation of the radar image, in order to verify the matching accuracy of the calibrated radar image and the new optical image, the average Euclidean distance error of the coastline feature points of the two images needs to be calculated. This process is realized through multi-level confidence control, which not only connects the image calibration results in the previous text, but also provides quantitative basis for possible secondary compensation. It forms a closed-loop verification logic with the generation of global time-series monitoring data. The coastline edge pixel chain refers to the continuous pixel sequence along the coastline contour in the new optical image and the calibrated radar image. These pixels have significant differences in gray scale or scattering characteristics from the surrounding sea area. For example, the gray value of the coastline pixel in the optical image suddenly changes, and the radar image shows strong scattering signals. Together, they form a continuous pixel set reflecting the shape of the coastline. When the coastline edge pixel chain is sampled to generate a verification point set at equal arc length intervals, the complete pixel chain is first obtained through edge tracking algorithm. For example, a certain section of the coastline pixel chain contains 1000 continuous pixels. The total arc length of the pixel chain is calculated, and the pixels are converted to actual distance based on geographic coordinates. For example, the total arc length is 500 meters. The sampling interval is determined at equal arc length intervals (such as 10 meters). Starting from the starting point of the pixel chain, every 10 meters selects a pixel as a verification point. If the end point is less than 10 meters, it is combined into the last verification point. Finally, 50 verification points are generated to form a verification point set, ensuring that the sampling points uniformly cover the entire coastline and avoiding local dense or sparse errors. After taking the coordinates of each verification point on the new optical image as the reference, a pixel search window is established at the corresponding position of the radar image: First, through geographic coordinate conversion, the latitude and longitude coordinates of the optical image verification point are mapped to the pixel coordinate system of the radar image to obtain the initial corresponding position. Considering the possible residual deviation, a square search window of 11×11 pixels is established with this position as the center, that is, 5 pixels are extended to each side. The window size can be dynamically adjusted according to the accuracy of the previous calibration. For example, when the calibration error is large, it can be expanded to 15×15 pixels to ensure that the window contains the true matching point. When calculating the point with the highest spectral similarity with the reference point in the search window as the matching point, the optical reference point is extracted. The grayscale feature vectors of the 3×3 neighborhood, such as mean, variance, and gradient direction, are extracted for each pixel in the radar image search window. The corresponding dual-polarization scattering feature vector (such as VV / VH polarization ratio and scattering entropy value) is extracted. The matching degree of the two sets of feature vectors is calculated using the cosine similarity algorithm (the value range is 0-1, 1 indicates a complete match). The pixel with the highest similarity is selected as the matching point, and the Euclidean distance is calculated based on the geographic coordinates of the two. For example, the distance between the optical reference point coordinates (x1, y1) and the radar matching point coordinates (x2, y2) is √[(x1-x2) 2 +(y1-y2) 2 ], record the distance value as the error of the verification point, perform box plot analysis on the Euclidean distance values ​​of all verification points to eliminate outliers, first sort all distance values ​​from small to large, calculate the quartiles, Q1 is the 25th percentile, Q3 is the 75th percentile, determine the interquartile range IQR=Q3-Q1, set the outlier judgment threshold to Q1-1.5×IQR and Q3+1.5×IQR, all verification points with distance values ​​outside this range are judged as outliers, such as a verification point distance value of 15 meters, which is much higher than Q3+1.5×IQR of 8 meters, and is judged as an anomaly. After eliminating these outliers, the remaining distance values ​​are considered to be credible matching errors Finally, the remaining Euclidean distance values ​​are weighted averaged according to the coast type, and the coastline is divided into reef coast, beach coast, artificial coast and other types. The weights are assigned according to the proportion of each type in the total coastline, and the average distance value within each type is calculated. Then, the weighted sum is calculated according to the weight to generate the final average Euclidean distance error. For example, the average error of the reef coast is 2 meters, and that of the beach coast is 3 meters. The total error after weighting is 2×0.4+3×0.6=2.6 meters. The final average Euclidean distance error can accurately reflect the coastline matching accuracy of the two types of images, providing a scientific basis for judging whether secondary compensation is needed, and further ensuring the consistency of the global time series monitoring data.

[0021] The present invention obtains optical images and cloud data of target islands through optical and meteorological satellites, predicts the probability of cloud cover, dispatches radar satellites when the threshold is exceeded, compares the effective area of ​​the optical image with the benchmark requirement, triggers radar imaging when the difference reaches a preset ratio, extracts a set of coastline feature points from historical images, calculates the spatial distribution entropy to filter the subset, constructs an optimized affine transformation matrix to compensate the radar image, obtains a new optical image when the cloud cover is low, calculates the error of the coastline feature points between the calibration radar image and the new optical image, performs secondary compensation, generates global time-series monitoring data, improves monitoring continuity and data consistency in cloud cover scenarios, and is suitable for dynamic monitoring of uninhabited islands.

[0022] The above shows and describes the basic principles, main features, and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The above embodiments and descriptions are merely preferred examples of the present invention and are not intended to limit the present invention. Various changes and improvements may be made to the present invention without departing from the spirit and scope of the present invention. Such changes and improvements fall within the scope of the present invention. The scope of protection claimed in the present invention is defined by the appended claims and their equivalents.

Claims

1. A method for monitoring uninhabited islands using satellite remote sensing, characterized by: The following steps are involved: Build a three-dimensional cloud field dynamic model and output the cloud obscuration probability and the difference between the effective area of ​​the optical image and the reference required area. When the cloud obscuration probability and the difference reach the corresponding preset thresholds, the radar satellite is driven to acquire radar images. When the cloud cover probability and difference both exceed the corresponding preset thresholds: The historical optical image coastline feature point set and radar image feature point set are extracted from the historical optical image and radar image, respectively. The spatial distribution entropy values ​​of the historical optical image coastline feature point set and the radar image feature point set are calculated. A feature point subset is screened based on the spatial distribution entropy values. An inverse distance weight matrix is ​​constructed based on the feature point subset. The L1 norm is introduced to optimize the inverse distance weight matrix to obtain an affine transformation matrix. The affine transformation matrix is ​​applied to perform coordinate offset compensation on the radar satellite to obtain a calibrated radar image. When the cloud cover probability and difference are both lower than the corresponding preset thresholds: Reacquire the new optical image, calculate the average Euclidean distance error of the coastline feature points of the calibration radar image and the new optical image, perform secondary compensation calibration on the calibration radar image based on the average Euclidean distance error, and generate time series monitoring data for the entire uninhabited island area.

2. The method for monitoring uninhabited islands using satellite remote sensing according to claim 1, characterized in that: Predicting cloud cover probability includes the following steps: Synchronously access infrared cloud top temperature from meteorological satellites, visible light cloud albedo, and sea surface humidity from ground meteorological stations to construct a three-dimensional cloud field dynamic model and obtain real-time cloud field data; The monsoon vortex characteristic vector of the target island and sea area is extracted based on the historical cloud motion trajectory library. The real-time cloud field data is dynamically time-warped matched with the monsoon vortex characteristic vector to calculate the probability value of cloud coverage of the target island and sea area. When the probability value exceeds the preset probability threshold, a radar satellite dispatch instruction is generated, in which the monsoon vortex characteristic vector includes the angular velocity of the cloud movement direction, the cloud fragmentation index and the vertical development gradient.

3. The method for monitoring uninhabited islands using satellite remote sensing according to claim 1, characterized in that: The coastline feature point set of historical optical images is extracted using a tidal interference-resistant stability feature extraction algorithm, which specifically includes the following steps: Identify reef apex, sea cliff base turning points, and permanent artificial feature corners on historical optical images to form a tide-insensitive feature set; A spatial density clustering analysis is performed on each type of feature points in the tide-insensitive feature set. Feature points with spatial density lower than the preset threshold are eliminated, and the main directions of the remaining feature points are screened through the edge gradient direction consistency test to finally generate the historical optical image coastline feature point set.

4. The method for monitoring uninhabited islands using satellite remote sensing according to claim 3, characterized in that: The multi-polarization scattering feature enhancement method is used to extract the radar image feature point set, which specifically includes the following steps: Perform sea clutter suppression filtering on the VV polarization and VH polarization channels of the radar image to obtain a dual-polarization image of the radar image. On the dual-polarization image, locate the strong scattering target points and calculate the polarization scattering entropy and scattering angle of each strong scattering target point. Target points that meet the preset conditions of polarization scattering entropy value and scattering angle are selected as candidate feature points. The candidate feature points are projected into the geographic coordinate system and spatially matched with the historical optical image coastline feature point set. The candidate feature points within the spatial range matching area are retained to form a radar image feature point set.

5. The method for monitoring uninhabited islands using satellite remote sensing according to claim 4, characterized in that: The coastline morphological constraint factor is introduced to calculate the spatial distribution entropy value, which includes the following steps: The coastline where the historical optical image coastline feature point set and the radar image feature point set are located is divided into straight line segments and curved line segments according to the curvature change; The Euclidean space grid division method is used in the straight segment to calculate the Shannon entropy value of the feature points in each Euclidean space grid. The Fréchet distance matching method is used in the curve segment to calculate the morphological deviation entropy value between the feature point sequence and the standard curve template. Finally, the weighted sum of the Shannon entropy value and the morphological deviation entropy value is used as the spatial distribution entropy value.

6. The method for monitoring uninhabited islands using satellite remote sensing according to claim 5, characterized in that: The feature point subset is screened through bidirectional stability verification, which includes the following steps: Calculate the spatial distribution entropy difference between the historical optical image coastline feature point set and the radar image feature point set, retain the point pairs that meet the difference threshold condition, and perform geometric invariance verification on the point pairs, namely: Each point pair is used as a control point to calculate the local affine transformation residual. The point pairs with residuals greater than the preset pixels are eliminated, and time invariance verification is performed on the remaining point pairs. The coordinate drift of the point pair in the three historical images is compared, and the verification point pairs are judged based on the standard deviation of the coordinate drift. Finally, the verified point pairs constitute the feature point subset.

7. The method for monitoring uninhabited islands using satellite remote sensing according to claim 6, characterized in that: The terrain shielding coefficient is integrated to construct an inverse distance weight matrix, which includes the following steps: The digital elevation model data of the target island is obtained, and the terrain shielding angle in the direction of the connecting line of each feature point subset is calculated. The inverse of the spatial distance of the point pairs in the feature point subset is used as the initial weight, and multiplied by the inverse of the terrain shielding coefficient to perform weighted correction to obtain the distance inverse weight matrix, where the terrain shielding coefficient is determined by the altitude difference in the connecting line direction, the horizontal distance, and the terrain shielding angle.

8. The method for monitoring uninhabited islands using satellite remote sensing according to claim 7, characterized in that: When introducing the L1 norm to optimize the distance inverse weight matrix, an adaptive smoothing constraint mechanism is adopted, which specifically includes the following steps: Initialize the rotation, translation, and scaling parameters of the affine transformation matrix, use the inverse distance weight matrix as the residual weighting basis, and construct the objective function to minimize the sum of the absolute values ​​of the weighted displacement residuals; The smoothing constraint strength is dynamically adjusted according to the spatial distribution density of the feature point subset, and the optimal affine transformation matrix is ​​solved by iterative reweighted least squares method.

9. The method for monitoring uninhabited islands using satellite remote sensing according to claim 8, characterized in that: Gradient domain fusion compensation is used to perform coordinate offset compensation on radar satellites, which specifically includes the following steps: The affine transformation matrix is ​​decomposed into global affine components and local elastic components. Bilinear interpolation is used to resample the image grid for the global affine components to obtain the global displacement field. The gradient field of the radar image is extracted for the local elastic components. Poisson reconstruction is performed on the gradient field to generate the local displacement field. The global displacement field and the local displacement field are superimposed, and coordinate offset compensation is performed on the radar satellite to generate a calibrated radar image.

10. The method for monitoring uninhabited islands using satellite remote sensing according to claim 9, characterized in that: The average Euclidean distance error of the coastline feature points between the calibrated radar image and the new optical image is calculated by multi-level confidence control. The following steps are involved: Coastline edge pixel chains are extracted synchronously from the new optical image and the calibrated radar image. Verification point sets are generated by sampling along the coastline edge pixel chains at equal arc length intervals. The coordinates of each verification point set on the new optical image are used as a benchmark, and a pixel search window is established at the corresponding position of the radar image. The point with the highest spectral similarity with the benchmark point within the pixel search window is calculated as the matching point, and the Euclidean distance between the two is recorded. A box plot analysis is performed on the Euclidean distance values ​​of all verification points, outliers are eliminated, and the remaining Euclidean distance values ​​are weighted averaged by coast type partition to generate the average Euclidean distance error.

Images