Coastal zone human activity interference remote sensing quantification and ecological profit and loss evaluation method and device

Abstract

The invention discloses a coastal zone human activity interference remote sensing quantification and ecological profit and loss evaluation method and device, belongs to the technical field of remote sensing monitoring, and solves the problems of low efficiency due to dependence on field investigation and insufficient reliability caused by separation of ecological evaluation and remote sensing data in the prior art. The method is characterized by comprising the following steps: acquiring two-stage multispectral remote sensing images containing multiple wavebands such as blue light at an interval of 5-10 years and with spatial resolution of 10-30m, and synchronously integrating tide, elevation and geological background data; eliminating environmental errors of image radiometric calibration and atmospheric correction; calculating a normalized building index based on the short-wave and near-infrared band reflectivity, and detecting a quantized interference area in combination with land utilization classification change; and the change area and the corrected intensity change value are extracted and input into an ecological system service value model, and an evaluation result is output by using a unit area value equivalent factor method. According to the method, full-process remote sensing automatic processing is realized, field investigation is avoided, the difference adjustment cost is effectively reduced, and the pertinence of management decision is remarkably improved.

Description

TECHNICAL FIELD

[0001] The present application belongs to the technical field of remote sensing monitoring, and particularly relates to a method and device for quantifying human activity interference and evaluating ecological loss and gain of a coastal zone based on multispectral remote sensing physical property analysis. BACKGROUND

[0002] Existing quantification of human activity interference in the coastal zone mainly relies on field sampling and manual investigation. The coastal zone has a large spatial span and complex terrain, and on-site investigation requires a large amount of manpower and time cost. The accessibility of intertidal zones and nearshore waters is poor, resulting in spatial blind spots in data collection. In multi-period monitoring tasks, on-site work is difficult to achieve time synchronization, and different year data may have comparability errors due to sampling point deviation or standard deviation. In addition, on-site investigation can only obtain discrete point information, and it is difficult to completely depict the spatial heterogeneity characteristics of human activity interference.

[0003] More importantly, the unique tidal dynamic environment of the coastal zone poses a severe challenge to remote sensing monitoring. Conventional remote sensing land use classification methods are prone to misjudgment of periodically submerged areas (such as low-tide exposed beaches) as permanent land class changes (such as bare land or artificial ground), while the same area at high tide is classified as water body, resulting in a large number of false change patches in multi-period image comparison, which seriously interferes with the accurate extraction of the real interference area of human activities. At the same time, the problems of image data loss, reflection anomalies (such as tidal ditch mirror reflection), and noise caused by tidal phenomena also destroy the continuity and accuracy of the time series of key parameters such as vegetation index, further reducing the reliability of change detection.

[0004] Ecological loss and gain evaluation usually adopts the static value coefficient method, but the existing technology fails to establish a quantitative conversion mechanism from the physical quantity changes (such as normalized building index NDBI change, land use change area) obtained from remote sensing monitoring to the loss of ecosystem service value. Although remote sensing data can provide information on large-scale land changes, these physical quantity indicators (such as NDBI only reflect the physical change of building density) cannot automatically map to ecological value loss (such as loss of biodiversity, decrease of carbon sink capacity). Traditional methods require manual input of interference range data obtained from field statistics or simple interpretation into the evaluation model, which is prone to human error. This disconnection between remote sensing physical quantities and ecological value evaluation results in a lack of response to the ecological impact of differences in coastal zone development intensity (such as high-intensity reclamation and low-intensity facilities) and spatial patterns (such as centralized development and fragmentation development), making it difficult to reflect the nonlinear cumulative damage effect of high-intensity or sensitive area development on ecological function.

[0005] The difficulty of realizing automatic evaluation lies in: first, the above tidal dynamic influence makes the accuracy of conventional remote sensing classification and change detection result insufficient, which directly affects the reliability of interference area extraction. Second, the ecological effect of human activity interference has nonlinear characteristics, and simple area statistics cannot represent the cumulative impact of fragmentation development and other modes, and a quantitative correlation model of interference strength and ecological value decay needs to be established. The existing technology has not solved the cross-scale conversion problem from remote sensing physical quantity to ecological parameter. SUMMARY

[0006] In order to achieve these objects and other advantages and in accordance with the purpose of the application, a coastal zone human activity interference remote sensing quantification and ecological loss and gain evaluation method is provided, which comprises the following steps: S1: Obtain multi-spectral remote sensing images of the target coastal zone at a first time point and a second time point, the first time point is earlier than the second time point, the time interval between the first time point and the second time point is between 5 years and 10 years, the spatial resolution of the multi-spectral remote sensing image is between 10 meters and 30 meters, the multi-spectral remote sensing image includes blue light band, green light band, red light band, near-infrared band and short wave infrared band, and tidal model data, digital elevation model data and geological background data are obtained at the same time; S2: Preprocess the multi-spectral remote sensing image, the preprocessing includes radiation calibration and atmospheric correction, to eliminate sensor error and atmospheric scattering effect; S3: Calculate the human activity interference index based on the preprocessed multi-spectral remote sensing image, the human activity interference index includes normalized building index and land use type change detection, wherein the normalized building index is calculated using the reflectance ratio of short wave infrared band to near infrared band, and the land use type change detection is realized by comparing the classification results of the first time point and the second time point, and the land use type change detection needs to be dynamically marked in the temporary flooding area based on the tidal model data and the digital elevation model data, and the area is excluded from the statistics in the change detection; S4: Quantify the degree of human activity interference, including calculating the area of the changed area and the intensity change value, the area of the changed area is based on the land use type change detection output, and the intensity change value is calculated based on the difference value of the normalized building index, and the intensity change value is corrected by the interference enhancement model to output the final intensity change value; The interference enhancement model generates a regional correction coefficient matrix based on the fragmentation interference index and the geological ecological data, and spatially weights the intensity change value; S5: The area of the changed area and the intensity change value are input into the ecological system service value model, and the ecological system service value model uses the unit area value equivalent factor method; S6: The ecological system service value model outputs the ecological loss and gain evaluation results including spatial distribution map and numerical report.

[0007] The present application is directed to the problem that the current human activity interference quantification means of the coastal zone excessively relies on field research, resulting in low efficiency, and the ecological loss and gain evaluation cannot be directly fused with remote sensing data, thereby affecting the reliability of the evaluation. The present application obtains two periods of multi-spectral remote sensing images (interval 5-10 years), eliminates environmental errors through radiation calibration and atmospheric correction; calculates the normalized building index (NDBI) based on the short-wave infrared and near-infrared bands, quantifies the interference area in combination with the land use classification change detection; inputs the change area and intensity value into the ecosystem service value model (using the unit area value equivalent factor method), and finally outputs the spatial distribution clear ecological loss and gain evaluation result. The present application realizes the whole-process remote sensing automatic processing, avoids the cost of field investigation; establishes a quantitative conversion mechanism from physical quantity change (NDBI difference) to ecological value loss, and improves the reliability of the evaluation.

[0008] Preferably, the pretreatment specifically comprises: performing radiation calibration on the original digital quantization value of the multi-spectral remote sensing image to output the radiation brightness data at the sensor entrance pupil; and performing atmospheric correction on the radiation brightness data to output the ground reflectance data.

[0009] Preferably, the radiation calibration comprises the following steps: a) reading the radiation calibration coefficient file provided by the sensor manufacturer, and converting the original digital quantization value into radiation brightness data based on the gain coefficient and the offset coefficient; b) performing sensor non-linear response compensation on the converted radiation brightness data: when the radiation brightness value is lower than 20% of the saturation radiation brightness value of the sensor, a low radiation zone gain compensation coefficient is used to improve the signal-to-noise ratio; when the radiation brightness value is higher than 50% of the saturation radiation brightness value of the sensor, a high radiation zone gain compensation coefficient is used to suppress the saturation effect; c) separately adding a specular reflection suppression process to the near-infrared band radiation brightness data, and using a polarization filtering algorithm in combination with the solar elevation angle, the sensor observation angle, the red band reflectance characteristics, and the sea surface wind speed data to correct the glint interference.

[0010] The radiation dynamic range of the coastal zone scene is too large, resulting in sensor non-linear response error. In the radiation calibration of the present application, partition compensation (low radiation zone to improve signal-to-noise ratio, high radiation zone to suppress saturation) is added; a polarization filtering algorithm is separately used for the near-infrared band in combination with the solar elevation angle, sea surface wind speed and other multi-source data to suppress the specular reflection. The radiation brightness data precision is improved, and the ground object reflectance distortion is reduced.

[0011] Preferably, the atmospheric correction comprises the following steps: a) obtaining aerosol optical thickness data and relative humidity data matching the imaging time of the multispectral remote sensing image, the aerosol optical thickness data derived from satellite remote sensing inversion products, the spatial resolution consistent with the multispectral remote sensing image, and the relative humidity data derived from meteorological reanalysis data, the time resolution less than or equal to 6 hours; b) based on the aerosol optical thickness data and the relative humidity data, calculating the aerosol scattering characteristic parameter at a wavelength of 550 nanometers, and using an aerosol scattering model for calculation, the model containing a mixing ratio parameter of sea salt type, dust type and artificial pollution type aerosols, and the mixing ratio parameter being set according to historical observation data of the target coastal zone region; c) converting the radiometrically calibrated radiance data into apparent reflectance data, and the conversion formula being: ; wherein, is the radiance value, d is the distance correction factor, is the extraterrestrial solar irradiance, is the solar zenith angle; d) using a radiative transfer model for atmospheric correction, the input parameters of the radiative transfer model including the aerosol scattering characteristic parameter, the sensor imaging geometry parameter, the ground elevation data and the apparent reflectance; e) in the radiative transfer model, when the relative humidity data exceeds 60%, introducing a humidity correction module, which uses Mie scattering theory to calculate the correction factor of the water vapor effect on the particle size expansion of the aerosol, to improve the accurate simulation of the model on the influence of humidity.

[0012] The high humidity environment and complex aerosol composition of the coastal zone lead to insufficient precision of the traditional atmospheric correction model. The present application uses an aerosol type parameterization model (mixed sea salt / dust / human pollution type), and introduces a relative humidity driven Mie scattering correction factor to quantify the particle expansion effect. The aerosol scattering effect is accurately corrected, and the reliability of the surface reflectance data is improved.

[0013] Preferably, the calculation of the human activity interference index comprises the following steps: a) performing land use classification on the surface reflectance data at the first time point and the second time point, respectively, using an object-oriented segmentation algorithm, setting the segmentation scale parameter to be between 10 and 30 pixels, and selecting the features to include the normalized building index, the normalized vegetation index, the improved normalized water body index and the mean value of the surface reflectance band; normalized building index ; normalized vegetation index ; improved normalized water body index ; wherein, is the surface reflectance of the blue light band, is the surface reflectance of the green light band, is the surface reflectance of the short-wave infrared light band, is the surface reflectance of the near-infrared light band, is the surface reflectance of the red light band; b) Based on the area in the classification result where the improved normalized water body index value is greater than the water body determination threshold, the tidal height data calculated by the tidal model are used to perform: When the tidal height value is greater than 1 meter of the local average sea level, the area with an elevation lower than the tidal height value and classified as a water body is marked as a temporary flooded area; The elevation data are derived from a digital elevation model, and the temporary flooded area does not participate in land use type change monitoring; c) In the non-flooded area and non-water body area, the building shadow area is determined based on the spectral characteristics of the vegetation covered area, if and , the building shadow area is determined; The average reflectance value of the same type of ground object in the adjacent non-shadow area is used to replace the reflectance of the shadow area; d) Based on the reflectance data after shadow compensation, the NDBI is recalculated; e) Land use type change detection is realized by comparing the classification results of the first time point and the second time point, and only the pixels in the non-temporary flooded area are counted for type conversion.

[0014] The intertidal zone water level change leads to land use misjudgment and building shadow reduction index precision. The present application dynamically marks the temporary flooded area based on the tidal model and digital elevation; in the vegetation covered area, the building shadow is identified through spectral criteria (near-infrared / red light band ratio <1.2 and red-blue band reflectance characteristic <0.2), and the reflectance is replaced with the average value of the same type of ground object in the adjacent non-shadow area. By reducing the influence of tides and shadows on classification, the extraction precision of the interference area is improved.

[0015] Preferably, the quantification of the degree of human activity interference comprises the following steps: a) Based on the land use type change monitoring result, the artificial land expansion area is extracted, and the total area of the changed area is calculated; b) Patch analysis is performed on the artificial land expansion area, the number of patches and the total length of patch edges are counted, the patch density index is calculated as the number of patches divided by the total area, and the edge density index is calculated as the total length of patch edges divided by the total area; c) The patch density index and the edge density index are fused according to the weights of 30% to 40% and 60% to 70% respectively to obtain the fragmentation interference index; d) The difference between the normalized building index at the first time point and the second time point is calculated as the basic intensity change value; e) input the base intensity change value into the interference enhancement model, and the interference enhancement model performs: When the fragmentation interference index is greater than the threshold value 0.5, the output enhanced intensity change value is the base intensity change value multiplied by the correction coefficient; If the fragmentation interference index does not exceed the threshold value of 0.5, the output enhanced intensity change value will remain consistent with the base intensity change value; f) output the enhanced intensity change value as the final intensity change value.

[0016] The fragmentation development interference is not well represented in the base intensity value. The present application generates a fragmentation interference index based on the fusion of patch density and edge density; when the index > 0.5, the base intensity change value is enhanced by a correction coefficient. The interference intensity quantification characterization of the intensive distributed development activities is strengthened.

[0017] Preferably, the interference enhancement model is constructed including the following steps: a) Obtain the geological background data of the target coastal zone, including the bedrock shoreline length proportion and the spatial distribution of the beach sediment type, and extract the ecological sensitive area spatial data, including the mangrove distribution boundary, the seagrass bed core area and the national coastal protection area range; b) Generate a regional correction coefficient matrix based on the geological background data and the ecological sensitive area data , the value range is 0.8~1.3, which is generated by coupling the fragmentation interference index and the geological and ecological data; c) The enhanced intensity change value is weighted according to the regional interference weight matrix: ; d) output as the enhanced intensity change value.

[0018] The single correction coefficient cannot adapt to the spatial heterogeneity of the coastal zone geology and ecology. The present application generates a spatial weight matrix by coupling the geological background data (bedrock shoreline proportion, sediment type) and the ecological sensitive area data (mangrove, seagrass bed), and weights the enhanced intensity change value by region. The spatial adaptive correction of the interference intensity quantification result is realized.

[0019] Preferably, the ecosystem service value model is constructed including the following steps: a) Establish a base value equivalent table, including the unit area annual service value coefficient of five types of land, i.e. mangrove, salt marsh, beach, shallow sea and artificial surface , wherein i =1: mangrove, 2: salt marsh, 3: beach, 4: shallow sea, 5: artificial surface; the coefficient is derived from the "Ecosystem Service Value Accounting Specification" issued by the national ecological protection department; b) Collecting dynamic parameters of coastal zone geographical environment: obtaining the following spatialized parameters of the target coastal zone area: intertidal zone width data; aquatic biodiversity index, fish species abundance calculated based on fishery resource survey data; storm surge historical occurrence frequency data; c) Generating a coastal zone spatial correction coefficient matrix: based on the parameters collected in step , dynamically correcting the value coefficient , and calculating the spatially corrected unit area annual service value coefficient ; ; When the intertidal zone width increases by more than 500 meters from the reference value, the ; When the storm surge occurrence frequency exceeds 0.5 times per year, the ; ; For land types or parameters that do not meet the above correction conditions ; d) Threshold value decay calculation driven by disturbance intensity: for artificial land expansion areas, based on the enhanced intensity change value, calculate the reduction factor through the disturbance decay module : When the final intensity change value is greater than the threshold value 0.3, ; When the final intensity change value is less than or equal to the threshold value 0.3,

[0020] e) Calculate the ecological loss and gain value: calculate the ecological service value change amount according to the following formula: ; Wherein, i is the land type index, i =1: mangrove, 2: salt marsh, 3: tidal flat, 4: shallow sea, 5: artificial land surface; is the area of the i th land type at the first time point; is the area of the i th land type at the second time point; the unit area annual service value coefficient of the i th land type after spatial correction; is the decay factor, and the ecological service value change caused by artificial land expansion is negative.

[0021] In order to solve the problem that the static value coefficient cannot reflect the spatial dynamics of the coastal zone and the nonlinear ecological influence of high-intensity development, the present application dynamically corrects the basic value coefficient based on the intertidal zone width, biodiversity index and storm surge frequency; for the artificial land expansion area, an interference attenuation factor (0.7-0.9 when 0.3 is taken) is applied according to the final intensity change value, so as to improve the spatio-temporal adaptability of the calculation of the change amount of the ecosystem service value.

[0022] Preferably, the interference attenuation module comprises the following steps: a) identifying the type of development activity: based on the patch characteristics of the artificial land expansion area, the patch characteristics including patch edge density and patch number, the type of development activity is determined; if the patch edge density is greater than 100 meters per hectare and the patch number is greater than 10, it is determined as port construction; if the patch number is 1 and the area is greater than 1 hectare, it is determined as land reclamation; the rest is determined as tourism facilities; b) obtaining the duration data of the interference, obtaining the time span data from the start of the development activity to the point of time for evaluating the ecological loss and gain, denoted as duration years T ; c) determining the type-based attenuation coefficient: according to the determined type of development activity, the corresponding basic attenuation coefficient is selected from the preset type-based attenuation coefficient table ; Land reclamation type: ; Port construction type: ; Tourism facilities type: ; d) calculating the time correction factor: based on duration years T , the time correction factor is calculated: ; e) determining the sensitive area correction factor: obtaining the spatial data of the ecological sensitive area, including the distribution boundary of mangrove forest, the core area of seagrass bed and the range of national coastal protection zone, and determining whether the artificial land expansion area overlaps with any ecological sensitive area: if there is spatial overlap, the sensitive area correction factor is set; if there is no spatial overlap, the ; f) calculating the final attenuation factor: based on , , the final output attenuation factor of the interference attenuation module is calculated : .

[0023] The present application is based on patch feature recognition to identify development type (reclamation / port / tourism); combined with duration (time correction factor) to calculate time correction factor; according to whether it overlaps with the ecological sensitive area to set the sensitive area correction factor; multiply the development type, duration and sensitive factor, and constrain the result in the interval [0.7, 0.9] to obtain the output value of the final attenuation factor. More accurate quantification of the cumulative ecological effects of high-intensity development. T

[0024] Other advantages, objects, and features of the present application will be apparent to those skilled in the art from the following specification. DETAILED DESCRIPTION

[0025] The present application will be further described in detail below, so that those skilled in the art can implement it according to the description.

[0026] It should be understood that the terms such as "have", "contain" and "include" used herein do not exclude the presence or addition of one or more other elements or combinations thereof.

[0027] Embodiment 1 A method for quantifying human activity interference in coastal zone and assessing ecological loss and gain by remote sensing, comprising the following steps: S1: Obtain multispectral remote sensing images of the target coastal zone at a first time point and a second time point, wherein the first time point is earlier than the second time point, the time interval is 5 to 10 years, and the spatial resolution of the multispectral remote sensing images is 10 to 30 meters. The multispectral remote sensing images include blue, green, red, near-infrared and short-wave infrared bands, and also include tidal model data, digital elevation model data and geological background data; S2: Preprocess the multispectral remote sensing images, which includes radiation calibration and atmospheric correction to eliminate sensor errors and atmospheric scattering effects; S3: Calculate the human activity interference index based on the preprocessed multispectral remote sensing images, which includes normalized building index and land use type change detection, wherein the normalized building index is calculated using the ratio of short-wave infrared band to near-infrared band reflectance, and the land use type change detection is realized by comparing the classification results of the first time point and the second time point, which needs to dynamically mark the temporary flooded area based on the tidal model data and the digital elevation model data, and exclude this area from the statistics in the change detection; ​S4: Quantify the degree of human activity interference, including calculating the area of the changed region and the intensity change value, the area of the changed region is based on the output of the land use type change detection, and the intensity change value is calculated based on the difference value of the normalized building index, and the intensity change value is finally obtained through the fine correction of the interference enhancement model; S5: The area of the changed region and the intensity change value are input into the ecosystem service value model, and the unit area value equivalent factor method is used in the ecosystem service value model; S6: The ecosystem service value model outputs the ecological loss and benefit evaluation results including spatial distribution map and numerical report.

[0028] This embodiment takes a certain bay area as an example. First, Landsat multispectral remote sensing images of the area in 2000 and 2010 are obtained, with a spatial resolution of 30 meters, including blue, green, red, near-infrared and short-wave infrared bands. At the same time, the same period tidal model data, 30-meter precision digital elevation model and geological background vector data are collected.

[0029] The traditional method needs to organize field teams to arrange sampling points in the coastal zone, manually record the building expansion range and draw sketches, which takes about 3 months. In addition, due to the positioning error (±5 meters) of handheld GPS and the deviation of paper registration, the actual spatial positioning accuracy of the field investigation results is 52.3 meters (calculated by 100 verification points, 95% confidence interval error 45-60 meters), and only the statistical average value of county scale can be presented. The present application directly uses satellite images, and the data acquisition time is not more than 1 week, which significantly improves the efficiency and real-time of data collection; through the whole process of automatic processing of remote sensing, the ecological loss and benefit spatial distribution map with 30m*30m grid precision is directly generated, which significantly improves the pertinence of management decision.

[0030] The pre-processing of the two images is performed: based on the sensor parameter file, the original digital value is converted into radiation brightness by completing radiation calibration; the surface reflectance data is output by using the special atmospheric correction model of the coastal zone. In the traditional method, the field investigators need to carry a spectrometer to measure the reflectivity on site, and only 2-3 point sampling can be completed per day; the present application generates the global surface reflectance data at one time through the processing of remote sensing physical quantities.

[0031] The human activity interference index is calculated: object-oriented classification is respectively performed on the reflectivity data of two periods, the segmentation scale is 20 pixels, five types of ground objects are divided based on NDBI, NDVI, MNDWI and band mean value; the artificial ground surface expansion area is extracted by comparing the classification results, and the change area is counted. At the same time, the NDBI difference value of two periods is calculated as the basic intensity change value. The traditional method relies on manual visual interpretation of the change spot, which is subjective and difficult to realize accurate quantification of intensity change; the present application realizes automatic quantification through classification comparison and NDBI difference value. The classification process can integrate a computing model based on a convolutional neural network to replace or assist the traditional object-oriented segmentation algorithm, improve the accuracy of ground object boundary recognition and type differentiation; the intensity change value calculation can introduce a random forest regression model to correct the nonlinear interference effect; the convolutional neural network model is applied to the land use classification task, and the multispectral features are automatically extracted through the trained network to replace the manually designed features, and the classification accuracy is significantly improved; the random forest model is used to fuse multiple source indicators (such as patch density and edge length) to predict the interference intensity and reduce the subjectivity of manually setting weights.

[0032] The change area and intensity value are input into the ecosystem service value model: the unit area value coefficient in the "Ecosystem Service Value Accounting Specification" is used, combined with the regional intertidal zone width and fish species abundance dynamic correction coefficient value, to output the ecological loss and gain spatial distribution map and value change amount report. The traditional method needs to manually input the expansion range drawn on site into the evaluation software, which is easy to introduce boundary error; the present application automatically converts the remote sensing physical quantity to the ecological parameter to directly generate the quantitative evaluation result.

[0033] The final result shows: from 2000 to 2010, the expansion of artificial ground surface led to the loss of ecosystem service value, and the spatial distribution map clearly presents the high-loss area superimposed by the core development area and the ecological sensitive area. Due to the limitation of the traditional method, the evaluation report can only present the average value at the county scale.

[0034] Further, the preprocessing specifically comprises: performing radiation calibration on original digital quantization values of the multispectral remote sensing image to output radiation brightness data at the sensor entrance pupil; performing atmospheric correction on the radiation brightness data to output surface reflectivity data.

[0035] After obtaining the Landsat images of 2000 and 2010, the original digital quantization values are first input into the radiation calibration module. By calling the coefficient file provided by the sensor manufacturer and using the gain and offset calculation method, the original pixel value is accurately converted into radiation brightness data. In the traditional method, different teams may use self-defined calibration formulas, resulting in inconsistent radiation brightness units of the same period images. For example, a study uses linear conversion and ignores the nonlinear characteristics of the sensor, resulting in systematic deviation in the radiation order of magnitude of the coastal dark pixels (such as deep water area) and bright pixels (such as sandy beach).

[0036] The radiance data is then input into an atmospheric correction module. In combination with aerosol optical depth products and relative humidity data provided by weather stations during satellite overpass, surface reflectance results are generated. In prior art systems, some processes choose to skip the atmospheric correction step and use apparent reflectance directly, or use a simplified model that ignores the specific effects of humidity on aerosols. For example, a coastal zone study used apparent reflectance to calculate vegetation indices, resulting in a 10-15% overestimation of water reflectance in near-shore waters, and misidentifying intertidal mud areas as vegetation degradation regions.

[0037] The present embodiment strictly limits the pre-processing output to two types of normalized data: sensor entrance pupil radiance (unit: W m -2 ·sr -1 ·μm -1 ) and surface reflectance (unitless). Traditional methods do not specify intermediate data standards, making it difficult to integrate results from different projects. For example, in a ten-year evaluation report for a certain region, data before 2005 used apparent reflectance, and data after 2005 used surface reflectance, resulting in an artificial discontinuity in the NDBI time series. The present invention unifies data standards, ensuring that physical quantities from two periods of images can be directly compared.

[0038] The final output of surface reflectance data is used for subsequent land use classification. Compared with traditional methods, standardized preprocessing eliminates the classification result deviation caused by differences in data sources. In a case of a prior art, using different atmospheric correction models to process the same image, the classification area of salt marshes differed by 8.2%; while the present invention forces the output to be standard, so that the classification result only responds to actual surface changes.

[0039] Further, the radiometric calibration includes the following steps: a) reading the radiometric calibration coefficient file provided by the sensor manufacturer, converting the original digital quantization value to radiance data based on the gain coefficient and offset coefficient; b) compensating for the non-linear response of the sensor to the converted radiance data: when the radiance value is less than 20% of the saturation radiance value of the sensor, use the low-radiance area gain compensation coefficient to improve the signal-to-noise ratio; when the radiance value is higher than 50% of the saturation radiance value of the sensor, use the high-radiance area gain compensation coefficient to suppress the saturation effect; c) separately add a specular reflection suppression process to the near-infrared band radiance data, using a polarization filtering algorithm combined with the solar elevation angle, sensor observation angle, red band reflectance characteristics, and sea surface wind speed data to correct the glint interference.

[0040] The traditional method directly uses the linear calibration coefficient provided by the sensor manufacturer, ignoring the characteristics of large dynamic range of coastal radiation. For example, a prior art uses the same conversion parameter for the deep water area (low radiation) and the salt field area (high radiation) of the same image, resulting in insufficient signal-to-noise ratio in the deep water area and saturation effect in the salt field area. The specific performance is as follows: the radiation brightness difference between the mangrove shadow and the turbid water body is less than 5%, and the confusion rate is 23% when classified; the salt field edge loses texture details due to signal truncation, and is misjudged as building land.

[0041] The present scheme implements partition compensation: Low radiation area processing: identify pixels with radiation values below 20% of the saturation value (such as deep water area, mangrove canopy shadow), use 1.15 times gain compensation coefficient to improve signal-to-noise ratio, and expand the reflectivity difference between shadow water body and vegetation to a distinguishable range.

[0042] High radiation area processing: for pixels exceeding 50% of the saturation value (such as salt field, bare beach), apply 0.9 times gain suppression coefficient to prevent signal overflow and preserve salt field crystal texture features.

[0043] Near-infrared flare suppression: under the condition of 5.2 meters / second southeast wind in the gulf, combined with the solar elevation angle of 58 degrees and the sensor observation angle of 8 degrees, a polarization filtering algorithm is used to correct the specular reflection. Compared with the prior art which only uses red light band threshold method, it has the disadvantage of eliminating wind speed influence, and the flare residual problem often causes the nearshore water body to be incorrectly classified as bare land.

[0044] After improvement, the radiation response curve of deep water area and salt field is closer to the actual measured spectrum, and the misclassification rate of mangrove boundary is significantly reduced. The traditional method often needs manual intervention to correct the classification results, which takes about two weeks; the present invention directly generates reliable data for subsequent processes through adaptive compensation.

[0045] Further, the atmospheric correction includes the following steps: a) Obtain aerosol optical depth data and relative humidity data matched with the imaging time of the multi-spectral remote sensing image, wherein the aerosol optical depth data is derived from high-precision satellite remote sensing inversion products, and the spatial resolution is completely matched with the multi-spectral remote sensing image; the relative humidity data is derived from the latest meteorological reanalysis data, and the time resolution is accurate to within 6 hours; b) Based on the aerosol optical depth data and the relative humidity data, calculate the aerosol scattering characteristic parameter at a wavelength of 550 nanometers, the calculation uses an aerosol scattering model, which includes the mixing ratio parameters of sea salt type, dust type and artificial pollution type aerosols, and the mixing ratio parameters are set according to the historical observation data of the target coastal area; c) Convert the radiation brightness data after radiation calibration to apparent reflectance data, and the conversion formula is: ; wherein, is the radiance value, d is the earth-sun distance correction factor, is the extraterrestrial solar irradiance, is the solar zenith angle; d) performing atmospheric correction using a radiative transfer model, the input parameters of which include aerosol scattering property parameters, sensor imaging geometry parameters, ground elevation data, and apparent reflectance; e) in the radiative transfer model, when the relative humidity data exceeds 60%, a humidity correction module is introduced, which uses Mie scattering theory to calculate the correction factor of the water vapor expansion effect on aerosol particle size, to ensure the accuracy of the model.

[0046] The prior art generally uses a fixed type aerosol model, such as a study that uniformly selects marine aerosol parameters for coastal images. This simplified process does not take into account the actual situation of bays affected by industrial emissions, resulting in systematic underestimation of surface reflectivity in urban nearshore areas. Specifically, the proportion of artificial pollution components in aerosols over industrial areas exceeds 40%, but traditional models still treat them as pure sea salt, resulting in a building roof reflectivity error of 8% to 12%, which further affects the recognition accuracy of artificial ground expansion.

[0047] The present scheme implements fine correction: first, obtain the aerosol optical depth product on the day of satellite overpass, with a spatial resolution matching 30-meter image data. Simultaneously integrate the 6-hourly relative humidity field from meteorological reanalysis data to identify cloud-free areas with humidity higher than 75%. Based on historical observation data, set the aerosol mixing ratio - in the northern industrial area of the bay, use a combination of 20% sea salt type, 30% dust type, and 50% artificial pollution type; in the southern tourist area, increase the sea salt type proportion to 60%. When the humidity exceeds the 60% threshold, start the Mie scattering correction module to quantify the water vapor expansion effect. Compared to a prior art that ignores the impact of humidity changes and directly applies dry season correction results to rainy season images, the reflectivity of mangrove areas abnormally rises, which is then incorrectly determined as vegetation degradation.

[0048] The corrected data significantly improves the classification reliability. The traditional method produces 15% of false building expansion patches in the port area due to inaccurate aerosol modeling; the present invention adjusts the dynamic parameters, making the reflectivity spatial gradient of the industrial area and the ecological protection area more consistent with the field survey results, providing reliable input for interference quantification.

[0049] Further, calculating the human activity interference index includes the following steps: a) Land use classification is performed on the ground reflectance data of the first and second time points respectively, using an object-oriented segmentation algorithm with a segmentation scale parameter set between 10 and 30 pixels, and features including normalized building index, normalized vegetation index, improved normalized water index, and mean value of ground reflectance bands; normalized building index ; normalized vegetation index ; improved normalized water index ; wherein, is the blue band ground reflectance, is the green band ground reflectance, is the short-wave infrared band ground reflectance, is the near-infrared band ground reflectance, is the red band ground reflectance; b) Based on the areas in the classification results where the improved normalized water index value is greater than the water body determination threshold, and combined with the tide height data calculated by the tide model, the following is performed: When the tide height value is greater than 1 meter of the local average sea level, the area with an elevation lower than the tide height value and classified as water body is marked as a temporary flooded area; The elevation data is derived from a digital elevation model, and the temporary flooded area does not participate in land use type change monitoring; c) In the non-flooded area and non-water body area, the building shadow area is determined based on the spectral characteristics of the vegetation covered area. If the vegetation area is and , it is determined as a building shadow area; The average reflectance value of the same type of ground object in the adjacent non-shadow area is used to replace the reflectance of the shadow area; d) Based on the reflectance data after shadow compensation, the NDBI is recalculated; e) Land use type change detection is achieved by comparing the classification results of the first and second time points, and only the pixels in the non-temporary flooded area that have undergone type conversion are counted.

[0050] Traditional land use classification methods do not fully consider the influence of tidal dynamics, resulting in significant errors in the classification results. For example, a technique classifies the exposed beach at low tide as bare land, while at high tide it is determined as water body. This classification method produces up to 18.7% of false change patches in two years. At the same time, the shadow of the port building in the vegetation area is misjudged as a vegetation degradation area due to reflectance distortion, resulting in an increase in the interference range.

[0051] This scheme implements dynamic masking and shadow compensation: Tidal inundation area marking: combine the tide height data at satellite overpass time with digital elevation model, mark the water body with elevation lower than 1.2m above mean sea level as temporary inundation area when the tide height exceeds 1.2m. For example, 650 hectares of tidal creek are excluded from the changeable area to avoid the misconversion between bare land and water body at high tide.

[0052] Building shadow identification and repair: detect spectral anomalies in vegetation-covered areas - areas with a near-infrared / red band ratio less than 1.15 and a red-blue reflectance difference less than 0.19 are determined to be shadows. For the mangrove area behind the port, we fill in the average reflectance value of the adjacent non-shadow area of the mangrove to repair the detected shadow patches. Compared with a prior art method that only uses NDVI thresholding, it fails to distinguish between shadows and real vegetation degradation, resulting in 42 hectares of false degradation patches around the wharf.

[0053] After optimization, the accuracy of intertidal zone area change detection is significantly improved. Traditional methods require manual correction when dealing with tidal misjudgment areas, which takes up to three weeks; this invention uses automated dynamic masking to directly output reliable change area data for disturbance quantification.

[0054] Further, the quantification of the degree of human activity disturbance includes the following steps: a) Based on the land use type change monitoring results, extract the artificial land expansion area and calculate the total area of the change area; b) Perform patch analysis on the artificial land expansion area, count the number of patches and the total length of patch edges, and calculate the patch density index by dividing the number of patches by the total area, and the edge density index by dividing the total length of patch edges by the total area; c) Fuse the patch density index and the edge density index into a fragmentation disturbance index with a weight of 30% to 40% and 60% to 70%; d) Calculate the difference between the normalized building index at the first time point and the second time point as the basic intensity change value; e) Input the basic intensity change value into the disturbance enhancement model, and the disturbance enhancement model performs: When the fragmentation disturbance index is greater than the threshold value 0.5, the enhanced intensity change value is the basic intensity change value multiplied by the correction coefficient; When the fragmentation disturbance index is less than or equal to the threshold value 0.5, the enhanced intensity change value is equal to the basic intensity change value; f) Output the enhanced intensity change value as the final intensity change value.

[0055] The prior art only uses two-period NDBI difference to represent the development intensity, resulting in the fragmentation construction being underestimated. In a certain existing case, the scattered layout of holiday villa groups (single area 0.1-0.3 hectares) and the centralized port area (5 hectares) obtain the same intensity value, which cannot reflect the difference of ecological impact. The traditional method generally assigns a strength value of 0.4 to a villa group of 20 hectares and a port area of 15 hectares, but the actual monitoring results show that the villa group has a much greater impact on bird habitats than the port area.

[0056] The present scheme implements fragmentation index correction: Patch feature extraction: morphological analysis is performed on the artificial surface expansion area, and 128 patches are counted (121 of which are villa groups), with a total edge length of 18.6 kilometers. The patch density index is 6.4 per hectare, and the edge density index is 93 meters per hectare.

[0057] Index fusion and correction: fusion with 35% and 65% weight gets the fragmentation disturbance index 0.72, which exceeds the threshold value 0.5. Trigger the correction mechanism, multiply the basic NDBI difference 0.25 by 1.28 times the enhancement coefficient, output the final intensity value 0.32. Compared with a certain study in the prior art, which uses a fixed area weight method to set the intensity of the port area to 0.42, while the villa group is only 0.19, this setting is contrary to the disturbance level assessed by ecologists in the field.

[0058] The corrected intensity value is more consistent with the observation of ecological effects. In the quantitative results of the traditional method, the ecological disturbance intensity value of the port area is 2.2 times that of the villa group. However, after enhancing the representation of fragmented areas by the present invention, we found that the ecological disturbance intensity value of the villa group actually exceeded the port area by 12%, which more accurately reflects the cumulative damage effect of distributed development on habitat connectivity.

[0059] In addition, the construction process of the disturbance enhancement model includes the following steps: a) Obtain the geological background data of the target coastal area, including the proportion of bedrock coastline length and the spatial distribution of tidal flat sediments, and extract the spatial data of ecologically sensitive areas, including the boundaries of mangrove distribution, seagrass bed core area, and national coastal protection area; b) Generate a regional correction coefficient matrix based on geological background data and ecologically sensitive area data , with a value range of 0.8-1.3, which is generated by coupling the fragmentation disturbance index and geological and ecological data; c) Weight the enhanced intensity change value according to the regional disturbance weight matrix:

[0060] d) Output as the enhanced intensity change value.

[0061] The prior art uses a single correction coefficient to process the entire data, ignoring the spatial heterogeneity of geology and ecology. A certain existing case uniformly applies a 1.25 times enhancement coefficient to the bedrock shore industrial area and the mudflat tourist area, which raises two problems: the bedrock area, due to its geological stability, actually suffers less interference, resulting in an 18% higher corrected intensity value; on the contrary, the mangrove mudflat area, due to its high ecological sensitivity, the correction effort is insufficient, thus the ecological impact brought by the construction of a 30-hectare resort is severely underestimated.

[0062] The present scheme implements spatial weighted correction: Regional weight matrix construction: Extract the geological data of the 82% bedrock shoreline on the north shore, and generate the weight coefficient in combination with the distribution boundary of the mangrove forest. The bedrock area is given a weight of 0.85 due to its strong anti-interference, the core area of the mangrove forest on the south shore is given a weight of 1.25 due to its high ecological sensitivity, and the ordinary mudflat area is set to the baseline value of 1.0.

[0063] Intensity value dynamic modulation: The output villa group enhancement intensity value 0.32 is input into this module. After 0.85 weight adjustment, the resort intensity value (original value 0.32) in the bedrock area on the north shore drops to 0.27; while the viewing platform intensity value (original value 0.28) in the mangrove forest area on the south shore rises to 0.35 under the support of 1.25 weight. Compared with the fixed coefficient method of a certain existing technology, which sets the viewing platform intensity value to 0.35 and the resort to 0.40, it is contradictory to the evaluation conclusion of ecologists - field monitoring shows that the viewing platform causes the mortality rate of mangrove seedlings to reach 23%, which is much higher than the 9% of the resort.

[0064] The correction result is more in line with the spatial law of the coastal zone. The traditional method quantifies the north shore intensity as 1.14 times that of the south shore; the present invention, through geological-ecological coupled weighting, makes the intensity value of the sensitive area on the south shore exceed that of the north shore by 29%, accurately reflecting the doubled ecological risk of development in the mangrove forest area. The final intensity distribution map is verified by the marine management department and directly used for ecological restoration priority delineation.

[0065] This matrix is calculated by coupling the fragmentation interference index with geological and ecological data: Calculate the geological stability factor G = bedrock shoreline length ratio × 0.01 + sediment resistance coefficient (sandy beach = 0.8, muddy beach = 1.2, mixed type = 1.0); Extract the ecological sensitivity factor E: core area of mangrove forest or seagrass bed = 1.3, buffer zone of protected area = 1.1, other areas = 1.0; Based on the fragmentation interference index F index Calculate the dynamic modulation term ; Output .

[0066] The ecosystem service value model construction includes the following steps: a) Establish a basic value equivalent table, including the annual service value coefficient of unit area of five types of land: mangrove, salt marsh, beach, shallow sea and artificial surface , wherein i =1: mangrove, 2: salt marsh, 3: beach, 4: shallow sea, 5: artificial surface; b) Collect the dynamic parameters of the coastal zone geographical environment: obtain the following spatialized parameters of the target coastal zone area: tidal zone width data; Aquatic biodiversity index, fish species abundance calculated based on fishery resource survey data; Storm surge historical occurrence frequency data; c) Generate a coastal zone spatial correction coefficient matrix: based on the parameters collected in step, dynamically correct the value coefficient , and calculate the spatially corrected annual service value coefficient of unit area ; When the intertidal zone width is greater than 500 meters than the reference value, the ; When the fish species abundance increases by 10 species, the ; When the storm surge occurrence frequency is more than 0.5 times per year, the ; ; For the land or parameter that does not meet the above correction conditions ; d) Threshold attenuation calculation driven by disturbance intensity: for the artificial surface expansion area, based on the enhanced intensity change value, calculate the reduction factor through the disturbance attenuation module : When the final intensity change value is greater than the threshold value 0.3, ; When the final intensity change value is less than or equal to the threshold value 0.3,

[0067] e) Calculate the ecological loss and gain value: calculate the change amount of ecological service value according to the following formula: ; wherein, i is the land index, i =1: mangrove, 2: salt marsh, 3: beach, 4: shallow sea, 5: artificial surface; is the first time pointi Class area; The second time point is the i Class area; for the spatial correction is the i Class unit area annual service value coefficient; is the attenuation factor, the ecological service value change caused by artificial land expansion is negative.

[0068] The prior art directly applies the static coefficient of the "ecosystem service value accounting specification", ignoring the geographical dynamic characteristics of the coastal zone. A certain prior art uses fixed value coefficients for the entire region from 2000 to 2010: 2.8 million yuan / acre / year for the beach, and 1.2 million yuan / acre / year for the shallow sea. However, actual monitoring shows that the intertidal zone on the north coast has expanded by 600 meters due to siltation, resulting in improved quality of shellfish habitats; while the south coast has seen an increase in storm tide frequency to an average of 0.7 times per year, significantly improving the protective value of mangroves. The traditional method does not take these changes into account, resulting in the north coast ecological gain not being reflected and the south coast disaster prevention value being underestimated by about 30%.

[0069] The present scheme implements dynamic correction: Spatial parameter integration: based on hydrological survey data, it is confirmed that the width of the north coast intertidal zone increases by 520 meters, triggering the beach value coefficient to rise to 3.36 million yuan / acre / year; the storm tide frequency on the south coast reaches 0.7 times / year, and the coefficients of mangroves and salt marshes are simultaneously increased to 1.2 times the baseline value.

[0070] Interference attenuation mechanism: for the 15 hectares of artificial reclamation area on the south coast, the final intensity change value of 0.41 exceeds the threshold value of 0.3, and a 0.75 attenuation factor is applied. The north coast has 8 hectares of tourism facilities with an intensity value of 0.29, which does not reach the threshold value, and the attenuation factor is 1.0. In contrast, a certain prior art uniformly applies an attenuation factor of 0.8 to all artificial land surfaces, which underestimates the ecological loss of the reclamation area by 1.6 times that of the tourism area, while the ecological model shows that the actual loss ratio should be 2.3 times.

[0071] When calculating the value, first calculate the ecological value of the two periods according to the corrected coefficients: the north coast beach in 2010 has an area expansion and a coefficient increase, with a value of 8.6 million yuan more than in 2000; although the south coast mangrove area decreases by 5 hectares, the coefficient increase partially offsets the loss. Finally, combine the attenuation factor to output the net loss and gain: the south coast reclamation area has a high-intensity development superimposed on a sensitive location, with an actual loss of 1.8 times the value estimated by the traditional method.

[0072] The results have been confirmed by the local marine bureau: the traditional static assessment suggests that the north coast should be prioritized for restoration, while the present invention reveals that high-intensity reclamation on the south coast is the core contradiction. After the redistribution of restoration funds, the recovery efficiency of the south coast seagrass bed has significantly improved, with a 40% increase, further proving the scientificity and effectiveness of the dynamic correction mechanism in the ecological restoration of seagrass beds.

[0073] Further, the interference attenuation module is constructed comprising the following steps: a) Identifying the development activity type: Based on the patch characteristics of artificial land expansion area, the patch characteristics include patch edge density and patch number, the development activity type is determined; If the patch edge density > 100 meters / ha and the patch number > 10, it is determined as port construction; If the patch number = 1 and the area > 1 ha, it is determined as land reclamation; The rest is determined as tourism facilities; b) Obtain the interference duration data, obtain the time span data from the start of the development activity to the ecological loss and benefit evaluation point, recorded as duration years T ; c) Determine the type-based attenuation coefficient: According to the determined development activity type, select the corresponding basic attenuation coefficient from the preset type-based attenuation coefficient table ; Land reclamation type: ; Port construction type: ; Tourism facilities type: ; d) Calculate the time correction factor: Based on duration years T , calculate the time correction factor : ; e) Determine the sensitive area correction factor: Obtain the ecological sensitive area spatial data, including the distribution boundary of mangrove forest, the core area of seagrass bed and the scope of national coastal protection area, and judge whether the artificial land expansion area overlaps with any ecological sensitive area: If there is spatial overlap, set the sensitive area correction factor ; If there is no spatial overlap, set ; f) Calculate the final attenuation factor: Based on , , Calculate the final output of the interference attenuation module, the attenuation factor of the applicability of ecological loss and benefit value : .

[0074] The prior art only sets a fixed attenuation coefficient according to the development area, such as a study that uniformly adopts a 0.7 attenuation value for reclamation projects. This simplified processing method fails to fully consider the cumulative effects of different development stages and the protection priority of ecological sensitive areas, thus leading to the following two misjudgments: one is that a new port in operation for 3 years and an abandoned old wharf for 10 years obtain the same attenuation value, ignoring the natural recovery process; the other is that a resort occupying a core area of mangrove and a peripheral facility are treated equally, underestimating the loss of a biodiversity hotspot.

[0075] The present technology implements a differentiated correction process: Development type identification: Through analysis of the reclamation area patches, when the area of a single patch reaches 18 hectares and there is no secondary patch, it is determined to be a reclamation type, and the basic attenuation coefficient is set to 0.72. In contrast, a prior art mistakenly judges the adjacent 12-hectare tourist facility (actually composed of 6 scattered patches) as a reclamation type, and incorrectly assigns an attenuation value of 0.73 (which should be 0.88 for tourist type).

[0076] Duration correction: The project has lasted for 7 years from construction to evaluation, triggering the calculation of the time correction factor: since 5 < T ≤ 10, according to the formula 0.9 + 0.01 × (7-5) = 0.92. The traditional method does not include the time dimension and treats new ports and historical wharfs equally.

[0077] Sensitive area superposition detection: Spatial analysis shows that the reclamation area overlaps with the core protection area of mangrove by 9.2 hectares, and the sensitive area correction factor is set to 0.87. The prior art only roughly judges according to the buffer zone of the protected area boundary and fails to accurately identify the actual overlapping area.

[0078] Final attenuation factor synthesis: The intermediate value is 0.577 by multiplying the basic coefficient 0.72, the time factor 0.92, and the sensitive factor 0.87, and the final attenuation factor is 0.7 after being constrained in the interval [0.7, 0.9]. The traditional method directly takes the reclamation type standard value 0.7, although the result is the same, but it lacks the quantitative response to the 7-year recovery period and the core area overlap.

[0079] Re-estimation of ecological value loss: The original loss value calculated by the module is 38 million yuan, and after correction by the attenuation factor of 0.7, the loss value becomes 26.6 million yuan. The traditional method, which does not distinguish the development attributes, uses the tourist facility attenuation factor of 0.85 for the same area, underestimating the loss to 32.3 million yuan. The historical compensation records of the Marine Bureau show that the actual ecological restoration investment in this area is 29.1 million yuan, which is consistent with the evaluation results of the present invention. Practice has proved that the multi-dimensional attenuation mechanism significantly improves the quantitative precision of damage, providing a reliable basis for the formulation of ecological compensation standards.

[0080] Example 2 The device comprises six core modules connected to each other, and realizes full-process automatic processing through software and hardware cooperation. The data acquisition module is equipped with a satellite sensor including a photovoltaic conversion unit, which is used to acquire two-period multi-spectral remote sensing images (with a spatial resolution of 10-30 meters, including blue, green, red, near-infrared and short-wave infrared bands) of the target coastal zone region at intervals of 5-10 years, and simultaneously collect tidal model data, digital elevation model and geological background data, so as to ensure the integrity and timeliness of the source data.

[0081] After receiving the original image, the pretreatment module first executes radiation calibration through GPU parallel computing architecture: according to the nonlinear response characteristics of the sensor, a gain compensation coefficient of 1.15 times is applied to the low radiation area to enhance the signal-to-noise ratio, and a suppression coefficient of 0.9 times is applied to the high radiation area to avoid signal saturation, in addition, a polarization filtering algorithm is implemented on the near-infrared band to reduce the influence of sea surface glare. Then, the real-time aerosol optical depth and humidity data are integrated, and the atmospheric correction is completed by using the radiation transfer model, and the accurate ground reflectance data set is output, and the processing efficiency is improved by 15 times compared with the traditional CPU scheme.

[0082] The interference index calculation module is built-in U-Net deep learning model, which implements object-oriented segmentation strategy for the two-period reflectance data after pretreatment: 20 pixels are used as segmentation reference, NDBI, NDVI, MNDWI and band mean features are fused to realize accurate classification of land use; at the same time, with the help of digital elevation data and tidal model, the temporary flooded area is dynamically identified, and the spectral criterion (i.e. the ratio of near-infrared / red light band is less than 1.2 and the red / blue reflectivity characteristic is lower than 0.2) is used to intelligently identify the building shadow area, so as to replace the abnormal reflectivity with the mean value of the same ground object in the adjacent non-shadow area, and finally output the NDBI and land use change map of non-flooded area after shadow compensation.

[0083] The interference quantification module extracts the artificial ground expansion area, counts the change area and fuses the patch density index (weight 35%) and edge density index (weight 65%) to generate the fragmentation interference index. The interference quantification module can combine the random forest regression model to predict the intensity change value according to the patch characteristics and environmental parameters, and improve the quantification reliability. When the index is greater than 0.5, the interference enhancement model is started, which generates a regional correction coefficient matrix based on the geological background data (such as bedrock shoreline proportion) and the distribution of ecological sensitive area (such as mangrove boundary), and performs spatial weighting on the NDBI difference value, and outputs the final intensity change value.

[0084] The ecological assessment module uses the unit area value equivalent factor method to construct a dynamic assessment model: the model adjusts and optimizes the ecological value coefficient in real time according to key parameters such as intertidal zone width and fish species abundance, and calculates the attenuation factor of the development type (reclamation / port / tourism), duration and overlap state of the ecological sensitive area through the interference attenuation module. During the assessment process, a random forest-based calculation model can be introduced to optimize the value coefficient correction and attenuation factor calculation, and to enhance the adaptability of the model to nonlinear ecological responses. The assessment results are uploaded to the government coastal belt supervision platform in real time through a special interface, triggering the automatic accounting process of the ecological compensation fund.

[0085] The result output module can generate ecological loss and gain spatial distribution maps and numerical reports with a 30-meter grid precision, supporting three-dimensional dynamic risk simulation. By superimposing the geological sensitive area boundary and the tidal inundation model, the system can deduce long-term ecological losses under the scenario of sea level rise in time series, providing visual decision support for restoration projects. For example, according to the new mode of mangrove ecological reconstruction effect prediction research, the cumulative impact of reclamation areas on mangrove larval survival rate in the next 30 years can be simulated, and high-risk areas where the loss increase may exceed 40% are warned.

[0086] Although the embodiments of the present application have been disclosed as above, they are not limited to the use listed in the specification and embodiments, and can be fully applied to various fields suitable for the present application, and additional modifications can be easily implemented by those skilled in the art, and therefore the present application is not limited to specific details without departing from the general concept defined by the claims and equivalent scope.

Claims

1. A method for remote sensing quantification of human activity disturbance in coastal zones and assessment of ecological damage and benefit, characterized in that, The method comprises the following steps: S1: Obtain multispectral remote sensing images of a target coastal zone at a first time point and a second time point, the first time point being earlier than the second time point, the time interval between the first time point and the second time point being between 5 years and 10 years, the spatial resolution of the multispectral remote sensing images being between 10 meters and 30 meters, the multispectral remote sensing images comprising a blue light waveband, a green light waveband, a red light waveband, a near-infrared waveband and a short-wave infrared waveband, and simultaneously obtaining tide model data, digital elevation model data and geological background data; S2: Preprocess the multispectral remote sensing images, the preprocessing comprising radiation calibration and atmospheric correction to eliminate sensor errors and atmospheric scattering effects; S3: Calculate a human activity interference index based on the preprocessed multispectral remote sensing images, the human activity interference index comprising a normalized building index and land use type change detection, wherein the normalized building index is calculated using the ratio of short-wave infrared waveband reflectivity to near-infrared waveband reflectivity, and the land use type change detection is achieved by comparing the classification results at the first time point and the second time point, and the land use type change detection needs to be dynamically marked in a temporary flooding area based on the tide model data and the digital elevation model data, and the area is excluded from statistics in the change detection; S4: Quantify the degree of human activity interference, including calculating the area of the changed region and the intensity change value, the area of the changed region being based on the land use type change detection output, and the intensity change value being calculated based on the difference of the normalized building index, and the final intensity change value being output after the intensity change value is corrected by an interference enhancement model; The interference enhancement model generates a regional correction coefficient matrix based on a fragmentation interference index and geological and ecological data, and spatially weights the intensity change value; S5: The area of the changed region and the intensity change value are input into an ecosystem service value model, and the ecosystem service value model uses the unit area value equivalent factor method; S6: The ecosystem service value model outputs an ecological loss and benefit evaluation result comprising a spatial distribution map and a numerical report.

2. The method according to claim 1, wherein, The preprocessing specifically comprises: performing radiation calibration on the original digital quantization values of the multispectral remote sensing images to output radiation brightness data at the entrance pupil of the sensor; and performing atmospheric correction on the radiation brightness data to output surface reflectivity data. 3.The method according to claim 1, characterized in that, The radiation calibration comprises the following steps: a) Read the radiation calibration coefficient file provided by the sensor manufacturer, and convert the original digital quantization values into radiation brightness data based on the gain coefficient and the offset coefficient; b) Perform sensor non-linear response compensation on the converted radiation brightness data: when the radiation brightness value is lower than 20% of the saturation radiation brightness value of the sensor, use a low radiation area gain compensation coefficient to improve the signal-to-noise ratio; when the radiation brightness value is higher than 50% of the saturation radiation brightness value of the sensor, use a high radiation area gain compensation coefficient to suppress the saturation effect; c) Independently add a specular reflection suppression process to the near-infrared waveband radiation brightness data, and use a polarization filtering algorithm combined with the solar elevation angle, the sensor observation angle, the red light waveband reflectivity characteristics and the sea surface wind speed data to correct the glare interference.

4. The method according to claim 1, wherein, The atmospheric correction comprises the following steps: a) Obtain aerosol optical thickness data and relative humidity data matched with the imaging time of the multispectral remote sensing image, the aerosol optical thickness data is derived from satellite remote sensing inversion products, the spatial resolution is consistent with the multispectral remote sensing image, and the relative humidity data is derived from meteorological reanalysis data, the time resolution is less than or equal to 6 hours; b) Based on the aerosol optical thickness data and the relative humidity data, calculate the aerosol scattering characteristic parameter at a wavelength of 550 nanometers, and calculate using an aerosol scattering model, which includes a mixing ratio parameter of sea salt type, dust type and artificial pollution type aerosols, the mixing ratio parameter is set according to the historical observation data of the target coastal zone area; c) Convert the radiometrically calibrated radiance data to apparent reflectance data, and the conversion formula is: ; wherein is a radiance value, d is a solar distance correction factor, is the solar irradiance outside the atmosphere, is the solar zenith angle; d) Perform atmospheric correction using a radiative transfer model, the input parameters of the radiative transfer model include aerosol scattering characteristic parameters, sensor imaging geometry parameters, ground elevation data and apparent reflectance; e) In the radiative transfer model, when the relative humidity data is greater than 60%, introduce a humidity correction module, which uses Mie scattering theory to calculate the correction factor of the water vapor on the size expansion effect of aerosol particles, to improve the accurate simulation of the model on the influence of humidity.

5. The method according to claim 1, wherein, The calculation of the human activity interference index includes the following steps: a) Perform land use classification on the ground reflectance data at the first time point and the second time point respectively, use object-oriented segmentation algorithm during classification, set the segmentation scale parameter between 10 to 30 pixels, and select features including normalized building index, normalized vegetation index, improved normalized water index and mean value of ground reflectance band; Normalized building index ; Normalized difference vegetation index ; Improved normalized difference water index ; wherein, is the surface reflectance in the blue light band, is the surface reflectance in the green light band, is the surface reflectance in the short wave infrared band, is the surface reflectance in the near infrared band, is the surface reflectance in the red light band; b) Based on the area where the improved normalized water index value in the classification result is greater than the water body judgment threshold, combined with the tidal height data calculated by the tide model, execute: When the tidal height value is greater than 1 meter of the local average sea level, mark the area with elevation lower than the tidal height value and classified as water body as temporary flooded area; The elevation data is derived from digital elevation model, and the temporary flooded area does not participate in land use type change monitoring; c) in non-submerged area and non-water area, based on the spectral characteristics of the vegetation covered area, determine the building shadow area, if the vegetation area is and , determine the building shadow area; Use the average reflectance value of the same type of ground object in the adjacent non-shaded area to replace the reflectance of the shaded area; d) Recalculate NDBI based on the reflectance data after shadow compensation; e) Land use type change detection is realized by comparing the classification results at the first time point and the second time point, and only the pixels in the non-temporary flooded area are counted. 6.The method according to claim 1, characterized in that, The quantification of the degree of human activity interference includes the following steps: a) Based on the land use type change monitoring result, extract the artificial land expansion area, and calculate the total area of the changed area; b) Perform patch analysis on the artificial land expansion area, count the number of patches and the total length of patch edges, calculate the patch density index as the number of patches divided by the total area, and calculate the edge density index as the total length of patch edges divided by the total area; c) Fuse the patch density index and the edge density index according to the weights of 30% to 40% and 60% to 70% respectively to obtain the fragmentation interference index; d) Calculate the difference value of the normalized building index at the first time point and the second time point as the basic intensity change value; e) inputting the base intensity change value into the interference enhancement model, the interference enhancement model performing: when the fragmentation interference index is greater than the threshold value 0.5, outputting an enhanced intensity change value as the base intensity change value multiplied by a correction coefficient; if the fragmentation interference index does not exceed the threshold value of 0.5, the output enhanced intensity change value will remain consistent with the base intensity change value; f) outputting the enhanced intensity change value as the final intensity change value.

7. The method according to claim 6, wherein, The interference enhancement model construction includes the following steps: a) obtaining geological background data of the target coastal zone, including bedrock coastline length proportion and spatial distribution of tidal flat sediment type, and extracting ecological sensitive area spatial data, including mangrove distribution boundary, seagrass bed core area and national coastal protection area range; b) generating a regional correction coefficient matrix based on the geological background data and the ecological sensitive area data , the value range is 0.8~1.3, the matrix is generated by coupling the fragmentation interference index and the geological ecological data; c) the enhanced intensity change values are weighted by a region-wise interference weight matrix weighted by a region-wise interference weight matrix: ； d) output as an enhanced intensity change value. 8.The method according to claim 1, characterized in that, The ecosystem service value model construction includes the following steps: a) Establishing the basic value equivalent table, including the unit area annual service value coefficient of five types of land, i.e. mangrove forest, salt marsh, mudflat, shallow sea and artificial land surface wherein i = 1: mangrove forest, 2: salt marsh, 3: mudflat, 4: shallow sea, 5: artificial land surface; b) collecting coastal geographic environment dynamic parameters: obtaining the following spatialized parameters of the target coastal zone: tidal zone width data; water-borne biodiversity index, fish species abundance calculated based on fishery resource survey data; storm surge historical occurrence frequency data; c) generating a coastal zone space correction coefficient matrix: based on the parameters collected in step a), the value coefficient is dynamically corrected, and the unit area annual service value coefficient after space correction is calculated ; when the tidal zone width increases by more than 500 meters compared to the reference value, the tidal flat land class ； when the fish species abundance increases by 10 species, the shallow sea land class ； when the storm surge occurrence frequency exceeds 0.5 times per year, the mangrove land class and the salt marsh land class ； ； For the land class or parameter that does not meet the above correction condition ; d) Interference intensity driven threshold decay calculation: For artificial surface expansion areas, a reduction factor is calculated by the interference decay module based on the enhanced intensity change value : When the final strength change value is greater than the threshold value 0.3, ; When the final strength change value is less than or equal to the threshold value 0.3, ; e) calculate the ecological loss and gain value: calculate the ecological service value change amount according to the following formula: ; wherein, i is a land class index, i = 1: mangrove, 2: salt marsh, 3: mudflat, 4: shallow sea, 5: artificial surface; is the area of the first time point, i land class; is the area of the second time point, i land class; is the spatially corrected area of the first time point, i land class; is a decay factor, and the change in ecological service value caused by artificial surface expansion is negative.

9. The method according to claim 8, wherein, The interference attenuation module construction includes the following steps: a) identify the type of development activities: based on the patch characteristics of artificial surface expansion area, the patch characteristics include patch edge density and patch number, to determine the type of development activities; if the patch edge density > 100 meters / hectare and the patch number > 10, it is determined as port construction; if the patch number = 1 and the area > 1 hectare, it is determined as land reclamation; the rest is determined as tourism facilities; b) obtaining interference duration data, obtaining data of the time span from the start of the development activity to the point in time of the ecological benefit-loss assessment, denoted as duration years T ; c) determining a type-based decay coefficient: according to the determined type of development activity, selecting a corresponding type-based decay coefficient from a pre-defined table of type-based decay coefficients ; Reclamation type: ; Port construction type: ; Tourist facility type: ; d) Calculate time correction factor: based on duration years T , calculate time correction factor : ； e) determine the sensitive area correction factor: obtain the ecological sensitive area spatial data, including mangrove distribution boundary, seagrass bed core area and national coastal protection area range, and judge whether the artificial surface expansion area overlaps with any ecological sensitive area: If there is spatial overlap, set a sensitivity zone correction factor ; If there is no spatial overlap, set ; f) calculating a final attenuation factor: based on , , calculating a final output of the interference attenuation module : calculating an attenuation factor for the applicability of the ecological profit and loss value 。 10. A device for quantifying human activity interference in coastal zones by remote sensing and assessing ecological gains and losses, characterized in that it comprises: including: The data acquisition module is configured to obtain multispectral remote sensing images of the target coastal zone at a first time point and a second time point, the multispectral remote sensing images have a spatial resolution of 10-30 meters and contain blue, green, red, near-infrared and short-wave infrared bands, and simultaneously obtain tidal model data, digital elevation model data and geological background data; The preprocessing module is connected to the data acquisition module and is configured to perform radiation calibration and atmospheric correction on the multispectral remote sensing images; The radiation calibration unit compensates for the non-linear response of the sensor based on the sensor non-linear response compensation mechanism, uses gain compensation coefficients to improve the signal-to-noise ratio in low radiation areas, uses gain compensation coefficients to suppress saturation effects in high radiation areas, and performs polarization filtering algorithm-driven specular reflection suppression on the near-infrared band; The atmospheric correction unit integrates an aerosol scattering model and a humidity correction module, and the humidity correction module generates an aerosol particle swelling effect correction factor based on Mie scattering theory when the relative humidity is greater than 60%. The interference index calculation module is connected with the preprocessing module and is configured to: Land use classification is performed on two-phase surface reflectivity data by using an object-oriented segmentation algorithm, the segmentation scale is 10-30 pixels, and the classification features include NDBI, NDVI, MNDWI, and the mean value of band reflectivity; Temporary inundated areas are dynamically marked based on a tidal model and a digital elevation model, and spectral criteria are used to identify building shadow areas, and the reflectivity of the shadow areas is replaced by the mean value of the same type of surface features in adjacent non-shadow areas; The output land use type change detection result of the non-temporary inundated area and the NDBI after shadow compensation are outputted. The interference quantification module is connected with the interference index calculation module and is configured to: The area of artificial surface expansion is extracted, and a fragmentation interference index is generated by fusing patch density index and edge density index; The NDBI difference is corrected by an interference enhancement model: when the fragmentation interference index is greater than 0.5, the final intensity change value is weighted and outputted by a regional correction coefficient matrix, and the matrix is generated by coupling geological background data and ecological sensitive area data; The ecological evaluation module is connected with the interference quantification module and is configured to: An ecosystem service value model is constructed by using the unit area value equivalent factor method, and the basic value coefficient is spatially corrected by integrating dynamic parameters of the coastal zone; An attenuation factor is calculated based on an interference attenuation module, and the module generates a multi-dimensional correction factor by identifying the development type based on patch characteristics, combining the duration and the overlap state of the ecological sensitive area; According to the formula ; output ecological profit and loss value; The result output module is connected with the ecological evaluation module and is configured to generate a spatial distribution map and a numerical ecological loss and gain evaluation report.

Images