Method for monitoring uninhabited islands by satellite remote sensing

Abstract

The present application relates to the technical field of satellite remote sensing monitoring, in particular, the present application relates to a kind of method for monitoring uninhabited islands by satellite remote sensing, the present application obtains target island optical image and cloud data by optical and meteorological satellite, predicts cloud cover probability, dispatches radar satellite when threshold value, compares the effective area of optical image with benchmark demand, difference reaches preset proportion when triggering radar imaging, extract the coastline feature point set from historical image, calculate spatial distribution entropy to filter subset, build optimized affine transformation matrix to compensate radar image, obtain new optical image when cloud cover is low, calculate the coastline feature point error of calibration radar image and new optical image, carry out secondary compensation, generate global time series monitoring data, improve the monitoring continuity and data consistency under cloud cover scene, applicable to the dynamic monitoring of uninhabited islands.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of satellite remote sensing monitoring, in particular to a method for monitoring uninhabited islands by satellite remote sensing. BACKGROUND

[0002] In the field of satellite remote sensing monitoring of uninhabited islands, cloud cover is the core problem leading to the absence of optical image data. The existing technology mainly constructs a cloud coverage probability model through infrared and visible light band data of meteorological satellites, combines historical cloud movement trajectories to predict short-term cloud dynamics, and accordingly dispatches radar satellites as an alternative to optical images.

[0003] However, this method lacks dynamic adaptability to complex meteorological systems such as monsoon vortices, and the simple correlation model of cloud top temperature, cloud albedo, and sea surface humidity cannot accurately quantify the cloud fragmentation index and vertical development gradient, resulting in large prediction deviation of cloud cover probability, and further causing ineffective scheduling of radar satellites or data absence in critical periods.

[0004] Although radar images have the advantage of penetrating clouds, they face technical bottlenecks in monitoring uninhabited islands. Tidal changes cause dynamic shifts in coastal feature points such as reef apexes and sea cliff base turning points. Traditional feature extraction algorithms do not distinguish between tide-sensitive and tide-insensitive features, leading to inconsistent spatial distribution of feature point sets between historical optical images and radar images. This results in essential differences between the polarization scattering characteristics of radar images, such as the scattering entropy and scattering angle of VV / VH channels, and the texture features of optical images. Existing spatial distribution entropy calculations do not incorporate a coastline curvature segmentation constraint mechanism, causing feature point selection results to deviate from actual spatial distribution. This further causes the affine transformation optimization process to rely on a distance inversely proportional weight matrix, but does not introduce a digital elevation model terrain shadow angle correction, resulting in significant error amplification in steep slope and valley terrain mutation areas during displacement field reconstruction, further deriving a time series data continuity disruption problem. Unoptimized radar image coordinate offset compensation accumulates residual errors, causing the average Euclidean distance error between new optical images and radar images of coastal feature points in the secondary calibration stage to exceed the limit, ultimately affecting the accuracy of uninhabited island-wide change monitoring. To solve this technical problem, we provide a method for monitoring uninhabited islands by satellite remote sensing. SUMMARY

[0005] The present application aims to provide a method for monitoring uninhabited islands by satellite remote sensing to solve the problems raised in the background.

[0006] Since the prior art cloud obscuration probability prediction has large deviation, and is not suitable for complex weather, leading to ineffective radar scheduling or data loss, therefore, the case can improve the cloud obscuration prediction accuracy and reduce resource waste by constructing a three-dimensional cloud field model through multi-source weather data, and calculating the probability combined with the monsoon vortex feature vector matching to trigger radar scheduling.

[0007] Since the prior art feature point extraction does not distinguish tidal sensitivity and does not fuse terrain shielding correction, resulting in large image calibration error, therefore, the case can improve the image matching precision and ensure the coherence of time series data by extracting stable feature points through the anti-tidal interference algorithm, and fusing the terrain shielding coefficient to construct a weight matrix to optimize the affine transformation.

[0008] To achieve the above purpose, one of the purposes of the present application is to provide a method for monitoring uninhabited islands by satellite remote sensing, comprising the following steps:

[0009] A three-dimensional cloud field dynamic model is constructed, and the cloud obscuration probability and the difference between the effective area of the optical image and the reference required area are outputted, when the cloud obscuration probability and the difference respectively reach the corresponding preset threshold, the radar satellite is driven to obtain radar images;

[0010] When the cloud obscuration probability and the difference respectively exceed the corresponding preset threshold:

[0011] The historical optical image coastline feature point set and the radar image feature point set are extracted from the historical optical image and the radar image respectively, the spatial distribution entropy value of the historical optical image coastline feature point set and the radar image feature point set is calculated, the feature point subset is selected according to the spatial distribution entropy value, and the distance inverse weight matrix is constructed based on the feature point subset, the L1 norm is introduced to optimize and solve the distance inverse weight matrix, the affine transformation matrix is obtained, the coordinate offset compensation is performed on the radar satellite by applying the affine transformation matrix, and the calibrated radar image is obtained;

[0012] When the cloud obscuration probability and the difference respectively are lower than the corresponding preset threshold:

[0013] The new optical image is reacquired, the average Euclidean distance error of the coastline feature points of the calibrated radar image and the new optical image is calculated, the calibrated radar image is secondarily compensated and calibrated according to the average Euclidean distance error, and the time series monitoring data of the entire uninhabited island is generated.

[0014] Compared with the prior art, the present application has the following advantages:

[0015] The anti-tidal interference algorithm is used to extract the tidal insensitive feature points such as the top of the reef, and the stable point set is screened through the spatial density clustering and the gradient direction test, so that the feature point drift caused by the tide is avoided, the matching points are extracted through the dual polarization scattering feature enhancement of the radar image, and the spatial distribution entropy is calculated by combining the coastline curvature segmentation, so that the feature point screening is ensured to be in line with the actual coastline shape, meanwhile, the terrain shielding coefficient of the digital elevation model is fused to construct a distance inverse weight matrix, and the affine transformation matrix is solved through L1 norm optimization, so that the error of the steep slope, valley and other terrain mutation areas is effectively corrected, and the alignment accuracy of the coastline of the radar image and the optical image is improved.

[0016] The gradient domain fusion compensation is used to realize the fine correction of the coordinate offset of the radar image, and then the new optical image is used for secondary compensation and calibration, the average Euclidean distance error of the coastline feature points is calculated through multi-level confidence control, the residual deviation is dynamically corrected, the residual accumulation problem of the traditional method is solved, the long-term time sequence monitoring data is consistent, and the subtle changes of the island terrain and the coastline can be accurately reflected, so that reliable data support is provided for the ecological protection and right protection of the uninhabited island. BRIEF DESCRIPTION OF DRAWINGS

[0017] Fig. 1 The method steps of the present application are shown in the figure;

[0018] Fig. 2 The overall workflow of the present application is shown in the figure. DETAILED DESCRIPTION

[0019] The technical solutions in the embodiments of the present application will be clearly and completely described below with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are only part of the embodiments of the present application, not all. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor are within the scope of protection of the present application.

[0020] Please refer to Figs. 1-2 The present embodiment provides a method for monitoring uninhabited islands by satellite remote sensing, which comprises the following steps:

[0021] A three-dimensional cloud field dynamic model is constructed, and the cloud shielding probability and the difference between the effective area of the optical image and the reference required area are output. When the cloud shielding probability and the difference value respectively reach the corresponding preset threshold value, the radar satellite is driven to obtain the radar image;

[0022] The cloud shielding probability is predicted, specifically:

[0023] In the process of monitoring uninhabited islands by satellite remote sensing, the accurate prediction of cloud cover probability is the key prerequisite for determining whether to dispatch radar satellites. This link is realized through a multi-source meteorological data fusion engine, which integrates multi-dimensional meteorological observation data and combines historical cloud movement rules to improve prediction accuracy. It forms a coherent technical chain with the subsequent image acquisition and calibration process, and synchronously accesses the infrared band cloud top temperature of meteorological satellites, the visible band cloud albedo, and the sea surface humidity of ground meteorological stations. The infrared band cloud top temperature is inversely calculated from the infrared radiation intensity at the top of the cloud layer, with lower values indicating higher cloud height and greater thickness. The visible band cloud albedo refers to the proportion of visible light reflected by the cloud layer, with higher albedo indicating thicker and more substantial cloud layers. The sea surface humidity of ground meteorological stations reflects the water vapor content of the sea surface air, with higher humidity being more conducive to the formation and maintenance of cloud layers. Based on these data, a three-dimensional cloud field dynamic model is constructed, which gridded the cloud top temperature and albedo data according to spatial coordinates and calculates the cloud formation potential of each grid in combination with sea surface humidity. Then, through numerical simulation, the spatial distribution and movement trend of the cloud layer in the next 1-3 hours are predicted, and finally, real-time cloud field data including cloud position, thickness, and movement direction are outputted.

[0024] The monsoon vortex feature vector of the target island sea area is extracted based on the historical cloud movement trajectory library, that is, from the cloud movement trajectory data in the past 5 years, the monsoon vortex events affecting the sea area are screened out, and the typical features are extracted to form a vector, wherein the cloud movement direction angular velocity reflects the rotation or translation speed of the cloud cluster (unit: degree / hour, positive value indicates clockwise, negative value indicates counterclockwise), the cloud cluster fragmentation index is the degree of fragmentation of the cloud cluster (value 0-1, 0 indicates a complete cloud cluster, and 1 indicates complete fragmentation), and the vertical development gradient reflects the thickening speed of the cloud layer in the vertical direction (unit: kilometer / hour, positive value indicates upward development). The real-time cloud field data is dynamically time warping matched with the monsoon vortex feature vector, and the cloud layer covering the target island sea area is calculated. The probability value is divided into every 10 minutes according to the time sequence, and each point contains cloud center coordinates, moving speed, and albedo average value. The monsoon vortex feature vector is also interpolated according to the same time interval to form a comparable feature sequence. The sea area weight factor is introduced, and the cloud cluster features close to the target island are given higher weight, such as 1.2 within 50 kilometers of the island, 0.8 within 50-100 kilometers. The dynamic time warping algorithm is used to find the matching path with the smallest feature difference between the two sequences, for example, the real-time cloud field cloud layer moving speed is 5 kilometers / hour, which is matched with the 4-6 kilometers / hour period in the historical vortex feature. The Euclidean distance of the two is calculated and accumulated, and the accumulated distance is converted into similarity. Combined with the cloud layer existence probability of the target island grid in the real-time cloud field, the cloud layer existence probability is calculated by the albedo and humidity, and the final cloud shielding probability value is obtained by Bayes formula fusion. For example, the real-time similarity is 0.8, and the grid cloud layer existence probability is 0.7, so the final shielding probability is 0.8*0.7+(1-0.8)*0.3=0.62. When the probability value exceeds the preset probability threshold, the radar satellite scheduling instruction is generated. This process combines multi-source data fusion and historical rule matching, avoids the deviation of single data prediction, and strengthens the pertinence to the target island through sea area weight adjustment, providing accurate basis for timely scheduling of radar satellites in cloud shielding, and ensuring the continuity of monitoring data.

[0025] When the cloud shielding probability and the difference value both exceed the corresponding preset threshold:

[0026] The historical optical image coastline feature point set and the radar image feature point set are extracted from the historical optical image and the radar image respectively, the spatial distribution entropy value of the historical optical image coastline feature point set and the radar image feature point set is calculated, the feature point subset is selected according to the spatial distribution entropy value, and the distance inverse weight matrix is constructed based on the feature point subset. The affine transformation matrix is obtained by introducing L1 norm to optimize and solve the distance inverse weight matrix. The affine transformation matrix is applied to the radar satellite to perform coordinate offset compensation, and the calibrated radar image is obtained.

[0027] The stability feature extraction algorithm against tidal interference is adopted to extract the feature point set of the historical optical image coastline, specifically:

[0028] After completing the cloud shielding probability prediction and determining whether to dispatch the radar satellite, in order to ensure the accuracy of subsequent image calibration, it is necessary to extract a stable coastline feature point set from historical optical images. This process uses the stability feature extraction algorithm against tidal interference to filter out feature points that are not affected by tidal fluctuations, providing a reliable reference for radar image and optical image matching. This 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:

[0029] 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 cliff's vertical height 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. Permanent artificial feature corners are 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 artificial objects 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 a 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:

[0030] The coastline is first divided into multiple sub-regions with a length of 500 meters. In each sub-region, the same type of feature points, such as reef points, are subjected to density clustering. A neighborhood with a radius of 30 meters is set around each feature point, and the number of feature points in the neighborhood is counted based on the complexity of the island terrain. If the number is greater than or equal to 5, it is marked as a high-density clustering core. Then, feature points within a 30-meter range from the core point are included in the same cluster to form multiple clustering units. For example, in a certain sub-region, there are 12 reef points, which are clustered into 2 clusters containing 7 and 5 points respectively. Subsequently, feature points with a low spatial density below a preset threshold are removed. The threshold is set based on the average point density in the clustering cluster, such as an average of 6 points per cluster in the entire region, and the threshold is set to 3. For clustering clusters with fewer than 3 points, they are determined as isolated points or noise points and are removed. The main direction of the remaining feature points is selected through edge gradient direction consistency test, which is a test to verify whether the gradient directions of the edge pixels around the feature point are concentrated around a dominant direction to determine whether the feature point is located on the same structure edge. The main direction of the feature point is the main extension direction of the edge where the feature point is located, which is used to ensure the consistency of the direction of the feature points in image matching. The specific selection process is as follows:

[0031] For each retained feature point, the edge gradient directions in the 3x3 pixel window around it are extracted, and the distribution histogram of these directions is counted. If the proportion of pixels in a certain direction in the histogram is more than 60%, i.e. most of the gradient directions are consistent, then the direction is determined as the main direction of the feature point. If the proportion is less than 60%, it indicates that the edge around the point is chaotic and may be a noise point, which is removed. For example, the gradient directions around a certain reef point are mostly directed northwest-southeast, with a proportion of 70%, so the main direction is northwest-southeast, and the point is retained. If the gradient directions around a certain point are scattered, with a maximum proportion of only 40%, then it is removed. This lays a high-quality foundation for subsequent matching with the radar image feature point set, making the calibration of cross-sensor images more accurate, and thus ensuring the continuity of the non-resident island monitoring data in time series.

[0032] The radar image feature point set is extracted using a multi-polarization scattering feature enhancement method, which is as follows:

[0033] In the extraction of radar image feature point set, the multi-polarization scattering feature enhancement method is adopted, which is complementary to the extraction of historical optical image coastline feature point set in the previous text. By strengthening the scattering characteristic difference between radar wave and island ground object, stable features are provided for cross-image matching, and the subsequent spatial distribution entropy calculation link is closely linked. The VV polarization channel of radar image is perpendicular to the emission and reception polarization mode, which has strong response to the scattering signal of sea roughness and metal artificial ground object. The VV polarization channel is perpendicular to the emission and reception polarization mode, which is more easy to capture the scattering characteristics of vegetation, reef and other complex structures. The combination of the two can fully reflect the scattering characteristics of the ground object. When the VV and VH polarization channels are filtered respectively, the polarization adaptive threshold filtering mechanism is adopted:

[0034] First, calculate the statistical characteristics of sea clutter for VV channel, set dynamic threshold, signals higher than 3 times the mean value of variance are judged as target, and strong scattering signals are reserved. For VV channel, the polarization ratio is introduced, and the power ratio of VV and VH is used to assist filtering. When the polarization ratio is greater than 2, it indicates that the ground object is more likely to be rigid structure, and the filtering threshold is reduced to retain more weak scattering features. For example, the mean value of VV channel sea clutter in a certain area is 20, the variance is 5, the dynamic threshold is set to 35, the pixels with signal intensity greater than 35 are reserved, and the polarization ratio of the area is calculated to be 3. The threshold of VV channel is reduced by 20% to ensure that the weak scattering targets such as reefs are not filtered out. After filtering, the effective signals of the two channels are fused according to the pixel position to form a dual polarization image, which not only suppresses more than 90% of the sea clutter, but also retains the polarization scattering difference of the ground object. In the dual polarization image, strong scattering target points are located by multi-window sliding detection:

[0035] 3x3, 5x5 window is used to slide through the image, the average scattering intensity in each window is calculated, when the average intensity of the two windows is higher than 2 times the mean value of the whole image, the center pixel of the window is marked as a strong scattering target point. For example, a certain pixel has an average intensity of 80 (the mean value of the whole image is 30) in a 3x3 window and 75 in a 5x5 window, which is located as a strong scattering target point. When calculating the polarization scattering entropy value, the polarization scattering matrix of the point needs to be obtained first to describe the scattering power relationship of VV, VH and other polarization channels. Through eigenvalue decomposition of the matrix, 3 eigenvalues are obtained, and then the entropy value is calculated according to the formula, the value is 0-1, 0 represents single scattering mechanism, and 1 represents multiple scattering mixture. For example, the top of the reef is mainly surface scattering, the entropy value is about 0.3, and the vegetation area has high body scattering proportion, the entropy value can reach 0.8. Scattering angle is the angle between radar incident direction and ground surface normal, 0-90 degrees, which is calculated by radar satellite orbit parameters, such as incident angle and target point elevation data. The scattering angle of nearshore reef is mostly between 30-45 degrees. When screening candidate feature points, the preset conditions are:

[0036] Polarization scattering entropy value ≤ 0.5, to ensure that the surface scattering is given priority to, corresponding to the rigid structure of the reef, artificial features, and scattering angle between 20-60 degrees, excluding the incident angle too steep or too slow caused by signal distortion, for example, a strong scattering target point entropy value 0.4, scattering angle 35 degrees, meet the conditions to be selected as candidate points, entropy value 0.7 points because of the complex scattering mechanism is excluded, the candidate feature points are projected into the geographic coordinate system, and the historical optical image coastline feature point set is matched in space range, and a buffer zone level matching strategy is adopted:

[0037] First, the minimum circumscribed rectangle of the historical optical feature point set is taken as the benchmark, and a 500-meter buffer zone is formed, covering the islands and nearshore areas, and the points in the candidate points that fall within the area are preliminarily retained. Then, the land use types in the first buffer zone are subdivided, and a 200-meter radius second buffer zone is established for each optical feature point. The spatial coincidence degree of the candidate points with the same type of optical feature points is calculated. For example, if a candidate point falls within the second buffer zone of three reef-type optical feature points, the coincidence degree is 3, and the candidate points with a coincidence degree greater than or equal to 1 are retained. For example, the reef vertexes in the historical optical image are concentrated near east longitude 121°30', north latitude 30°20', and their second buffer zone is 200 meters. The radar candidate points that meet the scattering characteristics of the reef attribute within the range are retained, and the final radar image feature point set is constructed, providing a high-quality matching basis for subsequent calculation of spatial distribution entropy and construction of affine transformation matrix, making the calibration of radar image and optical image more reliable.

[0038] The coastline shape constraint factor is introduced to calculate the spatial distribution entropy value, which is:

[0039] After obtaining the historical optical image coastline feature point set and the radar image feature point set, in order to accurately select the feature point subset that reflects the true shape of the coastline, the spatial distribution entropy values of the two sets need to be calculated. The introduction of the coastline shape constraint factor can make the entropy value calculation more consistent with the shape characteristics of different coastlines. This process builds on the results of the previous feature point extraction and lays a quantitative foundation for subsequent feature point matching. When dividing the coastline into straight segments and curved segments according to the curvature change, the sliding window curvature detection method is adopted:

[0040] First, the coastline feature points are connected in spatial order into a polyline, and the average curvature of the polyline in each 5 consecutive points window is calculated. The greater the curvature value, the more curved the line segment. A curvature threshold (such as 0.01 rad / m) is set. When the average curvature in the window is less than the threshold, the segment is determined to be a straight line segment. If the average curvature of the last three windows is greater than the threshold, it is marked as a curve segment, and the starting point of the first window and the ending point of the last window are taken as the boundary of the curve segment. For example, the 5-point window average curvature of a certain segment of the coastline is 0.008, 0.007, and 0.009, all of which are less than the threshold, and it is divided into a straight line segment. Another segment has a window curvature of 0.012, 0.015, and 0.013, all of which exceed the threshold, and it is divided into a curve segment. Through this dynamic segmentation method, the transition from flat to curved of the coastline 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 of the straight line segment into equal grids, so that each grid corresponds to a fixed area in the actual geographical space. The Shannon entropy value is used to measure the dispersion degree of the feature points in the grid. The greater the value, the more dispersed the distribution, and the smaller the value, the more concentrated the distribution. The specific calculation process is as follows:

[0041] A number of 10m x 10m grids are divided along the extension direction of the straight line segment, and the number of feature points in each grid is counted. According to the Shannon entropy formula, the proportion of the number of points in each grid to the total number of points is taken as the probability value, and the sum of the product of all grid probability values and probability logarithm is calculated. The negative value of the sum is the Shannon entropy value. For example, a straight line segment has a total of 50 feature points, distributed in 10 grids, of which 3 grids each contain 8 points, 5 grids each contain 6 points, and 2 grids each contain 2 points. The calculated Shannon entropy value is 1.8, reflecting that the feature point distribution of this segment is relatively uniform. In the curve segment, the Fréchet distance matching method is used to calculate the shape deviation entropy value. The Fréchet distance matching method is a method for measuring the similarity of two curves. By calculating the minimum maximum distance between points on two curves, the closeness of the curve shape is determined. The shape deviation entropy value is based on the Fréchet distance to quantify the shape difference between the feature point sequence and the standard curve template. The greater the value, the more significant the deviation, and the smaller the value, the more consistent the shape. The specific calculation process is as follows:

[0042] According to the curvature characteristics of the curve segment, a standard curve template is selected, the characteristic point sequence is connected into an actual curve, the Fréchet distance between the actual curve and the standard template is calculated, the ratio of the distance to the length of the curve segment is taken as the deviation probability, and then the deviation probability of multiple segments is converted into an entropy value by combining the calculation logic of Shannon entropy to obtain a morphological deviation entropy value. For example, a curve segment is divided into 3 segments, and the ratio of the Fréchet distance to the length is 0.1, 0.15, and 0.08, respectively. The morphological deviation entropy value is calculated as 0.6, indicating that the morphological deviation of the curve segment from the standard template is large. Finally, the weighted sum of the Shannon entropy value and the morphological deviation entropy value is taken as the spatial distribution entropy value, in which the weight of the Shannon entropy value of the straight line segment is 0.6, and the weight of the morphological deviation entropy value of the curve segment is 0.4. Because the morphological characteristics of the curve segment are more critical for coastline matching, for example, the Shannon entropy value of a straight line segment is 1.8, and the morphological deviation entropy value of a curve segment is 0.6, the comprehensive spatial distribution entropy value is 1.8*0.6+0.6*0.4=1.32. Through this entropy value calculation method combined with the morphology of the coastline, the spatial distribution characteristics of the feature point set can be more accurately reflected, providing a scientific basis for subsequent screening of the feature point subset, and making the feature matching of the optical image and the radar image more consistent with the actual coastline morphology.

[0043] The feature point subset is screened through bidirectional stability verification, specifically as follows:

[0044] After the spatial distribution entropy value of the feature point set of the historical optical image and the radar image is calculated, the feature point subset is screened through bidirectional stability verification. This process is a key link to ensure the accuracy of the subsequent affine transformation matrix, which not only connects the entropy value quantization results of the previous text, but also strengthens the stability of the feature points through the dual verification of the geometric and time dimensions, providing reliable control points for image calibration. The point pair is a combination of feature points with similar spatial positions selected from the feature point set of the historical optical image coastline and the feature point set of the radar image. Each point pair contains one optical feature point and one radar feature point, which correspond to the same ground object in geographical space, such as the same reef vertex, and is the basic unit for cross-sensor image matching. After calculating the spatial distribution entropy difference value of the feature point set of the historical optical image and the radar image, the point pairs with a difference value within a threshold range (such as a difference value ≤0.3) are retained. These point pairs have closer spatial distribution characteristics and have the basis for matching. Then, the point pairs are verified for geometric invariance by using a local grid residual analysis method:

[0045] A local grid of 50m x 50m is drawn around each point pair, and the optical feature point coordinates are taken as the reference. Through a preliminary affine transformation, including rotation and translation parameters, the predicted coordinates of the radar feature points in the optical coordinate system are calculated. The pixel difference between the two is the local affine transformation residual. One pixel corresponds to an actual distance of 1 meter. The preset residual threshold is 3 pixels, i.e., 3 meters. If the residual of a certain point pair is 4 meters, it is excluded if it exceeds the threshold. If the residual is 2 meters, the point pair is retained. This step quantifies the geometric deviation to ensure the consistency of the point pairs in the spatial transformation. When performing time invariance verification on the remaining point pairs, a three-phase image drift tracking mechanism is introduced:

[0046] The coordinates of the point pair in the past three historical images, spaced 4 months apart, are extracted. The coordinate drift between each two periods is calculated, such as 0.5 meters from the first to the second period and 0.3 meters from the second to the third period. The standard deviation of the three-period drift is then calculated. The standard deviation threshold is set to 0.4 meters. If the drift standard deviation of a certain point pair is 0.2 meters, it indicates that the position is stable, and it passes the verification. If the standard deviation is 0.5 meters, the position fluctuates greatly, and it is determined as an unstable point pair and is excluded. For example, the point pair of a reef vertex has minimal coordinate change in the three images, with a standard deviation of 0.15 meters, which is significantly lower than the threshold and is retained. The point pair of a nearshore sandy area is excluded due to erosion, resulting in a drift standard deviation of 0.6 meters. Finally, the point pairs that pass both geometric invariance and time invariance verification form a feature point subset. For example, 100 point pairs are initially selected, 70 are retained after entropy difference screening, 50 remain after geometric residual verification, and 35 stable point pairs are finally retained after time drift verification, forming a feature point subset. This process ensures the spatial matching accuracy of the point pairs and their stability in the time dimension through layer-by-layer screening. It provides high-quality control points for subsequent affine transformation matrix construction, making the calibration of radar images and optical images more consistent with the actual changes in island topography, and thus improving the consistency of global time series monitoring data.

[0047] A distance inverse weight matrix is constructed by fusing the terrain shielding coefficient, specifically:

[0048] After screening out the feature point subset, in order to construct the distance inverse weight matrix more in line with the island terrain, the terrain shielding coefficient needs to be fused for weighted correction. This process not only connects the achievements of the feature point stability verification in the previous text, but also improves the rationality of weight distribution by introducing the terrain factor, providing accurate weight basis for the subsequent optimization solution of affine transformation matrix. The digital elevation model data is the target island terrain elevation data obtained by satellite remote sensing or aerial survey. It records the elevation of each point in the form of a grid, such as every 10 meters of a grid point, recording the vertical distance of the point above sea level. It can intuitively reflect the terrain features of the island, such as the ups and downs of the terrain and the direction of the mountains. When calculating the terrain shielding angle of each feature point subset connecting line direction, first extract the two end points of the point pair in the digital elevation model, the elevations of the optical feature point and the radar feature point and the elevation of the midpoint of the connecting line, for example, the elevation of point A is 50 meters, the elevation of point B is 60 meters, and the elevation of the midpoint of the connecting line is 45 meters. Then take the midpoint of the connecting line as the vertex, calculate the included angle formed by the two end points and the midpoint of the connecting line. Specifically, through the elevation difference and horizontal distance between two points, the slope angles of point A to the midpoint and point B to the midpoint are calculated using the inverse tangent function. The difference between the two slope angles is the terrain shielding angle, for example, the slope angle of point A to the midpoint is 10°, and the slope angle of point B to the midpoint is 8°. The terrain shielding angle is 2°. This angle reflects the shielding degree of the terrain undulation to the connecting line of the two points. The larger the angle, the more significant the shielding. The acquisition of the terrain shielding coefficient considers the elevation difference, horizontal distance and terrain shielding angle of the connecting line direction:

[0049] First, the height difference between two points is calculated, such as the height difference between points A and B is 10 meters and the horizontal distance is calculated by geographic coordinates, such as 100 meters, to get the relative height difference ratio (10 / 100=0.1), and then combined with the terrain shielding angle (such as 2°), the terrain shielding coefficient is calculated by the formula: 0.5 x relative height difference ratio + 0.5 x (terrain shielding angle / 90°), for example, the coefficient in the above example is 0.5 x 0.1 + 0.5 x (2° / 90°)≈0.05 + 0.011=0.061, the larger the coefficient value, the stronger the influence of terrain on the shielding of the line between the two points, when constructing the distance inverse weight matrix, first take the reciprocal of the spatial distance of the point pair in the feature point subset as the initial weight, such as the horizontal distance between two points is 100 meters, the initial weight is 1 / 100=0.01, then multiply the initial weight by the inverse of the terrain shielding coefficient (1 / 0.061≈16.39) for weighted correction, the modified weight is 0.01 x 16.39≈0.164. Perform the same calculation on all point pairs in the feature point subset, and arrange the modified weights in order of point pairs to form a distance inverse weight matrix, for example, a feature point subset contains 3 pairs of points, the modified weights are 0.164, 0.21, and 0.18 respectively, the matrix is constructed with these three values as elements, ensuring that point pairs with strong terrain shielding effects obtain higher weights in subsequent matrix solving, making the affine transformation matrix more suitable for image matching under actual terrain conditions, providing a more accurate weight basis for subsequent optimization of the affine transformation matrix, making the compensation and calibration of radar images more in line with the complex terrain features of islands, and thus improving the accuracy of cross-sensor image matching.

[0050] When introducing the L1 norm to optimize the distance inverse weight matrix, an adaptive smoothing constraint mechanism is adopted, specifically:

[0051] After constructing the distance inverse weight matrix of the terrain shielding coefficient, the L1 norm is 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 the adaptive smoothing constraint mechanism, which not only connects the quantitative results of the weight matrix, but also provides a mathematical basis for the accurate compensation of the radar image, and forms a close link with the subsequent image calibration link. The L1 norm refers to the sum of the absolute values of each element in the vector, which is used to measure the overall deviation of the weighted displacement residual in this method. 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 is used to describe the geometric transformation relationship between the radar image and the optical image. Its source is based on the coordinate correspondence of the point pairs in the feature point subset. Through mathematical modeling, the spatial alignment of the two types of images is realized. The matrix contains three types of parameters: rotation, translation and scaling. The rotation parameter reflects the angle of rotation of the image around the origin (in degrees), the translation parameter represents the offset of the image in the x and y directions (in meters), and the scaling parameter represents the scaling ratio of the image in the horizontal and vertical directions (unitless, 1 represents equal scaling). When initializing these parameters, the coarse alignment + empirical value correction strategy is adopted:

[0052] 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. According to the overall trend of the island coastline, such as the east-west trend, 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 the optical image and 2 meters / pixel for the radar image. The scaling parameter is initialized to 0.5. The distance inverse weight matrix is used as the residual weighting basis because it has integrated the terrain shielding and spatial distance factors, which can reflect the matching reliability of different point pairs. The higher the weight of the point pair, such as the point pair with small terrain shielding and short distance, the higher the matching accuracy, which should be given more weight in the residual calculation. The residual weighting basis refers to weighting the displacement residual of each point pair with the elements in the weight matrix, so that the residual of reliable point pairs accounts for a higher proportion in the objective function. When constructing the objective function, a dynamic residual threshold is introduced:

[0053] For each point pair, the difference between the predicted coordinates of the radar feature points after affine transformation and the actual coordinates of the optical feature points, i.e. the displacement residual, is calculated, the residual is multiplied by the weight of the corresponding point pair in the distance inverse 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, the weight of a certain point pair is 0.2, the residual is 3 meters, the weighted residual is 0.6, and the residual of another point pair with a weight of 0.1 is 5 meters, the weighted residual is 0.5, and the objective function will preferentially reduce the residual of the high-weight point pair. Due to the possible unevenness of the spatial distribution density of the feature point subset, the smoothing constraint strength needs to be dynamically adjusted according to the density, and this step is realized through the density-constraint mapping:

[0054] The smoothing constraint strength refers to the coefficient that controls the degree of change in the affine transformation parameters. The larger the value, the more gradual the parameter change, avoiding over-fitting to local point pairs. First, calculate the point pair density in each 100m x 100m grid. For example, a certain grid contains 8 point pairs, with a density of 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 in dense areas and prevents over-distortion in sparse areas, naturally linking the construction and solution process of the objective function. Finally, solve the optimal affine transformation matrix through the iterative reweighted least squares method:

[0055] First, calculate the weighted residual of each point pair with the initial parameters, and give smaller iteration weights to high-weight residual points. For example, point pairs with a residual greater than twice the average value have a weight 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, and check the change in the objective function value after each iteration. Stop iteration 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.

[0056] The gradient domain fusion compensation is used to perform coordinate offset compensation on the radar satellite, including the following steps:

[0057] After obtaining the optimal affine transformation matrix, the radar image needs to be calibrated through gradient domain fusion compensation. This process is the key link of converting the mathematical transformation model into actual image coordinate correction, which not only connects the solution of the affine transformation matrix in the foregoing, but also improves the matching accuracy of the radar image and the optical image through global and local dual compensation, laying a foundation for subsequent secondary compensation calibration. When the affine transformation matrix is decomposed into a global affine component and a local elastic component, a matrix hierarchical analysis method is adopted:

[0058] The global affine component corresponds to the linear transformation part in 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 nonlinear deviation part remaining after the global affine component is removed from the affine transformation matrix, which reflects the local subtle deviation caused by terrain undulation and imaging distortion. When decomposing, the rotation angle, translation amount and scaling factor in the matrix are first extracted as the global affine component, and then the deviation matrix of the local elastic component is obtained through matrix subtraction to realize the separation of the two. The global affine component is resampled by bilinear interpolation to obtain the global displacement field, as follows:

[0059] First, the radar image is divided into a 10x10 pixel grid. According to the parameters of the global affine component, the target coordinates of each grid vertex in the optical image coordinate system are calculated. For example, a certain 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 by bilinear interpolation algorithm, i.e. according to the target coordinates of the four vertices of the grid, the displacement value of each pixel in the grid is calculated by weighting according to its relative position 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, making the image roughly aligned with the optical image. For the local elastic component, the gradient field of the radar image is extracted and Poisson reconstructed to generate the local displacement field, as follows:

[0060] First, the gradient values of the radar image in the x and y directions are calculated by the Sobel operator, which reflects the rate and direction of pixel gray scale change, forming a gradient field. Then, the gradient field is Poisson reconstructed with the deviation matrix of the local elastic component as the constraint, i.e. by solving the Poisson equation, the gradient field is converted into a continuous displacement distribution, so that the reconstructed displacement field not only satisfies the gradient constraint, but also fits the nonlinear deviation of the local elastic component. For example, in the edge area of the reef, the gradient field shows a sharp change, and the Poisson reconstruction generates a local stretching displacement to correct the edge distortion of the radar image caused by the imaging angle.

[0061] Firstly, according to the terrain complexity of the image area, dynamic weights are assigned to the global and local displacement fields. For 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. For complex terrain areas (such as multi-reef areas), the local displacement field weight is increased to 0.7, and the global displacement field weight is 0.3, to strengthen local detail correction. After superposition, the comprehensive displacement field is obtained, and then the coordinate offset compensation is performed for each pixel of the radar image. First, the target coordinates of the pixel are calculated according to the comprehensive displacement field, and then the sub-pixel level offset + edge sharpening strategy is used. For the offset pixel, the details are preserved through sub-pixel interpolation, and at the same time, the key edge areas such as the coastline are sharpened to avoid edge blur after compensation. For example, a reef pixel needs to be offset by 2.3 pixels according to the comprehensive displacement field calculation, and through sub-pixel interpolation, it is accurately positioned to the target position, and the edge contrast with the sea water is strengthened. The final generated calibrated radar image can be highly matched with the optical image in both overall and local, and the alignment error of the coastline feature points can be controlled within 2 pixels, providing high-quality basic data for subsequent secondary compensation calibration of new optical images, effectively improving the data consistency of cross-sensor images in the monitoring of uninhabited islands.

[0062] When the cloud cover probability and the difference value are both lower than the corresponding preset threshold values:

[0063] The new optical image is reacquired, and the average Euclidean distance error of the coastline feature points of the calibrated radar image and the new optical image is calculated. According to the average Euclidean distance error, the calibrated radar image is judged for secondary compensation calibration, and the time-series monitoring data of the entire uninhabited island is generated.

[0064] The average Euclidean distance error of the coastline feature points of the calculated calibrated radar image and the new optical image is calculated through multi-level confidence control, specifically:

[0065] After the gradient domain fusion compensation of the radar image is completed, 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, and 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 the possible secondary compensation, and 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 level or scattering characteristics from the surrounding sea area, such as sudden change in gray level of the coastline pixels in the optical image, and strong scattering signal in the radar image, which together form a continuous pixel set reflecting the coastline morphology. When the coastline edge pixel chain is sampled to generate a verification point set with equal arc length intervals, the total arc length of the pixel chain is calculated first, and then the pixels are converted to actual distances based on geographic coordinates. For example, the total arc length is 500 meters, and the sampling interval is determined based on equal arc length (e.g. 10 meters). Starting from the beginning of the pixel chain, every 10 meters, a pixel is selected as a verification point. If the end point is less than 10 meters, it is combined with the last verification point. Finally, 50 verification points are generated to form a verification point set, ensuring that the sampling points evenly cover the entire coastline and avoiding local dense or sparse errors. After taking the coordinates of each verification point in the new optical image as the reference, a pixel search window is established at the corresponding position in the radar image:

[0066] First, the latitude and longitude coordinates of the optical image verification point are mapped to the pixel coordinate system of the radar image through geographic coordinate conversion to obtain the initial corresponding position. Considering the possible residual error, a 11x11 pixel square search window is established around the center of the position, i.e. expanding 5 pixels in all directions. The window size can be dynamically adjusted based on the calibration accuracy, such as expanding to 15x15 pixels if the calibration error is large, to ensure that the window contains the true matching point. When calculating the pixel with the highest spectral similarity in the search window as the matching point, the 3x3 neighborhood gray feature vector of the optical reference point is extracted, such as mean, variance, gradient direction. The dual-polarization scattering feature vector (such as VV / VH polarization ratio, scattering entropy value) of each pixel in the search window of the radar image is extracted, and the matching degree of the two feature vectors is calculated by cosine similarity algorithm (value 0-1, 1 represents complete matching). The pixel with the highest similarity is selected as the matching point. Then the Euclidean distance between the two is calculated based on their geographic coordinates, such as the coordinates of the optical reference point (x1, y1) and the coordinates of the radar matching point (x2, y2), the distance is 2 +(y1-y2) 2The distance value is recorded as the error of the verification point, when the Euclidean distance values of all verification points are analyzed by box plot to eliminate abnormal values, first, all distance values are sorted from small to large, the quartiles are calculated, Q1 is the 25th percentile, Q3 is the 75th percentile, the interquartile range IQR = Q3-Q1 is determined, the abnormal value judgment threshold is Q1-1.5*IQR and Q3+1.5*IQR, and the verification point whose distance value is out of the range is determined as an abnormal value, for example, if the distance value of a certain verification point is 15 meters, which is much higher than 8 meters of Q3+1.5*IQR, it is determined as an abnormal value, after eliminating these abnormal values, the remaining distance values are considered as reliable matching errors, finally, the remaining Euclidean distance values are weighted and averaged according to the type of the coast, the coastline is divided into types such as rock coast, sandy beach coast and artificial coast, the weight is distributed according to the proportion of each type in the total coastline, the average value of the distance value in each type is calculated, and then the sum is weighted to generate the final average Euclidean distance error, for example, the average error of the rock coast is 2 meters, the average error of the sandy beach coast is 3 meters, and 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, which provides a scientific basis for judging whether secondary compensation is needed, and further guarantees the consistency of the global time series monitoring data.

[0067] The application obtains target island optical images and cloud data through optical and meteorological satellites, predicts cloud shielding probability, dispatches radar satellites when the threshold is exceeded, compares the effective area of the optical image with the reference requirement, triggers radar imaging when the difference reaches a preset proportion, extracts a set of coastline feature points from historical images, calculates spatial distribution entropy to screen a subset, constructs an optimized affine transformation matrix to compensate the radar image, obtains new optical images when the cloud cover is low, calculates the coastline feature point error between the calibrated radar image and the new optical image, performs secondary compensation, and generates global time series monitoring data, thereby improving the monitoring continuity and data consistency in the cloud shielding scenario, and being suitable for dynamic monitoring of uninhabited islands.

[0068] The basic principles, main features and advantages of the present application are shown and described above. Those skilled in the art should understand that the present application is not limited by the above examples, the above examples and descriptions in the specification are only preferred examples of the present application, and are not intended to limit the present application, various changes and improvements can be made to the present application without departing from the spirit and scope of the present application, and these changes and improvements all fall within the scope of the claimed present application. The scope of protection of the present application is defined by the appended claims and their equivalents.

Claims

1. A method for monitoring uninhabited islands by satellite remote sensing, characterized in that: The method comprises the following steps: The three-dimensional cloud field dynamic model is constructed, and the cloud shielding probability and the difference between the effective area of the optical image and the reference required area are outputted. When the cloud shielding probability and the difference respectively reach the corresponding preset threshold, the radar satellite is driven to obtain the radar image; When the cloud shielding probability and the difference respectively exceed the corresponding preset threshold: The historical optical image coastline feature point set and the radar image feature point set are respectively extracted from the historical optical image and the radar image. The spatial distribution entropy values of the historical optical image coastline feature point set and the radar image feature point set are calculated. The feature point subset is selected according to the spatial distribution entropy values. The distance inverse weight matrix is constructed based on the feature point subset. The L1 norm is introduced to optimize and solve the distance inverse weight matrix. The affine transformation matrix is obtained. The coordinate offset compensation is performed on the radar satellite by using the affine transformation matrix. The calibrated radar image is obtained. When the cloud shielding probability and the difference respectively are lower than the corresponding preset threshold: The new optical image is reacquired. The average Euclidean distance error of the coastline feature points of the calibrated radar image and the new optical image is calculated. The calibrated radar image is secondarily compensated and calibrated according to the average Euclidean distance error. The time sequence monitoring data of the entire area of the uninhabited island is generated. 2.The method for monitoring uninhabited islands by satellite remote sensing according to claim 1, characterized in that: The cloud shielding probability is predicted, specifically comprising the following steps: The infrared band cloud top temperature of the meteorological satellite, the visible band cloud layer albedo and the sea surface humidity of the ground meteorological station are synchronously accessed. The three-dimensional cloud field dynamic model is constructed. The real-time cloud field data are obtained. The monsoon vortex feature vector of the target island sea area is extracted based on the historical cloud motion trajectory library. The real-time cloud field data are dynamically time warping matched with the monsoon vortex feature vector. The probability value of the cloud layer covering the target island sea area is calculated. When the probability value exceeds the preset probability threshold, the radar satellite scheduling instruction is generated. The monsoon vortex feature vector includes the cloud layer moving direction angular velocity, the cloud cluster fragmentation index and the vertical development gradient. 3.The method for monitoring uninhabited islands by satellite remote sensing according to claim 1, characterized in that: The tidal interference-resistant stability feature extraction algorithm is used to extract the historical optical image coastline feature point set, specifically comprising the following steps: The reef vertex, the sea cliff base turning point and the permanent artificial feature corner on the historical optical image are identified. The tidal insensitive feature set is formed. The spatial density clustering analysis is performed on each type of feature point in the tidal insensitive feature set. The feature points with the spatial density lower than the preset threshold are removed. The main direction of the remaining feature points is selected through the edge gradient direction consistency test. Finally, the historical optical image coastline feature point set is generated.

4. The method for monitoring uninhabited islands by 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, specifically comprising the following steps: The sea clutter suppression filter is performed on the VV polarization and VH polarization channels of the radar image respectively. The dual-polarization image of the radar image is obtained. The strong scattering target points are located on the dual-polarization image. The polarization scattering entropy value and the scattering angle of each strong scattering target point are calculated. The target points meeting the preset conditions of the polarization scattering entropy value and the scattering angle are selected as candidate feature points. The candidate feature points are projected to the geographic coordinate system and matched with the historical optical image coastline feature point set in the spatial range. The candidate feature points in the spatial range matching area are retained and constitute the radar image feature point set.

5. The method for monitoring uninhabited islands by satellite remote sensing according to claim 4, characterized in that: The spatial distribution entropy value is calculated by introducing the coastline shape constraint factor, including 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 curve segments according to the curvature change; The Euclidean space grid division method is used in the straight line segment to calculate the Shannon entropy value of the feature points in each Euclidean space grid, and the Frechet distance matching method is used in the curve segment to calculate the shape deviation entropy value of the feature point sequence and the standard curve template, and finally the weighted sum of the Shannon entropy value and the shape deviation entropy value is taken as the spatial distribution entropy value.

6. The method for monitoring uninhabited islands by satellite remote sensing according to claim 5, characterized in that: The feature point subset is screened through bidirectional stability verification, including the following steps: The spatial distribution entropy difference value of the historical optical image coastline feature point set and the radar image feature point set is calculated, and the point pairs that meet the difference threshold condition are retained, and the geometric invariance verification is performed on the point pairs, that is: Each point pair is taken as a control point to calculate the local affine transformation residual, and the point pairs with residual greater than the preset pixel are removed, and the time invariance verification is performed on the remaining point pairs, the coordinate drift amount of the point pair in the historical three-period image is compared, and the verification point pair is determined according to the standard deviation of the coordinate drift amount, and finally the point pairs that pass the verification constitute the feature point subset.

7. The method for monitoring uninhabited islands by satellite remote sensing according to claim 6, characterized in that: The distance inverse ratio weight matrix is constructed by fusing the terrain shielding coefficient, including the following steps: The digital elevation model data of the target island is obtained, and the terrain shielding angle of the feature point subset is calculated, the reciprocal of the spatial distance of the point pair in the feature point subset is taken as the initial weight, and the inverse ratio of the terrain shielding coefficient is multiplied to correct the weight, and the distance inverse ratio weight matrix is obtained, wherein the terrain shielding coefficient is determined by the elevation difference, the horizontal distance and the terrain shielding angle in the connection direction.

8. The method for monitoring uninhabited islands by satellite remote sensing according to claim 7, characterized in that: When the L1 norm is introduced to optimize and solve the distance inverse ratio weight matrix, an adaptive smoothing constraint mechanism is adopted, including the following steps: The rotation, translation and scaling parameters of the affine transformation matrix are initialized, the distance inverse ratio weight matrix is taken as the residual weighted basis, and the target function is constructed to minimize the sum of the weighted displacement residual absolute values; 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 by satellite remote sensing according to claim 8, characterized in that: The gradient domain fusion compensation is used to perform coordinate offset compensation on the radar satellite, including the following steps: The affine transformation matrix is decomposed into global affine component and local elastic component, the global affine component is resampled by bilinear interpolation to obtain the global displacement field, the gradient field of the radar image is extracted, the local displacement field is generated by Poisson reconstruction in the gradient field, the global displacement field and the local displacement field are superimposed, and the coordinate offset compensation is performed on the radar satellite to generate the calibrated radar image.

10. The method for monitoring uninhabited islands by satellite remote sensing according to claim 9, characterized in that: The average Euclidean distance error of the coastline feature points of the calibrated radar image and the new optical image is calculated through multi-level confidence control, including the following steps: ​ The coastline edge pixel chain is extracted on the new optical image and the calibration radar image, the verification point set is generated by sampling the coastline edge pixel chain with equal arc length interval, the coordinate of each verification point set on the new optical image is taken as the reference, the pixel search window is established at the corresponding position of the radar image, the point with the highest spectral similarity to the reference point in the pixel search window is taken as the matching point, the Euclidean distance between the two points is recorded, the box plot analysis is performed on the Euclidean distance values of all verification points, the abnormal values are removed, and the average Euclidean distance error is generated by weighted averaging the remaining Euclidean distance values according to the coast type.

Images