Sea area monitoring method and system based on satellite remote sensing technology

Abstract

The invention discloses a sea area monitoring method and system based on a satellite remote sensing technology, and relates to the technical field of remote sensing, and the method comprises the steps: firstly, generating multi-stage and multi-band remote sensing raster data through a resource set, and secondly, obtaining joint sample vectors at different times through monitoring external environment data in a target monitoring region, a multivariable kernel density estimation algorithm is adopted to obtain a four-dimensional joint probability density distribution function, and an abnormal coverage area is obtained in combination with multi-stage and multi-band remote sensing raster data; then, based on the abnormal coverage area, identifying an enhanced abnormal unit, and performing dynamic correction operation; and finally, after dynamic correction, obtaining an abnormal connected region, analyzing an abnormal unit with an unstable form, determining an abnormal evolution sensitive region, and realizing a local adaptive optimization means based on the change of the abnormal evolution sensitive region.

Description

Technical Field

[0001] The present invention relates to the field of remote sensing technology, and in particular to a sea area monitoring method and system based on satellite remote sensing technology. Background Art

[0002] This invention relates to the field of remote sensing technology, specifically to a method for marine monitoring based on satellite remote sensing technology, primarily used for automated remote sensing monitoring of oil spills, illegal discharges of pollutants, and marine ecological anomalies. Within the broad technical field of marine information perception and monitoring, satellite remote sensing, due to its wide-range, periodic, and multi-band acquisition of physical parameters, has become an important technical means for monitoring marine pollution and abnormal ecological processes. Within specific branches of satellite remote sensing technology, multi-source remote sensing imagery, such as synthetic aperture radar (SAR), visible multispectral, and thermal infrared, is widely used to identify and track pollution incidents such as oil spills and sewage discharges.

[0003] Currently, monitoring of abnormal targets in these waters primarily relies on single-period synthetic aperture radar (SAR) imagery for anomaly identification. While SAR offers the advantage of all-weather, all-day imaging, complex wind and wave conditions often mask backscatter variations caused by abnormal oil slicks, minor pollution discharges, or incipient eutrophication. For example, at high wind speeds (greater than 10 m / s), surface spray and foam significantly enhance the scattering signal, preventing thin oil slicks from forming the typical "black spot" characteristic in SAR images, leading to missed detections. Conversely, at low wind speeds (less than 3 m / s), the mirror-like reflection effect of a calm sea surface can also create dark spots, which can be confused with oil spills and lead to misjudgments. This monitoring approach, which relies on a single radar band, has significant limitations in practical implementation and ecological risk management. Summary of the Invention

[0004] In view of the shortcomings of the existing technology, the present invention provides a sea area monitoring method and system based on satellite remote sensing technology, which solves the problems in the above-mentioned background technology.

[0005] To achieve the above objectives, the present invention is implemented through the following technical solutions: A method for monitoring sea areas based on satellite remote sensing technology, comprising the following steps: S1: Generate multi-period and multi-band remote sensing raster data by real-time monitoring of resource sets in the sea area; S2: By monitoring the external environmental data within the target monitoring area, joint sample vectors at different times are obtained. A multivariate kernel density estimation algorithm is used to obtain a four-dimensional joint probability density distribution function. Combined with multi-period and multi-band remote sensing raster data, the abnormal coverage area is obtained. S3: Based on the abnormal coverage area, identify and enhance abnormal units and perform dynamic correction operations; S4: After dynamic correction, the abnormal connected area is obtained, and the abnormal units with unstable morphology are analyzed to determine the abnormal evolution sensitive area. Based on the changes in the abnormal evolution sensitive area, local adaptive optimization means are implemented.

[0006] Preferably, S1 includes: S11: Using orbit prediction tools to obtain all satellite resources used for sea area monitoring within a predetermined monitoring period, and obtain a resource set; S12: Based on the available band capabilities and orbit transit timing of each satellite in the resource set, combined with the geometric coverage relationship of the entire sea area, greedy scheduling is performed through the time-sorted queue to generate a multi-satellite multi-mode schedule list; S13: Using a multi-satellite multi-mode scheduling list, SAR images, visible multispectral images, and thermal infrared images acquired from different satellites at different times are acquired to generate an image library. The images in the image library are then geometrically registered. Specifically, the image geometric registration tool is used to unify the images in the image library into the same map projection coordinate system and resample them to a uniform spatial resolution. The size of the sea area grid is determined based on the spatial resolution. S14: After completing the geometric registration, the pixel values ​​of the same location at different times and in different bands are organized into a four-dimensional matrix to form multi-period multi-band remote sensing raster data; Among them, the time sorting queue is a time window formed by sorting the satellites from early to late according to their orbital transit timing.

[0007] Preferably, S2 includes: S21: Obtain external environmental data within the target monitoring area through the buoy array, including wind speed, wind direction, wave height, and shortwave spectrum peak energy. Perform spatial interpolation on the external environmental data according to the sea area grid points to align it with the multi-period multi-band remote sensing grid data. Simultaneously, obtain a multi-period external environmental sample set, wherein the multi-period external environmental sample set includes joint sample vectors at different times. S22: Based on a multi-period external environment sample set, a multivariate kernel density estimation algorithm is used to perform kernel density estimation on the joint sample vector to obtain a four-dimensional joint probability density distribution function. The four-dimensional joint probability density distribution function is used to quantify the probability density value of any joint sample vector appearing simultaneously in historical samples.

[0008] Preferably, S2 further includes: S23: According to the four-dimensional joint probability density distribution function, the sum of all historical probability density values ​​lower than the current joint sample vector is counted to obtain a cumulative distribution function value. The cumulative distribution function value is used to quantify the cumulative occurrence probability of the current joint sample vector in the historical samples; S24: using the cumulative distribution function value, obtaining the abnormal sensitivity value based on wind field driving corresponding to the current monitoring time and spatial position, the abnormal sensitivity value is used to reflect the rarity of the currently monitored joint sample vector at the corresponding geographical location; S25: Determine the SAR backscatter coefficient corresponding to each pixel point in the SAR image; S26: At each sea area grid point, the local mean and standard deviation of the SAR backscatter coefficient within a sliding window centered at each sea area grid point are calculated, and the backscatter coefficient of each sea area grid point is normalized to obtain a scattering normalized outlier value; S27: Couple the anomaly sensitivity value with the scattering normalized anomaly value to calculate the scattering anomaly confidence value at each sea area grid point, specifically: ,in, is the confidence value of scattering anomaly, is the abnormal sensitivity value, is the scatter normalized outlier; S28: Construct a scattering anomaly confidence surface based on the scattering anomaly confidence values ​​at each sea area grid point; S29: Generate a dynamic local threshold surface based on the adaptive quantile method. By comparing the scattering anomaly confidence surface with the dynamic local threshold surface, an anomaly mask is formed. Specifically: ,in, For the location The abnormal mask on For the location The confidence value of the scattering anomaly on For the location Dynamic local threshold on ; S291: The sea area grid points with an anomaly mask of 1 are regarded as anomaly units, and the anomaly coverage area is obtained through statistics.

[0009] Preferably, S3 includes: S31: Based on the scattering anomaly confidence value, calculate the inverse sensitivity weight of each sea area grid point in the target monitoring area; S32: Determine weakly significant units based on the abnormal units, and extract the abnormal values ​​of the sea water color index, eutrophication index and thermal infrared brightness temperature within the weakly significant units; S33: Using the inverse sensitivity weight as the weighting factor, the local multivariate regression model is fitted by the weighted least squares method to establish a dual-domain photothermal prediction model based on the sea water color index, eutrophication index and thermal infrared brightness temperature anomaly values ​​to calculate the photothermal residual value on each weakly significant unit.

[0010] Preferably, S3 further includes: S34: At each sea area grid point, based on the light and thermal residual value, the 95th percentile value of the light and thermal residual value in the sliding window centered on the corresponding sea area grid point is used as the local dynamic threshold. If the light and thermal residual value exceeds the local dynamic threshold, the corresponding sea area grid point is marked as an enhanced anomaly unit; S35: Dynamically modify the abnormal coverage area by enhancing the abnormal unit to obtain the final abnormality detection mask in the target monitoring area.

[0011] Preferably, S4 includes: S41: Obtain multiple abnormal connected areas according to the final anomaly detection mask in the target monitoring area; S42: extracting the area and perimeter of each abnormal connected region and calculating the compactness, which is used to locate abnormal unit regions with unstable morphology; S43: When the compactness is lower than a preset compactness threshold, the corresponding abnormally connected area is marked as an abnormal evolution sensitive area, otherwise no marking is performed.

[0012] Preferably, S4 further includes: S44: After monitoring, obtain the abnormal evolution sensitive area in the next monitoring cycle, compare the deviation values ​​of the abnormal units in the abnormal evolution sensitive area during each monitoring cycle, and if the deviation value exceeds the preset deviation threshold, automatically repeat S2 to correct the abnormal sensitivity value driven by the wind field and dynamically adjust the scattering anomaly confidence value and inverse sensitivity weight.

[0013] A sea area monitoring system based on satellite remote sensing technology, including: The data module is used to generate multi-period and multi-band remote sensing raster data by real-time monitoring of resource sets in the sea area; The preliminary identification module is used to obtain joint sample vectors at different times by monitoring the external environmental data within the target monitoring area. The multivariate kernel density estimation algorithm is used to obtain the four-dimensional joint probability density distribution function. Combined with multi-period and multi-band remote sensing raster data, the abnormal coverage area is obtained. The secondary identification module is used to identify and enhance abnormal units based on abnormal coverage areas and perform dynamic correction operations; The local optimization module is used to obtain abnormal connected areas after dynamic correction, analyze abnormal units with unstable morphology, determine the abnormal evolution sensitive areas, and implement local adaptive optimization methods based on the changes in the abnormal evolution sensitive areas.

[0014] The present invention provides a method and system for monitoring sea areas based on satellite remote sensing technology, which has the following beneficial effects: (1) Through step S1, multiple satellite resources are used to obtain multiple remote sensing images at different times and different bands, unify the geometric registration and resolution, and form multi-period and multi-band remote sensing raster data, so that the monitoring area is fully covered in terms of time, band and spatial distribution, significantly enhancing the traceability and multi-dimensional comparison capabilities of abnormal water changes. Through step S2, based on the collection of multiple external environmental data (such as wind speed, wind direction, wave height and shortwave spectrum peak energy) and the acquisition of joint sample vectors at different times, a multivariate kernel density estimation algorithm is used to establish a four-dimensional joint probability density distribution function, effectively characterizing the rarity of the current monitoring environment in the historical data distribution, and then combining remote sensing data to form an abnormal coverage area. This method can adaptively identify potential anomalies under different wind and wave fields, and improve the detection accuracy under complex environmental conditions. Through step S3, based on the abnormal coverage area, the abnormal unit is further identified and enhanced, and the optical and thermal multi-domain residual method is used for dynamic correction, reducing the SAR scattering blind area problem caused by high wind speed or surge, thereby achieving effective detection of hidden oil film or weak pollution. Through step S4, the dynamically corrected abnormal connected areas are further analyzed, and their morphological parameters such as area, perimeter, and compactness are combined to identify morphologically unstable abnormal unit areas and determine the sensitive areas for abnormal evolution. Based on the periodic changes in the sensitive areas, the anomaly detection parameters are dynamically updated to complete the local adaptive optimization of the anomaly sensitivity and threshold model, effectively overcoming the problem of frequent fluctuations in anomaly detection boundaries in long-term series monitoring. In summary, this method forms a closed-loop monitoring system from data acquisition to model evolution through multi-period and multi-band multimodal data fusion, probability density-based environment-driven adaptive detection, and a local optimization mechanism with dynamic feedback. This improves the detection accuracy and spatial and temporal stability of abnormal water phenomena such as oil spills and eutrophication.

[0015] (2) In steps S23 and S24, the cumulative distribution function value CDF is obtained by integrating the four-dimensional joint probability density distribution function to quantify the cumulative occurrence probability of the current joint sample vector (i.e., the combination of wind speed, wind direction, wave height and shortwave spectrum peak energy) in the historical samples, and further inverted into an anomaly sensitivity value. The anomaly sensitivity value reflects the rarity of the current environmental combination, so that subsequent anomaly detection automatically improves the sensitivity under rare environmental conditions and reduces the missed detection rate under special meteorological-sea condition coupling conditions. Through S25 to S29, the SAR backscatter coefficient is first extracted at each sea area grid point, and then normalized by the mean and standard deviation of the local sliding window in S26 to obtain a scattering normalized anomaly value that can measure the degree of deviation between the water surface and the local background. Then, in S27, the anomaly value is coupled with the anomaly sensitivity value to calculate the scattering anomaly confidence value. In S28, a global scattering anomaly confidence surface is formed, and compared with the dynamic local threshold surface generated by the adaptive quantile method in S29 to form an anomaly mask. Finally, in S291, the anomaly unit is extracted and the anomaly coverage area is statistically obtained. Through this multi-stage sensitive detection mechanism that couples environmental rarity and local anomaly intensity, the ability to capture anomalies (such as weak oil film, slight pollution discharge or early eutrophication) under complex wind and wave conditions is significantly improved. In summary, through the multi-stage collaborative processing of S21 to S291 of the above-mentioned S2: dynamic adaptive anomaly sensitivity generation driven by the environment is realized, while maintaining the detection sensitivity, the spatial consistency of the anomaly detection results with changes in wind field, wave height and shortwave spectrum is significantly enhanced, so that the final detection results can maintain high stability and reliability under variable sea conditions, providing a solid foundation for subsequent tracking of the stability of abnormal unit morphology and identification of sensitive areas of abnormal evolution, and effectively breaking through the technical bottleneck of high false alarms and high missed detection rates of traditional methods in high dynamic environments.

[0016] (3) Through the inverse sensitivity weight guidance of S31 to S35, the residual analysis of the optical and thermal dual-domain fitting, and the dynamic quantile enhancement detection, the method of the present invention realizes the active attention allocation to the low-confidence units, the secondary deep excavation of the SAR weak anomaly area in the optical and thermal domain, and the adaptive correction of the anomaly coverage area, thereby improving the detection capability of hidden oil film, minor pollution discharge and early eutrophication in complex wind fields and multi-period sea conditions, and enhancing the refinement of the spatial distribution of the monitoring results and the detection sensitivity.

[0017] (4) In step S43, when the compactness of a certain abnormal connected area is detected to be lower than the preset compactness threshold, it can be determined that the boundary morphology of the area is loose and complex, and it is easy to have drastic changes in morphology under different monitoring cycles, and produce significant offsets with wind and wave disturbances. At this time, the connected area is marked as an abnormal evolution sensitive area, otherwise it is not marked. In this way, active spatial focusing on the morphologically changeable areas in the abnormal detection results is achieved, which is convenient for key tracking and analysis in subsequent cycles. In step S44, after obtaining the data of the next monitoring cycle, for these abnormal evolution sensitive areas, the deviation value is calculated by comparing the spatial coverage change of the abnormal unit in the current cycle detection result with the historical average detection result; if the deviation value exceeds the preset deviation threshold, it means that the detection result of the area is significantly inconsistent with the historical pattern, which may be due to the fact that the environmental statistical model or sensitivity calculation has not fully learned the long-term wind and wave-anomaly relationship here. At this time, S2 is automatically repeated, and the four-dimensional joint probability density distribution function is locally optimized based on the latest external environmental data and detection results in the area, and the abnormal sensitivity value, scattering anomaly confidence value and inverse sensitivity weight are dynamically updated to form a closed-loop adaptive adjustment mechanism. In summary, through S41 to S44, it is possible to infer the abnormally volatile areas from the morphological characteristics of the detection results, and automatically lock the areas with significant spatial fluctuations in the time series dimension, forming a key dynamic monitoring mechanism. It can continuously self-learn and gradually optimize the wind field-driven anomaly detection model in the long-term monitoring sequence, effectively avoiding the problems of frequent jumps in anomaly detection boundaries, false alarms or omissions caused by fixed thresholds and lack of spatiotemporal adaptability in traditional monitoring. BRIEF DESCRIPTION OF THE DRAWINGS

[0018] Figure 1 This is a flow chart of a method for sea area monitoring based on satellite remote sensing technology according to the present invention; Figure 2 This is a partial logic diagram of a sea area monitoring method based on satellite remote sensing technology of the present invention; Figure 3 This is a partial logic diagram of a sea area monitoring method based on satellite remote sensing technology of the present invention; Figure 4 This is a block diagram of a sea area monitoring system based on satellite remote sensing technology in the present invention. DETAILED DESCRIPTION

[0019] 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.

[0020] Embodiment 1: See also Figures 1 to 3The present invention provides a sea area monitoring method based on satellite remote sensing technology, comprising the following steps: S1: Generate multi-period and multi-band remote sensing raster data by real-time monitoring of resource sets in the sea area; S2: By monitoring the external environmental data within the target monitoring area, joint sample vectors at different times are obtained. A multivariate kernel density estimation algorithm is used to obtain a four-dimensional joint probability density distribution function. Combined with multi-period and multi-band remote sensing raster data, the abnormal coverage area is obtained. S3: Based on the abnormal coverage area, identify and enhance abnormal units and perform dynamic correction operations; S4: After dynamic correction, the abnormal connected area is obtained, and the abnormal units with unstable morphology are analyzed to determine the abnormal evolution sensitive area. Based on the changes in the abnormal evolution sensitive area, local adaptive optimization means are implemented.

[0021] Existing SAR (Synthetic Aperture Radar) satellite-based ocean monitoring is commonly used for oil spill detection and ice floe boundary identification. However, because SAR detection relies on the scattering properties of rough, shortwave units on the sea surface, low (<3 m / s) or high (>10 m / s) wind speeds can result in false positives or negatives in the image, distorting monitoring conclusions. This problem is particularly evident in the following ways: at low wind speeds, the sea surface is naturally smooth, making it easy to misinterpret oil slicks; at high wind speeds, the oil slick's wave suppression effect is lost, making it difficult to identify and potentially leading to missed detections. Therefore, a method for ocean monitoring based on satellite remote sensing technology is proposed.

[0022] In this embodiment, the present invention provides a method for monitoring sea areas based on satellite remote sensing technology. By leveraging deep coupled analysis of multi-period, multi-band remote sensing data and external driving data such as wind fields, this method effectively addresses the challenges of existing single-period, single-source monitoring models, such as poor detection stability, high false alarm rates, and susceptibility to masking by short-period wind disturbances, when faced with complex wind and wave conditions, variable oil slicks, and abnormally nutrient-rich water bodies. Specific beneficial effects are reflected in the following aspects: In step S1, the orbital timing, band capabilities, and sea area coverage geometry of multiple satellites with concentrated resources are dynamically scheduled. For example, within a 5-day monitoring cycle, a greedy scheduling algorithm is used to schedule multi-satellite transit missions such as Sentinel-1 / 2, GF-3, and Landsat-8, so that a distributed multi-period data grid sequence is formed in the time dimension for the target sea area.

[0023] Subsequently, through image geometric registration and unified spatial resolution (e.g., 30m grid) processing, the observation values ​​of different times and different sensors (such as SAR polarimetric, visible light multispectral, thermal infrared) are organized into a four-dimensional matrix to form multi-period and multi-band remote sensing raster data.

[0024] Step S2 acquires a joint sample vector and extracts the anomaly coverage area based on a four-dimensional probability distribution. This step uses a buoy array and a high-resolution numerical model to obtain external environmental data for the target monitoring sea area, including wind speed, wind direction, wave height, and shortwave spectrum peak energy. For example, during a certain monitoring period, an abnormal increase in wave height and shortwave spectrum energy obtained by multiple buoys indicates abnormally active surface dynamics in the sea area. This environmental data is interpolated onto a grid system corresponding to the remote sensing grid at the sea area grid point, forming a joint sample vector at the same spatial location but at different times.

[0025] Multivariate Gaussian kernel density estimation is used and the bandwidth matrix is ​​determined based on the Silverman rule to establish a four-dimensional joint probability density distribution function, thereby obtaining the cumulative distribution function value CDF, which is further converted into the anomaly sensitivity value S=1-CDF.

[0026] The scattering anomaly confidence value C is obtained by coupling the S value with the local SAR scattering anomaly normalization value N, which helps to identify grid cells with significant anomaly deviations in multi-period and multi-band data and form a preliminary anomaly coverage area.

[0027] S3 is based on enhanced abnormal unit detection and dynamic correction in the abnormal coverage area. On the basis of the abnormal coverage area, it further introduces inverse sensitivity weights to perform residual analysis of the light and heat dual domains (such as water color index, eutrophication index and thermal infrared brightness temperature anomalies) on weakly significant units other than the preliminary abnormal units. For example, in a certain monitoring, a decrease in water color index, an increase in eutrophication index and a significantly low brightness temperature all indicate potential oil film.

[0028] By combining the weighted least squares fitting model with the inverse sensitivity weight, the optical and thermal residual values ​​are calculated, and the enhanced abnormal units are determined based on the local dynamic quantile threshold (such as P95) to form a more sensitive abnormal layer. This not only avoids omissions caused by the lack of visibility of single SAR scattering, but also significantly improves the sensitivity of detection for thin or early oil films under high wind speed surge conditions.

[0029] In step S4, the abnormal connected regions formed by abnormal units are extracted through connectivity analysis, and their area, perimeter, and compactness are calculated (for example, when the compactness is less than 0.2, it indicates that the regional morphology is highly irregular). This can effectively identify abnormal regions with unstable morphology and prone to drastic boundary fluctuations, and determine them as sensitive areas for abnormal evolution.

[0030] In subsequent monitoring cycles, by comparing the spatial consistency deviation value of anomaly detection in the sensitive area (for example, the coverage deviation exceeds 0.25), the anomaly sensitivity S of the environmental conditions corresponding to the area is automatically triggered to be re-estimated, thereby dynamically adjusting the calculation link of the inverse sensitivity weight and the scattering anomaly confidence value to achieve local adaptive optimization.

[0031] In summary, through the gradual progression from S1 to S4, the present invention has formed a full-process technology chain from the spatiotemporal integration of multi-period and multi-band data, anomaly probability mapping driven by wind fields, enhanced detection of photothermal domains, and anomaly morphology feedback closed loop, which has improved the accuracy and stability of marine oil film, sewage discharge and eutrophication anomaly detection under different wind and wave conditions, and effectively avoided the misjudgment or missed judgment problems of traditional single-period SAR detection under high wind speed surge masking and low wind speed mirror confusion.

[0032] Example 2: Please refer to Figure 1 , specifically: S1 includes: S11: Use orbit prediction tools (e.g., professional orbit prediction software STK (Systems Tool Kit), GMAT (General Mission Analysis Tool)) to obtain all satellite resources used for sea area monitoring within the scheduled monitoring period and obtain the resource set; S12: Based on the available band capabilities and orbit transit timing of each satellite in the resource set, combined with the geometric coverage relationship of the entire sea area, greedy scheduling is performed through the time-sorted queue to generate a multi-satellite multi-mode schedule list; When developing a monitoring schedule, it is necessary to comprehensively consider which bands each satellite can observe (available band capacity), when it will pass over a specific place (transit sequence), and which part of the ocean it can see from orbit (geometric coverage). This will enable a reasonable multi-band, multi-period, and multi-regional observation plan to be arranged. The available band capability refers to the wavelengths at which it can obtain reflection or scattering information. For example, Sentinel-2 cannot obtain SAR, and Sentinel-1 cannot obtain thermal infrared. Greedy scheduling always chooses the optimal arrangement under current conditions without considering the future global optimum. In other words, it only pursues the current local optimum. Unlike genetic scheduling, it does not try many combinations to find the global optimum. Instead, it prioritizes the earliest opportunity to complete the task and quickly meet demand. For example, if Sentinel-1 can take SAR images on the morning of the first day, it will be scheduled immediately. If Sentinel-2 can take spectral images on the afternoon of the first day, it will be scheduled immediately.

[0033] S13: Using a multi-satellite multi-mode schedule, acquire SAR images, visible multispectral images, and thermal infrared images from different satellites at different times to generate an image library. Geometric registration is then performed on the images in the library. Specifically, the images in the library are unified into the same map projection coordinate system (e.g., WGS84 / UTM Zone) using an image geometric registration tool, and resampled to a uniform spatial resolution (e.g., 10m or 30m grid). The size of the sea area grid is determined based on the spatial resolution. Spatial resolution refers to the area represented by each pixel on the ground. The set spatial resolution of 10m or 30m determines the geographical range of each grid point, that is, the actual area represented by each sea grid point.

[0034] Resampling to a uniform spatial resolution means adjusting the spatial resolution of the image to a uniform scale to ensure that the actual ground coverage of each sea area grid point is consistent.

[0035] SAR images are used for monitoring dark spot anomalies, visible light multispectral images are used for monitoring turbidity and algal bloom indicators, and thermal infrared images are used for detecting temperature and sewage anomalies. S14: After completing the geometric registration, the pixel values ​​of the same location at different times and in different bands are organized into a four-dimensional matrix to form multi-period multi-band remote sensing raster data for subsequent arbitrary slice queries (for example, to view the changes in the same location over time, or the responses of different bands at the same time); Among them, the time sorting queue is a time window formed by sorting the satellites from early to late according to their orbital transit timing; The multi-satellite multi-mode scheduling list clearly stipulates which satellite will acquire which type of image data at different time nodes, thus providing a scheduling basis for subsequent time-series multi-band acquisition.

[0036] The four-dimensional matrix includes spatial dimension, time dimension and band dimension, wherein the spatial dimension includes longitude and latitude; Multi-period and multi-band remote sensing raster data is a pixel matrix data containing multiple bands obtained by remote sensing sensors at multiple different time points in the same monitoring area. Each pixel records the observation value of its corresponding geographical location in each band and is accompanied by a time label. It can be directly used for time series analysis, spectral comparison, or machine learning.

[0037] In this embodiment, through S11, using orbital prediction tools (such as GMAT, STK or self-developed satellite orbit simulation scheduling engine), the trajectories, transit times and visible windows of all satellite resources that can be used for sea area monitoring can be calculated in advance within the predetermined monitoring period, and a resource set can be quickly formed. This resource set generation based on orbital prediction avoids the inefficient method of manual satellite-by-satellite retrieval and static scheduling, and provides the necessary time series and available band capability input for subsequent scheduling.

[0038] In S12, based on the band capabilities of each satellite in the resource pool (for example, Sentinel-1 has C-band SAR, Landsat-8 has multi-spectral visible light and thermal infrared), combined with their orbital transit timing and the geometric coverage relationship of the entire ocean area, a multi-satellite multi-mode schedule list is quickly generated by constructing a time-sorted queue and executing greedy scheduling.

[0039] Greedy scheduling prioritizes satellite missions that cover the largest ocean area or acquire missing bands within each time window, thereby maximizing coverage of multiple spatial and temporal bands within limited observation opportunities. For example, if Sentinel-1, GF-3, and Landsat-8 pass through in sequence during a monitoring period, the system can prioritize Sentinel-1 for SAR, GF-3 for multi-polarization, and Landsat-8 for thermal infrared and visible light. This avoids band and time redundancy and enables efficient organization of cross-satellite resources.

[0040] In S13, SAR, visible multispectral, and thermal infrared images acquired at different times and by different satellites are acquired sequentially, following the scheduling sequence determined by the multi-satellite, multi-mode schedule list, to form an image library. Subsequently, image geometric registration tools (such as control point-based affine registration or automatic block matching algorithms) are used to unify all images in the library to the same map projection coordinate system (such as WGS84 / UTM Zone) and resample them to a uniform spatial resolution (such as a 10m grid). This step standardizes the spatial resolution and defines the grid size of the ocean area, enabling seamless overlay of satellite images from different sources and times within the same geographic reference frame. For example, a 30m-resolution Landsat-8 thermal infrared image and a 10m-resolution Sentinel-2 multispectral image can be unified to 10m resolution through bilinear interpolation, laying the geometric and resolution foundation for subsequent pixel-level multi-time series and cross-band joint analysis.

[0041] Through S14, the geometrically registered and resampled images are organized into a four-dimensional matrix based on spatial dimensions (longitude and latitude), temporal dimensions (acquisition time), and band dimensions (SAR polarimetric, visible multispectral, and thermal infrared), forming multi-period and multi-band remote sensing raster data. This allows for subsequent direct slice queries against any time, band, or location, facilitating rapid time series analysis, spectral combination comparison, or machine learning training. For example, one can directly retrieve "all SAR VH polarimetric backscatter sequences for the same sea area coordinate point over the past three months" or "spectral feature vectors in different bands on the same day," significantly reducing the complexity of subsequent cross-time and cross-band analysis and greatly improving the processing efficiency of multi-stage algorithms such as anomaly detection, photothermal residual analysis, and ecological prediction.

[0042] In summary, S1 and its substeps S11 to S14 have formed an efficient, automated, multi-period, multi-band, multi-satellite remote sensing data acquisition and standardization process in the entire satellite remote sensing-based sea area monitoring program, which not only improves the utilization efficiency of multi-source data, but also provides stable and consistent spatiotemporal data input for subsequent joint probability density, scattering anomaly confidence detection, photothermal residual analysis and dynamic threshold feedback, effectively supporting the precise detection and tracking of sea oil spills, illegal discharge of pollutants and eutrophication anomalies in complex environments.

[0043] SAR images are data products similar to black and white grayscale images generated by using synthetic aperture radar to transmit microwaves from satellites or aircraft, then receiving the signals reflected from the ground.

[0044] Example 3: Please refer to Figure 1 , specifically: S2 includes: S21: Obtain external environmental data within the target monitoring area through the buoy array, including wind speed, wind direction, wave height, and shortwave spectrum peak energy. Perform spatial interpolation on the external environmental data according to the sea area grid points to align it with the multi-period multi-band remote sensing grid data. Simultaneously, obtain a multi-period external environmental sample set, wherein the multi-period external environmental sample set includes joint sample vectors at different times. Spatial interpolation refers to the process of inferring the values ​​of other unknown data points based on the values ​​of known data points. Common interpolation methods include nearest neighbor interpolation, bilinear interpolation, or Kriging interpolation. Wind speed refers to the speed of air movement, typically measured at 10 meters above the ground. It reflects wind intensity and significantly influences sea surface roughness, wave formation, and variability. In SAR remote sensing, high wind speeds create a rough sea surface, resulting in stronger SAR backscatter and brighter images. Low wind speeds create a calmer sea surface, resulting in weaker SAR backscatter and darker images. Wind speed can be obtained from numerical meteorological models, such as ERA5 or GFS, which provide global and local wind speed forecasts, or from actual wind speed monitoring by in-situ buoys or radar. Wind direction refers to the angle of the wind, that is, where the wind is blowing from, usually relative to true north. Combined with wind speed, wind direction can affect the direction and height of ocean surface waves. In SAR imagery, wind direction can change the direction of microwave scattering from the sea surface, affecting the wind field's appearance in SAR images. Wind direction deviations can cause variations in wave patterns in different areas, affecting the reflection characteristics of SAR images. Wind direction can be measured by buoys: buoys on the sea surface measure wind direction using sensors, or wind direction data can be obtained through weather radar. Wave height represents the average wave height in the ocean, that is, the height of the largest 1 / 3 of the waves. It is a commonly used indicator to describe sea surface fluctuations. Wave height has a significant impact on sea surface reflection. High waves: large waves, rougher sea surface, strong SAR echo reflection, and brighter image. Low waves: relatively smooth sea surface, weak SAR echo reflection, and darker image. Wave height data is directly obtained by measuring sea surface fluctuations through buoys deployed on the sea surface. Shortwave spectrum peak energy is the peak energy in the shortwave band of the wave energy spectrum (often used to describe the response of short-period waves to microwave scattering). It reflects the wavelength and frequency of the waves and directly affects the scattering characteristics of microwaves and SAR. High shortwave energy usually indicates strong fluctuations on the sea surface, which has a strong impact on SAR backscattering. Low shortwave energy indicates a calm sea surface with weak scattering. Shortwave spectrum peak energy can be used to collect field data through buoy arrays. The joint sample vector is the state where wind speed, wind direction, wave height, and shortwave spectrum energy co-occur in the same space and time; S22: Based on a multi-period external environment sample set, a multivariate kernel density estimation algorithm is used to perform kernel density estimation on the joint sample vector to obtain a four-dimensional joint probability density distribution function. The four-dimensional joint probability density distribution function is used to quantify the probability density value of any joint sample vector (i.e., environmental combination) that appears simultaneously in historical samples.

[0045] The specific expression of the four-dimensional joint probability density distribution function is: ; in, is a four-dimensional joint probability density distribution function, which indicates the probability density value of the joint sample vector (i.e., environmental combination) appearing simultaneously in the historical samples; U is the wind speed, D is the wind direction, Hs is the wave height, Spk is the shortwave spectrum peak energy, m is the total number of historical samples, i is the historical sample index, and H is a four-dimensional bandwidth matrix (4×4 positive definite matrix), which describes the smoothness and mutual relationship between the four variables. is the determinant of H, which is used to normalize the multidimensional volume. Represents the inverse of the square root of the matrix H, which is used to normalize the original difference vector to the scale required by the kernel function. is the wind speed of the i-th sample, is the wind direction of the i-th sample, is the wave height of the i-th sample, is the shortwave spectrum peak energy of the i-th sample, is the current environment combination point (i.e., the current joint sample vector) for which the probability density is to be estimated, It is a four-dimensional kernel function. Here we generally use the four-dimensional Gaussian kernel function (Gaussian distribution function), which is often written as: , where, because it is four-dimensional, the exponential term is the square norm, is the kernel function value. In kernel density estimation, It is used to measure the local density contribution at a given normalized distance z. z represents the normalized offset vector after scaling by the bandwidth matrix H. It scales the difference between the four variables between each sample point and the target point to a uniform kernel scale. It is a four-dimensional column vector. Represents the square norm of the vector, that is, the square of the Euclidean distance of the z vector. It represents the square of the distance from the center in the Gaussian distribution exponential term and is used to control density attenuation. e is the natural base, and its value is approximately 2.71828; It is a normalization factor used to ensure that the integral of the kernel density function in the entire four-dimensional space is always 1; It means that the difference vector is scaled by H, which is equivalent to converting the original four-dimensional coordinates to the kernel function coordinate system with a variance of H. is the contribution of the i-th sample to the current environment combination; The multivariate kernel density estimation adopts a Gaussian kernel function and determines the bandwidth matrix based on the Silverman rule; The multivariate kernel density estimation algorithm is used to simultaneously estimate the joint probability density of multidimensional random variables.

[0046] S2 also includes: S23: According to the four-dimensional joint probability density distribution function, the sum of all historical probability density values ​​lower than the current joint sample vector (i.e., the environment combination) is counted to obtain a cumulative distribution function value. The cumulative distribution function value is used to quantify the cumulative occurrence probability of the current joint sample vector in the historical samples; The cumulative distribution function value is integrated under the corresponding four-dimensional density to obtain the cumulative probability of the current combination appearing, that is, the probability of the combination appearing in the past is less than or equal to the current combination. What is the cumulative proportion of samples? S24: Using the cumulative distribution function value, obtain the abnormal sensitivity value based on wind field drive corresponding to the current monitoring time and spatial position, specifically: , where S is the abnormal sensitivity value and CDF is the cumulative distribution function value; the abnormal sensitivity value is used to reflect the rarity of the currently monitored joint sample vector (i.e., environmental combination) at the corresponding geographical location; If the CDF value is large (close to 1), it means that the current environmental combination appears frequently in historical data and is normal. If the CDF value is small (close to 0), it means that it is relatively rare. Therefore, the abnormal sensitivity value indicates the degree to which the current environmental combination is rare in history. S25: Determine the SAR backscatter coefficient corresponding to each pixel point in the SAR image; In synthetic aperture radar (SAR) images, the value of each pixel is actually the backscatter coefficient (usually expressed in dB), which represents the power intensity of the radar wave reflected back to the radar receiver after hitting the sea surface. It is affected by sea surface roughness, wind speed, oil film, wave height, etc.

[0047] S26: At each sea area grid point, the local mean and standard deviation of the SAR backscatter coefficient within a sliding window centered at each sea area grid point are calculated, and the backscatter coefficient of each sea area grid point is normalized to obtain a scattering normalized outlier value; Step S26 is equivalent to doing a Z-score, which indicates how many standard deviations the backscatter of the corresponding point is higher than the surrounding average. This can better detect local anomalies. If the oil film causes a significant decrease in scattering, the scattering normalized anomaly value will be significantly lower than 0, otherwise, it will be significantly higher than 0. When performing local statistics (such as local mean and standard deviation), a small rectangular area (such as 3×3, 5×5, or 7×7 grid points) is formed with a certain sea area grid point as the center, and the statistical values ​​are collected in this area.

[0048] The scatter normalized outlier is used to measure the degree of deviation between the corresponding sea area grid point and the local normal water surface; S27: Couple the anomaly sensitivity value with the scattering normalized anomaly value to calculate the scattering anomaly confidence value at each sea area grid point, specifically: ,in, is the confidence value of scattering anomaly, is the abnormal sensitivity value, is the scattering normalized anomaly value, which is the degree of anomaly of the SAR backscatter coefficient of the corresponding sea area grid point relative to the local area; The scattering anomaly confidence value is used to form a scattering anomaly confidence surface in the monitoring area, providing a basis for the subsequent generation of anomaly detection thresholds based on local dynamic quantiles and the extraction of suspected anomaly masks; S28: Construct a scattering anomaly confidence surface based on the scattering anomaly confidence values ​​at each sea area grid point; The scattering anomaly confidence surface refers to a confidence map of anomaly detection results calculated at each grid point in the sea area by combining external environmental driving factors such as wind field (wind speed and direction) and waves (wave height and shortwave spectrum peak energy) with the backscattering coefficient of the SAR image. This map is used to characterize scattering anomaly areas on the sea surface and can identify anomalies such as oil spills or floating ice.

[0049] S29: Generate a dynamic local threshold surface based on the adaptive quantile method. By comparing the scattering anomaly confidence surface with the dynamic local threshold surface, an anomaly mask is formed. Specifically: ,in, For the location The abnormal mask on For the location The confidence value of the scattering anomaly on For the location Dynamic local threshold on Indicates the position coordinates in the monitoring space; The adaptive quantile method analysis is as follows: statistics of all positions The quantile values ​​of the scattering anomaly confidence value on the scattering anomaly (such as P90, P95, P98) are then adaptively set for different sea environments: ,in, Indicates the location The quantile adjustment coefficient is dynamically set according to environmental characteristics such as local wind speed, wind direction, wave height, and shortwave spectrum energy. It usually fluctuates within a certain range (such as 0.85 to 1.15) to adapt to sensitivity adjustment under different sea conditions. It represents the 95th percentile confidence value of the scattering anomaly obtained in the target monitoring area, which is used to provide the global reference percentile benchmark value for the entire domain.

[0050] S291: The sea area grid points with an anomaly mask of 1 are regarded as anomaly units, and the anomaly coverage area is obtained through statistics.

[0051] The cumulative distribution function value CDF is achieved by integrating the multivariate kernel density function in four-dimensional space and is used to quantify the rarity of the current environmental conditions in the monitoring area.

[0052] In this embodiment, S21 utilizes a numerical model (e.g., a wind and wave field model based on coupled atmosphere-ocean dynamics) and a buoy array deployed in the target sea area to obtain external environmental data at different time points, including wind speed, wind direction, wave height, and shortwave spectrum peak energy. For example, when monitoring the waters near the Zhoushan Islands, the wind speed field at a height of 10 meters can be derived in real time from the numerical model, and the corresponding measured wave height can be obtained from the buoy data.

[0053] The shortwave spectrum peak energy is a physical quantity used to characterize the energy of the shortwave band of the sea surface, which directly affects the imaging sensitivity of SAR to fine oil films; In step S21, these data are spatially interpolated according to the sea area grid points (i.e., the grid positions fixed after geometric registration) to ensure that different data sources are strictly aligned spatially with the multi-period and multi-band remote sensing raster data, laying a unified foundation for subsequent pixel-level multi-source coupling.

[0054] Based on the multi-period external environment sample set generated in step S21, a multivariate kernel density estimation algorithm (using a Gaussian kernel function and adaptively calculating the bandwidth matrix based on the Silverman rule) statistically models the joint sample vector—that is, the combined state of wind speed, wind direction, wave height, and shortwave spectral peak energy in the same space and time—to obtain a four-dimensional joint probability density distribution function. For example, if the frequency of wind speed of 6 m / s, wind direction of 120°, wave height of 1.5 m, and shortwave spectral peak energy of 50 is low in history, the corresponding joint probability density value will be small. This distribution function is used to comprehensively quantify the probability density of any marine environment combination in historical statistics and is an important basis for subsequent rarity detection.

[0055] At S23, the four-dimensional joint probability density function is integrated over the four-dimensional space, and the proportion of all historical cases where the probability density value is lower than the current joint sample vector is counted to obtain the cumulative distribution function (CDF). For example, when CDF = 0.8, it means that the current ocean environment is more common than 80% of the cases in the historical samples.

[0056] S24 uses the formula S = 1 − CDF to define anomaly sensitivity, S, to quantify the rarity of the current environment. An S = 0.2 indicates that the environment is relatively common, while an S = 0.95 indicates that the environment is extremely rare in history, worthy of higher vigilance for any remote sensing anomalies detected.

[0057] In S25, the SAR backscatter coefficient is extracted for each pixel based on the SAR image to characterize the radar scattering characteristics of that point. A lower backscatter coefficient often indicates a smoother sea surface, which may be caused by oil film, sewage discharge, etc.

[0058] At each sea area grid point, S26 calculates the local mean and standard deviation within a sliding window centered at that point, achieving local normalization. This results in a scatter normalized outlier value, N, which characterizes the degree of deviation of that pixel from the surrounding normal water surface. For example, if N = 2, this means that the pixel is twice the standard deviation higher than its surroundings.

[0059] In S27, the anomaly sensitivity value S is coupled with the scattering normalized anomaly value N to calculate the scattering anomaly confidence value. If a point is located in an extremely rare wind and wave combination (S is high) and its scattering anomaly is significant (N is high), the scattering anomaly confidence value C will be higher, indicating that the confidence in the anomaly is very high.

[0060] In S28, a scattering anomaly confidence surface is constructed across the entire monitoring area based on the scattering anomaly confidence value C. In S29, a dynamic threshold surface is generated in the sliding window using local dynamic quantiles (e.g., the 95th percentile value) and compared with the scattering anomaly confidence surface C(x,y) to form an anomaly mask. In S291, the sea area grid points marked as 1 in the mask are further treated as anomaly units, and statistics are used to form an anomaly coverage area.

[0061] In summary, S21 spatially aligns and temporally accumulates external environmental samples; S22-S24 quantify the rarity of environmental combinations in probability space and converts it into sensitivity S; S25-S27 normalize backscatter in physical scattering space and couple it with environmental rarity to form confidence; S28-S291 dynamically generate thresholds in the local statistical domain, extract anomaly units, and form anomaly coverage areas. This not only improves the adaptability to different environmental conditions and avoids the inherent defects of the single-threshold strategy, but also effectively identifies oil film, sewage discharge, and water eutrophication anomalies that may be hidden in high or low wind speeds through wind-wave drive and radar scattering coupling, thereby improving the overall detection accuracy and stability.

[0062] Example 4: Please refer to Figure 1 Specifically: S3 includes: S31: Based on the scattering anomaly confidence value, calculate the inverse sensitivity weight of each sea area grid point in the target monitoring area; The inverse sensitivity weight is obtained as follows: A=1-C, where A is the inverse sensitivity weight, which is used to form spatial attention guidance for SAR low-confidence units; S32: Determine weakly significant units based on the abnormal units, and extract the abnormal values ​​of the sea water color index, eutrophication index and thermal infrared brightness temperature within the weakly significant units; By extracting the water color index, eutrophication index and thermal infrared brightness temperature anomaly values ​​within the weakly significant units, a dual-domain optical and thermal supplement to the SAR detection blind spots is formed, thereby improving the spatial coverage and sensitivity of the overall sea area anomaly detection.

[0063] Weakly significant units refer to units other than abnormal units; The sea color index is calculated based on the ratio of the near infrared to red bands of the multispectral pixels. The specific calculation method is: ,in, is the sea color index, It is the near-infrared band reflectance pixel value at the monitoring sea area grid point, which is usually used to capture the near-infrared absorption and reflection characteristics of algae, phytoplankton and suspended particles in the water body; The red band reflectance pixel value at the monitoring sea area grid point is used for sensitivity analysis of water bodies, chlorophyll, etc. to red light absorption; The near-infrared band reflectance pixel value refers to the value of the near-infrared band (usually )'s photosensitive element and is radiation-corrected for the apparent surface reflectance.

[0064] The red band reflectance pixel value refers to the surface reflectance obtained by the satellite sensor in the red light band (usually 0.63–0.69 μm) at the grid point of the monitored sea area.

[0065] In seawater, pure water reflects almost no near-infrared light (absorbs very strongly) and strongly absorbs red light, resulting in a low water color index. However, the presence of large amounts of algae and suspended particles alters the relative reflectance of near-infrared and red light, causing the water color index to rise. Therefore, the water color index is used to distinguish abnormal water bodies (algal blooms, muddy water) from normal water bodies in the spectral dimension.

[0066] The eutrophication index is calculated based on the ratio of the green and red bands of the multispectral pixels. The specific calculation method is: ,in, is the eutrophication index, It is the green band reflectance pixel value at the grid point of the monitored sea area, which is often used to detect eutrophication indicators such as phytoplankton chlorophyll, because green light can reflect the scattering information of water pigments and algae; The green band reflectance pixel value refers to the surface reflectance obtained by the satellite sensor in the green light band (usually 0.52–0.60 μm) at the grid point of the monitored sea area.

[0067] Generally, green light can penetrate the upper water layer better and be scattered by phytoplankton chlorophyll, while red light is more strongly absorbed by algae chlorophyll. Therefore, the eutrophication index can significantly increase the sensitivity to eutrophication of water bodies (massive algae reproduction) and is used to detect eutrophication risk areas caused by early algae reproduction.

[0068] The thermal infrared brightness temperature anomaly is the deviation between the thermal infrared brightness temperature value of the corresponding sea grid point and the mean of its neighborhood sliding window. The specific calculation method is: ,in, is the abnormal value of thermal infrared brightness temperature, It is the thermal infrared brightness temperature (surface radiation temperature) at the grid points in the monitored sea area, which is used to reflect the spatial distribution of sea surface temperature; is the local mean of the thermal infrared brightness temperature within a sliding window centered at the grid point, which is used to calculate the temperature deviation relative to the neighborhood background; If there is an oil film on the surface of the water body, it will block evaporative cooling, resulting in slightly higher local temperatures. Or if there is abnormal sewage discharge, the temperature of the discharged wastewater is different, which will also cause local brightness temperature anomalies. Therefore, the thermal infrared brightness temperature anomaly value is used to detect potential abnormal areas with temperature disturbances compared with the surrounding normal waters. S33: Using the inverse sensitivity weight as the weighting factor, the local multivariate regression model is fitted by the weighted least squares method to establish a dual-domain photothermal prediction model based on the sea water color index, eutrophication index and thermal infrared brightness temperature anomaly values ​​to calculate the photothermal residual value on each weakly significant unit.

[0069] Fit a local multiple regression model using the inverse sensitivity weight as the weighting coefficient in the entire domain: ,in, It is a dual-domain light and heat prediction model, which is used to output regression prediction values, indicating the expected observed light and heat values ​​after multivariate fitting of light and heat characteristics at the sea grid point. is the intercept term of the regression model; 、 、 are all regression coefficients, indicating the influence of NDVI, FUI, and ΔTIR on the predicted value; they are obtained by fitting the global data using the weighted least squares method; The inverse sensitivity weighted least squares is used to make the SAR more sensitive to the photothermal residual in the low confidence region.

[0070] Calculate the photothermal residual value: ;in, is the photothermal residual value, For the location The right to reverse sensitivity, It is the actual light and thermal index value obtained directly from remote sensing images (or multi-band combinations), representing real observations and directly representing the actual spectral-thermal infrared joint anomaly extracted from satellite image pixels; It is a dual-domain light and heat prediction model, which is used to output regression prediction values, representing the expected observed light and heat values ​​after multivariate fitting of light and heat characteristics at the sea grid points; In this way, photothermal anomalies will be amplified in areas where SAR is not significant.

[0071] S3 also includes: S34: At each sea area grid point, based on the light and thermal residual value, the 95th percentile value of the light and thermal residual value in the sliding window centered on the corresponding sea area grid point is used as the local dynamic threshold. If the light and thermal residual value exceeds the local dynamic threshold, the corresponding sea area grid point is marked as an enhanced anomaly unit; S35: By enhancing the anomaly cells, the anomaly coverage area is dynamically corrected to obtain the final anomaly detection mask within the target monitoring area, thereby significantly improving the detection rate of hidden oil films or weak pollution discharges under high wind speed and surge conditions.

[0072] The final anomaly detection mask is based on the abnormal unit, and the enhanced abnormal unit is also marked as 1; Enhanced abnormal units are weakly significant units; In this embodiment, after obtaining the scattering anomaly confidence value in step S31, the present invention calculates an inverse sensitivity weight A = 1 - C for each sea area grid point, which is used to guide subsequent attention to low-confidence cells. For example, if the scattering anomaly confidence value C at some grid points has reached 0.9, indicating that the anomaly is highly reliable, the inverse sensitivity weight A = 0.1, and the subsequent analysis weight is reduced. For grid points with C = 0.2, the inverse sensitivity weight A = 0.8 is used to emphasize the enhanced detection sensitivity of these weakly significant areas in the photothermal analysis, thereby effectively focusing on potential areas that are easily missed.

[0073] Through step S32, the sea water color index is extracted from the weakly significant units except the abnormal units, and is calculated by the ratio of the near-infrared and red bands of the multispectral pixels. For example, when the sea water color index value increases significantly, it may indicate a large-scale reproduction of plankton in the water body; the eutrophication index is calculated by the ratio of the green and red bands. For example, when the green band is enhanced, it indicates the occurrence of algal bloom; the thermal infrared brightness temperature anomaly is the deviation between the thermal infrared brightness temperature value and the neighborhood sliding mean. For example, water temperature anomaly may be associated with pollution emissions. These optical and thermal multi-domain parameters can reveal surface organisms or temperature perturbations that are difficult to detect with SAR.

[0074] In step S33, the inverse sensitivity weight is used as a weighting factor to fit a local multivariate regression model using the sea color index, eutrophication index, and thermal infrared brightness temperature anomaly values ​​as inputs. This model is used to predict the photothermal response of weakly significant units under normal conditions. For example, if the photothermal fitting residual increases abnormally at a location with a high inverse sensitivity weight (e.g., A = 0.85), it means that even if the SAR confidence is low, there may be an anomaly in the photothermal coupling in that area, thus compensating for the limitations of SAR in detecting thin oil films or early pollution discharges.

[0075] In step S34, the 95th percentile of a sliding window centered at each sea area grid point is used as a local dynamic threshold based on the residual values ​​of light and heat. For example, in a eutrophic water area, if most residuals are concentrated between 1.0 and 1.5, and the 95th percentile reaches 1.7, a grid point with a residual value of 2.2 will be marked as an enhanced anomaly unit, avoiding the problem of local feature fluctuations being masked by the use of a global fixed threshold.

[0076] In step S35, the enhanced anomaly cells are superimposed onto the original anomaly coverage area to form the final anomaly detection mask. This means that, based on the original anomaly cells, these areas with significant optical and thermal residual anomalies are also marked as 1. For example, in high wind speeds (>10m / s), the original SAR anomalies are difficult to detect. However, this method can still accurately identify weak sewage discharge areas or oil slicks in areas of eutrophication or temperature anomalies, significantly improving detection rates and monitoring stability in complex sea conditions.

[0077] In summary, the present invention uses inverse sensitivity weights to adaptively guide attention in space, integrates multi-dimensional information of water color index, eutrophication and thermal infrared residual to construct a light-thermal dual-domain prediction, and performs local threshold discrimination based on the dynamic quantile method. It effectively overcomes the detection blind spots of a single SAR or a global fixed threshold in high waves and complex waters, and improves the detection rate and stability of abnormal phenomena such as weak oil spills, minor sewage discharges and early eutrophication.

[0078] Example 5: Please refer to Figure 1 , specifically: S4 includes: S41: Obtain multiple abnormal connected areas according to the final anomaly detection mask in the target monitoring area; S42: extracting the area (i.e., number of pixels) and perimeter (i.e., boundary length) of each abnormal connected region and calculating its compactness. The compactness is used to locate abnormal unit regions with unstable morphology and is an important geometric indicator for detecting abnormal morphological stability. The firmness is obtained by the following formula: ,in, For firmness, is the area, is the perimeter, is pi; An abnormally connected region refers to a region where all adjacent regions belong to the anomaly detection mask; S43: When the compactness is lower than a preset compactness threshold, the corresponding abnormally connected area is marked as an abnormal evolution sensitive area, otherwise no marking is performed.

[0079] The S4 also includes: S44: After monitoring, obtain the abnormal evolution sensitive area in the next monitoring cycle, compare the deviation values ​​of the abnormal units in the abnormal evolution sensitive area during each monitoring cycle, and if the deviation value exceeds the preset deviation threshold, automatically repeat S2 to correct the abnormal sensitivity value driven by the wind field and dynamically adjust the scattering anomaly confidence value and inverse sensitivity weight.

[0080] The deviation value is the difference in the distribution of abnormal units between the current cycle detection result and the historical average detection result in the abnormal evolution sensitive area, which can be expressed as the change in spatial coverage; In this embodiment, in step S41, multiple abnormal connected regions (closed regions where all adjacent pixels are abnormal units) are extracted based on the final anomaly detection mask within the target monitoring area (i.e., the anomaly marker layer after integrating SAR scattering and photothermal enhancement). For example, if a single monitoring event produces several dark oil slicks curving along the coast and scattered areas of foam disturbance on the sea surface, these different abnormal patches can be identified through regional connectivity.

[0081] In S42, the area (number of pixels) and perimeter (boundary length) of each anomalously connected region are further extracted, and the compactness is calculated. Compactness is a key indicator of the geometric compactness of the region and is used to quantify the shape regularity of the anomalous patch. For example, a nearly circular oil film patch, with a large area and a relatively short perimeter, has a compactness close to 1; whereas a long, slender, serpentine drainage strip, with a smaller area and a longer perimeter, has a compactness far below 1. This step automatically locates anomalous units with fragmented morphology, complex boundaries, and the tendency to split or shrink over multiple periods.

[0082] In S43, if the compactness of an abnormally connected region is detected to be below a preset compactness threshold (e.g., 0.5), the region is marked as a sensitive area for abnormal evolution. For example, if a narrow, bifurcated pollution belt is repeatedly detected near a port and its compactness is below the threshold, the area will be automatically marked as a sensitive area for tracking its future morphological evolution. Conversely, large abnormal patches with rounded shapes and stable boundaries will not be included in the sensitive area marking to avoid over-adjusting naturally stable abnormal areas.

[0083] In S44, after entering the next monitoring cycle, the distribution of abnormal units in the abnormal evolution sensitive area in the current cycle is compared with the historical cycle average. The deviation value is the difference between the spatial coverage rate of the abnormal units in the current cycle and the historical average (for example, if the historical average coverage rate is 30% and it suddenly increases to 50% in the current cycle, the deviation value is 20%). When the deviation value exceeds the preset deviation threshold (such as 15%), the S2 process is automatically triggered: the four-dimensional joint probability density distribution function is re-updated based on multi-period environmental samples, the abnormal sensitivity value is dynamically corrected, and the scattering anomaly confidence value and inverse sensitivity weight are further adjusted. For example, in multi-period monitoring, if the sewage discharge distribution in a narrow and long sensitive area deviates significantly from the historical normal due to tidal changes, the system will retrain the environmental probability density model in the sensitive area, so that the sensitivity value S under the same wind speed, wind direction, wave height, and shortwave spectrum combination is automatically increased, thereby improving the abnormality detection capability in future cycles.

[0084] In summary, through the cascade operation of S41 to S44, the present invention can not only quantify the stability of abnormal units based on geometric morphology within a single cycle, but also automatically feedback and correct the detection model through the deviation of the sensitive area of ​​abnormal evolution in multiple cycles, forming a closed loop of local dynamic optimization. In this way, in long-term sea area monitoring tasks, especially in sea areas affected by tides, local wind fields or periodic human pollution discharge, the frequent jitter of detection boundaries and the abnormality false detection rate can be significantly reduced, thereby improving the spatial consistency of the entire monitoring system in time series and the stability of long-term anomaly detection.

[0085] Abnormally connected regions are used to identify spatially continuous abnormal unit patches.

[0086] Compactness is used to quantify the geometric complexity of abnormal units and detect whether the shape is changeable.

[0087] The key dynamic tracking areas of the abnormal evolution sensitive areas are screened by compactness and used for adaptive adjustment of subsequent multi-cycle monitoring.

[0088] The deviation value is used to measure the significant difference in anomaly distribution in the anomaly evolution sensitive area between different cycles and to determine whether the detection model needs to be locally updated.

[0089] Example 6: Please refer to Figure 4 ,Specifically: A sea area monitoring system based on satellite remote sensing technology, including, The data module is used to generate multi-period and multi-band remote sensing raster data by real-time monitoring of resource sets in the sea area; The preliminary identification module is used to obtain joint sample vectors at different times by monitoring the external environmental data within the target monitoring area. The multivariate kernel density estimation algorithm is used to obtain the four-dimensional joint probability density distribution function. Combined with multi-period and multi-band remote sensing raster data, the abnormal coverage area is obtained. The secondary identification module is used to identify and enhance abnormal units based on abnormal coverage areas and perform dynamic correction operations; The local optimization module is used to obtain abnormal connected areas after dynamic correction, analyze abnormal units with unstable morphology, determine the abnormal evolution sensitive areas, and implement local adaptive optimization methods based on the changes in the abnormal evolution sensitive areas.

[0090] While embodiments of the present invention have been shown and described, it will be appreciated by those skilled in the art that various changes, modifications, substitutions, and variations may be made to these embodiments without departing from the principles and spirit of the invention, and that the scope of the invention is defined by the appended claims and their equivalents.

Claims

1. A method for sea area monitoring based on satellite remote sensing technology, characterized by: The following steps are included: S1: Generate multi-period and multi-band remote sensing raster data by real-time monitoring of resource sets in the sea area; S2: By monitoring the external environmental data within the target monitoring area, joint sample vectors at different times are obtained. A multivariate kernel density estimation algorithm is used to obtain a four-dimensional joint probability density distribution function. Combined with multi-period and multi-band remote sensing raster data, the abnormal coverage area is obtained. S3: Based on the abnormal coverage area, identify and enhance abnormal units and perform dynamic correction operations; S4: After dynamic correction, the abnormal connected area is obtained, and the abnormal units with unstable morphology are analyzed to determine the abnormal evolution sensitive area. Based on the changes in the abnormal evolution sensitive area, local adaptive optimization means are implemented.

2. The method for sea area monitoring based on satellite remote sensing technology according to claim 1, characterized in that: S1 includes: S11: Using orbit prediction tools to obtain all satellite resources used for sea area monitoring within a predetermined monitoring period, and obtain a resource set; S12: Based on the available band capabilities and orbit transit timing of each satellite in the resource set, combined with the geometric coverage relationship of the entire sea area, greedy scheduling is performed through the time-sorted queue to generate a multi-satellite multi-mode schedule list; S13: Using a multi-satellite multi-mode scheduling list, SAR images, visible multispectral images, and thermal infrared images acquired from different satellites at different times are acquired to generate an image library. The images in the image library are then geometrically registered. Specifically, the image geometric registration tool is used to unify the images in the image library into the same map projection coordinate system and resample them to a uniform spatial resolution. The size of the sea area grid is determined based on the spatial resolution. S14: After completing the geometric registration, the pixel values ​​of the same location at different times and in different bands are organized into a four-dimensional matrix to form multi-period multi-band remote sensing raster data; Among them, the time sorting queue is a time window formed by sorting the satellites from early to late according to their orbital transit timing.

3. The method for sea area monitoring based on satellite remote sensing technology according to claim 2, characterized in that: S2 includes: S21: Obtain external environmental data within the target monitoring area through the buoy array, including wind speed, wind direction, wave height, and shortwave spectrum peak energy. Perform spatial interpolation on the external environmental data according to the sea area grid points to align it with the multi-period multi-band remote sensing grid data. Simultaneously, obtain a multi-period external environmental sample set, wherein the multi-period external environmental sample set includes joint sample vectors at different times. S22: Based on a multi-period external environment sample set, a multivariate kernel density estimation algorithm is used to perform kernel density estimation on the joint sample vector to obtain a four-dimensional joint probability density distribution function. The four-dimensional joint probability density distribution function is used to quantify the probability density value of any joint sample vector appearing simultaneously in historical samples.

4. The method for sea area monitoring based on satellite remote sensing technology according to claim 3, characterized in that: S2 also includes: S23: According to the four-dimensional joint probability density distribution function, the sum of all historical probability density values ​​lower than the current joint sample vector is counted to obtain a cumulative distribution function value. The cumulative distribution function value is used to quantify the cumulative occurrence probability of the current joint sample vector in the historical samples; S24: using the cumulative distribution function value, obtaining the abnormal sensitivity value based on wind field driving corresponding to the current monitoring time and spatial position, the abnormal sensitivity value is used to reflect the rarity of the currently monitored joint sample vector at the corresponding geographical location; S25: Determine the SAR backscatter coefficient corresponding to each pixel point in the SAR image; S26: At each sea area grid point, the local mean and standard deviation of the SAR backscatter coefficient within a sliding window centered at each sea area grid point are calculated, and the backscatter coefficient of each sea area grid point is normalized to obtain a scattering normalized outlier value; S27: Couple the anomaly sensitivity value with the scattering normalized anomaly value to calculate the scattering anomaly confidence value at each sea area grid point, specifically: ,in, is the confidence value of scattering anomaly, is the abnormal sensitivity value, is the scatter normalized outlier; S28: Construct a scattering anomaly confidence surface based on the scattering anomaly confidence values ​​at each sea area grid point; S29: Generate a dynamic local threshold surface based on the adaptive quantile method. By comparing the scattering anomaly confidence surface with the dynamic local threshold surface, an anomaly mask is formed. Specifically: ,in, For the location The abnormal mask on For the location The confidence value of the scattering anomaly on For the location Dynamic local threshold on ; S291: The sea area grid points with an anomaly mask of 1 are regarded as anomaly units, and the anomaly coverage area is obtained through statistics.

5. The method for sea area monitoring based on satellite remote sensing technology according to claim 4, characterized in that: S3 includes: S31: Based on the scattering anomaly confidence value, calculate the inverse sensitivity weight of each sea area grid point in the target monitoring area; S32: Determine weakly significant units based on the abnormal units, and extract the abnormal values ​​of the sea water color index, eutrophication index and thermal infrared brightness temperature within the weakly significant units; S33: Using the inverse sensitivity weight as the weighting factor, the local multivariate regression model is fitted by the weighted least squares method to establish a dual-domain photothermal prediction model based on the sea water color index, eutrophication index and thermal infrared brightness temperature anomaly values ​​to calculate the photothermal residual value on each weakly significant unit.

6. The method for sea area monitoring based on satellite remote sensing technology according to claim 5, characterized in that: S3 also includes: S34: At each sea area grid point, based on the light and thermal residual value, the 95th percentile value of the light and thermal residual value in the sliding window centered on the corresponding sea area grid point is used as the local dynamic threshold. If the light and thermal residual value exceeds the local dynamic threshold, the corresponding sea area grid point is marked as an enhanced anomaly unit; S35: Dynamically modify the abnormal coverage area by enhancing the abnormal unit to obtain the final abnormality detection mask in the target monitoring area.

7. The method for sea area monitoring based on satellite remote sensing technology according to claim 6, characterized in that: S4 includes: S41: Obtain multiple abnormal connected areas according to the final anomaly detection mask in the target monitoring area; S42: extracting the area and perimeter of each abnormal connected region and calculating the compactness, which is used to locate abnormal unit regions with unstable morphology; S43: When the compactness is lower than a preset compactness threshold, the corresponding abnormally connected area is marked as an abnormal evolution sensitive area, otherwise no marking is performed.

8. The method for sea area monitoring based on satellite remote sensing technology according to claim 7, characterized in that: The S4 also includes: S44: After monitoring, obtain the abnormal evolution sensitive area in the next monitoring cycle, compare the deviation values ​​of the abnormal units in the abnormal evolution sensitive area during each monitoring cycle, and if the deviation value exceeds the preset deviation threshold, automatically repeat S2 to correct the abnormal sensitivity value driven by the wind field and dynamically adjust the scattering anomaly confidence value and inverse sensitivity weight.

9. A sea area monitoring system based on satellite remote sensing technology, used to implement the sea area monitoring method based on satellite remote sensing technology as described in any one of claims 1 to 8, characterized in that: include, The data module is used to generate multi-period and multi-band remote sensing raster data by real-time monitoring of resource sets in the sea area; The preliminary identification module is used to obtain joint sample vectors at different times by monitoring the external environmental data within the target monitoring area. The multivariate kernel density estimation algorithm is used to obtain the four-dimensional joint probability density distribution function. Combined with multi-period and multi-band remote sensing raster data, the abnormal coverage area is obtained. The secondary identification module is used to identify and enhance abnormal units based on abnormal coverage areas and perform dynamic correction operations; The local optimization module is used to obtain abnormal connected areas after dynamic correction, analyze abnormal units with unstable morphology, determine the abnormal evolution sensitive areas, and implement local adaptive optimization methods based on the changes in the abnormal evolution sensitive areas.

Images