Carbon sink baseline authenticity verification method and system based on multi-source remote sensing time series reconstruction

Abstract

The present application relates to remote sensing data processing and carbon sink accounting verification technical field, disclose a kind of carbon sink baseline authenticity verification method and system based on multi-source remote sensing time series reconstruction.The method extracts baseline year, surface type, project boundary and baseline accounting parameter from methodological declaration document;Obtain the multi-source satellite image in target time window, generate normalized time series data after effectiveness screening, pseudo-invariant feature extraction and cross-sensor normalization;Abnormal pixel is removed to time series data, climate zone adaptive smoothing and observation synthesis or missing reconstruction, obtain baseline year NDVI grid and plot baseline year NDVI;NDVI-land type mapping determined to infer surface type after climate zone calibration, and calculate baseline accounting parameter in combination with project type and regional parameter;Calculate land type consistency deviation, parameter deviation and baseline to the change consistency deviation of present situation, determine baseline authenticity verification grade according to preset grading rule, and output corresponding state marker signal and verification control instruction.

Description

Technical Field

[0001] This invention relates to the field of remote sensing data processing and carbon sequestration verification technology, and more specifically, to a method and system for verifying the authenticity of carbon sequestration baselines based on multi-source remote sensing time-series reconstruction. Background Technology

[0002] The baseline status of carbon sink projects is a crucial foundation for subsequent carbon sink increment accounting. For carbon sink projects targeting vegetation restoration, forest management, and grassland restoration, the baseline year's land cover status, vegetation growth status, and baseline accounting parameters such as biomass and carbon density directly affect the calculation results of subsequent carbon sink changes. Therefore, during project application or methodology filing, objectively verifying the baseline year, baseline land surface type, and related baseline accounting parameters recorded in the application materials is fundamental to ensuring the reliability of subsequent carbon sink accounting results.

[0003] In existing technologies, vegetation indices can be obtained from remote sensing images such as Landsat, Sentinel, and MODIS, and the Normalized Difference Vegetation Index (NDVI) can be used to reflect surface vegetation cover and growth status. However, different remote sensing satellites have differences in operating time, spatial resolution, revisit period, and sensor band response. Especially for earlier baseline years, it is often difficult to obtain temporally continuous and high-quality observation data from a single satellite data source. When observation data from different sensors are directly combined, differences in radiometric response between different sensors can also lead to systematic shifts in the NDVI time series, making it difficult to directly compare data from different sources.

[0004] On the other hand, remote sensing images are easily affected by factors such as clouds, cloud shadows, sensor malfunctions, and differences in observation periods. Within the target time window, some pixels in the same area may have insufficient observations or even be completely missing. Existing time-series smoothing or interpolation methods usually focus on generating continuous remote sensing images or analyzing long-term trends. When a uniform processing method is used for different climate zones, different effective observation numbers, and different degrees of missing data, it can easily lead to over-smoothing, outlier propagation, or unstable reconstruction results of missing pixels, thereby affecting the reliability of the baseline year's NDVI values.

[0005] Furthermore, existing NDVI data obtained through remote sensing processing are typically used primarily for vegetation monitoring, change detection, or mapping. The results lack a unified and automated conversion and cross-verification mechanism with parameters such as baseline land surface type, biomass, and carbon density in carbon sequestration project application materials. Vegetation cover status corresponding to the same NDVI may differ across different climate zones, and the carbon sequestration parameter conversion models and regional parameters used for different project types and ecological regions also vary. Therefore, directly using fixed NDVI intervals or unified conversion models for determination can easily lead to deviations in land surface type inference and baseline calculation parameter estimation.

[0006] Furthermore, existing carbon sink baseline verification typically relies on application materials, historical survey data, or manual sample plot data. There is a lack of a unified, mechanized cross-verification process between remote sensing reconstruction results and the declared land type, the calculated parameters of the declared baseline, and the carbon sink change trend from the baseline to the current status. It is difficult to form a complete technical link from historical remote sensing data reconstruction and parameter extrapolation to the determination of authenticity level.

[0007] In view of this, a method and system for verifying the authenticity of carbon sink baselines based on multi-source remote sensing time-series reconstruction is proposed. Summary of the Invention

[0008] This application proposes a method and system for verifying the authenticity of carbon sink baselines based on multi-source remote sensing time-series reconstruction, in order to solve one of the aforementioned technical problems.

[0009] The technical solution adopted in this application is as follows: In a first aspect, embodiments of this application provide a method for verifying the authenticity of carbon sink baselines based on multi-source remote sensing time-series reconstruction, including: Receive the methodology application document, extract the baseline year, baseline surface type, project vector boundary and baseline calculation parameters from the methodology application document, and generate a structured basic parameter object; Based on the baseline year and project vector boundary, multi-source satellite image metadata within the target time window is obtained. According to the image cloud coverage, the effective coverage ratio within the project vector boundary, the matching degree between the image acquisition time and the target time window, and the availability of sensor data, each screening index is compared with the corresponding preset validity conditions, and multi-source satellite images that meet the preset validity conditions are retained. Pseudo-invariant feature pixels are extracted from the retained images. Based on the pseudo-invariant feature pixels, a linear transformation relationship between different sensors is established, and the pixel observation values ​​of non-reference sensors are normalized to the same numerical reference to generate normalized time series data. Abnormal pixel identification and mask removal are performed on the normalized time series data. The time series smoothing parameters are determined from the preset parameter set according to the climate zone where the project is located. Based on the number of effective observations of the target pixel within the target time window, the corresponding observation value synthesis or missing pixel reconstruction is performed to generate the baseline year NDVI raster image. The plot-scale statistics of the baseline year NDVI raster image are performed on the project vector boundary to generate the plot baseline year NDVI. The preset NDVI-land type mapping interval is calibrated according to the climate zone where the project is located. The inferred land type is determined based on the calibrated NDVI-land type mapping interval and the baseline year NDVI raster image. The corresponding carbon sink parameter conversion model and regional parameters are determined based on the project type information obtained from the methodology application document and / or the inferred land type. The baseline accounting parameters are calculated based on the baseline year NDVI of the land parcel. The land use consistency deviation is calculated based on the NDVI of the baseline year of the land parcel and the NDVI interval corresponding to the surface type of the declared baseline. The parameter deviation is calculated based on the calculation parameters of the estimated baseline and the calculation parameters of the declared baseline. The carbon sink growth parameters are obtained from the methodological application document. Multi-source satellite images of the current year are obtained. The current year NDVI of the land parcel is obtained according to the same normalization, reconstruction and land parcel scale statistical rules as the baseline year. The change consistency deviation is determined based on the NDVI of the baseline year of the land parcel, the current year NDVI of the land parcel and the declared carbon sink growth parameters. The baseline authenticity verification level is determined according to the preset grading rules corresponding to each deviation, and the status mark signal and verification control command corresponding to the verification level are output.

[0010] Preferably, the baseline year, baseline surface type, project vector boundary, and baseline calculation parameters are extracted from the methodology application document, including: The method application document is checked to see if it contains a structured attachment that matches the preset data structure. If it does, the corresponding fields are extracted from the structured attachment according to the preset field mapping relationship. If the structured attachment does not exist or the structured attachment fails to be parsed, the layout of the methodological application document is analyzed and paragraphs are segmented. The corresponding fields are located by using preset regular expressions, keyword dictionaries and named entity recognition models, and the extraction results and extraction confidence of each field are output. For fields directly extracted from the structured attachments, they are used as fields to be verified; for fields automatically extracted from the main text of the methodology application document, they are used as fields to be verified, and their extraction confidence is compared with a preset extraction confidence threshold. When the extraction confidence of any field is lower than the preset extraction confidence threshold, a corresponding data correction request signal is generated and the returned data is received. When each field to be verified meets any of the following conditions: direct extraction of structured attachments, extraction confidence not lower than the preset extraction confidence threshold, or receiving feedback data that meets the preset field format for low confidence fields, the year range, numerical unit, geometric boundary, and parameter validity are verified for each current field. After the verification is passed, the basic parameter object is generated.

[0011] Preferably, the extraction of pseudo-invariant feature pixels includes: Candidate pixels are determined within the extended buffer of the project vector boundary, and time-series observation data of each candidate pixel are obtained within a preset cross-year time range including the baseline year. For each candidate pixel, calculate the temporal standard deviation σ_NIR of the near-infrared band reflectance, the difference between the temporal maximum and minimum NDVI values ​​Δ_NDVI, and the 95th percentile of the NDVI difference between adjacent time steps. Here, σ_NIR represents the degree of dispersion of the near-infrared band reflectance of the candidate pixel over time, and Δ_NDVI represents the magnitude of the NDVI change of the candidate pixel within the preset year-long time range. When σ_NIR < 0.015, Δ_NDVI < 0.10, and the 95th percentile < 0.05, the corresponding candidate pixel is determined as a pseudo-invariant feature pixel. When the number of pseudo-invariant feature pixels obtained is lower than the preset number of samples, the filtering thresholds corresponding to σ_NIR and Δ_NDVI are increased and the filtering is re-executed; when the number of pseudo-invariant feature pixels after re-filtering is still lower than the preset number of samples, a pseudo-invariant feature sample shortage status signal is output.

[0012] Preferably, the linear conversion relationship between different sensors and the supplementary auxiliary data include: Linear regression is performed based on paired observations of the same pseudo-invariant feature pixel in the reference sensor and the sensor to be converted to determine the first conversion coefficient k and the second conversion coefficient c, and a linear conversion relationship is established: NDVI baseline = k × NDVI to be converted + c; Wherein, NDVI reference represents the NDVI value after normalization to the reference sensor value reference, NDVI to be converted represents the NDVI observation value corresponding to the sensor to be converted, k represents the first conversion coefficient used to adjust the proportional relationship of NDVI values ​​of different sensors, and c represents the second conversion coefficient used to adjust the offset relationship of NDVI values ​​of different sensors. The linear transformation relationship that satisfies the preset fit condition and sample size condition is used to transform the time series observations of the corresponding sensor to be transformed, so as to generate time series data under a unified numerical benchmark. Before time series reconstruction, when the candidate images of the main data that meet the validity screening conditions are empty, or when the number of valid observations of the main sensor within the preset time range is lower than the preset number of scenes threshold, low spatial resolution auxiliary remote sensing data is acquired, and the auxiliary remote sensing data is supplemented into the time series data by directly supplementing or by fusing according to the preset fusion weight, depending on the triggering reason. After the initial temporal reconstruction is completed, when the proportion of pixels marked as unreconstructable exceeds a preset proportion threshold, the auxiliary remote sensing data is invoked to supplement the temporal data corresponding to the unreconstructable pixels, so that the reconstruction of the corresponding pixels can be performed again.

[0013] Preferably, time-series smoothing parameters are determined from a preset parameter set based on the climate zone where the project is located, including: The climate type of the project is determined based on the project vector boundary. A discrete candidate set of smoothing parameters λ corresponding to the climate type is obtained from the preset climate type-smoothing parameter mapping relationship. Here, λ represents the smoothing adjustment parameter that controls the relationship between the fidelity of time series observations and the smoothness of curves. Calculate the noise level of the target pixel temporal data. When the noise level is not lower than a preset high noise threshold, remove the minimum λ value from the discrete candidate set; when the noise level is not higher than a preset low noise threshold, remove the maximum λ value from the discrete candidate set. Weighted least squares temporal smoothing is performed using the remaining λ values, and the corresponding root mean square error of reconstruction is calculated through leave-one cross-validation. The λ value that minimizes the root mean square error of reconstruction is selected as the temporal smoothing parameter for the target pixel.

[0014] Preferably, the corresponding observation synthesis or missing pixel reconstruction is performed based on the number of effective observations of the target pixel within the target time window, including: The number of valid observations N for the target pixel within the baseline annual growing season window is counted, where N represents the number of valid NDVI observations of the target pixel after invalid observations are removed within the baseline annual growing season window. When N≥5, the median of the effective observations is taken as the baseline annual NDVI of the target pixel; When 3 ≤ N < 5, the arithmetic mean of the effective observations is taken as the baseline annual NDVI of the target pixel; When N=2, the difference between two valid observations is calculated. When the difference is not greater than a preset difference threshold, the average of the two valid observations is taken as the baseline annual NDVI of the target pixel. When the difference is greater than the preset difference threshold, the missing pixel reconstruction is initiated. When N=1, the single valid observation is used as the baseline annual NDVI of the target pixel, and a single observation synthesis status marker is set. When N=0, proceed to the missing pixel reconstruction stage; In the missing pixel reconstruction, effective neighboring pixels that meet the constraints of uncertainty, interpolation state, and NDVI difference are selected from the spatial neighborhood of the target pixel. When the number of effective neighboring pixels is not less than a preset number, the reconstructed value of the target pixel is calculated using an inverse distance weighting method. When the number of effective neighboring pixels is greater than zero but less than the preset number, the effective temporal observations of the years before and after the baseline year of the target pixel and the effective neighboring pixels are used as training data, and the reconstructed value of the target pixel is obtained through Gaussian process regression. When there are no effective neighboring pixels that meet the conditions, temporal smoothing extrapolation is performed based on the effective temporal observations of the years before and after the baseline year of the target pixel. When there are still no effective observations of the target pixel within a preset historical time range, the target pixel is set to an unreconstructable state.

[0015] Preferably, the preset NDVI-land type mapping interval is calibrated according to the climate zone where the project is located, and the inferred land type is determined, including: The climate zone where the project is located is determined based on the project vector boundary and the preset climate zone data, and the lower boundary correction coefficient and upper boundary correction coefficient of the corresponding NDVI interval are obtained. The lower and upper boundaries of the preset NDVI-land category mapping interval are calculated with the corresponding correction coefficients to form the regionally calibrated NDVI-land category mapping interval. The NDVI of each pixel in the baseline annual NDVI raster image is input into the NDVI-land class mapping interval after regional calibration for pixel-level classification. When the NDVI of the target pixel is within the boundary blur zone of the adjacent mapping interval, neighborhood determination and mixed pixel decomposition are performed on the target pixel to obtain at least one land surface type and its area ratio corresponding to the target pixel. Based on the classification results of each pixel within the project vector boundary, the area proportion of each surface type is calculated, and the land parcel-level inferred surface type is determined based on the area proportion.

[0016] Preferably, the corresponding carbon sink parameter conversion model and regional parameters are determined based on the project type information obtained from the methodology application document and / or the inferred land surface type, including: The project type is determined based on the project type keywords in the methodological application document; when no project type is identified, the project type is determined based on the area proportion of the inferred land surface type. When the project type corresponds to forestry carbon sinks, a spatial intersection operation is performed between the project vector boundary and the preset ecological forest area vector boundary. The corresponding ecological forest area is determined based on the intersection result, and the first biomass model parameter a and the second biomass model parameter b corresponding to the ecological forest area are obtained from the preset ecological forest area parameter set, according to: Estimated aboveground biomass = exp(a + b × plot baseline year NDVI) Calculate the extrapolated aboveground biomass, where a and b are regional model parameters determined based on the ecological forest area and used to characterize the mapping relationship between the baseline year NDVI and aboveground biomass of the plot; The estimated aboveground biomass is converted based on the carbon content corresponding to the inferred land surface type to obtain the estimated carbon density. When the project type corresponds to a carbon sink type using other biomass conversion models, the conversion model and regional parameters corresponding to the project type are retrieved from the preset project type parameter set, and the corresponding extrapolated baseline accounting parameters are calculated based on the baseline year NDVI of the plot. When the baseline calculation parameter to be verified includes soil carbon density, the estimated soil carbon density is determined from the preset soil carbon density parameter set based on the project surface type, climate zone, and available soil type information, and the corresponding confidence level is set according to the completeness of the input data.

[0017] Preferably, determining the baseline authenticity verification level includes: The land use consistency deviation δ_land is calculated using the following formula: δ_land=|NDVI_parcel-NDVI_mid| / W_land; Wherein, NDVI_parcel is the annual NDVI of the land parcel baseline, NDVI_mid is the midpoint of the NDVI interval after regional calibration corresponding to the declared baseline surface type, and W_land is the width of the NDVI interval after regional calibration. The units and calculation methods of the estimated baseline calculation parameters and the declared baseline calculation parameters are standardized, and the parameter deviation δ_param is calculated according to the following formula when the declared baseline calculation parameter is greater than zero: δ_param=|P_est-P_decl| / P_decl; Where P_est is the calculated baseline accounting parameter, and P_decl is the declared baseline accounting parameter corresponding to the calculated baseline accounting parameter; when P_decl is zero, the absolute deviation is compared with the corresponding preset absolute threshold. Based on the baseline year NDVI and the current year NDVI of the land parcel, the corresponding carbon density or carbon sink change parameters are calculated and converted into the same change indicators as the declared carbon sink growth parameters, resulting in the remote sensing change rate r_rs and the declared change rate r_decl. The change consistency deviation δ_change is then calculated according to the following formula: δ_change=|r_rs-r_decl| / |r_decl|; Wherein, r_rs represents the remote sensing carbon sink change rate determined by the conversion result of carbon sink parameters corresponding to the baseline year NDVI and the current year NDVI of the land parcel, and r_decl represents the declared carbon sink change rate determined by the declared carbon sink growth parameters in the methodology declaration document; when r_decl is zero, |r_rs-r_decl| is used as the absolute change deviation, and the absolute change deviation is compared with the corresponding preset absolute change threshold; The individual level of each verification dimension is determined based on the preset grading thresholds corresponding to land type consistency deviation, parameter deviation, and change consistency deviation; when the direction of remote sensing change is opposite to the direction of carbon sink growth reported, the change consistency dimension is determined as the highest anomaly level. The comprehensive verification level is determined according to the principle of prioritizing the highest anomaly level; when at least two verification dimensions reach the medium level, the comprehensive verification level is upgraded to the high level; when there are fewer than two valid verification dimensions, the comprehensive verification status is set to undetermined. Based on the comprehensive verification level, output at least one of the following: an abnormal status flag, a data correction request signal, a route review instruction, a verification termination signal, or a verification pass signal.

[0018] Secondly, embodiments of this application provide a carbon sink baseline authenticity verification system based on multi-source remote sensing time-series reconstruction to implement the method described in any of the above, comprising: The basic parameter extraction module is used to receive methodology application documents, extract baseline year, baseline surface type, project vector boundary and baseline calculation parameters from the methodology application documents, and generate a structured basic parameter object; The multi-source image normalization module is used to acquire multi-source satellite image metadata within the target time window based on the baseline year and project vector boundary. It compares each screening index with corresponding preset validity conditions based on image cloud coverage, effective coverage ratio within the project vector boundary, the matching degree between image acquisition time and the target time window, and sensor data availability. Multi-source satellite images that meet the preset validity conditions are retained. Pseudo-invariant feature pixels are extracted from the retained images. A linear transformation relationship between different sensors is established based on these pseudo-invariant feature pixels. Pixel observations from non-reference sensors are normalized to the same numerical reference, generating normalized time-series data. Furthermore, when performing change consistency verification, the current year's multi-source satellite images are processed according to the same validity screening and normalization rules to generate current year's normalized time-series data. The temporal reconstruction module is used to perform abnormal pixel identification and mask removal on the normalized temporal data, determine the temporal smoothing parameters according to the climate zone where the project is located, and perform corresponding observation value synthesis or missing pixel reconstruction based on the effective number of observations of the target pixel within the target time window to generate a baseline year NDVI raster image. It also performs plot-scale statistics on the baseline year NDVI raster image according to the project vector boundary to generate the plot baseline year NDVI; and processes the current year normalized temporal data according to the same temporal reconstruction and plot-scale statistics rules as the baseline year to generate the plot current year NDVI. The parameter conversion module is used to calibrate the preset NDVI-land type mapping interval according to the climate zone where the project is located, determine the inferred land type based on the calibrated NDVI-land type mapping interval and the baseline year NDVI raster image, and determine the corresponding carbon sink parameter conversion model and regional parameters based on the project type information and / or inferred land type obtained from the methodology application document, and calculate the baseline accounting parameters based on the baseline year NDVI of the plot. The authenticity verification module is used to determine the land type consistency deviation based on the annual NDVI of the land parcel baseline and the corresponding NDVI interval of the land surface type of the declared baseline, to determine the parameter deviation based on the calculated parameters of the estimated baseline and the calculated parameters of the declared baseline, and to determine the change consistency deviation based on the annual NDVI of the land parcel baseline, the annual NDVI of the current status of the land parcel, and the declared carbon sink growth parameters; to determine the baseline authenticity verification level according to the preset grading rules corresponding to each deviation, and to output the status flag signal and verification control command corresponding to the verification level.

[0019] This application discloses a method and system for verifying the authenticity of carbon sink baselines based on multi-source remote sensing time-series reconstruction. Compared with the prior art, this solution can produce at least the following technical effects.

[0020] First, this scheme acquires multi-source satellite imagery within the target time window based on the baseline year and project vector boundaries. Then, pseudo-invariant feature pixels are extracted through effectiveness screening. Based on these pseudo-invariant feature pixels, a linear transformation relationship is established between different sensors, normalizing pixel observations from non-reference sensors to the same numerical benchmark. This reduces the systematic NDVI shift caused by differences in band response and radiometric calibration between different satellite sensors, ensuring that observational data acquired from different periods and sensors have a unified numerical scale. This improves the comparability and combinatorial utilization of multi-source remote sensing data from historical baseline years.

[0021] Secondly, after generating normalized time-series data, this scheme identifies and masks out anomaly pixels, determines time-series smoothing parameters based on the climate zone of the project, and selects the corresponding observation synthesis or missing pixel reconstruction method based on the number of effective observations for the target pixel within the target time window. Therefore, compared to using uniform smoothing parameters or a single interpolation algorithm for all projects and pixels, this scheme can adjust the time-series processing strategy according to climate characteristics and the actual data completeness of pixels, reducing the impact of noise, cloud pollution, observational faults, and outliers on the baseline NDVI reconstruction results, and improving the stability of missing pixel recovery and baseline annual NDVI reconstruction.

[0022] Based on this, this scheme performs plot-scale statistics on the baseline year NDVI raster image according to the project vector boundary to form the plot baseline year NDVI, so that the pixel-level remote sensing reconstruction results can be further converted into quantitative characterization values ​​corresponding to the project plots, providing a unified data foundation for subsequent land use inference, carbon sink parameter conversion and authenticity verification.

[0023] Furthermore, this scheme calibrates the preset NDVI-land cover mapping interval based on the climate zone where the project is located, and determines the inferred land cover type based on the calibrated mapping interval. This can reduce the land cover identification bias caused by different vegetation cover states corresponding to the same NDVI due to different regional climate conditions. Simultaneously, based on the project type information and / or inferred land cover type obtained from the methodology application documents, the corresponding carbon sink parameter conversion model and regional parameters are selected. Baseline accounting parameters are calculated based on the baseline year NDVI of the land parcel, enabling the NDVI remote sensing observation results to be further converted into objective quantitative parameters corresponding to carbon sink accounting indicators such as land cover type, biomass, and carbon density, thereby improving the comparability between remote sensing data and carbon sink application parameters.

[0024] Furthermore, this scheme not only performs remote sensing reconstruction of the baseline year but also acquires multi-source satellite imagery of the current year and obtains the current year's NDVI according to the same normalization and reconstruction rules as the baseline year. Adopting unified data processing rules can reduce methodological biases caused by using different processing methods for the baseline and current years, enabling comparisons of changes in remote sensing results from the two periods under the same numerical scale and processing conditions.

[0025] Finally, this scheme calculates the land use consistency deviation between the baseline year NDVI and the declared baseline land surface type, the parameter deviation between the estimated baseline calculation parameters and the declared baseline calculation parameters, and the consistency deviation between the remote sensing change results of the baseline year and the current year relative to the declared carbon sink growth parameters. Based on preset grading rules, the baseline authenticity verification level is determined. This transforms the baseline authenticity judgment process, which previously relied primarily on manual data verification, into an objective calculation process based on multi-source remote sensing observations, time-series reconstruction results, and multi-dimensional parameter deviations. The system outputs machine-readable verification results through status marker signals and verification control commands, thus forming a complete technical closed loop: "multi-source remote sensing data normalization—baseline time-series reconstruction—land use and carbon sink parameter conversion—baseline and current status cross-verification—authenticity level output." Attached Figure Description

[0026] Various other advantages and benefits will become apparent to those skilled in the art upon reading the following detailed description of preferred embodiments. The accompanying drawings are for illustrative purposes only and are not intended to limit the invention. Furthermore, the same reference numerals denote the same parts throughout the drawings. In the drawings: Figure 1 A flowchart illustrating a method for verifying the authenticity of carbon sink baselines based on multi-source remote sensing temporal reconstruction, provided in an embodiment of the present invention. Figure 2 This is a schematic diagram of the structure of a carbon sink baseline authenticity verification system based on multi-source remote sensing time-series reconstruction, provided in an embodiment of the present invention. Detailed Implementation

[0027] Exemplary embodiments of the present disclosure will now be described in more detail with reference to the accompanying drawings. While exemplary embodiments of the present disclosure are shown in the drawings, it should be understood that the present disclosure may be implemented in various forms and should not be limited to the embodiments set forth herein. Rather, these embodiments are provided to enable a more thorough understanding of the present disclosure and to fully convey the scope of the disclosure to those skilled in the art. It should be noted that, unless otherwise specified, embodiments and features in the embodiments of the present invention can be combined with each other. The present invention will now be described in detail with reference to the accompanying drawings and embodiments.

[0028] like Figure 1 As shown in some embodiments of this application, this embodiment provides a method for verifying the authenticity of carbon sink baselines based on multi-source remote sensing time-series reconstruction. Specifically, the method includes the following steps: Step S101: Receive the methodology application document, extract the baseline year, baseline surface type, project vector boundary and baseline calculation parameters from the methodology application document, and generate a structured basic parameter object.

[0029] Specifically, in this step, the methodology application document serves as the data input for subsequent remote sensing data acquisition and baseline authenticity verification. The document text and attachments are parsed to obtain basic information for defining the baseline time, project spatial scope, and application baseline status, and these are then used to form basic parameter objects according to a pre-defined data structure.

[0030] The baseline year is used to determine the time range for retrieving subsequent historical remote sensing data; the baseline surface type is used to compare with the remotely inferred surface type; the project vector boundary is used to limit the spatial range for satellite image acquisition, pixel processing, and land parcel statistics; and the baseline calculation parameters are used to compare their consistency with the corresponding parameters obtained from subsequent remote sensing calculations. The generated basic parameter object serves as the data input for steps S102 to S105. The specific extraction and verification methods for the above-mentioned parameters will be further explained in subsequent implementation methods.

[0031] For example, after the above processing, a project obtains a baseline year of 2010, a baseline surface type of grassland, and obtains the corresponding project vector boundary and declared baseline calculation parameters, thus forming the basic parameter object of the project.

[0032] Step S102: Obtain multi-source satellite image metadata within the target time window based on the baseline year and project vector boundary. Compare each screening index with the corresponding preset validity conditions based on image cloud coverage, effective coverage ratio within the project vector boundary, matching degree between image acquisition time and target time window, and sensor data availability. Retain multi-source satellite images that meet the preset validity conditions. Extract pseudo-invariant feature pixels from the retained images. Establish linear transformation relationships between different sensors based on the pseudo-invariant feature pixels. Normalize the pixel observation values ​​of non-reference sensors to the same numerical reference to generate normalized time-series data.

[0033] Specifically, in this step, the target time window and corresponding satellite data sources that can provide effective observations are determined based on the baseline year, and candidate images intersecting with the project area are retrieved based on the project vector boundary. The candidate images are then screened for effectiveness, with screening indicators including image cloud coverage, effective coverage ratio within the project vector boundary, time distance between the image acquisition time and the baseline year, and sensor data availability status. Specifically, image cloud coverage is determined based on the cloud coverage field in the corresponding satellite image metadata or the corresponding quality control data; the effective coverage ratio within the project vector boundary refers to the ratio of the number of effective pixels remaining in the candidate image within the project vector boundary after removing clouds, cloud shadows, sensor anomalous stripes, and other invalid pixels to the total number of pixels within the project vector boundary; the matching degree between the image acquisition time and the target time window is determined based on the time distance between the image acquisition date and the baseline year or the target time window; and the sensor data availability status is used to determine whether the corresponding sensor has data available for subsequent processing in the image corresponding to the target time. In one implementation, the target time window is set as a range of two years before and after the baseline year. In this case, if the cloud coverage of the candidate image is not higher than 30%, the effective coverage ratio within the project vector boundary is not lower than 0.40, the image acquisition time is within the target time window (i.e., the time distance between the image acquisition time and the baseline year is not more than 730 days), and the corresponding sensor is available, the candidate image is determined to be an image that meets the validity conditions; if any condition is not met, the corresponding candidate image is removed from the current candidate set.

[0034] For the selected and retained multi-source images, pseudo-invariant feature pixels with relatively stable spectral responses across time periods are selected from the project area and its surroundings. A numerical transformation relationship between sensors is established using the corresponding observations of the same pseudo-invariant feature pixel in different sensors, and this transformation relationship is used to unify NDVI observations from different sources to the numerical scale of a preset reference sensor.

[0035] The observations that have undergone numerical normalization are organized according to the image acquisition time to obtain normalized time-series data, which is then used as input for anomaly handling, temporal smoothing, and missing data reconstruction in step S103. Specific selection criteria for pseudo-invariant feature pixels, cross-sensor conversion, and methods for supplementing auxiliary data when master data is insufficient will be further explained in subsequent embodiments.

[0036] Step S103: Perform abnormal pixel identification and mask removal on the normalized time series data, determine the time series smoothing parameters from the preset parameter set according to the climate zone where the project is located, and perform corresponding observation value synthesis or missing pixel reconstruction based on the effective observation quantity of the target pixel within the target time window to generate a baseline year NDVI raster image. Perform plot-scale statistics on the baseline year NDVI raster image according to the project vector boundary to generate the plot baseline year NDVI.

[0037] Specifically, this step first identifies clouds, cloud shadows, sensor anomalies, and other invalid observations in the normalized time series data, and removes these observations from the valid time series. Then, based on the project's climate zone and target pixel time series characteristics, the corresponding time series smoothing parameters are determined, and the remaining valid observations are time series smoothed.

[0038] For each target pixel within the project scope, the processing method is determined based on the number of valid observations within the target time window of the baseline year. When valid observations can directly characterize the vegetation status of the baseline year, statistical synthesis is performed on the valid observations; when valid observations are insufficient or cannot form reliable synthesis results, missing pixels are reconstructed by combining spatial neighborhood information and / or temporal information of adjacent years; pixels that cannot be reconstructed reliably are marked with corresponding status labels.

[0039] After processing each pixel, the obtained baseline year NDVI is used to form a baseline year NDVI raster image according to its spatial location. Then, plot-scale statistics are performed on the effective pixels using the project vector boundary to obtain the plot baseline year NDVI. The specific method for determining the temporal smoothing parameters and the synthesis and reconstruction rules corresponding to different numbers of effective observations will be further explained in subsequent implementation methods.

[0040] For example, after time-series processing, a project generates a baseline annual NDVI raster for 2010. The effective pixels within the project boundary are statistically analyzed, and the baseline annual NDVI of the land parcel is found to be 0.41. This result serves as the remote sensing input for subsequent land cover inference and baseline calculation parameter estimation.

[0041] Step S104: Calibrate the preset NDVI-land type mapping interval according to the climate zone where the project is located, determine the inferred land type based on the calibrated NDVI-land type mapping interval and the baseline year NDVI raster image; determine the corresponding carbon sink parameter conversion model and regional parameters based on the project type information obtained from the methodology application document and / or the inferred land type, and calculate the baseline accounting parameters based on the baseline year NDVI of the land parcel.

[0042] Specifically, vegetation NDVI levels vary across different climate zones. Therefore, the preset NDVI-land type mapping interval is modified according to the climate zone where the project is located. The effective pixels in the baseline year NDVI raster are matched with the modified mapping interval to form pixel-level land type results. Then, the spatial distribution of different land types within the project vector boundary is statistically analyzed to obtain the plot-level inferred land type.

[0043] Based on the methodological application information and inferred land surface type, the applicable carbon sink parameter conversion model is determined for the project. Regional parameters corresponding to this model are obtained based on the project's spatial location. The annual NDVI of the land parcel baseline is then converted into corresponding baseline biomass, carbon density, or other baseline accounting parameters requiring verification. For different types of carbon sink projects, parameter conversion relationships matching their project type are adopted to avoid using a uniform model to handle different vegetation types.

[0044] The specific climate zone calibration methods, boundary pixel processing, land parcel classification rules, and conversion models and regional parameters corresponding to different carbon sink projects will be further explained in subsequent implementation methods.

[0045] For example, if a project determines that a plot of land is mainly grassland based on the calibrated NDVI-land type mapping results, and further identifies it as a grassland carbon sink project, then the corresponding grassland carbon sink parameter conversion model and regional parameters are called to convert the plot's baseline year NDVI 0.41 into the corresponding extrapolated baseline accounting parameters.

[0046] Step S105: Calculate the land use consistency deviation based on the NDVI of the baseline year of the land parcel and the NDVI interval corresponding to the land surface type of the declared baseline, and calculate the parameter deviation based on the calculated parameters of the estimated baseline and the calculated parameters of the declared baseline; obtain the declared carbon sink growth parameters from the methodology application document, obtain multi-source satellite images of the current year, and obtain the current year NDVI of the land parcel according to the same normalization, reconstruction, and land parcel scale statistical rules as the baseline year; determine the change consistency deviation based on the NDVI of the baseline year of the land parcel, the current year NDVI of the land parcel, and the declared carbon sink growth parameters; determine the baseline authenticity verification level according to the preset grading rules corresponding to each deviation, and output the status marker signal and verification control command corresponding to the verification level.

[0047] Specifically, land cover consistency bias is used to evaluate the degree of consistency between the baseline vegetation status obtained from remote sensing reconstruction and the land surface type of the declared baseline; parameter bias is used to evaluate the numerical difference between the baseline calculation parameters obtained from remote sensing and the calculation parameters of the declared baseline.

[0048] For consistency verification, multi-source satellite imagery of the current year is acquired, and the current year's remote sensing results are obtained according to the same image normalization, temporal processing, and land parcel statistics rules as the baseline year. Remote sensing carbon sink change indicators are determined based on the remote sensing changes between the baseline year and the current year, and compared with the declared carbon sink growth parameters after conversion to the same caliber to obtain the consistency deviation.

[0049] Individual verification results are determined according to preset grading rules corresponding to the three verification dimensions, and the baseline authenticity verification level is formed by combining all effective verification dimensions. Based on this level, corresponding status marker signals and verification control commands are output to complete the processing loop from historical remote sensing reconstruction results to baseline authenticity verification results. The specific calculation methods for each deviation, grading thresholds, and comprehensive level determination rules are further explained in subsequent implementation methods.

[0050] For example, after completing the above three types of consistency checks for a project, land category consistency deviation, parameter deviation and change consistency deviation are generated respectively. Based on the corresponding classification rules, a comprehensive verification level is formed, and the system outputs corresponding control signals such as verification passed, data correction, review or termination of verification.

[0051] In some embodiments of this application, in step S101, to improve the adaptability of parameter extraction from methodological application materials of different formats, a hierarchical extraction method combining structured attachment parsing, text information recognition, and low-confidence correction can be adopted: This involves extracting baseline year, baseline surface type, project vector boundary, and baseline calculation parameters from the methodological application document, including: The method application document is checked to see if it contains a structured attachment that matches the preset data structure. If it does, the corresponding fields are extracted from the structured attachment according to the preset field mapping relationship. If the structured attachment does not exist or the structured attachment fails to be parsed, the layout of the methodological application document is analyzed and paragraphs are segmented. The corresponding fields are located by using preset regular expressions, keyword dictionaries and named entity recognition models, and the extraction results and extraction confidence of each field are output. For fields directly extracted from the structured attachments, they are used as fields to be verified; for fields automatically extracted from the main text of the methodology application document, they are used as fields to be verified, and their extraction confidence is compared with a preset extraction confidence threshold. When the extraction confidence of any field is lower than the preset extraction confidence threshold, a corresponding data correction request signal is generated and the returned data is received. When each field to be verified meets any of the following conditions: direct extraction of structured attachments, extraction confidence not lower than the preset extraction confidence threshold, or receiving feedback data that meets the preset field format for low confidence fields, the year range, numerical unit, geometric boundary, and parameter validity are verified for each current field. After the verification is passed, the basic parameter object is generated.

[0052] As described above, upon receiving the methodology application document, the system first checks whether the document includes structured attachments that match the preset data structure. Structured attachments can be data files with fixed field names and field hierarchies. When a corresponding attachment is detected, the attachment field names are matched against the preset field mapping relationship, and the field values ​​corresponding to the baseline year, baseline surface type, project vector boundary, and baseline calculation parameters are read respectively. Direct extraction using structured attachments can reduce the uncertainty in information recognition caused by the multiple expressions in natural language.

[0053] When no structured attachments are available, or when the structured attachments fail to be parsed, the layout of the methodological application document is analyzed and paragraphs are segmented to distinguish between titles, body text, tables, and paragraphs with different semantic meanings. Then, the candidate positions and candidate values ​​of the four types of target fields are determined by combining preset regular expressions, keyword dictionaries, and named entity recognition models, and each extraction result is assigned a corresponding extraction confidence level.

[0054] For the baseline year, priority should be given to identifying the four-digit year near keywords such as "baseline year," "benchmark year," and "benchmark year." When multiple years appear simultaneously in the text, the target value should be determined based on the semantic and positional relationships between the candidate years and the baseline year keywords, avoiding mistaking other years such as the project inclusion period or the application year as the baseline year. For example, if the text simultaneously states "the baseline year is 2010" and "the project inclusion period begins in 2015," then 2010 should be determined as the baseline year.

[0055] For baseline land surface types, the land category names in the text are matched against a pre-defined land category keyword dictionary, and different expressions are grouped into pre-defined land category categories. For example, for information with state or degree modifiers such as "degraded grassland" or "sparse woodland," the original description can be retained, and the corresponding standard land category can be determined based on the core land category terms. When the application materials contain multiple land surface types, the information for each land category is saved separately for subsequent NDVI interval matching based on the area proportion of each land category.

[0056] For project vector boundaries, when the attachments contain vector data that directly represents the project boundaries, the corresponding geometric objects are read; when the application materials only record the administrative location or latitude and longitude boundaries in the main text, the geographic coordinates are extracted from the corresponding location description to form the project spatial scope. After boundary extraction, the closure, area, and coordinate range of the geometric objects are verified, and the project area is converted to a unified area unit. For example, square kilometers can be converted to hectares to meet the unified calculation requirements for subsequent plot-scale statistics.

[0057] For baseline accounting parameters, the values ​​and units of measurement for parameters such as baseline biomass, baseline carbon density, and baseline soil carbon density are identified, and the values ​​are associated with and saved with their corresponding units. For parameters with the same meaning but different units of measurement, they are converted into a unified representation according to preset unit conversion relationships. For cases of abnormal values, unidentifiable units, or missing key parameters, corresponding abnormal statuses are set to prevent abnormal data from directly participating in subsequent deviation calculations.

[0058] After candidate fields are obtained through text recognition, the extraction confidence level of each field is compared with a preset extraction confidence threshold. The preset extraction confidence threshold can be set to, for example, 0.85. When the extraction confidence level of any key field is lower than the preset extraction confidence threshold, it indicates that the automatically extracted results are insufficient as reliable input for subsequent calculations. A data correction request signal corresponding to the low-confidence field is generated, and the corrected data returned by the reporting terminal for the corresponding field is received.

[0059] For structured corrected data returned using a preset field format, the corrected data is matched with the corresponding low-confidence field based on the field identifier, and the original extraction result is replaced with the corrected data. Since this corrected data is not obtained through automatic text recognition, no extraction confidence is set for it; instead, it directly enters the data validity verification process. For corrected data returned in unstructured text form that still requires automatic information extraction, the text recognition of the corresponding field is re-executed, and the extraction confidence of the field is recalculated. When the recalculated extraction confidence is still lower than the preset extraction confidence threshold, a corresponding data correction request signal is generated again.

[0060] When all fields meet the conditions for entering data validity verification, data validity verification is uniformly performed on all fields. Specifically, for fields directly extracted from structured attachments or structured correction fields returned using a preset field format, no extraction confidence level is set; for fields automatically identified through the main text, their extraction confidence level should not be lower than the preset extraction confidence threshold. For low-confidence fields returned as correction data in unstructured text form, the condition for entering data validity verification is that the extraction confidence level after re-identification is not lower than the preset extraction confidence threshold. Specifically, the baseline year is checked to ensure it falls within the allowed historical remote sensing data time range; the project boundaries are checked for geometric closure, area, and coordinate range; the baseline accounting parameters are checked for numerical symbols, units of measurement, and data integrity; and unit conversion is performed on data requiring unified measurement standards. If any verification fails, a corresponding data correction request is generated; after all verifications pass, all fields are organized into basic parameter objects.

[0061] The generated basic parameter object includes at least the baseline year, standardized baseline surface type, project boundary and area information, and the values ​​and units of each baseline calculation parameter. It can also associate data source information for each field, and for fields obtained through automatic identification, it associates corresponding extraction confidence levels. For fields directly extracted from structured attachments or obtained through structured correction, it associates the corresponding data source and correction status, facilitating data traceability in subsequent processing. This basic parameter object serves as the unified input for subsequent multi-source satellite image acquisition and baseline authenticity verification.

[0062] In some embodiments of this application, a further implementation of the pseudo-invariant feature pixel extraction in step S102 is provided: the extraction of pseudo-invariant feature pixels includes: Candidate pixels are determined within the extended buffer of the project vector boundary, and time-series observation data of each candidate pixel are obtained within a preset cross-year time range including the baseline year. For each candidate pixel, calculate the temporal standard deviation σ_NIR of the near-infrared band reflectance, the difference between the temporal maximum and minimum NDVI values ​​Δ_NDVI, and the 95th percentile of the NDVI difference between adjacent time steps. Here, σ_NIR represents the degree of dispersion of the near-infrared band reflectance of the candidate pixel over time, and Δ_NDVI represents the magnitude of the NDVI change of the candidate pixel within the preset year-long time range. When σ_NIR < 0.015, Δ_NDVI < 0.10, and the 95th percentile < 0.05, the corresponding candidate pixel is determined as a pseudo-invariant feature pixel. When the number of pseudo-invariant feature pixels obtained is lower than the preset number of samples, the filtering thresholds corresponding to σ_NIR and Δ_NDVI are increased and the filtering is re-executed; when the number of pseudo-invariant feature pixels after re-filtering is still lower than the preset number of samples, a pseudo-invariant feature sample shortage status signal is output.

[0063] As described above, firstly, an outer buffer is formed by extending a preset distance outward from the project vector boundary to filter candidate pixels, and the effective observations of each pixel within the buffer are statistically analyzed within a preset year-long time range. For any candidate pixel, the total number of observations covering that pixel within the preset year-long time range is denoted as N_total, and the number of observations that still yield effective near-infrared reflectance and NDVI values ​​after removing clouds, cloud shadows, sensor anomalies, and invalid band values ​​is denoted as N_valid. When N_total = 0, the corresponding pixel is directly removed from the candidate pixel set; when N_total > 0, the effective observation ratio R_cov of that pixel is calculated according to R_cov = N_valid / N_total. When R_cov is not lower than a preset effective observation ratio threshold, the corresponding pixel is retained as a candidate pixel; when R_cov is lower than the threshold, the corresponding pixel is removed from the candidate set. This condition avoids the inability to reliably calculate the temporal standard deviation of near-infrared reflectance, the long-term variation amplitude of NDVI, and the quantile of the NDVI difference between adjacent time steps due to an insufficient number of effective temporal observations. In one implementation, the project boundary can be extended outward by 5km, and R_cov≥0.80 can be used as a valid observation requirement for candidate pixels.

[0064] For each candidate pixel, cross-year time-series observation data including the baseline year is acquired. In one specific approach, the observation timeframe is extended forward and backward by 3 years, centered on the baseline year, thus forming a total observation timeframe of 7 years. Using cross-year data, rather than relying solely on the baseline year's single-year data, allows for the determination of whether candidate land features remain stable over a longer period, avoiding the misidentification of temporarily stable vegetation or other dynamic land features as pseudo-invariant feature pixels.

[0065] For each candidate pixel, calculate the temporal standard deviation σ_NIR of the near-infrared band reflectivity, the difference between the maximum and minimum temporal values ​​of NDVI Δ_NDVI, and the 95th percentile of the NDVI difference between adjacent time steps.

[0066] σ_NIR is used to reflect the dispersion of the near-infrared reflectance of a candidate pixel over time. The smaller the σ_NIR, the more stable the near-infrared reflectance of the pixel is at different observation times, making it more suitable as a numerical reference between different sensors.

[0067] Δ_NDVI reflects the overall change in NDVI of a candidate pixel over the entire yearly time span, and it can be determined by the difference between the maximum and minimum temporal NDVI values ​​of that pixel. If Δ_NDVI is small, it indicates that the candidate pixel has not undergone significant vegetation cover changes over the years; if Δ_NDVI is large, it indicates that the location may belong to a land feature with significant vegetation growth changes, land cover changes, or significant seasonal fluctuations, and is not suitable as a stable reference for cross-sensor normalization.

[0068] The 95th percentile of the NDVI difference between adjacent time steps is used to evaluate whether candidate pixels exhibit significant abrupt changes within a short period. Specifically, the absolute difference in NDVI between two adjacent time steps is calculated, and all adjacent differences are sorted to obtain the 95th percentile as an indicator of the degree of abrupt change. Compared to simply calculating long-term maximum and minimum values, this indicator can further identify short-term abrupt changes such as rapid changes in farmland phenology and intermittent water body changes.

[0069] After calculating the above three indicators, a joint screening process is performed on candidate pixels. When σ_NIR < 0.015, Δ_NDVI < 0.10, and the 95th percentile of the NDVI difference between adjacent time steps < 0.05, the candidate pixel is identified as a pseudo-invariant feature pixel. If all three conditions are met simultaneously, it indicates that the pixel possesses relatively stable near-infrared reflectance characteristics and does not exhibit significant multi-year NDVI variations or short-term NDVI abrupt changes; therefore, it can serve as a reference sample for establishing correspondences between observations from different sensors. If any condition is not met, the candidate pixel is removed from the pseudo-invariant feature sample set.

[0070] After screening all candidate pixels, the number of pseudo-invariant feature pixels is further counted and compared with a preset sample size. The preset sample size is used to ensure that there is sufficient data for subsequent cross-sensor linear relationships; for example, it can be set to 100.

[0071] When the number of pseudo-invariant feature pixels obtained from the initial screening reaches the preset sample size, the resulting pixel set is directly provided to the subsequent cross-sensor linear relationship establishment process. When the sample size is insufficient, the screening thresholds for some stability indicators are appropriately relaxed before re-screening. For example, the screening threshold for σ_NIR can be increased from 0.015 to 0.025, and the screening threshold for Δ_NDVI can be increased from 0.10 to 0.15, allowing some pixels with slightly lower stability but still reference value to enter the sample set.

[0072] After relaxing the threshold, the number of pseudo-invariant feature pixels is counted again. If the number of samples obtained from the re-screening reaches the preset number of samples, the subsequent cross-sensor normalization is performed using the newly formed sample set; if it is still lower than the preset number of samples, a pseudo-invariant feature sample shortage status signal is output, indicating that there are not enough reference pixels around the current project to establish a stable cross-sensor conversion relationship, and the conversion relationship is not directly forced to be established with insufficient samples.

[0073] Therefore, through the processing steps of "candidate region limitation - cross-year stability analysis - multi-index joint screening - sample quantity judgment - threshold adjustment - insufficient state output", a set of pseudo-invariant feature pixels that meet the requirements of subsequent processing in terms of both quantity and stability can be formed. This provides reference data for establishing linear transformation relationships between different satellite sensors based on paired pixels and reduces cross-sensor normalization errors introduced by changes in the reference ground objects themselves.

[0074] In some embodiments of this application, a further implementation of the cross-sensor conversion and auxiliary data supplementation in step S102 includes: the linear conversion relationship between different sensors and the auxiliary data supplementation. Linear regression is performed based on paired observations of the same pseudo-invariant feature pixel in the reference sensor and the sensor to be converted to determine the first conversion coefficient k and the second conversion coefficient c, and a linear conversion relationship is established: NDVI baseline = k × NDVI to be converted + c; Wherein, NDVI reference represents the NDVI value after normalization to the reference sensor value reference, NDVI to be converted represents the NDVI observation value corresponding to the sensor to be converted, k represents the first conversion coefficient used to adjust the proportional relationship of NDVI values ​​of different sensors, and c represents the second conversion coefficient used to adjust the offset relationship of NDVI values ​​of different sensors. The linear transformation relationship that satisfies the preset fit condition and sample size condition is used to transform the time series observations of the corresponding sensor to be transformed, so as to generate time series data under a unified numerical benchmark. Before time series reconstruction, when the candidate images of the main data that meet the validity screening conditions are empty, or when the number of valid observations of the main sensor within the preset time range is lower than the preset number of scenes threshold, low spatial resolution auxiliary remote sensing data is acquired, and the auxiliary remote sensing data is supplemented into the time series data by directly supplementing or by fusing according to the preset fusion weight, depending on the triggering reason. After the initial temporal reconstruction is completed, when the proportion of pixels marked as unreconstructable exceeds a preset proportion threshold, the auxiliary remote sensing data is invoked to supplement the temporal data corresponding to the unreconstructable pixels, so that the reconstruction of the corresponding pixels can be performed again.

[0075] As described above, from the candidate sensors that have passed the validity screening, the master sensor is determined based on spatial resolution, number of valid observations, and temporal continuity, and the valid image corresponding to the master sensor is used as the master data. In one embodiment, after obtaining pseudo-invariant feature pixels that meet the requirements, for the reference sensor and each sensor to be converted, the paired NDVI observation values ​​of the same pseudo-invariant feature pixel in both sensors are extracted, and linear regression is performed using multiple paired samples to determine the first conversion coefficient k and the second conversion coefficient c, establishing the following linear conversion relationship: NDVI reference = k × NDVI to be converted + c. Wherein, NDVI reference represents the NDVI value expressed according to the numerical scale of the reference sensor, and NDVI to be converted represents the original NDVI observation value obtained by the sensor to be converted; the first conversion coefficient k is used to correct the proportional difference between the NDVI values ​​of different sensors, and the second conversion coefficient c is used to correct the overall numerical offset between different sensors. Paired observation of the same pseudo-invariant feature pixel refers to valid observations obtained by the reference sensor and the sensor to be converted for the same spatial location, and the imaging time difference meets the preset time matching condition. By limiting the time difference of paired observations, the influence of changes in the real surface between two observations on the linear conversion parameters can be reduced.

[0076] Pseudo-invariant feature pixels are used as regression samples because these pixels exhibit relatively stable spectral responses over a multi-year time span. Therefore, the numerical differences of the same pixel across different sensors primarily reflect the differences in the sensor's own observational response, rather than actual changes in the Earth's surface. Regression based on a large number of such stable pixels yields transformation parameters used to correct for differences across sensor systems.

[0077] In one specific approach, the NDVI numerical scale corresponding to Landsat-8 / 9 can be used as a unified benchmark, and corresponding transformation relationships can be established for Landsat-5, Landsat-7, Sentinel-2, and MODIS used as auxiliary data. For data that uses the same numerical benchmark as the benchmark sensor, no additional transformation is required. When there are no direct observations meeting the pairing requirements between the sensor to be transformed and the benchmark sensor within the available time frame, a step-by-step transformation relationship can be established with intermediate sensors that have effective overlapping observations with both the sensor to be transformed and the benchmark sensor. Each level of transformation relationship uses paired observations of corresponding pseudo-invariant feature pixels for regression, satisfying preset goodness-of-fit and sample size conditions. When all levels of transformation relationships are valid, the observations of the sensor to be transformed are sequentially transformed to the intermediate sensor numerical benchmark, and then to the unified numerical benchmark of the final benchmark sensor. Paired observations are preferably from the same pseudo-invariant feature pixel location, and observation data with similar imaging times and that have passed validity screening are selected to reduce the impact of actual surface changes on the transformation parameters.

[0078] After establishing the linear transformation relationship, it is necessary to evaluate whether the transformation relationship has sufficient reliability. Specifically, this can be done by checking the goodness of fit of the linear regression and the number of paired samples of pseudo-invariant features involved in the regression. Only when both the goodness of fit and the number of samples meet the preset conditions will the current transformation relationship be used to process the time series observations of the sensor to be transformed.

[0079] For example, the regression determination coefficient R can be used 2 A minimum score of 0.85 and at least 100 valid paired samples are considered validity criteria. If the fit is insufficient or the sample size is inadequate, it indicates that the current transformation relationship cannot stably reflect the numerical correspondence between the two types of sensors. In this case, the pseudo-invariant feature candidate region can be expanded to obtain new samples, or an auxiliary data supplementation process can be initiated.

[0080] When the linear transformation relationship meets the requirements, the NDVI observations obtained by the sensor to be transformed within the target time range are substituted into the corresponding transformation relationship to unify the observation results of different sensors to the numerical scale of the reference sensor, and then rearranged according to the image acquisition time to form normalized time-series data. In this way, subsequent time-series smoothing and missing pixel reconstruction processing deal with NDVI observation sequences with consistent numerical calibers, rather than original multi-source observations with different systematic biases.

[0081] In one implementation, a master sensor is determined from candidate sensors that have passed validity screening. The master sensor can be determined based on a comprehensive assessment of the spatial resolution, number of effective observations, and temporal continuity of each candidate sensor within the project area, and the effective imagery corresponding to the master sensor is used as the master data. Imagery from other available sensors is used to fill temporal gaps in the master data or participate in cross-sensor normalization. When multiple candidate sensors meet the criteria, the sensor with a larger number of effective observations and higher spatial resolution can be preferentially selected as the master sensor. Furthermore, before time-series reconstruction, it is determined whether the master data can meet the basic observation requirements within the target time range. When the master data candidate images that meet the validity screening criteria are empty, or the number of effective observations of the master sensor within a preset time range is lower than a preset scene count threshold, low spatial resolution auxiliary remote sensing data is acquired to supplement the original time-series data.

[0082] Here, "empty candidate images for master data" means that no major satellite images that pass the validity screening were obtained within a preset time range near the baseline year; "insufficient number of effective observations from the master sensor" means that although some master data can be obtained, its temporal coverage is too sparse to support subsequent temporal smoothing and missing reconstruction. In one implementation, auxiliary data supplementation can be triggered when there is no master data that meets the conditions within two years before and after the baseline year, or when there are fewer than four effective observations from the master sensor within one year before and after the baseline year.

[0083] Auxiliary remote sensing data can be remote sensing products with good temporal continuity but lower spatial resolution than the master data, such as MODIS NDVI products. Since auxiliary data differs from master data such as Landsat in terms of spatial resolution and observation characteristics, it is also necessary to convert it to a unified NDVI numerical benchmark and perform spatial and temporal alignment before adding it to the time series.

[0084] For example, auxiliary data with a spatial resolution of 250m can be resampled to the same 30m spatial grid as the main data, and the auxiliary observations can be matched with the temporal positions of the main data according to the image time distance.

[0085] To reduce the substitutive impact of low spatial resolution auxiliary data on high spatial resolution master data, different data supplementation methods can be adopted depending on the reason for data insufficiency. When there is no high spatial resolution master data that meets the validity conditions within the target time range, auxiliary remote sensing data that has been numerically normalized and spatially and temporally aligned can be directly supplemented to the corresponding time series position to establish a base time series that can be used for subsequent reconstruction.

[0086] When the number of effective observations from the main sensor is insufficient, but other high spatial resolution images that meet the validity criteria still exist, auxiliary remote sensing data can be fused with existing high spatial resolution data according to preset fusion weights. For example, the weight of auxiliary remote sensing data can be set to 0.3 and the weight of high spatial resolution data to 0.7, so that high spatial resolution data remains the primary data source. For situations where the contribution of auxiliary data needs to be increased, the corresponding weights can be adjusted according to the data completeness. After completing the initial time-series reconstruction, the proportion of pixels marked as unreconstructable within the project area can be statistically analyzed. If the proportion of unreconstructable pixels is not higher than a preset threshold, the initial reconstruction result is used; if the proportion of unreconstructable pixels is higher than this threshold, it indicates that the existing high spatial resolution data is still insufficient to form reliable reconstruction results in some areas. In this case, auxiliary remote sensing data is called again, supplementing only the time series corresponding to the unreconstructable pixels.

[0087] In one specific approach, a 20% proportion of unreconstructable pixels can be set as the trigger condition. When this proportion exceeds 20%, auxiliary remote sensing observations for the corresponding pixel locations are extracted, and after numerical normalization and spatial and temporal alignment, they are added to the corresponding time series. Then, missing reconstruction is re-performed for these pixels without having to repeatedly introduce auxiliary data for pixels that can already be reliably reconstructed.

[0088] Therefore, through the processing of "pseudo-invariant feature pairing - linear regression - transformation quality judgment - unified numerical benchmark - data shortage triggering auxiliary supplementation - secondary supplementation of unreconstructable pixels", the numerical consistency and temporal integrity of multi-source satellite time series data can be improved while maintaining the priority use of high spatial resolution master data, providing stable data input for subsequent NDVI time series smoothing and baseline annual missing pixel reconstruction.

[0089] In some embodiments of this application, a further implementation of the determination of time-series smoothing parameters in step S103 includes: determining time-series smoothing parameters from a preset parameter set based on the climate zone where the project is located, including: The climate type of the project is determined based on the project vector boundary. A discrete candidate set of smoothing parameters λ corresponding to the climate type is obtained from the preset climate type-smoothing parameter mapping relationship. Here, λ represents the smoothing adjustment parameter that controls the relationship between the fidelity of time series observations and the smoothness of curves. Calculate the noise level of the target pixel temporal data. When the noise level is not lower than a preset high noise threshold, remove the minimum λ value from the discrete candidate set; when the noise level is not higher than a preset low noise threshold, remove the maximum λ value from the discrete candidate set. Weighted least squares temporal smoothing is performed using the remaining λ values, and the corresponding root mean square error of reconstruction is calculated through leave-one cross-validation. The λ value that minimizes the root mean square error of reconstruction is selected as the temporal smoothing parameter for the target pixel.

[0090] As described above, in one implementation, the climate type of the project is first determined based on the project vector boundary, and a candidate set of smoothing parameters λ is obtained from a pre-established climate type-smoothing parameter mapping relationship. λ is used to adjust the relationship between the fidelity of the original observations and the smoothness of the curve during the time series smoothing process: when λ increases, the suppression effect on time series fluctuations is enhanced; when λ decreases, more local variation characteristics in the original NDVI time series are preserved.

[0091] The climate type-smoothing parameter mapping relationship is pre-established based on the typical variation characteristics of vegetation NDVI time series under different climate types. Specifically, representative historical NDVI time series samples can be obtained for each climate type. The samples are preprocessed according to the same outlier removal and data normalization rules as in this application, and Whittaker time series smoothing is performed on multiple discrete λ candidate values. Leave-one-out cross-validation is performed using valid observations in the samples to calculate the root mean square error of reconstruction corresponding to different λ values, thereby obtaining the preferred λ values ​​for each representative sample. The preferred λ distribution of multiple representative samples under the same climate type is statistically analyzed, and several discrete λ values ​​that can cover the main preferred range of the climate type are selected to form a λ candidate set corresponding to the climate type. Different climate types and their corresponding λ candidate sets are associated and stored to form a pre-defined climate type-smoothing parameter mapping relationship.

[0092] The candidate λ values ​​for different climate types are set based on the difference between the actual temporal variation of vegetation and short-term random fluctuations. A larger λ value enhances the suppression of short-term fluctuations, while a smaller λ value preserves more local temporal variations. Therefore, for arid regions with relatively sparse vegetation cover and significant short-term observation noise, a larger overall λ candidate value can be configured; for humid forest regions with high vegetation cover and where phenological characteristics need to be preserved, a smaller overall λ candidate value can be configured. For other climate types, the corresponding discrete candidate sets can be determined using the same historical sample evaluation method. For example, in one implementation, the λ candidate set for arid regions can be set to {50, 100, 200}, and the λ candidate set for humid forest regions can be set to {1, 2, 5, 10}. This mapping relationship is only used for preliminary limitation of the λ search range at the climate type level; the final λ value used for the target pixel is still determined through subsequent noise level screening and leave-one-out cross-validation.

[0093] After obtaining the initial candidate set corresponding to the climate type, the noise level of the target pixel's NDVI time series data is further calculated. In one implementation, the effective NDVI observations of the target pixel are sorted according to the observation time, and an initial time series trend is formed using a sliding median process independent of the selected smoothing parameter λ. For the i-th effective observation in the middle of the time series, the median of the (i-1), i-th, and (i+1)-th effective NDVI observations is used as the trend value corresponding to the i-th observation time; for the first effective observation at the beginning of the time series, the arithmetic mean of the first and second effective NDVI observations is used as the corresponding trend value; for the last effective observation at the end of the time series, the arithmetic mean of the last and second-to-last effective NDVI observations is used as the corresponding trend value. Thus, the initial trend value corresponding to each effective observation time is obtained. For each effective observation, the difference between its NDVI observation value and the corresponding initial trend value is calculated to obtain the time series residual, and the standard deviation of all time series residuals is calculated as the time series noise level of the target pixel. This initial trend is only used to evaluate short-term random fluctuations in the observation sequence and is not used as the final NDVI smoothing result, nor is it involved in the determination of the final smoothing parameter λ. Noise level can be used to reflect the degree of fluctuation of actual observed data relative to the overall time-series change trend, and serves as a basis for further adjusting the smoothing strength. When the noise level is not lower than the preset high noise threshold, it indicates that the current pixel contains significant short-term random fluctuations. In this case, the smallest λ value is removed from the candidate set to avoid selecting parameters with weak smoothing ability. When the noise level is not higher than the preset low noise threshold, it indicates that the observation sequence itself is relatively stable. In this case, the largest λ value is removed to reduce the possibility of over-smoothing weakening the details of real vegetation changes. For example, when the time-series noise level is ≥0.08, it is determined to be a high noise case; when the time-series noise level is ≤0.03, it is determined to be a low noise case; when it is between the two, the original λ candidate set is maintained.

[0094] For example, a lower λ candidate set {1, 2, 5, 10} can be set for humid forest areas; for arid areas, a λ candidate set {50, 100, 200} can be set. Different candidate ranges are used for different climatic zones to adapt to the differences in the temporal variation of vegetation and noise characteristics in different regions. For target pixels in arid areas, the initial λ candidate set is {50, 100, 200}. If the calculated temporal noise level reaches the high noise condition, the minimum value of 50 is removed, and {100, 200} is retained as subsequent candidates; if it is in the low noise case, the maximum value of 200 is removed, and {50, 100} is retained. For ordinary noise cases between the high and low noise thresholds, the original candidate set can be kept unchanged.

[0095] After noise level screening, weighted least-squares smoothing is performed on the valid NDVI temporal observations of the target pixel using the remaining λ values, and leave-one-out cross-validation is used to evaluate the reconstruction capability corresponding to different parameters. The weight of each temporal observation is determined based on the validity status and data quality of the corresponding observation. Actual observations that pass the validity screening are assigned non-zero weights, while invalid observations removed by clouds, cloud shadows, or sensor anomaly masks are assigned zero weights. For valid observations with different quality levels, different non-zero weights can be configured according to their quality levels. Specifically, each time, one observation is temporarily removed from the existing valid observations, temporal smoothing is performed using the remaining observations, and the NDVI of the removed position is predicted using the smoothed result. The predicted value is then compared with the original valid observation value. This process is repeated for each valid observation, and the root mean square error of reconstruction under the corresponding λ is calculated.

[0096] By comparing the root mean square error of reconstruction corresponding to each candidate λ, the candidate value with the smallest error is selected as the final temporal smoothing parameter used for the current target pixel. Therefore, the final parameter is not merely fixed in advance based on the project's climate zone, but is further determined by combining the data noise characteristics of the target pixel itself and the actual reconstruction error. When the reconstruction errors corresponding to multiple candidate λ are relatively close, a smaller λ can be preferentially used to preserve more details of vegetation temporal changes while meeting reconstruction accuracy requirements; in the example, when the difference in reconstruction errors among multiple candidate parameters is less than 5%, a smaller λ is preferred.

[0097] Through the above processing, each target pixel can determine an appropriate smoothing parameter according to its climate region and its own time series data quality, thereby reducing the mismatch caused by fixed smoothing parameters for pixels in different regions and with different data quality, improving the ability of the smoothed NDVI time series to maintain the true vegetation change trend, and providing a stable data foundation for subsequent observation synthesis and missing pixel reconstruction.

[0098] In some embodiments of this application, a further implementation of observation synthesis and missing pixel reconstruction in step S103 includes: performing corresponding observation synthesis or missing pixel reconstruction based on the number of effective observations of the target pixel within the target time window, including: The number of valid observations N for the target pixel within the baseline annual growing season window is counted, where N represents the number of valid NDVI observations of the target pixel after invalid observations are removed within the baseline annual growing season window. When N≥5, the median of the effective observations is taken as the baseline annual NDVI of the target pixel; When 3 ≤ N < 5, the arithmetic mean of the effective observations is taken as the baseline annual NDVI of the target pixel; When N=2, the difference between two valid observations is calculated. When the difference is not greater than a preset difference threshold, the average of the two valid observations is taken as the baseline annual NDVI of the target pixel. When the difference is greater than the preset difference threshold, the missing pixel reconstruction is initiated. When N=1, the single valid observation is used as the baseline annual NDVI of the target pixel, and a single observation synthesis status marker is set. When N=0, proceed to the missing pixel reconstruction stage; In the missing pixel reconstruction, effective neighboring pixels that meet the constraints of uncertainty, interpolation state, and NDVI difference are selected from the spatial neighborhood of the target pixel. When the number of effective neighboring pixels is not less than a preset number, the reconstructed value of the target pixel is calculated using an inverse distance weighting method. When the number of effective neighboring pixels is greater than zero but less than the preset number, the effective temporal observations of the years before and after the baseline year of the target pixel and the effective neighboring pixels are used as training data, and the reconstructed value of the target pixel is obtained through Gaussian process regression. When there are no effective neighboring pixels that meet the conditions, temporal smoothing extrapolation is performed based on the effective temporal observations of the years before and after the baseline year of the target pixel. When there are still no effective observations of the target pixel within a preset historical time range, the target pixel is set to an unreconstructable state.

[0099] The baseline year growing season window is used to define the effective observation time range for NDVI synthesis in the baseline year. In one implementation, the baseline monitoring period or vegetation growing season explicitly specified in the methodology application document is preferentially used as the baseline year growing season window; when the methodology application document does not specify the corresponding time range, the corresponding growing season time range is obtained from the preset climate type-growing season window mapping relationship according to the climate zone where the project is located. For the current year, a time window that is the same as or in the same vegetation phenological stage as the baseline year is used to ensure the comparability of remote sensing results between the baseline year and the current year. As mentioned above, in one implementation, the growing season window corresponding to the baseline year is first determined, and the number N of valid NDVI observations remaining within the window after removing clouds, cloud shadows, sensor failures, and other invalid observations is counted. N is used to reflect the data completeness of the current pixel in the main vegetation growth stage of the baseline year and serves as the basis for selecting the synthesis method or reconstruction method.

[0100] When N≥5, it indicates that there are sufficient effective observations during the baseline year's growing season, and the median of all effective NDVI observations can be taken as the baseline year's NDVI for the target pixel. Using the median can reduce the impact of a small number of residual outlier observations on the final result.

[0101] When 3 ≤ N < 5, although the number of effective observations is relatively small, it can still reflect the overall vegetation status during the growing season. Therefore, the arithmetic mean of each effective observation is taken as the baseline annual NDVI of the target pixel. If necessary, corresponding confidence levels can be set for this type of small sample synthesis results to identify the reliability of the data during subsequent statistical analysis or authenticity verification.

[0102] When N=2, the absolute difference between two valid NDVI observations is further calculated. If the difference is not greater than a preset difference threshold, it indicates that there is no significant conflict in vegetation status between the two observations, and the average of the two can be taken as the baseline annual NDVI. If the difference is greater than the preset difference threshold, it may be that the two observations are in different phenological stages or that one of the observations is not representative enough. Therefore, the average is not taken directly, and the missing pixel reconstruction process is initiated. In one specific implementation, the preset difference threshold can be 0.15.

[0103] When N=1, since there is still one valid observation, this observation value can be directly used as the baseline year NDVI. At the same time, a single observation composite status marker is set to indicate that the baseline year value of this pixel is formed by only one valid observation. This preserves the real observation information while enabling subsequent plot statistics processes to identify plots with relatively weak data foundations.

[0104] When N=0, it means that there are no valid NDVI observations that can be used directly within the baseline annual growing season window, so the missing pixel reconstruction process is initiated.

[0105] In the missing pixel reconstruction process, a spatial neighborhood is first established centered on the target pixel, and valid pixels that can be used as the basis for reconstruction are selected from the neighborhood. The selection process comprehensively considers the uncertainty of neighboring pixels, their synthesis or interpolation status, and the NDVI difference between them and the target pixel. The uncertainty of neighboring pixels represents the reliability of the baseline year NDVI result for that pixel, and can be determined based on the dispersion of valid observations used in the baseline year NDVI synthesis or the uncertainty index output by the corresponding reconstruction model. For pixels obtained by direct statistical synthesis from multiple actual valid observations, the uncertainty can be calculated based on the dispersion of the corresponding valid NDVI observations; for pixels formed by the reconstruction model, the prediction uncertainty output by the corresponding model can be used as their uncertainty. The interpolation status is used to identify whether the baseline year NDVI of neighboring pixels was obtained from missing pixel reconstruction. To avoid the continuous propagation of reconstruction errors in the spatial range, non-interpolated pixels obtained directly from actual valid observations or statistically synthesized from valid observations are preferentially selected as valid neighboring pixels; pixels already obtained from missing pixel reconstruction can be marked with an interpolation status and excluded during neighborhood sample selection. Specifically, pixels with high uncertainty, pixels whose values ​​have already been obtained through interpolation, and pixels whose NDVI features differ significantly from the target pixel can be removed to prevent low-quality data or duplicate interpolation results from continuing to participate in new pixel reconstruction. In one implementation, the uncertainty of neighboring pixels is required to be less than 0.08, and the NDVI difference between the neighboring pixels and the center pixel is required to be no greater than 0.20.

[0106] When the number of effective neighboring pixels obtained through screening is not less than a preset number, the reconstructed value of the target pixel is calculated using an inverse distance weighting method. Inverse distance weighting sets weights based on the spatial distance between neighboring pixels and the target pixel; the closer the neighboring pixels are, the greater their contribution to the reconstruction result, thus utilizing the spatial continuity of the local surface state to recover missing pixels. In one specific implementation, the preset number can be three; inverse distance weighting can be used when there are at least three effective neighboring pixels.

[0107] When the number of effective neighboring pixels is greater than zero but less than a preset number, it is difficult to obtain a stable inverse distance weighted result based on only a small number of spatial samples. Therefore, time-dimensional data is further introduced. Specifically, effective time-series observations of the target pixel in the years before and after the baseline year can be obtained and used together with the currently obtained effective neighboring pixels as training data. The predicted NDVI of the target pixel in the baseline year is obtained through Gaussian process regression.

[0108] Specifically, when using Gaussian process regression to reconstruct missing pixels, the temporal observation information of the target pixel and the spatial observation information of its effective neighboring pixels are organized into training data. The temporal training samples include effective NDVI time-series observations of the target pixel in the years immediately before and after the baseline year; each temporal training sample includes at least the observation time and the corresponding NDVI value. The spatial training samples include effective neighboring pixels selected through the aforementioned uncertainty, interpolation state, and NDVI difference constraints; each spatial training sample includes at least the spatial position of the neighboring pixel relative to the target pixel and the corresponding baseline year's NDVI. In one implementation, effective temporal data within a one-year range before and after the baseline year are preferentially used. When the effective time samples within this range are insufficient, the time sample range can be further expanded.

[0109] Gaussian process regression uses a covariance kernel function to represent the correlation between training samples. When prediction mainly relies on temporal observations of target pixels, the Matérn 3 / 2 kernel function can be used to describe the temporal correlation between NDVI observations. When both effective temporal samples and spatial neighborhood samples are available, a spatiotemporal additive kernel function containing temporal and spatial correlation terms can be used to allow the model to consider both the temporal continuity and local spatial continuity of NDVI. For items with significant periodic phenological changes, a periodic correlation term can also be added to the kernel function.

[0110] The hyperparameters corresponding to the kernel function include at least length scale parameters and observation noise parameters. During model training, the initial values ​​of the corresponding parameters can be set according to the time interval, spatial distance, and NDVI dispersion of the training samples, and the kernel function parameters corresponding to the current training data can be determined through a preset parameter optimization process; for example, the parameters can be iteratively adjusted using the Gaussian process marginal likelihood of the training data as the optimization objective. After the model training is completed, the temporal and spatial features corresponding to the target pixel at the baseline year target time position are input into the model to obtain the corresponding NDVI prediction mean and prediction standard deviation σ_gpr. The prediction mean is used as the reconstructed value of the target pixel, and σ_gpr is used to characterize the uncertainty of the reconstructed value. For example, when σ_gpr < 0.05, it can be set to a high confidence state; when 0.05 ≤ σ_gpr < 0.10, it can be set to a medium confidence state; when 0.10 ≤ σ_gpr < 0.15, it can be set to a low confidence state; when σ_gpr ≥ 0.15, it can be set to an extremely low confidence state, and the corresponding reconstructed value can be excluded from the subsequent plot representative value statistics.

[0111] When there are no valid neighboring pixels that meet the screening requirements around the target pixel, spatial information is no longer used. Instead, temporal smoothing extrapolation is performed based on the valid temporal observations of the pixel in the years before and after the baseline year to predict the NDVI corresponding to the target time window of the baseline year based on the changing trend of the historical time series.

[0112] If no valid observations are found for the target pixel within a preset historical timeframe, the existing remote sensing data is insufficient to generate a reliable reconstruction result. In this case, the target pixel is set to an unreconstructable state, rather than being artificially assigned an overly speculative NDVI value. Unreconstructable pixels do not participate in the subsequent calculation of the representative NDVI value for the plot, but their quantity and area information can be retained to evaluate the effective coverage of the baseline year's remote sensing reconstruction results.

[0113] Through the hierarchical processing described above, which involves "branching and synthesizing based on the number of effective observations, spatial neighborhood filtering, inverse distance weighting, spatiotemporal joint regression, temporal extrapolation, and non-reconstructable markers," we can prioritize the use of real and effective observations and only invoke reconstruction methods step by step when actual data is insufficient. This approach balances the data authenticity and spatial integrity of the baseline year NDVI results, providing a stable pixel-level data foundation for the subsequent formation of baseline year NDVI raster images and plot baseline year NDVI.

[0114] In some embodiments of this application, a further implementation of the NDVI-land cover mapping calibration and land surface type inference in step S104 includes: calibrating a preset NDVI-land cover mapping interval based on the climate zone where the project is located and determining the inferred land surface type, including: The climate zone where the project is located is determined based on the project vector boundary and the preset climate zone data, and the lower boundary correction coefficient and upper boundary correction coefficient of the corresponding NDVI interval are obtained. The lower and upper boundaries of the preset NDVI-land category mapping interval are calculated with the corresponding correction coefficients to form the regionally calibrated NDVI-land category mapping interval. The NDVI of each pixel in the baseline annual NDVI raster image is input into the NDVI-land class mapping interval after regional calibration for pixel-level classification. When the NDVI of the target pixel is within the boundary blur zone of the adjacent mapping interval, neighborhood determination and mixed pixel decomposition are performed on the target pixel to obtain at least one land surface type and its area ratio corresponding to the target pixel. Based on the classification results of each pixel within the project vector boundary, the area proportion of each surface type is calculated, and the land parcel-level inferred surface type is determined based on the area proportion.

[0115] As described above, in one implementation, the climate zone where the project is located is first determined based on the project vector boundary and preset climate zoning data. Specifically, the project vector boundary can be spatially matched with the preset climate zoning boundary or climate zoning raster to obtain the climate type corresponding to the project. When the entire project is located in a single climate zone, the calibration parameters corresponding to that climate zone are directly used; when the project vector boundary crosses multiple climate zones, the project area is divided into corresponding climate sub-regions based on the intersection results of the project vector boundary and the spatial boundaries of each climate zone; the pixel-level land surface type classification within the corresponding sub-region is completed using the lower boundary correction coefficient and upper boundary correction coefficient of the NDVI interval corresponding to each climate sub-region, and then the classification results of each climate sub-region are summarized to obtain the area proportion of each land surface type within the entire project area. In one implementation, it is also stated that when the methodology has already clearly defined the project climate type, the corresponding information can be directly used; when only point coordinates are available, location matching can be performed through the climate zoning raster.

[0116] After determining the climate zone, the lower and upper boundary correction coefficients of the corresponding NDVI intervals for that climate zone are obtained from a preset parameter set and applied to the lower and upper boundaries of each preset NDVI-land type mapping interval to form a regionally calibrated mapping interval. This process is used to reflect the differences in the overall NDVI level of the same type of vegetation under different climatic conditions. In one implementation, the preset parameter set pre-stores the basic NDVI intervals corresponding to each standard land surface type, for example, setting the lower and upper boundaries of basic NDVI for water bodies, bare land, grassland, shrubs, sparse forest, and woodland respectively. Each basic interval is the basic mapping relationship used before regional climate correction. After the climate zone to which the project belongs is determined, the basic intervals are then regionally calibrated using the corresponding lower and upper boundary correction coefficients. In the example, the upper and lower boundary correction coefficients are both 1.00 for the humid zone, 0.92 and 0.95 for the semi-humid zone, 0.80 and 0.85 for the semi-arid zone, 0.65 and 0.72 for the arid zone, and 0.75 and 0.80 for the plateau climate zone. In one implementation, the lower boundary of the preset NDVI-land type mapping interval is multiplied by the lower boundary correction factor, and the upper boundary is multiplied by the upper boundary correction factor to obtain the lower and upper boundaries of the regionally calibrated NDVI interval, thus forming the regionally calibrated NDVI-land type mapping interval. Specifically, if the preset NDVI interval corresponding to a certain land surface type is [L, U], and the lower boundary correction factor corresponding to the climate zone where the project is located is α_L and the upper boundary correction factor is α_U, then the calibrated NDVI interval can be represented as [α_L×L, α_U×U].

[0117] After completing the regional calibration, each valid pixel in the baseline year NDVI raster image located within the project vector boundary is matched one by one with the calibrated NDVI-land type mapping interval. When a pixel's NDVI is within a certain mapping interval and far from the boundary of adjacent intervals, the pixel can be directly identified as the corresponding land type. In this way, continuous NDVI values ​​are converted into pixel-level classification results with actual land cover meanings, such as water bodies, bare soil, grassland, shrubs, sparse forest, and woodland.

[0118] For the classification boundary between two adjacent NDVI-land type mapping intervals, a boundary blurring band can be set according to the corresponding interval width. In one embodiment, based on the preset width corresponding to the adjacent mapping intervals, an NDVI range equivalent to 5% of the interval width is set on both sides of the classification boundary as a boundary blurring band. When the NDVI of a target pixel falls into this range, the land type is not directly determined according to a single threshold, but neighborhood determination and hybrid pixel decomposition are further performed.

[0119] When a target pixel may contain two land surface types simultaneously, a hybrid pixel decomposition can be performed. For binary hybrid cases, the composition ratio of the two land surface types can be estimated based on the relative position of the pixel's NDVI between the feature values ​​of the two candidate land surface types. For example, if the first candidate land surface type represents NDVI_1, the second candidate land surface type represents NDVI_2, the target pixel is NDVI_p, and the area ratio of the first candidate land surface type is q, then: NDVI_p = q × NDVI_1 + (1-q) × NDVI_2. For three or more candidate types, the two types closest to the current NDVI can be selected for decomposition. In the example, when the estimated proportion of a certain land surface type is less than 10%, it can be merged into the dominant land surface type; when the proportion is not less than 10%, the corresponding land surface type is retained, and its respective area ratio is recorded.

[0120] After classifying all valid pixels within the project area, the number of pixels corresponding to different land surface types is counted based on the project vector boundaries. This number is then converted into the area of ​​each land surface type and its proportion of the project's effective area based on the raster spatial resolution. Land surface types are then inferred at the plot level based on this area proportion.

[0121] In one specific approach, when a single land surface type accounts for 80% of the area, the project site can be identified as that single dominant land type; when the combined area of ​​two main land types accounts for 80%, a mixed land parcel type consisting of the two land types can be output; when there are three or more main land types and no single type accounts for less than 50%, the types with higher area proportions and their corresponding proportions can be used as the land parcel-level inference results. When the area proportion of non-vegetation types is high, corresponding abnormal status information can be output.

[0122] Through the above-mentioned processing of "climate zone matching, mapping interval calibration, pixel-level classification, boundary blurring, hybrid pixel decomposition, and plot area summarization", the baseline annual NDVI raster image can be transformed from continuous numerical data into spatial classification results with land surface type meaning. Furthermore, it can form plot-level inferred land surface types and the area proportion of each type, providing a data foundation for subsequent selection of carbon sink parameter conversion models, determination of regional parameters, and verification of the consistency of land types in the declared baseline.

[0123] In some embodiments of this application, a further implementation of the determination of the carbon sink parameter conversion model and regional parameters in step S104 includes: determining the corresponding carbon sink parameter conversion model and regional parameters based on the project type information obtained from the methodology application document and / or the inferred land surface type, including: The project type is determined based on the project type keywords in the methodological application document; when no project type is identified, the project type is determined based on the area proportion of the inferred land surface type. When the project type corresponds to forestry carbon sinks, a spatial intersection operation is performed between the project vector boundary and the preset ecological forest area vector boundary. The corresponding ecological forest area is determined based on the intersection result, and the first biomass model parameter a and the second biomass model parameter b corresponding to the ecological forest area are obtained from the preset ecological forest area parameter set, according to: Estimated aboveground biomass = exp(a + b × plot baseline year NDVI) Calculate the extrapolated aboveground biomass, where a and b are regional model parameters determined based on the ecological forest area and used to characterize the mapping relationship between the baseline year NDVI and aboveground biomass of the plot; The estimated aboveground biomass is converted based on the carbon content corresponding to the inferred land surface type to obtain the estimated carbon density. When the project type corresponds to a carbon sink type using other biomass conversion models, the conversion model and regional parameters corresponding to the project type are retrieved from the preset project type parameter set, and the corresponding extrapolated baseline accounting parameters are calculated based on the baseline year NDVI of the plot. When the baseline calculation parameter to be verified includes soil carbon density, the estimated soil carbon density is determined from the preset soil carbon density parameter set based on the project surface type, climate zone, and available soil type information, and the corresponding confidence level is set according to the completeness of the input data.

[0124] As described above, in one implementation, keywords characterizing carbon sink project types are first identified from the methodological application documents, and the project type corresponding to the current project is determined accordingly. For example, when information such as "afforestation" or "forest management" appears in the application materials, the project can be identified as a forestry carbon sink type; when information such as "grassland" or "terrace" appears, it can be identified as a grassland carbon sink type; when information such as "wetland" or "mangrove" appears, it can be identified as a wetland carbon sink type; and when information such as "farmland" or "conservation tillage" appears, it can be identified as a farmland carbon sink type.

[0125] When the project type cannot be determined based on keywords in the methodology application document, the inferred land cover types and their area proportions obtained from the aforementioned steps can be used for supplementary judgment. For example, when forest types such as arbor forests and plantations occupy the majority of the project area, the project can be matched to the forestry carbon sink type; when grassland types occupy the majority of the area, it can be matched to the grassland carbon sink type. Thus, the determination of the project type not only relies on the textual descriptions in the application document but also utilizes land cover composition inferred from remote sensing as a supplementary basis.

[0126] When the project type is determined to be forestry carbon sink, the corresponding ecological forest area is further determined based on the project vector boundary. Specifically, spatial intersection operation is performed between the project vector boundary and the pre-stored ecological forest area vector boundary, and the ecological forest area to which the project belongs is determined based on the intersection area. When the entire project is located in a single ecological forest area, the regional parameters corresponding to that ecological forest area are directly called; when the project vector boundary crosses multiple ecological forest areas, the project is divided into multiple ecological forest area sub-regions based on the spatial intersection results between the project vector boundary and each ecological forest area, and the biomass model parameters corresponding to each sub-region are obtained. For the i-th ecological forest area sub-region, the estimated aboveground biomass B_i of the sub-region is calculated based on the plot or pixel statistical NDVI of the sub-region and the corresponding first biomass model parameter a_i and second biomass model parameter b_i, B_i=exp(a_i+b_i×NDVI_i), where NDVI_i represents the statistical NDVI of the i-th ecological forest area sub-region; then, B_i is weighted and summarized according to the proportion of the area of ​​each ecological forest area sub-region to the effective area of ​​the project to obtain the estimated aboveground biomass of the project area. Therefore, the regional model parameters of each ecological forest area are applied to the NDVI of the corresponding area, instead of simply weighting the nonlinear model parameters a and b themselves first.

[0127] For projects located entirely within a single ecological forest area, after identifying the ecological forest area, the first biomass model parameter 'a' and the second biomass model parameter 'b' are obtained from a pre-set ecological forest area parameter set. The baseline year NDVI of the plot is then input into the corresponding biomass conversion relationship to obtain the estimated aboveground biomass. For the forestry index model given in the example, the calculation can be performed as: Estimated aboveground biomass = exp(a + b × plot baseline year NDVI). Here, 'a' and 'b' are regional model parameters used to characterize the correspondence between the plot baseline year NDVI and aboveground biomass under specific ecological forest area conditions. Different 'a' and 'b' parameters are used for different ecological forest areas to adapt to differences in vegetation communities, growth conditions, and biomass levels in different regions.

[0128] After obtaining the estimated aboveground biomass, the carbon content is converted based on the carbon content corresponding to the estimated land surface type to obtain the estimated carbon density. For example, in one parameter configuration, the carbon content of a arbor forest can be taken as 0.50, a bamboo forest as 0.48, and a shrub forest as 0.45. Therefore, the estimated aboveground biomass can be multiplied by the corresponding carbon content to obtain the estimated carbon density per unit area.

[0129] For carbon sink projects that do not use the aforementioned forestry index model, the corresponding conversion model and regional parameters are retrieved from the preset project type parameter set according to the project type. One example clearly distinguishes the model formats for different project types: forestry arbor forests use an index model, bamboo forests use a corresponding linear model, grassland carbon sinks use a linear model corresponding to grassland regional parameters, wetlands or mangroves use the corresponding index model, and farmland carbon sinks do not directly utilize NDVI to invert biomass, but primarily estimate soil carbon parameters. For other carbon sink project types, the preset project type parameter set must at least store the project type identifier, the corresponding conversion model, the regional parameters required by the model, and the type of output calculation parameters. Once the project type is determined, the corresponding model and regional parameters are indexed according to the project type, and the baseline year NDVI of the plot or other necessary remote sensing parameters corresponding to this project type are input into the model to obtain the extrapolated parameters of the same type as the declared baseline calculation parameters.

[0130] Therefore, during actual equipment operation, the corresponding model can be called based on the determined project type, instead of directly applying the regional parameters of forestry projects to grasslands, wetlands, or farmland. For example, when the project type is grassland carbon sink, the conversion model and regional parameters of the corresponding grassland area are obtained from the preset grassland project parameter set, and then the baseline year NDVI of the plot is input into the model to obtain the corresponding extrapolated baseline calculation parameters; when the project type is wetland carbon sink, the parameter model that matches wetlands or mangroves is called.

[0131] When a project encompasses multiple major land surface types, the corresponding conversion relationships can be applied based on the spatial area of ​​each land surface type to obtain the estimated parameters for different areas. These parameters are then summarized based on the area proportion of each land surface type. This approach can reduce the errors caused by using a single model parameter for mixed-use land parcels.

[0132] When the baseline calculation parameters to be verified also include soil carbon density, instead of directly establishing a single soil carbon density conversion function based on NDVI, the corresponding reference value is determined from the preset soil carbon density parameter set based on the project's surface type, climate zone, and available soil type information.

[0133] In one specific approach, when measured soil data independent of the parameters to be verified in the declaration can be obtained, this independent measured data can be used as high-confidence reference data, and the estimated soil carbon density or the preset soil carbon density parameter can be determined accordingly. When the obtained measured soil data is the declared soil carbon density to be verified in the methodology declaration document, this data is only used as the declaration baseline calculation parameter and not as the basis for determining the estimated soil carbon density. When measured data is lacking but soil type can be determined, parameter matching is performed based on three levels of conditions: surface type, climate zone, and soil type, and the data is set to a medium confidence level. When only surface type and climate zone can be determined, the corresponding regional average soil carbon density reference value can be used and set to a low confidence level. When the existing input information is insufficient to form a valid estimate, the soil carbon density is set to an incalculable state and does not participate in the quantitative verification of the corresponding parameter dimension.

[0134] Therefore, through the processing of "project type identification - ecological region matching - model and regional parameter calling - biomass conversion - carbon content conversion - model adaptation to other project types - soil carbon density classification estimation", the annual NDVI of the plot baseline can be further converted into the calculation parameters of the inferred baseline that match the project type and regional conditions, and the corresponding model, regional parameters and confidence status are recorded, providing unified remote sensing inferred data with clear sources for subsequent calculation of parameter consistency deviation.

[0135] In some embodiments of this application, a further implementation of determining the authenticity verification level in step S105 includes: determining the baseline authenticity verification level includes: The land use consistency deviation δ_land is calculated using the following formula: δ_land=|NDVI_parcel-NDVI_mid| / W_land; Wherein, NDVI_parcel is the annual NDVI of the land parcel baseline, NDVI_mid is the midpoint of the NDVI interval after regional calibration corresponding to the declared baseline surface type, and W_land is the width of the NDVI interval after regional calibration. The units and calculation methods of the estimated baseline calculation parameters and the declared baseline calculation parameters are standardized, and the parameter deviation δ_param is calculated according to the following formula when the declared baseline calculation parameter is greater than zero: δ_param=|P_est-P_decl| / P_decl; Where P_est is the calculated baseline accounting parameter, and P_decl is the declared baseline accounting parameter corresponding to the calculated baseline accounting parameter; when P_decl is zero, the absolute deviation is compared with the corresponding preset absolute threshold. Based on the baseline year NDVI and the current year NDVI of the land parcel, the corresponding carbon density or carbon sink change parameters are calculated and converted into the same change indicators as the declared carbon sink growth parameters, resulting in the remote sensing change rate r_rs and the declared change rate r_decl. The change consistency deviation δ_change is then calculated according to the following formula: δ_change=|r_rs-r_decl| / |r_decl|; Wherein, r_rs represents the remote sensing carbon sink change rate determined by the conversion result of carbon sink parameters corresponding to the baseline year NDVI and the current year NDVI of the land parcel, and r_decl represents the declared carbon sink change rate determined by the declared carbon sink growth parameters in the methodology declaration document; when r_decl is zero, |r_rs-r_decl| is used as the absolute change deviation, and the absolute change deviation is compared with the corresponding preset absolute change threshold; The individual level of each verification dimension is determined based on the preset grading thresholds corresponding to land type consistency deviation, parameter deviation, and change consistency deviation; when the direction of remote sensing change is opposite to the direction of carbon sink growth reported, the change consistency dimension is determined as the highest anomaly level. The comprehensive verification level is determined according to the principle of prioritizing the highest anomaly level; when at least two verification dimensions reach the medium level, the comprehensive verification level is upgraded to the high level; when there are fewer than two valid verification dimensions, the comprehensive verification status is set to undetermined. Based on the comprehensive verification level, output at least one of the following: an abnormal status flag, a data correction request signal, a route review instruction, a verification termination signal, or a verification pass signal.

[0136] As described above, in one implementation, a land use consistency check is first performed between the declared baseline surface type and the remote sensing reconstruction results. Based on the NDVI-land use mapping relationship determined in step S104 after regional calibration, the NDVI interval corresponding to the declared baseline surface type is obtained, and the midpoint NDVI_mid and interval width W_land of this interval are calculated. The annual NDVI of the parcel baseline obtained in step S103 is denoted as NDVI_parcel, and the land use consistency deviation is calculated according to the following formula: δ_land = |NDVI_parcel - NDVI_mid| / W_land.

[0137] Here, NDVI_mid represents the representative location of the NDVI interval corresponding to the declared land surface type, and W_land is used to normalize the absolute difference between the NDVI of the plot and the midpoint of the interval. By using the interval width as the denominator, the impact of the difference in NDVI interval width between different land types on the deviation evaluation can be reduced, making the deviations corresponding to different land types comparable.

[0138] For example, if the NDVI range for the area corresponding to the declared baseline surface type is 0.20 to 0.50 after calibration, then NDVI_mid is 0.35 and W_land is 0.30; if the annual NDVI of the plot baseline is 0.41, then the land type consistency deviation is |0.41-0.35| / 0.30=0.20, or 20%.

[0139] When the methodology application document contains only one primary baseline surface type, the land use consistency deviation corresponding to that surface type is directly used as the land use verification result. When the application project contains multiple baseline surface types, the land use consistency deviation can be calculated separately for each applied surface type and its corresponding project area, and the final level of the land use consistency dimension can be determined based on the individual level corresponding to each land use type. In one implementation, the highest anomaly level can be used as the final level of this dimension to avoid significant local land use conflicts being offset by other normal land use results.

[0140] Subsequently, the consistency between the estimated baseline accounting parameters and the declared baseline accounting parameters was verified. Before calculating parameter deviations, the units of execution and accounting caliber of the two comparison parameters were standardized to ensure that they have the same physical meaning and measurement basis. For example, the unit area parameter can be standardized as tC / ha; when the original declared data uses units such as square kilometers, mu, kilograms, and carbon dioxide equivalents, they are converted to a unified unit according to preset conversion relationships; for different accounting objects such as aboveground biomass, total biomass, and carbon density, they are compared according to their corresponding calibers, and parameters with different physical meanings are not directly mixed for calculation.

[0141] When the baseline calculation parameter P_decl is greater than zero, the parameter deviation is calculated using the following formula: δ_param=|P_est-P_decl| / P_decl. Where P_est represents the estimated baseline calculation parameter obtained from the baseline year's NDVI and the corresponding carbon sink conversion model, and P_decl represents the declared value in the application materials that belongs to the same parameter category as P_est. This deviation reflects the relative difference between the remote sensing calculation result and the declared value.

[0142] When P_decl is zero, the relative deviation formula cannot be used. Therefore, the absolute difference between P_est and P_decl is calculated instead, and this absolute difference is compared with the preset absolute threshold of the corresponding parameter. For cases where the declared parameters are missing, the units are unclear, or a consistent accounting standard cannot be formed, the parameter dimension can be set to have no valid quantitative results and is not required to participate in the comprehensive scoring.

[0143] When the baseline calculation parameters to be verified include multiple different physical quantities, the parameter deviation is calculated for each calculation parameter that can form a consistent estimated value and declared value, and the individual level of each parameter is determined. Instead of directly averaging the numerical deviations between different physical quantities, the final level of the parameter consistency dimension is determined based on the individual level of each parameter. In one implementation, the highest anomaly level can be used as the final level of the parameter consistency dimension.

[0144] Further verification of the consistency between baseline and current status was performed. First, the current year's NDVI was obtained using the same multi-source image normalization and temporal reconstruction rules as the baseline year, ensuring that the remote sensing results for both periods maintained the same numerical benchmark and processing caliber. Subsequently, based on the carbon sink parameter conversion model corresponding to the project type, the baseline year NDVI and the current year NDVI of the land parcel were converted into corresponding carbon density or other carbon sink change parameters, and a change index was formed at the remote sensing end based on the time span from the baseline to the current year.

[0145] Simultaneously, the carbon sink growth parameters declared in the methodology application documents are standardized and converted into the same change indicators as those in the remote sensing data. For example, they can be uniformly converted into the annual average carbon density growth rate per unit area, thereby obtaining the remote sensing change rate r_rs and the declared change rate r_decl, respectively.

[0146] When comparing relative change rates, the consistency deviation is calculated using the following formula: δ_change = |r_rs - r_decl| / |r_decl|. Where r_rs represents the carbon sink change rate calculated from remote sensing results of the baseline year and the current year, and r_decl represents the carbon sink growth change rate corresponding to the reported data. This deviation is used to quantify the difference between the actual change magnitude reflected by remote sensing observations and the reported growth magnitude.

[0147] For special projects where relative change rates cannot be used, other comparison methods can be selected based on the available data. For example, absolute growth can be used for comparisons in cases such as afforestation on barren land where the baseline value is zero; for projects lacking a reliable biomass conversion model, the NDVI change rate can be used for auxiliary judgment; when the data is insufficient to support quantitative comparisons, the method can be reduced to judging the consistency of the direction of change.

[0148] After obtaining the deviations corresponding to the three verification dimensions, the individual level is determined according to the preset grading thresholds corresponding to each dimension. In one implementation, a land type consistency deviation δ_land ≤ 30% is determined to be a level with no obvious contradiction, 30% < δ_land ≤ 50% is determined to be a medium level, and δ_land > 50% is determined to be a high level; a parameter deviation δ_param ≤ 15% is determined to be a level with no obvious contradiction, 15% < δ_param ≤ 30% is determined to be a medium level, and δ_param > 30% is determined to be a high level. For change consistency deviation δ_change, δ_change ≤ 30% is determined to be a level with no obvious contradiction, 30% < δ_change ≤ 50% is determined to be a medium level, 50% < δ_change ≤ 80% is determined to be a high level, and δ_change > 80% is determined to be the highest anomaly level.

[0149] In addition to numerical deviation, it is also necessary to determine whether the direction of change in remote sensing is consistent with the direction of carbon sink growth reported. When remote sensing results indicate a decrease in the project's carbon density or vegetation status relative to the baseline year, while the reporting materials show a continuous increase in carbon sink during the same period, the two directly conflict in terms of the direction of change. In this case, the consistency dimension of change can be directly determined as the highest level of anomaly, instead of judging solely based on the magnitude of relative deviation.

[0150] After obtaining the individual anomaly levels, a comprehensive verification level is formed according to the principle of prioritizing the highest anomaly level. That is, the level with the highest degree of anomaly is selected from all valid verification dimensions as the initial result of the comprehensive verification level. For example, if the land category dimension and the change dimension are in a state of no obvious contradiction, but the parameter dimension reaches a high level, the comprehensive result is still determined to be a high level, so as to avoid serious individual anomalies being offset by the average of other normal dimensions.

[0151] Based on this, further multi-dimensional anomaly adjustments are performed. When at least two valid verification dimensions reach the medium level, it indicates that multiple independent verification results show significant deviations simultaneously, thus the overall verification level is raised to the high level.

[0152] Among them, a valid verification dimension refers to a verification dimension that possesses the required declaration information and remote sensing inference results, and can complete the corresponding deviation calculation or change direction consistency judgment according to a unified unit and accounting standard. For situations where declaration parameters are missing, units cannot be unified, remote sensing data cannot form valid inference results, or other situations where corresponding comparisons cannot be completed, the corresponding verification dimension will be set as an invalid verification dimension. When there are fewer than two verification dimensions that can form valid results, the existing information is insufficient to support a stable comprehensive authenticity judgment, the comprehensive verification status will be set to "unable to determine," and data supplementation processing will begin, rather than directly forming a final authenticity conclusion based on a single dimension.

[0153] Finally, corresponding device status flags and control signals are generated based on the comprehensive verification level. For the highest anomaly level, a highest anomaly status flag and a control signal to terminate the current verification process can be output; for a high level, an anomaly status flag and a data correction request signal can be output; for a medium level, a review routing instruction can be output, transferring the current task to the corresponding review process; for states without obvious contradictions, a verification pass signal can be output; for states that cannot be determined, a data correction request signal is output, along with the parameters or missing dimensions that need to be supplemented. By using status signals and control instructions to represent the verification results, administrative actions such as "rejection" and "manual review" can be avoided as direct method steps, ensuring that the process remains a data processing flow that computer equipment can automatically execute.

[0154] Through the above-mentioned processing steps of "quantification of three types of deviations - determination of single-item level - identification of directional anomalies - priority of the highest level - escalation of multi-dimensional anomalies - judgment of validity - output of status and control signals", the remote sensing reconstruction and carbon sink parameter estimation results obtained in the previous steps can be converted into baseline authenticity verification results that are quantifiable, traceable and can be automatically processed by the equipment.

[0155] Compared with existing technologies, this application discloses a method for verifying the authenticity of carbon sink baselines based on multi-source remote sensing time-series reconstruction. It establishes a cross-sensor normalization relationship through pseudo-invariant features to reduce the systematic differences in NDVI from multiple satellites. It adaptively determines the time-series smoothing and reconstruction methods based on climate zones and the number of effective observations, improving the stability of annual NDVI reconstruction of the baseline. It forms comparable baseline accounting parameters through regionally calibrated NDVI-land type mapping and a carbon sink parameter conversion model matched to the project type. Furthermore, it establishes an authenticity verification level based on three dimensions: land type, accounting parameters, and changes from baseline to current status, achieving closed-loop processing from multi-source remote sensing reconstruction to verification signal output.

[0156] Based on the same inventive concept as the methods described above, this application also proposes a carbon sink baseline authenticity verification system based on multi-source remote sensing time-series reconstruction, such as... Figure 2The diagram shown is a structural schematic of a carbon sink baseline authenticity verification system based on multi-source remote sensing time-series reconstruction. The system includes: The basic parameter extraction module is used to receive methodology application documents, extract baseline year, baseline surface type, project vector boundary and baseline calculation parameters from the methodology application documents, and generate a structured basic parameter object; The multi-source image normalization module is used to acquire multi-source satellite image metadata within the target time window based on the baseline year and project vector boundary. It compares each screening index with corresponding preset validity conditions based on image cloud coverage, effective coverage ratio within the project vector boundary, the matching degree between image acquisition time and the target time window, and sensor data availability. Multi-source satellite images that meet the preset validity conditions are retained. Pseudo-invariant feature pixels are extracted from the retained images. A linear transformation relationship between different sensors is established based on these pseudo-invariant feature pixels. Pixel observations from non-reference sensors are normalized to the same numerical reference, generating normalized time-series data. Furthermore, when performing change consistency verification, the current year's multi-source satellite images are processed according to the same validity screening and normalization rules to generate current year's normalized time-series data. The temporal reconstruction module is used to perform abnormal pixel identification and mask removal on the normalized temporal data, determine the temporal smoothing parameters according to the climate zone where the project is located, and perform corresponding observation value synthesis or missing pixel reconstruction based on the effective number of observations of the target pixel within the target time window to generate a baseline year NDVI raster image. It also performs plot-scale statistics on the baseline year NDVI raster image according to the project vector boundary to generate the plot baseline year NDVI; and processes the current year normalized temporal data according to the same temporal reconstruction and plot-scale statistics rules as the baseline year to generate the plot current year NDVI. The parameter conversion module is used to calibrate the preset NDVI-land type mapping interval according to the climate zone where the project is located, determine the inferred land type based on the calibrated NDVI-land type mapping interval and the baseline year NDVI raster image, and determine the corresponding carbon sink parameter conversion model and regional parameters based on the project type information and / or inferred land type obtained from the methodology application document, and calculate the baseline accounting parameters based on the baseline year NDVI of the plot. The authenticity verification module is used to determine the land type consistency deviation based on the annual NDVI of the land parcel baseline and the corresponding NDVI interval of the land surface type of the declared baseline, to determine the parameter deviation based on the calculated parameters of the estimated baseline and the calculated parameters of the declared baseline, and to determine the change consistency deviation based on the annual NDVI of the land parcel baseline, the annual NDVI of the current status of the land parcel, and the declared carbon sink growth parameters; to determine the baseline authenticity verification level according to the preset grading rules corresponding to each deviation, and to output the status flag signal and verification control command corresponding to the verification level.

[0157] Those skilled in the art will understand that embodiments of this application can be provided as methods, systems, or computer program products. Therefore, this application can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, this application can take the form of a computer program product embodied on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0158] This application is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of this application. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0159] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0160] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0161] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit it. Although the present invention has been described in detail with reference to the above embodiments, those skilled in the art should understand that modifications or equivalent substitutions can still be made to the specific implementation of the present invention. Any modifications or equivalent substitutions that do not depart from the spirit and scope of the present invention should be covered within the scope of protection of the claims of the present invention.

Claims

1. A method for verifying the authenticity of carbon sink baselines based on multi-source remote sensing time-series reconstruction, characterized in that, include: Receive the methodology application document, extract the baseline year, baseline surface type, project vector boundary and baseline calculation parameters from the methodology application document, and generate a structured basic parameter object; Based on the baseline year and project vector boundary, multi-source satellite image metadata within the target time window is obtained. According to the image cloud coverage, the effective coverage ratio within the project vector boundary, the matching degree between the image acquisition time and the target time window, and the availability of sensor data, each screening index is compared with the corresponding preset validity conditions, and multi-source satellite images that meet the preset validity conditions are retained. Pseudo-invariant feature pixels are extracted from the retained images. Based on the pseudo-invariant feature pixels, a linear transformation relationship between different sensors is established, and the pixel observation values ​​of non-reference sensors are normalized to the same numerical reference to generate normalized time series data. Abnormal pixel identification and mask removal are performed on the normalized time series data. The time series smoothing parameters are determined from the preset parameter set according to the climate zone where the project is located. Based on the number of effective observations of the target pixel within the target time window, the corresponding observation value synthesis or missing pixel reconstruction is performed to generate the baseline year NDVI raster image. The plot-scale statistics of the baseline year NDVI raster image are performed on the project vector boundary to generate the plot baseline year NDVI. The preset NDVI-land type mapping interval is calibrated according to the climate zone where the project is located. The inferred land type is determined based on the calibrated NDVI-land type mapping interval and the baseline year NDVI raster image. The corresponding carbon sink parameter conversion model and regional parameters are determined based on the project type information obtained from the methodology application document and / or the inferred land type. The baseline accounting parameters are calculated based on the baseline year NDVI of the land parcel. The land use consistency deviation is calculated based on the NDVI of the baseline year of the land parcel and the NDVI interval corresponding to the surface type of the declared baseline. The parameter deviation is calculated based on the calculation parameters of the estimated baseline and the calculation parameters of the declared baseline. The carbon sink growth parameters are obtained from the methodological application document. Multi-source satellite images of the current year are obtained. The current year NDVI of the land parcel is obtained according to the same normalization, reconstruction and land parcel scale statistical rules as the baseline year. The change consistency deviation is determined based on the NDVI of the baseline year of the land parcel, the current year NDVI of the land parcel and the declared carbon sink growth parameters. The baseline authenticity verification level is determined according to the preset grading rules corresponding to each deviation, and the status mark signal and verification control command corresponding to the verification level are output.

2. The method according to claim 1, characterized in that, Extract baseline year, baseline surface type, project vector boundary, and baseline calculation parameters from the methodology application documents, including: The method application document is checked to see if it contains a structured attachment that matches the preset data structure. If it does, the corresponding fields are extracted from the structured attachment according to the preset field mapping relationship. If the structured attachment does not exist or the structured attachment fails to be parsed, the layout of the methodological application document is analyzed and paragraphs are segmented. The corresponding fields are located by using preset regular expressions, keyword dictionaries and named entity recognition models, and the extraction results and extraction confidence of each field are output. For fields directly extracted from the structured attachments, they are used as fields to be verified; for fields automatically extracted from the main text of the methodology application document, they are used as fields to be verified, and their extraction confidence is compared with a preset extraction confidence threshold. When the extraction confidence of any field is lower than the preset extraction confidence threshold, a corresponding data correction request signal is generated and the returned data is received. When each field to be verified meets any of the following conditions: direct extraction of structured attachments, extraction confidence not lower than the preset extraction confidence threshold, or receiving feedback data that meets the preset field format for low confidence fields, the year range, numerical unit, geometric boundary, and parameter validity are verified for each current field. After the verification is passed, the basic parameter object is generated.

3. The method according to claim 1, characterized in that, The extraction of pseudo-invariant feature pixels includes: Candidate pixels are determined within the extended buffer of the project vector boundary, and time-series observation data of each candidate pixel are obtained within a preset cross-year time range including the baseline year. For each candidate pixel, calculate the temporal standard deviation σ_NIR of the near-infrared band reflectance, the difference between the temporal maximum and minimum NDVI values ​​Δ_NDVI, and the 95th percentile of the NDVI difference between adjacent time steps. Here, σ_NIR represents the degree of dispersion of the near-infrared band reflectance of the candidate pixel over time, and Δ_NDVI represents the magnitude of the NDVI change of the candidate pixel within the preset year-long time range. When σ_NIR < 0.015, Δ_NDVI < 0.10, and the 95th percentile < 0.05, the corresponding candidate pixel is determined as a pseudo-invariant feature pixel. When the number of pseudo-invariant feature pixels obtained is lower than the preset number of samples, the filtering thresholds corresponding to σ_NIR and Δ_NDVI are increased and the filtering is re-executed; when the number of pseudo-invariant feature pixels after re-filtering is still lower than the preset number of samples, a pseudo-invariant feature sample shortage status signal is output.

4. The method according to claim 3, characterized in that, The linear conversion relationship between different sensors and supplementary auxiliary data include: Linear regression is performed based on paired observations of the same pseudo-invariant feature pixel in the reference sensor and the sensor to be converted to determine the first conversion coefficient k and the second conversion coefficient c, and a linear conversion relationship is established: NDVI baseline = k × NDVI to be converted + c; Wherein, NDVI reference represents the NDVI value after normalization to the reference sensor value reference, NDVI to be converted represents the NDVI observation value corresponding to the sensor to be converted, k represents the first conversion coefficient used to adjust the proportional relationship of NDVI values ​​of different sensors, and c represents the second conversion coefficient used to adjust the offset relationship of NDVI values ​​of different sensors. The linear transformation relationship that satisfies the preset fit condition and sample size condition is used to transform the time series observations of the corresponding sensor to be transformed, so as to generate time series data under a unified numerical benchmark. Before time series reconstruction, when the candidate images of the main data that meet the validity screening conditions are empty, or when the number of valid observations of the main sensor within the preset time range is lower than the preset number of scenes threshold, low spatial resolution auxiliary remote sensing data is acquired, and the auxiliary remote sensing data is supplemented into the time series data by directly supplementing or by fusing according to the preset fusion weight, depending on the triggering reason. After the initial temporal reconstruction is completed, when the proportion of pixels marked as unreconstructable exceeds a preset proportion threshold, the auxiliary remote sensing data is invoked to supplement the temporal data corresponding to the unreconstructable pixels, so that the reconstruction of the corresponding pixels can be performed again.

5. The method according to claim 1, characterized in that, The step of determining the time-series smoothing parameters from a preset parameter set based on the climate zone of the project includes: The climate type of the project is determined based on the project vector boundary. A discrete candidate set of smoothing parameters λ corresponding to the climate type is obtained from the preset climate type-smoothing parameter mapping relationship. Here, λ represents the smoothing adjustment parameter that controls the relationship between the fidelity of time series observations and the smoothness of curves. Calculate the noise level of the target pixel temporal data. When the noise level is not lower than a preset high noise threshold, remove the minimum λ value from the discrete candidate set; when the noise level is not higher than a preset low noise threshold, remove the maximum λ value from the discrete candidate set. Weighted least squares temporal smoothing is performed using the remaining λ values, and the corresponding root mean square error of reconstruction is calculated through leave-one cross-validation. The λ value that minimizes the root mean square error of reconstruction is selected as the temporal smoothing parameter for the target pixel.

6. The method according to claim 1, characterized in that, The step of performing corresponding observation value synthesis or missing pixel reconstruction based on the number of effective observations of the target pixel within the target time window includes: The number of valid observations N for the target pixel within the baseline annual growing season window is counted, where N represents the number of valid NDVI observations of the target pixel after invalid observations are removed within the baseline annual growing season window. When N≥5, the median of the effective observations is taken as the baseline annual NDVI of the target pixel; When 3 ≤ N < 5, the arithmetic mean of the effective observations is taken as the baseline annual NDVI of the target pixel; When N=2, the difference between two valid observations is calculated. When the difference is not greater than a preset difference threshold, the average of the two valid observations is taken as the baseline annual NDVI of the target pixel. When the difference is greater than the preset difference threshold, the missing pixel reconstruction is initiated. When N=1, a single valid observation is used as the baseline annual NDVI of the target pixel, and a single observation composite status marker is set. When N=0, proceed to the missing pixel reconstruction stage; In the missing pixel reconstruction, effective neighboring pixels that meet the constraints of uncertainty, interpolation state, and NDVI difference are selected from the spatial neighborhood of the target pixel. When the number of effective neighboring pixels is not less than a preset number, the reconstructed value of the target pixel is calculated using an inverse distance weighting method. When the number of effective neighboring pixels is greater than zero but less than the preset number, the effective temporal observations of the years before and after the baseline year of the target pixel and the effective neighboring pixels are used as training data, and the reconstructed value of the target pixel is obtained through Gaussian process regression. When there are no effective neighboring pixels that meet the conditions, temporal smoothing extrapolation is performed based on the effective temporal observations of the years before and after the baseline year of the target pixel. When there are still no effective observations of the target pixel within a preset historical time range, the target pixel is set to an unreconstructable state.

7. The method according to claim 1, characterized in that, The process of calibrating the preset NDVI-land type mapping interval based on the climate zone where the project is located and determining the inferred land type includes: The climate zone where the project is located is determined based on the project vector boundary and the preset climate zone data, and the lower boundary correction coefficient and upper boundary correction coefficient of the corresponding NDVI interval are obtained. The lower and upper boundaries of the preset NDVI-land category mapping interval are calculated with the corresponding correction coefficients to form the regionally calibrated NDVI-land category mapping interval. The NDVI of each pixel in the baseline year NDVI raster image is input into the NDVI-land class mapping interval after regional calibration for pixel-level classification. When the NDVI of the target pixel is within the boundary blur zone of the adjacent mapping interval, neighborhood determination and mixed pixel decomposition are performed on the target pixel to obtain at least one land surface type and its area ratio corresponding to the target pixel. Based on the classification results of each pixel within the project vector boundary, the area proportion of each surface type is calculated, and the land parcel-level inferred surface type is determined based on the area proportion.

8. The method according to claim 1, characterized in that, The step of determining the corresponding carbon sink parameter conversion model and regional parameters based on the project type information obtained from the methodology application document and / or the inferred land surface type includes: The project type is determined based on the project type keywords in the methodological application document; when no project type is identified, the project type is determined based on the area proportion of the inferred land surface type. When the project type corresponds to forestry carbon sinks, a spatial intersection operation is performed between the project vector boundary and the preset ecological forest area vector boundary. The corresponding ecological forest area is determined based on the intersection result, and the first biomass model parameter a and the second biomass model parameter b corresponding to the ecological forest area are obtained from the preset ecological forest area parameter set, according to: Estimated aboveground biomass = exp(a + b × plot baseline year NDVI) Calculate the extrapolated aboveground biomass, where a and b are regional model parameters determined based on the ecological forest area and used to characterize the mapping relationship between the baseline year NDVI and aboveground biomass of the plot; The estimated aboveground biomass is converted based on the carbon content corresponding to the inferred land surface type to obtain the estimated carbon density. When the project type corresponds to a carbon sink type using other biomass conversion models, the conversion model and regional parameters corresponding to the project type are retrieved from the preset project type parameter set, and the corresponding extrapolated baseline accounting parameters are calculated based on the baseline year NDVI of the plot. When the baseline calculation parameter to be verified includes soil carbon density, the estimated soil carbon density is determined from the preset soil carbon density parameter set based on the project surface type, climate zone, and available soil type information, and the corresponding confidence level is set according to the completeness of the input data.

9. The method according to claim 1, characterized in that, The determination of the baseline authenticity verification level includes: The land use consistency deviation δ_land is calculated using the following formula: δ_land=|NDVI_parcel-NDVI_mid| / W_land; Wherein, NDVI_parcel is the annual NDVI of the land parcel baseline, NDVI_mid is the midpoint of the NDVI interval after regional calibration corresponding to the declared baseline surface type, and W_land is the width of the NDVI interval after regional calibration. The units and calculation methods of the estimated baseline calculation parameters and the declared baseline calculation parameters are standardized, and the parameter deviation δ_param is calculated according to the following formula when the declared baseline calculation parameter is greater than zero: δ_param=|P_est-P_decl| / P_decl; Where P_est is the calculated baseline accounting parameter, and P_decl is the declared baseline accounting parameter corresponding to the calculated baseline accounting parameter; when P_decl is zero, the absolute deviation is compared with the corresponding preset absolute threshold. Based on the baseline year NDVI and the current year NDVI of the land parcel, the corresponding carbon density or carbon sink change parameters are calculated and converted into the same change indicators as the declared carbon sink growth parameters, resulting in the remote sensing change rate r_rs and the declared change rate r_decl. The change consistency deviation δ_change is then calculated according to the following formula: δ_change=|r_rs-r_decl| / |r_decl|; Wherein, r_rs represents the remote sensing carbon sink change rate determined by the conversion result of carbon sink parameters corresponding to the baseline year NDVI and the current year NDVI of the land parcel, and r_decl represents the declared carbon sink change rate determined by the declared carbon sink growth parameters in the methodology declaration document; when r_decl is zero, |r_rs-r_decl| is used as the absolute change deviation, and the absolute change deviation is compared with the corresponding preset absolute change threshold; The individual level of each verification dimension is determined based on the preset grading thresholds corresponding to land type consistency deviation, parameter deviation, and change consistency deviation; when the direction of remote sensing change is opposite to the direction of carbon sink growth reported, the change consistency dimension is determined as the highest anomaly level. The comprehensive verification level is determined according to the principle of prioritizing the highest anomaly level; when at least two verification dimensions reach the medium level, the comprehensive verification level is upgraded to the high level; when there are fewer than two valid verification dimensions, the comprehensive verification status is set to undetermined. Based on the comprehensive verification level, output at least one of the following: an abnormal status flag, a data correction request signal, a route review instruction, a verification termination signal, or a verification pass signal.

10. A carbon sink baseline authenticity verification system based on multi-source remote sensing time-series reconstruction, characterized in that, include: The basic parameter extraction module is used to receive methodology application documents, extract baseline year, baseline surface type, project vector boundary and baseline calculation parameters from the methodology application documents, and generate a structured basic parameter object; The multi-source image normalization module is used to acquire multi-source satellite image metadata within the target time window based on the baseline year and project vector boundary. According to the image cloud coverage, the effective coverage ratio within the project vector boundary, the matching degree between the image acquisition time and the target time window, and the availability status of sensor data, it compares each screening index with the corresponding preset validity conditions, retains multi-source satellite images that meet the preset validity conditions, extracts pseudo-invariant feature pixels from the retained images, establishes a linear transformation relationship between different sensors based on the pseudo-invariant feature pixels, normalizes the pixel observation values ​​of non-reference sensors to the same numerical reference, and generates normalized time series data. In addition, when performing change consistency verification, the current year multi-source satellite images are processed according to the same validity screening and normalization rules to generate current year normalized time series data; The temporal reconstruction module is used to perform abnormal pixel identification and mask removal on the normalized temporal data, determine the temporal smoothing parameters according to the climate zone where the project is located, and perform corresponding observation value synthesis or missing pixel reconstruction based on the effective number of observations of the target pixel within the target time window to generate a baseline year NDVI raster image. It also performs plot-scale statistics on the baseline year NDVI raster image according to the project vector boundary to generate the plot baseline year NDVI; and processes the current year normalized temporal data according to the same temporal reconstruction and plot-scale statistics rules as the baseline year to generate the plot current year NDVI. The parameter conversion module is used to calibrate the preset NDVI-land type mapping interval according to the climate zone where the project is located, determine the inferred land type based on the calibrated NDVI-land type mapping interval and the baseline year NDVI raster image, and determine the corresponding carbon sink parameter conversion model and regional parameters based on the project type information and / or inferred land type obtained from the methodology application document, and calculate the baseline accounting parameters based on the baseline year NDVI of the plot. The authenticity verification module is used to determine the land type consistency deviation based on the annual NDVI of the land parcel baseline and the corresponding NDVI interval of the land surface type of the declared baseline, to determine the parameter deviation based on the calculated parameters of the estimated baseline and the calculated parameters of the declared baseline, and to determine the change consistency deviation based on the annual NDVI of the land parcel baseline, the annual NDVI of the current status of the land parcel, and the declared carbon sink growth parameters; to determine the baseline authenticity verification level according to the preset grading rules corresponding to each deviation, and to output the status flag signal and verification control command corresponding to the verification level.

Images