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

Abstract

The application discloses a kind of coast zone human activity interference remote sensing quantification and ecological loss and gain evaluation method and device, belong to remote sensing monitoring technical field, solve the problem of low efficiency of existing technology depending on field investigation and ecological evaluation and remote sensing data disjointedness and reliability insufficient.The gist is as follows: obtain interval 5-10 years, spatial resolution 10-30 meters, two period multispectral remote sensing image containing blue light and other multiband, synchronously integrate tidal, elevation and geological background data;Image radiation calibration and atmospheric correction eliminate environmental error;Based on shortwave and near infrared band reflectivity algorithm normalized building index, combine land use classification change detection quantification interference area;Extract change area and corrected intensity change value input ecosystem service value model, output evaluation result using unit area value equivalent factor method.The method realizes full-process remote sensing automation processing, avoids field investigation, effectively reduces the cost of adjustment, significantly improves the pertinence of management decision.

Description

Technical Field

[0001] This invention belongs to the field of remote sensing monitoring technology, specifically relating to a method and device for quantifying human activity disturbance and assessing ecological damage and benefit in coastal zones based on multispectral remote sensing physical property analysis. Background Technology

[0002] Current methods for quantifying human-induced disturbances in coastal zones primarily rely on field sampling and manual surveys. Coastal zones span vast areas with complex topography, requiring significant manpower and time for on-site investigations. Poor accessibility in areas such as the intertidal zone and nearshore waters leads to spatial blind spots in data collection. In multi-phase monitoring missions, achieving temporal synchronization in fieldwork is challenging, and data from different years suffers from comparability errors due to sampling point offsets or standard differences. Furthermore, field surveys only acquire discrete point information, failing to fully characterize the spatial heterogeneity of human-induced disturbances.

[0003] More importantly, the unique tidal dynamics of the coastal zone pose a severe challenge to remote sensing monitoring. Conventional remote sensing land use classification methods are prone to misclassifying periodically inundated areas (such as bare mudflats at low tide) as permanent land use changes (such as bare land or artificial surfaces), while classifying the same area as water bodies at high tide. This results in numerous false change patches when comparing multiple images, severely interfering with the accurate extraction of areas truly affected by human activities. Simultaneously, issues such as missing image data, abnormal reflections (such as specular reflections in tidal channels), and noise caused by tidal phenomena also disrupt the continuity and accuracy of time series data for key parameters such as vegetation indices, further reducing the reliability of change detection.

[0004] Ecological loss and gain assessment typically employs the static value coefficient method, but current technologies have failed to establish a quantitative conversion mechanism from changes in physical quantities obtained from remote sensing monitoring (such as changes in the Normalized Difference Building Index (NDBI) and land use change area) to losses in ecosystem service value. While remote sensing data can provide information on large-scale surface changes, these physical indicators (such as NDBI, which only reflects changes in building density) cannot be automatically mapped to ecological value losses (such as biodiversity loss and reduced carbon sequestration capacity). Traditional methods require manually inputting disturbance range data obtained from field statistics or simple interpretation into the assessment model, which is prone to introducing human error. This disconnect between remote sensing physical quantities and ecological value assessment results in assessments lacking responsiveness to the ecological impacts of differences in coastal development intensity (such as high-intensity reclamation versus low-intensity infrastructure) and spatial patterns (such as concentrated development versus fragmented development), making it difficult to reflect the nonlinear cumulative destructive effects of high-intensity or sensitive area development on ecological functions.

[0005] The challenges in achieving automated assessment lie in two aspects: First, the aforementioned tidal dynamics render conventional remote sensing classification and change detection results inaccurate, directly impacting the reliability of disturbance area extraction. Second, the ecological effects of human activity disturbance exhibit non-linear characteristics; simple area statistics cannot characterize the cumulative impact of fragmentation and other development patterns, necessitating the establishment of a quantitative correlation model between disturbance intensity and ecological value decay. Current technologies have not yet resolved the cross-scale conversion problem from remote sensing physical quantities to ecological parameters. Summary of the Invention

[0006] To achieve these objectives and other advantages of the present invention, a method for remote sensing quantification of human activity disturbance in coastal zones and assessment of ecological losses and benefits is provided, comprising the following steps:

[0007] S1: Acquire multispectral remote sensing images of the target coastal area at a first time point and a second time point, wherein the first time point is earlier than the second time point, and the time interval between the first time point and the second time point is between 5 and 10 years. The spatial resolution of the multispectral remote sensing images is between 10 meters and 30 meters. The multispectral remote sensing images include blue light band, green light band, red light band, near-infrared band and short-wave infrared band. At the same time, acquire tidal model data, digital elevation model data and geological background data.

[0008] S2: Preprocess the multispectral remote sensing image, including radiometric calibration and atmospheric correction, to eliminate sensor errors and atmospheric scattering effects;

[0009] S3: Based on the preprocessed multispectral remote sensing image, calculate the human activity disturbance index, which includes the normalized building index and land use type change detection. The normalized building index is calculated using the ratio of reflectance of the shortwave infrared band to the near-infrared band. The land use type change detection is achieved by comparing the classification results of the first time point and the second time point. The land use type change detection needs to dynamically mark the temporary flooded area based on the tidal model data and digital elevation model data, and exclude the area from the statistics of change detection.

[0010] S4: Quantify the degree of human activity interference, including calculating the area and intensity change value of the changed area. The area of ​​the changed area is based on the output of land use type change detection, and the intensity change value is calculated based on the difference of the normalized building index. The intensity change value is corrected by the interference enhancement model and the final intensity change value is output. The interference enhancement model generates a regional correction coefficient matrix based on the fragmentation interference index and geological and ecological data, and spatially weights the intensity change value.

[0011] S5: The changes in area and intensity of the changed region are input into the ecosystem service value model, which uses the unit area value equivalent factor method;

[0012] S6: The output of the ecosystem service value model includes spatial distribution maps and numerical reports of ecological cost-benefit assessment results.

[0013] To address the current problem that quantifying human activity disturbances in coastal zones relies excessively on field surveys, leading to inefficiency, and that ecological loss assessments are not directly integrated with remote sensing data, thus affecting assessment reliability, this invention acquires two periods of multispectral remote sensing images (5-10 years apart), eliminates environmental errors through radiometric calibration and atmospheric correction; calculates the Normalized Difference Building Index (NDBI) based on shortwave and near-infrared bands, and quantifies the disturbance area by combining land use classification change detection; inputs the change area and intensity values ​​into an ecosystem service value model (using the unit area value equivalent factor method), and finally outputs a spatially clear ecological loss assessment result. This invention achieves fully automated remote sensing processing, avoiding the costs of field surveys; and establishes a quantitative conversion mechanism from physical quantity changes (NDBI difference) to ecological value loss, improving assessment reliability.

[0014] Preferably, the preprocessing specifically involves: radiometric calibration of the raw digital quantization values ​​of the multispectral remote sensing image to output radiance data at the sensor entrance pupil; atmospheric correction of the radiance data to output surface reflectance data.

[0015] Preferably, radiation calibration includes the following steps:

[0016] a) Read the radiometric calibration coefficient file provided by the sensor manufacturer, and convert the original digital quantization values ​​into radiance data based on the gain coefficient and offset coefficient;

[0017] b) Compensate the sensor nonlinear response for the converted radiance data: when the radiance value is less than 20% of the sensor saturation radiance value, use the low-radiance region gain compensation coefficient to improve the signal-to-noise ratio; when the radiance value is greater than 50% of the sensor saturation radiance value, use the high-radiance region gain compensation coefficient to suppress the saturation effect.

[0018] c) Add specular reflection suppression processing to the near-infrared band radiance data separately, and use polarization filtering algorithm combined with solar altitude angle, sensor observation angle, red band reflectivity characteristics and sea surface wind speed data to correct flare interference.

[0019] The excessively large dynamic range of radiation in coastal scenes leads to nonlinear response errors in sensors. This invention incorporates zonal compensation in its radiometric calibration (improving the signal-to-noise ratio in low-radiation areas and suppressing saturation in high-radiation areas); a separate polarization filtering algorithm is used for the near-infrared band, combined with multi-source data such as solar altitude angle and sea surface wind speed to suppress specular reflection. This improves the accuracy of radiance data and reduces distortion of ground object reflectivity.

[0020] Preferably, atmospheric correction includes the following steps:

[0021] a) Acquire aerosol optical thickness data and relative humidity data that match the imaging time of the multispectral remote sensing image. The aerosol optical thickness data is derived from satellite remote sensing inversion products, with a spatial resolution consistent with the multispectral remote sensing image. The relative humidity data is derived from meteorological reanalysis data, with a temporal resolution less than or equal to 6 hours.

[0022] b) Based on aerosol optical thickness data and relative humidity data, calculate the aerosol scattering characteristic parameters at a wavelength of 550 nm. The calculation uses an aerosol scattering model that includes the mixing ratio parameters of sea salt, dust and anthropogenic pollution aerosols. The mixing ratio parameters are set according to historical observation data of the target coastal area.

[0023] c) Convert the radiometrically calibrated radiance data into apparent reflectance data. The conversion formula is:

[0024] ;in, This is the radiance value. d This is the Earth-Sun distance correction factor. Solar irradiance outside the atmosphere. The solar zenith angle;

[0025] d) Atmospheric correction is performed using a radiative transfer model, the input parameters of which include aerosol scattering characteristics, sensor imaging geometry, surface elevation data, and apparent reflectivity.

[0026] e) In the radiative transfer model, when the relative humidity data exceeds 60%, a humidity correction module is introduced. This module uses Mie scattering theory to calculate the correction factor for the water vapor's effect on the aerosol particle-scale expansion effect, so as to improve the model's accurate simulation of the influence of humidity.

[0027] The high humidity and complex aerosol composition of coastal areas lead to insufficient accuracy in traditional atmospheric correction models. This invention employs an aerosol type parameterized model (mixed sea salt / dust / anthropogenic pollution type) and introduces a relative humidity-driven Mie scattering correction factor to quantify particle expansion effects. This accurately corrects for the influence of aerosol scattering and improves the reliability of surface reflectance data.

[0028] Preferably, calculating the human activity disturbance index includes the following steps:

[0029] a) Land use classification was performed on the surface reflectance data at the first and second time points. An object-oriented segmentation algorithm was used for classification, and the segmentation scale parameter was set between 10 and 30 pixels. The selected features included the normalized building index, the normalized vegetation index, the improved normalized water index, and the mean of the surface reflectance band.

[0030] Normalized Building Index ;

[0031] Normalized Difference Vegetation Index ;

[0032] Improved Normalized Water Index ;

[0033] in, The surface reflectance in the blue light band. The surface reflectance is in the green light band. For shortwave infrared band surface reflectance, Near-infrared surface reflectance, The surface reflectance is in the red band.

[0034] b) Based on the regions where the improved normalized water index value is greater than the water body determination threshold in the classification results, and combined with the tidal height data estimated by the tidal model, the following steps are performed:

[0035] When the tide height is 1 meter higher than the local mean sea level, the area with an elevation lower than the tide height and classified as a body of water is marked as a temporary flooding zone.

[0036] Elevation data are derived from digital elevation models; temporary inundation areas are not included in land use type change monitoring.

[0037] c) In non-submerged and non-water body areas, determine the building shadow area based on the spectral characteristics of the vegetation cover area. If the vegetation area is within... and At that time, it was determined to be a building shadow area;

[0038] The reflectance of the shaded area is replaced by the average reflectance value of the same type of ground features in adjacent non-shaded areas.

[0039] d) Recalculate NDBI based on the reflectance data after shading compensation;

[0040] e) Land use type change detection is achieved by comparing the classification results at the first and second time points, and only the pixels in non-temporary flooded areas that have undergone type transformation are counted.

[0041] Intertidal water level fluctuations lead to misjudgments of land use and reduced index accuracy due to building shadows. This invention uses a tidal model and digital elevation model to dynamically mark temporarily flooded areas. In vegetated areas, building shadows are identified using spectral criteria (near-infrared / red band ratio <1.2 and red-blue band reflectance characteristics <0.2), and the reflectance is replaced by the average value of similar land features in adjacent non-shaded areas. By reducing the impact of tides and shadows on classification, the accuracy of extracting interference areas is improved.

[0042] Preferably, quantifying the degree of human activity interference includes the following steps:

[0043] a) Based on the monitoring results of land use type change, extract the artificial surface expansion area and calculate the total area of ​​the changed area;

[0044] b) Conduct 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 by dividing the number of patches by the total area, and calculate the edge density index by dividing the total length of patch edges by the total area.

[0045] c) The patch density index and the edge density index are fused with weights of 30% to 40% and 60% to 70%, respectively, to obtain the fragmentation interference index;

[0046] d) Calculate the difference between the normalized building index at the first time point and the second time point as the base strength change value;

[0047] e) Input the base intensity change value into the disturbance enhancement model, and the disturbance enhancement model will execute:

[0048] When the fragmentation interference index is greater than the threshold of 0.5, the output enhancement intensity change value is multiplied by the basic intensity change value and the correction coefficient.

[0049] If the fragmentation interference index does not exceed the threshold of 0.5, the output enhancement intensity change value will be consistent with the basic intensity change value.

[0050] f) Output the enhanced strength change value as the final strength change value.

[0051] The representation of fragmentation development interference in the base intensity value is insufficient. This invention generates a fragmentation interference index based on the fusion of patch density and edge density; when the index > 0.5, a correction coefficient is used to enhance the base intensity variation value. This strengthens the quantitative representation of interference intensity from distributed development activities.

[0052] Preferably, the construction of the interference enhancement model includes the following steps:

[0053] a) Obtain geological background data of the target coastal zone, including the proportion of bedrock coastline length and the spatial distribution of tidal flat sediment types, and extract spatial data of ecologically sensitive areas, including mangrove distribution boundaries, seagrass bed core areas and the scope of national coastal protected areas;

[0054] b) Generate a regional correction coefficient matrix based on geological background data and ecologically sensitive area data. The value ranges from 0.8 to 1.3. This matrix is ​​generated by coupling the fragmentation disturbance index with geological and ecological data.

[0055] c) The change in reinforcement strength Weighted by the regional interference weight matrix:

[0056] ;

[0057] d) Output As a value representing the change in strength.

[0058] To address the issue that a single correction coefficient cannot adapt to the spatial heterogeneity of coastal geological and ecological environments, this invention generates a spatial weight matrix by coupling geological background data (bedrock shoreline ratio, sediment type) with ecologically sensitive area data (mangroves, seagrass beds), and then applies regional weights to the enhancement intensity variation values. This achieves spatial adaptive correction of the quantification results of interference intensity.

[0059] Preferably, the construction of the ecosystem service value model includes the following steps:

[0060] a) Establish a basic value equivalent table, including the annual service value coefficient per unit area for five land types: mangroves, salt marshes, tidal flats, shallow seas, and artificial surfaces. ,in i =1: Mangroves, 2: Salt marshes, 3: Mudflats, 4: Shallow seas, 5: Artificial surfaces; The coefficients are derived from the "Ecosystem Service Value Accounting Standard" issued by the national ecological protection authorities.

[0061] b) Collect dynamic parameters of the coastal zone geographic environment: Obtain the following spatial parameters for the target coastal zone region:

[0062] Intertidal zone width data;

[0063] Aquatic biodiversity index, which is the abundance of fish species calculated based on fishery resource survey data;

[0064] Historical frequency data of storm surges;

[0065] c) Generate the coastal zone spatial correction coefficient matrix: Based on the parameters collected in the steps, adjust the value coefficients... Dynamic adjustments are made to calculate the spatially adjusted annual service value coefficient per unit area. ;

[0066] When the intertidal zone width increases by more than 500 meters from the baseline value, the tidal flat type... ;

[0067] For every 10 fish species added to the abundance, the abundance of shallow marine species... ;

[0068] When storm surges occur more frequently than 0.5 times per year, mangrove and salt marshland types...

[0069] ;

[0070] ;

[0071] For land categories or parameters that do not meet the above correction conditions ;

[0072] d) Threshold attenuation calculation driven by interference intensity: For artificial land expansion areas, the reduction factor is calculated through the interference attenuation module based on the change in enhancement intensity. :

[0073] When the final intensity change value is greater than the threshold of 0.3, ;

[0074] When the final intensity change is less than or equal to the threshold of 0.3,

[0075] e) Calculate the ecological gain / loss value: Calculate the change in the value of ecosystem services using the following formula:

[0076] ;in, i For land category index, i =1: Mangrove forest, 2: Salt marsh, 3: Mudflats, 4: Shallow sea, 5: Artificial surface; For the first time point i Land type area; For the second time point i Land type area; the spatially corrected area of ​​the first type. i Annual service value coefficient per unit area for each land category; As a decay factor, the change in ecosystem service value caused by artificial land expansion is negative.

[0077] To address the issue that static value coefficients cannot reflect the spatial dynamics of coastal zones and the nonlinear ecological impacts of high-intensity development, this invention dynamically corrects the basic value coefficients based on intertidal width, biodiversity index, and storm surge frequency. For artificial land expansion zones, a disturbance attenuation factor (0.7-0.9 when the value is 0.3) is applied based on the final intensity change value. This improves the spatiotemporal adaptability of calculating changes in ecosystem service value.

[0078] Preferably, the interference attenuation module is constructed by including the following steps:

[0079] a) Identify the type of development activity: Based on the patch characteristics of artificial land expansion areas, including patch edge density and number of patches, determine the type of development activity;

[0080] If the density of patch edges is greater than 100 meters per hectare and the number of patches is greater than 10, it is determined to be port construction.

[0081] If the number of patches is 1 and the area is greater than 1 hectare, it is determined to be land reclamation.

[0082] The rest were classified as tourist facilities;

[0083] b) Obtain disturbance duration data, specifically the time span from the start of development activities to the ecological cost-benefit assessment point, denoted as duration in years. T ;

[0084] c) Determine the type-based attenuation coefficient: Based on the determined development activity type, select the corresponding basic attenuation coefficient from the preset type-based attenuation coefficient table. ;

[0085] Types of land reclamation: ;

[0086] Port construction types: ;

[0087] Types of tourist facilities: ;

[0088] d) Calculate the time correction factor: based on the duration in years. T Calculate the time correction factor :

[0089] ;

[0090] e) Determine the sensitive area correction factor: Obtain spatial data of ecologically sensitive areas, including mangrove distribution boundaries, seagrass bed core areas, and the scope of national coastal protected areas, and determine whether artificial land expansion areas spatially overlap with any ecologically sensitive area:

[0091] If spatial overlap exists, a sensitive area correction factor is set. ;

[0092] If there is no spatial overlap, then set ;

[0093] f) Calculate the final decay factor: based on , , The final output of the interference attenuation module is an attenuation factor used to calculate the applicability of ecological gain / loss values. :

[0094] .

[0095] This invention addresses the issue of interference attenuation factors failing to consider the cumulative effects of development type, duration, and ecologically sensitive areas. It identifies development type (reclamation / port / tourism) based on patch features; and combines this with duration (…). TCalculate the time correction factor; set the sensitivity zone correction factor based on whether it overlaps with ecologically sensitive areas; multiply the development type, duration, and sensitivity factor, and constrain the result to the interval [0.7, 0.9], using this as the output value of the final attenuation factor. This more accurately quantifies the cumulative ecological effects of high-intensity development.

[0096] Other advantages, objectives and features of the present invention will become apparent in part from the following description, and in part from those skilled in the art through study and practice of the invention. Detailed Implementation

[0097] The present invention will now be described in further detail so that those skilled in the art can implement it based on the description.

[0098] It should be understood that terms such as “having,” “comprising,” and “including” as used herein do not exclude the presence or addition of one or more other elements or combinations thereof.

[0099] Example 1

[0100] A remote sensing quantification method for human activity disturbance in coastal zones and an assessment of ecological losses and benefits, comprising the following steps:

[0101] S1: Acquire multispectral remote sensing images of the target coastal area at a first time point and a second time point, wherein the first time point is earlier than the second time point, with a time interval of 5 to 10 years, and the spatial resolution of the multispectral remote sensing images is between 10 meters and 30 meters. The multispectral remote sensing images include blue light band, green light band, red light band, near-infrared band, and shortwave infrared band, and simultaneously acquire tidal model data, digital elevation model data, and geological background data;

[0102] S2: Preprocess the multispectral remote sensing image, including radiometric calibration and atmospheric correction, to eliminate sensor errors and atmospheric scattering effects;

[0103] S3: Based on the preprocessed multispectral remote sensing image, calculate the human activity disturbance index, which includes the normalized building index and land use type change detection. The normalized building index is calculated using the ratio of reflectance of the shortwave infrared band to the near-infrared band. The land use type change detection is achieved by comparing the classification results of the first time point and the second time point. The land use type change detection needs to dynamically mark the temporary flooded area based on the tidal model data and digital elevation model data, and exclude the area from the statistics of change detection.

[0104] S4: Quantify the degree of human activity interference, including calculating the area and intensity change value of the changed area. The area of ​​the changed area is based on the output of land use type change detection, and the intensity change value is calculated based on the difference of normalized building index. Through fine correction of the interference enhancement model, the intensity change value is finally obtained.

[0105] S5: The changes in area and intensity of the changed region are input into the ecosystem service value model, which is based on the unit area value equivalent factor method.

[0106] S6: The output of the ecosystem service value model includes spatial distribution maps and numerical reports of ecological cost-benefit assessment results.

[0107] This implementation method uses a bay area as an example. First, Landsat multispectral remote sensing images of the area from 2000 and 2010 were acquired, with a spatial resolution of 30 meters, including blue, green, red, near-infrared, and shortwave infrared bands. Simultaneously, contemporaneous tidal model data, a 30-meter precision digital elevation model, and geological background vector data were collected.

[0108] Traditional methods require organizing field teams to set up sampling points along the coast, manually recording the extent of building expansion and drawing sketches, which takes about three months. Furthermore, limited by handheld GPS positioning errors (±5 meters) and map registration deviations, the actual spatial positioning accuracy of the field survey results is 52.3 meters (calculated at 100 verification points, with a 95% confidence interval error of 45-60 meters), and can only present statistical averages at the county scale. In contrast, this invention directly uses satellite imagery, with data acquisition taking no more than one week, significantly improving the efficiency and real-time nature of data collection; through fully automated remote sensing processing, it directly generates a spatial distribution map of ecological damage and loss with a 30m × 30m grid precision, significantly improving the targeting of management decisions.

[0109] Preprocessing was performed on the two imagery phases: radiometric calibration was completed based on sensor parameter files, converting the raw digital values ​​into radiance; surface reflectance data was output using a coastal zone-specific atmospheric correction model. In traditional methods, field surveyors need to carry spectrometers to measure reflectance in the field, and can only complete sampling at 2-3 locations per day; this invention generates full-area surface reflectance data in a single step through remote sensing physical quantity processing.

[0110] The human activity interference index is calculated as follows: Object-oriented classification is performed on two periods of reflectance data, with a segmentation scale of 20 pixels. Five land cover categories are defined based on NDBI, NDVI, MNDWI, and band mean values. The classification results are compared to extract artificial land expansion areas, and the changed area is statistically analyzed. Simultaneously, the difference in NDBI between the two periods is calculated as the baseline intensity change value. Traditional methods rely on manual visual interpretation of changed patches, a process that is highly subjective and difficult to accurately quantify intensity changes. This invention achieves automated quantification through classification comparison and NDBI difference. The classification process can integrate a computational model based on convolutional neural networks to replace or assist traditional object-oriented segmentation algorithms, improving the accuracy of land cover boundary identification and type differentiation. A random forest regression model can be introduced into the intensity change value calculation to correct for nonlinear interference effects. The convolutional neural network model is applied to the land use classification task, automatically extracting multispectral features through the trained network, replacing manually designed features, significantly improving classification accuracy. The random forest model is used to fuse multi-source indicators (such as patch density and edge length) to predict interference intensity, reducing the subjectivity of manually set weights.

[0111] The changed area and intensity values ​​are input into the ecosystem service value model: using the unit area value coefficient in the "Ecosystem Service Value Accounting Standard," combined with the dynamic correction coefficient values ​​of regional intertidal width and fish species abundance, an ecological loss and gain spatial distribution map and a value change report are output. Traditional methods require manually inputting the expanded range drawn on-site into the assessment software, which is prone to introducing boundary errors; this invention directly generates quantitative assessment results through the automatic conversion of remote sensing physical quantities to ecological parameters.

[0112] The final results show that artificial land expansion from 2000 to 2010 led to a loss of ecosystem service value, and the spatial distribution map clearly shows the high-loss areas where core development areas and ecologically sensitive areas overlap. Due to the limitations of traditional methods in terms of field sampling density, their assessment reports often only present average values ​​at the county level.

[0113] Furthermore, the preprocessing specifically involves: radiometric calibration of the raw digital quantization values ​​of the multispectral remote sensing image to output radiance data at the sensor entrance pupil; atmospheric correction of the radiance data to output surface reflectance data.

[0114] After acquiring Landsat imagery from 2000 and 2010, the raw digital quantization values ​​are first input into the radiometric calibration module. By calling the coefficient files provided by the sensor manufacturer and using gain and offset calculation methods, the raw pixel values ​​are accurately converted into radiometric data. In traditional methods, different teams may use custom calibration formulas, resulting in inconsistent radiometric units in imagery from the same period. For example, one study used linear transformation while ignoring the nonlinear characteristics of the sensor, causing a systematic deviation in the radiometric magnitudes of dark pixels (such as deep water areas) and bright pixels (such as beaches) in the coastal zone.

[0115] The radiance data is then input into the atmospheric correction module. Combined with aerosol optical thickness products from satellite transit and relative humidity data from meteorological stations, surface reflectance results are generated. In existing technologies, some processes skip the atmospheric correction step and directly use apparent reflectance, or employ simplified models that ignore the specific impact of humidity on aerosols. For example, a coastal zone study directly used apparent reflectance to calculate vegetation indices, leading to an overestimation of nearshore water reflectance by 10%–15%, mistakenly identifying intertidal silt areas as vegetation-degraded zones.

[0116] This implementation strictly limits the preprocessing output to two types of standardized data: sensor entrance pupil radiance (unit: W·m²). -2 ·sr -1 ·μm -1 The present invention uses apparent reflectance (dimensionless) and surface reflectance. Traditional methods, lacking standardized intermediate data, struggle to integrate results from different projects. For example, in a ten-year assessment report for a certain region, apparent reflectance was used for data before 2005, while surface reflectance was used after 2005, resulting in artificially created breaks in the NDBI time series. This invention, however, unifies the data benchmark, ensuring direct comparison of physical quantities between two image periods.

[0117] The final output surface reflectance data is used for subsequent land use classification. Compared with traditional methods, standardized preprocessing eliminates classification bias caused by differences in data sources. In a prior art case, using different atmospheric correction models to process the same image resulted in a difference of 8.2% in the classified area of ​​salt marshes; however, this invention, by forcing a standard output, ensures that the classification results only reflect actual surface changes.

[0118] Furthermore, radiation calibration includes the following steps:

[0119] a) Read the radiometric calibration coefficient file provided by the sensor manufacturer, and convert the original digital quantization values ​​into radiance data based on the gain coefficient and offset coefficient;

[0120] b) Compensate the sensor nonlinear response for the converted radiance data: when the radiance value is less than 20% of the sensor saturation radiance value, use the low-radiance region gain compensation coefficient to improve the signal-to-noise ratio; when the radiance value is greater than 50% of the sensor saturation radiance value, use the high-radiance region gain compensation coefficient to suppress the saturation effect.

[0121] c) Add specular reflection suppression processing to the near-infrared band radiance data separately, and use polarization filtering algorithm combined with solar altitude angle, sensor observation angle, red band reflectivity characteristics and sea surface wind speed data to correct flare interference.

[0122] Traditional methods directly use linear calibration coefficients provided by sensor manufacturers, ignoring the large dynamic range of coastal radiation. For example, a certain existing technology uses the same transformation parameters for deep water (low radiation) and salt pan (high radiation) areas in the same image, resulting in insufficient signal-to-noise ratio in deep water and saturation effect in salt pan areas. Specifically, the difference in radiance between mangrove shadows and turbid water is less than 5%, but the confusion rate during classification reaches 23%; the edges of salt pans lose texture details due to signal truncation and are misclassified as built-up areas.

[0123] This plan implements zoned compensation:

[0124] Low-radiation area processing: Identify pixels with radiation values ​​below 20% of the saturation value (such as deep water areas and mangrove canopy shadows), and use a 1.15x gain compensation coefficient to improve the signal-to-noise ratio, thus expanding the difference in reflectivity between shaded water bodies and vegetation to a distinguishable range.

[0125] High radiation area processing: For pixels exceeding 50% of the saturation value (such as salt fields and bare beaches), apply a gain suppression coefficient of 0.9 to prevent signal overflow and preserve the crystal texture characteristics of salt fields.

[0126] Near-infrared flare suppression: Targeting the common southeasterly wind of 5.2 m / s in the bay, and combining a solar altitude angle of 58 degrees with a sensor observation angle of 8 degrees, a polarization filtering algorithm is employed to correct specular reflection. Compared to existing technologies that only use red-band thresholding, this method is insufficient in eliminating the impact of wind speed, and the residual flare problem often leads to nearshore waters being incorrectly classified as bare land.

[0127] After the improvement, the radiation response curves of deep water areas and salt fields are closer to the measured spectra, and the misclassification rate of mangrove boundaries is significantly reduced. Traditional methods often require manual intervention to correct classification results due to radiation distortion, a process that takes up to two weeks; this invention, through adaptive compensation, directly generates reliable data for subsequent processes.

[0128] Furthermore, atmospheric correction includes the following steps:

[0129] a) Acquire aerosol optical thickness data and relative humidity data that match the imaging time of the multispectral remote sensing image. The aerosol optical thickness data is derived from high-precision satellite remote sensing inversion products, and its spatial resolution is perfectly matched with the multispectral remote sensing image. The relative humidity data is derived from the latest meteorological reanalysis data, with a temporal resolution accurate to within 6 hours.

[0130] b) Based on aerosol optical thickness data and relative humidity data, calculate the aerosol scattering characteristic parameters at a wavelength of 550 nm. The calculation uses an aerosol scattering model that includes the mixing ratio parameters of sea salt type, dust type and anthropogenic pollution type aerosols. The mixing ratio parameters are set according to historical observation data of the target coastal area.

[0131] c) Convert the radiometrically calibrated radiance data into apparent reflectance data. The conversion formula is:

[0132] ;in, This is the radiance value. d This is the Earth-Sun distance correction factor. Solar irradiance outside the atmosphere. The solar zenith angle;

[0133] d) Atmospheric correction is performed using a radiative transfer model, the input parameters of which include aerosol scattering characteristics, sensor imaging geometry, surface elevation data, and apparent reflectivity.

[0134] e) In the radiative transfer model, when the relative humidity data exceeds 60%, a humidity correction module is introduced. This module uses Mie scattering theory to calculate the correction factor for the water vapor's effect on the aerosol particle size expansion effect, in order to ensure the accuracy of the model.

[0135] Existing technologies typically employ fixed-type aerosol models; for example, one study uniformly used marine aerosol parameters for coastal imagery. This simplistic approach fails to consider the actual impact of industrial emissions on bays, leading to a systematic underestimation of surface reflectance in near-shore urban areas. Specifically, anthropogenic pollutants account for over 40% of aerosols above industrial areas, yet traditional models still treat them as pure sea salt models, resulting in an 8%–12% error in building roof reflectance, which in turn affects the accuracy of identifying artificial foundation spread.

[0136] This solution implements refined correction: First, it acquires aerosol optical thickness products for the day the satellite passes overhead, matching the spatial resolution to 30-meter image data. Simultaneously, it integrates relative humidity fields at 6-hour intervals from meteorological reanalysis data to identify under-cloud areas with humidity exceeding 75%. Based on historical observation data, it sets aerosol mixing ratios—20% sea salt type, 30% dust type, and 50% anthropogenic pollution type for the northern industrial area of ​​the Gulf; and increases the sea salt type ratio to 60% for the southern tourist area. When humidity exceeds the 60% threshold, the Mie scattering correction module is activated to quantify the water vapor expansion effect. This approach differs from a certain existing technology that ignores the impact of humidity changes and directly applies dry season correction results to rainy season images, leading to an abnormal increase in reflectivity in mangrove areas, which is then incorrectly identified as vegetation degradation.

[0137] The corrected data significantly improves classification reliability. Traditional methods generate 15% false building expansion patches in port areas due to inaccurate aerosol modeling; this invention, through dynamic parameter adjustment, makes the spatial gradient of reflectance between industrial areas and ecological protection areas more consistent with the field survey results, providing a reliable input for interference quantification.

[0138] Furthermore, calculating the human activity disturbance index includes the following steps:

[0139] a) Land use classification was performed on the surface reflectance data at the first and second time points. An object-oriented segmentation algorithm was used for classification, and the segmentation scale parameter was set between 10 and 30 pixels. The selected features included the normalized building index, the normalized vegetation index, the improved normalized water index, and the mean of the surface reflectance band.

[0140] Normalized Building Index ;

[0141] Normalized Difference Vegetation Index ;

[0142] Improved Normalized Water Index ;

[0143] in, The surface reflectance in the blue light band. The surface reflectance is in the green light band. For shortwave infrared band surface reflectance, Near-infrared surface reflectance, The surface reflectance is in the red band.

[0144] b) Based on the regions where the improved normalized water index value is greater than the water body determination threshold in the classification results, and combined with the tidal height data estimated by the tidal model, the following steps are performed:

[0145] When the tide height is 1 meter higher than the local mean sea level, the area with an elevation lower than the tide height and classified as a body of water is marked as a temporary flooding zone.

[0146] Elevation data are derived from digital elevation models; temporary inundation areas are not included in land use type change monitoring.

[0147] c) In non-submerged and non-water body areas, determine the building shadow area based on the spectral characteristics of the vegetation cover area. If the vegetation area is within... and At that time, it was determined to be a building shadow area;

[0148] The reflectance of the shaded area is replaced by the average reflectance value of the same type of ground features in adjacent non-shaded areas.

[0149] d) Recalculate NDBI based on the reflectance data after shading compensation;

[0150] e) Land use type change detection is achieved by comparing the classification results at the first and second time points, and only the pixels in non-temporary flooded areas that have undergone type transformation are counted.

[0151] Traditional land use classification methods fail to adequately account for the impact of tidal dynamics, leading to significant errors in classification results. For example, one technique classified exposed mudflats as bare land during low tide and as water bodies during high tide, resulting in a false change patch rate of up to 18.7% over two years. Furthermore, shadows cast by port structures in vegetated areas are misclassified as vegetation degradation zones due to reflectivity distortion, thus artificially inflating the extent of disturbance.

[0152] This solution implements dynamic masking and shadow compensation:

[0153] Tidal inundation zone marking: Combining satellite transit tidal height data with digital elevation models, when the tide level exceeds the mean sea level by 1.2 meters, classified water bodies with elevations below this value are marked as temporary inundation zones. For example, 650 hectares of tidal channels are excluded from unpredictable tidal ranges at high tide, avoiding misclassification between bare land and water bodies.

[0154] Building Shadow Identification and Restoration: In vegetation-covered areas, spectral anomalies are detected—regions with a near-infrared / red band ratio below 1.15 and a red-blue reflectance difference of less than 0.19 are identified as shadows. For the mangrove area behind the port, we filled in the detected shadow patches using the average reflectance value of adjacent non-shaded mangrove areas. In contrast, a certain existing technology using only the NDVI thresholding method failed to distinguish between shadows and true vegetation degradation, generating 42 hectares of false degradation patches around the dock.

[0155] After optimization, the accuracy of intertidal zone change detection is significantly improved. Traditional methods require manual correction when dealing with tidal misjudgment areas, a process that can take up to three weeks. This invention, through automated dynamic masking, directly outputs reliable change area data for interference quantification.

[0156] Furthermore, quantifying the degree of human interference includes the following steps:

[0157] a) Based on the monitoring results of land use type change, extract the artificial surface expansion area and calculate the total area of ​​the changed area;

[0158] b) Perform patch analysis on the artificial land expansion area, count the number of patches and the total length of patch edges. The patch density index is calculated by dividing the number of patches by the total area, while the edge density index is obtained by dividing the total length of patch edges by the total area.

[0159] c) The patch density index and the edge density index are combined with weights of 30% to 40% and 60% to 70% respectively to form a fragmentation disturbance index;

[0160] d) Calculate the difference between the normalized building index at the first time point and the second time point as the base strength change value;

[0161] e) Input the base intensity change value into the disturbance enhancement model, and the disturbance enhancement model will execute:

[0162] When the fragmentation interference index is greater than the threshold of 0.5, the output enhancement intensity change value is multiplied by the basic intensity change value and the correction coefficient.

[0163] When the fragmentation interference index is less than or equal to the threshold of 0.5, the output enhancement intensity change value is equal to the basic intensity change value;

[0164] f) Output the enhanced strength change value as the final strength change value.

[0165] Current technology uses only the difference between two NDBI values ​​to characterize development intensity, leading to an underestimation of fragmented development. In one existing case, a scattered cluster of holiday villas (each villa ranging from 0.1 to 0.3 hectares) and a centralized marina area (5 hectares) received the same intensity value, failing to reflect the differences in ecological impact. Traditional methods generally assign an intensity value of 0.4 to both a 20-hectare villa cluster and a 15-hectare marina area; however, actual monitoring results show that the villa cluster's fragmentation effect on bird habitats is far greater than that of the marina area.

[0166] This plan implements fragmentation index correction:

[0167] Patch feature extraction: Morphological analysis of the artificial land expansion area identified 128 patches (121 of which were villa clusters), with a total edge length of 18.6 kilometers. The calculated patch density index was 6.4 patches / hectare, and the edge density index was 93 meters / hectare.

[0168] Index fusion and correction: The fragmentation disturbance index, fused with weights of 35% and 65%, is 0.72, exceeding the threshold of 0.5. This triggers a correction mechanism, multiplying the base NDBI difference of 0.25 by an enhancement factor of 1.28, resulting in a final intensity value of 0.32. This contrasts with a previous study that used a fixed area weighting method to set the intensity for the dock area at 0.42, while the intensity for the villa complex was only 0.19. This setting contradicts the disturbance levels assessed by ecologists in the field.

[0169] The corrected intensity values ​​better align with ecological effect observations. In the quantification results of traditional methods, the ecological disturbance intensity value of the dock area is 2.2 times that of the villa cluster. However, by enhancing the characterization of fragmented areas using this invention, we found that the ecological disturbance intensity value of the villa cluster actually exceeds that of the dock area by 12%. This finding more accurately reflects the cumulative destructive effect of distributed development on habitat connectivity.

[0170] In addition, the process of constructing the interference enhancement model includes the following steps:

[0171] a) Obtain geological background data of the target coastal zone, including the proportion of bedrock coastline length and the spatial distribution of tidal flat sediment types, and extract spatial data of ecologically sensitive areas, including mangrove distribution boundaries, seagrass bed core areas and the scope of national coastal protected areas;

[0172] b) Generate a regional correction coefficient matrix based on geological background data and ecologically sensitive area data. The value ranges from 0.8 to 1.3. This matrix is ​​generated by coupling the fragmentation disturbance index with geological and ecological data.

[0173] c) The change in reinforcement strength Weighted by the regional interference weight matrix:

[0174]

[0175] d) Output As a value representing the change in strength.

[0176] Current technologies use a single correction coefficient to process data across the entire area, ignoring the spatial heterogeneity of geology and ecology. One existing case applied a uniform enhancement coefficient of 1.25 to both the bedrock shore industrial zone and the tidal flat tourist zone. This led to two major problems: the bedrock zone, due to its geological stability, experienced less actual disturbance, resulting in a 18% overestimation of the corrected intensity value; conversely, the mangrove tidal flat zone, due to its high ecological sensitivity, received insufficient correction, thus severely underestimating the ecological impact of the 30-hectare resort construction.

[0177] This scheme implements spatial weighted correction:

[0178] Regional weight matrix construction: Geological data of 82% of the bedrock coastline on the north bank were extracted and weighted by combining the mangrove distribution boundary. The bedrock area was assigned a weight of 0.85 due to its strong resistance to interference, the core mangrove area on the south bank was assigned a weight of 1.25 due to its high ecological sensitivity, and the ordinary tidal flat area was set as the baseline value of 1.0.

[0179] Dynamic Intensity Modulation: The output villa complex enhancement intensity value of 0.32 is input into this module. After a weight adjustment of 0.85, the intensity value of the resort in the north bedrock area (original value 0.32) decreases to 0.27; while the intensity value of the viewing platform in the south mangrove area (original value 0.28) increases to 0.35 with a weight of 1.25. This contrasts with a certain existing fixed-coefficient method, which sets the intensity value of the viewing platform at 0.35 and the resort at 0.40, contradicting the ecologist's assessment—field monitoring shows that the viewing platform caused a 23% mortality rate for mangrove larvae, far exceeding the 9% mortality rate of the resort.

[0180] The correction results better align with the spatial patterns of the coastal zone. Traditional methods quantify the intensity on the north shore as 1.14 times that on the south shore; this invention, through geological-ecological coupling weighting, makes the intensity value of sensitive areas on the south shore exceed that on the north shore by 29%, accurately reflecting the exponential ecological risk from mangrove development. The final intensity distribution map has been verified by marine management departments and is directly used for ecological restoration priority delineation.

[0181] This matrix was calculated by coupling the fragmentation disturbance index with geological and ecological data:

[0182] The geological stability factor G is calculated as follows: G = bedrock shoreline length ratio × 0.01 + sediment disturbance resistance coefficient (sandy tidal flats = 0.8, muddy tidal flats = 1.2, mixed type = 1.0).

[0183] Ecological sensitivity factor E was extracted: core area of ​​mangrove or seagrass bed = 1.3, buffer zone of protected area = 1.1, other areas = 1.0;

[0184] Based on fragmentation interference index F index Calculate the dynamic modulation term ;

[0185] Output .

[0186] The construction of the ecosystem service value model includes the following steps:

[0187] a) Establish a basic value equivalent table, including the annual service value coefficient per unit area for five land types: mangroves, salt marshes, tidal flats, shallow seas, and artificial surfaces. ,in i =1: Mangrove forest, 2: Salt marsh, 3: Mudflats, 4: Shallow sea, 5: Artificial surface;

[0188] b) Collect dynamic parameters of the coastal zone geographic environment: Obtain the following spatial parameters for the target coastal zone region:

[0189] Intertidal zone width data;

[0190] Aquatic biodiversity index, which is the abundance of fish species calculated based on fishery resource survey data;

[0191] Historical frequency data of storm surges;

[0192] c) Generate the coastal zone spatial correction coefficient matrix: Based on the parameters collected in the steps, adjust the value coefficients... Dynamic adjustments are made to calculate the spatially adjusted annual service value coefficient per unit area. ;

[0193] When the intertidal zone width increases by more than 500 meters from the baseline value, the tidal flat type... ;

[0194] For every 10 fish species added to the abundance, the abundance of shallow marine species... ;

[0195] When storm surges occur more frequently than 0.5 times per year, mangrove and salt marshland types...

[0196] ;

[0197] ;

[0198] For land categories or parameters that do not meet the above correction conditions ;

[0199] d) Threshold attenuation calculation driven by interference intensity: For artificial land expansion areas, the reduction factor is calculated through the interference attenuation module based on the change in enhancement intensity. :

[0200] When the final intensity change value is greater than the threshold of 0.3, ;

[0201] When the final intensity change is less than or equal to the threshold of 0.3,

[0202] e) Calculate the ecological gain / loss value: Calculate the change in the value of ecosystem services using the following formula:

[0203] ;in, i For land category index, i =1: Mangrove forest, 2: Salt marsh, 3: Mudflats, 4: Shallow sea, 5: Artificial surface; For the first time point i Land type area; For the second time point i Land type area; the spatially corrected area of ​​the first type. i Annual service value coefficient per unit area for each land category; As a decay factor, the change in ecosystem service value caused by artificial land expansion is negative.

[0204] Current technologies directly apply static coefficients from the "Ecosystem Service Value Accounting Standard," ignoring the dynamic geographical characteristics of the coastal zone. One existing case used fixed value coefficients for the entire area from 2000 to 2010: a uniform 28,000 yuan / hectare / year for tidal flats and 12,000 yuan / hectare / year for shallow seas. However, actual monitoring showed that the intertidal zone on the north shore widened by 600 meters due to siltation, improving the quality of shellfish habitat; while on the south shore, the increased frequency of storm surges to an average of 0.7 times per year significantly enhanced the protective value of mangroves. Traditional methods failed to incorporate these changes, resulting in the ecological gains on the north shore not being reflected, and the disaster prevention value on the south shore being underestimated by approximately 30%.

[0205] This plan will be dynamically revised:

[0206] Spatial parameter integration: Based on hydrological mapping data, it was confirmed that the intertidal zone width on the north bank increased by 520 meters, triggering the mudflat value coefficient to rise to 33,600 yuan / hectare / year; the storm surge frequency on the south bank reached 0.7 times / year, and the mangrove and salt marsh coefficients simultaneously increased to 1.2 times the benchmark value.

[0207] Disturbance attenuation mechanism: For the 15-hectare reclaimed area on the south bank, the final intensity change value of 0.41 exceeded the threshold of 0.3, so an attenuation factor of 0.75 was applied. The intensity value of the 8-hectare tourist facilities on the north bank was 0.29, which did not reach the threshold, so an attenuation factor of 1.0 was applied. In contrast, a certain existing technology uniformly applies an attenuation factor of 0.8 to all artificial surfaces, which underestimates the ecological loss of the reclaimed area to 1.6 times that of the tourist area, while the ecological model shows that the actual loss should be 2.3 times.

[0208] In calculating the value, the ecological value for the two periods was first calculated using the adjusted coefficients: in 2010, the value of the northern shore mudflats increased by 8.6 million yuan compared to 2000 due to the expansion of the area and the increase in the coefficients; although the area of ​​the southern shore mangrove forest decreased by 5 hectares, the loss was partially offset by the increase in the coefficients. Finally, the net profit and loss were output by combining the attenuation factor: due to high-intensity development and the sensitive location, the actual loss of the southern shore reclamation area reached 1.8 times the value assessed by the traditional method.

[0209] The findings, confirmed by the local marine bureau, indicate that traditional static assessments prioritized the restoration of the north shore, while this invention reveals that the core issue lies in the high-intensity reclamation on the south shore. Following the redistribution of restoration funds, the recovery efficiency of seagrass beds on the south shore significantly improved, reaching a 40% increase. This further demonstrates the scientific validity and effectiveness of the dynamic correction mechanism in seagrass bed ecological restoration.

[0210] Furthermore, the construction of the interference attenuation module includes the following steps:

[0211] a) Identify the type of development activity: Based on the patch characteristics of artificial land expansion areas, including patch edge density and number of patches, determine the type of development activity;

[0212] If the density of patch edges is greater than 100 meters per hectare and the number of patches is greater than 10, it is determined to be port construction.

[0213] If the number of patches is 1 and the area is greater than 1 hectare, it is determined to be land reclamation.

[0214] The rest were classified as tourist facilities;

[0215] b) Obtain disturbance duration data, specifically the time span from the start of development activities to the ecological cost-benefit assessment point, denoted as duration in years. T ;

[0216] c) Determine the type-based attenuation coefficient: Based on the determined development activity type, select the corresponding basic attenuation coefficient from the preset type-based attenuation coefficient table. ;

[0217] Types of land reclamation: ;

[0218] Port construction types: ;

[0219] Types of tourist facilities: ;

[0220] d) Calculate the time correction factor: based on the duration in years. T Calculate the time correction factor :

[0221] ;

[0222] e) Determine the sensitive area correction factor: Obtain spatial data of ecologically sensitive areas, including mangrove distribution boundaries, seagrass bed core areas, and the scope of national coastal protected areas, and determine whether artificial land expansion areas spatially overlap with any ecologically sensitive area:

[0223] If spatial overlap exists, a sensitive area correction factor is set. ;

[0224] If there is no spatial overlap, then set ;

[0225] f) Calculate the final decay factor: based on , , The final output of the interference attenuation module is an attenuation factor used to calculate the applicability of ecological gain / loss values. :

[0226] .

[0227] The prior art only sets a fixed attenuation coefficient based on the developed area. For example, a certain study uniformly uses an attenuation value of 0.7 for reclamation projects. This simplified processing method fails to fully consider the cumulative effects of different development stages and the protection priorities of ecological sensitive areas, resulting in the following two misjudgments: First, a new port in operation for 3 years and an old dock abandoned for 10 years obtain the same attenuation value, ignoring the natural recovery process; Second, a resort occupying the core area of mangroves and peripheral facilities are treated equally, underestimating the loss in biodiversity hotspots.

[0228] This technology implements a differential correction process:

[0229] Identification of development type: By analyzing the patches in the reclamation area, when the area of a single patch reaches 18 hectares and there are no secondary patches, it is determined as the reclamation type, and its basic attenuation coefficient is set to 0.72. In contrast, a certain prior art misjudged adjacent 12 - hectare tourism facilities (actually composed of 6 scattered patches) as the reclamation type and wrongly assigned an attenuation value of 0.73 (it should be 0.88 for tourism).

[0230] Duration correction: The project lasted 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 did not incorporate the time dimension and treated a newly built port and a historical dock equally.

[0231] Detection of sensitive area superposition: Spatial analysis shows that the reclamation area overlaps with the core protected area of mangroves by 9.2 hectares, and the sensitive area correction factor is set to 0.87. The prior art only roughly judges based on the buffer zone of the protected area boundary and fails to accurately identify the actual overlapping area.

[0232] Synthesis of the final attenuation factor: Multiply the basic coefficient 0.72, the time factor 0.92, and the sensitive factor 0.87 to obtain an intermediate value of 0.577. After being constrained by the [0.7, 0.9] interval, the final attenuation factor 0.7 is output. If the traditional method directly takes the standard value of 0.7 for the reclamation type, although the result value is the same, it lacks the quantitative response to the 7 - year recovery period and the core area overlap.

[0233] Revaluation of ecological value loss: The original loss value calculated by the module is 38 million yuan. After being corrected by the attenuation factor of 0.7, the loss value becomes 26.6 million yuan. Due to the failure to distinguish development attributes, the traditional method uses an attenuation factor of 0.85 for tourism facilities in the same area, underestimating the loss to 32.3 million yuan. The historical compensation records of the Ocean Bureau show that the actual ecological restoration investment in this area is 29.1 million yuan, which is consistent with the evaluation result of this invention. Practice has proved that the multi - dimensional attenuation mechanism significantly improves the accuracy of damage quantification and provides a reliable basis for formulating ecological compensation standards.

[0234] Example 2

[0235] A remote sensing quantification and ecological loss assessment device for coastal zone human activity disturbance is disclosed. This device comprises six interconnected core modules, achieving fully automated processing through hardware and software collaboration. The data acquisition module is equipped with a satellite sensor containing a photovoltaic conversion unit to acquire two-stage multispectral remote sensing images (spatial resolution 10–30 meters, including blue, green, red, near-infrared, and shortwave infrared bands) of the target coastal zone at 5–10-year intervals. Simultaneously, it collects tidal model data, digital elevation model data, and geological background data to ensure the integrity and timeliness of the source data.

[0236] After receiving the raw image, the preprocessing module first accelerates radiometric calibration using a GPU parallel computing architecture. Based on the sensor's nonlinear response characteristics, a gain compensation coefficient of 1.15 is applied to low-radiance regions to enhance the signal-to-noise ratio, while a suppression coefficient of 0.9 is used to avoid signal saturation in high-radiance regions. Furthermore, a polarization filtering algorithm is implemented in the near-infrared band to reduce the impact of sea surface flares. Subsequently, real-time aerosol optical thickness and humidity data are integrated, and atmospheric correction is performed using a radiative transfer model, outputting an accurate surface reflectance dataset. The processing efficiency is 15 times higher than that of traditional CPU solutions.

[0237] The interference index calculation module incorporates a U-Net deep learning model and implements an object-oriented segmentation strategy for the preprocessed reflectance data from the two periods: using 20 pixels as the segmentation benchmark, it integrates NDBI, NDVI, MNDWI, and band mean features to achieve accurate land use classification; simultaneously, it dynamically identifies temporary flooded areas using digital elevation data and a tidal model, and intelligently identifies building shadow areas using spectral criteria (i.e., near-infrared / red band ratio less than 1.2 and red-blue reflectance characteristics less than 0.2), replacing abnormal reflectance with the mean of similar land features in adjacent non-shadowed areas, and finally outputting the shadow-compensated NDBI and land use change map of the non-flooded area.

[0238] The disturbance quantification module extracts artificial land expansion areas, statistically analyzes the changed area, and merges the patch density index (weight 35%) and edge density index (weight 65%) to generate a fragmentation disturbance index. This module can be combined with a random forest regression model to predict intensity changes based on patch characteristics and environmental parameters, improving quantification reliability. When the index > 0.5, the disturbance enhancement model is activated—based on geological background data (such as the proportion of bedrock shoreline) and the distribution of ecologically sensitive areas (such as mangrove boundaries), a regional correction coefficient matrix is ​​generated, spatially weighting the NDBI difference to output the final intensity change value.

[0239] The ecological assessment module employs the unit area value equivalent factor method to construct a dynamic assessment model. This model adjusts and optimizes the ecological value coefficient in real time based on key parameters such as intertidal zone width and fish species abundance. It also calculates attenuation factors linked to development type (reclamation / port / tourism), duration, and overlapping of ecologically sensitive areas through a disturbance attenuation module. A random forest-based computational model can be introduced during the assessment process to optimize value coefficient correction and attenuation factor calculation, enhancing the model's adaptability to nonlinear ecological responses. Assessment results are uploaded in real time to the government's coastal zone monitoring platform via a dedicated interface, triggering an automatic ecological compensation calculation process.

[0240] The output module can generate spatial distribution maps and numerical reports of ecological gains and losses with a 30-meter grid precision, supporting three-dimensional dynamic risk simulation. By overlaying geologically sensitive area boundaries with tidal inundation models, this system can extrapolate long-term ecological losses under sea-level rise scenarios in a time series manner, providing visualized decision support for restoration projects. For example, based on a new model for predicting the effects of mangrove ecological reconstruction, the cumulative impact of reclamation areas on the survival rate of mangrove larvae over the next 30 years can be simulated, and early warnings can be issued for high-risk areas where the increase in losses may exceed 40%.

[0241] Although the embodiments of the present invention have been disclosed above, they are not limited to the applications listed in the specification and embodiments. They can be applied to various fields suitable for the present invention. For those skilled in the art, other modifications can be easily made. Therefore, without departing from the general concept defined by the claims and their equivalents, the present invention is not limited to the specific details.

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, Includes the following steps: S1: Acquire multispectral remote sensing images of the target coastal area at a first time point and a second time point, wherein the first time point is earlier than the second time point, and the time interval between the first time point and the second time point is between 5 and 10 years. The spatial resolution of the multispectral remote sensing images is between 10 meters and 30 meters. The multispectral remote sensing images include blue light band, green light band, red light band, near-infrared band and short-wave infrared band. At the same time, acquire tidal model data, digital elevation model data and geological background data. S2: Preprocess the multispectral remote sensing image, including radiometric calibration and atmospheric correction, to eliminate sensor errors and atmospheric scattering effects; S3: Based on the preprocessed multispectral remote sensing image, calculate the human activity disturbance index, which includes the normalized building index and land use type change detection. The normalized building index is calculated using the ratio of reflectance of the shortwave infrared band to the near-infrared band. The land use type change detection is achieved by comparing the classification results of the first time point and the second time point. The land use type change detection needs to dynamically mark the temporary flooded area based on the tidal model data and digital elevation model data, and exclude the area from the statistics of change detection. S4: Quantify the degree of human activity interference, including calculating the area and intensity change value of the changed area. The area of ​​the changed area is based on the output of land use type change detection, and the intensity change value is calculated based on the difference of the normalized building index. The intensity change value is corrected by the interference enhancement model and the final intensity change value is output. The interference enhancement model is based on the fragmentation interference index and the geological and ecological data to generate a regional correction coefficient matrix, and spatially weights the intensity change values. S5: The changes in area and intensity of the changed region are input into the ecosystem service value model, which uses the unit area value equivalent factor method; S6: The output of the ecosystem service value model includes spatial distribution maps and numerical reports of ecological cost-benefit assessment results.

2. The method for remote sensing quantification of human activity disturbance and ecological loss assessment in coastal zones according to claim 1, characterized in that, The preprocessing specifically involves: radiometric calibration of the raw digital quantization values ​​of the multispectral remote sensing image to output radiance data at the sensor entrance pupil; atmospheric correction of the radiance data to output surface reflectance data.

3. The method for remote sensing quantification of human activity disturbance and ecological loss assessment in coastal zones according to claim 1, characterized in that, Radiation calibration includes the following steps: a) Read the radiometric calibration coefficient file provided by the sensor manufacturer, and convert the original digital quantization values ​​into radiance data based on the gain coefficient and offset coefficient; b) Compensate the sensor nonlinear response for the converted radiance data: when the radiance value is less than 20% of the sensor saturation radiance value, use the low-radiance region gain compensation coefficient to improve the signal-to-noise ratio; when the radiance value is greater than 50% of the sensor saturation radiance value, use the high-radiance region gain compensation coefficient to suppress the saturation effect. c) Add specular reflection suppression processing to the near-infrared band radiance data separately, and use polarization filtering algorithm combined with solar altitude angle, sensor observation angle, red band reflectivity characteristics and sea surface wind speed data to correct flare interference.

4. The method for remote sensing quantification of human activity disturbance and ecological loss assessment in coastal zones according to claim 1, characterized in that, Atmospheric correction includes the following steps: a) Acquire aerosol optical thickness data and relative humidity data that match the imaging time of the multispectral remote sensing image. The aerosol optical thickness data is derived from satellite remote sensing inversion products, with a spatial resolution consistent with the multispectral remote sensing image. The relative humidity data is derived from meteorological reanalysis data, with a temporal resolution less than or equal to 6 hours. b) Based on aerosol optical thickness data and relative humidity data, calculate the aerosol scattering characteristic parameters at a wavelength of 550 nm. The calculation uses an aerosol scattering model that includes the mixing ratio parameters of sea salt, dust and anthropogenic pollution aerosols. The mixing ratio parameters are set according to historical observation data of the target coastal area. c) Convert the radiometrically calibrated radiance data into apparent reflectance data. The conversion formula is: ;in, This is the radiance value. d This is the Earth-Sun distance correction factor. Solar irradiance outside the atmosphere. The solar zenith angle; d) Atmospheric correction is performed using a radiative transfer model, the input parameters of which include aerosol scattering characteristics, sensor imaging geometry, surface elevation data, and apparent reflectivity. e) In the radiative transfer model, when the relative humidity data exceeds 60%, a humidity correction module is introduced. This module uses Mie scattering theory to calculate the correction factor for the water vapor's effect on the aerosol particle-scale expansion effect, so as to improve the model's accurate simulation of the influence of humidity.

5. The method for remote sensing quantification of human activity disturbance and ecological loss assessment in coastal zones according to claim 1, characterized in that, Calculating the human activity disturbance index includes the following steps: a) Land use classification was performed on the surface reflectance data at the first and second time points. An object-oriented segmentation algorithm was used for classification, and the segmentation scale parameter was set between 10 and 30 pixels. The selected features included the normalized building index, the normalized vegetation index, the improved normalized water index, and the mean of the surface reflectance band. Normalized Building Index ; Normalized Difference Vegetation Index ; Improved Normalized Water Index ; in, The surface reflectance in the blue light band. The surface reflectance is in the green light band. For shortwave infrared band surface reflectance, Near-infrared surface reflectance, The surface reflectance is in the red band. b) Based on the regions where the improved normalized water index value is greater than the water body determination threshold in the classification results, and combined with the tidal height data estimated by the tidal model, the following steps are performed: When the tide height is 1 meter higher than the local mean sea level, the area with an elevation lower than the tide height and classified as a body of water is marked as a temporary flooding zone. Elevation data are derived from digital elevation models; temporary inundation areas are not included in land use type change monitoring. c) In non-submerged and non-water body areas, determine the building shadow area based on the spectral characteristics of the vegetation cover area. If the vegetation area is within... and At that time, it was determined to be a building shadow area; The reflectance of the shaded area is replaced by the average reflectance value of the same type of ground features in adjacent non-shaded areas. d) Recalculate NDBI based on the reflectance data after shading compensation; e) Land use type change detection is achieved by comparing the classification results at the first and second time points, and only the pixels in non-temporary flooded areas that have undergone type transformation are counted.

6. The method for remote sensing quantification of human activity disturbance and ecological loss assessment in coastal zones according to claim 1, characterized in that, Quantifying the degree of human interference includes the following steps: a) Based on the monitoring results of land use type change, extract the artificial surface expansion area and calculate the total area of ​​the changed area; b) Conduct 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 by dividing the number of patches by the total area, and calculate the edge density index by dividing the total length of patch edges by the total area. c) The patch density index and the edge density index are fused with weights of 30% to 40% and 60% to 70%, respectively, to obtain the fragmentation interference index; d) Calculate the difference between the normalized building index at the first time point and the second time point as the base strength change value; e) Input the base intensity change value into the disturbance enhancement model, and the disturbance enhancement model will execute: When the fragmentation interference index is greater than the threshold of 0.5, the output enhancement intensity change value is multiplied by the basic intensity change value and the correction coefficient. If the fragmentation interference index does not exceed the threshold of 0.5, the output enhancement intensity change value will be consistent with the basic intensity change value. f) Output the enhanced strength change value as the final strength change value.

7. The method for remote sensing quantification of human activity disturbance and ecological loss assessment in coastal zones according to claim 6, characterized in that, The construction of the interference enhancement model includes the following steps: a) Obtain geological background data of the target coastal zone, including the proportion of bedrock coastline length and the spatial distribution of tidal flat sediment types, and extract spatial data of ecologically sensitive areas, including mangrove distribution boundaries, seagrass bed core areas and the scope of national coastal protected areas; b) Generate a regional correction coefficient matrix based on geological background data and ecologically sensitive area data. The value ranges from 0.8 to 1.

3. This matrix is ​​generated by coupling the fragmentation disturbance index with geological and ecological data. c) The change in reinforcement strength Weighted by the regional interference weight matrix: ； d) Output As the final intensity change value.

8. The method for remote sensing quantification of human activity disturbance and ecological loss assessment in coastal zones according to claim 1, characterized in that, The construction of the ecosystem service value model includes the following steps: a) Establish a basic value equivalent table, including the annual service value coefficient per unit area for five land types: mangroves, salt marshes, tidal flats, shallow seas, and artificial surfaces. ,in i For land category index, i =1 represents mangrove forest. i =2 indicates a salt marsh. i =3 indicates mudflats. i =4 indicates shallow sea. i =5 indicates an artificial surface; b) Collect dynamic parameters of the coastal zone geographic environment: Obtain the following spatial parameters for the target coastal zone region: Intertidal zone width data; Aquatic biodiversity index, which is the abundance of fish species calculated based on fishery resource survey data; Historical frequency data of storm surges; c) Generate the coastal zone spatial correction coefficient matrix: Based on the parameters collected in step b), perform... Perform dynamic correction and calculate the spatially corrected first... i Land type annual service value coefficient per unit area ; When the intertidal zone width increases by more than 500 meters from the baseline value, the tidal flat type... ; For every 10 fish species added to the abundance, the abundance of shallow marine species... ; When storm surges occur more than 0.5 times per year, mangrove and salt marshland types... ； ； For land categories that do not meet the correction conditions, ; d) Threshold attenuation calculation driven by interference intensity: For artificial land expansion areas, the attenuation factor is calculated using the interference attenuation module based on the change in enhancement intensity. : When the final intensity change value is greater than the threshold of 0.3, ; When the final intensity change is less than or equal to the threshold of 0.3, ; e) Calculate the ecological gain / loss value: Calculate the change in the value of ecosystem services using the following formula: ;in, i For land category index, i =1: Mangrove forest, 2: Salt marsh, 3: Mudflats, 4: Shallow sea, 5: Artificial surface; For the first time point i Land type area; For the second time point i Land type area; For the spatially corrected first i Annual service value coefficient per unit area for each land category; As a decay factor, the change in ecosystem service value caused by artificial land expansion is negative.

9. The method for remote sensing quantification of human activity disturbance in coastal zones and assessment of ecological losses and benefits according to claim 8, characterized in that, The construction of the interference attenuation module includes the following steps: a) Identify the type of development activity: Based on the patch characteristics of artificial land expansion areas, including patch edge density and number of patches, determine the type of development activity; If the density of patch edges is greater than 100 meters per hectare and the number of patches is greater than 10, it is determined to be port construction. If the number of patches is 1 and the area is greater than 1 hectare, it is determined to be land reclamation. The rest were classified as tourist facilities; b) Obtain disturbance duration data, specifically the time span from the start of development activities to the ecological cost-benefit assessment point, denoted as duration in years. T ; c) Determine the type-based attenuation coefficient: Based on the determined development activity type, select the corresponding basic attenuation coefficient from the preset type-based attenuation coefficient table. ; Types of land reclamation: ; Port construction types: ; Types of tourist facilities: ; d) Calculate the time correction factor: based on the duration in years. T Calculate the time correction factor : ； e) Determine the sensitive area correction factor: Obtain spatial data of ecologically sensitive areas, including mangrove distribution boundaries, seagrass bed core areas, and the scope of national coastal protected areas, and determine whether artificial land expansion areas spatially overlap with any ecologically sensitive area: If spatial overlap exists, a sensitive area correction factor is set. ; If there is no spatial overlap, then set ; f) Calculate the final decay factor: based on , , The final output of the interference attenuation module is an attenuation factor used to calculate the applicability of ecological gain / loss values. : 。 10. A remote sensing quantification and ecological loss assessment device for human activity disturbance in coastal zones, characterized in that, include: The data acquisition module is configured to acquire multispectral remote sensing images of the target coastal zone at the first and second time points. The spatial resolution of the multispectral remote sensing images is between 10 and 30 meters and includes blue light band, green light band, red light band, near-infrared band and short-wave infrared band. Simultaneously, it acquires tidal model data, digital elevation model data and geological background data. The preprocessing module, connected to the data acquisition module, is configured to perform radiometric calibration and atmospheric correction on the multispectral remote sensing image; Among them, the radiation calibration unit is based on the sensor nonlinear response compensation mechanism. It uses a gain compensation coefficient to improve the signal-to-noise ratio in the low radiation region, a gain compensation coefficient to suppress the saturation effect in the high radiation region, and performs specular reflection suppression driven by polarization filtering algorithm in the near-infrared band. The atmospheric correction unit integrates an aerosol scattering model and a humidity correction module. The humidity correction module generates an aerosol particle expansion effect correction factor based on Mie scattering theory when the relative humidity is >60%. The interference index calculation module, connected to the preprocessing module, is configured as follows: An object-oriented segmentation algorithm was used to classify land use based on two periods of surface reflectance data. The segmentation scale was 10-30 pixels, and the classification features included NDBI, NDVI, MNDWI and mean band reflectance. Temporary flooding areas are dynamically marked based on tidal models and digital elevation models, and building shadow areas are identified using spectral criteria. The reflectance of shadow areas is replaced by the average value of similar features in adjacent non-shaded areas. Output the land use type change detection results and the NDBI after shading compensation for non-temporary flooding areas; The interference quantization module, connected to the interference index calculation module, is configured as follows: The area of ​​artificially expanded land is extracted, and a fragmentation disturbance index is generated by fusing the patch density index and the edge density index. The NDBI difference is corrected by an interference enhancement model: when the fragmentation interference index > 0.5, the final intensity change value is output by weighting the regional correction coefficient matrix, which is generated by coupling geological background data and ecologically sensitive area data; The ecological assessment module, connected to the interference quantification module, is configured as follows: An ecosystem service value model was constructed using the unit area value equivalent factor method, and the basic value coefficient was spatially corrected by integrating dynamic parameters of the coastal zone. The attenuation factor is calculated based on the interference attenuation module. The module generates multiple maintenance positive factors by identifying the development type through patch features and combining the duration and the overlapping state of ecologically sensitive areas. According to the formula Output ecological gain / loss value; For the first time point i Land type area; For the second time point i Land type area; For the spatially corrected first i Annual service value coefficient per unit area for each land category; It is the attenuation factor; The results output module is connected to the ecological assessment module and is configured to generate spatial distribution maps and numerical ecological gain and loss assessment reports.

Images