Artificial forest and natural forest identification method based on multi-source remote sensing and spatio-temporal joint optimization

Abstract

The application discloses a kind of plantation and natural forest identification method and system based on multi-source remote sensing and space-time joint optimization, belong to remote sensing image processing, earth observation and forest ecological dynamic monitoring technical field.For the problems of large label noise and time sequence logic disorder in prior art, the application obtains high-purity training samples by fusing multi-stage iterative purification of CCDC mutation history and deep semantic features;Independent initial separation is carried out according to geographical division to obtain pixel level confidence;The HMM model is constructed by fusing ecological prior constraint and dynamic adaptive emission probability to carry out Viterbi decoding;Restricted spatial repair with cross-dimensional time sequence security check is performed.The application effectively suppresses time sequence false change and salt and pepper noise, realizes time logic without paradox and spatial patch regular annual high-resolution thematic mapping, and is suitable for fine classification of tropical and subtropical region forests and zero deforestation bill compliance monitoring.

Description

Technical Field

[0001] This invention belongs to the fields of remote sensing image processing, earth observation and forest ecological dynamic monitoring technology, and in particular relates to a method for identifying artificial forests and natural forests based on multi-source remote sensing and spatiotemporal joint optimization. Background Technology

[0002] In the grand context of global efforts to address climate change and biodiversity loss, accurate monitoring, reporting, and verification (MRV) of forest dynamics in time and space is crucial for achieving carbon neutrality and sustainable development goals. Tropical forests, as key regulators of the global carbon cycle and biodiversity, have received unprecedented attention. Southeast Asia, for example, possesses nearly 15% of the world's tropical forests but is also one of the regions with the highest deforestation rates globally. Unlike other regions (such as South America, primarily driven by livestock farming, and Africa, primarily by agricultural expansion), forest loss in Southeast Asia and elsewhere is largely driven by the expansion of economic plantations fueled by global commodity demand (such as oil palm and rubber). This process is characterized by the large-scale clearing of biodiversity-rich, high-carbon-storage natural forests and their conversion into intensively managed plantations. Plantations and natural forests differ fundamentally in their ecological functions, carbon storage capacity, and biodiversity carrying capacity: natural forests have complex structures, strong ecological stability, and persistent carbon sequestration, while plantations often have short rotation periods and limited species diversity. However, existing global land cover data products (such as MODIS VCF, ESRI 10m Land Cover, ESA WorldCover, etc.) include the "Forest / Tree" category, but none of them further distinguish between plantations and natural forests at the remote sensing observation scale. This mixed classification masks the hidden deforestation process of natural forests degrading into plantations, resulting in distorted carbon emission accounting and making it difficult to meet the urgent needs of refined forest management and zero deforestation legislation (such as EUDR) at the national and continental levels.

[0003] For distinguishing between plantations and natural forests and for dynamic monitoring, existing remote sensing mapping techniques mainly rely on large-scale static forest remote sensing mapping methods with medium to low spatial resolution. These techniques are primarily based on optical satellite imagery with a spatial resolution of 500 meters (e.g., MODIS data) or 30 meters (e.g., Landsat data). By extracting vegetation indices and multispectral features for a specific year (e.g., 2020), they utilize traditional machine learning algorithms such as Support Vector Machines (SVM) or Random Forests (RF) to perform large-scale classification mapping, generating static snapshot-style forest type distribution maps, such as the Global Forest Type Map (GFT2020) published by the European Union Joint Research Centre. However, the 30-meter or 500-meter resolution results in a large number of mixed pixels in the imagery. In highly fragmented landscapes caused by agricultural expansion, it is impossible to effectively capture the boundaries of small-scale economic forest plantations (usually less than 30 meters) under a small-scale farming system, leading to serious distortions in forest area calculations. In addition, due to the limitations of cloud and rain obstruction and algorithm computing power, the above methods usually only provide a single static snapshot of a specific year, which cannot reflect the continuous spatiotemporal evolution of the "game of advance and retreat" of natural forests being cut down and gradually converted into plantations, and it is difficult to capture the precise time nodes of deforestation.

[0004] To compensate for insufficient resolution, another approach is to employ spatial vector database construction technology based on multi-source data compilation. This technology does not directly perform low-level remote sensing pixel classification. Instead, it collects official forestry survey reports from various countries, concession maps of large agricultural companies, and sampling survey data from local areas. Using Geographic Information Systems (GIS), it stitches and merges these heterogeneous polygonal vector data from multiple sources to form a large-scale spatially distributed plantation forest database, such as the Spatial Distributed Plantation Forest Database (SDPT) published by the World Resources Institute (WRI). However, different countries have different legal definitions of "plantation forest" and "agricultural land" (for example, some countries define oil palm and rubber as purely agricultural crops rather than forests). This leads to abrupt differences in statistical caliber and patch density when the database crosses national borders, making it impossible to establish a unified mapping standard. More seriously, because the data source contains a large amount of historical outdated data and approximate boundaries, if the coordinate points in this database are directly extracted as training samples in the remote sensing algorithm, a large number of pixels that are actually other land features (or whose boundaries are drawn incorrectly) will be incorrectly labeled as artificial forests, introducing systematic error samples into the model, leading to severe overfitting in the subsequent classification model.

[0005] With the widespread availability of high-resolution satellite data, direct classification techniques based on single-temporal phases of high-resolution remote sensing have become a new research direction. This type of technology uses 10-meter high-resolution satellite data (such as Sentinel-2 satellite data) to extract multispectral bands and texture features, directly inputting them into deep learning semantic segmentation networks (such as U-Net and DeepLab) or machine learning models for end-to-end pixel-level direct classification, outputting a forest category distribution map for the current year. However, at 10-meter resolution, mature plantations (such as high-density rubber or oil palm forests) have extremely closed canopies, and their spectral characteristics highly overlap with natural forests. Single-temporal direct classification easily generates a large number of fragmented misclassification noise points (i.e., salt-and-pepper noise) within homogeneous forest patches, resulting in severe spatial fragmentation and disrupting the patch continuity of the real landscape. In addition, tropical regions are often cloudy and rainy, resulting in a severe lack of high-quality optical images. Combined with phenological fluctuations such as leaf fall, high-resolution classification in a single year is prone to "flickering" phenomena in multi-year series, such as "natural forest → plantation forest → natural forest," which are logically incorrect in terms of physical reality and ecological succession.

[0006] To address the temporal fluctuations caused by single-phase classification, existing technologies have further introduced temporal optimization methods such as Conditional Random Fields (CRF), Long Short-Term Memory Networks (LSTM), or standard Hidden Markov Models (HMM) in the post-classification processing stage. Taking the standard HMM as an example, it constructs a fixed "state transition probability matrix" by statistically analyzing the macro-level land use conversion frequencies of historical years, while directly using the global confusion matrix obtained by the base classifier on the test set as the global probability matrix. Figure 1The process involves using an invariant "emission probability matrix" and then employing the Viterbi decoding algorithm to find the true land cover category sequence with the highest likelihood across the entire time series, thus achieving timeline smoothing. However, deep learning time series models such as LSTM are "black box" algorithms, making it extremely difficult to directly embed explicit hard facts of forestry ecology (such as plantations almost never regressing directly back to original closed natural forests within one or two years) as absolute limiting parameter constraints into the underlying decision engine. For standard HMMs, directly applying a global static confusion matrix as the invariant emission probability for all spatial pixels completely ignores the spatial spectral heterogeneity of remote sensing images. That is, the "confidence" of pixels located at the forest edge with severe spectral mixing (high classification uncertainty) and pixels located in the forest core with extremely pure features (low classification uncertainty) should be drastically different in their impact on the classifier's output. This rigid and inflexible emission probability can lead to either over-smoothing (masking true changes) or under-smoothing (preserving spurious changes). More importantly, traditional time-series smoothing models only perform independent optimization of a single pixel on the time axis (one-dimensional vertical), completely ignoring the "homogeneous continuity of landscape patches" that should exist between adjacent pixels on the two-dimensional plane of remote sensing images. As a result, even if the time axis is smoothed, the final classification map output is still full of isolated fragment spots in the spatial dimension, and the spatial filtering process may destroy the logical consistency that has been optimized on the time axis.

[0007] In summary, current technologies for effectively and accurately acquiring the evolution characteristics of large-scale natural and planted forests at a high spatial resolution of 10 meters and a long-term, yearly frequency are still lacking. Existing technologies mainly face key technical bottlenecks that urgently need to be overcome, including high proportions of label noise in sample acquisition, temporal logic errors caused by single-phase classification, the inability of traditional time-series models to perceive pixel-level classification uncertainties, and severe fragmentation in spatiotemporal dimension processing. Therefore, there is an urgent need for a novel method for identifying planted and natural forests that can integrate multi-source high-resolution remote sensing data, automatically purify training samples, and embed ecological prior constraints within a spatiotemporal joint framework. Summary of the Invention

[0008] To address the aforementioned shortcomings of existing technologies, the present invention aims to provide a method for identifying planted and natural forests based on multi-source remote sensing and spatiotemporal joint optimization, specifically solving the following four major technical problems:

[0009] (1) To address the problem of low resolution and inconsistent standards in existing technologies leading to high label noise: This invention aims to provide a "multi-stage iterative automated sample purification mechanism". By continuously changing detection of historical records and using Mahalanobis distance chi-square statistics for edge removal, erroneous samples are removed, and a high-purity large-scale training sample base can be provided without manual labeling.

[0010] (2) In view of the problems that existing technologies are prone to temporal “pseudo-changes” and that black-box models cannot forcibly embed professional forestry laws: This invention aims to forcibly embed the minimum penalty parameter representing “irreversible ecological succession” into the state transition matrix, so as to completely cut off unreasonable and illegal succession paths from the mathematical bottom.

[0011] (3) In view of the problem that the standard HMM emission probability is extremely rigid in the existing technology: The present invention aims to provide a "dynamic adaptive emission probability model that senses classification uncertainty", which converts the confidence of the spectral features of the pixel into adjustment weights, so that the algorithm can make intelligent decisions like an expert: trust the ecological history succession trajectory when the features are ambiguous, and trust the current remote sensing spectrum when the features are certain.

[0012] (4) In view of the problems of spatiotemporal fragmentation and severe spatial fragmentation in the existing technology, the present invention aims to provide a "two-dimensional spatial context adaptive repair mechanism with ecological boundary blocking pre-verification", which, while performing spatial smoothing to eliminate salt and pepper noise, wonderfully protects the real boundary of small economic forest farms from being swallowed up and accidentally damaged through the time dimension security check lock.

[0013] To achieve the above objectives, the first technical solution of this application discloses a method for identifying planted forests and natural forests based on multi-source remote sensing and spatiotemporal joint optimization, comprising the following steps:

[0014] Input steps: Obtain a multi-source remote sensing time-series dataset of the target monitoring area within a specified year range. The dataset contains deep semantic embedding feature vectors of high-resolution remote sensing images and pixel-by-pixel surface change history parameters generated by continuous change detection and classification algorithms. Define the spatial range of candidate samples based on the dataset and a reference forest distribution base map.

[0015] Purification steps: Perform multi-stage iterative purification on the deep semantic embedding feature vector and surface mutation history parameters within the candidate sample space to output a high-purity training sample set; the multi-stage iterative purification includes initial blind selection, time-series stability logic screening based on the surface mutation history parameters, extreme outlier and edge confusion sample reduction based on Mahalanobis distance and chi-square distribution probability in the feature space, and physical grid quota balancing.

[0016] Initial classification step: Apply the high-purity training sample set to the multi-source remote sensing time series dataset, perform geographically independent annual initial classification, and output the annual preliminary hard classification results and pixel-level classification confidence distribution;

[0017] The time-series optimization steps are as follows: Based on the preliminary hard classification results and pixel-level classification confidence distribution over the years, a state transition probability matrix and a dynamic adaptive emission probability matrix that integrate ecological prior constraints are constructed, and the time-series optimized classification results are output through Viterbi decoding; the ecological prior constraints include the limit penalty parameter for preventing the reverse succession of plantations into natural forests in the short term; the dynamic adaptive emission probability matrix performs pixel-level dynamic modulation on the static baseline emission probability based on the pixel-level classification confidence distribution;

[0018] Spatial restoration steps: Perform restricted two-dimensional spatial context restoration on the temporal optimized classification results, and before modifying the spatial state of any pixel, pre-test whether the modification will trigger a violation of the ecological prior constraints in the temporal dimension. If it is triggered, the modification will be revoked, and the annual classification thematic raster data of planted forests and natural forests will be output.

[0019] Furthermore, the pixel-by-pixel surface mutation history parameters include the year of occurrence of the surface mutation t_Break and the total number of historical cumulative mutations n_Breaks; the deep semantic embedding feature vector is a 64-dimensional dense feature representation extracted from a pre-trained deep learning surface basic model.

[0020] Furthermore, the temporal stability logic screening based on surface mutation history parameters includes: for plantation forest candidate samples, it is mandatory that the year t_Break in which the surface mutation occurs is later than the preset baseline year and the subsequent state remains stable; for natural forest candidate samples, it is mandatory that the total number of historical cumulative mutations n_Breaks is less than 2 during the entire monitoring period and the year of the last disturbance is earlier than the preset baseline year.

[0021] Furthermore, the extreme outlier and marginal confusion sample reduction based on Mahalanobis distance and chi-square distribution probability in the feature space includes: calculating the Mahalanobis distance between the deep semantic embedding feature vector of the candidate sample and the mean center of its provisional class; removing outlier samples with a significance level P value less than 0.001 corresponding to the Mahalanobis distance; and randomly reducing marginal confusion samples with a P value between 0.001 and 0.01 by 70%.

[0022] Furthermore, the pixel-level dynamic modulation of the static baseline emission probability based on the pixel-level classification confidence distribution includes: for any pixel at time t, the deep semantic embedding feature vector Calculate the Mahalanobis distance between the target class j and the center of the target class j and convert it into spectral classification confidence weights. The static reference emission probability is obtained by using the nonlinear amplitude modulation index parameter β. Synthesized into dynamic adaptive emission probability elements ;in,

[0023] , This represents the distance between the current pixel feature and the standard center point of the target category;

[0024] And β is set to 2-5.

[0025] Furthermore, the limiting penalty parameter for preventing the plantation from reverting to natural forest in the short term sets the probability of the plantation transitioning to natural forest in the state transition probability matrix to be less than or equal to 1×10. -9 The probabilities of the remaining rows in the matrix are adjusted accordingly using normalization.

[0026] Furthermore, the pre-simulation of whether the modification triggers a violation of the ecological prior constraints in the temporal dimension includes: after simulating the execution of the spatial state modification, checking whether the pixel has a violation of the succession path in which the artificial forest directly reverses into the natural forest within one year. If it does, a safety brake is triggered and the modification is revoked.

[0027] Furthermore, the restricted two-dimensional spatial context repair includes: calculating the mainstream class consistency ratio C of a cell within its connected neighborhood window. p Only if the ratio is within the preset middle range of 0.25 ≤ C p Spatial assimilation is triggered when C < 0.5; p ≥0.5 or C p When <0.25, the original category of the pixel remains unchanged.

[0028] The second technical solution of this application discloses a system for identifying planted and natural forests based on multi-source remote sensing and spatiotemporal joint optimization, including:

[0029] The data input interface is configured to receive a multi-source remote sensing time-series dataset of the target monitoring area within a specified year range. The dataset includes deep semantic embedding feature vectors of high-resolution remote sensing images and pixel-by-pixel surface mutation history parameters generated by continuous change detection and classification algorithms.

[0030] The multi-stage purification module is connected to the data input interface and is configured to perform temporal stability logic filtering, Mahalanobis distance chi-square statistical cleaning and grid quota balancing based on the surface mutation history parameters and deep semantic embedding feature vectors, and output a high-purity training sample set.

[0031] The geographic zoning initial classification module is connected to the multi-stage purification module and the data input interface. It is configured to perform annual initial classification on the multi-source remote sensing time series dataset based on the high-purity training sample set, and output annual preliminary hard classification results and pixel-level classification confidence distribution.

[0032] The adaptive temporal optimization module is connected to the geographic partitioning initial classification module and is configured to construct a state transition probability matrix and a dynamic adaptive emission probability matrix that integrate ecological prior constraints based on the annual preliminary hard classification results and the pixel-level classification confidence distribution, and output the temporal optimization classification results through Viterbi decoding.

[0033] A confined space repair module, connected to the adaptive temporal optimization module, is configured to perform two-dimensional spatial context repair on the temporal optimization classification results, and to pre-simulate the compliance of the modification in the temporal dimension before performing the spatial state modification; if a violation of ecological prior constraints is triggered, the modification is revoked; and

[0034] The data output interface is connected to the confined space restoration module and configured to output annual thematic raster data of planted forests and natural forests.

[0035] And a computer-readable storage medium having a computer program stored thereon, wherein when the computer program is executed by a processor, it implements the method for identifying artificial and natural forests based on multi-source remote sensing and spatiotemporal joint optimization as described in the first technical solution.

[0036] Compared with the prior art, the present invention has the following beneficial effects:

[0037] (1) By using a multi-stage iterative sample purification mechanism, label noise is significantly reduced and the model's generalization ability is improved:

[0038] By employing a multi-stage purification method that integrates CCDC time-series stability screening, Mahalanobis distance chi-square extreme statistical cleaning, and physical grid quota balancing, this invention addresses the high proportion of label noise caused by directly using open-source coarse-resolution data as training samples in existing technologies. This invention can automatically filter out outlier mislabeled samples and decision boundary confusion samples in the spectral space without manual annotation, significantly improving the purity of the core endmembers of the training samples. This effectively overcomes the overfitting defect of the ensemble tree model caused by spectral similarity and significantly improves the generalization performance of the classification model on cross-regional and cross-year data.

[0039] (2) By constructing a dynamic adaptive emission probability model that perceives pixel-level classification uncertainty, the flexibility and accuracy of time series analysis are optimized:

[0040] By employing a technique that calculates spectral classification confidence weights based on Mahalanobis distance in the feature space and introduces a nonlinear amplitude modulation index to dynamically modulate the static reference emission probability at the pixel level, this invention achieves a flexible judgment mechanism that "trusts current remote sensing observations when features are clear and follows historical succession trajectories when features are ambiguous." This effectively suppresses temporal pseudo-changes while avoiding the real change hysteresis error caused by traditional hard smoothing, significantly improving the timeliness accuracy of long-term annual monitoring.

[0041] (3) By incorporating state transition limit constraints based on ecological prior physical constraints, temporal logic paradoxes are blocked at their source:

[0042] Because a limit penalty parameter (≤1×10) representing "planted forests cannot be reversed to natural forests in the short term" is forcibly embedded in the state transition probability matrix of the Hidden Markov Model. -9 Compared to existing technologies where black-box models cannot embed hard forestry knowledge or standard HMMs lack succession constraints, this invention cuts off the illegal state transition path from "planted forest to natural forest" from the mathematical level. This ensures that the annual classification results output by the algorithm strictly follow the physical laws of forest ecological succession, effectively eliminating the "high-frequency flickering" type of temporal pseudo-changes caused by tropical cloud and rain obstruction and phenological fluctuations, and ensuring the logical consistency of long-term monitoring data.

[0043] (4) Achieving spatiotemporal consistency through a restricted two-dimensional spatial context repair mechanism with ecological boundary blocking pre-verification:

[0044] Before performing two-dimensional spatial context repair, an ecological security interception procedure based on temporal dimension pre-simulation is introduced. Compared with the "temporal-spatial paradox" problem caused by the separation of temporal optimization and spatial filtering in the existing technology (i.e., spatial smoothing destroys the temporal optimization results), this invention effectively removes "salt and pepper noise" inside homogeneous patches, while accurately protecting the real boundaries of small economic forest farms in highly fragmented landscapes from being swallowed by excessive smoothing through a safety braking mechanism. This achieves an organic unity between the logical consistency of the temporal dimension and the regularity and coherence of the spatial dimension patches, and outputs a thematic classification product that combines high accuracy and landscape ecological rationality.

[0045] (5) Enhance adaptability to large-scale heterogeneous environments through a modular architecture with independent geographical partitions:

[0046] By employing a strategy of pre-dividing a large-scale study area into multiple computational partitions based on geographic and ecological homogeneity and independently performing supervised ensemble classification, this invention effectively mitigates the negative impact of spectral heterogeneity across climate zones and geomorphic units on classification accuracy, compared to the poor regional adaptability caused by the use of a globally unified model in existing technologies. This ensures balanced and high accuracy across the entire region, from tropical to subtropical and from plains to mountains, and significantly improves the transferability and robustness of the method in different ecological regions around the world. Attached Figure Description

[0047] To more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings used in the description of the embodiments of the present invention will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0048] Figure 1 This is a flowchart of the construction process for a multi-source remote sensing time-series dataset;

[0049] Figure 2 This is a flowchart of the sample purification process;

[0050] Figure 3 This is a schematic diagram of the construction of an adaptive hidden Markov model;

[0051] Figure 4 This is a diagram illustrating the process of finding the globally optimal path based on the Viterbi dynamic programming algorithm.

[0052] Figure 5 This is a schematic diagram of the global optimal path solution process based on the Viterbi dynamic programming algorithm, showing the transition relationship between three states (planted forest PF, natural forest NF, and non-forest land Other) in the time series from 2017 to 2024;

[0053] Figure 6 This is a schematic diagram of the spatial distribution of planted and natural forests in 2024, extracted from examples in 11 Southeast Asian countries based on this invention patent.

[0054] Figure 7 This is a comparison chart of the verification results of user accuracy (UA) and producer accuracy (PA) for three categories: plantations, natural forests and non-forest land from 2017 to 2024. UA OF PF represents user accuracy for plantations, PA OF PF represents producer accuracy for plantations, UA OF NF represents user accuracy for natural forests, PA OF NF represents producer accuracy for natural forests, UA OF Others represents user accuracy for non-forest land, and PA OF Others represents producer accuracy for non-forest land. Detailed Implementation

[0055] To make the technical problems solved, the technical solutions, and the beneficial effects of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and are not intended to limit the invention.

[0056] First, to explain the implementation scheme of this application, the following technical terms are defined:

[0057] The term "PF (Plantation Forest)" refers to planted forests (such as economic plantations like oil palm, rubber, and eucalyptus).

[0058] The term "NF (Natural Forest)" refers to a natural forest.

[0059] The term "RF (Random Forest)" refers to a random forest.

[0060] The term "HMM (Hidden Markov Model)" refers to a hidden Markov model.

[0061] The term "GEE (Google Earth Engine)" refers to the Google Earth Engine cloud computing platform.

[0062] The term "FAO (Food and Agriculture Organization)" refers to the Food and Agriculture Organization of the United Nations.

[0063] The term "CCDC (Continuous Change Detection and Classification)" refers to a continuous change detection and classification model.

[0064] The term "SDPT (Spatially Distributed Plantation Trees)" refers to spatially distributed plantation data.

[0065] The term "EUDR (European Union Deforestation Regulation)" refers to the European Union's zero-forestation legislation.

[0066] The term "OA (Overall Accuracy)" refers to overall accuracy.

[0067] The term "UA (User's Accuracy)" refers to user accuracy.

[0068] The term "PA (Producer's Accuracy)" refers to producer accuracy.

[0069] The first embodiment of this application discloses a method for identifying planted forests and natural forests based on multi-source remote sensing and spatiotemporal joint optimization, including the following steps:

[0070] Input steps: Obtain a multi-source remote sensing time-series dataset of the target monitoring area within a specified year range. The dataset contains deep semantic embedding feature vectors of high-resolution remote sensing images and pixel-by-pixel surface change history parameters generated by continuous change detection and classification algorithms. Define the spatial range of candidate samples based on the dataset and a reference forest distribution base map.

[0071] like Figure 1 As shown, the multi-source remote sensing time-series dataset obtained in the above steps includes deep semantic embedding feature vectors of high-resolution remote sensing images (multi-source high-resolution remote sensing time-series feature dataset) and pixel-by-pixel surface mutation history parameters generated by continuous change detection and classification algorithms.

[0072] The deep semantic embedding feature vector of the high-resolution remote sensing image is a 64-dimensional dense feature representation extracted by a pre-trained deep learning surface model. Specifically, it is generated by the system acquiring high-resolution remote sensing observation sequences of a specified year range (e.g., 2017-2024) of the target monitoring area, and processing them with a deep learning surface model (e.g., Google AlphaEarth large model) to produce annual cloudless mosaic image features (this processing can overcome the expression bottleneck of traditional multispectral limited bands). The system directly extracts its output dense multidimensional temporal semantic embedding feature vector (EmbeddingVector), so that each pixel with a spatial resolution of 10 meters has a 64-dimensional high-level composite feature vector (represented as vector x, with a dimension of 64), thereby capturing the comprehensive implicit features of spectrum, physical texture and annual growth phenology.

[0073] The pixel-by-pixel surface mutation history parameters are mutation detection products generated by the Continuous Change Detection and Classification (CCDC) algorithm, including the year of occurrence of surface mutation t_Break and the total number of historical cumulative mutations n_Breaks.

[0074] It should be noted that the parameter set also includes existing forest spatial distribution base maps (such as SDPT plantation vector polygons, 30-meter natural forest / plantation product NPF30), annual basic land cover category mask data (such as ESRI 10-meter land cover data), etc., which, together with the pixel-by-pixel land cover abrupt change history parameters, constitute the auxiliary parameter set. Figure 1 Furthermore, the system performs equal-area projection conversion, precise spatial registration, and unified resampling alignment to 10 meters for all the aforementioned heterogeneous data.

[0075] Purification steps: Perform multi-stage iterative purification on the deep semantic embedding feature vector and surface mutation history parameters within the candidate sample space to output a high-purity training sample set; the multi-stage iterative purification includes initial blind selection, time-series stability logic screening based on the surface mutation history parameters, reduction of extreme outliers and edge confusion samples based on Mahalanobis distance and chi-square distribution probability in the feature space, and physical grid quota balancing.

[0076] like Figure 2 The diagram shows the sample purification steps, where:

[0077] The initial blind selection specifically involves: within the reference range defined by auxiliary data, and in conjunction with the "forest / non-forest" mask of land cover products, randomly and discretely generating initial candidate sample point groups of three categories—planted forests, natural forests, and non-forest land (other land features)—within the study area at a preset sparse density (e.g., 1 point per square kilometer).

[0078] The temporal stability logical screening based on the surface mutation history parameters (CCDC temporal succession stability and logical blockade preliminary screening) includes extracting the mutation year parameter t_Break and the total number of mutations n_Breaks for each candidate pixel over the past twenty years. For "planted forest" candidate samples: it is mandatory that a significant surface mutation marking the clearing and reforestation process must have occurred after the selected reference year (e.g., 2010) (i.e., satisfying t_Break > 2010 and subsequent stability). For "natural forest" candidate samples: it is mandatory that the original surface state remains extremely stable for a long period of time, limiting the cumulative number of mutations to very low over the entire monitoring period (i.e., satisfying n_Breaks < 2), and the year of the last occasional minor disturbance must be old (earlier than 2010). Candidate points that violate the above logic are directly eliminated.

[0079] The extreme outlier and marginal confusion sample reduction based on Mahalanobis distance and chi-square distribution probability in the feature space includes: calculating the Mahalanobis distance between the deep semantic embedding feature vector of the candidate sample and the mean center of its provisional class; removing outlier samples with a significance level P-value less than 0.001 corresponding to the Mahalanobis distance; and randomly reducing marginal confusion samples with P-values ​​between 0.001 and 0.01 by 70%; specifically:

[0080] (1) For the remaining candidate samples, extract the 64-dimensional feature vector x obtained in step one, and calculate the mean center of each pixel point to its provisional class. The spatial Mahalanobis distance (DM) is calculated using the following formula: .

[0081] Where Σc is the covariance matrix of the samples in this category. Since the squared Mahalanobis distance follows a chi-square distribution, the system calculates the significance level P-value and performs double washing:

[0082] ① Absolute exclusion: Remove outlier mislabeled samples whose significance P-value falls below the minimum red line threshold (e.g., P < 0.001).

[0083] ② Reduce unstable samples: For suspected confusing samples located in the confusion transition zone (e.g., the set interval is 0.001 ≤ P < 0.01), the system forcibly applies a large proportion (e.g., a reduction and discard rate of up to 70%) of random blind box discarding actions. This step greatly compresses the fuzzy zone of the decision surface by sacrificing the sample size, forcing the remaining samples to be the core endmembers with the purest features.

[0084] The physical grid quota balancing specifically involves forced stratified quota sampling of the cleaned sample group within a preset latitude and longitude physical boundary (such as defining a large-scale grid array of 0.5 degrees x 0.5 degrees). The system forces the proportion of plantation forest, natural forest, and non-forest land samples remaining in each grid to converge to an absolutely equal 1:1:1 ratio. This operation aims to reduce spatial autocorrelation overfitting and ultimately output a high-purity training set.

[0085] Initial classification step: Apply the high-purity training sample set to the multi-source remote sensing time series dataset, perform geographically independent annual initial classification, and output the annual preliminary hard classification results and pixel-level classification confidence distribution.

[0086] In this step, considering the vastly different geographical environments and characteristics of planted and natural forests across the country and even globally, this patent employs a classification technique based on geographical zoning. The relevant steps are as follows:

[0087] (1) In order to digest the climate and landform heterogeneity in large-scale mapping, this implementation method uses the large study area to be pre-divided into multiple relatively isomorphic geoecological computing zones (e.g., divided into 10 zones) based on latitude and natural environment.

[0088] (2) Within each partition, the supervised ensemble network model for machine learning is trained independently, relying entirely on the high-quality training sample set output from step two. This scheme preferably uses a random forest classification array containing no less than 300 parallel decision trees.

[0089] (3) The training and calibrated model array is used to extrapolate all feature pixels in the target time series year by year. The system output includes: the preliminary hard classification map of the current year, the "predicted probability tensor distribution" of each category, and the "static global confusion matrix" obtained by the classifier on the validation set.

[0090] The time-series optimization steps are as follows: Based on the preliminary hard classification results and pixel-level classification confidence distribution over the years, a state transition probability matrix and a dynamic adaptive emission probability matrix that integrate ecological prior constraints are constructed, and the time-series optimized classification results are output through Viterbi decoding; the ecological prior constraints include the limit penalty parameter for preventing the reverse succession of plantations into natural forests in the short term; the dynamic adaptive emission probability matrix performs pixel-level dynamic modulation on the static baseline emission probability based on the pixel-level classification confidence distribution;

[0091] Satellite images are often obscured by clouds and rain, causing the system to classify the same area inconsistencies between this year and last year (e.g., it was woodland last year, grassland this year, and woodland again next year). This step involves establishing a "timeline optimization engine" to analyze the timeline of each pixel and correct these "flickering pseudo-changes" caused by weather or algorithmic errors. The specific technical process is as follows: Figure 3 As shown.

[0092] (1) Construct the state transition probability matrix (T matrix) of ecological physical constraints

[0093] Based on macro-level land evolution area sampling, an interannual transfer probability matrix is ​​established, comprising three states: planted forest (state 1), natural forest (state 2), and non-forest land (state 3). Subsequently, the system enforces mandatory human intervention, applying forestry physical laws: highly artificially intervened planted forests, without a long period of natural secondary succession, cannot directly revert to their original, closed natural state within one year. That is, the probability of transferring from "planted forest" to "natural forest" is mathematically suppressed at the underlying level, with the state transition parameter suppressed and locked into an extremely small, non-zero limit penalty value. .

[0094] The probabilities of the remaining pairs in the matrix are then normalized and adjusted accordingly. This measure fundamentally prevents the high-frequency logical jump pseudo-change of "natural forest → plantation forest → natural forest".

[0095] (2) Reconstructing the dynamic adaptive emission probability matrix of "perceptual pixel-level classification uncertainty" ( matrix)

[0096] Traditional models treat all pixels equally, but in reality, some pixels are located at forest edges or partially obscured by clouds, resulting in blurred features. This patent abandons the approach of applying a single static confusion matrix to the entire image. For any pixel p at time t, the 64-dimensional features... Calculate the Mahalanobis distance (DM) to the center of the target class j, and convert it into a single-pixel independent spectral classification confidence representation weight. :

[0097] ;

[0098] in, The distance represents the "distance" between the current pixel feature and the standard center point of the target category. The larger the distance, the more "non-standard" the pixel is. Adding 1 to the denominator is to prevent the denominator from being 0, thus mapping the infinite distance inversely to a weight between 0 and 1. .

[0099] Next, the static baseline emission probabilities obtained from the preliminary classification will be... By introducing a nonlinear amplitude modulation index parameter (Set to 2-5, preferably 2), synthesize dynamic adaptive emission probability elements specific to the current period of this pixel. :

[0100] ;

[0101] The molecule utilizes weights of The power represents the original fixed baseline probability. Dynamically zoom in or out. The denominator ( The summation is a standard normalization operation, ensuring that the sum of the adjusted probabilities still equals 1 (i.e., 100%). The dynamic network mechanism works as follows: if a pixel's features are extremely prominent and pure (confidence weight w is extremely large), the system fully trusts the remote sensing image from that year; if the pixel's features are blurry, it falls into a transition zone (weight w is extremely small), and the system abandons the image from that year and instead trusts the historical record from the previous year. This gives the model a "flexible judgment" ability similar to that of human experts, preventing over-smoothing.

[0102] In this embodiment, the nonlinear amplitude modulation index parameter β is set to a value of 2-5, preferably 2:

[0103] When β=0, the model degenerates into the case of traditional HMM using a global static confusion matrix, losing its adaptive ability;

[0104] When β=1, the linear case occurs, where the weights have a linear effect on the probability, and the adjustment range is relatively small.

[0105] Nonlinear amplification when β>1 (β=2 is preferred) will amplify this advantage when the weight w is large (features are pure and confidence is high), making the model more confident in the current observation; when the weight w is small (features are fuzzy and confidence is low), it will suppress this disadvantage, making the model more confident in the historical transfer prior; β=2 provides a moderate nonlinear adjustment intensity, which can effectively distinguish between high and low confidence pixels without being overly sensitive.

[0106] When β > 5, it is meaningless.

[0107] (3) Global optimal path solution based on Viterbi dynamic programming algorithm

[0108] Suppose we want to monitor the actual land cover type of a piece of land over eight consecutive years (2017-2024). Since there are three possible states each year (planted forest, natural forest, non-forest land), theoretically, every pixel must contain land cover. A complete historical evolution path. The introduction of the Viterbi algorithm is precisely to efficiently and accurately find the most probable "correct" path from a massive number of possibilities. The Viterbi algorithm does not employ exhaustive search, but rather a strategy of "deriving the global optimum from the local optimum," such as... Figure 4 As shown: It specifically includes the following three steps:

[0109] ① Year-by-year local optimization and memory: When the algorithm extrapolates to a certain state in year t (e.g., assessing the probability that the year is a "planted forest"), it does not look back at earlier years, but only multiplies the scores of the three states in year t-1 to reach this point by the corresponding state transition moment T and the emission probability matrix of that year. After calculation, the algorithm retains only the access path with the highest score and decisively discards the other poorer paths.

[0110] ②Physical pruning: In the aforementioned "year-by-year competition," the present invention sets forth the following... The penalty parameter plays a crucial role. If a path attempts to traverse "from last year's planted forest to this year's natural forest," the system inputs this minimum value during multiplication, and the total score of this violating path instantly drops to zero. Therefore, it is directly eliminated during the year's "local optimization" phase (i.e., algorithmic "pruning"), and cannot qualify for the next year's iteration.

[0111] ③ Backtracking Search: When the simulation progresses to the final year (e.g., 2024), each of the three final state nodes—planted forest, natural forest, and non-forest land—holds a tightly bound "optimal historical long-term line" leading to its own outcome. At this point, the system only needs to compare the cumulative total scores of these three final nodes, select the one with the highest score, and then trace back to accurately extract the most reasonable and complete trajectory of change over these eight years, free from any physical breaks or ecological violations.

[0112] Spatial restoration steps: Perform restricted two-dimensional spatial context restoration on the temporal optimized classification results, and before modifying the spatial state of any pixel, pre-test whether the modification will trigger a violation of the ecological prior constraints in the temporal dimension. If it is triggered, the modification will be revoked, and the annual classification thematic raster data of planted forests and natural forests will be output.

[0113] HMM only addresses jumps along the time axis. To clean up isolated "salt and pepper" fragments remaining in the two-dimensional planar space, a confined space repair strategy is set.

[0114] (1) After time-series clearing, calculate the coverage consistency ratio Cp of the majority of land cover categories within each pixel p in its 8 connected or extended 24 connected topological physical window Np (i.e., the number of mainstream category pixels in the neighborhood divided by the total number of pixels in the neighborhood).

[0115] (2) Multi-branch confined space repair and reset rules: If C p ≥0.5 (located in a broad, homogeneous, safe core area) or C p <0.25 (in an extremely heterogeneous, interlocking, fragmented zone), the system is immune to smoothing intervention, forcibly maintaining the status quo to prevent the true boundaries of micro-farms and plantations from being blindly erased.

[0116] If 0.25≤C p <0.5 (characterized by a typical abrupt, discrete cluster of noise spots), triggering a spatial assimilation mechanism, which proposes to reset the identity of the weak pixel in the center and replace it with the dominant category in the surrounding area, as shown in Table 1.

[0117] Table 1

[0118] .

[0119] (3) Execute the pre-simulated ecological security check and interception procedure:

[0120] Before officially executing the physical overwrite replacement command, the system will perform a preliminary simulation in the background to determine whether the two-dimensional modification of the space will cause the pixel to trigger a violation of the red line set in step four, namely "artificial forests are prohibited from being converted into natural forests within one year," on the "one-dimensional timeline." If the simulation touches the red line of historical succession violations, the system will immediately trigger a red safety brake, forcibly canceling the current spatial smoothing command and upholding the original legal category of the pixel.

[0121] Finally, through the above five steps, the output results are automatically and seamlessly stitched together to produce a high-resolution, continuous thematic raster dataset of plantations and natural forests that combines logical inconsistencies in time history with regular and coherent spatial patch patterns.

[0122] The second embodiment of this application discloses a system for identifying planted and natural forests based on multi-source remote sensing and spatiotemporal joint optimization, including:

[0123] The data input interface is configured to receive a multi-source remote sensing time-series dataset of the target monitoring area within a specified year range. The dataset includes deep semantic embedding feature vectors of high-resolution remote sensing images and pixel-by-pixel surface mutation history parameters (mutation year t_Break, mutation number n_Breaks) generated by a continuous change detection and classification algorithm.

[0124] A multi-stage purification module, connected to the data input interface, is configured to perform temporal stability logic filtering (e.g., if there have been fewer than 2 land cover type changes since 2010: t_Break>2010 and n_Breaks<2), Mahalanobis distance chi-square statistical cleaning (P<0.001 is removed, 0.001≤P<0.01 is reduced by 70%), and grid quota balancing (1:1:1 quota within a 0.5°×0.5° grid) based on the land cover mutation history parameters and deep semantic embedding feature vectors, and output a high-purity training sample set;

[0125] The geographic partitioning module is connected to the multi-stage purification module and the data input interface. It is configured to perform annual initial classification on the multi-source remote sensing time series dataset based on the high-purity training sample set (preferably 10 partitions, and the classifier is preferably a random forest or LightGBM with 300 trees), and output the annual preliminary hard classification results and pixel-level classification confidence distribution.

[0126] The adaptive temporal optimization module, connected to the geographic partitioning initial classification module, is configured to construct a state transition probability matrix and a dynamic adaptive emission probability matrix that integrate ecological prior constraints based on the annual preliminary hard classification results and pixel-level classification confidence distribution. The state transition probability matrix and dynamic adaptive emission probability matrix (β=2) are used to output the time-optimized classification results via Viterbi decoding;

[0127] A confined space restoration module, connected to the adaptive temporal optimization module, is configured to perform two-dimensional spatial context restoration based on 8 / 24 connected neighborhoods (threshold 0.25≤Cp<0.5) on the temporal optimization classification results, and to pre-simulate the compliance of the modification in the temporal dimension (whether it triggers a violation of converting artificial forests to natural forests within one year) before performing the spatial state modification; if a violation of ecological prior constraints is triggered, the modification is revoked; and

[0128] The data output interface is connected to the confined space restoration module and configured to output annual thematic raster data of planted forests and natural forests.

[0129] Example 1

[0130] Based on the aforementioned multi-stage iterative sample purification mechanism, and compared to the high-proportion label noise problem caused by directly using open-source coarse-resolution data as training samples in existing technologies, this invention covers the Southeast Asian region, including 11 countries: Myanmar, Thailand, Laos, Vietnam, Cambodia, Malaysia, Indonesia, the Philippines, Brunei, East Timor, and Singapore. This region has complex terrain and diverse climates, possesses abundant tropical rainforest resources, and is also one of the regions with the highest deforestation rates globally, exhibiting typical highly fragmented and heterogeneous landscape characteristics, making it an ideal test area for verifying the effectiveness of the method of this invention. In the field measurements in Southeast Asia (see the 2024 distribution map of planted and natural forests in Southeast Asia),... Figure 6 As shown in the figure, the original candidate samples were subjected to CCDC time-series stability screening, Mahalanobis distance chi-square statistical cleaning (100% removal of extreme outliers at significance level P<0.001, and 70% random reduction of marginal confusion samples at P∈[0.001,0.01]), and 1:1:1 quota balancing of a 0.5°×0.5° physical grid, resulting in approximately 600,000 high-purity core endmember samples. Experiments show (experimental results are shown in the figure)... Figure 7 The comparison of user accuracy (UA) and producer accuracy (PA) between the initial RF classification and the HMM-optimized data from 8 periods of 2017-2024 is shown. This sample purification strategy effectively overcomes the overfitting defect of the ensemble tree model caused by spectral similarity.

[0131] Based on the aforementioned dynamic adaptive emission probability model and Viterbi decoding optimization (e.g., for perceptual pixel-level classification uncertainty), Figure 5 As shown, in independent validation of the present invention using eight periods of data from Southeast Asia spanning 2017-2024 (approximately 120,000 retained samples), the initial average overall classification accuracy was 94.9%, which improved to 98.0% (an improvement of approximately 3.1 percentage points) after HMM time-series optimization. Specifically, the accuracy for plantation forest users improved from 91.6% to 96.5%–97.8%, and the accuracy for natural forest users improved from 94.0% to approximately 97.5%. Experimental results demonstrate that this dynamic adaptive emission probability model effectively achieves a flexible judgment mechanism of "trusting current remote sensing observations when features are clear, and obeying historical succession trajectories when features are ambiguous." This effectively suppresses spurious temporal changes while avoiding the hysteresis error of real changes caused by traditional hard smoothing processing.

[0132] In a comprehensive test covering 11 Southeast Asian countries, the study area was divided into 10 geoecological zones and a random forest model was trained independently in the case study. This effectively mitigated the spectral heterogeneity across climate zones and geomorphic units. The output data showed a correlation coefficient of 0.99 with the natural forest area statistics from FAO and a correlation coefficient of 0.69 with the plantation forest area.

[0133] The above-described embodiments are only used to illustrate the technical solutions of the present invention, and are not intended to limit it. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention, and should all be included within the protection scope of the present invention.

Claims

1. A method for identifying planted forests and natural forests based on multi-source remote sensing and spatiotemporal joint optimization, characterized in that, Includes the following steps: Input steps: Obtain a multi-source remote sensing time-series dataset of the target monitoring area within a specified year range. The dataset contains deep semantic embedding feature vectors of high-resolution remote sensing images and pixel-by-pixel surface change history parameters generated by continuous change detection and classification algorithms. Define the spatial range of candidate samples based on the dataset and a reference forest distribution base map. Purification steps: Within the candidate sample space, perform multi-stage iterative purification on the deep semantic embedding feature vector and surface mutation history parameters to output a high-purity training sample set; the multi-stage iterative purification includes initial blind selection, time-series stability logic screening based on the surface mutation history parameters, reduction of extreme outliers and marginal confusion samples based on Mahalanobis distance and chi-square distribution probability in the feature space, and physical grid quota balancing. The physical grid quota balancing is to perform forced stratified quota sampling on the cleaned sample group within the preset latitude and longitude physical boundaries, so that the proportion of artificial forest, natural forest and non-forest land samples remaining in each grid converges to an absolutely equal 1:1:1 situation; Initial classification step: Apply the high-purity training sample set to the multi-source remote sensing time series dataset, perform geographically independent annual initial classification, and output the annual preliminary hard classification results and pixel-level classification confidence distribution; Temporal optimization steps: Based on the preliminary hard classification results and pixel-level classification confidence distribution over the years, construct a state transition probability matrix and a dynamic adaptive emission probability matrix that integrate ecological prior constraints, and output the temporal optimization classification results through Viterbi decoding; the ecological prior constraints include the limit penalty parameter for preventing the reverse succession of plantations into natural forests in the short term; The dynamic adaptive emission probability matrix performs pixel-level dynamic modulation on the static baseline emission probability based on the pixel-level classification confidence distribution, specifically including: For any pixel at time t, the deep semantic embedding feature vector Calculate the Mahalanobis distance between the target class j and the center of the target class j, and convert it into spectral classification confidence weights. The static reference emission probability is obtained by using the nonlinear amplitude modulation index parameter β. Synthesized into dynamic adaptive emission probability elements ;in, , This represents the distance between the current pixel feature and the standard center point of the target category; And β is set to 2-5; Spatial restoration steps: Perform restricted two-dimensional spatial context restoration on the temporal optimized classification results, and before modifying the spatial state of any pixel, pre-test whether the modification will trigger a violation of the ecological prior constraints in the temporal dimension. If it is triggered, the modification will be revoked, and the annual classification thematic raster data of planted forests and natural forests will be output.

2. The method according to claim 1, characterized in that, The pixel-by-pixel surface mutation history parameters include the year of occurrence of the surface mutation t_Break and the total number of historical mutations n_Breaks; the deep semantic embedding feature vector is a 64-dimensional dense feature representation extracted from a pre-trained deep learning surface basic model.

3. The method according to claim 2, characterized in that, The temporal stability logic screening based on the surface mutation history parameters includes: for plantation forest candidate samples, it is mandatory that the year t_Break in which the surface mutation occurs is later than the preset baseline year and the subsequent state remains stable; for natural forest candidate samples, it is mandatory that the total number of historical cumulative mutations n_Breaks is less than 2 during the entire monitoring period and the year of the last disturbance is earlier than the preset baseline year.

4. The method according to claim 1, characterized in that, The extreme outlier and marginal confusion sample reduction based on Mahalanobis distance and chi-square distribution probability in feature space includes: calculating the Mahalanobis distance between the deep semantic embedding feature vector of the candidate sample and the mean center of its provisional class; removing outlier samples with a significance level P value less than 0.001 corresponding to the Mahalanobis distance; and randomly reducing marginal confusion samples with a P value between 0.001 and 0.01 by 70%.

5. The method according to claim 1, characterized in that, The limiting penalty parameter for the short-term reverse succession of plantations to natural forests sets the probability of plantation transitioning to natural forest in the state transition probability matrix to be less than or equal to 1×10. -9 The probabilities of the remaining rows in the matrix are adjusted accordingly using normalization.

6. The method according to claim 1, characterized in that, The pre-simulation of whether the modification triggers a violation of the ecological prior constraints in the temporal dimension includes: after simulating the execution of the spatial state modification, checking whether the pixel has a violation of the succession path in which the artificial forest directly reverses into the natural forest within one year. If it does, a safety brake is triggered and the modification is revoked.

7. The method according to claim 1, characterized in that, The restricted two-dimensional spatial context repair includes: calculating the mainstream class consistency ratio C of a cell within its connected neighborhood window. p Only if the ratio is within the preset middle range of 0.25 ≤ C p Spatial assimilation is triggered when C < 0.5; p ≥0.5 or C p When <0.25, the original category of the pixel remains unchanged.

8. A system for identifying planted and natural forests based on multi-source remote sensing and spatiotemporal joint optimization, used to implement the method for identifying planted and natural forests based on multi-source remote sensing and spatiotemporal joint optimization as described in claim 1, characterized in that... include: The data input interface is configured to receive a multi-source remote sensing time-series dataset of the target monitoring area within a specified year range. The dataset includes deep semantic embedding feature vectors of high-resolution remote sensing images and pixel-by-pixel surface mutation history parameters generated by continuous change detection and classification algorithms. The multi-stage purification module is connected to the data input interface and is configured to perform temporal stability logic filtering, Mahalanobis distance chi-square statistical cleaning, and physical grid quota balancing based on the surface mutation history parameters and deep semantic embedding feature vectors, and output a high-purity training sample set. The geographic zoning initial classification module is connected to the multi-stage purification module and the data input interface. It is configured to perform annual initial classification on the multi-source remote sensing time series dataset based on the high-purity training sample set, and output annual preliminary hard classification results and pixel-level classification confidence distribution. The adaptive temporal optimization module is connected to the geographic partitioning initial classification module and is configured to construct a state transition probability matrix and a dynamic adaptive emission probability matrix that integrate ecological prior constraints based on the annual preliminary hard classification results and the pixel-level classification confidence distribution, and output the temporal optimization classification results through Viterbi decoding. The confined space repair module is connected to the adaptive temporal optimization module and is configured to perform two-dimensional spatial context repair on the temporal optimization classification result. Before performing spatial state modification, it rehearses the compliance of the modification in the temporal dimension. If the violation of ecological prior constraints is triggered, the modification is revoked. as well as The data output interface is connected to the confined space restoration module and configured to output annual thematic raster data of planted forests and natural forests.

9. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by the processor, it implements the method for identifying artificial and natural forests based on multi-source remote sensing and spatiotemporal joint optimization as described in claim 1.

Images