A global sustainable supply chain certification method and apparatus

Abstract

The application discloses a global sustainable supply chain certification method and device, and belongs to the technical field of data management. The method comprises the following steps: performing GIS positioning analysis on a production site to generate a geographic fence area; calculating sustainability quantitative index data of the geographic fence area; obtaining activity data and emission factors with a credibility level from a plurality of data collection channels classified by credibility; calculating the life cycle carbon emission of a product; determining the probability distribution model and the uncertainty parameter followed by the activity data and the emission factors according to a preset mapping relationship between the credibility level and the probability distribution model to generate a confidence interval of the life cycle carbon emission; and verifying the sustainability quantitative index data and the life cycle carbon emission according to the confidence interval. The technical problems of insufficient geographic accuracy, inaccurate carbon emission calculation and unquantifiable result credibility in traditional certification are solved.

Description

Technical Field

[0001] This application relates to the fields of carbon emission accounting and data governance, and data management technology, and in particular to a global sustainable supply chain certification method and equipment. Background Technology

[0002] With increasing global focus on climate change and sustainable development, various sustainability certification standards (such as ISCC and RSPO) play a crucial role in supply chain management. International sustainability and carbon certification systems have become one of the most widely recognized standards in the biomass and bioenergy sector.

[0003] However, traditional certification methods have the following technical limitations: Geographic Information Systems are used for supply chain traceability, but are mostly limited to simple path visualization; geospatial analysis is not accurate enough to achieve precise location of production sites; and life cycle carbon accounting often uses regional or national averages, which cannot accurately reflect the real environmental impact of a specific production site.

[0004] In addition, current accounting practices lack transparency regarding data sources, boundary assumptions, and processing methods. Different institutions have vastly different accounting results for the same product, which provides room for selectively using favorable data to "greenwash" practices. Summary of the Invention

[0005] To address the aforementioned issues, this application proposes a global sustainable supply chain certification method, comprising: performing GIS location analysis on the production location of products within the supply chain to generate geofenced areas; calculating sustainability quantitative index data for the geofenced areas based on multi-source remote sensing data and GIS-LCA spatial analysis methods; acquiring activity data and emission factors with credibility level identifiers from multiple data collection channels with credibility levels; calculating the life cycle carbon emissions of the products using an LCA carbon accounting model; determining the probability distribution models and uncertainty parameters followed by the activity data and emission factors respectively, based on the credibility levels corresponding to the activity data and emission factors, according to a preset mapping relationship between credibility levels and probability distribution models, to generate confidence intervals for life cycle carbon emissions through Monte Carlo simulation; and verifying the sustainability quantitative index data and life cycle carbon emissions for compliance with sustainability standards based on the confidence intervals.

[0006] On the other hand, embodiments of this application provide a global sustainable supply chain certification device, including: at least one processor; and a memory communicatively connected to the at least one processor; wherein the memory stores instructions executable by the at least one processor, the instructions being executed by the at least one processor to enable the at least one processor to perform a global sustainable supply chain certification method as described above.

[0007] The above-described technical solutions adopted in the embodiments of this application can achieve the following beneficial effects:

[0008] By performing GIS location analysis on production sites and generating geofenced areas, the analysis unit is refined from a macro-level region to a micro-level geographic grid. This allows subsequent sustainable quantitative analysis and the extraction of activity data and emission factors to be precisely spatially linked to specific production activities, eliminating errors caused by overgeneralization of geographical location ranges at the source.

[0009] A credibility-tiered data collection channel has been established. The system can intelligently select the optimal available data source and assign a clear credibility level label to the data. This mechanism acknowledges and systematically manages the objectively existing data quality differences in the global supply chain. Under the reality of incomplete data, it produces the most credible carbon footprint accounting results possible through a layered, degradable, and transparently quantifiable data governance and computing framework, and significantly improves the transparency of the accounting process and the traceability of the data. Attached Figure Description

[0010] To more clearly illustrate the technical solution of this application, some embodiments of this application will be described in detail below with reference to the accompanying drawings, in which:

[0011] Figure 1 This is a flowchart illustrating a global sustainable supply chain certification method provided in an embodiment of this application. Detailed Implementation

[0012] To make the objectives, technical solutions, and advantages of this application clearer, the technical solutions of this application will be clearly and completely described below in conjunction with specific embodiments and corresponding drawings. Obviously, the described embodiments are only a part of the embodiments of this application, and not all of them. Based on the embodiments in this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0013] Some embodiments of this application will now be described in detail with reference to the accompanying drawings.

[0014] Figure 1This is a schematic diagram of a global sustainable supply chain certification method provided in an embodiment of this application. Certain input parameters or intermediate results in this process can be manually adjusted to help improve accuracy.

[0015] Figure 1 The process includes the following steps:

[0016] S101: Perform GIS location analysis on the production location of products in the supply chain to generate geofence areas.

[0017] In this example, a precise, dynamically adjustable spatial boundary (geofence) is defined for the product's production location in order to accurately bind the product to its specific geographical environment and production activities.

[0018] For example, the Voronoi diagram algorithm can be used to generate polygons, and geofenced areas of production sites can be determined based on polygon regions.

[0019] In summary, by expanding the point location of the production site into a geographically significant area, the scope of subsequent refined calculations can be defined.

[0020] S102: Calculate the sustainability quantitative index data of the geofence area based on multi-source remote sensing data and GIS-LCA spatial analysis method.

[0021] In this example, multi-source remote sensing data refers to Earth observation data acquired from various types of remote sensing platforms or sensors, obtained through open and standardized international data platform APIs. Using the geofence boundary coordinates (typically polygons in the WGS84 coordinate system) and a preset time range (such as the past 5 years for historical analysis) as input parameters, the system automatically retrieves and downloads raw data products covering that region and time period.

[0022] For example, data sources include: Optical remote sensing satellites (such as Landsat and Sentinel-2): providing high-resolution multispectral imagery for calculating the Normalized Difference Vegetation Index (NDVI), and serving as core data for monitoring forest cover, vegetation health, and land use change. Radar satellites (such as Sentinel-1): capable of penetrating clouds and some vegetation, unaffected by day / night cycles or weather, and can be used to supplement optical data, monitoring surface moisture, topography, and certain vegetation structures. Meteorological satellites: providing climate data such as temperature, precipitation, and solar radiation, used to assess climate risk and water resource pressure.

[0023] It should be noted that all the remote sensing data mentioned above were stored and transmitted in GeoTIFF format. This is a TIFF image format that includes embedded geographic coordinate information, ensuring that each pixel accurately corresponds to a specific geographic location on the Earth's surface.

[0024] Optical remote sensing data can originate from USGS's Landsat series satellites (such as Landsat 8 / 9 OLI) and ESA's Sentinel-2 satellites under the Copernicus program. Acquisition can be achieved through the USGS EarthExplorer API or ESA's Copernicus Open Access Hub API, using the boundary coordinates (WGS84 coordinate system) of the geofenced area and a time range (such as a preset historical duration) as query criteria to programmatically download Level-2 scientific products (such as Landsat's Surface Reflectance products and Sentinel-2's L2A atmospheric bottom reflectance products) for the corresponding area. The data format of optical remote sensing data can be multi-band GeoTIFF image files, each containing reflectance information for a specific band and complete geographic metadata (projection information, pixel size, coverage area).

[0025] Radar remote sensing data can be sourced from ESA's Sentinel-1 satellite, acquired via the Copernicus Open Access Hub API by downloading ground distance detection (GRD) products in interferometric wide swath (IW) mode as needed. These products are radiometrically calibrated and topographically corrected. Data with an appropriate polarization scheme (e.g., VV+VH) is selected based on the geofenced area and date range. The data is in GeoTIFF format image files, including backscattering coefficient information for polarization channels such as vertical transmission / vertical reception (VV) and vertical transmission / horizontal reception (VH).

[0026] Weather satellite data can be obtained from the climate dataset integrated into the Google Earth Engine (GEE) platform, which provides gridded climate data corresponding to the geofenced area, including but not limited to daily / monthly average temperature (e.g., at a height of 2 meters), cumulative precipitation, and surface solar radiation flux.

[0027] Specifically, the data structure and content can be as follows: Optical Imagery: Each COG file is a multi-band two-dimensional array. For example, the Sentinel-2 L2A product contains 13 bands, and the pixel value in each band array represents the surface reflectance of that pixel at a specific wavelength (the unit is a dimensionless scaling factor, typically ranging from 0 to 1 or scaled by 10000). These reflectance values ​​are the physical basis for calculating derived indicators such as the Normalized Difference Vegetation Index (NDVI), calculated using the formula NDVI = (B8 - B4) / (B8 + B4), where B8 is the near-infrared reflectance and B4 is the red band reflectance. Radar Imagery: Sentinel-1 GRD products are typically dual-band (VV, VH) COG files. Each pixel value represents the backscattering coefficient, usually expressed in decibels (dB). This data is used to monitor changes in surface moisture or identify specific surface structures under cloudy conditions.

[0028] Meteorological data: ERA5 and other data are provided in the form of univariate, multi-time-slice gridded data. Each COG file represents the geospatial distribution of a climate variable (such as temperature) at a specific time, with pixel values ​​representing physical quantities (such as Kelvin temperature, wind speed in meters per second, and precipitation in millimeters). Metadata: Each data scene is accompanied by complete metadata, clearly recording its spatial reference system, geographic coordinate boundaries of the coverage area, pixel spatial resolution, timestamp of data acquisition or simulation, and data quality identifiers (such as cloud mask files). This metadata is the key basis for the system to automatically perform spatiotemporal matching and data fusion.

[0029] In summary, the data format is scientific-grade raster data (GeoTIFF) with precise geocoding, and the data content consists of physically calibrated quantitative measurements of surface attributes (such as reflectance, backscattering coefficient, and radiance). Those skilled in the art can reliably acquire and read this data using existing Geographic Information Science (GIS) and remote sensing data processing tools (such as the GDAL library, QGIS, and Google Earth Engine SDK) based on the publicly available data sources, product levels, and formats described above, providing input for subsequent quantitative indicator calculations.

[0030] The GIS-LCA spatial analysis method is an analytical framework that deeply integrates the spatial data processing, analysis, and visualization capabilities of Geographic Information Systems (GIS) with Life Cycle Assessment (LCA) models. Its core objective is to address the problem of traditional LCA using global average data and neglecting spatial heterogeneity. It utilizes GIS tools to spatially process, extract features, and model data from remote sensing and other sources, and then uses the analysis results to serve international sustainability certification standards. Specifically, GIS provides the necessary spatial data processing, analysis, and visualization capabilities for computation (e.g., fence delineation, remote sensing image processing, kernel density, hotspot analysis), while LCA sets the computational objectives and framework (whether it meets the requirements of standards such as ISCC for environmental and social sustainable development), and transforms the results of GIS analysis into quantitative indicators that can be used for certification decisions.

[0031] Among them, the quantitative indicators of sustainability can include historical data on deforestation, social risk data, etc.

[0032] Deforestation history reflects environmental baselines, ensuring that economic activity does not come at the expense of critical ecosystems. Social risks reflect social baselines. Climate risks reflect governance and long-term sustainability, assessing the future resilience of supply chains.

[0033] Forest deforestation history analysis refers to determining whether deforestation or ecosystem degradation has occurred within a designated geographic fenced area over a certain period (e.g., 5 years). This means that the area cannot be newly planted after deforestation. If large-scale, continuous vegetation loss is found within the geographic fence, it is considered a risk of deforestation. Time-series satellite remote sensing imagery, such as Landsat (30-meter resolution), Sentinel-2 (10-meter resolution), or even higher resolution commercial satellite data, can be used. Interannual variations can be analyzed using NDVI (Normalized Difference Vegetation Index) or more specialized forest cover change maps.

[0034] Social risk analysis refers to assessing whether there are serious social risk factors in a geofenced area.

[0035] For example, data from multiple sources can be collected, including population distribution within geofenced areas, news and public opinion data, and historical events. Machine learning models (such as classifiers) can then be used to integrate this multi-source data to generate a social risk heatmap. Hotspot areas on the social risk heatmap can represent high-risk areas.

[0036] In conclusion, by utilizing GIS and big data analytics, the abstract rules in international certification standards can be transformed into quantifiable spatial judgments, thereby greatly improving the reliability, efficiency, and credibility of certification.

[0037] S103: Obtain activity data and emission factors with credibility level identifiers from multiple data collection channels with credibility level classification.

[0038] In this example, the theoretical accounting model requires high-precision input data, but data availability, accuracy, and monitoring levels vary greatly. For instance, large, modern factories may be equipped with complete IoT sensors and automatic metering systems, while scattered small and medium-sized farmers or collection points may rely solely on manual records or even lack systematic records altogether. Using a single, coarse estimate would sacrifice the reliability of the results. Therefore, a hierarchical data input system was constructed, and a multi-data collection channel with tiered reliability was proposed, enabling reliable accounting results even under conditions of incomplete data.

[0039] It's important to note that different data collection channels correspond to different credibility levels. In other words, determining the specific credibility level of activity data or emission factors depends on the data collection channel from which the data or emission factors originate. Data classification primarily relies on the verifiability of its source, the objectivity of the collection method, and the inherent error level. In non-transportation stages, geographical location is a key index for triggering and selecting appropriate data channels. Therefore, by assigning clear credibility level labels to data of different qualities, subsequent accounting can recognize and quantify the uncertainty of the input data, thereby outputting highly transparent accounting results with confidence intervals.

[0040] In life cycle assessment, the emission factor is a core parameter used to represent the amount of greenhouse gas emissions corresponding to a unit of activity data (such as consuming 1 kilowatt-hour of electricity or burning 1 liter of diesel). Given the same activity data, using a higher emission factor will result in higher greenhouse gas emissions, indicating that the activity has a greater potential environmental impact on climate change.

[0041] It should be noted that activity data includes material data and energy consumption data. Material data refers to the physical quantities of raw materials, intermediate products, and the products themselves involved in the production, processing, and transportation stages. Examples include raw material input (e.g., the tons of waste oil consumed in producing a batch of biodiesel, the kilograms of seeds consumed in planting crops), product output (e.g., the liters of biodiesel produced in a factory, the tons of grain harvested from a farm), and material transfer / loss (e.g., material loss during transportation or storage, the amount of by-products generated during processing). Energy data refers to the quantity of energy and resources directly consumed to complete production, processing, and transportation activities. Examples include energy consumption (e.g., electricity, natural gas, and coal consumed in factory production), resource consumption (e.g., water resources consumed in agricultural irrigation, steam consumed in industrial production), and transportation fuel consumption (e.g., diesel fuel and marine fuel oil consumed by trucks and ships transporting materials).

[0042] In this example, the activity data and emission factors each have their own data collection channel classification.

[0043] For example, the data collection channel for the first credibility level of activity data refers to the data that is real-time monitored and directly uploaded through Internet of Things sensors and automatic metering devices (such as weighbridges and flow meters). The data collection channel for the second credibility level refers to the operational business data from enterprise resource planning (ERP), production management systems, or third-party certified operation ledgers (such as purchase invoices, logistics bills of lading, energy bills). The data collection channel for the third credibility level refers to the estimation based on the geographical fence area of the production site, by invoking geographical information data (such as the cultivated area evaluated by satellite images and the catering scale estimated by POI data), and through a built-in model. It can be seen that the data collection channel for the third credibility level is used for supplementing data-poor scenarios. It should be noted that for the transportation link, there is no such data collection channel for the third credibility level.

[0044] For example, for non-transportation links, the data collection channel for the first credibility level refers to querying and matching the localized, specific facility or park emission factor list officially released by the country or local area based on the precise geographical coordinates of the production site. For example, matching the specific steam emission factor of a chemical industrial park where a biodiesel plant is located. The data collection channel for the second credibility level refers to obtaining the baseline factor of the region (such as the provincial power grid average factor) according to the administrative division of the province or city where the production site is located, and further using the characteristic parameters of the region (such as the municipal power supply structure and industrial composition), and making a localized adjustment through a preset standardized correction model to generate a corrected factor. The data collection channel for the third credibility level refers to using the average factor in a national or internationally recognized database as the default value.

[0045] For the transportation link, the data collection channel for the first credibility level refers to querying and matching the vehicle carbon emission factor list officially released based on the specific vehicle model. The data collection channel for the third credibility level refers to the average factor based on the vehicle type. It should be noted that there is no data collection channel for the second credibility level.

[0046] Based on this, for the data collection channels of activity data and emission factors, the first credibility level is higher than the second credibility level, and the second credibility level is higher than the third credibility level. That is, the first credibility level is equivalent to high credibility, the second credibility level is equivalent to medium credibility, and the third credibility level is equivalent to reference credibility.

[0047] It should be noted that if activity data or emission factors can be collected from the first level of confidence data collection channel, the collection action from the second level of confidence data collection channel will not be performed. If data can be collected from the second level of confidence data collection channel, the collection action from the third level of confidence data collection channel will not be performed.

[0048] In this example, each piece of activity data and emission factor acquired carries an identifier of its level, which allows for uncertainty quantification and propagation (such as Monte Carlo simulation).

[0049] It should be noted that different confidence levels are associated with preset uncertainty parameters (such as error ranges). For example, the uncertainty of data at the first confidence level might be preset to ±2%, while at the third confidence level it might be ±25% or higher. This allows the final carbon footprint results to be presented in the form of a numerical ± range, and also supports the generation of a comprehensive data quality score.

[0050] Based on this, the tiered data acquisition mechanism, especially in the production and raw material collection stages, enables precise accounting tailored to local conditions and data availability. It can produce high-precision results when data conditions are favorable, and also provide estimates with a clear range of uncertainty when data is scarce, guiding enterprises to improve to higher data levels and enhancing the operability and credibility of carbon accounting in the global supply chain.

[0051] S104: Calculate the life cycle carbon emissions of the product using the LCA carbon accounting model.

[0052] In this example, the LCA carbon accounting model is an application in the carbon emission calculation scenario. Its core function is to systematically quantify the carbon footprint of a product from its birth to its demise.

[0053] Lifecycle scope can be flexibly defined depending on the certification standards and assessment objectives: some standards (such as the complete product carbon footprint standard) require calculation of emissions across the entire chain, covering all stages including production, transportation, use, and disposal; while other standards or specific application scenarios may only focus on calculating emissions during the production stage up to the point where the product leaves the factory. However, regardless of the lifecycle scope adopted, carbon emissions during the production stage are a core component that must be calculated.

[0054] It should be noted that for each consumption or output, carbon emissions are calculated as the product of the activity data and the emission factor. The emission factor represents the carbon emission intensity corresponding to a unit of activity data.

[0055] For the production process, activity data are specific quantitative values ​​of material or energy flows that are linked to geofences (such as cubic meters of natural gas consumed, kilowatt-hours of electricity generated, and tons of products produced).

[0056] Based on this, each activity data (such as a factory consuming 1,500 cubic meters of natural gas) is matched with an emission factor that best represents the geographical or technical characteristics (such as the specific emission factor of the natural gas pipeline in the city where the factory is located).

[0057] Based on this, the above multiplication calculation is performed on all relevant activities in the corresponding stage, and then the carbon emissions of all activities are summed to obtain the total carbon emissions of that stage. Over the entire lifecycle of the supply chain, the total carbon emissions of all stages are summed to obtain the lifecycle carbon emissions.

[0058] In summary, combining the globally accepted LCA carbon accounting framework with a localized and differentiated data input system based on high-precision geospatial analysis enables carbon footprint accounting to move from macro-average to micro-precision.

[0059] S105: For activity data and emission factors with confidence level identifiers, determine the probability distribution models followed by the activity data and emission factors and the corresponding uncertainty parameters according to the preset mapping relationship between confidence level and probability distribution model, so as to generate confidence intervals for life cycle carbon emissions through Monte Carlo simulation.

[0060] In this example, instead of simply outputting a single carbon emission estimate, the uncertainty of the input data is quantified to ultimately present a statistical range that reflects the reliability of the results. This provides a scientific basis for high-standard certifications such as ISCC, ensuring that risks are known and decisions are robust.

[0061] Based on this, a corresponding uncertainty parameter is preset for each confidence level (e.g., Level 1, Level 2, and Level 3) of the activity data and emission factors. This parameter quantitatively describes the possible error range or probability distribution characteristics of the data at that level. The uncertainty parameter for Level 1 is smaller than that for Level 2, and the uncertainty parameter for Level 2 is smaller than that for Level 3.

[0062] For example, the first confidence level corresponds to data from direct monitoring or precise geographic matching, with low uncertainty, possibly preset to ±2% (the standard deviation of a normal distribution). The second confidence level corresponds to verified business records or model-corrected data, with moderate uncertainty, possibly preset to ±10%. Data corresponding to general average factors or geographically based estimations has higher uncertainty, possibly preset to ±25%. This mapping relationship is based on a comprehensive preset of historical data verification, measurement technology error analysis, and industry standards (such as GB / T32150), forming the basis for uncertainty quantification.

[0063] Specifically, high-confidence data at the first confidence level (such as directly monitored data) has small errors and is usually symmetrically distributed around the true value, so a normal distribution is suitable. Medium-confidence data at the second confidence level (such as business record data) has some errors but still shows a central tendency, so a normal distribution is also suitable, but with a larger standard deviation. Low-confidence data (such as estimated data) has a large error range and may not have a clear central value; therefore, a uniform distribution (the most conservative assumption) is suitable.

[0064] The pre-defined mapping relationship between confidence levels and probability distribution models represents the probability distribution model corresponding to the confidence level and its uncertainty parameters (such as standard deviation). For example, the first confidence level corresponds to a normal distribution with a pre-defined standard deviation of 2%. The second confidence level corresponds to a normal distribution with a pre-defined standard deviation of 10%. The third confidence level corresponds to a uniform distribution with a pre-defined range of ±25%.

[0065] In this example, the process of simulating and generating confidence intervals for lifecycle carbon emissions may include the following steps:

[0066] First, based on the confidence level of each data point, the corresponding probability distribution model and corresponding uncertainty parameters are found in the mapping table.

[0067] Then, based on the activity data values, emission factor values, and corresponding uncertainty parameters, probability distributions corresponding to the activity data and emission factors are constructed respectively.

[0068] For example, construct a corresponding probability distribution model, such as a normal or log-normal distribution centered on the reported value and with the uncertainty parameter as the standard deviation. Or, center on the baseline value of the data (for a uniform distribution).

[0069] Specifically, for calculating the distribution parameters, for a normal distribution, the mean equals the reported value, and the standard deviation equals the reported value multiplied by a preset percentage. For a uniform distribution: the lower limit equals the reported value multiplied by (1 - preset percentage), and the upper limit equals the reported value multiplied by (1 + preset percentage).

[0070] Then, multiple random samples are taken from the probability distribution, and the activity data values ​​and emission factor values ​​of each sample are calculated to obtain the probability distribution of life-cycle carbon emissions. For example, Monte Carlo simulation is used as the core engine to perform uncertainty propagation, as follows:

[0071] Sampling: From the probability distribution of each activity data and emission factor constructed above, a large number (e.g., tens of thousands) of random independent samples are taken. Each sampling forms a set of possible combinations of input data.

[0072] Calculation: For each combination of data obtained from sampling, the life cycle carbon emissions are recalculated according to the carbon emission formula to obtain a possible result.

[0073] Aggregation: After repeating the above process thousands of times, thousands of possible results will be obtained, which together constitute the probability distribution of life cycle carbon emissions.

[0074] Finally, based on the probability distribution of life-cycle carbon emissions, the corresponding confidence intervals are determined.

[0075] In this example, based on the carbon emission probability distribution generated by the Monte Carlo simulation above, the system performs statistical analysis and generates a specified confidence interval.

[0076] For example, a 95% confidence interval means that there is a 95% probability that the actual carbon emissions fall within this interval. The calculation method typically involves finding the 2.5% and 97.5% quantiles in the distribution of simulation results, which are used as the lower and upper limits of the interval, respectively.

[0077] S106: Based on the confidence interval, verify the sustainability quantitative indicator data and life cycle carbon emissions for compliance with sustainability standards.

[0078] Among them, the sustainability standard compliance verification can be in accordance with the ISCC standard for international sustainability and carbon certification.

[0079] This example changes the traditional compliance verification model that relies on a single, definite numerical value for a binary yes / no judgment. By introducing a probability-based decision framework, the inherent uncertainty of the data is made explicit and incorporated into the certification logic, thereby greatly improving the scientific rigor, robustness, and credibility of ISCC compliance verification.

[0080] Since the confidence interval for carbon emissions is a quantitative range generated through Monte Carlo simulation, it represents a probabilistic understanding of the true value of carbon emissions given the current data quality. Therefore, risk-based compliance decisions can be implemented based on the relative position of the confidence interval to a preset threshold (completely below, exceeding, or completely above), allowing for expedited approval, focused review, rejection, or requirements for improvement. For example, the specific logic for risk-based compliance assessment of carbon emissions is as follows:

[0081] The upper limit of the confidence interval is lower than the preset threshold. This indicates that even under the most unfavorable combination of errors (i.e., carbon emissions reaching their maximum), the product still meets carbon emission requirements. The system determines this dimension to be low-risk compliance and can quickly proceed to subsequent comprehensive judgment. This greatly improves the efficiency of certifying a high-quality data supply chain.

[0082] The lower limit of the confidence interval is higher than the preset threshold, indicating that even under the most favorable error combination (i.e., minimum carbon emissions), the product is almost impossible to meet the standards. This dimension is deemed a high-risk non-compliance, allowing for a direct rejection or a requirement for fundamental process modifications. This avoids wasting audit resources on clearly non-compliant projects.

[0083] The confidence interval crosses a preset threshold (i.e., the lower limit is below the threshold, but the upper limit is above the threshold). This indicates that the uncertainty of the calculation results may lead to different conclusions. The system marks this as a critical case requiring close review. One or more of the following mechanisms will be automatically triggered:

[0084] Request higher quality data: prompt companies to provide activity data or emission factor evidence at a higher level (e.g., from the third confidence level to the second confidence level).

[0085] Initiate a special on-site verification: Mark this node in the certification process and arrange on-site verification of the authenticity of key data.

[0086] Use conservative estimates for provisional judgments: When new data cannot be obtained immediately, the upper limit of the confidence interval can be used for judgment (i.e., consider the worst case).

[0087] Based on this, the certification resources were precisely allocated, manual review was focused on the areas where questions truly arose, and clear data improvement paths were provided for enterprises.

[0088] After completing the risk classification assessment for carbon emissions, and if the risk assessment result is low-risk compliance, the system will perform the final joint verification of ISCC compliance.

[0089] pass Figure 1 Traditional methods, which use administrative divisions or rough coordinates, cannot accurately link production activities with specific ecological environments. This solution uses GIS location analysis of production sites and generates geofenced areas, refining the analysis unit from a macro-regional level to a micro-geographic grid. This allows subsequent sustainable quantitative analysis and the extraction of activity data and emission factors to be precisely spatially linked to specific production activities, eliminating errors caused by overgeneralization of geographical location ranges at the source, and laying an immutable spatial foundation for accurate accounting and verification.

[0090] Traditional carbon footprint accounting faces a conflict between ideal models and a lack of real-world data. This solution establishes a tiered data collection channel based on credibility. The system intelligently selects the optimal available data source and assigns a clear credibility level to the data. This mechanism acknowledges and systematically manages the objectively existing data quality differences in the global supply chain. Under conditions of incomplete data, it produces the most credible carbon footprint accounting results possible through a layered, degradable, and transparently quantifiable data governance and computing framework. This ensures that accounting can be initiated under any data conditions and significantly improves the transparency of the accounting process and the traceability of the data. It does not require complex models (such as machine learning) to construct mathematical relationships between emission factors and does not rely on large amounts of high-quality training data, making it easy to implement in real-world scenarios with weak data foundations.

[0091] Even with the existence of GIS-based supply chain traceability systems, most are limited to simple geographic location marking and fail to achieve deep coupling between GIS data and LCA models. Traditional LCA calculations often use national or regional average emission factors. This solution moves beyond simple table lookups in determining emission factors. Especially for the production and processing stages, a localized decision path can be triggered based on the production location coordinates, dynamically coupling localized features extracted based on geofencing with the LCA model. Furthermore, activity data can be integrated with local geographic features for estimation or calibration. This allows the final calculated carbon emissions from the production stage to truly reflect the unique environmental impact of a specific production location due to differences in resource endowment, energy structure, and technological conditions, resulting in a qualitative leap in calculation accuracy.

[0092] Traditional certification methods rely on a single numerical value to make a binary judgment of pass / fail, which fails to reflect the reliability of the data itself. This solution uses Monte Carlo simulation based on confidence levels to generate confidence intervals for carbon emissions, explicitly and quantifying the uncertainty of data quality and conveying it to the final result.

[0093] Finally, combining sustainability (risk screening) and carbon footprint accounting for certification reveals the intrinsic link between risk and carbon emissions and improves the accuracy of compliance verification. This collaborative analysis supports the development of comprehensive emission reduction and risk avoidance strategies. Through multi-dimensional compliance verification, automated certification decisions are achieved, and reliable certification results are output.

[0094] It should be noted that the scope of this solution covers multiple stages of the product lifecycle, and the activity data and emission factors for each stage can be obtained in layers and the uncertainty can be quantified.

[0095] based on Figure 1 In addition to the method described herein, this application also provides some specific implementation schemes and extended schemes of the method, which will be further described below.

[0096] This example addresses the shortcomings of traditional supply chain traceability systems, which are susceptible to data tampering and document loss. It deeply integrates key data reported by each participating node in the supply chain with a blockchain-based immutable evidence storage system, achieving full-chain, verifiable, and tamper-proof data traceability from raw material sources to the final product. Based on this, the specific process may include the following steps:

[0097] First, digital identities are created and registered for each participating node in the supply chain. Then, material data digitally signed by each participating node is received, and the verified data is recorded as a transaction in a distributed ledger to generate a material balance data chain.

[0098] It should be noted that within a specific system (such as a factory or warehouse), the correspondence and balance between the quantity of material inputs and outputs is key to verifying the consistency and authenticity of the data.

[0099] Based on this, by continuously and immutably recording these material data, a material balance data chain that runs through the entire chain can be automatically constructed and verified, ensuring that any claimed sustainable product (such as biodiesel) can be matched in physical quantity with the source of its sustainable raw materials (such as waste oil), fundamentally eliminating the risk of data falsification and document loss, and providing atomic-level credible evidence for compliance verification of standards such as ISCC.

[0100] Then, at least one of the following is identified as a sustainability document: sustainable quantitative indicator data, life cycle carbon emission data, and geographic information of the geofence area. The first smart contract deployed on the blockchain is invoked. The first smart contract is configured to automatically verify the validity of the digital signature and the compliance of the data format of the sustainability document. After the verification is passed, the hash value and metadata of the document are recorded on the blockchain.

[0101] Then, after the product's sustainability standard compliance verification is passed, a second smart contract deployed on the blockchain is invoked. The second smart contract is configured to: collect on-chain evidence records related to product certification, including material balance chain records, sustainability document hashes, compliance verification conclusions, and generate a unique certification evidence chain identifier and a certification digital certificate including the certification evidence chain identifier, product information, and certification conclusions.

[0102] Finally, the authentication evidence chain identifier and the hash value of the authentication digital certificate are submitted to an external trusted evidence storage system for storage.

[0103] Specifically, after compliance verification is passed, all key data and results are automatically generated into a data package and written to the blockchain via a smart contract, achieving end-to-end traceability and tamper-proof evidence storage. For example, the generated digital certificate includes a unique certificate ID, product / batch number, manufacturer information, certification standard (ISCC), validity period, and final certification conclusion (pass / fail). The digital certificate is primarily proof of the certification conclusion and is an output of blockchain-based evidence storage.

[0104] Furthermore, lifecycle carbon emission results can be pushed to a dynamic carbon footprint traceability dashboard, providing visualized decision support for users in different roles. For example, Sankey diagrams or stacked bar charts can be used to display the carbon emission percentage of a single product or batch at each stage of its lifecycle (production, transportation, processing, and disposal). Users can view the detailed calculation process and activity data for emissions at each stage, as well as the credibility level of emission factors.

[0105] In one example, performing GIS location analysis on the production location of products in the supply chain to generate geofenced areas includes the following steps:

[0106] Voronoi diagrams are a spatial partitioning method whose core rule is that every location within a point on a plane is closer to that point than to any other point. Voronoi diagrams solve the problem of how a set of discrete points can divide an entire geographic area into several regions, naturally and fairly assigning each certified entity a unique geographic area of ​​influence.

[0107] Specifically, obtain the geometric center coordinates or precise geographic coordinates of the production site. Use these coordinates as the seed point for the Voronoi algorithm.

[0108] It should be noted that, due to the diversity of geographical information about production sites, for production sites with clear and unique geographical coordinates (such as the place of registration or the geographical location of the factory), the precise geographical coordinates can be obtained directly. This point itself is clear and unique.

[0109] For geographically dispersed production areas without a single precise coordinate (e.g., planting areas composed of numerous farmer plots), where a single point cannot accurately represent the location of all scattered plots, a surrogate strategy is employed to determine the seed point. Specifically, the geometric center of the bounding polygon of all scattered plots can be calculated, and the coordinates of this center are used as the seed point representing that production area. This allows for the allocation of a dedicated influence area and data analysis scope (i.e., the subsequently generated geofence area) for that production area in spatial location analysis.

[0110] It should be noted that when the production scale (area, annual output, etc.) of a production site exceeds a preset scale threshold, it is necessary to first find precise geographic coordinates. This is because entities with larger production scales have a greater environmental impact and require higher data accuracy; therefore, it is essential to perform the operation of obtaining precise geographic coordinates to generate the most representative geofence, ensuring the highest credibility of subsequent analysis and authentication results. However, performing the operation of obtaining precise geographic coordinates does not guarantee that precise geographic coordinates will be obtained.

[0111] Then, the production area is divided into Voronoi polygons to obtain the initial geofence area.

[0112] It should be noted that the boundaries of the polygon are determined by other geographical features near the production area or by a predefined calculation range. Other geographical features include neighboring farms, towns, or natural boundaries.

[0113] For example, the seed point set formed by combining other geographical features is ( ), The coordinates of the production location are given. Based on this, the following calculation is performed on each location (which can be considered as a pixel or an infinitesimally small point) within the target area (e.g., a rectangular area surrounding these seed points):

[0114] First, calculate the straight-line distance (or spherical distance) from the location to the four seed points. Then, compare these distances and find the seed point that is closest to the location. Finally, assign the location to that closest seed point.

[0115] Repeat the comparison process until all locations within the target area have been assigned. All locations assigned... The position constitutes The Voronoi elements (i.e., polygons). Similarly, they are assigned to... Their positions constitute their respective units.

[0116] Based on this, the seed points are extracted. All the locations, when connected, form a closed polygon, which is the initial geographic fence area of ​​the production site.

[0117] Then, the hexagonal grid is overlaid onto the initial geofence area, and the intersection area ratio of each hexagonal grid with the polygon is calculated.

[0118] It should be noted that geographic grids are spatial data structures based on hierarchical partitioning of rules. They achieve spatial positioning and data organization by discretizing the Earth's surface using polygonal grid cells. In geography, based on the discretization of geographic information, geographic space is divided into hexagonal grids, with each hexagonal cell representing a geographic area within a certain range.

[0119] Finally, hexagonal grids with significant spatial overlap with the polygon are selected to determine the geofencing area; significant spatial overlap is determined by comparing the ratio of the intersection area to a preset ratio threshold. It should be noted that this preset ratio threshold ensures that the selected grids effectively represent the original polygon area. Selecting hexagonal grids with significant spatial overlap with the polygon means selecting hexagonal grids whose intersection area ratio is higher than the preset ratio threshold.

[0120] For example, with a total grid area of ​​0.5 square kilometers, hexagonal grid A has an overlap area of ​​0.45 square kilometers with the polygon, resulting in an intersection ratio of 90%. Hexagonal grid B has an overlap area of ​​0.2 square kilometers, resulting in an intersection ratio of 40%. Hexagonal grid C is completely outside the polygon, resulting in an intersection ratio of 0%. If the preset ratio threshold is 50%, then grid A is included in the geofence area. All selected hexagonal grids together constitute a continuous, regular spatial range, which is the precise geofence of the farm.

[0121] In other words, when the overlap between the hexagonal grid and the Voronoi polygon is large enough, the hexagonal grid is included in the final geofence area. A geofence is a clearly defined spatial area that serves as the baseline geographic unit for all subsequent environmental data analysis, sustainability quantification calculations, and localized carbon emission accounting. It should be noted that it is also possible to select hexagonal grid units whose center points all fall within the Voronoi polygon.

[0122] In summary, polygon division is entirely determined by the product's geographical location, avoiding the unfairness or inaccuracy of subjective human division. Furthermore, a fast query index for the nearest point from different points is established. To find the affiliation of a point, the system does not need to calculate the distances to all points; it only needs to perform a rapid spatial relationship judgment. Directly performing statistical calculations on an irregular polygon is inefficient and inaccurate. Therefore, based on the above, spatial analysis is used to select hexagonal grid cell sets that can effectively represent Voronoi polygon regions.

[0123] In one example, calculating quantitative indicators of sustainability for geofenced areas may include the following steps:

[0124] In this example, when the sustainability quantification indicator is the history of deforestation, the sustainability quantification indicator data for the geofenced area is calculated, specifically including:

[0125] First, the normalized vegetation index (NDVI) of each pixel in the satellite imagery of the geofenced area within a preset historical time period and the current time imagery is calculated to obtain the NDVI distribution map corresponding to each image.

[0126] The calculation formula is: NDVI=(B8-B4) / (B8+B4), where NDVI is the normalized vegetation index, B8 is the near-infrared reflectance, and B4 is the red reflectance.

[0127] It should be noted that satellite imagery covering geofenced areas can be retrieved from the satellite data center. For example, Sentinel-2 L2A-level surface reflectance images of the same period each year within a preset historical timeframe (e.g., 5 years) covering the geofenced area can be obtained. Each image can be radiometrically calibrated and atmospherically corrected to eliminate the effects of atmospheric, illumination, and sensor differences, ensuring fair data comparison.

[0128] Then, for each specified pixel in the NDVI distribution map of the current time image, the NDVI value from each historical image is extracted to obtain the NDVI time series for each specified pixel in the NDVI distribution map of the current time image. For example, an NDVI time series vector of length n is generated for each specified pixel. In other words, what is obtained for each specified pixel is an NDVI curve that fluctuates over time.

[0129] Then, the NDVI time series of each specified pixel is analyzed according to the preset deforestation detection network model to obtain the deforestation labeling result of each specified pixel.

[0130] It should be noted that deforestation data can include whether a pixel was deforestation, the time of deforestation, and the extent of deforestation. For example, the time of deforestation was July 2022, and the magnitude of change was a decrease of approximately 0.6 in NDVI. The state after the change was from forest to non-forest.

[0131] A sudden and permanent drop in NDVI values ​​typically occurs when a forest transitions to a non-forested state. For example, we can look for a point in time after which the overall NDVI level experienced a precipitous decline, and whether this decline was permanent or only temporary (e.g., recovery after harvesting or fire).

[0132] The pre-defined forest deforestation detection network model can distinguish between normal seasonal fluctuations (such as natural yellowing of vegetation during the dry season) and abnormal permanent declines (such as deforestation). It can capture abrupt patterns of precipitous drops in vegetation cover and accurately pinpoint the time when the changes occur.

[0133] The training process for a deforestation detection network model can be as follows: using the NDVI time series of each sample pixel as input and deforestation data as sample labels, train the neural network model to obtain the deforestation detection network model. For example, adopting a Transformer network model architecture can capture long-distance dependencies through a self-attention mechanism, improving the accuracy of change detection. This could be a temporal classification model based on a convolutional neural network and an attention mechanism. This model sequentially includes: an input layer: receiving the NDVI time series; a one-dimensional convolutional layer: using multiple convolutional kernels of different sizes (e.g., lengths 3 and 5) to extract local temporal features; a temporal attention layer: calculating the importance weights of features at different time points in the sequence, enabling the model to focus on key time regions where abrupt changes occur; and a fully connected classification layer: the final output layer uses a Sigmoid activation function (corresponding to probability output) or a Softmax function (corresponding to binary classification). If a probability output is used, a classification probability threshold is set (e.g., a classification probability of 0.7). When the classification probability threshold is exceeded, the pixel is determined to be a deforestation pixel; otherwise, it is a non-deforestation pixel. If a direct classification method is used, then the label 1 output by the model is used directly as the determination of the deforestation pixel.

[0134] Then, if the ratio of the area marked as a deforestation pixel to the total area of ​​the geofenced region is greater than a preset area ratio threshold, the deforestation index data is determined to be deforestation. For example, the preset area ratio threshold is 5%.

[0135] If the ratio of the area marked as a deforestation pixel to the total area of ​​the geofence region is less than or equal to a preset area ratio threshold, the forest deforestation index data is determined to be non-deforested.

[0136] In summary, time series analysis can filter out seasonal and temporary interference, thereby improving the accuracy of logging verification results.

[0137] In this example, when the sustainability quantification indicator is the history of deforestation, the calculation of the sustainability quantification indicator data for the geofenced area can also specifically include:

[0138] First, the normalized difference in vegetation index (NDVI) of each pixel in satellite imagery of the geofenced area at a preset historical time and in the current time satellite imagery are calculated to obtain NDVI distribution maps for the two images respectively. Different satellite imagery can be from the same month; for example, if the current time is January of year A, January of any of the past 5 years can be selected.

[0139] Then, change detection and analysis are performed on the NDVI distribution maps corresponding to the two images to obtain NDVI change maps.

[0140] Then, the change magnitude of each pixel in the NDVI change map is compared with a preset pixel change threshold, and pixels exceeding the preset pixel change threshold are marked. The ratio of the geographic area corresponding to the pixels exceeding the preset pixel change threshold to the total area of ​​the geographic fence area is calculated to obtain the vegetation degradation area ratio.

[0141] If the proportion of vegetation degradation area is greater than the preset area proportion threshold (e.g., 5%), the forest deforestation index data is determined to be deforestation; otherwise, it is determined to be no deforestation.

[0142] In this example, when the sustainability quantification indicator is social risk, the sustainability quantification indicator data for the geofenced area is calculated, specifically including:

[0143] Step 1: Using kernel density estimation, perform smoothing analysis on the social risk data of the smallest administrative unit where the geofenced area is located to generate a surface distribution map of social risk density.

[0144] Social risk data can include population and labor data, news and public opinion data, and historical risk event data.

[0145] For example, data on population density, poverty rate, and unemployment rate for the smallest administrative unit can be obtained from publicly available government data. For instance, the World Bank's sub-national poverty database provides regional poverty rate data; International Labour Organization (ILO) statistics provide regional employment rates and informal employment rates; and the World Population Dataset provides high-resolution raster data on population density distribution.

[0146] Historical risk event information for the smallest administrative unit (e.g., a rally of agricultural workers demanding higher wages, previously reported in the media) is obtained from the government's publicly available risk event database, and data with coordinates is exported. Verified labor, environmental, or community conflict events and their geographical locations are extracted from publicly available reports.

[0147] Through the API interface, local news information from news websites and social media for the smallest administrative unit over a preset period of time (e.g., 6 months) can be retrieved. For example, keyword lists (such as strikes, conflicts, pollution protests, land disputes) and geographic bounding boxes can be configured for periodic retrieval.

[0148] It's important to note that each news item or warning message includes geographic coordinates (e.g., the location of the report, social media location), allowing us to identify spatial clustering patterns of negative public opinion. Areas with higher density of negative public opinion are at higher risk of social instability. Historical risk event data is recorded as points with location information, thus identifying hotspots where social risk events have repeatedly occurred in the past; these areas have a higher probability of experiencing risk events again in the future.

[0149] Based on this, all acquired event points (news, history) are uniformly converted into point vector layers with latitude and longitude coordinates. Attributes can be assigned to each point, such as event type, occurrence time, and severity (e.g., based on ACLED event type or GDELT Goldstein score).

[0150] Kernel density estimation is a method that transforms discrete points into a continuous, smooth surface density map. It calculates the influence density of the surrounding region centered on each event point. The closer to the event point, the greater the influence and the higher the density value.

[0151] Based on this, the specific process of smoothing analysis includes the following steps:

[0152] Step 11: Standardize the population and labor data, and then perform a weighted summation of the standardized population and labor data to obtain a composite index.

[0153] Specifically, firstly, the acquired statistical data of administrative units are interpolated into a spatial raster layer with a resolution consistent with the target analysis area using methods such as area weighting or inverse distance weighting. Then, the raster layer of each indicator (such as poverty rate and unemployment rate) is subjected to min-max standardization or Z-score standardization to eliminate dimensions.

[0154] It should be noted that the weight of each indicator in each human and labor data point can be predetermined using methods such as expert scoring (AHP) or principal component analysis. Specifically, for example, a judgment matrix can be constructed using expert questionnaires to calculate the weights of poverty rate (0.4), unemployment rate (0.35), and population rate (0.25), and a consistency test (CR < 0.1) can also be performed.

[0155] Based on this, the standardized indicator grids are weighted and summed at the pixel level to generate a raster map of the population and labor force composite index.

[0156] Step 12: Use the comprehensive index as the weight of the administrative unit center point, and use the administrative unit center point as the center of the kernel function to perform weighted kernel density estimation to obtain the population and labor risk density surface map.

[0157] Bandwidth is a core parameter in kernel density estimation, controlling the smoothness. Those skilled in the art can determine it empirically: a geographic distance is preset based on the analytical scale and data characteristics. For example, 50 kilometers can be used for provincial / state-level analysis; 10-20 kilometers can be used for city / county-level analysis. This is based on empirical judgment of the scope of social risk impact.

[0158] Specifically, the center point of the administrative unit is used as the center point of the kernel function. The weight of this point is the comprehensive index value of that administrative unit. Using the kernel density analysis tool in GIS, the search radius (bandwidth) is set, and the weight field is selected as the comprehensive index value field.

[0159] Step 13: Perform weighted kernel density estimation with the location of each negative news event as the center of the kernel function to obtain the public opinion risk density surface map.

[0160] Specifically, each negative news event is used as the kernel center. The weight field can be selected as a severity score; otherwise, the default weight is 1. Kernel density analysis is performed using the same bandwidth.

[0161] It should be noted that all points are typically assigned the same weight (or different weights may be assigned based on the severity of the news). The kernel function represents that the closer one is to the location of a negative news event, the higher the probability of social instability caused by it. The public opinion risk density surface map reflects the outbreak points and clusters of social tensions in the near future, and has a stronger timeliness.

[0162] Step 14: Perform weighted kernel density estimation with the location of each historical risk event as the center of the kernel function to obtain the historical risk density surface map.

[0163] Specifically, each historical risk event point serves as the kernel center. The weight field can be selected as a severity score; otherwise, the default weight is 1. Kernel density analysis is performed using the same bandwidth.

[0164] It should be noted that the uniform smoothing process also relies on the distance decay principle to identify areas where problems repeatedly occur. The historical risk density surface map reveals long-standing, structural areas of social conflict.

[0165] Step 15: Standardize the pixel values ​​of the population and labor risk density surface map, the public opinion risk density surface map, and the historical risk density surface map to a preset range (e.g., between 0 and 1) to eliminate the influence of dimensions and facilitate comparison and integration.

[0166] Specifically, each kernel density analysis generates a continuous floating-point raster layer, where the cell value represents the risk density at that location. The common standardization method is min-max standardization, which can be performed using a GIS raster calculator.

[0167] Step 16: Weighted overlay the standardized population and labor risk density surface map, public opinion risk density surface map, and historical risk density surface map to generate a social risk density surface distribution map.

[0168] For example, the Analytic Hierarchy Process (AHP) can be used to construct an importance comparison matrix of population and labor force, news and public opinion, and historical events. For example, if the long-term structural impact of historical events is considered the most important, followed by recent public opinion, and then basic population background, weights can be assigned such as [0.3, 0.4, 0.3].

[0169] It should be noted that either two values ​​can be weighted and superimposed, or either value can be used as the surface distribution map of social risk density.

[0170] For example, taking a palm oil plantation in City A as an example, City A has a poverty rate of 8.2%, a youth unemployment rate of 22.1%, and a population density of 1200 people / km². Min-max standardization is performed on each indicator, mapping the values ​​to between 0 and 1: assuming the poverty rate ranges from 5% to 25%, the poverty rate standardization is 0.16. The unemployment rate ranges from 15% to 35%, so the unemployment rate standardization is 0.355. The population density ranges from 50 to 1500 people / km², so the population density standardization is 0.793. The weights of the three indicators are determined using the expert AHP method: poverty rate weight 0.45, unemployment rate weight 0.45, and population density weight 0.10. The comprehensive index of the administrative unit is 0.311. The geometric center point coordinates of the administrative unit are determined: longitude 109.342°, latitude -0.026°. The kernel density analysis tool of GIS software (such as ArcGIS or QGIS) is used: Input data: a vector file containing the above point coordinates and weight values. Weight field: Select the field containing the weight values. Search radius (bandwidth): Set to 20 kilometers (based on the provincial analysis scale), Output cell size: Set to 1000 meters (1 square kilometer), Output: Generate a continuous risk density raster map with a value range of approximately 0.02-0.65.

[0171] Negative news related to palm oil cultivation over the past six months was retrieved via a news API, including: palm oil worker protests against low wages (location: longitude 109.521°, latitude -0.482°), land disputes between communities and plantations (location: longitude 109.612°, latitude -0.415°), and protests sparked by plantation environmental pollution (location: longitude 109.488°, latitude -0.551°). A total of 35 relevant negative news events were collected. Each news event was assigned a weight (all weights were 1 if no severity information was available). Kernel density analysis was performed using GIS software (such as ArcGIS or QGIS): Input data: Vector file containing the coordinates of the 35 news events. Weight field: Weight field (all 1), search radius: 20 km, Output: Generated public opinion risk density raster map with values ​​ranging from approximately 0.01 to 0.45.

[0172] Social risk events from the past 5 years were retrieved from a historical event database: land conflict events (15 points), labor rights events (20 points), and community protest events (25 points), totaling 60 historical risk event points, each with specific geographic coordinates. Different weights were assigned to different event types (reflecting event severity): labor rights events: weight 1.0, community protest events: weight 1.2, land conflict events: weight 1.5. The same kernel density analysis parameters were used.

[0173] Search radius: 20 km, output cell size: 1000 m, output: generated historical risk density raster map, value range approximately 0.03-0.70.

[0174] The three risk density maps were subjected to min-max normalization, mapping the original density values ​​to a range of 0-1 based on their respective value ranges. After normalization, the value range of all three maps is 0-1, with a value of 1 indicating that the risk at that location in that dimension reaches the highest level in the region.

[0175] The relative importance weights of the three dimensions were determined using the Expert AHP method: Historical risk (reflecting long-term structural contradictions): most important, weight 0.55; Public opinion risk (reflecting current dynamics): second most important, weight 0.30; Population and labor risk (reflecting background conditions): relatively minor, weight 0.15. Weighted overlay was performed using a raster calculator in GIS, generating a comprehensive social risk density surface distribution map with values ​​ranging from 0 to 1. Specifically: high value area (0.7-1.0): extremely high social risk, potentially indicating long-standing unresolved social contradictions; median value area (0.3-0.7): moderate social risk, requiring attention to potential risks; low value area (0-0.3): relatively low social risk, indicating relative social stability.

[0176] Step 2: Perform spatial correlation analysis on the surface distribution map of social risk density to identify statistically significant hotspot areas and generate a social risk heat map.

[0177] It should be noted that the risk surface map generated by kernel density estimation shows the risk concentration, but it is unclear whether it is statistically significant or merely a random distribution. That is, spatial correlation analysis is needed to determine whether high-risk areas appear randomly or form meaningful clusters, and whether the clustering patterns are statistically significant. Spatial correlation analysis can be performed using the Moran index or the Getis-Ord Gi statistic.

[0178] When using the Getis-Ord Gi statistic, pixels are classified as follows: Hotspots: high Z-score and small P-value (statistically significant); Colds: low Z-score and small P-value (statistically significant); Not significant. Significance can be determined based on the Z-score. For example, a significant hotspot is defined as Z > 1.96 and P < 0.05, indicating that the risk density of this area and its surroundings is significantly higher than the global average, forming a statistically significant high-risk cluster. A significant coldspot is defined as Z < -1.96 and P < 0.05, indicating that the risk density of this area and its surroundings is significantly lower than the global average, forming a statistically significant low-risk cluster. If these conditions are not met, the pixel is not significant.

[0179] This analysis can be performed using ArcGIS's hotspot analysis tool (Getis-Ord Gi*), GeoDa, or Python's PySAL library (esda.G_Local function), and will output a result layer with GiZScore and GiPValue fields.

[0180] For example, prominent hot spots are represented in red, cold spots in blue, and insignificant areas in gray or no color, thus generating a social risk heat map.

[0181] Step 3: Analyze the geofenced area based on the social risk heat map to obtain the social risk level of the production site.

[0182] Among them, the hot spots, cold spots or inconspicuous areas where the geofenced area is located correspond to different social risk compliance levels, and the levels increase progressively.

[0183] For example, prominent hotspots indicate a high-risk cluster, while prominent cold spots indicate a low-risk cluster. Insignificant hotspots indicate a risk-free cluster. If more than 30% of the geofenced area falls within a prominent hotspot, it is classified as high-risk. If the area within a prominent hotspot is greater than 10% and less than or equal to 30%, it is classified as medium-risk; the remainder is classified as low-risk.

[0184] In this example, performing spatial correlation analysis on the surface distribution map of social risk density may include the following steps:

[0185] The surface distribution map of social risk density is converted into grid data. Each cell has a risk density value, meaning that the value of each cell represents the risk density at that location.

[0186] It should be noted that a continuous risk density surface map is generated by kernel density estimation, representing the risk density at each geographic location. However, computer data is typically processed discretely. Therefore, the continuous surface needs to be discretized into grid data, that is, divided into regularly arranged cells (pixels). Each cell represents a grid area and is assigned a risk density value. This value typically represents the risk density at the center of the cell, or the average risk density within the cell area.

[0187] Then, based on the preset grid distance or adjacency relationship, the spatial weight matrix corresponding to the grid data is determined.

[0188] When based on distance, units that are closer together should have a greater mutual influence, while units that are farther apart should have a smaller influence. The weight is the reciprocal of the distance. When based on adjacency, units that are adjacent vertically, horizontally, or vertically have a weight of 1, otherwise a weight of 0.

[0189] For example, after grid transformation of the risk density surface map, there are 4 grid cells. The risk value of grid cell 1 is 0.92, the risk value of grid cell 2 is 0.65, the risk value of grid cell 3 is 0.78, and the risk value of grid cell 4 is 0.18. A distance threshold of 20km is set, and the distance between each grid cell and other grid cells, as well as the weight of each grid cell relative to its surrounding grid cells, are calculated, as shown in Table 1.

[0190] Table 1 Grid Weights

[0191]

[0192] It should be noted that grid 1 and grid 2 are very close (5km), with a weight of 0.8 (strong influence); grid 1 and grid 3 are moderately close (15km), with a weight of 0.3 (moderate influence); and grid 1 and grid 4 are very far apart (50km), with a weight of 0 (no influence).

[0193] Then, based on the spatial weight matrix, the Getis-Ord Gi statistic for each grid is calculated to obtain the Z score and P value for each grid.

[0194] The statistical expression is as follows:

[0195]

[0196] in, It is a statistic. It is the risk density value at pixel location j (grid location). It is the spatial weight between pixel position i and pixel position j, where n is the total number of pixels. S is the average risk density, and S is the standard deviation of the risk density.

[0197] Finally, multiple grids were classified based on Z-scores and P-values ​​to identify significant hot and cold areas.

[0198] In one example, an automated decision-making logic tightly coupled with clear priorities and linked to geographic location was constructed to determine activity data outside of transportation processes. This logic may include the following steps:

[0199] Step 1: Obtain direct monitoring data of the geographical location, and when the direct monitoring data meets the first preset valid conditions, obtain activity data with a first credibility level identifier.

[0200] In this example, direct monitoring data for geofenced areas specifically refers to physical quantity data that is automatically and transmitted in real time through calibrated IoT sensors or smart meters (such as mass flow meters, smart meters, and weighbridges) installed on production lines, tanks, pipelines, or critical equipment. This type of data is considered the highest quality data source due to its automated collection, real-time nature, and minimal human intervention.

[0201] In this case, the activity data can be acquired as follows: Smart sensors on the production line (such as mass flow meters and smart meters) transmit real-time readings to the enterprise edge gateway via industrial IoT protocols (such as MQTT and OPC UA). The edge gateway performs preliminary range verification on the data (such as whether the instantaneous flow rate is between 5% and 95% of the device's capacity), and then encrypts and transmits the data to the authentication service platform.

[0202] The first preset validity conditions may include the following: Data validity: The value is within a reasonable physical range (e.g., the flow rate is not negative). Temporal continuity: There are no large gaps or abnormal interruptions in the data within the calculation period. Equipment status identification: The sensor or instrument from which the data comes reports a normal status and is not in calibration, fault, or maintenance period. Communication integrity: The data upload link is reliable, and the timestamps are complete and continuous. For example, the numerical validity is 100% (no abnormal values ​​exceeding the physical range, such as the flow rate not being negative), the temporal integrity is ≥98%, the equipment status identification is continuously normal operation, and the timestamps are continuous and completely match the calculation period.

[0203] Based on this, when the activity data meets the above conditions, the system will automatically label the activity data as the first confidence level.

[0204] Step 2: If the first preset valid condition is not met, obtain the operational business data of the geographical location.

[0205] In this example, when the direct monitoring data verification finds that the first preset valid condition is not met (such as the discovery of a long-term data gap or instrument fault indicators), the system automatically switches to the next data channel.

[0206] Operational business data refers to digitized data recorded from a company's formal operational management system. Examples include production work orders, material input records, and production reports from Enterprise Resource Planning (ERP) or Manufacturing Execution System (MES); material quantity information on purchase invoices, sales invoices, and logistics waybills that have been signed and confirmed by both parties; and consumption readings on monthly water, electricity, and gas bills provided by utility companies.

[0207] Step 3: When the operational business data meets the second preset validity conditions, activity data with a second credibility level identifier is obtained.

[0208] At this point, activity data can be acquired by periodically retrieving key business tables from the ERP or MES system, such as production work orders, material movement vouchers, and daily energy consumption reports, through secure APIs provided by the enterprise (such as RESTful APIs) or encrypted file transfers (such as SFTP). For external documents (such as invoices and waybills), OCR (Optical Character Recognition) technology is used to identify key fields (material name, quantity, date) and associate them with structured data.

[0209] The second preset validity condition may include the following: Multi-source consistency verification: For example, logically balancing the raw material purchase quantity on the purchase invoice, the material input quantity on the production work order, and the finished product warehousing quantity. If they can form a closed balance within a reasonable error range, they are considered consistent. Due to the possibility of data entry errors or human adjustments, the system introduces cross-validation logic. Time range matching: Ensure that the period of the ledger records matches the product batch period to be calculated. For example, if the multi-source consistency error is ≤5%, the data is considered logically consistent (such as the logical balance error of the purchase quantity, material input quantity, and output). Time range matching degree 100% (the business record period is completely consistent with the product batch period). The document or record has a verifiable unique identifier (such as invoice number, work order number) and is not marked as abnormal.

[0210] It should be noted that the formula for calculating the balance error is as follows:

[0211]

[0212] in, To balance the error, Beginning inventory, For this period's investment, For this period's output, Ending inventory.

[0213] Step 4: If the second preset valid condition is not met, acquire geographic image data and geographic feature data related to the production space within the geographic fence area.

[0214] In this example, geographic feature data within the geofence is retrieved from a map service platform (e.g., via API call). Satellite remote sensing or aerial imagery is automatically retrieved to obtain recent (e.g., within 30 days prior to the current accounting date) high-resolution images (e.g., images with less than 20% cloud cover) to obtain geographic imagery data.

[0215] Geographic feature data includes POIs, the number, size, or density of specific types of facilities. For example, for factories, this might include the identification and counting of chimneys, storage tanks, and factory building outlines. For farms, it might include the number and distribution of greenhouses and irrigation facilities. For collection points, it might include the number and size of catering businesses within their service area.

[0216] For example, POI data is obtained by calling commercial or open-source map service APIs (such as the map's surrounding search interface or OpenStreetMap's Overpass API), using the geometric center or boundary of the geofenced area as the query range and a preset industry classification code as the keyword (e.g., for waste oil collection points, the keywords are Chinese restaurants, fast food restaurants, hotel restaurants, etc. under the catering service category). This allows for the batch acquisition of POI coordinates, names, types, and size levels (if any). The API typically returns data in JSON format. The acquired POI data is then deduplicated (based on coordinates and names), corrected (by converting map coordinates to a unified WGS84 coordinate system), and classified statistically analyzed. For example, the number of large (revenue > X million yuan), medium, and small catering businesses within the geofenced area can be counted, and their density (number per square kilometer) can be calculated.

[0217] Step 5: Assess the production activity of the product within the geofenced area based on the geographic imagery data.

[0218] Geographic imagery data is used to assist in assessing regional productivity. By analyzing geographic imagery (such as nighttime light data, traffic density, and new building identification), the system uses image recognition or change detection models to dynamically assess the productivity of the region. For example, an area with frequent recent construction activity and bright nighttime lights would be rated as having high productivity.

[0219] In this example, the evaluation method may include the following steps:

[0220] First, based on the boundary coordinates of the geofenced area, geographic imagery data covering a predetermined historical period within the geofenced area is acquired. This geographic imagery data includes at least two of the following: optical imagery, nighttime light imagery, and synthetic aperture radar (SAR) imagery. To ensure the comprehensiveness and robustness of the assessment, the system avoids reliance on a single data source and instead employs a multi-source remote sensing data fusion strategy. Different types of data reflect activities from different physical dimensions (visible light, nighttime light, microwave), and their complementarity overcomes the limitations of single data sources (such as cloud cover obstructing optical imagery and limited resolution of nighttime light).

[0221] It should be noted that the preset historical time period can be set according to actual needs, such as the most recent six months or one year. Setting at least two types of image data ensures minimal data fusion, thereby enabling assessments to be initiated under various weather and data availability conditions.

[0222] The geographic imagery data is then processed to generate one or more feature layers to characterize regional dynamic changes. These feature layers may include at least one of the following: a construction land change layer reflecting changes in production facilities, a nighttime light intensity change layer reflecting the intensity of economic activity, and a surface activity index layer reflecting surface disturbance.

[0223] Among these methods, land cover classification algorithms (such as semantic segmentation models based on deep learning) are applied to temporal optical imagery to identify land types such as industrial production facilities, warehousing land, and construction sites. By comparing the classification results of different time phases, the system can automatically detect changes in facilities, such as new construction, expansion, and demolition, forming change detection patches and generating vector or raster change layers.

[0224] For time-series nighttime light imagery, calculate statistical values ​​(such as average light intensity and total light volume) within the geofenced area. By analyzing the time series of these statistical values ​​(e.g., creating monthly variation curves), features such as brightness growth trends and seasonal fluctuations can be extracted to form a layer representing the dynamics of economic activity.

[0225] For time-series synthetic aperture radar (SAR) images, techniques such as interferometric processing or temporal coherence analysis are employed. For example, continuous material stacking, equipment movement, or ground compaction can cause regular changes in radar echo signals. Based on this, the system can identify active work areas and generate index layers reflecting the intensity and extent of surface disturbance.

[0226] Then, based on one or more feature layers, one or more corresponding activity metrics are calculated.

[0227] Specifically, for the construction land change layer, the percentage of newly added facilities area (new area / total area) or the monthly average number of change patches can be calculated. For the nighttime light intensity change layer, the annual average growth rate of light intensity or the comparison index between recent light brightness and the same period in history can be calculated. For the surface activity index layer, the proportion of highly disturbed areas or the average surface deformation rate can be calculated. Based on this, these indicators are all quantified into numerical values, forming a feature vector that describes the regional activity from different perspectives.

[0228] Finally, by integrating one or more activity metrics, the overall production activity coefficient of the geofenced area is determined.

[0229] This process involves normalizing various activity metrics of different dimensions to a uniform range (e.g., 0-1) using methods such as max-min scaling. Then, a weighted linear combination or a simple machine learning model (e.g., training a classifier) ​​is used to merge the normalized metrics into a comprehensive production activity coefficient. The weights can be set based on expert knowledge or determined through historical data regression analysis. It should be noted that the final output can be a continuous coefficient (e.g., 0.85, indicating very high activity).

[0230] For example, weights can be determined using regression analysis based on historical data: a sample of geofenced areas with known levels of actual production activity (such as historical output or energy consumption data obtained through direct monitoring) is collected. For each sample area, its various activity metrics are calculated ( Intensity of construction land expansion Nighttime light growth rate (Ratio of surface disturbance area). Using the actual level of production activity as the dependent variable Y, and [...] , , Using as the independent variable, construct a multiple linear regression model: The regression coefficients were obtained by using the least squares method. , , After normalization (making the sum of the three equal to 1), they can be used as weights in the calculation of the comprehensive production activity coefficient.

[0231] In summary, unstructured geospatial data such as multi-source, multi-temporal remote sensing images can be transformed into structured indicators that characterize the intensity of regional production activities through a series of automated processing and analysis, thereby providing objective and dynamic input parameters for subsequent material or energy consumption estimation models.

[0232] For example, to assess the production activity of a palm oil plantation (geofence area) in Southeast Asia over the past year to help estimate its fresh fruit bunch harvest, the system automatically invoked remote sensing data services to acquire the following imagery data covering the geofenced area over the most recent 13 months (e.g., June 2023 to June 2024):

[0233] Monthly Sentinel-2 optical imagery (10-meter resolution) is used to observe vegetation and surface changes. Monthly VIIRS nighttime light imagery is used to observe nighttime operational lighting. Sentinel-1 radar imagery (C-band) is provided every 12 days for all-weather monitoring of micro-surface changes.

[0234] The system performs automated preprocessing (correction and registration) on the above images, and then runs the built-in model to generate feature layers, as detailed below:

[0235] Construction land change layer generation: Land classification was performed on each period of Sentinel-2 imagery to identify mature palm groves, newly cleared / sapling forests, roads, processing plant areas, etc. Furthermore, by comparing the classification results from June 2023 and June 2024, the system detected approximately 20 hectares of land in the southwest corner of the region that had transitioned from secondary forest to newly cleared land. This change was generated as a vector patch, forming the facility expansion feature layer.

[0236] Nighttime light intensity variation layer generation: Extracting monthly nighttime light brightness values ​​within a 500-meter radius of the factory area's coordinates. The system plotted brightness curves for 13 months, revealing a stable upward trend in brightness values ​​since November 2023, peaking in May 2024, representing an average increase of approximately 40% compared to the same period in 2023. This trend was quantified into a time-series layer.

[0237] Surface activity index layer generation: Interferometric coherence analysis was performed on time-series Sentinel-1 radar imagery. The analysis revealed that the main road area connecting the plantation hinterland and the processing plant, as well as the raw material storage area within the processing plant, maintained extremely low coherence (<0.2) throughout the year, indicating that the surface in these areas is continuously disturbed by frequent human activities or vehicle traffic. Based on this, the system generated a raster layer representing a high-disturbance area.

[0238] Based on the above layers, the system calculates three core quantitative indicators as follows: Expansion intensity indicator: 2%; Newly reclaimed area (20 hectares) / Total fenced area of ​​the plantation (assuming 1000 hectares). Nighttime operation indicator: Average annual growth rate of nighttime light intensity in the processing plant area is 40%. Logistics disturbance indicator: 5%; Area of ​​high disturbance zone (roads and storage yards, assuming 50 hectares) / Total fenced area of ​​the plantation.

[0239] Based on this, the system normalizes the three indicator values ​​(2%, 40%, 5%) to the range of 0-1 according to the historical data range, assuming the values ​​to be (0.4, 0.8, 0.5). If the weight of nighttime lighting is 0.5, the expansion intensity is 0.3, and the logistics disturbance is 0.2, then the overall production activity level obtained through weighted fusion is 0.62.

[0240] Step 6: Input the production activity and the geographic feature data into the preset estimation model to obtain activity data with a third confidence level identifier.

[0241] In this example, the preset estimation model is as follows:

[0242]

[0243] in, For activity data, To increase production activity, Geographic feature data, This refers to the industry average capacity coefficient or the regional statistical benchmark value.

[0244] For example, to estimate the monthly collection volume of waste oil at a collection point without a weighbridge, the geographical feature data could be the number of catering businesses. This represents the average amount of waste cooking oil generated by businesses in the region.

[0245] In one example, for non-transportation links, the factors that best fit the local characteristics of the specific production location are prioritized. When these factors are unavailable, local characteristic parameters are used to scientifically correct the higher-level factors, and finally, a general factor is used as a safety net. This process ensures that the accounting results achieve the highest possible localization accuracy while guaranteeing reliable execution of the accounting under any circumstances. It includes the following steps:

[0246] Step 1: Match the emission factor type of the geofence area in the localized emission factor database of the administrative region or localized geographical unit to which the production site belongs.

[0247] The matching process is as follows: The precise geographic coordinates of the production site and the required emission factor type are used as query keywords to prioritize searching the localized emission factor database. For example, a spatial query is performed in the localized emission factor database, filtering out the first set of records whose spatial range includes the geographic coordinates. If the first set of records is not empty, a second set of records that exactly matches the required emission factor type is filtered out. If the second set of records is not empty, and there are multiple second records (e.g., matching both industrial park and city-level factors), the most granular record is selected according to a preset spatial granularity priority (e.g., specific facilities > industrial parks > districts / counties > cities) to determine the corresponding emission factor value. If there is only one second record, the emission factor value corresponding to that second record is adopted.

[0248] The localized emission factor database includes: detailed emission inventories issued by specific industrial parks or economic development zones, localized guideline emission inventories issued by the ecological and environmental departments of the directly administered regions, and verified emission factor inventories for specific large-scale facilities or joint ventures.

[0249] In this example, the system automatically determines which specific region, which has been compiled and included in the emissions factor database, belongs to the production site based on its geographical coordinates. Specifically:

[0250] Legally defined administrative regions: such as municipal, county, and district-level administrative regions, whose ecological and environmental departments may publish localized greenhouse gas emission factor lists.

[0251] Localized geographical units refer to any spatial unit with independent emission characteristics, delineated based on a shared energy structure, industrial characteristics, or management policies. Examples include economic and technological development zones, high-tech industrial parks, eco-industrial parks, free trade zones, and key river basin management areas. These units may have specific emission factors that distinguish them from higher-level regions due to the use of centralized heating, biomass energy, or special waste treatment methods.

[0252] It should be noted that the localized emission factor database is a spatial-attribute database, and each record contains at least three key fields: the applicable geographic unit identifier, the emission factor type (such as grid electricity, pipeline natural gas, and industrial steam), the specific value of the factor type, and its unit of measurement.

[0253] It should be noted that this rating factor comes from the most direct local production and sales or actual measurement data, and has the smallest spatial representativeness error.

[0254] Based on this, the errors caused by directly using national or provincial average factors are avoided. By matching localized geographical units, the emission factors can reflect the unique energy structure, industrial policies, and technological level of the region (for example, a park that uses a large amount of wind power has a grid factor that is significantly lower than the national average). This solves the problem of insufficient spatial accuracy.

[0255] Step 2: If the match is successful, obtain the emission factor with the first confidence level identifier.

[0256] That is, the emission factor values ​​obtained from the data collection channels in step 1 are determined as the corresponding first confidence level.

[0257] Step 3: If the matching fails, obtain the characteristic parameters of the administrative region to which the production site belongs, and correct them according to the standard emission factor of the next higher level administrative region to which the production site belongs; the characteristic parameters include at least one of energy structure, industrial composition, and climate zone type.

[0258] In this example, instead of directly reverting to a broader regional average, a localization correction engine is activated. The purpose is to use the characteristic information of the administrative unit to which the production location belongs (such as a city or county) to fine-tune the standard factors of the next higher level of administrative region (such as a province or country) to make them closer to the local reality.

[0259] It should be noted that energy structure parameters include the proportion of local renewable energy generation and the ratio of coal-fired to gas-fired power generation, obtained from local statistical yearbooks or power balance sheets. If a city's clean energy share is significantly higher than the provincial average, its grid emission factor should be adjusted downwards from the provincial factor. Industrial composition parameters include the proportion of added value from high-energy-consuming industries or the type of dominant industries in the region. For example, the industrial heat emission factor of a region dominated by steel may differ from that of a region dominated by light industry. Climate zone type parameters include the building thermal design zone or the number of heating degree-days. This directly affects building energy consumption-related emission factors; heating emissions in colder regions are typically higher.

[0260] In this example, the system includes pre-built standardized correction models for different types of emission factors (such as electricity, heat, and transportation fuels). These models are typically based on historical statistical data and establish a quantitative relationship between characteristic parameters and emission intensity (such as regression models).

[0261] Based on this, the correction process may include the following steps:

[0262] First, quantitative data of at least one characteristic parameter of the administrative region to which the production site belongs is obtained.

[0263] At least one of these refers to the system's ability to be flexibly configured, allowing for the selection of one or more of the most representative parameters for correction based on the sensitivity of different emission factor types to various parameters, in order to balance accuracy and data availability.

[0264] Then, based on the quantified data of the at least one characteristic parameter, a localized correction coefficient for the standard emission factor of the next higher administrative region is calculated using a pre-set correction model that reflects the relationship between the characteristic parameter and emission intensity.

[0265] The core of establishing a quantitative mapping relationship from characteristic parameters to emission intensity adjustment ratios is calculating a correction coefficient. This coefficient reflects the expected deviation of local emission intensity from the average level of the higher-level region.

[0266] Before system deployment, one or more modified models were trained and solidified for each type of emission factor (such as grid power, heat supply, and diesel combustion) through historical big data analysis and machine learning methods.

[0267] For example, a typical simplified correction model for the grid emission factor might look like this:

[0268]

[0269] in, This refers to the emission factor of the municipal power grid. The proportion of renewable energy power generation in the city, It is the average proportion of renewable energy power generation in the province.

[0270] For example, if the city's renewable energy share is 40% and the provincial average is 30%, then the correction factor is 0.857, which means that the city's grid factor should theoretically be 85.7% of the provincial factor.

[0271] Finally, the localized correction factor is combined with the standard emission factor to generate the corrected emission factor.

[0272] Among them, the product of the localization correction factor and the standard emission factor of the next higher level of administrative region to which the production site belongs can be used as the corrected emission factor.

[0273] Step 4: If the correction is successful, the emission factor with the second confidence level identifier is obtained.

[0274] That is, the necessary feature parameters have been obtained, and the correction model is running normally and outputting reasonable coefficients (e.g., 0.7-1.3). If any of the core feature parameters required by the correction model cannot be obtained (e.g., there is neither energy structure data nor industry composition data) or the correction coefficients output by the correction model are outside the reasonable range, then the correction is considered to have failed.

[0275] Step 5: If the correction fails, use the average factor from the general standard database to obtain the emission factor with a third confidence level identifier.

[0276] Explicit sources of common standard databases can include international databases: the default factor in the Intergovernmental Panel on Climate Change (IPCC) National Greenhouse Gas Inventory Guidelines, the European Commission's ELCD Core Database, and the Ecoinvent Database.

[0277] National / Regional Databases: Average Emission Factors of Chinese Provincial Power Grids (released by the Ministry of Ecology and Environment), and the Emission Factors Center (EFDB) of the U.S. Environmental Protection Agency (EPA).

[0278] This system can include a built-in factor mapping table. For example, it can map emission factor types from production locations to query keys in one or more backup databases. The system automatically performs the lookup in the mapping table. The database selection priority can be: national factors > regional (e.g., Asian) average factors > internationally accepted default factors.

[0279] For example, for the power grid factor of a certain city in China, if the first and second levels fail, the latest national average emission factor value of the power grid published by the Ministry of Ecology and Environment in that year will be directly used as the emission factor value of the third confidence level.

[0280] In summary, dynamic coupling between emission factors and geospatial data has been achieved: the selection of emission factors has been transformed from a static table lookup operation into a dynamic decision-making and calculation process triggered by geographic coordinates, enabling carbon accounting results to truly reflect the differences in environmental impacts caused by where production takes place.

[0281] By incorporating local characteristic parameters for correction, the level of granularity in accounting is significantly improved, even in the absence of an existing local factor list. Combined with the quantification of different levels of uncertainty, the confidence interval for the final carbon footprint result can be scientifically calculated, resulting in a highly transparent report.

[0282] It should be noted that for other stages in the product lifecycle (such as raw material transportation, product distribution, usage, and waste disposal / recycling), carbon emission calculations can be based on existing and mature universal LCA accounting models, standard databases, and calculation methods within the industry. This is because the carbon emission intensity of these stages primarily depends on general technical parameters, fuel type, transportation mode, or average processing technology, rather than the unique geographical characteristics of a specific production location. Therefore, using industry-standard data and methods is reasonable, reliable, and efficient. However, other stages can also be calculated by collecting activity data and emission factors from multiple data collection channels with varying levels of credibility, and then combining this with universally accepted carbon emission accounting models. For example, emission factors in the transportation stage (such as carbon emissions per unit distance for trucks, ships, and airplanes) have internationally recognized databases (such as IPE, Ecoinvent, and GREET models); emission factors in the waste disposal stage (such as landfill, incineration, and recycling) have standard values ​​in environmental management databases of various countries. The calculation methods are mature: the accounting models for these links (such as transportation carbon emissions being the product of transportation distance, transportation volume, and transportation mode emission factors) have been widely adopted and verified by international standards such as ISO 14040 series standards and PAS 2050.

[0283] For example, the carbon emission calculation process for the transportation segment can also be as follows: A method combining route optimization and uncertain emission factors is used to calculate carbon emissions for the transportation segment. Based on the starting and ending points of the transportation route, an optimization algorithm (such as a multi-objective genetic algorithm) is used to comprehensively consider constraints such as road type, altitude changes, and real-time traffic conditions to simulate and calculate a route with the lowest carbon emissions. Then, based on data such as the distance and estimated fuel consumption of this route, combined with emission factors for the transportation vehicles that have a credibility level indicator, the carbon emissions for the transportation segment are accurately calculated.

[0284] In one example, for the production stage, after obtaining activity data and emission factors, the emission factors can be compensated by incorporating the localized environmental characteristics of the geofenced area. The carbon emission calculation process for the production stage is as follows:

[0285] First, retrieve the environmental characteristic weights of the production site from the environmental characteristic weight database.

[0286] The environmental characteristic weights can be determined through the following methods: multiple regression analysis, which is based on the statistical relationship between a large amount of field-measured carbon emission data and geographical parameters; or by combining the experience and judgment of LCA and ecological experts. This approach simultaneously considers three key dimensions: carbon, water, and ecology, thus achieving a refined integration across multiple dimensions.

[0287] Then, based on the environmental characteristic weights, the environmental characteristic data are weighted and calculated to obtain the regional environmental characteristic coefficient of the production site.

[0288] In this example, to reflect the potential impact of environmental differences in different geographical regions on carbon emissions, this application introduces a regional environmental characteristic coefficient (denoted as ). This coefficient is used to localize emission factors. It is calculated by integrating key regional ecological and environmental parameters, and its expression is as follows:

[0289]

[0290] in, For soil carbon storage, As the characteristic weight of soil carbon storage, Due to water resource pressure, As the weight of water resource pressure characteristics, The vegetation cover index, The feature weights of the vegetation cover index are denoted as .

[0291] It should be noted that high soil carbon storage represents significant carbon sink potential, but once damaged (e.g., by land-use change), it can also become a source of high carbon emission risk. Therefore, this item is presented as a positive term to quantify this potential risk. Greater water resource pressure typically means higher energy inputs required for agricultural activities such as irrigation, indirectly leading to increased carbon emissions. Therefore, this item is presented as a positive term to reflect the indirect emission costs resulting from water scarcity. Healthy vegetation possesses strong carbon sequestration capabilities and ecological regulation services, reducing the net environmental impact of the region. Therefore, this item is presented as a negative term to reflect the carbon sink service value of the ecosystem.

[0292] Based on this, this application constructs a dynamically adjusted emission calculation model by introducing regional environmental characteristic coefficients. This model enables the final product carbon footprint accounting results to reflect the specific geographical and environmental characteristics of the production location, thereby achieving localized dynamic adaptability and higher spatial accuracy of the accounting results.

[0293] Finally, based on the emission factors of the activity data and the regional environmental characteristic coefficients, carbon calculations are performed on the activity data to obtain the carbon emissions of the product during the production process.

[0294] In this example, the calculation method is based on the coefficient adjustment method, and the expression is as follows:

[0295]

[0296] in, express Type of activity data, express The emission factor for the type of activity data, where G represents carbon emissions.

[0297] The above expression treats the entire production process as a whole and adjusts it uniformly. That is, it calculates the standard total carbon emissions based on general data and multiplies this standard total carbon emissions by a regional environmental characteristic coefficient. This regional environmental characteristic coefficient acts as a global regulator. If the regional environmental characteristic coefficient is greater than 1, it indicates that the environmental characteristics of the production site provide positive ecosystem services, and the total emissions are reduced. If the regional environmental characteristic coefficient is less than 1, it indicates that the local environment provides positive ecosystem services, and the total emissions are reduced.

[0298] In summary, the calculation is efficient and the concept is intuitive, summarizing the complex effects of geographical heterogeneity into a single overall adjustment factor.

[0299] In one example, uncertainty analysis can be performed again on carbon emissions from the production process based on the characteristics of geographical parameters.

[0300] Because there are inherent spatial correlations among geographical parameters—for example, areas with high soil carbon storage often also have good vegetation cover—independent random sampling may produce unrealistic parameter combinations such as high soil carbon storage and extremely low vegetation cover, which are rarely seen in reality, thus distorting uncertain results. To address this issue, the Copula function is introduced. Based on this, the parameter correlation structure is embedded into Latin hypercube sampling using the Copula function. Specifically, samples can first be generated in the correlation space and then mapped back to the actual distribution of each parameter.

[0301] It should be noted that a spatial correlation covariance matrix (the spatial correlation structure between parameters) is established in advance through historical data analysis.

[0302] This correlation structure effectively quantifies the pairwise correlation strength between various geographic feature parameters (such as soil carbon, water resource pressure, and NDVI). Then, during sampling, the Copula function ensures that the generated parameter sample set not only conforms to the probability distribution (marginal distribution) of each parameter but also maintains the true spatial correlation structure between them. This guarantees that the parameter combinations used in each Monte Carlo simulation are geographically consistent and realistic.

[0303] Drawing a large number of samples from the aforementioned correlated multidimensional probability distributions would be inefficient through simple random sampling, potentially requiring an extremely large sample size to cover all possible parameter combinations. Therefore, Latin hypercube sampling, a stratified random sampling strategy, is employed.

[0304] In this example, the process of obtaining the target sampling point set includes the following steps:

[0305] First, determine the marginal probability distribution and spatial correlation structure between the multiple geographic feature parameters to be sampled.

[0306] Then, based on the Latin hypercube sampling strategy, an initial set of sampling points is generated within the unit hypercube to ensure uniform projection distribution across each parameter dimension.

[0307] The process involves uniformly dividing the value range of each parameter into N equally probable intervals, where N is the number of sampling attempts. A value is randomly selected from each interval to form the sampling sequence for each parameter. Finally, the sampling sequences of each parameter are randomly arranged and combined to form the initial set of sampling points.

[0308] For example, for the soil carbon storage parameter, the probability distribution from 0% to 100% is divided into 1000 intervals (if N=1000). In the soil carbon storage dimension, 1000 sample points will evenly cover all possible ranges from the lowest to the highest value, with no blank areas. This ensures a uniform projected distribution.

[0309] Then, based on the Copula function, the initial set of sampling points is transformed into an initial set of associated sampling points with a target correlation structure.

[0310] For example, the Copula function is Gaussian Copula, which transforms the initial sample point set into the correlated normal space through the quantile function of the standard normal distribution, decomposes the spatial covariance matrix using Cholesky decomposition, applies a correlation structure, and transforms the samples in the correlated normal space back into the unit hypercube space through the cumulative distribution function of the standard normal distribution, thus obtaining the initial correlated sample point set.

[0311] Then, based on the iterative optimization algorithm, the initial set of associated sampling points is optimized to obtain the final set of associated sampling points.

[0312] Then, each point in the associated sampling point set is mapped to the actual parameter space through the inverse cumulative distribution function corresponding to the marginal probability distribution of the corresponding parameters, thus obtaining the target sampling point set.

[0313] Finally, Monte Carlo simulations were performed on the target sampling point set to calculate carbon emissions and obtain the confidence interval for carbon emissions in the production process.

[0314] In summary, this method can efficiently and uniformly explore the entire high-dimensional parameter space with far fewer samples than pure random sampling, avoiding sampling blank areas, thus significantly improving the efficiency of Monte Carlo simulation.

[0315] In this example, to further improve the quality of the Latin hypercube sampled samples, an iterative optimization algorithm (such as exchange-based stochastic optimization) can be used to process the associated sample point set initially generated by the Copula function.

[0316] Optimization Objective: The optimization objective is to minimize the difference between the empirical correlation matrix of the sampled point set and the target correlation matrix defined by the Copula function. Commonly used measures of difference include the Frobenius norm of the matrix or the sum of squared errors of the Pearson correlation coefficient. The specific optimization steps are as follows:

[0317] 1. Calculate the current difference: Calculate the empirical correlation matrix of the initial associated sampling point set and evaluate the difference between the empirical correlation matrix and the target correlation matrix.

[0318] It's important to note that the empirical correlation matrix calculates the correlation between multiple parameters represented by all points in the current initial correlation sampling set. It's not a matrix of a single point or coordinate, but rather a statistical measure calculated based on data from all points in the entire set, describing the degree of correlation between these parameters under the current sampling results. The target correlation matrix is ​​a pre-defined, desired ideal correlation structure (e.g., a positive correlation between soil carbon storage and vegetation index). The difference values ​​compare the gap between the correlation presented by the current sampling set and the target correlation.

[0319] For example, having Several parameters (e.g., soil carbon storage, water resource pressure, vegetation cover index) were used to obtain a current associated sampling point set containing n sampling points. This point set can be represented as a... The matrix X, Indicates the first The sampling point of the nth sampling point The values ​​of each parameter are determined. Based on this, the process of calculating the empirical correlation matrix is ​​as follows:

[0320] The objective is to calculate a A symmetric matrix R, whose elements Indicates the first The empirical (sample) Pearson correlation coefficient between the qth parameter and the qth parameter.

[0321] For any two parameters p and q (i.e., the p-th and q-th columns of matrix X), the correlation coefficient

[0322] The calculation formula is:

[0323]

[0324] in, and These are the values ​​of the i-th sampling point in terms of parameters p and q, respectively. and are the average values ​​of all sampling points for parameters p and q, respectively. The numerator is the covariance of the two parameters. The denominator is the product of the standard deviations of the two parameters, used to standardize the correlation coefficient to the interval [-1, +1].

[0325] Based on this, all empirical correlation coefficients are represented as an empirical correlation matrix. It should be noted that the values ​​on the diagonal of the matrix are always 1, because the parameters are perfectly positively correlated with themselves.

[0326] 2. Generate new samples: Generate a new set of associated sampling points by randomly swapping the coordinates of the initial set of associated sampling points on a certain parameter dimension.

[0327] For example, two sampling points (e.g., point A and point B) are randomly selected, and their coordinates on a specific parameter (e.g., soil carbon storage) are swapped. The values ​​of other parameters (e.g., water resource pressure, vegetation index) remain unchanged. This operation maintains the uniformity of the marginal distribution of each parameter, slightly perturbs the joint distribution structure between parameters, and thus attempts to make the empirical correlation matrix of the point set closer to the target matrix.

[0328] 3. Evaluation and Replacement: Calculate the difference between the empirical correlation matrix and the target correlation matrix of the new correlation sampling point set. If the new difference is less than the current difference, then the new correlation sampling point set replaces the current initial correlation sampling point set.

[0329] Iteration: Repeat steps 2 and 3 until the convergence condition is met (e.g., the difference value no longer decreases significantly) or the maximum number of iterations is reached. The resulting set of sample points is the final set of associated sample points.

[0330] In summary, through iterative optimization, while retaining the advantages of Latin hypercube sampling stratification, it is possible to ensure that the joint distribution of the sample point set in the multidimensional parameter space more accurately reflects the true geographic correlation pattern. This compensates for the potential approximation limitations of the Copula function in high-dimensional cases, generating statistically more consistent and realistic parameter samples. This provides higher-quality input for subsequent Monte Carlo simulations and enhances the credibility of the uncertainty quantification results.

[0331] In one example, verifying sustainability and carbon emissions compliance with sustainability standards includes the following steps:

[0332] First, when the upper limit of the confidence interval for life-cycle carbon emissions is less than a preset emission threshold, climate prediction data for the geofenced area over a future timeframe is extracted from global climate change prediction models. For example, if the climate change prediction model is CMIP6 with a future timeframe of 30 years, the future climate prediction data includes forecasts for indicators such as precipitation and temperature. It should be noted that climate risk analysis refers to assessing the potential physical and climatic risks that the geofenced area may face in the future, primarily extreme weather events such as floods, droughts, extreme heat, and wildfires. This is crucial for the long-term stability of the supply chain.

[0333] For example, the system will call the prediction data of indicators such as precipitation and temperature in CMIP6, and use specialized disaster models (such as hydrological models and drought index models) to calculate the probability of severe flooding or drought in the region in the next 30 years.

[0334] Then, based on climate prediction data, the climate risk level of the geofenced area is predicted over a future period.

[0335] This involves analyzing future climate forecast data to predict the drought and flood risk levels of the geographically fenced area over a future period. For example, this analysis can be based on hydrological models and drought index models. Then, based on the drought and flood risk levels, the climate risk level within the geographically fenced area is determined.

[0336] It should be noted that the higher the drought risk level or the higher the flood risk level, the higher the climate risk level.

[0337] Then, based on the climate risk level, the preset emission threshold is compensated. The climate risk level is negatively correlated with the compensation coefficient.

[0338] Specifically, the compensation mapping table retrieves the compensation coefficient for each climate risk level, and then multiplies this coefficient by a preset emission threshold. The higher the climate risk level, the smaller the compensation coefficient; the compensation coefficient is greater than 0 and less than 1. This is equivalent to lowering the emission threshold.

[0339] Finally, if there is no history of deforestation, carbon emissions are less than the preset compensation emission threshold, and the social risk level is lower than the preset level, the product's ISCC compliance verification is deemed successful.

[0340] It should be noted that if a product has a history of deforestation, carbon emissions are not less than the preset compensation emission threshold, or the social risk level is not lower than the preset level, the product's compliance with sustainability standards will fail the verification.

[0341] Based on this, a history of no deforestation and compliance with lifecycle carbon emission standards are established as mandatory environmental and carbon performance thresholds. This ensures that certified products meet the most basic sustainable development red lines. Subsequently, a social risk level assessment is introduced, incorporating social dimensions into the considerations to prevent green production at the expense of social interests. This sequential screening mechanism, filtering at each level, ensures that the supply chain entering the final certification stage is reliable in the three basic dimensions of environment, carbon, and society.

[0342] By triggering the climate risk prediction module, the long-term physical resilience of the supply chain is assessed by analyzing the probability of future extreme weather events such as droughts and floods.

[0343] Based on climate risk levels, carbon emission thresholds are dynamically adjusted downwards using preset compensation coefficients. In regions facing higher future climate risks (such as severe drought), production activities (such as agricultural irrigation) may be forced to consume more energy to sustain operations, indirectly increasing the potential risk of future carbon emissions. Therefore, current carbon emission requirements should be more stringent to incentivize these regions to either invest in climate adaptation measures (such as efficient water-saving irrigation) in advance or further reduce their current carbon footprint, thereby enhancing their ability to cope with future climate pressures.

[0344] Furthermore, it embodies the essence of climate resilience certification. Certification standards are no longer rigid but dynamically adjusted based on the future physical climate threats faced by the supply chain's location. This enhances the entire supply chain's ability to adapt to climate change. It establishes a risk-driven, threshold-fluctuating resilient decision-making framework, improving the scientific rigor and guidance of the certification process.

[0345] This application details an intelligent assessment method for international sustainable carbon certification. This method integrates geographic information systems, life cycle assessment, and uncertainty quantification technologies to achieve accurate carbon footprint accounting and robust compliance decisions for the supply chain. This application is particularly applicable to solving the certification challenges of exporting waste oil biofuels to the EU, providing reliable technical support for production enterprises.

[0346] In one example, this application provides a generated traceability digital certificate that includes not only the certification conclusion but also a detailed record of the scientific basis for carbon emission accounting, with a focus on the uncertainty information of carbon emissions. The traceability digital certificate includes: product identification information, a central estimate of life-cycle carbon emissions, a confidence interval characterizing the uncertainty of the central estimate, a data quality score, the basis for uncertainty analysis, sustainability quantitative indicators, and compliance verification conclusions.

[0347] Product identification information refers to a set of data used to uniquely identify a specific product batch or product, including but not limited to product batch code, manufacturer / supplier identifier, production date / time period, product category / model, geofence area identifier, etc.

[0348] The central estimate of life-cycle carbon emissions is the mean or median of the carbon emission probability distribution generated through Monte Carlo simulations. Assume 10,000 Monte Carlo simulations yield 10,000 possible carbon emission outcomes. These outcomes constitute a probability distribution. The central estimate might be the average of these 10,000 outcomes, for example, 152 tons. The confidence interval is the range of data calculated from this probability distribution based on a specified confidence level (e.g., 95%). For example, after sorting, the values ​​at the 250th and 9750th positions constitute a 95% confidence interval. For example, 140 tons to 165 tons.

[0349] Confidence intervals representing the uncertainty of the center estimate can be generated through the following steps:

[0350] Assign confidence levels to activity data and emission factors from different data collection channels; determine uncertainty parameters based on the probability distribution model followed by each confidence level; and generate confidence intervals by performing uncertainty propagation analysis through Monte Carlo simulation.

[0351] The data quality score is a comprehensive score calculated based on the credibility level of the input data.

[0352] Uncertainty analysis is based on the probability distribution model type (such as normal distribution or uniform distribution) used for the confidence level distribution of each input data, the number of iterations of the Monte Carlo simulation, and the convergence criterion.

[0353] Sustainability quantitative indicators: historical deforestation assessment results, social risk level, climate risk level, and threshold compensation coefficient.

[0354] Compliance verification conclusion: final certification result (pass / fail), certification evidence chain identifier, and blockchain evidence hash value.

[0355] It should be noted that the traceable digital certificate uses a machine-readable structured format (such as JSON-LD, XML), supporting automated verification and third-party auditing. The certificate's confidence interval generation process can be traced back to the specific input data's confidence level and uncertainty propagation analysis, ensuring the scientific transparency of the authentication results.

[0356] This application provides a global sustainable supply chain certification device, comprising:

[0357] At least one processor; and,

[0358] A memory communicatively connected to the at least one processor; wherein,

[0359] The memory stores instructions that can be executed by the at least one processor to enable the at least one processor to perform any of the above-described global sustainable supply chain certification methods.

[0360] The above description is merely an embodiment of this application and is not intended to limit this application. Various modifications and variations can be made to this application by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the technical principles of this application should fall within the protection scope of this application.

Claims

1. A global sustainable supply chain certification method, characterized in that, The method includes: Perform GIS location analysis on the production location of products in the supply chain to generate geofenced areas; Based on multi-source remote sensing data and GIS-LCA spatial analysis method, the sustainability quantitative index data of the geofence area are calculated. Obtain activity data and emission factors with credibility level identifiers from multiple data collection channels with credibility level classification; The life cycle carbon emissions of the product are calculated using the LCA carbon accounting model. For activity data and emission factors with confidence level labels, the probability distribution models followed by the activity data and emission factors and the corresponding uncertainty parameters are determined according to the preset mapping relationship between confidence level and probability distribution model, so as to generate confidence intervals for life cycle carbon emissions through Monte Carlo simulation. Based on the confidence interval, the sustainability quantitative indicator data and life cycle carbon emissions are verified for compliance with sustainability standards. In non-transportation phases, activity data with credibility level identifiers is obtained from multiple data collection channels with credibility grading, specifically including: Obtain direct monitoring data of geographical location, and when the direct monitoring data meets the first preset valid conditions, obtain activity data with a first credibility level identifier; If the first preset valid condition is not met, obtain the operational business data of the geographical location; When the operational business data meets the second preset validity conditions, activity data with a second credibility level identifier is obtained; If the second preset valid condition is not met, acquire geographic image data and geographic feature data related to the production space within the geographic fence area; The production activity of the product within the geographic fenced area is assessed based on the aforementioned geographic imagery data. Based on the production activity level and the geographical feature data, activity data with a third confidence level identifier is obtained; Based on the confidence interval, the sustainability quantitative indicator data and life cycle carbon emissions are verified for compliance with sustainability standards, specifically including: When the upper limit of the confidence interval for life cycle carbon emissions is less than a preset threshold, climate prediction data for the geofenced area over the future time period is extracted from the global climate change prediction model. Based on climate prediction data, predict the climate risk level of the geofenced area over a future period. Compensation is applied to preset emission thresholds based on climate risk levels; the climate risk level and the compensation coefficient are negatively correlated. When there is no history of deforestation, carbon emissions are less than the preset compensation emission threshold, and the social risk level is lower than the preset level, the product sustainability standard compliance verification is deemed successful.

2. The method according to claim 1, characterized in that, GIS location analysis is performed on the production location of products in the supply chain to generate geofenced areas, specifically including: Obtain the geometric center coordinates or precise geographic coordinates of the production site, and use these coordinates as seed points for the Voronoi algorithm. The production site is divided into Voronoi polygons to obtain the initial geofence area; Overlay hexagonal grids onto the initial geofence area and calculate the intersection area ratio of each hexagonal grid with the polygon; A hexagonal grid with significant spatial overlap with the polygon is selected to determine the geofence area; the significant spatial overlap is determined by comparing the ratio of the intersecting areas with a preset ratio threshold.

3. The method according to claim 1, characterized in that, When the sustainability quantification indicator is the history of deforestation, the sustainability quantification indicator data for the geographically fenced area is calculated, specifically including: The normalized vegetation index (NDVI) of each pixel in satellite images of the geofenced area within a preset historical time period and in the current time image is calculated to obtain the NDVI distribution map corresponding to each image. For each specified pixel in the NDVI distribution map of the current time image, extract the NDVI value from each historical image to obtain the NDVI time series of each specified pixel in the NDVI distribution map of the current time image. The NDVI time series of each specified pixel is analyzed based on the preset deforestation detection network model to obtain the deforestation labeling result for each specified pixel. If the ratio of the area marked as a deforestation pixel to the total area of ​​the geofence region is greater than a preset area ratio threshold, the deforestation index data is determined to be deforestation. If the ratio of the area marked as a deforestation pixel to the total area of ​​the geofence region is less than or equal to a preset area ratio threshold, the forest deforestation index data is determined to be non-deforested.

4. The method according to claim 1, characterized in that, When the sustainability quantification indicator is social risk, the sustainability quantification indicator data for the geofenced area is calculated, specifically including: Using kernel density estimation, a smoothing analysis is performed on the social risk data of the smallest administrative unit where the geographic fence area is located, generating a surface distribution map of social risk density. Spatial correlation analysis was performed on the surface distribution map of social risk density to identify statistically significant hotspot areas and generate a social risk heat map. The social risk level of the production site is obtained by analyzing the geographically fenced area based on the social risk heat map.

5. The method according to claim 4, characterized in that, Using kernel density estimation, a smoothing analysis is performed on the social risk data of the smallest administrative unit where the geofenced area is located, generating a surface distribution map of social risk density, specifically including: The population and labor index data are weighted and summed to obtain a composite index; the composite index is used as the weight of the center point of the smallest administrative unit where the geographic fence area is located; the center point of the administrative unit is used as the center of the kernel function to perform weighted kernel density estimation, and the population and labor risk density surface map is obtained. We estimate the weighted kernel density by taking the location of each negative news event in the news and public opinion data as the center of the kernel function, and obtain the public opinion risk density surface map. By taking the location of each historical risk event as the center of the kernel function, a weighted kernel density estimation is performed to obtain a historical risk density surface map; Generate a social risk density surface distribution map based on at least one of the population and labor risk density surface map, the public opinion risk density surface map, and the historical risk density surface map.

6. The method according to claim 4, characterized in that, The spatial correlation analysis of the social risk density surface distribution map specifically includes: Convert the surface distribution map of social risk density into grid data; Based on preset grid distances or adjacency relationships, determine the spatial weight matrix corresponding to the grid data; Based on the spatial weight matrix, calculate the Getis-Ord Gi statistic for each grid to obtain the Z score and P value for each grid; Multiple grids were classified based on Z-scores and P-values ​​to identify significant hot and cold areas.

7. The method according to claim 1, characterized in that, In non-transportation phases, emission factors with credibility level identifiers are obtained from multiple data collection channels with credibility grading, specifically including: In the local emission factor database of the administrative region or localized geographical unit to which the production site belongs, the emission factor type of the geographically fenced area is matched; If a match is successful, obtain the emission factor with the first confidence level identifier; If the matching fails, the characteristic parameters of the administrative region to which the production site belongs are obtained, and the standard emission factor of the next higher level administrative region to which the production site belongs is corrected according to the characteristic parameters; the characteristic parameters include at least one of energy structure, industrial composition, and climate zone type; If the correction is successful, an emission factor with a second confidence level identifier will be obtained; If the correction fails, the average factor in the general standard database is used to obtain the emission factor with a third confidence level identifier.

8. The method according to claim 1, characterized in that, The method further includes: Create and register digital identities for each participant in the supply chain; Receive material data digitally signed by each participating node and record the verified data as a transaction in the distributed ledger to generate a material balance data chain; At least one of the following is identified as a sustainability document: sustainable quantitative indicator data, life cycle carbon emission data, and geographic information of a geofence area. A first smart contract deployed on the blockchain is invoked. The first smart contract is configured to automatically verify the validity of the digital signature and the compliance of the data format of the sustainability document. After verification, the hash value and metadata of the document are recorded on the blockchain. After the product's sustainability standard compliance verification is passed, a second smart contract deployed on the blockchain is invoked. The second smart contract is configured to: collect on-chain evidence records related to product certification, including material balance chain records, sustainability document hashes, compliance verification conclusions, and generate a unique certification evidence chain identifier and a certification digital certificate including the certification evidence chain identifier, product information, and certification conclusions. The authentication evidence chain identifier and the hash value of the authentication digital certificate are submitted to an external trusted evidence storage system for storage.

9. The method according to any one of claims 1-8, characterized in that, The method further includes: generating a traceability digital certificate, the traceability digital certificate comprising: Product identification information is a set of data used to uniquely identify a specific product batch or product. The central estimate of life-cycle carbon emissions is the mean or median of the carbon emission probability distribution generated by Monte Carlo simulation. The confidence interval characterizing the uncertainty of the central estimate is generated through the following steps: Assign confidence levels to activity data and emission factors from different data collection channels; Uncertainty parameters are determined based on the probability distribution model followed by each confidence level; Uncertainty propagation analysis was performed using Monte Carlo simulation to generate confidence intervals.

10. A global sustainable supply chain certification device, characterized in that, include: At least one processor; as well as, A memory communicatively connected to the at least one processor; wherein, The memory stores instructions that can be executed by the at least one processor to enable the at least one processor to perform a global sustainable supply chain certification method as described in any one of claims 1-9.

Images