A remote sensing-based method for estimating carbon emissions

Abstract

The application discloses a kind of carbon emission estimation methods based on remote sensing, it is related to remote sensing technical field.The method includes: step 1: reading target area boundary, obtains a single time phase orthographic optical remote sensing image and is cut out to form basic image;The basic image is divided into several carbon activity units according to uniform grid;Step 2: on the scale of carbon activity unit, construct observation field atlas, process field atlas and income and expenditure field atlas, and align three field atlas with uniform index;By inversion selector, determine the only growth trajectory in each carbon activity unit and the discrete grade of the artificial carbon emission corresponding to the only growth trajectory;Step 3: on the scale of carbon activity unit, generate pixel level artificial carbon emission estimation graph and grid list.The present application can improve the accuracy of carbon source identification while maintaining spatial resolution, and enhance the stability and rationality of estimation results through multi-domain coupling constraint and scenario coverage.

Description

Technical Field

[0001] This invention relates to the field of remote sensing technology, and in particular to a method for estimating carbon emissions based on remote sensing. Background Technology

[0002] For a long time, the estimation of carbon emissions has mainly relied on statistical data and ground-based observation stations. Statistical data is usually compiled by governments or industry departments at various levels and can reflect the overall emission level of a region or industry over a certain time scale, but its spatial resolution is low, making it difficult to reveal differences in carbon emissions at finer-grained scales. Ground-based observation methods can provide continuous, fixed-point monitoring data with high accuracy at a single station or within a limited area, but due to limited deployment density, their coverage is significantly insufficient, making them unsuitable for large-scale and dynamic monitoring needs.

[0003] In recent years, with the development of remote sensing technology, researchers have attempted to use multi-source remote sensing data for carbon emission monitoring. For example, satellite-based carbon dioxide concentration retrieval methods can provide relatively uniform monitoring results over a large area. Some studies have achieved high spatial resolution carbon dioxide concentration distribution maps by fusing data from multiple satellites and using atmospheric transport models and Bayesian inversion algorithms. These methods can accurately identify emission characteristics in key areas on a global scale, but they still suffer from significant errors and susceptibility to cloud and aerosol interference in regional and urban scale applications. Furthermore, the long acquisition cycle of satellite data is not conducive to the formation of real-time or near-real-time dynamic monitoring of carbon emissions.

[0004] Another type of research attempts to improve the accuracy of emission factors by deploying continuous emission monitoring systems and combining them with artificial intelligence models to correct data from the perspective of industrial facilities and energy consumption. These methods have significant advantages in point source monitoring, capable of capturing emission surges on short timescales, but their application is limited by the cost and deployment conditions of monitoring equipment, making large-scale deployment difficult. At the urban scale, some studies combine low-cost sensor networks, traffic data, and remote sensing data to generate high spatiotemporal resolution carbon emission maps through data fusion methods. This type of research has made breakthroughs in spatial coverage, but the limited accuracy of low-cost sensors and high data noise lead to overall estimation biases. Summary of the Invention

[0005] In view of this, the present invention provides a remote sensing-based carbon emission estimation method that can improve the accuracy of carbon source identification while maintaining spatial resolution, and enhance the stability and rationality of the estimation results through multi-domain coupling constraints and scenario coverage. This effectively overcomes the shortcomings of existing technologies in terms of data fragmentation, model limitations, and insufficient utilization of spatial features, and has the value of wide application in carbon emission monitoring and policy support.

[0006] The technical solution adopted in this invention is as follows:

[0007] A remote sensing-based method for estimating carbon emissions, comprising:

[0008] Step 1: Read the boundary of the target area, acquire a single-temporal orthophoto remote sensing image and crop it to form a base image; divide the base image into several carbon activity units according to a uniform grid;

[0009] Step 2: Construct observation domain maps, process domain maps, and budget domain maps at the carbon activity unit scale, and align the three domain maps with a unified index; call the FORMIND model to generate a set of candidate growth trajectories covering multiple scenarios; establish a shading propagation zone using buildings and roads within the carbon activity unit as shading sources, and write the attributes of the shading propagation zone into the edge attributes of the three domain maps; apply observation consistency constraints, process consistency constraints, and budget consistency constraints in a fixed order, and update the candidate growth trajectory set and the shading propagation zone in conjunction with residual propagation; determine a unique growth trajectory and the discrete level of anthropogenic carbon emissions corresponding to the unique growth trajectory within each carbon activity unit using an inversion selector;

[0010] Step 3: Generate a pixel-level anthropogenic carbon emission estimation map and a gridded inventory at the carbon activity unit scale.

[0011] Furthermore, in step 1, the geographic reference point at the upper left corner of the base image is used as the starting point of the unified grid, and the unified grid is generated in a way that is strictly aligned with the pixel grid of the base image; the side length of the unified grid is set to an integer multiple of the side length of the base image pixel to avoid spatial misalignment and sampling deviation caused by cross-pixel cutting; the range of the unified grid is the smallest rectangular range that extends to completely cover the boundary of the target area.

[0012] Furthermore, in step 2, all generated carbon activity units are treated as a node set; adjacency edges are established according to the boundary adjacency relationships between carbon activity units to form a map of the observation domain; in the map of the observation domain, a node attribute area is established for each carbon activity unit, which includes: fragment position, texture description, edge description, and color combination description, for subsequent consistency verification; based on the geometric morphological cues of the base image within the carbon activity unit, object-level recognition is performed on the pixel set of each carbon activity unit through morphological segmentation and linear feature extraction, and the output object label is architecture, road, non-architecture, or non-road, and the object label is written into the node attributes of the map of the observation domain.

[0013] Furthermore, in step 2, in the process domain, for carbon activity units labeled as neither buildings nor roads, a process domain node is established for each such carbon activity unit, and process domain edges are established based on possible community succession relationships to form a graph of the process domain; in the budget domain, with all carbon activity units as nodes, budget domain edges are established according to spatial adjacency and potential source-sink flow directions to form a graph of the budget domain; a unified index lookup table is established to record the correspondence of each carbon activity unit in the observation domain node, process domain node, and budget domain node, as well as the correspondence of edges in the three-domain graphs, to achieve one-to-one alignment of the observation domain graph, process domain graph, and budget domain graph.

[0014] Furthermore, in step 2, for carbon activity units whose object labels are neither buildings nor roads, three types of driving information—site category, canopy layer, and void ratio—are inferred based on their texture description, edge description, and color combination description, and multi-scenario combinations are constructed. For each such carbon activity unit, the driving information and scenario combination of the carbon activity unit are input into the FORMIND model one by one, and the FORMIND model is run to obtain several time-ordered growth trajectories. All growth trajectories corresponding to the same carbon activity unit are aggregated into a candidate growth trajectory set for that carbon activity unit. For carbon activity units whose object labels are buildings or roads, no growth trajectory is generated, and only the source attribute is retained in the budget domain map. The number of the candidate growth trajectory set is linked to the unique identifier of the carbon activity unit in the unified index lookup table to ensure that the candidate growth trajectory set and the nodes of the observation domain map, process domain map, and budget domain map can be bidirectionally retrieved.

[0015] Furthermore, in step 2, the step of establishing a shading propagation band using buildings and roads in the carbon activity unit as shading sources and writing the attributes of the shading propagation band into the edge attributes of the three-domain map includes: in the observation domain map, taking the carbon activity unit with the object label as a road as the center, extracting the road skeleton direction, and generating fan-shaped or strip-shaped shading propagation bands on both sides along the skeleton direction; taking the large-volume continuous carbon activity unit with the object label as a building as the impedance center, generating a ring-shaped shading propagation band around the impedance center; for carbon activity units that simultaneously receive shading effects from roads and buildings, synthesizing shading propagation bands in the order of road guidance priority followed by building continuous inhibition; projecting the synthesized shading propagation bands onto the edges of the observation domain map, process domain map, and budget domain map, and writing the visibility level, succession velocity level, and flux impedance level into the edge attributes of the three-domain map, respectively.

[0016] Furthermore, in step 2, the loop order is set as follows: observation consistency constraints first, process consistency constraints second, and budget consistency constraints last. In each loop, observation consistency constraints are first applied to each carbon activity unit on the map of the observation domain. The texture description, edge description, and color combination description of the base image are compared with the visibility status of the candidate growth trajectory set in the corresponding time period. Growth trajectories that do not meet the visibility level requirements are deleted or the visibility level of the corresponding edge is reduced. Then, process consistency constraints are applied to the map of the process domain. The candidate growth trajectory set is screened according to the structural stability cues output by the FORMIND model. Growth trajectories that are unstable on the void structure or succession path are deleted or the succession velocity level of the corresponding edge is reduced. Finally, budget consistency constraints are applied to the map of the budget domain. The flux impedance level is locally adjusted according to the source-sink balance between adjacent carbon activity units. The edges corresponding to the candidate growth trajectories that do not meet the balance requirements are marked as pending repair.

[0017] Furthermore, in step 2, after a deletion, downgrade, or repair mark occurs in the graph of any domain, a change record for that domain's graph is generated. Residual propagation is triggered through the unified index lookup table, and the change record is mapped to the relevant nodes and edges of the other two domain graphs. The visibility level, succession rate level, and flux impedance level of the candidate growth trajectory set and the corresponding shading propagation band are adjusted accordingly. When all three types of consistency constraints are completed in this round, the next cycle begins until all carbon activity units have no repair marks in the three domain graphs and the size of the candidate growth trajectory set no longer changes.

[0018] Furthermore, in step 2, in the three-domain map state after the loop ends, for each carbon activity unit, the inversion selector is invoked to make decisions according to a fixed priority rule. Priority is given to resolving the remaining source-sink conflicts in the budget domain map, followed by resolving the remaining observation inconsistencies in the observation domain map, and then resolving the remaining structural instability in the process domain map. If multiple candidate growth trajectories simultaneously satisfy the constraints of the three-domain map, the growth trajectory with the highest consistency score is adopted as the unique growth trajectory of the carbon activity unit. After the unique growth trajectory is determined, based on the source-sink performance of the growth trajectory in the budget domain map and the degree of influence of the shading propagation band in the observation domain map and the process domain map, the anthropogenic carbon emission dispersion level of the carbon activity unit is determined by looking up a table, and the unique growth trajectory number and the anthropogenic carbon emission dispersion level are written back to the record corresponding to the carbon activity unit in the unified index lookup table.

[0019] By employing the above technical solutions, this invention achieves the following beneficial effects: Firstly, by constructing observation domain maps, process domain maps, and budget domain maps at the carbon activity unit scale and aligning them using a unified index, structured fusion of multi-source information is realized. This allows observational data from different sources and at different scales to be compared and updated within the same reference system, avoiding the problems of information fragmentation and scale mismatch in traditional methods. Secondly, this invention introduces the FORMIND model to generate a set of candidate growth trajectories covering multiple scenarios. Trajectories are then progressively screened and corrected through the cyclical application of observation consistency constraints, process consistency constraints, and budget consistency constraints. This approach ensures that the estimation results conform to direct observational evidence from remote sensing imagery, satisfy the rationality of ecological succession, and maintain a relative balance of source-sink relationships within the neighborhood, thereby improving the accuracy and stability of the estimation. Furthermore, by establishing shading propagation zones using buildings and roads as shading sources and incorporating them into the edge attributes of the three-domain maps, the shading effect can be finely characterized along directional paths, thus more realistically reflecting the impact of human activities on the visibility, succession rate, and flux exchange of the surrounding area, avoiding the loss of directional information caused by node averaging. Ultimately, at the carbon activity unit scale, pixel-level anthropogenic carbon emission estimation maps and gridded inventories are output. This allows the results to be used for high-resolution spatial visualization and to be integrated with statistical accounting systems, facilitating regulatory and policy applications. In summary, this invention, through the coupling constraints of multi-domain maps, the scenario coverage of trajectory sets, and the introduction of occlusion propagation mechanisms, achieves spatial refinement, process rationalization, and operational feasibility in carbon emission estimation. It effectively overcomes the shortcomings of existing technologies, such as insufficient spatial resolution, inaccurate carbon source identification, and insufficient integration of estimation models with spatial features. Attached Figure Description

[0020] Figure 1 This is a schematic diagram of the method flow for the carbon emission estimation method based on remote sensing in an embodiment of the present invention;

[0021] Figure 2 This is a schematic diagram illustrating the changes in carbon emission coefficients for different land use types in embodiments of the present invention;

[0022] Figure 3 This is a schematic diagram illustrating the relationship between remote sensing image resolution and recognition accuracy in an embodiment of the present invention;

[0023] Figure 4 This is a schematic diagram illustrating the impact of seasonal variations on vegetation carbon storage estimation in an embodiment of the present invention.

[0024] Figure 5 This is a schematic diagram of the standard deviation analysis of the spatiotemporal distribution of carbon emissions in an embodiment of the present invention. Detailed Implementation

[0025] All features disclosed in this specification, or all steps in all disclosed methods or processes, may be combined in any way, except for mutually exclusive features and / or steps.

[0026] Any feature disclosed in this specification, unless otherwise stated, may be replaced by other equivalent or similar features. That is, unless otherwise stated, each feature is merely one example of a series of equivalent or similar features.

[0027] refer to Figure 1 Step 1: Read the boundary of the target area, acquire a single-temporal orthophoto remote sensing image and crop it to form a base image; divide the base image into several carbon activity units according to a uniform grid.

[0028] Specifically, step 1 first reads the target region boundary as input vector surface data. This data should consist of one or more closed polygons, carrying a clear coordinate reference system and unit information. After reading, geometric validity checks and topology repair are performed, processing self-intersections, duplicate vertices, dangling edges, fragment loops, and non-closed loops in sequence. The direct benefit of this processing is that subsequent clipping will not produce narrow gaps or duplicate coverage, avoiding double counting or omissions at the boundary of the base image. When the target region boundary is composed of multiple non-adjacent closed polygons, its multi-part structure is maintained, and a spatial index is established to support block clipping and querying. The value of retaining the multi-part structure is that it can completely correspond to the actual management partitions, avoiding cross-regional interference caused by manual merging.

[0029] Select a single-temporal orthorectified optical remote sensing image from a public or authorized data source. Single-temporal means the acquisition time is concentrated on the same date or within the same observation period, not spanning multiple imaging dates. Choosing a single-temporal image reduces the differences in illumination, shadow, vegetation status, and human activities of surface features at different times, reducing uncertainty in subsequent estimations from the outset. Select an image with completed georegistration, radiometric correction, and geometric correction. Using orthorectified images allows for direct overlay on the same Earth reference plane as the target area boundary, reducing interpolation errors introduced by resampling and ensuring the base image retains the geometric and radiometric characteristics of the original image. When multiple candidate images exist for the same coverage area, prioritize the one with lower cloud cover, weaker haze influence, and significantly less stripe noise. This priority setting is based on the intuitive signal-to-noise ratio principle: clouds and haze smooth out the texture and color combinations of ground features, and stripe noise introduces false edges in the raster direction, both of which reduce the discriminability and repeatability of subsequent carbon activity units. If the coordinate reference system of the target area boundary is inconsistent with that of the single-temporal orthophoto image, the target area boundary is reprojected to a coordinate reference system and units completely consistent with the single-temporal orthophoto image, and closure and geometric validity checks are performed again. This ensures that the raster rows and columns correspond to the boundary coordinates in the same dimension during cropping, avoiding the accumulation of sub-pixel misalignment at the edges.

[0030] Using the target area boundary as a clipping mask, the single-temporal orthorectified optical remote sensing image is clipped and masked, outputting a base image that only covers the target area boundary. Pixels outside the boundary are assigned invalid values, while the original radiometric values ​​and georeferenced information are retained within the boundary. The direct advantage of using mask clipping instead of rectangular clipping is that the base image is strictly consistent with the management area, and the subsequently generated gridded list will not contain pixel records unrelated to the management area. When the target area boundary contains multiple non-adjacent closed polygons, each closed polygon is clipped separately, and a unified georeferenced image is maintained in the output, ensuring a clear one-to-many relationship between the base image and the target area boundary in space. This approach allows for parallel summarization of components and overall data at the carbon activity unit scale, facilitating verification and recalculation.

[0031] The starting point for the unified grid is determined to be the georeference point at the top left corner of the base image, and the unified grid is generated by strictly aligning it with the base image's cell grid. Choosing the top left corner as the starting point aligns with the storage order of most raster data from top to bottom and from left to right, which helps maintain integer alignment of rows and columns during writing and reading, avoiding row and column drift caused by cumulative rounding. The side length of the unified grid is set to an integer multiple of the side length of the base image's cells. This setting ensures that the boundary of the unified grid falls exactly on the cell boundary, rather than cutting off the cell, thus avoiding information pollution caused by mixed cells. For subsequent statistical and visualization operations at the carbon activity unit level, this integer multiple alignment significantly improves the mapping stability between the cell level and the unit level. The range of the unified grid is extended to the smallest rectangular area that completely covers the boundary of the target area. The significance of extending the range is that, regardless of the complexity of the target area boundary shape, the unified grid can form a closed row and column index, ensuring that the carbon activity units around the boundary are clearly defined and can be located.

[0032] A uniform grid is laid out within the base image, and the spatial relationship between the uniform grid cell and the target region boundary is determined cell by cell. Any uniform grid cell that intersects the target region boundary with any positive area is defined as a carbon active cell, and the intersection geometry of the uniform grid cell and the target region boundary is used as the effective spatial range of the carbon active cell. Retaining the strategy of arbitrary positive area intersection maximizes the preservation of the true details of the boundary and avoids "gaps" near the boundary that could cause discontinuities in subsequent statistical analysis. Uniform grid cells that only contact the target region boundary at the edge, forming a narrow geometry, are also defined as carbon active cells, and their true geometry is fully recorded. This process maintains the continuity of the adjacency relationships of carbon active cells in the observation domain map and the budget domain map, avoiding topological discontinuities caused by isolated pixels at the boundary. A unique identifier is generated for each carbon active cell, recording its row number, column number, starting vertex coordinates, center point coordinates, and effective area, and a correspondence is established between the carbon active cell and the rows and columns of the base image. This correspondence eliminates the need for resampling when performing bidirectional retrieval between the carbon active cell scale and the pixel scale, reducing the accumulation of bias caused by repeated interpolation. A consistency check is performed on the carbon active unit set. The check includes: no overlap between carbon active units, no holes within the target area boundary, and a one-to-one correspondence between the effective spatial extent of each carbon active unit and the rows and columns of the base image. This consistency check allows for timely detection of mismatches between boundary data and image data during implementation, preventing the amplification of errors in steps 2 and 3. The base image and the carbon active unit set are used as the final outputs of step 1. The former is used to maintain the complete radiometric and geometric information of the single-temporal orthophoto remote sensing image within the target area boundary, while the latter serves as the sole spatial carrier for subsequent processing and statistical analysis at the carbon active unit scale.

[0033] Choosing single-temporal imagery instead of multi-temporal fusion is to maintain the consistency of surface conditions. While multi-temporal fusion can cover a larger area or reduce cloud cover, it introduces shading and land cover variations from different dates, leading to unexplained temporal differences within the same carbon activity unit and increasing the burden of subsequent consistency constraints. Employing a strategy of strict alignment with base image pixels in unified grid generation allows for the direct inheritance of the base image's row and column structure, thus achieving lossless mapping between the pixel level and the carbon activity unit level. This lossless mapping eliminates the need for spatial resampling when generating the pixel-level anthropogenic carbon emission estimation map in step 3, reducing jagged edges and patch effects at boundaries. Unified grid units with only minimal intersection areas are still defined as carbon activity units, preserving the geometric details of the target area boundaries and ensuring that the topological relationships between the observation domain map and the budget domain map at the edges remain consistent with the true geographic boundaries. This directly reduces the ambiguity of adjacency determination in step 2. Using mask clipping instead of bounding rectangle clipping in the clipping stage can prevent background pixels that do not belong to the target area boundary from entering the subsequent processing flow, reduce the dilution effect of invalid data on statistical results, and improve the readability and audit efficiency of the gridded list in step 3.

[0034] When the uniform grid's side length is set to an integer multiple of the base image pixel side length, options such as 1x, 2x, and 4x can be selected. 1x provides higher spatial detail, suitable for urban built-up areas; 2x and 4x effectively reduce the number of carbon activity units, suitable for large-scale processing and rapid statistical analysis. Integer multiples maintain the overlap between pixel boundaries and uniform grid boundaries, avoiding pixel mixing. Using the top-left geographic reference point of the base image as the starting point, the starting point can be shifted to the right or downward by less than one pixel side length to minimize the generation of numerous elongated carbon activity units at the boundaries. This fine-tuning maintains row and column alignment while making the uniform grid more closely fit the boundary shape, thus reducing fragmentation of subsequent boundary units. For uniform grid units that only form elongated geometries at the boundaries, they can be unconditionally retained to maintain continuous adjacency between the observation domain map and the budget domain map; alternatively, a minimum width threshold can be set, and when the width of an elongated unit is below this threshold and does not contain a valid pixel center point, it is merged into an adjacent carbon activity unit, reducing noise in subsequent statistical analysis. When merging, maintain a consistent grid index to avoid duplicate identifiers.

[0035] Step 2: Construct observation domain maps, process domain maps, and budget domain maps at the carbon activity unit scale, and align the three domain maps with a unified index; call the FORMIND model to generate a set of candidate growth trajectories covering multiple scenarios; establish a shading propagation zone using buildings and roads within the carbon activity unit as shading sources, and write the attributes of the shading propagation zone into the edge attributes of the three domain maps; apply observation consistency constraints, process consistency constraints, and budget consistency constraints in a fixed order, and update the candidate growth trajectory set and the shading propagation zone in conjunction with residual propagation; determine a unique growth trajectory and the discrete level of anthropogenic carbon emissions corresponding to the unique growth trajectory within each carbon activity unit using an inversion selector;

[0036] In the specific implementation, the edge construction method of the observation domain map can be selected as either using only shared boundaries or using both shared boundaries and shared vertices. Using only shared boundaries can reduce long-distance interference introduced by higher-order adjacencies, making it suitable for urban areas with complex boundaries; while using shared vertices can improve connectivity, making it suitable for near-natural areas dominated by continuous surfaces. The succession sequence of the process domain map can be determined by either determining the direction based on the texture similarity of adjacent carbon activity units or by determining the direction based on the void ratio gradient of carbon activity units. The former is more sensitive to image quality and is suitable for high-quality base images; the latter is more robust to noise and is suitable for scenes with light fog. The shape of the shading propagation band can be selected as a bidirectional strip along the road framework, a fan-shaped extension along the road framework, or a star-shaped extension superimposed at intersections. The bidirectional strip is more in line with linear corridors, the fan shape is more suitable for networks with main and branch road levels, and the star-shaped extension can emphasize the multi-directional influence of intersections. The shading propagation band of buildings can be selected as a single ring or multiple rings; multiple rings are beneficial for expressing the differentiated influence of large buildings at different distances.

[0037] The edge attribute update strategy for the three-domain graph can be either conservative or aggressive. Conservative updates only reduce visibility, succession rate, or increase flux impedance when evidence is sufficient; their advantage is stability, but their disadvantage is slower convergence. Aggressive updates make adjustments in advance when evidence shows a clear bias; their advantage is faster convergence, but their disadvantage is the potential need for backtracking in subsequent rounds. The residual propagation trigger strategy can be either event-by-event immediate propagation or timed propagation after batch merging. Immediate propagation ensures rapid consistency of states between domains and is suitable for real-time processing; timed propagation reduces communication overhead and is suitable for batch processing. Both strategies use a unified index for location, resulting in consistent outcomes. The parallel candidate processing for the inversion selector can be either using the candidate with the highest consistency score or prioritizing those that have consistently appeared in historical rounds. The former emphasizes current evidence, while the latter emphasizes temporal stability. For stronger repeatability, a deterministic order of row and column numbers can be introduced when scores are the same to ensure that the same input yields the same output.

[0038] Specifically, firstly, the carbon active units obtained in step 1 are used as the sole source of nodes. To avoid misconnections caused by crossing gaps, each carbon active unit only establishes an adjacency relationship with adjacent carbon active units that have a real boundary contact on the base image. The reason for using boundary contact instead of distance threshold is that boundary contact directly corresponds to connected surface objects, which can stably carry the continuous texture and edge changes from the base image, reducing the risk of incorrect coupling of incoherent areas. The construction of the observation domain map includes two parts: nodes and edges. Nodes establish node attribute areas on a carbon active unit basis. The attribute area includes at least fragment location, texture description, edge description, and color combination description. Fragment location is used to stably locate in the base image to ensure repeatability of subsequent comparisons; texture description and edge description reflect the surface roughness and contour intensity of ground features, which can distinguish paving, vegetation, and structures; color combination description is used to capture the long-term stable characteristics of spectral combinations, which are still distinguishable under different illumination levels. Edges are created based on adjacency relationships. When two units share a real boundary, an undirected edge is established, and a visibility level field is reserved on the edge to accommodate the writing of occlusion propagation bands. Writing visibility on edges instead of nodes allows us to express the differences in visibility of the same unit in different directions, avoiding the loss of directional information due to averaging.

[0039] The process domain map establishes nodes only for carbon activity units whose object labels are neither buildings nor roads. This is because these units carry the temporal evolution of vegetation communities and are well-suited for aligning with the growth trajectories output by the FORMIND model. Edge placement follows a combination of adjacency and succession order: when two carbon activity units are adjacent and their texture and color combinations exhibit a gradual change from sparse to dense or from light to dark, directed edges are established to represent possible succession directions. Directional connections help extend the temporal sequence to spatial paths, creating coherent chains of candidate growth trajectories between adjacent units and preventing the same continuous land surface from being fragmented into contradictory successional states. The budget domain map uses all carbon activity units as nodes, with edges established based on spatial adjacency and potential source-sink flow directions. The determination of potential source-sink flow directions depends on object labels and neighborhood patterns; for example, buildings and roads tend to act as sources, while non-building and non-road units are more likely to reflect sink responses. Encoding this information into edges allows for targeted repair of local imbalances during subsequent consistency constraint phases. It prioritizes propagation and adjustment along adjacent directions with actual exchange channels, rather than spreading it evenly across the entire graph, thus preserving hotspots and gradients.

[0040] After the three-domain graph is established, a unified index is created. The unified index uses the unique identifier of each carbon activity unit as the primary key, recording the node number of each carbon activity unit in the observation domain graph, process domain graph, and budget domain graph, and maintaining an edge mapping table to establish a one-to-one correspondence between edges pointing to the same pair of carbon activity units in the observation domain graph, process domain graph, and budget domain graph. The value of the unified index lies in the fact that any addition, deletion, or downgrading operation in one domain can be deterministically mapped to the other two domains, ensuring that residual propagation is neither lost nor duplicated, and that the timing of the same event is consistent across different domains. To eliminate uncertainties in implementation, the unified index requires generating numbers in a fixed order of row number priority followed by column number, and restoring this order after any parallel processing or block merging, thereby ensuring that the same input yields the same output. This determinism not only facilitates verification but also reduces the probability of oscillations during the cyclic constraint phase, because each state update returns to the same reference frame. In multi-component target region boundary scenarios, the unified index uses hierarchical encoding for component numbers, row numbers, and column numbers, maintaining global uniqueness across components while allowing for fast subset retrieval within a single component. For cross-tile processing, the unified index employs a shared edge rule at tile boundaries: when carbon activity units on two tiles share the same real boundary, only one edge record is retained and uniformly assigned during the merging phase, avoiding degree inflation caused by duplicate edge connections. These measures ensure that the topology of the three-domain graph remains consistent at any scale, thus providing a stable platform for subsequent candidate growth trajectory set writing, occlusion propagation band writing, and the three types of consistency constraints.

[0041] Based on the scale response of both color combination and texture descriptions, this carbon activity unit was classified into three site categories: herbaceous, shrubby, and tree-like. Color combination descriptions reflect the relative differences in leaf area density and water content, while texture descriptions are sensitive to the fine structure of leaves and branches at small scales and to the aggregation of canopies and patches at medium scales. Herbaceous tendencies typically exhibit relatively uniform color combinations and low-amplitude texture fluctuations; shrubby tendencies show blocky textures and moderate undulations; and tree-like tendencies have distinct patch layers and stronger texture undulations. The advantage of this approach is that stable classifications can be obtained in a single time phase, without relying on external surveys, and the correspondence with community types in the FORMIND model is clear. The vertical structure of this carbon activity unit was graded using the density and crossover of edge descriptions at multiple scales. Single-layer canopies are common in herbaceous tendencies or low shrublands, with sparse and uniformly oriented edge lines; double-layer canopies are common in shrublands covered by sparse trees or young forests, with two dominant edge directions; and multi-layer canopies are common in mature forests, with dense and interwoven edges. The reason for using the "density and intersection" of edges as clues is that the complexity of vertical structures is transformed into the superposition and occlusion of contours in planar images. The more complex the vertical structure, the more edges are generated in the same area and the more dispersed the directions are. This mapping is stable and reproducible.

[0042] Within the carbon activity unit, the area ratio of connected slices with low texture, low edge intensity, and bright color combinations is statistically analyzed as a measure of void ratio. A higher void ratio indicates more update space, lower competition intensity, and a more active regeneration process in the FORMIND model; a lower void ratio indicates canopy closure, full resource utilization, and a more slow succession progression. Using connected slices instead of scattered pixels helps eliminate the interference of noise and shadow points on the results, making this indicator more stable in reflecting the real structure. Carbon activity units near buildings and roads are susceptible to occlusion and mixing. First, based on the object label and fragment position, the core region is shrunken by a certain pixel distance at the edge. The texture description, edge description, and color combination description are then recalculated, and the results of the core region are used as the standard, with the outer ring region only used as a reference. If there is a significant discrepancy between the core and the outer ring, the confidence correction for the carbon activity unit is performed using the majority consensus conclusion of adjacent carbon activity units. This approach can reduce misjudgments caused by boundary mixing while maintaining spatial continuity. The site category, canopy level, and void ratio are encoded as discrete levels and written into a unified index along with the unique identifier of the carbon activity unit, ensuring that subsequent multi-scenario combinations and FORMIND model inputs can be directly retrieved.

[0043] Without altering single-phase observations, this study uses discrete binning and combination of three types of driving information to cover the possible structural states and achievable successional paths of a carbon activity unit. The starting point for setting up multiple scenarios is to avoid interpreting a single observation as the only true state, allowing subsequent consistency constraints to "select the best rather than be forced" within a wider candidate space, thereby reducing cascading errors caused by early assumption errors. Based on the overall distribution of the project area, site categories are divided into fixed set of categories, canopy layers are divided into a sequence of levels from simple to complex, and void proportions are divided into a sequence of levels from low to high. Using intra-regional distribution for binning instead of applying external thresholds can automatically adapt to contrast differences between different cities or forest types, avoiding systematic bias caused by overall brightness or darkness. A Cartesian combination of the three discrete levels yields an initial set of multiple scenarios. Two screening processes are then applied. The first is neighborhood consistency screening: if a scenario is completely opposite to the mainstream level of an adjacent carbon activity unit and has no transitional unit support, it is temporarily not adopted to avoid creating spatially isolated scenarios. The second step is self-consistency screening: if a scenario is internally contradictory, such as binding herbaceous tendencies to multi-layered canopies, it is eliminated. The purpose of these two screenings is to reduce "theoretical feasibility" to "spatial feasibility and structural feasibility," thereby reducing unnecessary subsequent calculations.

[0044] For each selected scenario, the initial community and driving configuration required for the FORMIND model are constructed. Site category determines the species pool and physiological constraints, canopy hierarchy determines the initial hierarchical structure and canopy cluster distribution, and void proportion determines the initial void network and regeneration intensity. The configuration determined by these three factors is used as input for an independent run, ensuring consistency between the model evolution path and image observation support. The FORMIND model is run for each scenario, recording the time series from the initial to mature stage. To facilitate subsequent comparisons, only slices from several representative periods are retained, such as the initial rapid recovery stage, the mid-term structural differentiation stage, and the late stable stage. Using "representative slices" significantly reduces storage and comparison costs without sacrificing diversity, while preserving key resolution of successional stages. For each time-ordered growth trajectory, a description consistent with observations is generated, including structural hierarchy, canopy cluster connectivity, void network morphology, and expected visibility labels. Converting the trajectory into an observational description allows for item-by-item verification during the observation consistency constraint stage, avoiding dimensional mismatches caused by directly comparing the model's internal state to the original image values. All time-ordered growth trajectories generated by all scenarios within the same carbon activity unit are compiled into a candidate growth trajectory set using a unified index. Each growth trajectory is assigned a unique number, and its source scenario, representative time period sequence, and observational description summary are recorded. This organization allows for precise location at the trajectory level for subsequent deletion, downgrading, rollback, and write-back. If a growth trajectory's observational description across all representative time periods deviates entirely from the fragment position, texture description, edge description, and color combination description of that carbon activity unit, and there are no similar paths in adjacent carbon activity units, it is removed as an unreliable candidate. This is a preliminary health check, aimed at eliminating paths that clearly conflict with the current single-phase evidence in advance, reducing the burden of subsequent consistency constraints, and preventing extreme values ​​from dominating residual propagation.

[0045] Employing three types of cues—color combination description, texture description, and edge description—explicits the classic correspondence between "spectrum to matter," "texture to structure," and "edge to morphology." Under single-temporal conditions, relying on seasonal curves is insufficient; stable signals must be sought within spatial distribution. These three cues stably correspond to water content and greenness characteristics, layering and roughness, and contour and occlusion, respectively, and their combination provides sufficient differentiation of community states. Using scenario combinations instead of a single interpretation shifts uncertainty forward. Multiple possibilities are initially allowed to coexist, then eliminated through subsequent consistency constraints. Compared to directly selecting a single interpretation, this approach significantly reduces the irreversibility of early errors, giving the system room to maneuver when facing local anomalies. Replacing complete continuous sequences with representative time-slices strikes a balance between interpretability and resource overhead. Subsequent constraints primarily rely on alignment checks with image evidence, focusing on whether the stage structure matches, rather than meticulously tracking every minute time step. Retaining only a few slices is sufficient to support robust comparisons.

[0046] In addition to edge density and intersection conditions, multi-scale self-similarity can be introduced as an auxiliary clue for canopy hierarchy inference. Multi-layered canopies exhibit similar texture patterns at different scales, resulting in high self-similarity; single-layered canopies show rapid texture decay with scale changes, leading to low self-similarity. This clue is only used for parallel judgment and does not alter the original hierarchy. When the proportion of carbon activity units labeled as neither buildings nor roads is too low, a neighborhood migration strategy can be enabled: for a small number of isolated carbon activity units of this type, the distribution of the three types of driving information from neighboring similar carbon activity units is referenced to smooth the data, generating a more robust scenario set to avoid excessive sparsity in the candidate growth trajectory set affecting subsequent convergence.

[0047] If the object is labeled as a carbon activity unit of a road, its skeleton direction is extracted based on edge direction statistics and thinning within the fragment location. Thinning can converge the road surface connectivity structure to a single-pixel-wide centerline, while edge direction statistics are used to stabilize the main direction in the presence of noise. The main direction is determined first because traffic corridors have a significant directional impact on visibility and process propagation. Masking and interference along the corridor direction are continuous and far-reaching, with faster lateral attenuation. If the direction is not distinguished, linear propagation will be incorrectly averaged across the circumference, reducing the resolution of subsequent constraints. Strip-shaped or fan-shaped masking propagation bands are generated on both sides of the road along the skeleton direction. The strip shape is used for a single road segment, while the fan shape is used when the skeleton extends in multiple directions near intersections. The bandwidth and angle are determined based on the combined strength of fragment location, texture description, and edge description: the stronger the texture and edge, the more obvious the road surface material, traffic traces, and hard boundaries, resulting in a relatively wider range of influence. Using this adaptive range allows wide and narrow roads to show differences in the masking propagation bands, avoiding over-masking or under-masking caused by a uniform constant.

[0048] For carbon activity units labeled as buildings, large-scale continuous carbon activity units are first identified in the observation domain map through connectivity determination. Their outer boundaries and approximate center points are calculated, and a ring-shaped shading propagation zone is generated centered on this continuum. Buildings are used in a ring shape rather than directional expansion because large entities tend to provide balanced circumferential shading. The three-dimensional blocks create stable line-of-sight obstructions and airflow disturbances in all directions, exhibiting strong near-distance effects that gradually decrease with distance. Describing this pattern as a ring shape provides directional, unbiased, and continuous suppression without introducing complex geometry. Within the same or adjacent carbon activity units, if road shading propagation zones overlap with building shading propagation zones, synthesis is performed in order of priority: road guidance followed by continuous building suppression. Specifically, the dominant strip or fan-shaped morphology along the road framework is preserved, and the building ring-shaped influence is superimposed as a baseline to enhance the overall circumferential shading level. The advantage of this approach is that it preserves the directionality of the corridor while maintaining the uniform suppression of the volumetric entities, resulting in both "distant guidance" and "near equilibrium" after stacking, which is beneficial for subsequently providing differentiated levels in different side directions.

[0049] When the road framework has small breaks or offsets, linear interpolation is performed on the shading propagation zone to restore connectivity. When narrow gaps exist between building continuums but appear as shadowed areas of the same volume in the image, the annular shading propagation zone is expanded and the contact rings are merged. This avoids abrupt changes in shading levels caused by minor discontinuities, maintaining consistency with the real scene. The synthesized shading propagation zone is projected onto the edges of the three-domain map. For the observation domain map, the spatial relationship between each edge and the shading propagation zone is calculated and written into the visibility level; for the process domain map, the succession velocity level is written; and for the budget domain map, the flux impedance level is written. Edges are written instead of nodes because the strength of the shading effect varies when the same carbon activity unit faces different adjacent directions. Edge-level representation can preserve directional differences, ensuring that subsequent constraints take effect on the actual propagation path, rather than being averaged out at nodes. The levels are discretely graded, combined with distance and morphology. The closer an edge is to the core area of ​​the shading propagation zone, the higher its level. If the shape is located in the main direction of the road framework, it has a higher visibility level in the observation domain map, a lower succession rate level in the process domain map, and a higher flux impedance level in the budget domain map. If it is located in a direction orthogonal to the road, the level is lowered. This directional differentiated grading reflects the common characteristics of "linear reach and rapid lateral decay" in the edge attributes, improving the selectivity and stability of subsequent cyclic constraints. When the same edge is covered by multiple shading propagation zones, it is first determined whether it includes the direct direction of the road framework. If it does, the road shading propagation zone takes precedence, and the baseline enhancement of the building ring shading is superimposed. If it does not, the building ring shading takes precedence. This order makes the road network structure clearly expressed in the map, while preserving the long-term suppression effect of the built environment, making it easier to obtain an intuitive level distribution in complex blocks.

[0050] The shape of road shading propagation zones can be strip-shaped or fan-shaped, with a star-shaped expansion at intersections to enhance multi-directional guidance. Building shading propagation zones can be single-ring or multi-ring, with multi-ring used for particularly large continuums to distinguish near, middle, and far-reaching influence zones. The grading system can use three or five levels; three levels facilitate rapid processing, while five levels allow for fine-tuning. In areas with a high proportion of natural environment, the priority of road guidance can be appropriately reduced, relatively increasing the contribution of buildings and terrain continuums to the grading. For projects processing across tiles, a half-bandwidth buffer can be added at tile boundaries. Shading propagation zones are first calculated within the buffer and deduplicated during the merging phase to avoid grading breaks at boundaries. For long-running systems, the edge attribute grading from the previous round can be recorded in a unified index, with change amplitude limits applied during new rounds of writing to prevent frequent grading oscillations from affecting subsequent convergence.

[0051] In the observation domain map, carbon activity units labeled as roads and carbon activity units labeled as buildings form a set of shading sources. All are bound by a unique carbon activity unit identifier and a unified index to ensure consistency in subsequent write-backs and backtracking. For carbon activity units labeled as roads, the road skeleton direction is determined based on the statistical and fine-line results of edge orientation within the segment location. Determining the main direction first is valuable because the influence of traffic corridors exhibits significant anisotropy, extending far along the corridor and attenuating rapidly laterally. If the direction is ignored, shading will be improperly distributed circumferentially, reducing the sensitivity of subsequent constraints to the actual propagation path. Strip-shaped or fan-shaped shading propagation bands are generated on both sides of the road along the skeleton direction. Strips are used for continuous straight sections, and fans are used near intersections and bends. The bandwidth and angle are adaptively set based on the strength of texture and edge descriptions: the more stable and intense the texture and edges, the harder and more accessible the corridor, and the greater the influence range can be. This allows main roads and branch roads to show differences in shading, avoiding excessive or insufficient shading caused by a uniform width. For carbon activity units labeled as "buildings," large-scale continuous carbon activity units are first merged based on connectivity. The outer boundary and approximate center point are then calculated, generating a ring-shaped shielding propagation zone around the continuum. A ring-shaped rather than directional expansion is used because the circumferential line-of-sight obstruction and near-ground disturbance of the mass entity are relatively uniform, smoothly attenuating with distance; a ring-shaped representation better reflects this continuous and nearly uniform characteristic in all directions. If road shielding propagation zones overlap with building shielding propagation zones within the same domain, synthesis is performed prioritizing road guidance and followed by continuous building suppression: the dominant strip or fan-shaped form is retained along the road framework direction, while the circumferential horizontal is raised using the building ring-shaped shielding as a baseline. This preserves the reachability of linear guidance without sacrificing the near-domain consistency of mass suppression.

[0052] Using each edge connecting adjacent carbon activity units in the three-domain map as a carrier, the spatial and directional relationships between the edge and the shading propagation zone are calculated. Visibility level is written into the edge attributes of the observation domain map, succession velocity level into the edge attributes of the process domain map, and flux impedance level into the edge attributes of the budget domain map. Edge-level writing is used because the same carbon activity unit is affected differently when facing different adjacent directions; edges can preserve directional differences and avoid information loss caused by node averaging. When the edge direction is approximately in the same direction as the road framework, the visibility level is increased in the observation domain map, the succession velocity level is decreased in the process domain map, and the flux impedance level is increased in the budget domain map; when the edge direction is approximately orthogonal to the road framework, the above levels are appropriately lowered. For building rings, levels are set segment by segment according to the distance to the continuum boundary, with higher levels for closer distances. This approach accurately maps the commonalities of "linear reach, rapid lateral decay, and strong near-term intensity" to edge attributes. If the same edge falls into multiple shading propagation zones simultaneously, priority is given to determining whether it includes the direct direction of the road framework. If it does, road shading is prioritized, with building baselines superimposed; otherwise, building rings are prioritized. If parallel zones still exist, priority is given to those closer to the center of the shading source, followed by the fixed order of carbon activity unit row and column numbers to ensure the same input yields the same output. For each carbon activity unit, the outgoing edge level is checked to ensure it changes smoothly radially or in a strip pattern along space. If a checkerboard-like jump occurs, backtracking is performed to check for breaks or mis-superimposed shading propagation zones and corrections are made. After passing these checks, the level and source shading propagation zone identifiers are written into the three-domain map and the writing time is registered in a unified index for subsequent residual propagation and historical comparison.

[0053] The algorithm prioritizes observation consistency constraints, followed by process consistency constraints, and then expenditure consistency constraints. Prioritizing observation consistency constraints allows for early narrowing of the candidate space based on the closest evidence, reducing subsequent computational burden. Prioritizing expenditure consistency constraints allows for minimal local balancing after preserving more observational and structural information, preventing premature smoothing that masks true differences. For each edge, the visibility of the candidate growth trajectory set is compared item by item with the corresponding time period using fragment location, texture description, edge description, and color combination description. When an occlusion propagation band indicates that strong textures or clear edges should not appear in that direction, but the base image provides contradictory evidence, the visibility level of that edge is reduced, and the related growth trajectory is marked for deletion. When an occlusion propagation band indicates that the direction should be visible but image evidence is insufficient, the priority of that growth trajectory is reduced. This prioritizes paths compatible with the image, reducing false consistency caused by assumptions. Along the directional edges of the process domain map, the structural continuity of the same candidate growth trajectory and the sequential progression between adjacent units are examined over representative time periods. Once an unreasonable stage jump or a long-term conflict with the overall succession direction of the neighborhood occurs, the succession velocity level of that edge is reduced, and the relevant growth trajectory is marked as to be deleted or downgraded. Using structural continuity signals, unreliable paths can be eliminated without introducing external data. The closure of source-sink exchanges is checked within the local neighborhood of the feed-in domain map. When persistent unidirectional accumulation occurs and is difficult to explain by observation and process, the flux impedance level of the relevant edge is increased to limit rapid leveling across that edge, and the growth trajectory involved in the accumulation is marked as to be repaired, allowing subsequent rounds to complete adjustments along the direction of lower impedance. Adopting neighborhood-by-neighbor correction rather than a global one-time balancing preserves hotspots and gradients, preventing the smoothing out of true peak values. After any domain is deleted, downgraded, or marked as to be repaired, a change record is generated and mapped to the other two domains through a unified index. The affected edges and growth trajectories are added to the update queue and processed in queue order, avoiding a full graph recalculation. This ensures synchronization across the three domains, controls computational load, and provides a traceable cross-domain response for each modification. After completing one round of three-class consistency constraints, if there are no markers to be repaired within the three-domain graph, or if the size of the candidate growth trajectory set remains unchanged for several consecutive rounds, the loop terminates. If necessary, the maximum number of rounds can be increased as an upper limit to ensure that the processing can be completed controllably in complex scenarios.

[0054] In the three-domain map state after the cycle terminates, the inversion selector is invoked for each carbon activity unit to make a decision. The inversion selector executes according to a fixed priority: first, it resolves residual source-sink conflicts in the budget domain map; second, it resolves residual observation inconsistencies in the observation domain map; and third, it resolves residual structural instability in the process domain map. This order is considered because the budget domain map directly affects the local closure of the discrete levels of anthropogenic carbon emissions. If image or structural details are resolved first, it may lead to a large-scale backtracking requirement in the end. When multiple candidate growth trajectories simultaneously satisfy the constraints of the three-domain map, the one with the highest consistency in the degree of influence of the determined trajectory in the shading propagation zone within the neighborhood is selected as the unique growth trajectory. Selecting the trajectory with the highest neighborhood consistency allows for a continuous spatial interpretation of the results, avoiding the alternating chessboard pattern. After determining the unique growth trajectory, based on its source-sink representation in the budget domain map and the influence of shading propagation zones in the observation and process domain maps, the discrete level of anthropogenic carbon emissions for that carbon activity unit is obtained by looking up a table. The unique growth trajectory number and the discrete level of anthropogenic carbon emissions are then written back to the unified index. Using a table lookup method ensures comparability of the same level across different regions, facilitating the direct generation of pixel-level anthropogenic carbon emission estimation maps and gridded inventories in step 3. When small-scale discontinuities exist in the road framework or building continuums are interrupted by narrow gaps, resulting in a consistent shadow on the image, linear interpolation or ring merging of shading propagation zones is performed to avoid non-physical level jumps and ensure spatial coherence and interpretability.

[0055] Step 3: Generate a pixel-level anthropogenic carbon emission estimation map and a gridded inventory at the carbon activity unit scale.

[0056] In the specific implementation, a new blank raster is created as a pixel-level anthropogenic carbon emission estimation map. Its row count, column count, pixel size, and spatial origin are completely consistent with the base image, and it inherits the georeference of the base image. A specific invalid value code is set specifically for pixels outside the target area boundary to distinguish between unassigned areas and true zero-value areas. For each carbon activity unit, its minimum row number, maximum row number, minimum column number, and maximum column number are first obtained from a unified index to form a bounding window. Then, the center point is determined pixel by pixel within the window: pixels whose center point is within the effective spatial range of the carbon activity unit are included in the pixel set of that carbon activity unit. Using center point determination instead of area ratio allocation has three advantages: first, the determination result is not sensitive to minor boundary fluctuations, and repeated runs can obtain consistent results; second, it avoids the generation of fractional weights, and both assignment and statistics are integer-based, facilitating auditing; third, a one-to-one association is formed between pixels and carbon activity units, making subsequent tracing simple and clear. Individual pixel center points may fall on shared boundaries due to numerical rounding errors. A fixed and publicly known decision-making order is applied: carbon activity units with smaller row numbers are selected first, and if the row numbers are the same, carbon activity units with smaller column numbers are selected. This order is consistent with the unified index, ensuring that the same input yields the same output and facilitating reproducibility under different environments. After assignment, the discrete level of anthropogenic carbon emissions corresponding to the carbon activity unit is written to that cell; cells outside the target area boundary are uniformly assigned invalid value codes.

[0057] When a carbon activity unit exhibits a narrow, elongated shape at its boundary, resulting in single-cell holes in the pixel set, a combination of center point determination and single-layer adjacency filling is used for repair within that carbon activity unit. The reason for limiting it to single-layer adjacency is to only fill the discrete gaps caused by rasterization, without crossing the true boundary and thus altering the extent of ground features. After repair, the pixel values ​​within the carbon activity unit are checked again to ensure they all equal the discrete levels of its anthropogenic carbon emissions. The number of pixels and the number of pixel value types are counted for each carbon activity unit. If the number of pixel value types is greater than 1, it indicates an anomaly in boundary attribution or repair, and the process is reverted to the mapping and assignment steps. If the number of pixels does not match the bounding window and effective area calculation recorded in the unified index, the implementer is prompted to review the alignment of the target area boundary with the base image. A full-map check ensures that invalid value codes only appear outside the target area boundary, thus confirming clear clipping and filling boundaries. A gridded list is generated using a unified grid as the recording unit, with each row corresponding to one carbon activity unit. To ensure a closed loop for verification and traceability, the gridded list must contain at least the following fields and be stored in a fixed order: unique identifier for carbon activity unit, row number, column number, unique growth trajectory number, anthropogenic carbon emission dispersion level, effective area, number of pixels, minimum row number, maximum row number, minimum column number, maximum column number, and coordinate reference system identifier. These fields are set up because: the unique identifier, row number, and column number ensure direct connection to the unified index; the unique growth trajectory number and anthropogenic carbon emission dispersion level support the decision results from step 2; the effective area and number of pixels provide quantitative basis for subsequent accounting and spot checks; row and column boundaries facilitate rapid location; and the coordinate reference system identifier ensures no misreading occurs when sharing across systems. The gridded list is verified row by row to ensure that the number of pixels for each carbon activity unit is consistent with the pixel count within the same spatial range in the pixel-level anthropogenic carbon emission estimation map, and that the anthropogenic carbon emission dispersion level is consistent with the pixel values ​​within that range. Any inconsistent entry is marked and backtracked to the mapping or assignment step for correction until both are consistent. This approach adheres to the principle of "two-way verification," ensuring that tables can be derived from images while images can be reconstructed from tables, thus establishing a mutual verification relationship at the data level.

[0058] Figure 2 The graph shows the differences in carbon emission coefficients for six major land use types. As can be seen from the graph, construction land has the highest carbon emission coefficient, reaching [value missing]. They exhibit significant carbon source characteristics. Farmland and grassland show relatively low positive emission coefficients, respectively. In contrast, forest land exhibited the strongest carbon sequestration function, with a carbon emission coefficient of [missing value]. Wetlands also exhibit significant carbon sequestration characteristics, with a coefficient of [missing value]. The carbon emission coefficient of unused land is close to neutral, at only [value missing]. This result indicates that different land use types contribute significantly to the carbon cycle, providing an important parameter basis for the classification and management of carbon activity units and emission estimation.

[0059] Figure 3 This study reveals the impact of remote sensing image spatial resolution on the accuracy of ground feature identification. As image resolution decreases from 0.5m to 100m, the accuracy rate for building identification drops sharply from 96.8% to 38.2%, and for roads, it drops from 92.3% to 24.8%. Within the high-resolution range (0.5-2m), the accuracy rate remains above 90%. When the resolution is coarsened to 5-10m, a significant inflection point in accuracy is observed. Building identification is slightly less sensitive to resolution than road identification, which is related to the relatively large geometric features of building targets. This relationship curve provides a quantitative basis for the selection of remote sensing data sources, indicating that achieving high-precision carbon activity unit identification requires remote sensing image data with a resolution better than 10m.

[0060] Figure 4 The seasonal variation patterns of carbon storage in broadleaf and coniferous forests are described. Broadleaf forest carbon storage exhibits a clear seasonal fluctuation, gradually increasing from 145.2 tC / ha in January, peaking at 241.2 tC / ha in August, and then declining to 152.4 tC / ha in December, with an annual variation of 96.0 tC / ha. Coniferous forests show a relatively gentler seasonal variation, with carbon storage fluctuating between 182.5 and 218.9 tC / ha, and an annual variation of only 36.4 tC / ha. Both forest types reach their peak carbon storage in summer, but the seasonal sensitivity of broadleaf forests is significantly higher than that of coniferous forests. This variation pattern reflects the differences in phenological characteristics between different vegetation types and provides an important reference for carbon storage estimation from time-series remote sensing data and seasonal correction of FORMIND model parameters.

[0061] Figure 5 The variability of carbon emissions in the spatiotemporal distribution was quantified. Spatially, the standard deviation of carbon emissions was largest in winter, reaching [value missing]. Spring and autumn are respectively The minimum in summer is The standard deviation over time is generally lower than the standard deviation over space, and is higher in winter. Autumn is Springtime is The lowest in summer is The results indicate that the spatial heterogeneity of carbon emissions is stronger than their temporal variability, with both spatial and temporal variability reaching their maximum in winter. This variability pattern reflects the combined effects of human activities such as heating and natural processes such as vegetation dormancy, providing a statistical basis for parameter setting of the residual propagation mechanism and the construction of the budget domain map.

[0062] This invention is not limited to the specific embodiments described above. The invention extends to any new feature or combination disclosed in this specification, as well as any new method or process step or combination disclosed herein.

Claims

1. A remote sensing-based method for estimating carbon emissions, characterized in that, The method includes: Step 1: Read the boundary of the target area, acquire a single-temporal orthophoto remote sensing image and crop it to form a base image; divide the base image into several carbon activity units according to a uniform grid; Step 2: Construct observation domain maps, process domain maps, and budget domain maps at the carbon activity unit scale, and align the three domain maps with a unified index; call the FORMIND model to generate a set of candidate growth trajectories covering multiple scenarios; establish shading propagation bands using buildings and roads within the carbon activity unit as shading sources, and write the attributes of the shading propagation bands into the edge attributes of the three domain maps; apply observation consistency constraints, process consistency constraints, and budget consistency constraints in a fixed order, and update the candidate growth trajectory set and the shading propagation band in conjunction with residual propagation; determine a unique growth trajectory within each carbon activity unit using an inversion selector; after the unique growth trajectory is determined, determine the anthropogenic carbon emission dispersion level of the carbon activity unit by looking up a table based on the source and sink performance of the growth trajectory in the budget domain map and the degree of influence of the shading propagation band in the observation domain map and process domain map; Step 3: Generate a pixel-level anthropogenic carbon emission estimation map and a gridded inventory at the carbon activity unit scale; In step 1, the geographic reference point at the top left corner of the base image is used as the starting point of the unified grid, and the unified grid is generated by strictly aligning it with the pixel grid of the base image. The side length of the unified grid is set to an integer multiple of the side length of the base image pixel to avoid spatial misalignment and sampling deviation caused by cross-pixel cutting. The range of the unified grid is the smallest rectangular range that extends to completely cover the boundary of the target area. In step 2, all generated carbon activity units are treated as a node set; adjacency edges are established according to the boundary adjacency relationships between carbon activity units to form a map of the observation domain; in the map of the observation domain, a node attribute area is established for each carbon activity unit, which includes: fragment position, texture description, edge description, and color combination description, for subsequent consistency verification; based on the geometric morphological cues of the base image within the carbon activity units, object-level recognition is performed on the pixel set of each carbon activity unit through morphological segmentation and linear feature extraction, and the output object label is either building, road, non-building, or non-road, and the object label is written into the node attributes of the map of the observation domain; In the process domain, for carbon activity units labeled as neither buildings nor roads, a process domain node is established for each such carbon activity unit, and process domain edges are established based on possible community succession relationships to form a graph of the process domain. In the budget domain, with all carbon activity units as nodes, budget domain edges are established according to spatial adjacency and potential source-sink flow directions to form a graph of the budget domain. A unified index lookup table is established to record the correspondence of each carbon activity unit in the observation domain node, process domain node, and budget domain node, as well as the correspondence of edges in the three-domain graphs, to achieve one-to-one alignment of the observation domain graph, process domain graph, and budget domain graph. The steps for establishing a shading propagation band using buildings and roads in carbon activity units as shading sources and writing the attributes of the shading propagation band into the edge attributes of the three-domain map include: in the observation domain map, taking the carbon activity unit with the object label as a road as the center, extracting the road skeleton direction, and generating fan-shaped or strip-shaped shading propagation bands on both sides along the skeleton direction; taking the large-volume continuous carbon activity unit with the object label as a building as the impedance center, generating a ring-shaped shading propagation band around the impedance center; for carbon activity units that simultaneously receive shading effects from roads and buildings, synthesizing shading propagation bands in the order of road guidance priority followed by building continuous inhibition; projecting the synthesized shading propagation bands onto the edges of the observation domain map, process domain map, and budget domain map, and writing the visibility level, succession velocity level, and flux impedance level into the edge attributes of the three-domain map respectively; The iteration order is set as follows: observation consistency constraints first, process consistency constraints second, and budget consistency constraints last. In each iteration, observation consistency constraints are first applied to each carbon activity unit on the map of the observation domain. The texture description, edge description, and color combination description of the base image are compared with the visibility status of the candidate growth trajectory set in the corresponding time period. Growth trajectories that do not meet the visibility level requirements are deleted or the visibility level of the corresponding edge is reduced. Then, process consistency constraints are applied to the map of the process domain. The candidate growth trajectory set is screened according to the structural stability cues output by the FORMIND model. Growth trajectories that are unstable in the void structure or succession path are deleted or the succession velocity level of the corresponding edge is reduced. Finally, budget consistency constraints are applied to the map of the budget domain. The flux impedance level is locally adjusted according to the source-sink balance between adjacent carbon activity units. The edges corresponding to the candidate growth trajectories that do not meet the balance requirements are marked as pending repair. After a deletion, downgrade, or repair mark occurs in any domain's graph, a change record for that domain's graph is generated. Residual propagation is triggered through the unified index lookup table, and the change record is mapped to the relevant nodes and edges in the other two domains' graphs. The visibility level, succession rate level, and flux impedance level of the candidate growth trajectory set and the corresponding shading propagation band are adjusted accordingly. When all three types of consistency constraints are completed in this round, the next cycle begins until all carbon activity units have no repair marks in the three domains' graphs and the size of the candidate growth trajectory set no longer changes. In the three-domain graph state after the cycle ends, for each carbon activity unit, the inversion selector is invoked to make decisions according to a fixed priority rule. Priority is given to resolving residual source-sink conflicts in the budget domain graph, followed by resolving residual observation inconsistencies in the observation domain graph, and then resolving residual structural instability in the process domain graph. If multiple candidate growth trajectories simultaneously satisfy the constraints of the three-domain graph, the growth trajectory with the highest consistency score is adopted as the unique growth trajectory for that carbon activity unit. After the unique growth trajectory is determined, based on the source-sink performance of the growth trajectory in the budget domain graph and the degree of influence of the shading propagation band in the observation domain graph and process domain graph, the anthropogenic carbon emission dispersion level of the carbon activity unit is determined by looking up a table, and the unique growth trajectory number and the anthropogenic carbon emission dispersion level are written back to the record corresponding to the carbon activity unit in the unified index lookup table.

2. The method as described in claim 1, characterized in that, In step 2, for carbon activity units whose object labels are neither buildings nor roads, three types of driving information—site category, canopy layer, and void ratio—are inferred based on their texture description, edge description, and color combination description, and a multi-scenario combination is constructed. For each such carbon activity unit, the driving information and scenario combination of the carbon activity unit are input into the FORMIND model one by one, and the FORMIND model is run to obtain several time-ordered growth trajectories. All growth trajectories corresponding to the same carbon activity unit are aggregated into a candidate growth trajectory set for that carbon activity unit. For carbon activity units whose object labels are buildings or roads, no growth trajectory is generated, but the source attribute is only retained in the budget domain map. The number of the candidate growth trajectory set is linked to the unique identifier of the carbon activity unit in the unified index lookup table to ensure that the candidate growth trajectory set and the nodes of the observation domain map, process domain map, and budget domain map can be bidirectionally retrieved.

Images