A method for remote sensing mapping of plantations by fusing TVSSI indexes

Abstract

The application discloses a kind of artificial forest remote sensing mapping methods of fusing TVSSI index, rely on GEE platform pretreatment Landsat time series image, fit NDVI time series and calculate SAVI, after building vegetation soil collaborative index NVSSI, get hyperbolic tangent vegetation soil collaborative index TVSSI by Tanh transformation;Using continuous change detection and classification algorithm CCDC generates artificial forest and natural forest sample;Through simple non-iterative clustering SNIC superpixel segmentation realizes object-oriented feature extraction;Fusion spectrum, time series, topography and TVSSI, combined with random forest algorithm RF completes artificial forest classification;Carries out accuracy verification and dataset comparison, obtains artificial forest spatial distribution result.The application can effectively overcome soil background interference and spectral saturation problem, significantly improve the distinguishing ability of artificial forest and natural forest, improve the accuracy and stability of artificial forest remote sensing mapping under complex terrain conditions.

Description

Technical Field

[0001] This invention relates to the fields of dynamic monitoring of forest resources, remote sensing information extraction, and geospatial data analysis. In particular, it relates to a method for high-precision remote sensing mapping of plantations based on hyperbolic tangent transform to construct the TVSSI vegetation-soil synergy index, and relying on the Google Earth Engine (GEE) cloud computing platform to integrate Continuous Change Detection and Classification (CCDC), Simple Non-Iterative Clustering (SNIC), and Random Forest (RF) algorithms. This method can be widely used in the verification of the effectiveness of shelterbelts, regional ecological engineering monitoring, forest land resource surveys, and carbon storage assessments. Background Technology

[0002] As a key component of terrestrial ecosystems, planted forests play an irreplaceable ecological role in water conservation, soil and water conservation, carbon sequestration and oxygen release, and regional soil erosion control. Accurately distinguishing the spatial distribution of planted forests from natural forests and understanding the annual distribution status of forest land are prerequisites for scientifically evaluating afforestation effectiveness, optimizing forest land management policies, and quantifying regional ecological benefits.

[0003] Traditional methods of spatial surveying of plantations mainly rely on field surveys and remote sensing visual interpretation. These methods require a large amount of manpower and resources for on-site sampling, resulting in long field operation cycles, high survey costs, and difficulty in implementing large-scale, comprehensive surveys. Furthermore, manual visual interpretation is affected by the experience of the interpreters, leading to significant subjective errors and making it impossible to achieve routine, dynamic monitoring year after year.

[0004] Currently, mainstream remote sensing identification methods generally rely on vegetation indices such as the Normalized Differential Vegetation Index (NDVI) and the Soil Adjusted Vegetation Index (SAVI) for forest land classification. However, these indices have two key drawbacks when identifying and extracting vegetation from plantations: First, areas with high vegetation cover exhibit significant spectral saturation, causing the index values ​​to flatten out and making it difficult to capture subtle spectral differences between well-ordered plantations and disordered natural forests. Second, when the proportion of young plantations, forest gaps, and bare soil in gullies is high, soil reflection signals severely interfere with the true value of the vegetation index, resulting in large-scale misclassification and omissions. In addition, existing technologies mostly rely on manual selection of classification samples, which involves a large sampling workload and makes it difficult to ensure sample representativeness. Conventional single-pixel classification is prone to salt-and-pepper noise, resulting in poor final mapping accuracy and interannual stability in complex geomorphic environments with broken gullies and dramatic topographic relief. Summary of the Invention

[0005] This invention provides a remote sensing mapping method for plantations that integrates the Hyperbolic Tangent Vegetation Soil Synergy Index (TVSSI). The aim is to overcome existing shortcomings such as spectral saturation of vegetation indices, severe interference from bare soil background, and low efficiency of manual sample selection. Utilizing GEE cloud processing of massive Landsat remote sensing data, a comprehensive technical solution is employed, including harmonic fitting, hierarchical index construction, automatic temporal perturbation sampling, superpixel object-oriented processing, multi-feature fusion machine learning classification, multi-dimensional index evaluation, and joint verification with field images. This achieves accurate spatial mapping of plantations at a 30m spatial resolution.

[0006] Therefore, the present invention adopts the following technical solution: A remote sensing mapping method for plantations that integrates the TVSSI index, the method comprising the following steps in sequence: S1: Remote sensing image acquisition and preprocessing to obtain a well-organized and usable remote sensing dataset of the study area; S2: Based on the preprocessed Landsat multiband surface reflectance data, perform NDVI time-series harmonic fitting (Harmonic Analysis, HA) and SAVI calculation. S3: Based on NDVI harmonic analysis fitting, NDVI_fitted and SAVI are used to construct the Normalized Vegetation Soil Synergy Index (NVSSI). S4: Perform hyperbolic tangent transform (Tangent, Tan) on NVSSI to construct TVSSI; S5: Automatically generate training samples for planted forests and natural forests based on the CCDC temporal disturbance detection algorithm; S6: Build an independent validation sample dataset based on high-resolution imagery and field measurement data; S7: Use the SNIC algorithm to perform superpixel segmentation and complete object-oriented feature extraction; S8: Integrates multiple indicators such as spectral, temporal, topographic, and TVSSI to construct a multi-dimensional fusion feature set; S9: After the samples are stratified, the random forest algorithm is used to classify and identify planted forests and natural forests. S10: Conduct multi-parameter quantitative evaluation from two dimensions: classification accuracy and index discrimination performance; S11: Generate 30m resolution spatial raster data of plantations and complete the resulting mapping; S12: Compare and verify the mapping results by combining field photos, high-resolution satellite imagery, and existing datasets.

[0007] Specifically, S1 involves: using the GEE geographic big data cloud computing platform to select and call all available Landsat 9 Collection 2 Level-2 surface reflectance images for 2025, with the selection criterion being that the cloud cover of a single image is controlled at ≤20%; using the JavaScript programming language to batch process the original remote sensing images in sequence, performing atmospheric correction processing, cloud masking processing, cloud shadow masking processing, removal of invalid snow and ice pixels, annual temporal median synthesis, and vector clipping of the administrative boundary of Dingxi City. After clipping, only the effective remote sensing data of the study area is retained.

[0008] S2 specifically involves: calculating NDVI based on the preprocessed Landsat multi-band reflectivity data; in addition, using a harmonic analysis model to fit and denoise the annual continuous NDVI time series to obtain NDVI_fitted time series data that eliminates cloud and fog interference and random noise interference; and simultaneously calculating the soil-regulated vegetation index (SAVI).

[0009] The formulas for calculating NDVI, NDVI_fitted, and SAVI are as follows: NDVI = (NIR - RED) / (NIR + RED); SAVI= (NIR−RED)×(1+L) / (NIR+RED+L); NDVI_fitted (t)=a+b×t+∑(k=1~n)[α k cos (2πkt / T)+β k sin (2πkt / T)] In the formula: NIR is the near-infrared band, RED is the red band, t represents the number of days in the time series, T is taken as the annual phenological cycle of vegetation (365 days), a and b are the overall trend fitting coefficients, and α k β k is the harmonic oscillation coefficient, n is the harmonic fitting order, and the soil adjustment coefficient L is fixed at 0.5; Specifically, S3 involves: using the ratio of the harmonic-fitted NDVI_fitted sequence to the soil-corrected SAVI data to construct the intermediate synergistic index NVSSI, calculated as: NVSSI = NDVI_fitted / SAVI. This amplifies the difference in the original spectra of plantations and natural forests by leveraging the inverse response law between the two types of indices.

[0010] S4 specifically involves performing a Tanh hyperbolic tangent nonlinear transform on the obtained NVSSI raster data to construct the final TVSSI index. Leveraging the S-curve characteristics of the Tanh function, the characteristic difference is sensitively amplified within the intermediate range of the index, while convergent extreme noise is compressed at the maximum and minimum points of the index. This effectively overcomes the spectral saturation defects of traditional linear vegetation indices under high vegetation cover conditions. The calculation formula is as follows: TVSSI = Tanh (NVSSI) = (e^NVSSI − e^(−NVSSI)) / (e^NVSSI + e^(−NVSSI)).

[0011] Specifically, S5 involves: employing the Continuous Change Detection and Classification (CCDC) algorithm, selecting time-series reflectance data from six bands in Landsat: blue, green, red, near-infrared, shortwave infrared 1, and shortwave infrared 2. A robust iterative weighted least squares method is used to perform time-series regression fitting pixel-by-pixel, automatically identifying abrupt breakpoints in time-series disturbances for each pixel. The multi-year disturbance frequency (FD) of a single pixel is statistically analyzed, and pixel regions with FD ≥ 3 are designated as training samples for planted forests. Furthermore, regions with FD = 0 and unchanged forest attributes in the China Land Cover Dataset (CLCD) are designated as training samples for natural forests, thus achieving automated extraction of samples from the entire region.

[0012] Specifically, S6 involves: using high-resolution historical images from Google Earth Pro, photos of locations collected through field surveys, and field investigation records, manually setting up evenly distributed sampling points, and separately constructing two independent verification sample sets: artificial forests and natural forests. These samples do not participate in model training but are only used for subsequent classification accuracy calculations and TVSSI index discrimination capability verification.

[0013] Specifically, S7 involves running the SNIC superpixel segmentation algorithm in the GEE cloud computing environment. The algorithm parameters are set as follows: seed spacing 15, compactness coefficient 4, and neighborhood search size 128. The generated superpixel objects are used as the basic classification units. The average spectral features and spatial texture features inside the objects are extracted in a unified manner to eliminate salt-and-pepper noise caused by single pixel classification and improve the spatial integrity of the classified patches.

[0014] Specifically, S8 involves integrating multi-source information to construct a full-dimensional fusion feature set. The features include: Landsat original B1~B7 single-band reflection features, NDVI_fitted, SAVI, TVSSI index features, NDVI time-series derived features, and DEM topographic features. The NDVI time-series derived features include seven types of indicators: annual maximum value, annual minimum value, coefficient of variation, variance, amplitude, phase, and fitting residual. The DEM topographic features include three topographic parameters: elevation, slope, and aspect.

[0015] All training samples obtained from S5 are randomly split in a 7:3 ratio. 70% of the samples are used as the training subset for fitting the random forest model, and the remaining 30% of the samples are used as the internal test subset. The random forest classifier is set with 100 decision trees. Based on the constructed multidimensional feature set, binary classification of planted forests and natural forests is performed, and the classification raster results are output.

[0016] S10 specifically involves: on the one hand, constructing a confusion matrix based on validation samples, and calculating the overall accuracy OA and Kappa coefficient to evaluate the overall classification effect; on the other hand, selecting four quantitative indicators—inter-class difference (Cohen's d), comprehensive discrimination ability (FisherRatio), difference significance (T Statistic), and overlap coefficient—to horizontally compare the advantages and disadvantages of the four indices NDVI_fitted, SAVI, NVSSI, and TVSSI in distinguishing between plantations and natural forests, and quantitatively verifying the anti-interference and discrimination performance of the TVSSI index.

[0017] Specifically, S11 involves outputting 30m spatial resolution raster data of planted forests in Dingxi City in 2025 based on the optimal classification model, forming a standardized spatial distribution map of planted forests.

[0018] S12 specifically involves randomly selecting multiple typical verification plots in different landforms and forest stand types within the study area, and verifying the reliability of the final results from the spatial distribution by comparing them with field photos taken at the corresponding locations, concurrent high-resolution satellite imagery, and the existing 2021 dataset, Global Natural and Planted Forests Mapping (GNPF) and Spatial Database of Planted Trees (SPDT).

[0019] Technical principle of the invention: This invention addresses a series of problems with traditional vegetation indices such as NDVI and SAVI, including susceptibility to bare soil background interference, spectral saturation in high-vegetation areas, low efficiency of manual sampling, salt-and-pepper noise in single-pixel classification, and insufficient classification accuracy in complex terrain. It employs a multi-level technology combination scheme to optimize the entire process from mechanism to application: First, a harmonic analysis model is used to fit and denoise the original NDVI time series, resulting in NDVI_fitted, which removes random noise from clouds, atmospheric disturbances, and short-term bare soil patches, preserving the stable phenological and spectral characteristics of forest land. This is combined with the SAVI index, which weakens soil reflection signals, to preemptively reduce background interference from forest clearings, young forests, and bare gullies. Then, a transition index, NVSSI, is constructed through the ratio calculation of the two indices. Based on the different spectral response patterns of the two indices to plantations and natural forests, this further amplifies the inherent differences in the characteristics of the two types of forest land and continuously suppresses soil noise. Finally, hyperbolic tangent Tanh is applied to NVSSI. The TVSSI index is obtained through nonlinear transformation. Leveraging the S-curve characteristic of this function, subtle differences in forest stands are amplified in the middle range and outliers are converged in the extreme range, fundamentally solving the spectral saturation problem of traditional linear vegetation indices and enhancing land cover differentiation capabilities. In the sample acquisition stage, the CCDC continuous change detection and classification algorithm is employed. Based on multi-band temporal reflectance data, pixel-by-pixel fitting is used to identify temporal disturbance breakpoints and count disturbance frequencies. Training samples are automatically divided into plantation and natural forests based on disturbance frequencies, overcoming the drawbacks of large workload, strong subjectivity, and insufficient sample representativeness of manual sampling, achieving fully automated and efficient sample extraction. In the image segmentation and feature extraction stage, the SNIC simple non-iterative clustering algorithm is used for superpixel segmentation. Consecutive superpixel objects replace scattered single pixels as classification units, simultaneously extracting spectral and texture features within the units. This effectively eliminates salt-and-pepper noise caused by single-pixel classification, improving the spatial integrity and visual effect of the classified patches. Finally, a fusion feature set is constructed by integrating multi-dimensional features such as spectral, temporal, topographic, and TVSSI data. The RF algorithm is introduced to build a classification model to complete the binary classification of planted forests and natural forests. The synergistic effect of multiple technologies greatly improves the automation, anti-interference ability and classification accuracy of the whole method, which can efficiently support the routine annual survey of regional planted forests and the quantitative evaluation of the effectiveness of ecological engineering.

[0020] The beneficial effects of this invention are as follows: This invention systematically addresses soil background noise interference from bare soil and gully clearings in forests by progressively constructing the NVSSI transition index and the novel TVSSI index, combined with harmonic time-series fitting noise reduction, soil disturbance correction, and hyperbolic tangent nonlinear transformation processing. This significantly widens the phenological spectral differences between well-ordered plantation stands and disordered natural forest stands, enhancing the ability to distinguish and identify plantation stands. Furthermore, this invention introduces CCDC continuous change detection and classification time-series disturbance detection algorithms, automatically identifying forest disturbance breakpoints based on long-term remote sensing images and automatically generating data based on the frequency of disturbances over many years. The training samples for planted and natural forests overcome the significant drawbacks of traditional manual visual sampling methods, such as massive workload, long sampling cycles, strong human subjectivity, and insufficient spatial representativeness, greatly reducing the cost of manual interpretation and field sampling. Simultaneously, SNIC is used for superpixel segmentation, employing contiguous objects as the basic classification unit to achieve object-oriented feature extraction, effectively eliminating salt-and-pepper noise commonly found in single-pixel classification, optimizing the spatial integrity of the classification results, and then integrating spectral, temporal phenological, topographic factors, and TVSSI multidimensional features to collaboratively build a random forest classification model for plantation identification. This entire technical solution can efficiently support various forestry operations in complex terrain areas, including routine annual surveys of plantation resources, verification of the effectiveness of shelterbelt construction, dynamic monitoring of regional ecological restoration projects, and quantitative assessment of ecological benefits. Attached Figure Description

[0021] Figure 1 This invention provides a remote sensing mapping process for plantations that integrates the TVSSI index. Figure 2 This is a Landsat-9 false-color image from an embodiment of the present invention; Figure 3 The embodiment of this invention is a Landsat-9 true-color image; Figure 4 This is an altitude distribution map for an embodiment of the present invention; Figure 5 This is a slope aspect distribution diagram for an embodiment of the present invention; Figure 6 This is a slope distribution map for an embodiment of the present invention; Figure 7 This is a SAVI spectral index distribution diagram as an embodiment of the present invention; Figure 8 This is a spatial distribution diagram of NDVI_fitted after harmonic fitting in an embodiment of the present invention; Figure 9 The embodiment of this invention is the distribution of the annual average NDVI (NDVI_amp); Figure 10 This is a spatial distribution diagram of NDVI dispersion (NDVI_diss) in an embodiment of the present invention. Figure 11This is a spatial distribution diagram of the NDVI autocorrelation feature NDVI_asm in an embodiment of the present invention; Figure 12 This is a spatial distribution diagram of the NDVI contrast feature NDVI_contrast in an embodiment of the present invention; Figure 13 This is a spatial distribution map of NDVI_idm, an NDVI information metric, as described in this embodiment of the invention. Figure 14 This is a spatial distribution diagram of the annual maximum NDVI value, NDVI_max, according to an embodiment of the present invention. Figure 15 This is a spatial distribution diagram of the annual minimum NDVI value, NDVI_min, according to an embodiment of the present invention. Figure 16 This is a grid map of the NDVI edge features of the Canny operator in an embodiment of the present invention. Figure 17 This is a spatial distribution diagram of the mean square error NDVI_mse for NDVI fitting in an embodiment of the present invention. Figure 18 This is a spatial distribution diagram of the root mean square NDVI_rmse residuals in an embodiment of the present invention. Figure 19 This is a spatial distribution map of the NDVI phenological phase NDVI_pha according to an embodiment of the present invention; Figure 20 Image showing the results of SNIC superpixel segmentation in the region; Figure 21 The pixel perturbation frequency FD spatial distribution map output by the CCDC algorithm; Figure 22 CCDC time series fitting curves for plantations and natural forests under different FD conditions; Figure 23 Spatial distribution map of plantation forests and natural training samples after screening by the CCDC algorithm; Figure 24 A raster distribution map of the TVSSI index; Figure 25 Separation test charts for different indicators; Figure 26 A graph for evaluating TVSSI separation capability indicators; Figure 27 This is a verification map comparing the distribution of planted forests with existing datasets and high-resolution images. Figure 28 This is a statistical chart showing the accuracy of each indicator. Detailed Implementation

[0022] The technical solution of the present invention will be further described below with reference to the accompanying drawings and specific embodiments: S1: Relying on the GEE cloud computing platform, Landsat 9 Collection 2 Level 2 surface reflection remote sensing products from 2025 were selected. The screening threshold was set to cloud cover ≤ 20%. JavaScript scripts were used to perform preprocessing operations such as atmospheric correction, cloud body and cloud shadow removal, removal of invalid snow and ice pixels, annual time series median synthesis, and administrative boundary vector clipping to obtain a regular and usable remote sensing dataset of the study area.

[0023] S2: The NDVI index is calculated based on the preprocessed Landsat multiband surface reflectance data. The Harmonic Analysis (HA) algorithm is used to fit and denoise the annual continuous NDVI time series to remove random noise caused by clouds, fog, and atmospheric disturbances, generating stationary NDVI_fitted time series data. SAVI is calculated simultaneously to reduce soil spectral interference caused by bare gullies, forest clearings, and bare ground in young forests. The NDVI_fitted and SAVI calculation results are used as the basic data for subsequent index construction.

[0024] S3: By using the ratio of harmonic-optimized NDVI_fitted to soil-corrected SAVI, an intermediate transition index, the vegetation-soil synergy index (NVSSI), is constructed. By utilizing the opposite spectral response characteristics of the two indices to plantations and natural forests, the numerical differences between the two types of forest land are initially amplified.

[0025] S4: Perform hyperbolic tangent Tanh nonlinear transformation on each pixel of the NVSSI raster data to generate the final hyperbolic tangent vegetation-soil synergy index (TVSSI). By leveraging the inherent mathematical characteristics of the S-curve of the Tanh function, subtle feature differences are amplified in the middle variation range of the index, and convergence abnormal noise is compressed in the high and low extreme value range, thus eliminating the saturation drawbacks of traditional linear exponents in principle.

[0026] S5: The CCDC continuous temporal disturbance detection algorithm is adopted, and six sets of band temporal reflectance data are selected: blue, green, red, near-infrared, shortwave infrared 1, and shortwave infrared 2. Robust iterative weighted least squares method is used for pixel-by-pixel temporal fitting to automatically identify annual spectral abrupt change breakpoints and count the per-pixel per-year disturbance frequency (FD). Frequently disturbed plots with FD≥3 are designated as artificial forest training samples, and plots with FD=0 and continuously stable as forest land in the China Land Cover Dataset (CLCD) are designated as natural forest training.

[0027] S6: Based on high-resolution historical satellite imagery from CCDC and Google Earth Pro, a training sample dataset of independent planted forests and natural forests was constructed. In addition, a validation dataset was constructed by combining field surveys. This dataset was not used in model training throughout the process and was only used for subsequent quantitative verification of accuracy.

[0028] S7: Deploy the SNIC simple non-iterative clustering algorithm on the GEE platform, set the seed spacing to 15, compactness coefficient to 4, and neighborhood search size to 128, and complete the superpixel segmentation of the entire study area. Use contiguous superpixel objects to replace scattered single pixels as the basic classification unit, extract the average spectral and spatial texture information of the objects, and suppress the salt-and-pepper noise in the classification.

[0029] S8: Integrates multi-source information to build a multi-dimensional fusion feature database. The feature components include Landsat original B1~B7 bands, NDVI_fitted, SAVI, TVSSI index, NDVI derived time series indicators (maximum value, minimum value, coefficient of variation, variance, amplitude, phase, fitting residual), and DEM derived topographic parameters (elevation, slope, aspect), realizing the collaborative input of multi-source information such as spectrum, phenological time series, and topography.

[0030] S9: Randomly split all training samples generated in S5 at a ratio of 7:3, with 70% of the samples used as the training subset and 30% as the internal test subset; the random forest classifier has a fixed number of decision trees of 100, and performs binary classification operations on planted forests and natural forests based on the multidimensional feature dataset, outputting classification raster results.

[0031] S10: Quantitative evaluation is carried out in two dimensions: The first dimension selects four quantitative indicators, namely Cohen's d, FisherRatio, T Statistic, and Overlap Coefficient, and compares the forest area classification performance of four indices, namely NDVI_fitted, SAVI, NVSSI, and TVSS. The second dimension relies on the validation samples to construct a confusion matrix, solve the overall accuracy of OA and the Kappa coefficient to evaluate the overall classification level of the model.

[0032] S11: Output a 30m resolution spatial raster of plantations in Dingxi City for 2025 based on the optimized model, and output a standardized spatial distribution map of plantations.

[0033] S12: Multiple verification sites were randomly selected from different afforestation plots, and spatial and land cover comparisons were conducted using real-world field photographs and concurrent high-resolution satellite imagery. Furthermore, the reliability of the results was verified by comparing them with existing datasets GNPF and SPDT.

[0034] Example 1: Remote Sensing Preprocessing and TVSSI Index Construction Perform steps S1-S4, using the GEE platform to screen valid Landsat 9 remote sensing images from 2025, limiting cloud cover to ≤20% for each scene. Then, use JavaScript to batch perform atmospheric correction, cloud and cloud shadow masking, removal of invalid snow and ice pixels, annual median composite, and regional vector clipping to obtain standardized reflectance data and remote sensing feature data for the study area (see...). Figure 2-19 The process involves calculating NDVI based on preprocessed bands, fitting the data using harmonic formulas to obtain NDVI_fitted, calculating SAVI with a soil coefficient L=0.5, performing a ratio calculation of NVSSI=NDVI_fitted / SAVI, and then obtaining TVSSI through pixel-by-pixel nonlinear transformation using the Tanh function. Finally, a TVSSI global spatial raster is output and visualized. The results are shown below. Figure 24 As shown in the figure. This embodiment overcomes the defects of conventional vegetation index saturation and bare soil interference by using the Tanh nonlinearity through a two-level exponential construction.

[0035] Example 2: CCDC Time Series Analysis and Automatic Extraction of Classification Samples Following step S5, Landsat time-series data in six bands (blue, green, red, near-infrared, SWIR1, and SWIR2) were selected. A robust iterative weighted least squares algorithm was used to perform full-time CCDC fitting, identifying spectral abrupt change breakpoints year by year and statistically analyzing the perturbation frequency (FD) of each pixel. Plots with FD ≥ 3 were designated as plantation training samples, and areas with FD = 0 and long-term stable CLCD forests were designated as natural forest samples. The spatial distribution results of the perturbation frequency across the entire region (…) Figure 21 ), and randomly selected three representative forest locations to plot the annual NDVI observation and fitting time series curves ( Figure 22 Finally, all sample points are summarized to generate a sample spatial distribution map. Figure 23 ).

[0036] Example 3: SNIC Superpixel Segmentation and Multidimensional Feature Construction Steps S7 and S8 were executed, and the SNIC superpixel segmentation algorithm was run on the GEE platform with fixed parameters: seed spacing of 15, compactness of 4, and neighborhood size of 128. Object-oriented block processing was applied to the entire study area, and the average spectral and texture features within each superpixel unit were extracted. The segmentation results are shown below. Figure 20 The original bands, NDVI_fitted, SAVI, TVSSI, NDVI time-series derived parameters, elevation, slope, and aspect topographic factors were summarized to construct a complete multidimensional feature dataset, providing multi-source feature support for subsequent random forest classification.

[0037] Example 4: Quantitative Verification of the Discriminating Performance of Four Types of Indices In the S10 index evaluation step, random samples were taken from the extracted plantation and natural forest samples. The sample value distribution of the four indices NDVI_fitted, SAVI, NVSSI, and TVSSI was statistically analyzed, and a scatter plot of the distribution was plotted for testing. Figure 25 Furthermore, based on the evaluation indicators, four quantitative indicators—Cohen's d, Fisher Ratio, T Statistic, and Overlap Coefficient—a bar chart was generated. Figure 26 This confirms the distinguishing advantage of TVSSI compared to traditional indices.

[0038] Example 5: Random Forest Classification and Model Accuracy Statistics Following the S9 and S10 classification steps, the samples were split in a 7:3 ratio, and 100 decision trees were set up in a random forest. Automatic classification of planted and natural forests across the entire region was completed based on a multidimensional feature set. A confusion matrix was constructed using independent validation samples, and five classification accuracy indicators (OA, Kappa, PA, UA, and F1) were calculated. The results were then summarized and a precision statistics chart was plotted. Figure 28 The reliability of the method of the present invention is verified by quantitative analysis of the overall classification level.

[0039] Example 6: Field and High-Resolution Image Comparison Verification The S12 verification step was performed, and verification points were randomly set up in various landforms in the southern gully area, northern hilly area, and plain area of ​​the region. Real-world forest photos were taken on-site, and corresponding high-resolution satellite imagery from Google Earth Pro was retrieved. The plantation distribution map output by this invention was compared one by one with the real-world locations, high-resolution image patch locations, and existing research datasets (SPDT, GNPF). The comparison results are as follows: Figure 27 As shown, the verification of the results has been completed.

Claims

1. A remote sensing mapping method for plantations integrating the TVSSI index, characterized in that, The method includes the following steps in sequence: S1: Remote sensing image acquisition and preprocessing to obtain a well-organized and usable remote sensing dataset of the study area; S2: Based on the preprocessed Landsat multiband surface reflectance data, the normalized vegetation index (NDVI) is calculated. The annual continuous NDVI time series is fitted and denoised using a harmonic analysis model. The soil-regulated vegetation index (SAVI) is calculated simultaneously. S3: Constructing the vegetation-soil synergistic index NVSSI based on NDVI_fitted by HA fitting and SAVI using NDVI harmonic analysis; NVSSI = NDVI_fitted / SAVI; S4: Perform hyperbolic tangent transform on the NVSSI to construct the nonlinear hyperbolic tangent vegetation-soil synergy index TVSSI; TVSSI = Tanh (NVSSI) = (e^NVSSI − e^(−NVSSI)) / (e^NVSSI + e^(−NVSSI)) S5: Automatically generate training samples for planted forests and natural forests based on the CCDC temporal disturbance detection algorithm for continuous change detection and classification; S6: Build an independent validation sample dataset based on high-resolution imagery and field measurement data; S7: Use the SNIC algorithm to perform superpixel segmentation and complete object-oriented feature extraction; S8: Integrates multiple indicators such as spectral, temporal, topographic, and TVSSI to construct a multi-dimensional fusion feature set; S9: After the samples are stratified, the random forest algorithm is used to classify and identify planted forests and natural forests. S10: Conduct multi-parameter quantitative evaluation from two dimensions: classification accuracy and index discrimination performance; conduct a horizontal comparison of the advantages and disadvantages of four indices, NDVI_fitted, SAVI, NVSSI, and TVSSI, in distinguishing between plantations and natural forests, and quantitatively verify the anti-interference and discrimination performance of the TVSSI index.

2. The method for remote sensing mapping of plantations by integrating the TVSSI index according to claim 1, characterized in that, Specifically, S1 involves: using the GEE geographic big data cloud computing platform to select and call Landsat 9 Collection 2 Level-2 surface reflectance images, with the selection criterion being that the cloud cover of a single image is controlled to be ≤20%; using the JavaScript programming language to batch process the original remote sensing images in sequence, performing atmospheric correction processing, cloud masking processing, cloud shadow masking processing, removal of invalid pixels due to ice and snow, annual time series median synthesis, and administrative boundary vector clipping, retaining only the effective remote sensing data of the study area after clipping.

3. The method for remote sensing mapping of plantations incorporating the TVSSI index according to claim 2, characterized in that, S2 specifically involves: calculating NDVI based on preprocessed Landsat multiband reflectance data; in addition, using a harmonic analysis model to fit and denoise the annual continuous NDVI time series to obtain NDVI_fitted time series data that eliminates cloud and fog and random noise interference; and simultaneously calculating the soil-regulated vegetation index (SAVI). The formulas for calculating NDVI, NDVI_fitted, and SAVI are as follows: NDVI = (NIR - RED) / (NIR + RED); SAVI= (NIR−RED)×(1+L) / (NIR+RED+L); NDVI_fitted (t)=a+b×t+∑(k=1~n)[α k cos (2πkt / T)+β k sin (2πkt / T)] In the formula: NIR is the near-infrared band, RED is the red band, t represents the number of days in the time series, T is taken as the annual phenological cycle of vegetation (365 days), a and b are the overall trend fitting coefficients, and α k β k is the harmonic oscillation coefficient, n is the harmonic fitting order, and the soil adjustment coefficient L is fixed at 0.

5.

4. The method for remote sensing mapping of plantations fused with TVSSI index according to claim 3, characterized in that: Specifically, S3 involves: using the ratio of the harmonic-fitted NDVI_fitted sequence to the soil-corrected SAVI data to construct the intermediate synergistic index NVSSI, calculated as: NVSSI = NDVI_fitted / SAVI. This amplifies the difference in the original spectra of plantations and natural forests by leveraging the inverse response law between the two types of indices.

5. The method for remote sensing mapping of plantations fused with TVSSI index according to claim 4, characterized in that: S4 specifically involves: performing a Tanh hyperbolic tangent nonlinear transform on the obtained NVSSI raster data to construct the final TVSSI index. Leveraging the S-curve characteristics of the Tanh function, the characteristic difference is sensitively amplified within the intermediate range of the index, and convergence extreme noise is compressed at the maximum and minimum points of the index to overcome the spectral saturation defect of traditional linear vegetation indices under high vegetation cover conditions. The calculation formula is as follows: TVSSI = Tanh (NVSSI) = (e^NVSSI − e^(−NVSSI)) / (e^NVSSI + e^(−NVSSI)).

6. The method for remote sensing mapping of plantations fused with TVSSI index according to claim 5, characterized in that, Specifically, S5 involves: employing the Continuous Change Detection and Classification (CCDC) algorithm, selecting time-series reflectance data from six Landsat bands (blue, green, red, near-infrared, shortwave infrared 1, and shortwave infrared 2), and using a robust iterative weighted least squares method to perform time-series regression fitting pixel by pixel, automatically identifying breakpoints of temporal disturbances in each pixel; statistically analyzing the multi-year disturbance frequency (FD) of a single pixel, defining pixel regions with FD≥3 as training samples for planted forests; and then defining regions with FD=0 and maintaining unchanged local forest attributes as training samples for natural forests, thereby achieving automated extraction of samples from the entire region.

7. The method for remote sensing mapping of plantations incorporating the TVSSI index according to claim 6, characterized in that, Specifically, S6 involves: using high-resolution historical images from Google Earth Pro, photos of locations collected through field surveys, and field investigation records, manually setting up evenly distributed sampling points, and separately constructing two independent verification sample sets: artificial forests and natural forests. These samples do not participate in model training but are only used for subsequent classification accuracy calculations and TVSSI index discrimination capability verification.

8. The method for remote sensing mapping of plantations fused with the TVSSI index according to claim 7, characterized in that, Specifically, S7 involves running the Simple Non-Iterative Clustering (SNIC) superpixel segmentation algorithm in the GEE cloud computing environment. The algorithm parameters are set as follows: seed spacing 15, compactness coefficient 4, and neighborhood search size 128. The generated superpixel objects are used as the basic classification units. The average spectral features and spatial texture features inside the objects are extracted in a unified manner to eliminate salt-and-pepper noise caused by single-pixel classification and improve the spatial integrity of the classified patches.

9. A remote sensing mapping method for plantations based on TVSSI index according to claim 8, wherein S8 specifically comprises: integrating multi-source information to construct a full-dimensional fusion feature set, the features specifically including: Landsat original B1~B7 single-band reflectance features, NDVI_fitted, SAVI, TVSSI index features, NDVI time-series derived features, and DEM topographic features; wherein the NDVI time-series derived features include seven types of indicators: annual maximum value, annual minimum value, coefficient of variation, variance, amplitude, phase, and fitting residual; and the DEM topographic features include three topographic parameters: elevation, slope, and aspect.

10. The remote sensing mapping method for planted forests fused with TVSSI index according to claim 9, wherein S9 specifically comprises: randomly splitting all training samples obtained in S5 according to a 7:3 ratio, using 70% of the samples as a training subset for fitting the random forest model, and using the remaining 30% of the samples as an internal test subset; setting the random forest classifier parameters to have 100 decision trees, performing binary classification operations on planted forests and natural forests based on the constructed multidimensional feature set, and outputting classification raster results; S10 specifically involves: on the one hand, constructing a confusion matrix based on validation samples, calculating the overall accuracy OA and Kappa coefficient to evaluate the overall classification effect; on the other hand, selecting evaluation indicators, including four quantitative indicators: inter-class difference, comprehensive discrimination ability, difference significance, and overlap coefficient, and comparing the advantages and disadvantages of the four indices NDVI_fitted, SAVI, NVSSI, and TVSSI in distinguishing between plantations and natural forests, and quantitatively verifying the anti-interference and discrimination performance of the TVSSI index.

Images