A remote sensing inversion method for soil organic matter layer thickness based on a hybrid model

Abstract

This invention discloses a remote sensing inversion method for soil organic matter layer thickness based on a hybrid model, including satellite-ground heterogeneous data fusion, extreme value-aware sample partitioning, quantile transform preprocessing, wavelet-PCA feature engineering, six-branch hybrid model construction and training, and remote sensing inversion steps. By constructing a satellite-ground heterogeneous fused dataset, after feature selection and preprocessing, it is input into a six-branch parallel CNN-BiLSTM-Transformer hybrid model. A multi-granularity spatiotemporal feature interactive attention mechanism is used to explicitly model the coupling relationship between "climate-vegetation" and "spectrum-topography." Five environmental factors are dynamically weighted through a feature group adaptive gating fusion mechanism. Finally, an uncertainty-aware dual-head prediction network synchronously outputs predicted values ​​and confidence intervals. This method addresses the problems of insufficient long-term dependency capture, a single multi-source feature fusion mechanism, smooth extreme value sample prediction, and lack of uncertainty quantification in existing methods. It achieves high inversion accuracy and strong robustness, providing data support for ecological protection and carbon cycle assessment.

Description

Technical Field

[0001] This invention relates to the fields of remote sensing inversion and soil detection technology, and to a remote sensing inversion method for soil organic matter layer thickness based on a hybrid model. Specifically, it relates to a remote sensing inversion method for soil organic matter layer thickness in typical areas of the Qinghai-Tibet Plateau based on a CNN-BiLSTM-Transformer hybrid model. This method is suitable for large-scale, high-precision soil organic matter layer thickness inversion in alpine meadow areas, providing data support for ecological protection, carbon cycle assessment, and climate change monitoring. Background Technology

[0002] As the "Roof of the World" and the "Water Tower of Asia," the Qinghai-Tibet Plateau's ecological environment plays a crucial indicative role in global climate change. The alpine meadow ecosystem of the Qinghai-Tibet Plateau possesses enormous carbon storage, primarily in the form of organic matter. Soil carbon storage, particularly the organic matter layer thickness (SOMT), directly reflects soil carbon sequestration capacity and ecosystem health. SOMT is a core carrier of vital ecological functions such as soil productivity, water conservation, and carbon fixation in alpine meadows, and is of great significance for maintaining plateau ecological stability and the well-being of herders. Therefore, accurately obtaining the distribution of SOMT over a large area is of significant practical importance for assessing regional carbon cycling, ecological protection, and climate change monitoring.

[0003] Traditional methods for obtaining organic matter layer thickness rely on field sampling and laboratory analysis. While these methods offer high accuracy, they are costly in terms of manpower and resources, resulting in poor economic efficiency and timeliness, and are insufficient to meet the needs of large-scale dynamic monitoring. In recent years, the combination of remote sensing technology and deep learning has provided a new approach for soil property inversion. However, existing methods still have significant drawbacks: they often employ single models, making it difficult to simultaneously capture the complex spatial heterogeneity of the alpine meadow surface and the long-term vegetation evolution patterns, thus limiting inversion accuracy; obtaining field samples on the Qinghai-Tibet Plateau is costly and difficult, and under small sample conditions, complex deep learning models are prone to overfitting, while simple machine learning models cannot fit the nonlinear relationship between high-dimensional remote sensing data and organic matter layer thickness; the fusion strategies for multi-source heterogeneous data (climate, vegetation, topography, soil properties, etc.) lack specificity, resulting in insufficient information utilization and susceptibility to noise interference.

[0004] In the existing technologies, the relevant inversion methods mainly include traditional machine learning, single convolutional neural network (CNN) spatial feature extraction, and hybrid models of Transformer and CNN, but none of them can solve the core problems of insufficient multi-scale spatiotemporal feature capture, poor generalization ability for small samples, and insufficient fusion of multi-source data, making it difficult to meet the high-precision inversion requirements of organic matter layer thickness in the alpine meadow region of the Qinghai-Tibet Plateau. Summary of the Invention

[0005] The purpose of this invention is to provide a remote sensing inversion method for soil organic matter layer thickness based on a hybrid model. It aims to provide a method that integrates multi-source heterogeneous data with a six-branch hybrid deep learning model, captures spatiotemporal features through a multi-granularity spatiotemporal feature interaction attention mechanism (MGSTFIA) and a feature group adaptive gating fusion mechanism (FGAGF), and employs an uncertainty-aware dual-head prediction network (UADPN) to achieve high-precision, large-scale remote sensing inversion of soil organic matter layer thickness in alpine meadow areas of typical Qinghai-Tibet Plateau regions.

[0006] According to the purpose of this invention, a remote sensing inversion method for soil organic matter layer thickness based on a hybrid model is provided, comprising the following steps: S1. Data Acquisition: A heterogeneous fusion dataset of ground and satellite data was constructed, including ground-based measured data and satellite remote sensing data. The ground-based measured data consisted of soil organic matter layer thickness and soil type data in the alpine meadow region of the Qinghai-Tibet Plateau. The satellite data included Landsat time-series spectral data, KNDVI vegetation index time-series data, climate time-series data, solar radiation data, and static environmental data. As shown in Figures 21-27, this invention constructed a KNDVI spatial feature matrix with 38 time dimensions from 1986 to 2023. This time-series GIS mapping not only reflects the growth status of alpine meadow vegetation over the past 40 years, but also provides physically meaningful dynamic evolution inputs for the time-series feature branches in the hybrid model through the coupling of spatial and temporal dimensions, effectively solving the problem of lack of spatiotemporal coherence of single-point data in traditional methods.

[0007] S2. Sample partitioning: The samples are stratified based on the extreme value-aware binning strategy, and the training set and test set are divided by stratified K-fold cross-validation. S3. Data preprocessing transformation: Quantile transformation is performed on the target variable, organic matter layer thickness (SOMT), robust standardization is performed on the features, and soil type is encoded and embedded. S4. Feature Engineering: Wavelet transform and principal component analysis were performed on Landsat time-series spectral data and KNDVI vegetation index time-series data for dimensionality reduction. A random forest algorithm was used to select key auxiliary features and group them according to physical attributes. S5. Hybrid Model Construction and Training: A six-branch parallel CNN-BiLSTM-Transformer hybrid model is constructed. Multi-source features are fused through the multi-granularity spatiotemporal feature interaction attention mechanism (MGSTFIA) and the feature group adaptive gating fusion mechanism (FGAGF). The uncertainty-aware dual-head prediction network (UADPN) is used to output the predicted value and confidence interval, and the model is trained. S6. Remote sensing inversion: Input the preprocessed multi-source data into the trained model and output the soil organic matter layer thickness inversion results.

[0008] Further, in step S1, the Landsat time-series spectral data is multi-band surface reflectance data for 39 years from 1985 to 2023; the KNDVI vegetation index time-series data is the annual composite product of nuclear normalized vegetation index for 38 years from 1986 to 2023; the climate time-series data includes a 12-month series composed of multi-year monthly averages of potential evapotranspiration from 1995 to 2024 and temperature and precipitation from 2001 to 2024; the solar radiation data is the annual total solar radiation data for 15 years from 2010 to 2024; the ground-based measured data is obtained by stratified random sampling, excavating soil profiles down to the mineral soil layer, measuring the vertical thickness of the organic matter layer, and recording environmental information.

[0009] Further, in step S2, the extreme value-aware binning strategy includes: calculating the extreme value threshold based on the interquartile range method, where the upper extreme value threshold is... The lower extreme threshold is The samples are divided into high extreme value samples, low extreme value samples, and normal samples. Normal samples are binned equally according to percentiles, and divided into 10 bins. Extreme value samples are evenly distributed to each bin through a round-robin allocation strategy. Based on bin coding, hierarchical 10-fold cross-validation is used to ensure that the distribution of the target variable in each fold is consistent with the overall dataset.

[0010] Further, in step S3, the quantile transformation uses QuantileTransformer to map the non-normally distributed target value of organic matter layer thickness (SOMT) to a standard normal distribution; the robust standardization process uses RobustScaler to standardize the eight feature groups based on the median and interquartile range; the encoding embedding is to encode the soil type into integers and then map it into a 16-dimensional continuous vector through the embedding layer; the data preprocessing transformation also includes calculating adaptive weights for samples, assigning dynamic weights of 1.5 to 5.0 times to extreme value samples.

[0011] Further, in step S4, the wavelet transform uses Morlet continuous wavelet transform to extract multi-scale time-frequency features; the principal component analysis dimensionality reduction compresses the Landsat wavelet coefficients to 15 dimensions and the KNDVI wavelet coefficients to 10 dimensions; the random forest feature selection constructs 200 decision trees with a maximum depth of 8, calculates the importance score of each feature based on variance reduction, selects the top 50 key features, and classifies them into 8 feature groups according to physical attributes: temperature, precipitation, evapotranspiration, solar radiation, topography, aridity, soil properties, and others.

[0012] Further, the six-branch parallel hybrid model described in step S5 includes: a climate time series branch, which uses a three-layer one-dimensional convolutional neural network to process monthly time series data of temperature, precipitation, and evapotranspiration, with channel changes from 3->32->64->64; a solar radiation branch, which uses a two-layer one-dimensional convolutional neural network to process 15-year solar radiation time series data, with channel changes from 1->32->64; a Landsat spectral branch, which uses a two-layer fully connected network (MLP) to map the 15-dimensional spectral features after wavelet-PCA dimensionality reduction to a 64-dimensional hidden space; a KNDVI vegetation branch, which uses a two-layer fully connected network (MLP) to map the 10-dimensional vegetation features after wavelet-PCA dimensionality reduction to a 64-dimensional hidden space; a static environment branch, which uses a two-layer fully connected network (MLP) to process topography, aridity, and soil property features to a 64-dimensional hidden space; and a soil type branch, which uses an embedding layer to map soil type integer encodings to 16-dimensional continuous vectors.

[0013] Furthermore, the BiLSTM module is a three-layer bidirectional long short-term memory network with a hidden layer dimension of 64. After bidirectional splicing, it becomes 128-dimensional, which is reduced to 64-dimensional by linear projection, and the dropout rate is 0.35. The Transformer encoder module contains a 3-layer encoder, with 4 attention heads in each layer. The feedforward network has a dimension of 256 and the activation function is GELU.

[0014] Further, in step S5, the Feature Group Adaptive Gated Fusion Mechanism (FGAGF) dynamically calculates the weights of five feature groups—climate, solar radiation, spectrum, vegetation, and static attributes—through a Softmax gated network. Each feature group is then transformed by an independent projection network and fused according to its weights. The Uncertainty Aware Dual-Head Prediction Network (UADPN) includes a mean prediction head and a log-variance prediction head, which simultaneously output the predicted mean and the predicted variance. The log-variance is pruned and limited to the range of -10 to 10.

[0015] Further, in step S5, the model training employs an extreme value-aware hybrid loss function, fusing Huber loss, negative log-likelihood loss, and quantile loss, with weights of 0.4, 0.3, and 0.3 respectively, and quantile loss quantiles of 0.1, 0.5, and 0.9; the AdamW optimizer is used with an initial learning rate of 0.0015 and a weight decay of 0.08; a course learning strategy is employed, gradually increasing the influence of extreme value samples in the first 200 training rounds; cosine annealing learning rate scheduling is used, with 40 warm-up rounds and a maximum training round count of 600; gradient norm pruning is employed with a threshold of 1.0; and an early stopping mechanism is used, stopping training if there is no improvement in R² on the validation set for 100 consecutive rounds.

[0016] Further, in step S6, the enhancement during testing employs a 7-cycle prediction. Each time, Gaussian noise with increasing magnitude is added to the climate, solar radiation, Landsat spectral features processed by wavelet transform and PCA dimensionality reduction, and KNDVI vegetation features. The noise scale is 0.003×i (i=0,1,...,6), and the mean of the 7 predictions is taken as the final normal space prediction value. The inverse transformation restores the normal space prediction value to the original scale through the inverse function of quantile transformation, and the value range of the prediction value is clipped. The inversion result is cross-validated with 10-fold cross-validation, and the model determination coefficient R² reaches 0.60, and the root mean square error RMSE reaches 4.91cm.

[0017] The beneficial effects of this invention are: This invention, through the fusion of heterogeneous satellite and ground data and an innovative six-branch hybrid model architecture, possesses significant advantages: a six-branch parallel feature extraction network extracts multiple types of temporal and static features; a three-layer BiLSTM captures bidirectional temporal dependencies; a three-layer Transformer encoder models global long-distance relationships; a multi-granularity spatiotemporal feature interaction attention mechanism (MGSTFIA) explicitly constructs the coupling relationship between "climate-vegetation" and "spectrum-topography"; a feature group adaptive gating fusion mechanism (FGAGF) dynamically weights five types of environmental factors to achieve "site-specific" feature fusion, comprehensively mining spatiotemporal information; Morlet wavelet transform and PCA dimensionality reduction effectively extract multi-scale time-frequency features from the temporal data; random forest feature selection and robust normalization preprocessing strategies effectively reduce noise interference; an extreme value-aware hybrid loss function combined with quantile transformation overcomes the "smoothing effect" of traditional models, significantly improving the accuracy of extreme value inversion; and an uncertainty-aware dual-head prediction network (UADPN) simultaneously outputs predicted values ​​and confidence intervals, achieving "prediction + risk assessment." After ten-fold cross-validation, the model R... 2 The accuracy reached 0.60, with an RMSE of 4.91 cm, demonstrating superior inversion accuracy compared to single models. This method is well-suited to the complex surface features of alpine meadows, exhibits strong robustness, and can achieve large-scale, high-precision organic matter layer thickness inversion, providing reliable data support for ecological protection and carbon cycle assessment. Attached Figure Description

[0018] Figure 1 This is a flowchart of the remote sensing inversion technology of the present invention; Figure 2 This is a sample point distribution diagram in an embodiment of the present invention; Figure 3 This is a schematic diagram of the extreme value sensing bin structure in an example of the present invention; Figure 4 This is a schematic diagram of the data preprocessing transformation structure in an example of the present invention; Figure 5 This is a schematic diagram of the Landsat spectral branching structure in an example of the present invention; Figure 6 This is a schematic diagram of the climate time series branch structure in an example of the present invention; Figure 7 This is a schematic diagram of the solar radiation branch structure in an example of the present invention; Figure 8 This is a schematic diagram of the static environment branch structure in an example of the present invention; Figure 9 This is a schematic diagram of the KNDVI vegetation branching structure in an example of the present invention; Figure 10 This is a schematic diagram of the soil type branch structure in an example of the present invention; Figure 11 This is a schematic diagram of the Multi-Granularity Spatiotemporal Feature Interactive Attention Mechanism (MGSTFIA) structure in an example of the present invention; Figure 12 This is a schematic diagram of the Bidirectional Long Short-Term Memory (BiLSTM) network structure in an example of the present invention; Figure 13 This is a schematic diagram of the Transformer encoder structure in an example of the present invention; Figure 14 This is a schematic diagram of the Feature Group Adaptive Gated Fusion (FGAGF) structure in an example of the present invention; Figure 15 This is a schematic diagram of the residual deep semantic fusion (RDSC) module structure in an example of the present invention; Figure 16 This is a schematic diagram of the uncertainty-aware dual-head prediction network (UADPN) module structure in an example of the present invention; Figure 17 This is a scatter plot of the 10-fold cross-validation of the CNN-BiLSTM-Transformer model in this invention example; Figure 18 This is a scatter plot of the 10-fold cross-validation of the CNN model in this invention example; Figure 19 This is a scatter plot of the 10-fold cross-validation of the BiLSTM model in this invention example; Figure 20 This is a scatter plot of the ten-fold cross-validation of the Transformer model in this invention example; Figure 21 This is a schematic diagram of the annual KNDVI time-series GIS mapping from 1986 to 1991 in an embodiment of the present invention; Figure 22 This is a schematic diagram of the annual KNDVI time-series GIS mapping from 1992 to 1997 in an embodiment of the present invention; Figure 23 This is a schematic diagram of the annual KNDVI time-series GIS mapping from 1998 to 2003 in an embodiment of the present invention; Figure 24This is a schematic diagram of the KNDVI time-series GIS mapping from 2004 to 2009 in an embodiment of the present invention. Figure 25 This is a schematic diagram of the KNDVI time-series GIS mapping from 2010 to 2015 in an embodiment of the present invention. Figure 26 This is a schematic diagram of the KNDVI time-series GIS mapping from 2016 to 2021 in an embodiment of the present invention. Figure 27 This is a schematic diagram of the KNDVI time-series GIS mapping from 2022 to 2023 in an embodiment of the present invention. Detailed Implementation

[0019] The specific embodiments of the present invention will be further described below. It should be noted that these descriptions are for the purpose of aiding understanding the present invention, but do not constitute a limitation thereof. Furthermore, the technical features involved in the various embodiments of the present invention described below can be combined with each other as long as they do not conflict with each other.

[0020] Example 1 like Figures 1-27 As shown, a remote sensing inversion method for soil organic matter layer thickness based on a hybrid model includes the following steps: S1. Data Acquisition: A heterogeneous fusion dataset of ground and satellite data was constructed, including ground-based measured data and satellite remote sensing data. The ground-based measured data consisted of soil organic matter layer thickness and soil type data in the alpine meadow region of the Qinghai-Tibet Plateau. The satellite data included Landsat time-series spectral data, KNDVI vegetation index time-series data, climate time-series data, solar radiation data, and static environmental data. As shown in Figures 21-27, this invention constructed a KNDVI spatial feature matrix with 38 time dimensions from 1986 to 2023. This time-series GIS mapping not only reflects the growth status of alpine meadow vegetation over the past 40 years, but also provides physically meaningful dynamic evolution inputs for the time-series feature branches in the hybrid model through the coupling of spatial and temporal dimensions, effectively solving the problem of lack of spatiotemporal coherence of single-point data in traditional methods.

[0021] S2. Sample partitioning: The samples are stratified based on the extreme value-aware binning strategy, and the training set and test set are divided by stratified K-fold cross-validation. S3. Data preprocessing transformation: Perform quantile transformation on the target variable, robust standardization on the features, and encoding and embedding of soil types; S4. Perform wavelet transform and principal component analysis to reduce the dimensionality of Landsat time-series spectral data and KNDVI vegetation index time-series data, and use the random forest algorithm to screen key auxiliary features and group them according to physical attributes. S5. Hybrid Model Construction and Training: A six-branch parallel CNN-BiLSTM-Transformer hybrid model is constructed. Multi-source features are fused through the multi-granularity spatiotemporal feature interaction attention mechanism (MGSTFIA) and the feature group adaptive gating fusion mechanism (FGAGF). The uncertainty-aware dual-head prediction network (UADPN) is used to output the predicted value and confidence interval, and the model is trained. S6. Remote sensing inversion: Input the preprocessed multi-source data into the trained model and output the soil organic matter layer thickness inversion results.

[0022] Further, in step S1, the Landsat time-series spectral data is multi-band surface reflectance data for 39 years from 1985 to 2023; the KNDVI vegetation index time-series data is the annual composite product of nuclear normalized vegetation index for 38 years from 1986 to 2023; the climate time-series data includes a 12-month series composed of multi-year monthly averages of potential evapotranspiration from 1995 to 2024 and temperature and precipitation from 2001 to 2024; the solar radiation data is the annual total solar radiation data for 15 years from 2010 to 2024; the ground-based measured data is obtained by stratified random sampling, excavating soil profiles down to the mineral soil layer, measuring the vertical thickness of the organic matter layer, and recording environmental information.

[0023] Further, in step S2, the extreme value-aware binning strategy includes: calculating the extreme value threshold based on the interquartile range method, where the upper extreme value threshold is... The lower extreme threshold is The samples are divided into high extreme value samples, low extreme value samples, and normal samples. Normal samples are binned equally according to percentiles, and divided into 10 bins. Extreme value samples are evenly distributed to each bin through a round-robin allocation strategy. Based on bin coding, hierarchical 10-fold cross-validation is used to ensure that the distribution of the target variable in each fold is consistent with the overall dataset.

[0024] Further, in step S3, the quantile transformation uses QuantileTransformer to map the non-normally distributed target value (organic matter layer thickness) to a standard normal distribution; the robust standardization process uses RobustScaler to standardize the eight feature groups based on the median and interquartile range; the encoding embedding is to encode the soil type into integers and then map it into a 16-dimensional continuous vector through the Embedding layer; the data preprocessing transformation also includes calculating adaptive weights for samples and assigning dynamic weights of 1.5 to 5.0 times to extreme value samples.

[0025] Further, in step S4, the wavelet transform uses Morlet continuous wavelet transform to extract multi-scale time-frequency features; the principal component analysis dimensionality reduction compresses the Landsat wavelet coefficients to 15 dimensions and the KNDVI wavelet coefficients to 10 dimensions; the random forest feature selection constructs 200 decision trees with a maximum depth of 8, calculates the importance score of each feature based on variance reduction, selects the top 50 key features, and classifies them into 8 feature groups according to physical attributes: temperature, precipitation, evapotranspiration, solar radiation, topography, aridity, soil properties, and others.

[0026] Further, the six-branch parallel hybrid model described in step S5 includes: a climate time series branch, which uses a three-layer one-dimensional convolutional neural network to process monthly time series data of temperature, precipitation, and evapotranspiration, with channel changes from 3->32->64->64; a solar radiation branch, which uses a two-layer one-dimensional convolutional neural network to process 15-year solar radiation time series data, with channel changes from 1->32->64; a Landsat spectral branch, which uses a two-layer fully connected network (MLP) to map the 15-dimensional spectral features after wavelet-PCA dimensionality reduction to a 64-dimensional hidden space; a KNDVI vegetation branch, which uses a two-layer fully connected network (MLP) to map the 10-dimensional vegetation features after wavelet-PCA dimensionality reduction to a 64-dimensional hidden space; a static environment branch, which uses a two-layer fully connected network (MLP) to process topography, aridity, and soil property features to a 64-dimensional hidden space; and a soil type branch, which uses an embedding layer to map soil type integer encodings to 16-dimensional continuous vectors.

[0027] Furthermore, the BiLSTM module is a three-layer bidirectional long short-term memory network with a hidden layer dimension of 64. After bidirectional splicing, it becomes 128-dimensional, which is reduced to 64-dimensional by linear projection, and the dropout rate is 0.35. The Transformer encoder module contains a 3-layer encoder, with 4 attention heads in each layer. The feedforward network has a dimension of 256 and the activation function is GELU.

[0028] Further, in step S5, the Feature Group Adaptive Gated Fusion Mechanism (FGAGF) dynamically calculates the weights of five feature groups—climate, solar radiation, spectrum, vegetation, and static attributes—through a Softmax gated network. Each feature group is then transformed by an independent projection network and fused according to its weights. The Uncertainty Aware Dual-Head Prediction Network (UADPN) includes a mean prediction head and a log-variance prediction head, which simultaneously output the predicted mean and the predicted variance. The log-variance is pruned and limited to the range of -10 to 10.

[0029] Further, in step S5, the model training employs an extreme value-aware hybrid loss function, fusing Huber loss, negative log-likelihood loss, and quantile loss, with weights of 0.4, 0.3, and 0.3 respectively, and quantile loss quantiles of 0.1, 0.5, and 0.9; the AdamW optimizer is used with an initial learning rate of 0.0015 and a weight decay of 0.08; a course learning strategy is employed, gradually increasing the influence of extreme value samples in the first 200 training rounds; cosine annealing learning rate scheduling is used, with 40 warm-up rounds and a maximum training round count of 600; gradient norm pruning is employed with a threshold of 1.0; and an early stopping mechanism is used, stopping training if there is no improvement in R² on the validation set for 100 consecutive rounds.

[0030] Further, in step S6, the enhancement during testing employs a 7-cycle prediction. Each time, Gaussian noise with increasing magnitude is added to the climate, solar radiation, Landsat spectral features processed by wavelet transform and PCA dimensionality reduction, and KNDVI vegetation features. The noise scale is 0.003×i (i=0,1,...,6), and the mean of the 7 predictions is taken as the final normal space prediction value. The inverse transformation restores the normal space prediction value to the original scale through the inverse function of quantile transformation, and the value range of the prediction value is clipped. The inversion result is cross-validated with 10-fold cross-validation, and the model determination coefficient R² reaches 0.60, and the root mean square error RMSE reaches 4.91cm.

[0031] Example 2 like Figures 1-4 As shown, a remote sensing inversion method for soil organic matter layer thickness based on a hybrid model includes the following steps in its data acquisition and preprocessing process: 1. Field measurement data collection: In the typical area of ​​the Qinghai-Tibet Plateau (geographical range: 96°50′-99°20′E, 33°50′-35°40′N), 150 sampling points were randomly distributed in layers according to vegetation type (alpine grassland, alpine meadow, marsh meadow), altitude gradient (3800-4500m), and climate zone (semi-arid, semi-humid). A 1m×0.5m×0.8m soil profile was excavated at each sampling point. The vertical thickness of the organic matter layer (from the surface vegetation residue layer to the top of the mineral soil layer) was measured using a tape measure with an accuracy of 0.1cm. Environmental information such as latitude, longitude, altitude, and vegetation cover of the sampling point was recorded.

[0032] 2. Collection of multi-source satellite remote sensing data: Climate time-series data: Multi-year monthly average data of potential evapotranspiration from 1995 to 2024, and temperature and precipitation from 2001 to 2024, were collected from the National Meteorological Science Data Center; Vegetation phenology data: Time series data of Normalized Difference Vegetation Index (KNDVI) from 1986 to 2023, derived from Landsat's annual synthetic product and synthesized from monthly averages over many years; Landsat time-series spectral data: composite images from the annual growing season (June-September) of Landsat 5 / 7 / 8 satellite from 1985 to 2023, from which six bands were extracted: blue, green, red, near-infrared, shortwave infrared 1, and shortwave infrared 2. Solar radiation data: Annual average solar radiation data from 2010 to 2024; Aridity data: Annual average aridity data from 1995 to 2024, sourced from the National Tibetan Plateau Scientific Data Center; Static environmental characteristics: Topographic data such as elevation, slope, and aspect are derived from ASTER GDEM V3; Soil property data (20 indicators including available nitrogen and soil bulk density) are sourced from the China Soil Database. Soil type data were obtained from field sampling in typical areas of the Qinghai-Tibet Plateau.

[0033] 3. Construction of a heterogeneous satellite-ground fusion dataset: Based on the latitude and longitude coordinates of ground sampling points, a spatial multi-value extraction method was used to extract corresponding pixel values ​​from the aforementioned multi-source satellite remote sensing dataset. Ground-measured soil organic matter layer thickness data was spatially matched with satellite remote sensing features to construct a satellite-ground heterogeneous fusion dataset. This dataset includes ground-measured soil organic matter layer thickness and soil type as target variables and category features, and Landsat time-series spectra, KNDVI vegetation index, climate time series, solar radiation, and static environmental features provided by satellite remote sensing data as input variables.

[0034] 4. Extreme value perception sample partitioning: The extreme value thresholds are calculated based on the interquartile range method, with the upper extreme value threshold being Q3 + 1.5 × IQR and the lower extreme value threshold being Q1 - 1.5 × IQR. The samples are divided into high extreme value samples, low extreme value samples, and normal samples. The normal samples are binned equally according to percentiles, resulting in 10 bins. The extreme value samples are evenly distributed into each bin using a round-robin allocation strategy. Hierarchical 10-fold cross-validation is used based on the bin coding. 5. Data preprocessing and transformation: Target variable transformation: The organic matter layer thickness data were processed using QuantileTransformer, and the absolute value of the skewness of the transformed data was less than 0.3; Sample weight calculation: Adaptive sample weights are calculated based on the degree of deviation of the sample from the extreme value threshold. The weight of normal samples is 1.0, and the weight of extreme value samples is 1.5 to 5.0. Feature standardization: Eight feature groups—temperature, precipitation, evapotranspiration, solar radiation, topography, aridity, soil properties, and others—were standardized using RobustScaler. The median and interquartile range of each feature were calculated, and then standardized according to the formula... Scaling; Soil type encoding: After integer encoding of soil type, it is mapped to a 16-dimensional vector through the Embedding layer.

[0035] Example 3 As shown in Figure 1, a remote sensing inversion method for soil organic matter layer thickness based on a hybrid model includes the following steps: Landsat time series data processing: Morlet continuous wavelet transform was performed on 39 years of Landsat multi-band time series data to extract multi-scale time-frequency features; RobustScaler robust standardization was performed on the wavelet coefficients; Principal component analysis was performed on the standardized wavelet coefficient matrix to reduce dimensionality, and the first 15 principal components were retained.

[0036] KNDVI time series data processing: The same processing flow as Landsat time series data is adopted, and Morlet wavelet transform, wavelet coefficient matrix construction, RobustScaler standardization, and PCA dimensionality reduction are performed in sequence; PCA dimensionality reduction retains the first 10 principal components.

[0037] Auxiliary feature processing: 200 decision trees were constructed using random forest feature selection, with a maximum depth of 8 and a minimum number of samples per leaf node of 4. The model was trained using SOMT as the target variable, and the importance score of each feature was calculated. The top 50 most important features were selected and categorized into 8 feature groups based on physical attributes: temperature, precipitation, evapotranspiration, solar radiation, topography, aridity, soil properties, and others.

[0038] Feature organization: The climate feature tensor stacks the three sets of features—temperature, precipitation, and evapotranspiration—into a three-channel time series tensor with a shape of (number of samples, 12, 3), where 12 represents the months and 3 represents the number of variables; Landsat wavelet-PCA features are organized into (number of samples, 15); KNDVI wavelet-PCA features are organized into (number of samples, 10); solar radiation features are organized into (number of samples, 15); static feature vectors concatenate static features such as topography, aridity, and soil properties; soil type embedding vectors are organized into (number of samples, 16).

[0039] Example 4 like Figure 1 , Figures 5-16 As shown, a remote sensing inversion method for soil organic matter layer thickness based on a hybrid model includes the following steps in its hybrid model construction and training process: 1. Model structure parameters: The six-branch parallel feature extraction network includes a climate time series branch that uses a three-layer one-dimensional convolutional neural network. The input shape (batch, 12, 3) is transposed to (batch, 3, 12), and the channel changes from 3->32->64->64. Each layer contains Conv1D, BatchNorm1D, GELU activation and Dropout. Finally, it outputs 64-dimensional features through adaptive average pooling. The solar radiation branch uses a two-layer one-dimensional convolutional neural network. The input shape (batch, 15) is expanded to (batch, 1, 15), the channels change from 1->32->64, and the output is 64-dimensional features. The Landsat spectral branch uses a two-layer fully connected network (MLP) with the structure Linear(15->64)->LayerNorm->GELU->Dropout->Linear(64->64)->LayerNorm->GELU; KNDVI vegetation branches adopt a two-layer fully connected network (MLP) with the structure Linear(10->64)->LayerNorm->GELU->Dropout->Linear(64->64)->LayerNorm->GELU; The static environment branch uses a two-layer fully connected network (MLP) and is mapped to a 64-dimensional hidden space; The soil type branch maps the soil type integer encoding to a 16-dimensional continuous vector through the Embedding layer; The BiLSTM module is a three-layer bidirectional long short-term memory network with a hidden layer dimension of 64. After bidirectional splicing, it becomes 128-dimensional, which is reduced to 64-dimensional by linear projection. The dropout rate is 0.35. The Transformer encoder module contains a 3-layer encoder, with 4 attention heads per layer. The feedforward network has a dimension of 256, the activation function is GELU, and the dropout rate is 0.35.

[0040] 2. Multi-granularity spatiotemporal feature interaction attention mechanism (MGSTFIA): Climate-vegetation interaction attention uses climate features as queries, concatenates climate and vegetation features as key-value pairs, and calculates interaction features through multi-head attention. The spectral-terrain interaction attention uses spectral features as queries, and concatenates spectral and static features as key-value pairs. Global temporal self-attention is used to perform self-attention calculations on four types of features: climate, vegetation, spectrum, and statics. The output is then processed by mean pooling. The outputs from the three channels are concatenated and then fused into an MLP (192→128→64 dimensions) to obtain 64-dimensional interactive features.

[0041] 3. Feature Group Adaptive Gated Fusion Mechanism (FGAGF): Five groups of 64-dimensional features, including climate, solar radiation, spectrum, vegetation, and static attributes, are concatenated into a 320-dimensional vector; normalized weights are generated by a gating network (320->64->5) and Softmax; each feature group is transformed by an independent projection network (64->64) and then fused according to the weights to obtain 64-dimensional gated features.

[0042] Deep Feature Fusion and Prediction: Residual Deep Semantic Fusion (RDSC) concatenates BiLSTM pooling features, interaction features, Transformer pooling features, gated features, and soil type embedding features (a total of 272 dimensions), and then fuses them through two fully connected network layers to obtain a 64-dimensional deep fusion representation; Uncertainty-Aware Dual-Head Prediction Network (UADPN) outputs the predicted mean from the mean prediction head and the predicted variance from the log-variance prediction head after the shared feature layer. The variance is cropped and limited to the range of -10 to 10.

[0043] Example 5 Figure 1 shows a remote sensing inversion method for soil organic matter layer thickness based on a hybrid model. The model inference process includes the following steps: 1. Test-Time Augmentation (TTA) Inference: Load optimal model weights; perform 7 iterations of prediction, with noise incrementing according to the following rule. = 0.003 × i (i=0,1,...,6); Gaussian noise with increasing magnitude is added to the climate time series features, solar radiation features, Landsat wavelet-PCA features, and KNDVI wavelet-PCA features; the mean of 7 predictions is taken as the final normal spatial prediction value.

[0044] 2. Inverse Transformation and Range Clipping: The predicted values ​​in the normal space are restored to their original scale using the inverse function of the quantile transformation; the range of the predicted values ​​is clipped to a reasonable range.

[0045] Example 6 As shown in Figure 1 and Figure 17- Figure 20 As shown, a remote sensing inversion method for soil organic matter layer thickness based on a hybrid model is described. Its model validation and inversion application include the following steps: 1. Model Validation: The results of the 10-fold cross-validation show that model R... 2 The RSI reached 0.60, and the RMSE was 4.91 cm, significantly outperforming a single CNN model (R²). 2 =0.53, RMSE=5.29cm), BiLSTM model (R 2 =0.39, RMSE=6.04cm) and Transformer model (R 2 =0.45, RMSE=5.71cm).

[0046] 2. Inversion Application: Preprocessed multi-source remote sensing data of a typical area of ​​the Qinghai-Tibet Plateau is input into the trained model to predict and obtain a 30m resolution map of soil organic matter layer thickness distribution across the entire region. This prediction result can be used to assess regional carbon storage, guide ecological protection and restoration, and monitor the impact of climate change on alpine meadow ecosystems.

[0047] As shown in Figure 17, the present invention uses the ten-fold cross-validation evaluation method to evaluate R. 2 The RMSE reached 4.91 (cm), with the horizontal axis representing the predicted value and the vertical axis representing the observed value. Different colors represent the verification results of different folds.

[0048] In summary, this invention has strong multi-scale feature capture capabilities: it extracts multiple types of temporal and static features through a six-branch parallel feature extraction network (including CNN, MLP and Embedding layers), captures bidirectional temporal dependencies through BiLSTM, and models global long-distance relationships through Transformer, thereby achieving comprehensive mining of multi-scale spatiotemporal features and adapting to the complex surface features of alpine meadows. This invention is based on explicit feature interaction of physical mechanisms: by using the multi-granularity spatiotemporal feature interaction attention mechanism (MGSTF), independent interaction channels of "climate-vegetation" and "spectrum-topography" are explicitly constructed, enabling the model to simulate the hydrothermal coupling effect and surface reflection mechanism in the soil development process, thereby improving the physical interpretability of features. This invention features a dynamic weight adaptive allocation of heterogeneous environmental factors: the Feature Group Adaptive Gated Fusion Mechanism (FGAGF) can dynamically calculate the confidence weights of five groups of features based on the habitat characteristics of the sample, achieving "locally adapted" feature fusion and effectively improving the model's generalization ability in complex spatial heterogeneity. This invention has strong extreme value sample prediction capability: the extreme value-aware hybrid loss function combined with the quantile transformation strategy overcomes the "smoothing effect" of traditional models and significantly improves the inversion accuracy for extreme high and low values; This invention is highly robust: the data preprocessing and transformation process and model structure design both take into account the effects of outliers and data skewness, adapt to the noise characteristics of remote sensing data, and the inversion results are stable and reliable.

[0049] Although embodiments of the present invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the present invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A remote sensing inversion method for soil organic matter layer thickness based on a hybrid model, characterized in that, Includes the following steps: S1. Data Acquisition: A heterogeneous fusion dataset of ground and satellite data is constructed, including ground-based measured data and satellite remote sensing data. The ground-based measured data includes soil organic matter layer thickness and soil type data in the alpine meadow region of the Qinghai-Tibet Plateau. The satellite data includes Landsat time-series spectral data, KNDVI vegetation index time-series data, climate time-series data, solar radiation data, and static environmental data. As shown in Figures 21-27, this invention constructs a KNDVI spatial feature matrix with 38 time dimensions from 1986 to 2023. This time-series GIS mapping not only reflects the growth status of alpine meadow vegetation over the past 40 years, but also provides physically meaningful dynamic evolution inputs for the time-series feature branches in the hybrid model through the coupling of spatial and temporal dimensions, effectively solving the problem of lack of spatiotemporal coherence of single-point data in traditional methods. S2. Sample partitioning: The samples are stratified based on the extreme value-aware binning strategy, and the training set and test set are divided by stratified K-fold cross-validation. S3. Data preprocessing transformation: Quantile transformation is performed on the target variable, organic matter layer thickness (SOMT), robust standardization is performed on the features, and soil type is encoded and embedded. S4. Feature Engineering: Wavelet transform and principal component analysis were performed on Landsat time-series spectral data and KNDVI vegetation index time-series data for dimensionality reduction. A random forest algorithm was used to select key auxiliary features and group them according to physical attributes. S5. Hybrid Model Construction and Training: A six-branch parallel CNN-BiLSTM-Transformer hybrid model is constructed. Multi-source features are fused through the multi-granularity spatiotemporal feature interaction attention mechanism (MGSTFIA) and the feature group adaptive gating fusion mechanism (FGAGF). The uncertainty-aware dual-head prediction network (UADPN) is used to output the predicted value and confidence interval, and the model is trained. S6. Remote sensing inversion: Input the preprocessed multi-source data into the trained model and output the soil organic matter layer thickness inversion results.

2. The remote sensing inversion method according to claim 1, characterized in that, In step S1, the Landsat time-series spectral data consists of multi-band surface reflectance data for 39 years from 1985 to 2023; the KNDVI vegetation index time-series data consists of annual composite nuclear normalized vegetation index products for 38 years from 1986 to 2023; the climate time-series data includes a 12-month series composed of multi-year monthly averages of potential evapotranspiration from 1995 to 2024 and temperature and precipitation from 2001 to 2024; the solar radiation data consists of annual total solar radiation data for 15 years from 2010 to 2024; and the ground-based measured data are obtained by stratified random sampling, excavating soil profiles down to the mineral soil layer, measuring the vertical thickness of the organic matter layer, and recording environmental information.

3. The remote sensing inversion method according to claim 1, characterized in that, In step S2, the extreme value-aware binning strategy includes: calculating the extreme value threshold based on the interquartile range method, where the upper extreme value threshold is... The lower extreme threshold is The samples are divided into high extreme value samples, low extreme value samples, and normal samples. Normal samples are binned equally according to percentiles, and divided into 10 bins. Extreme value samples are evenly distributed to each bin through a round-robin allocation strategy. Based on bin coding, hierarchical 10-fold cross-validation is used to ensure that the distribution of the target variable in each fold is consistent with the overall dataset.

4. The remote sensing inversion method according to claim 1, characterized in that, In step S3, the quantile transformation uses QuantileTransformer to map the non-normally distributed target value (organic matter layer thickness) to a standard normal distribution; the robust standardization process uses RobustScaler to standardize the eight feature groups based on the median and interquartile range; the encoding embedding encodes the soil type into integers and then maps it to a 16-dimensional continuous vector through an embedding layer; the data preprocessing transformation also includes calculating adaptive weights for samples, assigning dynamic weights of 1.5 to 5.0 times to extreme value samples.

5. The remote sensing inversion method according to claim 1, characterized in that, In step S4, the wavelet transform uses Morlet continuous wavelet transform to extract multi-scale time-frequency features; the principal component analysis dimensionality reduction compresses Landsat wavelet coefficients to 15 dimensions and KNDVI wavelet coefficients to 10 dimensions; the random forest feature selection constructs 200 decision trees with a maximum depth of 8, calculates the importance score of each feature based on variance reduction, selects the top 50 key features, and classifies them into 8 feature groups according to physical attributes: temperature, precipitation, evapotranspiration, solar radiation, topography, aridity, soil properties, and others.

6. The remote sensing inversion method according to claim 1, characterized in that, In step S5, the six-branch parallel hybrid model includes: a climate time series branch, which uses a three-layer one-dimensional convolutional neural network to process monthly time series data of temperature, precipitation, and evapotranspiration, with channel changes from 3->32->64->64; a solar radiation branch, which uses a two-layer one-dimensional convolutional neural network to process 15-year solar radiation time series data, with channel changes from 1->32->64; a Landsat spectral branch, which uses a two-layer fully connected network (MLP) to map 15-dimensional spectral features after wavelet-PCA dimensionality reduction to a 64-dimensional hidden space; a KNDVI vegetation branch, which uses a two-layer fully connected network (MLP) to map 10-dimensional vegetation features after wavelet-PCA dimensionality reduction to a 64-dimensional hidden space; a static environment branch, which uses a two-layer fully connected network (MLP) to process topography, aridity, and soil property features to a 64-dimensional hidden space; and a soil type branch, which uses an embedding layer to map soil type integer encodings to 16-dimensional continuous vectors.

7. The remote sensing inversion method according to claim 1, characterized in that, In step S5, the BiLSTM module is a three-layer bidirectional long short-term memory network with a hidden layer dimension of 64. After bidirectional splicing, it becomes 128-dimensional, which is reduced to 64-dimensional by linear projection, and the dropout rate is 0.

35. The Transformer encoder module contains a 3-layer encoder, with 4 attention heads in each layer. The feedforward network dimension is 256, and the activation function is GELU.

8. The remote sensing inversion method according to claim 1, characterized in that, In step S5, the Multi-Granularity Spatiotemporal Feature Interactive Attention Mechanism (MGSTFIA) includes three parallel channels: climate-vegetation interactive attention, which uses climate features as queries and concatenates climate and vegetation features as key-value pairs. Spectral-terrain interactive attention uses spectral features as queries and concatenates spectral and static features as key-value pairs. Global temporal self-attention is used to perform self-attention calculations on four types of features: climate, vegetation, spectrum, and statics. The outputs of the three channels are concatenated and then fused through a network to obtain 64-dimensional interactive features.

9. The remote sensing inversion method according to claim 1, characterized in that, In step S5, the Feature Group Adaptive Gated Fusion Mechanism (FGAGF) dynamically calculates the weights of five feature groups—climate, solar radiation, spectrum, vegetation, and static attributes—through a Softmax gated network. Each feature group is then transformed by an independent projection network and fused according to its weights. The Uncertainty Aware Dual-Head Prediction Network (UADPN) includes a mean prediction head and a log-variance prediction head, which simultaneously output the predicted mean and prediction variance. The log-variance is pruned and limited to the range of -10 to 10. In step S5, the model training employs an extreme value-aware hybrid loss function, fusing Huber loss, negative log-likelihood loss, and quantile loss, with weights of 0.4, 0.3, and 0.3, respectively, and quantile loss quantiles of 0.1, 0.5, and 0.

9. The AdamW optimizer is used with an initial learning rate of 0.0015 and a weight decay of 0.

08. A course learning strategy is employed, gradually increasing the weight influence of extreme value samples in the first 200 training rounds. Cosine annealing learning rate scheduling is used, with 40 warm-up rounds and a maximum training round count of 600. Gradient norm pruning is employed with a threshold of 1.

0. An early stopping mechanism is used, stopping training if there is no improvement in R² on the validation set for 100 consecutive rounds.

10. The remote sensing inversion method according to claim 1, characterized in that, In step S6, the enhancement during testing employs a 7-cycle prediction. Each time, Gaussian noise with increasing increments is added to the climate, solar radiation, Landsat spectral features processed by wavelet transform and PCA dimensionality reduction, and KNDVI vegetation features. The noise scale is 0.003×i (i=0,1,...,6), and the mean of the 7 predictions is taken as the final normal space prediction value. The inverse transformation restores the normal space prediction value to the original scale through the inverse function of quantile transformation and performs value range clipping on the prediction value. The inversion result is cross-validated with 10-fold cross-validation, and the model determination coefficient R² reaches 0.60, and the root mean square error RMSE reaches 4.91cm.

Images