A sea area monitoring method and system based on satellite remote sensing technology

Abstract

The application discloses a kind of sea area monitoring method and system based on satellite remote sensing technology, it is related to remote sensing technical field, first by resource set generation multi-period multi-band remote sensing grid data, second, by the external environment data in the monitoring target monitoring area, obtain the joint sample vector on different time, using multivariate kernel density estimation algorithm, obtain four-dimensional joint probability density distribution function, combined with multi-period multi-band remote sensing grid data, obtain abnormal coverage area;Then, based on abnormal coverage area, identify enhanced anomaly unit, and carry out dynamic correction operation;Finally, after dynamic correction, obtain abnormal connected region, and analyze the unstable abnormal unit of form, determine abnormal evolution sensitive area, based on the change of abnormal evolution sensitive area, realize local adaptive optimization means.

Description

TECHNICAL FIELD

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

[0002] The present application relates to the technical field of remote sensing, in particular to a sea area monitoring method based on satellite remote sensing technology, mainly applied to automatic remote sensing monitoring of oil spill, illegal pollution and marine ecological abnormalities. In the broad technical field of marine information perception and monitoring, satellite remote sensing has become an important technical means for monitoring sea pollution and abnormal ecological processes due to its characteristics of large range, periodicity and multi-band acquisition of physical parameters. In the specific branch of satellite remote sensing technology, synthetic aperture radar (SAR), visible light multispectral, thermal infrared and other multi-source remote sensing images are widely used for discrimination and tracking of pollution events such as oil spill and pollution.

[0003] At present, in the monitoring of the above-mentioned sea area abnormal targets, single-period synthetic aperture radar (SAR) images are mainly relied on for abnormal identification. Although SAR has the advantages of all-weather and all-day imaging, in complex wind and wave environment, the backscattering changes caused by abnormal oil film, slight pollution or initial eutrophication are often covered up. For example, when the wind speed is high (more than 10 m / s), the sea surface layer of spray and foam will significantly enhance the scattering signal, so that the thin oil film cannot form a typical "black spot" feature in the SAR image, resulting in missed detection. On the contrary, when the wind speed is too low (less than 3 m / s), the mirror reflection effect is produced on the calm sea surface, which may also appear dark spots, which are confused with oil spill dark spots, causing misjudgment. This single radar band dependent monitoring method has obvious limitations in actual implementation and ecological risk control. SUMMARY

[0004] In view of the deficiencies of the prior art, the present application provides a sea area monitoring method and system based on satellite remote sensing technology, which solves the problems in the above background art.

[0005] To achieve the above purpose, the present application is realized by the following technical scheme: a sea area monitoring method based on satellite remote sensing technology, comprising the following steps,

[0006] S1: generating multi-period multi-band remote sensing raster data by monitoring resources in the sea area in real time;

[0007] S2: obtaining joint sample vectors at different times by monitoring external environmental data in the target monitoring area, obtaining a four-dimensional joint probability density distribution function by using a multivariate kernel density estimation algorithm, and obtaining an abnormal coverage area in combination with the multi-period multi-band remote sensing raster data;

[0008] S3: Based on the abnormal coverage area, identify the enhanced anomaly unit and carry out dynamic correction work;

[0009] S4: After dynamic correction, obtain the abnormal connected region, analyze the unstable anomaly unit, determine the abnormal evolution sensitive area, and realize local adaptive optimization means based on the change of the abnormal evolution sensitive area.

[0010] Preferably, S1 includes:

[0011] S11: Obtain all satellite resources for sea area monitoring within a predetermined monitoring period using an orbit prediction tool, and obtain a resource set;

[0012] S12: According to the available band capability and orbit transit timing of each satellite in the resource set, combined with the geometric coverage relationship of the global sea area, execute the greedy scheduling through the time sorting queue to generate a multi-satellite multi-mode scheduling list;

[0013] S13: Through the multi-satellite multi-mode scheduling list, respectively obtain SAR images, visible light multispectral images and thermal infrared images from different times and different satellites to generate an image library, and geometrically register each image in the image library. Specifically, use an image geometric registration tool to unify each image in the image library to the same map projection coordinate system, and resample to a unified spatial resolution. According to the spatial resolution, the size of the sea area grid is determined;

[0014] S14: After completing the geometric registration, organize the pixel values of the same place at different times and different bands into a four-dimensional matrix to form multi-period multi-band remote sensing raster data;

[0015] The time sorting queue is a time window formed by sorting each satellite's orbit transit timing from early to late.

[0016] Preferably, S2 includes:

[0017] S21: Obtain external environmental data in the target monitoring area through the buoy array, including wind speed, wind direction, wave height and short-wave spectral peak energy, and spatially interpolate the external environmental data according to the sea area grid to align it with the multi-period multi-band remote sensing raster data, while obtaining a multi-period external environmental sample set, wherein the multi-period external environmental sample set includes joint sample vectors at different times;

[0018] S22: Based on the multi-period external environmental sample set, a multivariate kernel density estimation algorithm is used to perform kernel density estimation on the joint sample vectors to obtain a four-dimensional joint probability density distribution function, which is used to quantify the probability density value of any joint sample vector appearing simultaneously in the historical sample.

[0019] Preferably, S2 further includes:

[0020] 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, which is used to quantify the cumulative occurrence probability of the current joint sample vector in the historical samples;

[0021] S24: Using the cumulative distribution function value, an abnormal sensitivity value based on the wind field driving corresponding to the current monitoring time and spatial position is obtained, which is used to reflect the rarity of the current monitored joint sample vector at the corresponding geographic location;

[0022] S25: Determine the SAR backscattering coefficient corresponding to each pixel point in the SAR image;

[0023] S26: At each sea area grid point, the local mean and standard deviation of the SAR backscattering coefficient in the sliding window centered on each sea area grid point are calculated, and the backscattering coefficient of each sea area grid point is normalized to obtain a scattering normalized anomaly value;

[0024] S27: Coupling the abnormal sensitivity value and the scattering normalized anomaly value, the scattering anomaly confidence value at each sea area grid point is calculated, specifically: , wherein is the scattering anomaly confidence value, is the abnormal sensitivity value, is the scattering normalized anomaly value;

[0025] S28: Based on the scattering anomaly confidence value at each sea area grid point, a scattering anomaly confidence surface is constructed;

[0026] S29: According to the adaptive quantile method, a dynamic local threshold surface is generated, and by comparing the scattering anomaly confidence surface with the dynamic local threshold surface, an anomaly mask is formed, specifically: , wherein is the anomaly mask at position , is the scattering anomaly confidence value at position , is the dynamic local threshold value at position ;

[0027] S291: The sea area grid points with an anomaly mask of 1 are taken as anomaly units, and by counting, an anomaly coverage area is obtained.

[0028] Preferably, S3 comprises:

[0029] S31: Based on the scattering anomaly confidence value, the inverse sensitivity weight of each sea area grid point in the target monitoring area is calculated;

[0030] ​​S32: determining a weak significant unit according to the abnormal unit, and extracting the sea water color index, the eutrophication index and the thermal infrared brightness temperature abnormal value in the weak significant unit;

[0031] S33: taking the inverse sensitivity weight as a weighting factor, fitting a local multiple regression model by a weighted least square method, and establishing a light-heat dual-domain prediction model based on the sea water color index, the eutrophication index and the thermal infrared brightness temperature abnormal value, so as to calculate the light-heat residual value on each weak significant unit.

[0032] Preferably, S3 further comprises:

[0033] S34: based on the light-heat residual value on each sea area grid point, taking the 95th percentile value of the light-heat residual value in a sliding window with the corresponding sea area grid point as the center as a local dynamic threshold, and if the light-heat residual value exceeds the local dynamic threshold, the corresponding sea area grid point is marked as an enhanced abnormal unit.

[0034] S35: dynamically correcting the abnormal coverage area through the enhanced abnormal unit, so as to obtain a final abnormal detection mask in the target monitoring area.

[0035] Preferably, S4 comprises:

[0036] S41: obtaining a plurality of abnormal connected regions according to the final abnormal detection mask in the target monitoring area.

[0037] S42: extracting the area and the perimeter in each abnormal connected region, and calculating the compactness, which is used to locate the abnormal unit area with unstable shape.

[0038] S43: when the compactness is lower than a preset compactness threshold, marking the corresponding abnormal connected region as an abnormal evolution sensitive area, otherwise, no marking processing is performed.

[0039] Preferably, S4 further comprises:

[0040] S44: after monitoring, obtaining the abnormal evolution sensitive area in the next monitoring period, comparing the deviation value of the abnormal unit in the abnormal evolution sensitive area between monitoring periods, and if the deviation value exceeds a preset deviation threshold, automatically repeating S2 to correct the abnormal sensitivity value based on the wind field driving, and dynamically adjusting the scattering abnormal confidence value and the inverse sensitivity weight.

[0041] A sea area monitoring system based on satellite remote sensing technology, comprising,

[0042] The data module is configured to generate multi-period multi-band remote sensing raster data by monitoring resources in the sea area in real time.

[0043] The preliminary identification module is used for obtaining joint sample vectors at different times by monitoring external environment data in a target monitoring area, obtaining a four-dimensional joint probability density distribution function by using a multivariate kernel density estimation algorithm, and obtaining an abnormal coverage area in combination with multi-period and multi-band remote sensing grid data;

[0044] The secondary identification module is used for identifying enhanced abnormal units based on the abnormal coverage area and performing dynamic correction work.

[0045] The local optimization module is used for obtaining an abnormal connected region after dynamic correction, analyzing abnormal units with unstable shapes, determining an abnormal evolution sensitive area, and realizing local adaptive optimization means based on changes in the abnormal evolution sensitive area.

[0046] The application provides a sea area monitoring method and system based on satellite remote sensing technology, and has the following beneficial effects:

[0047] (1) By step S1, multi-satellite resources are used to obtain multi-period remote sensing images at different times and different bands, geometric registration and resolution are unified, and multi-period and multi-band remote sensing grid data are formed, so that the monitoring area is completely covered in time, band and spatial distribution, and the traceability and multi-dimensional comparison ability of abnormal water area changes are significantly enhanced. By step S2, on the basis of collecting multi-period external environment data (such as wind speed, wind direction, wave height and short-wave spectral peak energy) and obtaining joint sample vectors at different times, a four-dimensional joint probability density distribution function is established by using a multivariate kernel density estimation algorithm, which effectively describes the rarity of the current monitoring environment in the historical data distribution, and then the abnormal coverage area is formed in combination with the remote sensing data. This method can adaptively identify potential abnormalities under different wind and wave fields, and improve the detection accuracy under complex environmental conditions. By step S3, enhanced abnormal units are further identified based on the abnormal coverage area, and dynamic correction is performed by using a light-heat multi-domain residual method, thereby reducing the SAR scattering blind area problem caused by high wind speed or surges, so as to realize effective detection of hidden oil films or weak pollution. By step S4, the abnormal connected region after dynamic correction is further analyzed, the shape parameters such as area, perimeter and compactness of the abnormal connected region are combined, the abnormal unit area with unstable shape is identified, and the abnormal evolution sensitive area is determined, and then the abnormal detection parameters are dynamically updated based on the periodic changes of the sensitive area, the local adaptive optimization of the abnormal sensitivity and threshold model is completed, and the problem of frequent fluctuation of the abnormal detection boundary in long-time sequence monitoring is effectively overcome. In summary, the method realizes a closed-loop monitoring system from data acquisition to model evolution through multi-period and multi-band multi-modal data fusion, environment-driven adaptive detection based on probability density and dynamic feedback local optimization mechanism, and improves the detection accuracy and spatial and temporal stability of abnormal water area phenomena such as oil spill and eutrophication.

[0048] (2) In steps S23 and S24, the cumulative distribution function value CDF is calculated by integrating the four-dimensional joint probability density distribution function, so as to quantify the cumulative occurrence probability of the current joint sample vector (i.e. the combination of wind speed, wind direction, wave height and short wave spectral peak energy) in the historical sample, and further inverted into an anomaly sensitivity value, which reflects the rarity of the current environmental combination, so that the subsequent anomaly detection automatically improves the sensitivity under rare environmental conditions, and reduces the missed detection rate under special meteorological-sea state coupling conditions. Through S25 to S29, first, the SAR backscattering coefficient is extracted at each sea area grid point, and in S26, the mean value and standard deviation of the local sliding window are normalized to obtain the scattering anomaly value which can measure the deviation degree of the water surface from 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, the global scattering anomaly confidence surface is formed, and in S29, the dynamic local threshold surface generated by the adaptive quantile method is compared to form the anomaly mask, and finally in S291, the anomaly unit is extracted and the anomaly coverage area is obtained. Through the multi-stage sensitive detection mechanism coupled with environmental rarity and local anomaly intensity, the ability to capture abnormal phenomena (such as weak oil film, slight pollution or early eutrophication) under complex wind and wave conditions is significantly improved. In summary, through the multi-stage cooperative processing of S21 to S291 of S2 described above: the dynamic adaptive anomaly sensitivity generation under environmental driving is realized, while maintaining the detection sensitivity, the spatial consistency of the anomaly detection result with the change of wind field, wave height and short wave spectrum is significantly enhanced, so that the final detection result can maintain high stability and reliability under varying sea conditions, providing a solid foundation for subsequent anomaly unit shape stability tracking and identification of anomaly evolution sensitive area, effectively breaking through the technical bottleneck of high false alarm and high missed detection rate of traditional methods in high dynamic environment.

[0049] (3) Through the inverse sensitivity weight guidance, photothermal dual-domain fitting residual analysis and dynamic quantile enhanced detection of S31 to S35, the method of the present application realizes the active attention allocation to low confidence units, the secondary deep digging of SAR weak anomaly area in the photothermal domain, and the adaptive correction of anomaly coverage area, improves the detection ability of hidden oil film, slight pollution and early eutrophication under complex wind field and multi-period sea state, and enhances the refinement and detection acuteness of the monitoring result in spatial distribution.

[0050] (4) In step S43, when it is detected that the compactness of a certain abnormal connected region is lower than the preset compactness threshold, it is determined that the boundary form of the region is loose and complex, and is prone to dramatic changes in form under different monitoring periods, and significant deviations caused by wind and wave disturbance. At this time, the connected region is marked as an abnormal evolution sensitive area, otherwise it is not marked. In this way, active spatial focusing on the form variable region in the abnormal detection result is achieved, which facilitates key tracking and analysis in subsequent periods. In step S44, after obtaining the data of the next monitoring period, for these abnormal evolution sensitive areas, the spatial coverage rate change of the abnormal units in the current period detection result and the historical average detection result is compared, and the deviation value is calculated; if the deviation value exceeds the preset deviation threshold, it means that the detection result of the region is significantly inconsistent with the historical pattern, which may be caused by insufficient learning of the long-term wind and wave-abnormal relationship in this area by the environmental statistical model or sensitivity calculation. At this time, S2 is automatically repeated to re-optimize the four-dimensional joint probability density distribution function based on the latest external environmental data and detection results in the region, and dynamically update the abnormal sensitivity value, scattered abnormal confidence value and inverse sensitivity weight, forming a closed-loop adaptive adjustment mechanism. In summary, steps S41 to S44 achieve the reverse deduction of abnormal variable regions from the form characteristics of the detection results, and automatically lock the regions 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-driven abnormal detection model in the long-term monitoring sequence, effectively avoiding the frequent boundary jumping, false alarm or missed alarm problems caused by fixed threshold and lack of time and space self-adaptation in traditional monitoring. BRIEF DESCRIPTION OF DRAWINGS

[0051] Figure 1 A flowchart of a sea area monitoring method based on satellite remote sensing technology according to the present application;

[0052] Figure 2 A partial logic diagram of a sea area monitoring method based on satellite remote sensing technology according to the present application;

[0053] Figure 3 A partial logic diagram of a sea area monitoring method based on satellite remote sensing technology according to the present application;

[0054] Figure 4 A block diagram of a sea area monitoring system based on satellite remote sensing technology according to the present application. DETAILED DESCRIPTION

[0055] The technical solutions in the embodiments of the present application will be described clearly and completely 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 the embodiments. 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.

[0056] Embodiment 1:

[0057] Referring to Figures 1 to 3 The application provides a sea area monitoring method based on satellite remote sensing technology, comprising the following steps,

[0058] S1: generating multi-period multi-band remote sensing grid data by monitoring resources in the sea area in real time;

[0059] S2: obtaining joint sample vectors at different times by monitoring external environmental data in the target monitoring area, obtaining a four-dimensional joint probability density distribution function by using a multivariate kernel density estimation algorithm, and obtaining an abnormal coverage area in combination with the multi-period multi-band remote sensing grid data;

[0060] S3: identifying an enhanced abnormal unit based on the abnormal coverage area and performing dynamic correction work;

[0061] S4: obtaining an abnormal connected region after dynamic correction, analyzing an abnormal unit with unstable morphology, determining an abnormal evolution sensitive area, and realizing local adaptive optimization means based on changes in the abnormal evolution sensitive area.

[0062] In the existing sea area monitoring based on SAR (synthetic aperture radar) satellites, it is commonly used for oil spill monitoring and ice boundary identification, etc. However, since SAR detection relies on the scattering characteristics of sea surface short wave roughness units, when the sea surface wind speed is too low (<3 m / s) or too high (>10 m / s), false positives or false negatives may occur in the image, causing distortion of the monitoring conclusion. This problem is particularly prominent in that the low wind speed sea surface is naturally smooth, which is easy to misjudge as an oil film black spot; under high wind speed, the oil film fails to suppress waves, making it difficult to be identified, resulting in missed detection. Therefore, a sea area monitoring method based on satellite remote sensing technology is proposed.

[0063] In this embodiment, the sea area monitoring method based on satellite remote sensing technology provided by the application relies on deep coupling analysis of multi-period multi-band remote sensing data and external driving data such as wind fields, which can effectively solve the problems of poor detection stability, high false alarm rate and easy to be covered by short period wind field disturbance when facing complex wind and wave conditions, changeable oil film and eutrophic abnormal water body in the existing single-period single-source monitoring mode. The specific beneficial effects are as follows:

[0064] The S1 step dynamically schedules the orbit timing, band capacity and sea area coverage geometry of multiple satellites in the resource set, for example, in a 5-day monitoring period, a greedy scheduling algorithm is used to schedule Sentinel-1 / 2, GF-3, Landsat-8 and other multi-satellite transit tasks, so that the target sea area forms a distributed multi-period data grid sequence in the time dimension.

[0065] Subsequently, the observations of different time and different sensors (such as SAR polarization, visible light multi-spectrum, and thermal infrared) are organized into a four-dimensional matrix to form multi-period and multi-band remote sensing grid data through image geometric registration and uniform spatial resolution (for example, 30m grid) processing.

[0066] In the S2 step, the joint sample vector is obtained and the abnormal coverage area is extracted based on the four-dimensional probability distribution. The external environmental data of the target monitoring sea area are obtained through the buoy array and the high-resolution numerical model, including wind speed, wind direction, wave height, and short-wave spectral peak energy. For example, in a certain monitoring, the wave height and short-wave spectral energy obtained through multiple buoys are abnormally high, indicating that the sea surface dynamics are abnormally active. These environmental data are interpolated into the grid system corresponding to the remote sensing grid at the same spatial position and different times to form a joint sample vector at the same spatial position and different times.

[0067] The multivariate Gaussian kernel density estimation is used and the bandwidth matrix is determined based on the Silverman rule to establish the four-dimensional joint probability density distribution function, so as to obtain the cumulative distribution function value CDF, which is further converted into the abnormal sensitivity value S=1-CDF.

[0068] The S value is coupled with the normalized value N of the local SAR scattering anomaly to obtain the scattering anomaly confidence value C, which is helpful to identify the grid cells with significant abnormal deviation in the multi-period and multi-band data, and form the preliminary abnormal coverage area.

[0069] In the S3 step, the enhanced abnormal cell detection and dynamic correction are based on the abnormal coverage area. Based on the abnormal coverage area, the inverse sensitivity weight is further introduced, and the light-thermal dual-domain (such as water color index, eutrophication index, and thermal infrared brightness temperature anomaly) residual analysis is performed on the weak significant cells except the preliminary abnormal cells. For example, in a certain monitoring, the water color index decreases, the eutrophication index increases, and the brightness temperature is significantly low, which jointly indicates a potential oil film.

[0070] Through the weighted least squares fitting model combined with the inverse sensitivity weight, the light-thermal residual value is calculated, and the enhanced abnormal cells are determined based on the local dynamic quantile threshold (such as P95) to form a more sensitive abnormal layer. In this way, not only the omission caused by the non-visibility of single SAR scattering is avoided, but also the sensitivity to light and thin or early oil films is significantly improved under high wind speed and surge conditions.

[0071] In the S4 step, the abnormal connected region formed by the abnormal cells is extracted through connectivity analysis, and the area, perimeter, and compactness (for example, when the compactness is less than 0.2, it indicates that the region morphology is highly irregular) are calculated. In this way, the abnormal region with unstable morphology and prone to boundary dramatic fluctuation can be effectively identified, and the abnormal evolution sensitive area is determined.

[0072] In the subsequent monitoring period, by comparing the spatial consistency deviation value (for example, the coverage deviation exceeds 0.25) of the anomaly detection in the sensitive area, the re-estimation of the anomaly sensitivity S of the corresponding environmental condition of the area is automatically triggered, so as to dynamically adjust the calculation link of the inverse sensitivity weight and the scattering anomaly confidence value, and realize local adaptive optimization.

[0073] In summary, through the step-by-step progression of S1 to S4, the present application forms a full-process technical chain from the spatio-temporal integration of multi-period multi-band data, the anomaly probability mapping driven by the wind field, the light-heat domain enhanced detection, and the anomaly morphology feedback closed loop, which improves the accuracy and stability of the sea oil film, pollution discharge and eutrophication anomaly detection under different wind and wave conditions, and effectively avoids the misjudgment or omission problem of traditional single-period SAR detection under the cover of high wind speed and surge and the confusion of low wind speed mirror.

[0074] Embodiment 2:

[0075] Please refer to Figure 1 Specifically, S1 includes:

[0076] S11: Obtain all satellite resources for sea monitoring in a predetermined monitoring period by using an orbit prediction tool (for example, professional orbit prediction software STK (Systems Tool Kit), GMAT (General Mission Analysis Tool)), and obtain a resource set;

[0077] S12: According to the available band capability and orbit transit timing of each satellite in the resource set, combined with the geometric coverage relationship of the global sea area, execute greedy scheduling through time-sequenced queue to generate a multi-satellite multi-mode scheduling list;

[0078] When formulating the monitoring scheduling plan, it is necessary to comprehensively consider which bands each satellite can observe (available band capability), when it passes through a specific place (transit timing), and which part of the sea area it can see from the orbit (geometric coverage), so as to arrange a reasonable multi-band, multi-period and multi-area coverage observation plan;

[0079] Among them, the available band capability refers to which wavelength of reflection or scattering information it can obtain, for example, Sentinel-2 cannot obtain SAR, and Sentinel-1 cannot obtain thermal infrared;

[0080] The greedy scheduling is to choose the optimal arrangement under the current condition each time, without considering the global optimum in the future, that is, it only pursues the current local optimum each time, unlike genetic scheduling which tries many combinations to find the global optimum, but prefers to use the opportunity that can complete the task earliest to quickly cover the demand; for example, Sentinel-1 can take SAR in the morning of the first day, and Sentinel-2 can take spectrum in the afternoon of the first day.

[0081] S13: Obtain SAR images, visible multispectral images, and thermal infrared images acquired from different satellites at different times through multi-satellite multi-mode scheduling list to generate an image library, and geometrically register each image in the image library. Specifically, use image geometric registration tools to unify each image in the image library to the same map projection coordinate system (such as WGS84 / UTM Zone) and resample to a unified spatial resolution (e.g., 10m or 30m grid). According to the spatial resolution, determine the size of the sea area grid.

[0082] Spatial resolution refers to the size of the area represented by each pixel on the ground. The set spatial resolution of 10m or 30m determines the geographic range of each grid, i.e., the actual area size represented by each sea area grid.

[0083] Resampling to a unified spatial resolution means adjusting the spatial resolution of the image to a unified scale to ensure consistent actual ground coverage of each sea area grid.

[0084] SAR images are used for black spot anomaly monitoring, visible multispectral images are used for turbidity and algal bloom indication monitoring, and thermal infrared images are used for temperature and pollution discharge anomaly detection.

[0085] S14: After geometric registration, organize pixel values of the same location at different times and different bands into a four-dimensional matrix to form multi-period multi-band remote sensing grid data for subsequent arbitrary slice queries (e.g., viewing changes over time at the same location or different band responses at the same time).

[0086] The time ordering queue is a time window formed by sequentially ordering the orbits of each satellite from early to late.

[0087] The multi-satellite multi-mode scheduling list specifies which satellite to acquire which type of image data at different time nodes, providing scheduling basis for subsequent time sequence multi-band acquisition.

[0088] The four-dimensional matrix includes spatial dimensions, time dimensions, and band dimensions. The spatial dimensions include longitude and latitude.

[0089] Multi-period multi-band remote sensing grid data is pixel matrix data containing multiple bands acquired by remote sensing sensors at multiple different time points for the same monitoring area. Each pixel records the observation value of its corresponding geographic location at each band with a time label, which can be directly used for time series analysis, spectral comparison, or machine learning.

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

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

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

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

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

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

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

[0097] Example 3:

[0098] Please refer to Figure 1 , specifically: S2 includes:

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

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

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

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

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

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

[0105] 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;

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

[0107] The specific expression of the four-dimensional joint probability density distribution function is:

[0108] ;

[0109] where, is the four-dimensional joint probability density distribution function, representing the probability density value of the joint sample vector (i.e. the 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 short-crested wave peak energy, m is the total number of historical samples, i is the index of the historical sample, H is a four-dimensional bandwidth matrix (a 4x4 positive definite matrix) characterizing the smoothness and cross-correlation between the four variables, is the determinant of H, used for the normalization of the multi-dimensional volume, is the inverse of the matrix square root of H, 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 short-crested wave peak energy of the i-th sample, is the current environmental combination point (i.e. the current joint sample vector) whose probability density is to be estimated, is the four-dimensional kernel function, which is generally chosen to be the four-dimensional Gaussian kernel function (Gaussian distribution function) and is often written as: where the exponent term is the squared norm because it is four-dimensional, is the kernel function value, in the kernel density estimation, is used to measure the local density contribution at a given normalized distance z, z represents the normalized offset vector scaled by the bandwidth matrix H, which scales the difference between the four variables at each sample point and the target point to a unified kernel scale, and is a four-dimensional column vector, represents the squared norm of the vector, that is, the squared Euclidean distance of the z vector, which represents the squared distance from the center in the Gaussian distribution exponent term, used to control the density decay. e is the natural base number, with a value of about 2.71828; is the normalization factor, used to ensure that the kernel density function integrates to 1 over the entire four-dimensional space;

[0110] represents scaling the difference vector with 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 environmental combination;

[0111] The multivariate kernel density estimation uses a Gaussian kernel function and determines the bandwidth matrix based on the Silverman rule;

[0112] The multivariate kernel density estimation algorithm is used to estimate the joint probability density of multiple random variables simultaneously.

[0113] S2 further comprises:

[0114] S23: According to the four-dimensional joint probability density distribution function, the total 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, which is used to quantify the cumulative occurrence probability of the current joint sample vector in the historical samples;

[0115] The cumulative distribution function value is integrated under the corresponding four-dimensional density to obtain the cumulative probability of the occurrence of the current combination, that is, the cumulative proportion of samples in history that simultaneously satisfy ;

[0116] S24: Using the cumulative distribution function value, obtain the abnormal sensitivity value based on the wind field driving corresponding to the current monitoring time and spatial position, specifically: , wherein 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 current monitored joint sample vector (i.e. the environment combination) at the corresponding geographic location;

[0117] If the CDF value is large (close to 1), it means that the current environment combination appears frequently in the historical data and belongs to the normal state, and if the CDF value is small (close to 0), it means that it is less common, so the abnormal sensitivity value represents the degree of rarity of the current environment combination in history;

[0118] S25: Determine the SAR backscattering coefficient corresponding to each pixel point in the SAR image;

[0119] In a synthetic aperture radar (SAR) image, the value of each pixel is actually the backscattering 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, affected by sea surface roughness, wind speed, oil film, waves, etc.

[0120] S26: Based on the SAR backscattering coefficient in the sliding window centered on each sea area grid, calculate the local mean and standard deviation, and normalize the backscattering coefficient of each sea area grid to obtain the scattering normalized anomaly value;

[0121] The S26 step is equivalent to doing Z-score, which represents how many standard deviations the backscattering of the corresponding point is higher than the average value around it, which can better discover local anomalies. If the oil film causes the scattering to decrease significantly, the scattering normalized anomaly value will be significantly lower than 0, and vice versa, it will be significantly higher than 0;

[0122] When performing local statistics (such as local mean and standard deviation), a small rectangular area (such as 3x3, 5x5, 7x7 grids) is formed by expanding around a certain sea area grid, and the values in this area are counted.

[0123] The scattering normalized anomaly value is used to measure the deviation of the corresponding sea area grid point from the local normal water surface;

[0124] S27: Coupling the anomaly sensitivity value and the scattering normalized anomaly value, calculating the scattering anomaly confidence value on each sea area grid point, specifically: , wherein, is the scattering anomaly confidence value, is the anomaly sensitivity value, is the scattering normalized anomaly value, which is the abnormal degree of the corresponding sea area grid point SAR backscattering coefficient relative to the local;

[0125] The scattering anomaly confidence value is used to form a scattering anomaly confidence surface in the monitoring area, providing a basis for subsequent generation of local dynamic quantile-based anomaly detection threshold and extraction of suspected anomaly mask;

[0126] S28: Based on the scattering anomaly confidence value on each sea area grid point, a scattering anomaly confidence surface is constructed;

[0127] The scattering anomaly confidence surface refers to the confidence map of the anomaly detection result calculated by the wind field (wind speed, wind direction), wave (wave height, short wave spectral peak energy) and other external environmental driving factors, combined with the backscattering coefficient of SAR image, at each grid point of the sea area. This map is used to represent the scattering anomaly area on the sea surface, and can identify abnormal phenomena such as oil spill or ice floe.

[0128] S29: Generating a dynamic local threshold surface according to the adaptive quantile method, and forming an anomaly mask by comparing the scattering anomaly confidence surface with the dynamic local threshold surface, specifically: , wherein, is the anomaly mask at position , is the scattering anomaly confidence value at position , is the dynamic local threshold at position , is the position coordinate on the monitoring space;

[0129] The adaptive quantile method analysis specifically includes: calculating the quantile value (e.g. P90, P95, P98) of the scattering anomaly confidence value at all positions , and then adaptively setting a local dynamic threshold for different sea area environments:

[0130] , wherein, represents the position ​​​The quantile adjustment coefficient is dynamically set according to local wind speed, wind direction, wave height, short-wave spectral energy and other environmental characteristics, and usually fluctuates within a certain range (such as 0.85-1.15) to adapt to the sensitivity adjustment under different sea conditions. The 95th percentile scattering anomaly confidence value obtained by statistics in the target monitoring area is used to provide a global reference quantile benchmark value.

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

[0132] The cumulative distribution function value CDF is obtained 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.

[0133] In this embodiment, external environmental data at different time points, including wind speed, wind direction, wave height and short-wave spectral peak energy, are obtained by S21 using a numerical model (such as a wind wave field model based on atmospheric-ocean coupled dynamics) and a buoy array deployed in the target sea area. For example, when monitoring the sea area near the Zhoushan Islands, the wind speed field at a height of 10 m can be derived in real time from the numerical model, and the corresponding measured wave height can be obtained from the buoy data.

[0134] The short-wave spectral peak energy is a physical quantity used to describe the energy of the short-wave band on the sea surface, which directly affects the imaging sensitivity of SAR to fine oil films;

[0135] In the S21 step, these data are spatially interpolated according to the sea area grid (i.e. the fixed grid position after geometric registration) to ensure that different data sources are strictly aligned in space with multi-period and multi-band remote sensing grid data, laying a unified foundation for subsequent pixel-level multi-source coupling.

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

[0137] At S23, the cumulative distribution function value CDF is obtained by integrating the four-dimensional joint probability density function in the four-dimensional space to count the proportion of all historical probabilities that are lower than the current joint sample vector. For example, when CDF = 0.8, it means that the current marine environment is more common than 80% of the historical samples.

[0138] S24 defines the anomaly sensitivity value S using the formula S = 1 - CDF, which is used to quantify the rarity of the current environment. If S = 0.2, it means that the environment is relatively common, and if S = 0.95, it means that the environment is extremely rare in history, and more attention should be paid to the remote sensing anomalies monitored.

[0139] At S25, based on the SAR image, the SAR backscatter coefficient of each pixel point is extracted to characterize the radar scattering characteristics of the point. The lower the backscatter coefficient, the smoother the sea surface, which may be caused by oil film, pollution, etc.

[0140] S26 calculates the local mean and standard deviation in a sliding window centered on each sea area grid point to achieve local normalization and obtain the scattering normalized anomaly value N, which is used to characterize the deviation of the pixel point from the surrounding normal water surface. For example, if N = 2, it means that the pixel is two standard deviations higher than the surrounding area.

[0141] In S27, the anomaly sensitivity value S and the scattering normalized anomaly value N are coupled to calculate the scattering anomaly confidence value. If a point is in an extremely rare wind wave combination (S is high) and its scattering anomaly is significant (N is high), the scattering anomaly confidence value C is higher, representing a high degree of confidence in the anomaly.

[0142] In S28, a scattering anomaly confidence surface is constructed based on the scattering anomaly confidence value C in the entire monitoring area. In S29, a dynamic threshold surface is generated in a sliding window by local dynamic quantile (such as the 95th percentile), and the scattering anomaly confidence surface C(x, y) is compared 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 the anomaly coverage area is counted.

[0143] In summary, S21 aligns the external environment samples in space and accumulates them in time, S22-S24 quantifies the rarity of the environment combination in the probability space and converts it into sensitivity S, S25-S27 normalizes the backscatter in the physical scattering space and couples the environmental rarity to form the confidence, S28-S291 dynamically generates a threshold in the local statistical domain, extracts anomaly units and forms an anomaly coverage area. This not only improves the adaptability to different environmental conditions and avoids the inherent defects of single threshold strategy, but also effectively identifies oil film, pollution and water eutrophication anomalies that may be hidden under high or low wind speed through wind wave driving and radar scattering coupling, improving the accuracy and stability of the overall detection.

[0144] Embodiment 4:

[0145] Please refer to Figure 1 , specifically: S3 comprises:

[0146] S31: based on the scattering anomaly confidence value, calculating the inverse sensitivity weight of each sea area grid point in the target monitoring area;

[0147] The inverse sensitivity weight is obtained in the following manner: A=1-C, wherein A is the inverse sensitivity weight, which is used to form spatial attention guidance to SAR low confidence units;

[0148] S32: determining weak significant units according to the anomaly units, and extracting the sea water color index, eutrophication index and thermal infrared brightness temperature anomaly value in the weak significant units;

[0149] By extracting the sea water color index, eutrophication index and thermal infrared brightness temperature anomaly value in the weak significant units, the light-thermal dual-domain supplement to the SAR detection blind spot is formed, and the spatial coverage and sensitivity of the overall sea area anomaly detection are improved.

[0150] Weak significant units refer to units other than anomaly units;

[0151] The sea water color index is calculated according to the ratio of near-infrared to red band of multi-spectral pixels, and the specific calculation method is as follows: , wherein, is the sea water color index, is the near-infrared band reflectance pixel value at the monitoring sea area grid point, which is usually used to capture the absorption and reflection characteristics of algae, phytoplankton and suspended particles in water to near-infrared; is the red band reflectance pixel value at the monitoring sea area grid point, which is used for sensitive analysis of water, chlorophyll and other red light absorption;

[0152] The near-infrared band reflectance pixel value refers to the apparent surface reflectance directly received by the satellite sensor at the monitoring sea area grid point and after radiation correction in the near-infrared band (usually 0.76–0.90 μm).

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

[0154] In seawater, pure water has little reflection to near-infrared (strong absorption) and strong absorption to red light, so the sea water color index is low. If there are a large number of algae and suspended particles, the relative reflection of near-infrared and red light will be changed, resulting in an increase in the sea water color index value. Therefore, the sea water color index is used to distinguish abnormal water (algal blooms, silt turbidity) from normal water in the spectral dimension.​

[0155] 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;

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

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

[0158] 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;

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

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

[0161] 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 regression coefficients, representing the influence intensity of NDVI, FUI, ΔTIR on the predicted value; obtained by fitting the global data according to the weighted least squares method;

[0162] By inverse sensitivity weight weighted least squares, the light-heat residual is more sensitive to the low confidence area of SAR.

[0163] Calculate the light-heat residual value: ; wherein, is the light-heat residual value, is the inverse sensitivity weight at the position , is the actual light-heat index value directly obtained from the remote sensing image (or multi-band combination), representing the true observation, directly representing the actual spectral-thermal infrared combined anomaly extracted from the satellite image pixel; is the light-heat dual-domain prediction model, used to output the regression predicted value, representing the expected observed light-heat value at the sea area grid according to the light-heat feature multivariate fitting;

[0164] In this way, the light-heat anomaly in the SAR insignificant place will be amplified.

[0165] S3 also includes:

[0166] S34: At each sea area grid, based on the light-heat residual value, use the 95% quantile value of the light-heat residual value in the sliding window centered on the corresponding sea area grid as the local dynamic threshold, if the light-heat residual value exceeds the local dynamic threshold, the corresponding sea area grid is marked as an enhanced anomaly unit;

[0167] S35: Through the enhanced anomaly unit, dynamically correct the anomaly coverage area to obtain the final anomaly detection mask in the target monitoring area, thereby significantly improving the detection rate of hidden oil film or weak pollution under high wind speed and surge conditions.

[0168] The final anomaly detection mask is also marked as 1 at the enhanced anomaly unit based on the anomaly unit;

[0169] The enhanced anomaly unit belongs to a weak significant unit;

[0170] In this embodiment, after obtaining the scattering anomaly confidence value through step S31, the present application calculates the inverse sensitivity weight A = 1 - C of each sea area grid, which is used to guide the attention degree of the low confidence unit in the subsequent analysis. For example, the scattering anomaly confidence value C at some grid has reached 0.9, indicating that the anomaly here is highly reliable, then the inverse sensitivity weight A = 0.1, the subsequent analysis weight is reduced; while for the grid with C = 0.2, the inverse sensitivity weight A = 0.8, which is used to emphasize the detection sensitivity of this kind of weak significant area in the light-heat analysis, so as to effectively focus on the potential potential leakage detection area.

[0171] By step S32, the sea water color index is extracted in the weak significant unit except the abnormal unit, and the sea water color index value is calculated by the ratio of the near-infrared and red wave bands of the multi-spectral pixel, for example, when the sea water color index value is significantly increased, it may indicate that the water body has a large number of phytoplankton; the eutrophication index is calculated by the ratio of the green and red wave bands, for example, when the green wave band is enhanced, it indicates that the water bloom occurs; the thermal infrared brightness temperature abnormal value is the deviation of the thermal infrared brightness temperature value and the sliding mean value of the neighborhood, for example, the water temperature abnormality may be related to the pollution discharge, and these light-thermal multi-domain parameters can reveal the surface biological or temperature perturbation which is difficult to be detected by SAR.

[0172] In step S33, the inverse sensitivity weight is taken as a weighting factor, and a local multiple regression model with the sea water color index, the eutrophication index and the thermal infrared brightness temperature abnormal value as inputs is fitted by using the weighted least square method, which is used to predict the light-thermal response of the weak significant unit under normal conditions. For example, if the light-thermal fitting residual is abnormally increased at a position with high inverse sensitivity weight (such as A=0.85), it means that even if the SAR confidence is low, there may be light-thermal coupling abnormality in this region, thereby compensating for the deficiency of SAR in the detection of thin oil film or early pollution discharge.

[0173] In step S34, the light-thermal residual value is used as the basis, and the 95% quantile of the sliding window with each sea grid point as the center is used as the local dynamic threshold. For example, in a eutrophic water area, if most of the residuals are concentrated in 1.0-1.5, the 95% quantile reaches 1.7, and when the residual value of a certain grid point reaches 2.2, it can be marked as an enhanced abnormal unit, thereby avoiding the problem that the local feature fluctuation is covered due to the use of a global fixed threshold.

[0174] By step S35, the above enhanced abnormal unit is superimposed on the original abnormal coverage area to form a final abnormal detection mask, that is, on the basis of the original abnormal unit, the light-thermal residual significantly abnormal position is also marked as 1. For example, under the condition of high wind speed (>10 m / s), the original SAR abnormality is difficult to appear, but through this method, the weak pollution discharge area or oil film can still be accurately identified in the eutrophication abnormal or temperature abnormal area, and the detection rate and monitoring stability under complex sea conditions are significantly improved.

[0175] In summary, the present application uses the inverse sensitivity weight to adaptively guide attention in space, fuses the multi-dimensional information of the water color index, eutrophication and thermal infrared residual to construct the light-thermal dual-domain prediction, and performs local threshold discrimination based on the dynamic quantile method, thereby effectively overcoming the detection blind area of single SAR or global fixed threshold under high wind and wave, complex water area, and improving the detection rate and detection stability of weak oil spill, slight pollution and early eutrophication and other abnormal phenomena.

[0176] Example 5:

[0177] Please refer to Figure 1Specifically, S4 comprises:

[0178] S41: obtaining a plurality of abnormal connected regions according to the final abnormality detection mask in the target monitoring region;

[0179] S42: extracting the area (i.e. the number of pixels) and the perimeter (i.e. the boundary length) in each abnormal connected region, and calculating the compactness, which is used to locate the abnormal unit region with unstable shape and is an important geometric index for detecting the stability of abnormal shape;

[0180] The compactness is obtained by the following formula: wherein, is the compactness, is the area, is the perimeter, is the constant pi;

[0181] The abnormal connected region refers to a region in which all adjacent regions belong to the abnormality detection mask;

[0182] S43: when the compactness is lower than a preset compactness threshold, marking the corresponding abnormal connected region as an abnormal evolution sensitive area, otherwise, not marking.

[0183] S4 further comprises:

[0184] S44: after monitoring, obtaining the abnormal evolution sensitive area in the next monitoring period, comparing the deviation value of the abnormal unit in the abnormal evolution sensitive area between monitoring periods, and if the deviation value exceeds a preset deviation threshold, automatically repeating S2 to correct the abnormal sensitivity value based on the wind field driving, and dynamically adjusting the scattering abnormal confidence value and the inverse sensitivity weight.

[0185] The deviation value is the difference between the current period detection result and the historical average detection result in the abnormal unit distribution in the abnormal evolution sensitive area, which can be expressed by the change amount of spatial coverage;

[0186] In the embodiment, in step S41, a plurality of abnormal connected regions (i.e. closed regions composed of all adjacent pixels being abnormal units) are extracted according to the final abnormality detection mask (i.e. the abnormal marker layer integrated after SAR scattering and light-heat enhancement) in the target monitoring region. For example, if several oil film dark spots along the coast and scattered foam disturbance regions on the sea surface are formed in one monitoring, different abnormal patches can be identified by region connectivity.

[0187] In S42, the area (number of pixels) and the perimeter (boundary length) of each abnormal connected region are further extracted, and the compactness is calculated. The compactness is a key indicator representing the geometric compactness of the region, and is used to quantify the shape regularity of the abnormal patch. For example, an oil film patch close to a circle has a large area and a relatively short perimeter, and the compactness is close to 1; while a long and winding pollution belt has a small area but a long perimeter, and the compactness is much lower than 1. Through this step, those abnormal units with broken morphology, complex boundary, easy to split or shrink in multiple stages can be automatically located.

[0188] In S43, if the compactness of a certain abnormal connected region is lower than a preset tightness threshold (for example, 0.5), the region is marked as an abnormal evolution sensitive area. For example, if a long and branched pollution belt is repeatedly monitored near a port, and the compactness is lower than the threshold, the region will be automatically marked as a sensitive area for tracking its future morphological evolution. On the contrary, large-area abnormal patches with round shape and stable boundary will not be included in the sensitive area marking, so as to avoid excessive adjustment of natural stable abnormal regions.

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

[0190] In summary, through the cascade operation of S41 to S44, the present application can not only quantify the stability of abnormal units based on geometric morphology in a single period, but also automatically correct the detection model through the deviation of the abnormal evolution sensitive area in multiple periods, forming a local dynamic optimization closed loop. In long-term sea area monitoring tasks, especially in sea areas affected by tidal current, local wind field or periodic human pollution, the detection boundary can be significantly reduced, and the abnormal error rate can be reduced, so as to improve the spatial consistency of the whole monitoring system in time sequence and the stability of long-term abnormal discovery.

[0191] The abnormal connected region is used to identify spatially continuous abnormal unit patches.

[0192] Compactness is used to quantify the geometric complexity of abnormal units, and to detect whether the shape is variable.

[0193] Abnormal evolution sensitive area is the key dynamic tracking area screened out by compactness, and is used for adaptive adjustment of subsequent multi-period monitoring.

[0194] Deviation value is used to measure the significant difference of abnormal distribution in different periods in the abnormal evolution sensitive area, and to determine whether local re-update of the detection model is needed.

[0195] Embodiment 6:

[0196] Please refer to Figure 4 , Specifically: a sea area monitoring system based on satellite remote sensing technology, comprising,

[0197] The data module is used to generate multi-period multi-band remote sensing raster data by monitoring the resource set in the sea area in real time;

[0198] The preliminary identification module is used to obtain joint sample vectors at different times by monitoring the external environmental data in the target monitoring area, and to obtain a four-dimensional joint probability density distribution function by using a multivariate kernel density estimation algorithm, and to obtain an abnormal coverage area by combining the multi-period multi-band remote sensing raster data.

[0199] The secondary identification module is used to identify enhanced abnormal units based on the abnormal coverage area and to perform dynamic correction work.

[0200] The local optimization module is used to obtain abnormal connected regions after dynamic correction, to analyze abnormal units with unstable shapes, to determine abnormal evolution sensitive areas, and to realize local adaptive optimization means based on changes in the abnormal evolution sensitive areas.

[0201] Although the embodiments of the present application have been shown and described, it can be understood by those of ordinary skill in the art that various changes, modifications, replacements and variations can be made to these embodiments without departing from the principles and spirits of the present application, and the scope of the present application is defined by the appended claims and their equivalents.

Claims

1. A sea area monitoring method based on satellite remote sensing technology, characterized in that: The method comprises the following steps, S1: generating multi-period multi-band remote sensing grid data by monitoring resource sets in the sea area in real time; S2: obtaining joint sample vectors at different times by monitoring external environment data in the target monitoring area, obtaining a four-dimensional joint probability density distribution function by using a multivariate kernel density estimation algorithm, and obtaining an abnormal coverage area by combining the multi-period multi-band remote sensing grid data; S3: identifying enhanced abnormal units based on the abnormal coverage area and performing dynamic correction work; S4: obtaining an abnormal connected region after dynamic correction, analyzing abnormal units with unstable shapes, determining an abnormal evolution sensitive area, and realizing local adaptive optimization based on changes in the abnormal evolution sensitive area; S1 comprises: S11: obtaining all satellite resources for sea area monitoring in a predetermined monitoring period by using an orbit prediction tool to obtain a resource set; S12: generating a multi-satellite multi-mode scheduling list by using a time sorting queue to perform a greedy scheduling according to the available band capabilities and orbit transit timing of each satellite in the resource set and the geometric coverage relationship of the global sea area; S13: generating an image library by obtaining SAR images, visible light multispectral images and thermal infrared images from different times and different satellites through the multi-satellite multi-mode scheduling list, and geometrically registering each image in the image library, specifically: using an image geometric registration tool to unify each image in the image library to the same map projection coordinate system and resample to a unified spatial resolution, and determining the size of the sea area grid according to the spatial resolution; S14: after completing the geometric registration, organizing the pixel values of the same place at different times and different bands into a four-dimensional matrix to form multi-period multi-band remote sensing grid data; wherein the time sorting queue is a time window sorted from early to late according to the orbit transit timing of each satellite; S2 comprises: S21: obtaining external environment data in the target monitoring area by using a buoy array, including wind speed, wind direction, wave height and short-wave spectral peak energy, and spatially interpolating the external environment data according to the sea area grid to align it with the multi-period multi-band remote sensing grid data, while obtaining a multi-period external environment sample set, wherein the multi-period external environment sample set includes joint sample vectors at different times; S22: based on the multi-period external environment sample set, using a multivariate kernel density estimation algorithm to perform kernel density estimation on the joint sample vectors to obtain a four-dimensional joint probability density distribution function, which is used to quantify the probability density value of any joint sample vector appearing simultaneously in the historical samples; S2 further comprises: S23: according to the four-dimensional joint probability density distribution function, counting the total sum of all historical probability density values lower than the current joint sample vector to obtain a cumulative distribution function value, which 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 to obtain the abnormal sensitivity value based on the wind field at the corresponding monitoring time and spatial position, which is used to reflect the rarity of the currently monitored joint sample vector at the corresponding geographical location; S25: Determine the SAR backscattering coefficient corresponding to each pixel point in the SAR image; S26: Calculate the local mean and standard deviation of the SAR backscattering coefficient in the sliding window centered on each sea area grid point, and normalize the backscattering coefficient of each sea area grid point to obtain the scattering anomaly confidence value; S27: coupling the anomaly sensitivity value with the scattering normalized anomaly value, calculating the scattering anomaly confidence value at each grid point of the sea area, specifically: wherein, is the scattering anomaly confidence value, is the anomaly sensitivity value, is the scattering normalized anomaly value; S28: Construct a scattering anomaly confidence surface based on the scattering anomaly confidence values of each sea area grid point; S29: generating a dynamic local threshold surface according to an adaptive quantile method, and forming an anomaly mask by comparing the scattering anomaly confidence surface with the dynamic local threshold surface, specifically: wherein, is an anomaly mask at position , is a scattering anomaly confidence value at position , is a dynamic local threshold at position . S291: Take the sea area grid points with the anomaly mask value of 1 as the anomaly units, and obtain the abnormal coverage area by statistical analysis; S3 includes: S31: Calculate the inverse sensitivity weight of each sea area grid point in the target monitoring area based on the scattering anomaly confidence value; S32: Determine the weak significant unit based on the anomaly unit, and extract the sea color index, eutrophication index, and thermal infrared brightness temperature anomaly value in the weak significant unit; S33: Use the inverse sensitivity weight as a weighting factor to fit a local multiple regression model by weighted least squares method, and establish a light-heat dual-domain prediction model based on the sea color index, eutrophication index, and thermal infrared brightness temperature anomaly value to calculate the light-heat residual value in each weak significant unit; S3 further includes: S34: At each sea area grid point, use the 95th percentile value of the light-heat residual value in the sliding window centered on the corresponding sea area grid point as the local dynamic threshold based on the light-heat residual value, and if the light-heat residual value exceeds the local dynamic threshold, the corresponding sea area grid point is marked as an enhanced anomaly unit; S35: Dynamically correct the abnormal coverage area through the enhanced anomaly unit to obtain the final anomaly detection mask in the target monitoring area; S4 includes: S41: Obtain a plurality of abnormal connected regions according to the final anomaly detection mask in the target monitoring area; S42: Extract the area and perimeter of each abnormal connected region, calculate the compactness, and use the compactness to locate the abnormal unit region with unstable shape; S43: When the compactness is lower than the preset compactness threshold, mark the corresponding abnormal connected region as an abnormal evolution sensitive area, otherwise do not mark it; S4 further includes: S44: After monitoring, obtain the abnormal evolution sensitive area in the next monitoring period, compare the deviation values of the abnormal units in the abnormal evolution sensitive area between monitoring periods, and if the deviation value exceeds the preset deviation threshold, automatically repeat S2 to correct the abnormal sensitivity value based on the wind field and dynamically adjust the scattering anomaly confidence value and the inverse sensitivity weight.

2. A sea area monitoring system based on satellite remote sensing technology, used to realize the sea area monitoring method based on satellite remote sensing technology in claim 1, characterized in that: includes, The data module is configured to generate multi-period and multi-band remote sensing raster data by monitoring resources in the sea area in real time; The preliminary identification module is configured to obtain a joint sample vector at different times by monitoring external environmental data in the target monitoring area, obtain a four-dimensional joint probability density distribution function by using a multivariate kernel density estimation algorithm, and obtain an abnormal coverage area in combination with the multi-period and multi-band remote sensing raster data; The secondary identification module is configured to identify enhanced anomaly units based on the abnormal coverage area and perform dynamic correction work. The local optimization module is used to obtain abnormal connected regions after dynamic correction, analyze abnormal units with unstable morphology, determine abnormal evolution sensitive areas, and realize local adaptive optimization means based on the changes of the abnormal evolution sensitive areas.

Images