Crop phenology monitoring method combining sar backscatter and interferometric coherence

Abstract

The present application relates to a kind of crop phenology monitoring methods combined with SAR backscattering and interference coherence, belong to remote sensing image processing and agricultural information technology field.The existing method is mostly dependent on single parameter or simple linear weighting, lack of systematic quantification of the sensitivity difference of backscattering coefficient and interference coherence at different phenological stages, it is difficult to fully tap the complementary potential of both.This application obtains crop time series SAR image, extracts backscattering coefficient time series and interference coherence time series;Construct a joint sensitivity evaluation matrix containing time variation rate, coefficient of variation and independence measure;Calculate time domain weight and spatial domain weight;Introduce coupling coefficient for weighted fusion to generate fusion features;Combined with bare soil and lush vegetation reference backscattering coefficient, construct normalized SAR complex phenology index, determine crop phenological stage by its time series change.Mainly used for continuous monitoring and accurate identification of crop full phenological stage in cloudy and rainy areas.

Description

Technical Field

[0001] This invention belongs to the fields of remote sensing image processing and agricultural information technology, specifically relating to a method for monitoring crop phenology by combining SAR backscattering and interferometric coherence. Background Technology

[0002] Crop phenological monitoring is of great significance for agricultural production management, yield estimation, and food security. Traditional optical remote sensing technology can effectively monitor crop growth under clear sky conditions, but in cloudy and rainy areas, cloud cover severely limits the effective observation window, making it difficult to achieve continuous monitoring of the entire phenological period of crops.

[0003] Synthetic Aperture Radar (SAR), as an active microwave remote sensing technology, possesses all-weather, day-and-night imaging capabilities, providing a new technical approach for crop monitoring in cloudy and rainy areas. Among existing technologies, Nasirzadehdizaji et al., in their 2021 paper "Sentinel-1 interferometric coherence and backscattering analysis for crop monitoring" published in Volume 185 of *Computers and Electronics in Agriculture*, disclosed a method for extracting backscattering coefficients and interferometric coherence from Sentinel-1 satellite data to analyze the variation patterns of maize, sunflower, and wheat at different phenological stages. This study qualitatively indicated the potential for synergistic use of backscattering coefficients and interferometric coherence in crop monitoring by separately plotting time-series curves of these two parameters. Coherence helps estimate the main growth stages, while backscattering provides reliable information throughout the growing season. However, this study only focuses on the qualitative observation and separate description of the temporal variation trends of the two parameters, without providing specific technical means to numerically fuse the backscattering coefficients and interferometric coherence to generate a unified monitoring index.

[0004] As a result, the following three closely related technical problems have arisen in practical applications and have not yet been solved by existing technologies.

[0005] First, in the absence of a quantitative fusion framework, the "complementary" and "redundant" relationship between backscattering coefficients and interference coherence at different phenological stages cannot be systematically characterized and utilized. Their information contributions differ significantly at different crop growth stages; in some stages, backscattering dominates, while in others, coherence dominates, and sometimes the two highly overlap. Existing technologies offer no guidance on how to construct a quantitative fusion model that can adaptively balance this dynamic complementary-redundant relationship.

[0006] Second, due to the lack of demand for quantitative fusion, existing technologies have remained at a qualitative level in analyzing the differences in sensitivity of the two types of parameters at different phenological stages, judging which parameter is "more sensitive" in which stage solely based on the visual changes in time-series curves. This qualitative judgment cannot provide calculable numerical basis for any fusion process that requires weighted decision-making.

[0007] Third, even if an attempt is made to construct a fusion model, determining appropriate fusion weights in the temporal and spatial dimensions remains a key obstacle. In the temporal dimension, the two types of parameters carry varying degrees of information at different phenological stages; in the spatial dimension, the statistical relationship between the two types of parameters varies with local regions. Existing technologies completely fail to address how to adaptively determine fusion weights in these two dimensions. Summary of the Invention

[0008] One objective of this invention is to address the following problem: existing SAR crop phenology monitoring methods rely on a single parameter or only qualitatively compare backscattering coefficients and interferometric coherence, lacking technical means to quantitatively fuse the two numerically. At different phenological stages of crops, the information contributions of the two parameters fluctuate, sometimes complementing and sometimes redundant. Existing technologies cannot construct a quantitative fusion model that can adaptively balance this dynamic complementary-redundant relationship, resulting in incomplete phenological feature extraction and limited monitoring accuracy.

[0009] To achieve the above objectives, this invention provides a method for monitoring crop phenology by combining SAR backscattering and interferometric coherence, comprising the following steps: Step 1: Acquire time-series synthetic aperture radar (SAR) image data covering the entire growth cycle of crops, and extract the backscattering coefficient time series and interferometric coherence time series after preprocessing. Step 2: Based on the backscattering coefficient time series and the interference coherence time series, a joint sensitivity evaluation matrix is ​​constructed for each crop growth stage. The elements of the joint sensitivity evaluation matrix include the time change rate of the backscattering coefficient and the interference coherence, the coefficient of variation, and the independence measure that characterizes the degree of overlap between the two information. Step 3: Based on the joint sensitivity evaluation matrix, calculate the time domain weights and spatial domain weights. The time domain weights are determined by the product of the local information entropy of the backscattering coefficient or the time series of interference coherence and the spatial autocorrelation index. The spatial domain weights are determined by the ratio of the local variance to the covariance of the backscattering coefficient and the interference coherence. Step 4: Using time domain weights and spatial domain weights, the backscattering coefficient and interferometric coherence are weighted and fused, and a coupling coefficient is introduced during the fusion process to characterize the nonlinear interaction between the two, generating fusion features; Step 5: Using the fusion features, combined with the preset backscattering coefficients of bare soil and lush vegetation, a normalized SAR complex phenological index is constructed. By analyzing the temporal changes of this SAR complex phenological index, the crop phenological period is determined.

[0010] Preferably, in step 2, the rate of change of time is the partial derivative of the backscattering coefficient or interference coherence with respect to time during a specific crop growth stage; the coefficient of variation is the ratio of the standard deviation to the mean of the backscattering coefficient or interference coherence during a specific crop growth stage.

[0011] Preferably, in step 2, the independence measure is the difference between 1 and the absolute value of the Pearson correlation coefficient between the backscattering coefficient and the interference coherence.

[0012] Preferably, in step 3, the formula for calculating the local information entropy H is: H = -∑ p(x i log2p(x) i ); where x i p(x) represents the backscattering coefficient or interference coherence within the sliding window. i ) is x i The probability of it appearing within the sliding window.

[0013] Preferably, in step 3, the formula for calculating the spatial domain weight Ws is: Ws = |Var(σ°) - Var(γ)| / (|Cov(σ°, γ)| + ε); where Var(σ°) is the variance of the backscattering coefficient within the local spatial window, Var(γ) is the variance of the interference coherence within the local spatial window, Cov(σ°, γ) is the covariance of the backscattering coefficient and the interference coherence within the local spatial window, and ε is a preset positive constant.

[0014] Preferably, in step 4, the formula for calculating the fusion feature Ffusion is: Ffusion = Wt · σ° + Ws · γ + λ · σ° · γ; where Wt is the time domain weight, Ws is the spatial domain weight, σ° is the backscattering coefficient, γ is the interference coherence, and λ is the coupling coefficient.

[0015] Preferably, in step 5, the formula for constructing the normalized SAR complex climate index SCPI is: SCPI = α · Φstab + β · (σ° - σ°soil) / (σ°veg - σ°soil); where Φstab is the phase stability index in the fusion features, σ° is the current backscattering coefficient, σ°soil is the preset bare soil baseline backscattering coefficient, σ°veg is the preset lush vegetation baseline backscattering coefficient, and α and β are normalization coefficients, satisfying α + β = 1.

[0016] Preferably, before determining the crop phenological period in step 5, the method further includes: using a dynamic time warping algorithm to align and match the time series curve of the constructed SAR complex phenological index with the standard phenological period time series curve of the crop established based on historical observation data, so as to adaptively correct the time offset of the phenological period.

[0017] Preferably, after performing alignment matching using the dynamic time warping algorithm, the method further includes: calculating the time series similarity score after aligning the time series curve of the SAR complex phenological index with the time series curve of the standard phenological period; when the time series similarity score is higher than a preset threshold, outputting the corrected phenological period determination result.

[0018] Preferably, the crop is rice, and the synthetic aperture radar (SAR) image data is obtained from the Sentinel-1 satellite.

[0019] Preferably, in step 4, the value of the coupling coefficient λ is determined through the following iterative process: setting an initial value for λ and calculating the corresponding fusion feature Ffusion; calculating the first mutual information between the fusion feature and the backscattering coefficient time series, and the second mutual information between the fusion feature and the interference coherence time series; using the sum of the first mutual information and the second mutual information as the objective function, using an optimization algorithm to iteratively adjust λ until the objective function value is maximized, and using the λ at this time as the determined coupling coefficient.

[0020] Preferably, the calculation of the first mutual information and the second mutual information adopts an adaptive estimation algorithm based on k-nearest neighbors, specifically including: constructing a two-dimensional dataset by combining the time-series data points of the fusion feature Ffusion with the time-series data points of the backscattering coefficient, or by combining the time-series data points of the fusion feature Ffusion with the time-series data points of the interferometric coherence; calculating the Chebyshev distance between each data point in the two-dimensional dataset and its k-th nearest neighbor; counting the number of data points falling into a rectangular neighborhood centered on each data point and with the Chebyshev distance corresponding to that data point as half the side length; and calculating the first mutual information or the second mutual information using the number of data points falling into the rectangular neighborhood and the digamma function.

[0021] Preferably, the sum of the first mutual information and the second mutual information is used as the objective function, specifically including: calculating the first entropy of the backscattering coefficient time series and calculating the second entropy of the interference coherence time series; and using the ratio of the first mutual information to the first entropy plus the ratio of the second mutual information to the second entropy as the objective function.

[0022] Preferably, after extracting the interference coherence time series in step 1 and before step 2, the method further includes: obtaining the interference coherence of the flooded empty field period before rice transplanting as the reference water surface coherence; calculating the time series canopy coverage based on the backscattering coefficient time series, whereby the canopy coverage characterizes the proportion of the rice canopy obstructing the water surface; obtaining the canopy structure coherence time series through coherence decomposition using the reference water surface coherence, the time series canopy coverage, and the interference coherence time series; and using the canopy structure coherence time series as a substitute interference coherence time series for executing steps 2 to 5.

[0023] Preferably, the coherence time series of the canopy structure is obtained through coherence decomposition, specifically as follows: using a complex coherence model, with canopy coverage as the weight, the complex coherence coefficient corresponding to the interference coherence time series is expressed as a weighted sum of the complex coherence coefficient of the water surface and the complex coherence coefficient of the canopy; the complex coherence coefficient of the water surface corresponding to the reference water surface coherence is substituted as a known quantity, and the weighted sum is solved time-by-time using the time series canopy coverage to obtain the time series of the canopy complex coherence coefficient; the modulus of the canopy complex coherence coefficient time series is taken to obtain the time series of the canopy structure coherence.

[0024] The present invention has at least the following beneficial effects: First, by constructing a joint sensitivity evaluation matrix that includes time-varying rate of change, coefficient of variation, and independence measure, this study, for the first time, elevates the sensitivity differences of backscattering coefficient and interferometric coherence at different phenological stages from a qualitative description to a quantitative measure. The diagonal elements of the matrix directly provide comparable sensitivity values ​​for each parameter, while the off-diagonal independence measure indicates the degree of information overlap between parameters. This provides an objective numerical basis for subsequent fusion weight allocation, avoiding the blindness of existing technologies that rely on empirical judgment.

[0025] Second, by constructing a spatiotemporal adaptive weighted fusion model, the problem of existing technologies being unable to dynamically adjust fusion weights in both the temporal and spatial dimensions is solved. The temporal domain weights are determined by the product of local information entropy and spatial autocorrelation index, enabling the model to automatically identify and rely on parameters that are information-rich and spatially reliable within the current time window. The spatial domain weights are determined by the ratio of variance to covariance, enabling the model to automatically enhance fusion when complementarity is strong and automatically suppress repetition when redundancy is high. This dual adaptive mechanism systematically achieves a balance between the dynamic complementarity and redundancy relationship between backscattering coefficients and interferometric coherence.

[0026] Third, by introducing a coupling coefficient to characterize the nonlinear interaction effect between the backscattering coefficient and the interference coherence, the limitations of simple linear weighting in existing technologies are overcome. The product coupling term in the fusion feature enables the model to capture the complex physical relationship between the two parameters. Experimental results show that the coupling coefficient is not zero, verifying the objective existence of nonlinear interaction and the necessity of introducing this term.

[0027] Fourth, by constructing a normalized SAR complex phenological index (SCPI), the phase stability index and the normalized backscattering coefficient are organically integrated, achieving for the first time the synergistic expression of crop "structural morphology" changes and "physiological state" changes within the same index. The SCPI value ranges from [0,1], possessing clear threshold segmentation agronomic significance. Users do not need to manually switch between multiple curves for interpretation; they can directly determine the phenological stage of the crop based on the index value. Experiments show that the coefficient of determination R between the SCPI and the ground-measured leaf area index is [missing information]. 2 The value reached 0.87, which is lower than the single backscattering coefficient model (Ri). 2 =0.71) and the single coherence model (R 2 =0.65) represents an increase of approximately 22.5% and 33.8%, respectively.

[0028] Fifth, by introducing a dynamic time warping algorithm to correct phenological time shifts, the problem of fixed calendar determination failing when phenological periods are advanced or delayed due to climate anomalies is solved. Experiments show that when the heading period is delayed by about 10 days, the overall accuracy of phenological period identification can be improved from 78% to 92% after dynamic time warping alignment correction.

[0029] Other advantages, objectives and features of the present invention will become apparent in part from the following description, and in part from those skilled in the art through study and practice of the invention. Attached Figure Description

[0030] Figure 1 This is the overall flowchart of the crop phenology monitoring method combining SAR backscattering and interferometric coherence of the present invention; Figure 2 This is a schematic diagram illustrating the distribution pattern of SAR parameters of crops at different phenological stages in Embodiment 1 of the present invention; Figure 3 This is a schematic diagram illustrating the construction principle of the SAR complex phenotype index (SCPI) in Embodiment 1 of the present invention; Figure 4 This is a schematic diagram of the DTW timing alignment algorithm in Embodiment 1 of the present invention; Figure 5 This is a comparison chart of the accuracy of SCPI in Embodiment 1 of the present invention with that of traditional methods. Detailed Implementation

[0031] The present invention will be further described in detail below with reference to examples, so that those skilled in the art can implement it based on the description.

[0032] Example 1 This embodiment selects a rice-growing area in a cloudy and rainy region in southern my country as the study area. This region has an average annual precipitation exceeding 1500 mm, frequent cloud cover, and less than 30% of the time is effectively observed using traditional optical remote sensing, making it a typical area difficult to monitor with optical remote sensing. The rice-growing area in the study area is approximately 500 square kilometers, encompassing both early and late rice growing seasons.

[0033] Step 1: Multi-temporal SAR data preprocessing and basic parameter extraction This embodiment acquired Sentinel-1A / B time-series SAR images covering the entire rice growth cycle (from sowing to harvest, approximately 120 days), with time intervals of 6-12 days, for a total of 20 images. Figure 1 As shown, the data preprocessing process is as follows: (1) Multi-temporal precise registration: The intensity cross-correlation matching method is adopted, with the first image as the master image, and all slave images are registered at the sub-pixel level, with a registration accuracy better than 0.1 pixels. Precise registration is a key prerequisite for ensuring the accuracy of subsequent interferometric coherence calculations - registration error will directly lead to a lower coherence estimate, thereby affecting the reliability of the entire phenological monitoring chain.

[0034] (2) Backscattering coefficient extraction: Radiometric calibration and topographic correction were performed on the registered SAR images, and the digital quantization values ​​(DN values) were converted into backscattering coefficients σ° (unit: dB). Lee filtering (window size 7×7) was used to suppress speckle noise. Thus, the time series backscattering coefficients σ°(t) covering the complete growth cycle of rice were obtained, where t is the acquisition time of each image, for a total of 20 time points.

[0035] (3) Interference coherence calculation: Image pairs with a time baseline of less than 24 days are selected to generate interferometric pairs. The selection of the time baseline needs to balance two factors: an excessively long time baseline will lead to severe temporal incoherence and will not be able to effectively reflect crop structure changes; an excessively short time baseline will result in an insufficient number of interferometric pairs and incomplete temporal coverage. In this embodiment, 24 days is selected as the time baseline threshold, and a total of M interferometric pairs are generated. For each interferometric pair, Goldstein filtering (filter strength α=0.5) is used to suppress phase noise, and then the interference coherence γ is calculated through a 15×15 pixel sliding window. For the interferometric pair consisting of the master image S1 and the slave image S2, the formula for calculating its interference coherence γ is: γ = |Σ (S1,i · S2,i*)| / √( Σ|S1,i| 2 · Σ|S2,i| 2 (1); In the formula, S1,i and S2,i are the complex values ​​(including amplitude and phase information) of the i-th pixel in the master and slave images within the sliding window, respectively. * represents the complex conjugate operation, and Σ represents the summation over all N pixels in the sliding window. The range of the interferometric coherence γ is [0,1]. The closer γ is to 1, the higher the phase consistency between the two images, reflecting a more stable ground structure; the lower γ is, the more significant the changes in the ground surface (such as canopy structure reorganization caused by rapid crop growth).

[0036] (4) Generation of temporal coherence: such as Figure 2 As shown, the coherence calculation results for all interferometric pairs are time-stacked. For each interferometric pair, its effective time is recorded as the midpoint of the master-slave image time for that pair. This yields the interferometric coherence time series γ(t) corresponding to the backscattering coefficient time series on the time axis.

[0037] (5) Geocoding: SRTM 30m DEM is used for geocoding to convert SAR data from slant range coordinate system to WGS84 geographic coordinate system with an output resolution of 10m to ensure spatial registration with ground sampling points.

[0038] Existing techniques (such as Nasirzadehdizaji et al., 2021) can also perform the aforementioned data preprocessing and extraction of σ° and γ. This step shares similarities with existing techniques in that both use Sentinel-1 time-series data and extract backscattering coefficients and interference coherence as fundamental parameters. However, existing techniques, after obtaining the σ° and γ time series data, only perform descriptive statistics and qualitative comparisons, lacking a subsequent quantitative fusion mechanism.

[0039] Step 2: Quantification of SAR parameters sensitivity during crop phenological stages and construction of joint evaluation matrix Based on the rice phenological stage tags (including transplanting, tillering, jointing, heading, grain-filling, and maturity stages) obtained from ground surveys, the time-series data of γ(t) and σ°(t) corresponding to each phenological stage were extracted, and the following three levels of statistical analysis were performed: (1) Distribution pattern statistics: The mean, standard deviation, maximum and minimum values ​​of γ and σ° at each phenological stage were calculated. The results showed that σ° continued to increase from the tillering stage to the heading stage (from -15dB to -8dB), which was closely related to the increase in leaf biomass; γ showed a significant decrease at the heading stage (from 0.6 to 0.3), reflecting the structural instability caused by rapid stem growth.

[0040] (2) Sensitivity index calculation: For each crop growth stage (phenological stage), two basic sensitivity indices are calculated: Rate of change over time: This is the partial derivative of the backscattering coefficient σ° or the interference coherence γ with respect to time during a specific crop growth stage, i.e., ωγ / ωt and ωσ° / ωt. This index reflects the rate of change of the parameter during that phenological stage; the faster the change, the more sensitive the parameter is to the phenological transition at that stage.

[0041] Coefficient of variation: The ratio of the standard deviation to the mean of the backscattering coefficient σ° or interference coherence γ during a specific crop growth stage, i.e., Cv = σ / μ. This index reflects the dispersion of the parameter during phenological periods; a larger Cv indicates higher sensitivity.

[0042] (3) Construction of joint sensitivity evaluation matrix: Construct a 2×2 joint sensitivity evaluation matrix Msens: Msens = [ [ (ωσ° / ωt)·Cv,σ°, 1 - |r| ], [ 1 - |r|, (ωγ / ωt)·Cv,γ] ] (2) In the formula, r is the Pearson correlation coefficient between σ° and γ in the same phenological stage, and 1-|r| is the independence measure characterizing the degree of information overlap between the two. When |r| is close to 1, the independence measure approaches 0, indicating that the two parameters are highly correlated and their information overlaps significantly; when |r| is close to 0, the independence measure approaches 1, indicating that the two parameters are independent of each other and have good complementary potential.

[0043] The diagonal elements of the Msens matrix represent the combined sensitivity of σ° and γ during a specific phenological period (integrating both the rate of change and volatility dimensions), while the off-diagonal elements measure the independence between parameters. By comparing the values ​​of the diagonal elements, the sensitivity advantage of each parameter during a specific phenological period can be quantitatively determined.

[0044] The analysis results of this embodiment show that the diagonal element value corresponding to γ ​​in Msens during the heading stage is significantly higher than that of σ°, indicating that γ sensitivity is dominant; while the diagonal element value corresponding to σ° is dominant during the grain-filling stage. This quantitative conclusion is corroborated by the qualitative observation of the distribution pattern statistics.

[0045] Existing techniques (such as Nasirzadehdizaji et al., 2021) only qualitatively determine "which parameter changes significantly at which stage" by visually observing the fluctuations of time-series curves. They cannot provide comparable or calculable sensitivity values, nor can they distinguish between sensitivity types "dominated by rapid response" and "dominated by highly discrete factors." This step, by introducing two independent dimensions—the time partial derivative and the coefficient of variation—and constructing a joint sensitivity evaluation matrix that includes independence measures, elevates sensitivity analysis from "qualitative description" to "quantitative measurement" for the first time. This provides an objective and calculable numerical basis for determining the dynamic fusion weights in subsequent steps.

[0046] Step 3: Construction of Spatiotemporal Adaptive Weighted Feature Fusion Model Based on the principles of information entropy and spatial autocorrelation, a spatiotemporal adaptive weighted fusion model is constructed. The core innovation of this model lies in the fact that instead of using fixed weights, the model automatically adjusts the fusion weights dynamically according to the information richness, spatial reliability, and complementary-redundant relationship between σ° and γ at the current phenological stage.

[0047] (1) Calculation of time domain weights: First, calculate the local information entropy H of the σ° and γ time series data. For the set of SAR parameter time series values ​​X = {x1, x2, ..., x_n} within the sliding window, the formula for calculating its local information entropy is: H = - Σ p(x i log2p(x) i (3); In the formula, the numerical range of the backscattering coefficient or interference coherence within the sliding window is divided into several intervals, x i Let p(x) be the representative value of the i-th interval. i ) represents the frequency with which a value falls into the i-th interval within the sliding window (i.e., the proportion of samples in that interval to the total number of samples in the window). Information entropy H reflects the ability of a feature to distinguish phenological changes within a local time window: the higher the entropy value, the more dispersed the value of the parameter is within that time period, and the richer the phenological state information it contains.

[0048] Secondly, the spatial autocorrelation index—Moran's I index—is calculated. This index assesses the spatial clustering of features: a higher value indicates that the features exhibit a significant spatial clustering distribution, a low probability of local noise interference, and that the features are reliable.

[0049] The information entropy H is multiplied by Moran's I exponent and then normalized across parameters to obtain the time domain weight wt. In this embodiment, the H×Moran's I values ​​are calculated for σ° and γ respectively, and then normalized so that wt(σ°) + wt(γ) = 1. The physical meaning of the time domain weight is: within the current local time window, the model assigns a higher degree of confidence to the parameter that carries richer information and is more spatially reliable in the time dimension.

[0050] (2) Spatial domain weight calculation: Spatial domain weights ws are defined based on feature variance and covariance: ws = |Var(σ°) - Var(γ)| / ( |Cov(σ°, γ)| + ε ) (4); In the formula, Var(σ°) and Var(γ) are the normalized variances of the backscattering coefficient and interference coherence within the local spatial window, respectively; Cov(σ°, γ) is the normalized covariance of the two; and ε is a preset positive constant (in this embodiment, ε = 10). -6 To prevent the denominator from being zero, the standardization process converts both the backscattering coefficient and the interference coherence into dimensionless variables with a mean of 0 and a standard deviation of 1, thereby eliminating the influence of the difference in dimensions between the two on the comparability of their variances.

[0051] The design principle of spatial domain weights is as follows: when σ° and γ change asynchronously within the local spatial window (the absolute value of the covariance is small, that is, the two are highly independent and complementary), the numerator is relatively large and the denominator is small, ws increases, and the model automatically enhances the utilization of complementary information; when the two change highly synchronously (the absolute value of the covariance is large, that is, the redundancy is high), ws decreases, suppressing the repeated inclusion of redundant information.

[0052] (3) Calculation of phase stability index: To suppress the interference of interferometric phase noise on the fusion result, a phase stability index Φstab is introduced: Φstab = e^(-σφ) (5); In the formula, σφ represents the standard deviation of the interference phase within the local spatial window. The smaller σφ is (the more stable the phase), the closer Φstab is to 1; the larger σφ is (the more severe the phase noise), the closer Φstab is to 0. This index is used to characterize the structural stability of the crop canopy in the subsequent SCPI construction.

[0053] (4) Feature fusion generation: The time-domain weights, spatial-domain weights, and original parameters are weighted and fused together, and a coupling coefficient is introduced to characterize the nonlinear interaction between σ° and γ, generating the fused feature Ffusion: Ffusion = wt·σ° + ws·γ + λ·σ°·γ (6); In the formula, wt is the time domain weight, ws is the spatial domain weight, σ° is the backscattering coefficient, γ is the interference coherence, and λ is the coupling coefficient.

[0054] Formula (6) consists of three terms: the first term wt·σ° is the time-weighted backscattering contribution; the second term ws·γ is the spatially weighted coherence contribution; and the third term λ·σ°·γ is the product coupling term, which is used to characterize the special interaction effect produced when σ° and γ are at high or low values ​​at the same time, breaking through the limitations of simple linear weighting.

[0055] In this embodiment, the coupling coefficient λ is optimized by minimizing the mean square error (MSE) between the fused features and the measured phenological period labels using the gradient descent algorithm. The optimized result is λ=0.15. The fact that λ is not zero quantitatively proves that there is indeed a nonlinear interaction between σ° and γ, verifying the necessity of introducing the coupling coefficient.

[0056] Existing techniques, such as the work of Nasirzadehdizaji et al. (2021) published in Computers and Electronics in Agriculture (Vol. 185, document number 106118), entitled "Sentinel-1 interferometric coherence and backscattering analysis for crop monitoring" (hereinafter referred to as Nasirzadehdizaji et al., 2021), only plot the time-series curves of σ° and γ separately for qualitative comparison, without the concept of "weights," let alone adaptive weights and nonlinear coupling. Even if weighted fusion is envisioned, it cannot solve the core problem of "how to determine the weights." This step constructs time-domain weights through the product of information entropy × Moran's I, constructs spatial-domain weights through the variance / covariance ratio, and captures nonlinear interactions through coupling coefficients. For the first time, it establishes a complete, data-driven, adaptive fusion framework that does not require human pre-setting, systematically solving the fundamental technical problem of "how to adaptively balance complementary-redundant relationships."

[0057] Step 4: Construction of SAR Complex Phenotyping Index (SCPI) like Figure 3 As shown, based on the fusion model in step 3, the SAR complex phenological index SCPI is constructed. The design goal of SCPI is to fuse the phase stability index (reflecting structural dynamics) and the normalized backscattering coefficient (reflecting physiological dynamics) to achieve a synergistic characterization of the dual phenological response of "structure-physiology".

[0058] The construction formula for SCPI is as follows: SCPI = α·Φstab + β·(σ° - σ°soil) / (σ°veg - σ°soil) (7); The meanings of each term in the formula are as follows: Φstab is the phase stability index defined in equation (5) of step 3, which characterizes the stability of the crop canopy structure. During periods of drastic morphological changes in crops (such as rapid stem elongation during the heading stage), structural instability leads to increased interference phase noise, and Φstab decreases significantly.

[0059] σ° is the backscattering coefficient at the current moment.

[0060] σ°soil is the baseline backscattering coefficient of bare soil, which is the average backscattering coefficient of images taken from the flooded, empty field period before sowing. In this example, σ°soil = -16.5dB.

[0061] σ°veg is the baseline backscattering coefficient for lush vegetation, taken as the average backscattering coefficient of images from the heading stage (when canopy coverage is at its highest). In this example, σ°veg = -7.8dB.

[0062] (σ° - σ°soil) / (σ°veg - σ°soil) is the normalized backscattering coefficient. Range normalization is performed using the bare soil period and the lush vegetation period as dual benchmarks, mapping the original σ° to the [0,1] interval. This normalization method enables SCPI to have cross-temporal and spatial generalization capabilities—SCPI values ​​from different regions and different planting seasons are comparable.

[0063] α and β are normalization coefficients that satisfy α + β = 1, controlling the contribution ratios of the structural stability term and the physiological state term in SCPI, respectively. In this embodiment, sensitivity analysis determined α = 0.45 and β = 0.55, slightly favoring the physiological state term. This is consistent with the physical characteristics of rice growth, which exhibits a wide response range and high signal stability of the backscattering coefficient.

[0064] The physical significance of SCPI: The first term, α·Φstab, captures changes in crop "structural morphology"—such as canopy structure reorganization caused by jointing and heading; the second term, β·(normalized σ°), captures changes in crop "physiological characteristics"—such as increases in water content and biomass caused by leaf growth. Together, they achieve a complete expression of phenological information.

[0065] The SCPI value range is [0,1]. This embodiment establishes the following threshold segmentation standard based on ground phenological observation data: SCPI < 0.2: corresponds to bare soil or fallow period (before transplanting); 0.2 ≤ SCPI<0.5: Corresponds to the vegetative growth stage (from transplanting to tillering and jointing); 0.5 ≤ SCPI < 0.7: Corresponds to the reproductive growth stage (heading to grain filling); SCPI ≥ 0.7: corresponds to maturity stage.

[0066] Existing SAR phenological indices (such as Xin Bao, Rui Zhang, Xu He, Age Shama, Gaofei Yin, Jie Chen, Hongsheng Zhang, Guoxiang Liu, Xianjian Shi, A novel dual-polarization SAR vegetation index for crop phenology detection, Computers and Electronics in Agriculture, Volume 239, Part B, 2025) (referred to as the polarization-decomposition-based DRVIs index, 2025) propose a novel dual-polarization SAR vegetation index (DRVIs) that combines parameters only within a single physical domain of backscattering, without using interferometric coherence. While Nasirzadehdizaji et al. (2021) possess both σ° and γ curves, they did not construct any unified index. SCPI, for the first time, incorporates "structural stability" (derived from the phase information of γ) and "physiological state" (derived from the intensity information of σ°) into the same weighted summation index framework, achieving synergistic characterization across physical domains. The threshold segmentation standard gives SCPI a clear agronomical physical meaning, allowing it to be directly used for operational phenological period division.

[0067] Step 5: SCPI accuracy verification and timing alignment optimization (1) Accuracy verification: like Figure 4 and Figure 5 As shown, 50 ground sampling points were set up in the study area to simultaneously acquire rice leaf area index (LAI) and phenological observation data. A regression model between SCPI and LAI was established: LAI = a·e^(b·SCPI) + c (8); The regression results for this example are: a=2.35, b=1.82, c=0.28, and the coefficient of determination R0. 2 =0.87, RMSE=0.42.

[0068] For comparison, regression models (R0) were established for single σ° and LAI, respectively. 2 =0.71) and the regression model of single γ with LAI (R 2 =0.65). Compared to the single σ° model, R of SCPI... 2 It improved by approximately 22.5%; compared to the single γ model, R 2The improvement was approximately 33.8%. This significant improvement verifies the synergistic gain effect of fusing backscattering and interferometric coherence, and also demonstrates the substantial progress made by the spatiotemporal adaptive weighted fusion framework and SCPI index proposed in this invention in terms of phenological inversion accuracy.

[0069] (2) Timing alignment optimization: In response to the situation in the study area in 2023 where the heading period of rice was delayed by about 10 days due to continuous rainy weather, the Dynamic Time Warping (DTW) algorithm was used for time alignment correction.

[0070] First, standard phenological SCPI curves were constructed—based on SCPI time-series data of the same study area and the same rice variety from the previous five years (2018-2022), and synthesized by averaging according to phenological stages.

[0071] Then, the DTW distance matrix between the SCPI time series curve and the standard curve for the current season of 2023 is calculated. The DTW algorithm uses dynamic programming to find the optimal alignment path between the two curves, achieving flexible matching of the time axis—that is, allowing a certain time point on the current year's curve to be mapped to a different time point on the standard curve, thereby automatically correcting the overall shift in phenological periods caused by climate anomalies.

[0072] After alignment, calculate the temporal similarity score Dsim: Dsim = 1 - DTWdist / DTWmax (9); In the formula, DTWdist is the cumulative normalized distance between the two curves after DTW alignment, and DTWmax is the normalization factor (the product of the sum of the lengths of the two curves and their respective diagonal distances). The value of Dsim ranges from [0,1], and the closer it is to 1, the more similar the two curves are.

[0073] In this embodiment, Dsim=0.85 was calculated, which is higher than the preset threshold of 0.8. This indicates that the SCPI time series curve still closely matches the historical standard curve after DTW correction under the condition of phenological period shift, and the correction result has high reliability.

[0074] The core evaluation metric for DTW alignment effectiveness—phenological period identification accuracy—shows that before alignment, due to a delay of approximately 10 days in the heading stage, only 78% of ground sampling points were correctly identified if phenological periods were determined using a fixed calendar. After DTW alignment correction, the overall accuracy of phenological period identification improved to 92%, a net increase of 14 percentage points. This result demonstrates that the synergistic effect of DTW and SCPI effectively solves the problem of dynamic phenological period shifts that existing technologies cannot address.

[0075] Existing technologies (such as Nasirzadehdizaji et al., 2021) completely ignore the issue of phenological shift correction—they implicitly assume that phenological periods are fixed. Existing DTW remote sensing applications (such as rice planting area extraction studies based on optimized DTW) align the original backscattering coefficient curves for crop classification rather than phenological shift correction. This embodiment is the first to use DTW to align the fused SCPI time-series curves, correcting the prominent problem of "phenological period shift" in actual agricultural production. Because the SCPI curve incorporates complementary information from σ° and γ, its curve shape is more robust and its phenological inflection points are clearer than those of single-parameter curves, thus the DTW alignment effect is superior to aligning any single-parameter curve. The final recognition accuracy of 92% is significantly improved compared to traditional single-parameter methods (72%, 68%), fully verifying the overall advancement and practicality of the technical solution of this invention.

[0076] Example 2 This embodiment addresses the typical scenario of sparse ground phenological label data in the aforementioned study area by designing an unsupervised coupling coefficient optimization scheme. This scenario is common in the complex terrain and fragmented, cloudy and rainy areas of southern my country. Due to inconvenient transportation and high labor costs, ground phenological observation data can only be obtained from a very limited number of sample points, making it impossible to support the supervised optimization solution of λ in Embodiment 1.

[0077] I. In this embodiment, only 5 out of 50 ground sampling points in the study area (10%) have complete phenological period label data. The remaining 45 points are only used for final accuracy verification and do not participate in the determination of λ. This setting simulates the typical situation of extremely scarce ground data in actual operational monitoring.

[0078] In this scenario, the gradient descent + MSE supervised optimization method used in Example 1 faces serious difficulties due to the insufficient number of tag samples: with only tag information from 5 sampling points, the λ value obtained by MSE optimization is extremely unstable—the absence or anomalies of key sampling points during the heading stage cause λ to fluctuate significantly between 0.05 and 0.42, failing to converge to a reliable value. If this unstable λ value is forcibly adopted, the overall accuracy of subsequent SCPI index construction and phenological stage determination will significantly decrease. The performance degradation of this technique in the sparse tag scenario is not due to a defect in the method itself, but rather to the fact that its design premise (requiring sufficient ground tags) cannot be met in this scenario.

[0079] II. Comparison Scheme Setup To verify the effectiveness of the technical solution, the following comparison scheme is set up in this embodiment: Scheme A (Comparative Example 2-1): With a fixed empirical value of λ=0, this degenerates into a linear weighted fusion using only the weights, ignoring the nonlinear interaction term between σ° and γ. This scheme simulates the simplification strategy most likely adopted by those skilled in the art when faced with insufficient data—abandoning optimization and directly ignoring the nonlinear term.

[0080] Scheme B (Comparative Example 2-2): This scheme employs the supervised optimization method from Example 1, using only 5 labeled samples to minimize the MSE (Mean Sequence Size) λ. This scheme simulates the actual performance of forcibly using a supervised method under sparse label conditions, demonstrating the limitations of the supervised method in this scenario.

[0081] Scheme C (Example 2-1): Unsupervised optimization to determine λ is employed by maximizing mutual information. The objective function is the sum of the first mutual information between the fusion feature Ffusion and the backscattering coefficient time series, and the second mutual information between the fusion feature and the interference coherence time series. The initial value of λ is set to 0.10. A grid search algorithm is used to iteratively search within the interval [0, 1] with a step size of 0.01. In each iteration, the mutual information of Ffusion corresponding to the current λ and its two mutual information with σ° and γ time series is calculated, and the λ value that maximizes the sum of the two mutual information is selected.

[0082] In this embodiment, the mutual information estimation adopts an adaptive estimation algorithm based on k-nearest neighbors: A two-dimensional dataset is constructed by combining Flux time-series data points and σ° time-series data points (or γ time-series data points); the nearest neighbor number k=10, and the Chebyshev distance between each data point and its 10th nearest neighbor is calculated; the number of data points falling within a rectangular neighborhood centered on each data point and with its corresponding Chebyshev distance as half a side is counted; using this number and the digamma function, the mutual information value is calculated according to the k-nearest neighbor mutual information estimation formula. Specifically: Step 1: Construct a two-dimensional dataset. For two sets of time-series data for which mutual information needs to be calculated—the time-series values ​​of the fusion feature Ffusion {F1, F2, ..., F_n}, and the time-series values ​​of the backscattering coefficient σ° {σ°1, σ°2, ..., σ°_n} (or the time-series values ​​of the interference coherence γ {γ1, γ2, ..., γ_n})—pair their data points at the same time to construct a two-dimensional dataset Z = {z1, z2, ..., z_n}, where z_i = (F_i, σ°_i) (or z_i = (F_i, γ_i)), and n is the time-series length.

[0083] Step 2: Calculate the Chebyshev distance for each data point. Set the number of nearest neighbors k (k=10 in this example). For each data point z_i in the two-dimensional dataset Z, calculate its Chebyshev distance to all other data points z_j (j≠i) in the dataset. The Chebyshev distance between two points z_i = (x_i, y_i) and z_j = (x_j, y_j) is defined as: d_Chebyshev(z_i, z_j) = max( |x_i - x_j|, |y_i - y_j| ); Find the kth smallest distance value from z_i to all other data points, denoted as ε_x(i) (when calculating the first dimension x, i.e., the fusion direction) or ε_y(i) (when calculating the second dimension y, i.e., the σ° or γ direction).

[0084] Specifically, ε_x(i) is the projection component of the Chebyshev distance along the x-axis among the k nearest neighbors found with z_i as the center and d_Chebyshev(z_i, z_j) as the metric. Since the Chebyshev distance takes the maximum value of the difference between the two axes, ε_x(i) is the larger of the absolute value of the difference between the k nearest neighbor and z_i along the x-axis and the absolute value of the difference along the y-axis.

[0085] Step 3: Count the number of data points in the neighborhood. For each data point z_i, within the rectangular neighborhood centered on z_i and with half the length of ε_x(i) (i.e., the x-direction interval [x_i - ε_x(i), x_i + ε_x(i)], and the y-direction interval [y_i - ε_y(i), y_i + ε_y(i)]), count the number of data points n_x(i) and n_y(i) falling into this neighborhood, as well as the number of data points n_xy(i) falling into the intersection of the two directional neighborhoods.

[0086] In fact, according to the k-nearest neighbor-based estimation method: n_x(i) counts the number of points whose distance in the x-direction is less than ε_x(i), and n_y(i) counts the number of points whose distance in the y-direction is less than ε_y(i). Since ε_x(i) and ε_y(i) are both Chebyshev distances corresponding to the k-th nearest neighbor, under the definition of two-dimensional Chebyshev distance, ε_x(i) = ε_y(i) = ε(i).

[0087] Step 4: Calculate mutual information using the digamma function. The estimated value of mutual information is given by the following formula: MI(Ffusion, σ°) = ψ(k) - (1 / k) - (1 / n) · Σ_{i=1}^{n} [ψ(n_x(i)) + ψ(n_y(i))]+ ψ(n); In the formula, ψ(·) is the digamma function (i.e., the first derivative of the logarithm of the Gamma function), k is the number of nearest neighbors selected, n is the total number of time series data points, n_x(i) is the number of data points falling in the x-axis direction of the rectangular neighborhood centered at z_i and with side length ε(i) (excluding z_i itself), n_y(i) is the number of data points falling in the neighborhood in the y-axis direction (excluding z_i itself), and Σ represents the summation over all n data points.

[0088] Similarly, when calculating MI(Ffusion, γ), simply replace σ° with γ in the above formula.

[0089] The digamma function can be calculated using the following recursive relation: ψ(1) = -0.5772156649 (the negative value of Euler's constant γ), ψ(m+1) = ψ(m) + 1 / m, where m is a positive integer. For any positive real number, this recursive relation can be used in conjunction with the asymptotic expansion ψ(x) ≈ ln(x) - 1 / (2x) - 1 / (12x) 2 ) + 1 / (120x 4 The required accuracy is calculated.

[0090] In this embodiment, k=10 and n=20 (corresponding to 20 SAR image sampling times). Taking scheme C as an example, at λ=0.17, MI(Ffusion, σ°) = 1.847 nat and MI(Ffusion, γ) = 1.562 nat are calculated. The sum of the two is 3.409 nat, which is the maximum value among all candidate λ values. Therefore, the optimal λ=0.17 is determined.

[0091] Scheme D (Example 2-2): Based on Scheme C, further calculate the entropy of the σ° time series and the γ time series. The entropy is calculated using the same k-nearest neighbor estimation method as described above: H(σ°) = ψ(n) - ψ(k) + ln(c_d) + (d / n) ·Σ_{i=1}^{n} ln(ε(i)); where d is the data dimension (here d=1, because only the entropy of a single variable is calculated), c_d is the volume of a d-dimensional unit sphere (c_d=2 when d=1), and ε(i) is the distance from the σ° time series data point x_i to its k-th nearest neighbor. H(γ) is calculated similarly. In this embodiment, H(σ°) = 2.485 nat, H(γ) = 2.143 nat. The objective function of scheme D is revised as follows: Objective function = MI(Ffusion, σ°) / H(σ°) + MI(Ffusion, γ) / H(γ) = 1.847 / 2.485 + 1.562 / 2.143 = 0.743 + 0.729 = 1.472; this value reaches its maximum at λ=0.16, thus determining the optimal λ=0.16. The remaining steps are the same as those of scheme C.

[0092] III. Comparison of λ determination results The λ values ​​determined for each scheme are as follows: Table 1. Determined λ values The λ values ​​determined by schemes C and D are both within a reasonable neighborhood of the supervised optimization result (0.15) in Example 1, and the objective function exhibits a single-peak shape with a clear peak during the search process, indicating that the unsupervised optimization criterion has good convergence performance.

[0093] IV. Verification of Phenological Monitoring Accuracy Substituting the λ values ​​determined by each scheme into the fusion feature generation formula (6) in step 4, the SCPI index is constructed. Accuracy verification is performed using the same 50 ground sampling points (including 5 label points and 45 independent verification points) as in Example 1. The evaluation index is the coefficient of determination R between the SCPI and the ground-measured LAI. 2 The root mean square error (RMSE) and the overall accuracy of phenological period identification.

[0094] Table 2 Verification results of phenological monitoring accuracy Scheme A, which fixes λ=0, causes the fusion model to degenerate into a linear weighted model, completely ignoring the nonlinear interaction effect between σ° and γ, which has been verified by physical facts and Example 1 (supervised optimization results with λ=0.15 in Example 1). This simplification results in a significantly weaker response of SCPI during the heading stage—a stage where σ° and γ are both at medium to high values, and the product term contributes the most to fusion—leading to a decrease in phenological identification accuracy of approximately 8 percentage points compared to the scheme of this invention. This gap confirms the necessity of the nonlinear coupling term and also shows that simply ignoring this term cannot solve the problem of λ determination in the absence of labels.

[0095] In Scheme B, with only 5 labeled samples, the supervised optimization yields an unstable λ value. The λ value obtained by randomly selecting different subsets of 5 labels from the same data can fluctuate between 0.08 and 0.37. In this experiment, the λ=0.13 of Scheme B is only the median value of a single sampling. R 2 With a value of only 0.80 and a phenological recognition accuracy of 83%, it is actually lower than the 84% of scheme A, which directly ignores the nonlinear term—this indicates that the noise introduced by the unstable λ outweighs the gain brought by the nonlinear term. This confirms the limitations of supervised methods in sparse label scenarios.

[0096] R of scheme C 2 The accuracy reached 0.86, with a phenological identification accuracy of 91%, which is close to the effect of supervised optimization in Example 1 (Ri). 2 =0.87, accuracy 92%). Scheme D further refines the objective function through entropy normalization, ensuring a fair balance between the information contributions of σ° and γ, R 2 The accuracy was improved to 0.87, reaching 92%, which is on par with the supervised optimization in Example 1—but without relying on any ground phenology labels.

[0097] Existing technologies (such as Nasirzadehdizaji et al., 2021) do not involve the fusion of σ° and γ, therefore there is no problem in solving for λ. Example 1 provides a supervised optimization scheme for λ, but it relies on the construction of a joint sensitivity evaluation matrix and the support of ground phenology labels in step 2. When the acquisition of ground label data is limited, this method cannot be stably solved due to insufficient label samples.

[0098] This embodiment discloses for the first time a fully unsupervised scheme for determining λ: Starting from the physical intuition that "the fused feature should retain the effective information carried by each of the two original parameters to the greatest extent," it introduces mutual information (a standard tool in information theory for measuring the statistical dependency between two variables) as an unsupervised optimization criterion. This allows the model to autonomously perceive "what kind of λ allows the fused feature to most faithfully carry the information of σ° and γ" without any ground labels. Normalized mutual information (Scheme D) further standardizes the mutual information based on the entropy of each parameter itself, avoiding the optimization process being dominated by parameters with higher information content and achieving a fair trade-off.

[0099] The technical effect of this embodiment proves that in scenarios where ground tags are extremely scarce (e.g., only 5 tag samples), the unsupervised λ optimization scheme of the present invention can keep the phenological monitoring accuracy close to the level of supervised optimization in Example 1, while conventional simplified strategies (fixing λ=0) or forcibly using unreliable supervision methods cannot achieve acceptable accuracy.

[0100] Example 3 This embodiment addresses the issue of background coherence pollution caused by flooding when the method of the present invention is applied to rice, a specific crop, and employs a canopy-water surface coherence decomposition scheme. This embodiment uses the same study area and the same set of Sentinel-1 time-series SAR data as Embodiment 1, but focuses on the stage from the early transplanting stage to the tillering stage of rice, when the impact of flooding background is most significant.

[0101] I. Rice in the study area was transplanted using the flooding method. Before transplanting, the fields were plowed and leveled, and a water layer of 5-15 cm was maintained, with the water surface calm. From transplanting until the peak tillering period (about 30-40 days after transplanting), the water layer was maintained in the field, and the rice canopy grew from nothing to something, from small to large, above the water surface.

[0102] SAR interferometric coherence is highly sensitive to surface changes. Calm water surfaces, due to their extremely stable geometry and dielectric properties, exhibit almost no decoherence between interferometric pairs, with theoretical coherence approaching 1.0. This means that before and in the early stages of rice transplanting, the observed interferometric coherence γ_obs primarily reflects the contribution of the water surface rather than the canopy. As the rice canopy gradually unfolds and covers the water surface, the subtle movements of the canopy itself under the influence of wind, growth, and other factors begin to dominate the decoherence process within pixels. The observed coherence is actually a mixture of water surface and canopy contributions; its synthesis mechanism is not a simple intensity weighting but involves complex coherent superposition of the two scatterers in the complex domain.

[0103] The general processing flow of the existing technology or Example 1 does not perform water surface-canopy separation for observation coherence, and directly uses γ_obs as the interference coherence time series for subsequent sensitivity evaluation, fusion and SCPI construction. In the dryland crop scenario, the decoherence contribution of the soil background is relatively stable and low, and this simplification has a limited impact. However, in the rice flooding scenario, the pollution of the water surface with high coherence will systematically increase the γ value in the early stage of transplanting, causing the sensitivity evaluation matrix (2) to distort the sensitivity assessment of γ in this stage - manifested as underestimating the true response amplitude of γ in the key phenological transition from the absence of canopy to its presence, thereby affecting the allocation of spatiotemporal adaptive weights and the ability of the SCPI index to characterize the early phenological stage. This performance decline is not due to the error of the general method itself, but rather because the implicit assumption that "all pixel coherence comes from the canopy" on which it is based does not hold in the rice flooding scenario.

[0104] II. Comparison Scheme Setup To verify the effectiveness of the technical solution, the following comparison scheme is set up in this embodiment: Scheme A (Comparative Example 3-1): Without coherence decomposition, the original observational interferometric coherence γ_obs is directly used as the γ time series, and steps 2 to 5 of Example 1 are executed. This scheme simulates the direct application effect of the general processing flow of Example 1 in a rice flooding scenario.

[0105] Scheme B (Scheme 3 of Example): After extracting the interference coherence time series γ_obs(t) in step 1, the following coherence decomposition steps are performed to obtain the canopy structure coherence time series γ_canopy(t). Then, steps 2 to 5 are performed by replacing γ_obs(t) with γ_canopy(t): Step 1-1: Obtain the coherence of the reference water surface Two Sentinel-1 images were selected before rice transplanting, when the fields had been flooded but the seedlings had not yet been transplanted, with a time baseline of 12 days. During this period, the field surface was pure water (without any canopy obstruction), and the interference coherence of this image pair was calculated as the reference water surface coherence γ_water. In this embodiment, γ_water = 0.94 was measured on multiple pure water surface pixels.

[0106] Steps 1-2: Calculate time-series canopy coverage Based on the obtained backscattering coefficient time series σ°(t), the time series canopy cover f_v(t) is calculated. Canopy cover characterizes the proportion of water surface obstruction by the rice canopy, and its value ranges from [0,1]. This embodiment utilizes the high sensitivity of σ° to water surface and vegetation, and estimates it using the following formula: f_v(t) = (σ°(t) - σ°_water) / (σ°_full - σ°_water); In the formula, σ°_water is the average backscattering coefficient during the flooded empty field period before transplanting (-16.5 dB in this example), and σ°_full is the asymptotic value of the backscattering coefficient after the canopy completely covers the water surface (taking the average value of -7.8 dB at the heading stage). When σ°(t) approaches σ°_water, f_v(t) approaches 0, indicating that the water surface is completely exposed; when σ°(t) approaches σ°_full, f_v(t) approaches 1, indicating that the water surface is completely covered by the canopy. The canopy coverage f_v(t) shows an S-shaped curve that monotonically increases from 0 to approximately 0.95 as the rice grows.

[0107] Steps 1-3: Perform complex coherence decomposition to obtain the coherence time series of the canopy structure. For the interference coherence γ_obs(t) corresponding to each time t, it is first transformed into a complex coherence coefficient Γ_obs(t). Since the definition of interference coherence is known, the modulus of Γ_obs is γ_obs. This invention constructs the following complex coherence model: Γ_obs(t) = f_v(t) · Γ_canopy(t) + (1 - f_v(t)) · Γ_water (10); The physical meaning of this model is as follows: the observed complex coherence coefficient Γ_obs(t) is a complex weighted sum of the canopy complex coherence coefficient Γ_canopy(t) and the water surface complex coherence coefficient Γ_water, with the canopy coverage f_v(t) as the weight. The water surface complex coherence coefficient Γ_water is determined by the reference water surface coherence γ_water (the modulus of Γ_water is γ_water, and the phase is taken as the spatial mean of the water surface interference phase).

[0108] In formula (10), Γ_obs(t) comes from time-series observations, f_v(t) is obtained from steps 1-2, and Γ_water is the obtained baseline value. Substituting f_v(t), Γ_obs(t), and Γ_water into formula (10), the complex coherence coefficient Γ_canopy(t) of the canopy is obtained by solving time-by-time: Γ_canopy(t) = [Γ_obs(t) - (1 - f_v(t)) · Γ_water] / f_v(t) (11); Taking the modulus of Γ_canopy(t) yields the final desired canopy structure coherence time series: γ_canopy(t) = |Γ_canopy(t)| (12); When f_v(t) is extremely small (in the early stage of transplanting), formula (11) is more sensitive to small errors in Γ_water. In this embodiment, the threshold of f_v(t) is set to 0.05: when f_v(t) < 0.05, γ_canopy(t) is directly assigned to γ_obs(t) (because the contribution of the canopy is extremely small at this time, the decomposition is not very meaningful and the value is unstable); when f_v(t) ≥ 0.05, formulas (11) and (12) are executed normally.

[0109] Steps 1-4: Replace and execute Using γ_canopy(t) as a substitute for the interference coherence time series, input steps 2 to 5 of Example 1 to complete the subsequent sensitivity evaluation, spatiotemporal adaptive fusion, SCPI construction and phenological period determination.

[0110] III. Comparison of γ time series during key phenological stages To clearly demonstrate the effect of coherence decomposition, the table below shows a comparison of the mean interference coherence values ​​of scheme A (using γ_obs directly) and scheme B (using γ_canopy) in four typical phenological stages: Table 3 Mean values ​​of interference coherence at each stage The above comparison clearly demonstrates the quantitative impact of water surface pollution: during the critical phenological transition period from tillering to jointing, γ_obs was systematically overestimated by approximately 0.14-0.09 (relative deviation of approximately 24%-20%) due to the lifting effect of high water surface coherence. This stage is precisely the critical window through which the sensitivity matrix should capture the response of γ to changes in canopy structure; overestimation of γ directly leads to anomalies in the matrix's diagonal elements—severely weakening the characterization role of γ in the early vegetative growth stage.

[0111] IV. Verification of Phenological Monitoring Accuracy The SCPI indices obtained from schemes A and B were compared and verified with the measured LAI and phenological period labels from 50 ground sampling points: Table 4 Comparison of LAI and phenological label results Results Analysis: The SCPI of Scheme B and the R of LAI 2 The accuracy of phenological identification improved from 0.81 to 0.88, and from 87% to 93%. The gain mainly comes from the significant improvement in the early phenological stage (from transplanting to tillering).

[0112] The identification accuracy during the transplanting-tillering stage jumped from 79% to 92%, an improvement of 13 percentage points. This is because in Scheme A, γ_obs was contaminated by water surface coherence during this stage, causing γ to lose its inherent sensitivity to changes in canopy structure. After decomposition, γ_canopy regained this ability, allowing the sensitivity matrix to correctly assess the relative advantage of γ during this stage, thereby optimizing the allocation of spatiotemporal adaptive weights. This stage is the critical physiological period for rice to transition from "seedlings in water" to "independent plants," and accurate identification has significant guiding significance for field water and fertilizer management. The improvement during the jointing to maturity stage was within 3 percentage points. This is because as canopy coverage increases, the contribution of water surface naturally weakens, and γ_obs and γ_canopy tend to be consistent, so whether or not decomposition is performed has little impact on the results. This pattern is completely consistent with the physical explanation of the coherence mixing mechanism in this embodiment.

[0113] The general method in Example 1 performed well when processing dryland crops such as corn and wheat. Its scheme is based on the conventional physical assumption that "intra-pixel decoherence mainly comes from canopy changes," and its performance has been verified in dryland scenarios. However, the special farming method of rice flooding breaks this assumption. When the general method is directly applied, its characterization ability degrades in the critical early stage from transplanting to tillering due to the physical interference of high coherence on the water surface.

[0114] The coherence decomposition scheme disclosed in this embodiment constructs a binary complex coherence model of "canopy-water surface": using the dynamic estimation of canopy coverage driven by σ° time series as a bridge, and using a physical model of weighted summation of coherence coefficients in the complex domain, the coherence contributions of the water surface and the canopy are analyzed from the original observation γ, and the pure canopy structure coherence is used to replace the observation value and sent into the subsequent processing chain.

[0115] This scheme automatically satisfies boundary constraints under two extreme conditions of canopy cover: when canopy cover approaches 0 (full water surface), the decomposition result degenerates into the observed value, because the water surface is the only source of coherence at this time, and the observed value itself is a valid signal; when canopy cover approaches 1 (complete coverage), the decomposition result approaches the observed value, because the influence of the water surface has been naturally eliminated. This method does not make structural changes to the original process, but only performs pre-treatment physical purification, which can systematically eliminate the contamination of the phenological monitoring chain by the background of rice flooding, and improve the key phenological identification accuracy from transplanting to tillering stage by 13 percentage points.

[0116] Although embodiments of the present invention have been disclosed above, they are not limited to the applications listed in the specification and embodiments. It can be applied to various fields suitable for the present invention. Further modifications can be readily implemented by those skilled in the art.

Claims

1. A method for monitoring crop phenology by combining SAR backscattering and interferometric coherence, characterized in that, Includes the following steps: Step 1: Acquire time-series synthetic aperture radar (SAR) image data covering the entire growth cycle of crops, and extract the backscattering coefficient time series and interferometric coherence time series after preprocessing. Step 2: Based on the backscattering coefficient time series and the interference coherence time series, a joint sensitivity evaluation matrix is ​​constructed for each crop growth stage. The elements of the joint sensitivity evaluation matrix include the time change rate of the backscattering coefficient and the interference coherence, the coefficient of variation, and the independence measure that characterizes the degree of overlap between the two information. Step 3: Based on the joint sensitivity evaluation matrix, calculate the time domain weights and spatial domain weights. The time domain weights are determined by the product of the local information entropy of the backscattering coefficient or the time series of interference coherence and the spatial autocorrelation index. The spatial domain weights are determined by the ratio of the local variance to the covariance of the backscattering coefficient and the interference coherence. Step 4: Using time domain weights and spatial domain weights, the backscattering coefficient and interferometric coherence are weighted and fused, and a coupling coefficient is introduced during the fusion process to characterize the nonlinear interaction between the two, generating fusion features; Step 5: Using the fusion features, combined with the preset backscattering coefficients of bare soil and lush vegetation, a normalized SAR complex phenological index is constructed. By analyzing the temporal changes of this SAR complex phenological index, the crop phenological period is determined.

2. The method according to claim 1, characterized in that, In step 2, the rate of change over time is the partial derivative of the backscattering coefficient or interference coherence with respect to time during a specific crop growth stage; the coefficient of variation is the ratio of the standard deviation to the mean of the backscattering coefficient or interference coherence during a specific crop growth stage; in step 2, the independence measure is the difference between 1 and the absolute value of the Pearson correlation coefficient between the backscattering coefficient and the interference coherence.

3. The method according to claim 1, characterized in that, In step 3, the formula for calculating the local information entropy H is: H = -∑ p(x i log2 p(x) i ); where x i p(x) represents the backscattering coefficient or interference coherence within the sliding window. i ) is x i The probability of it appearing within the sliding window.

4. The method according to claim 1, characterized in that, In step 3, the formula for calculating the spatial domain weight Ws is: Ws = |Var(σ°) - Var(γ)| / ( |Cov(σ°, γ)| + ε); where Var(σ°) is the variance of the backscattering coefficient within the local spatial window, Var(γ) is the variance of the interference coherence within the local spatial window, Cov(σ°, γ) is the covariance of the backscattering coefficient and the interference coherence within the local spatial window, and ε is a preset positive constant.

5. The method according to claim 1, characterized in that, In step 4, the formula for calculating the fusion feature Ffusion is: Ffusion = Wt · σ° + Ws · γ + λ · σ° · γ; where Wt is the time domain weight, Ws is the spatial domain weight, σ° is the backscattering coefficient, γ is the interference coherence, and λ is the coupling coefficient.

6. The method according to claim 1, characterized in that, In step 5, the formula for constructing the normalized SAR complex climate index SCPI is: SCPI = α · Φstab + β · (σ° - σ°soil) / (σ°veg - σ°soil); where Φstab is the phase stability index in the fusion features, σ° is the current backscattering coefficient, σ°soil is the preset bare soil baseline backscattering coefficient, σ°veg is the preset lush vegetation baseline backscattering coefficient, and α and β are normalization coefficients, satisfying α + β = 1.

7. The method according to claim 1, characterized in that, Before determining the crop phenological period in step 5, the method further includes: using a dynamic time warping algorithm to align and match the time series curve of the constructed SAR complex phenological index with the standard phenological period time series curve of the crop established based on historical observation data, so as to adaptively correct the time offset of the phenological period. After using the dynamic time warping algorithm for alignment and matching, the algorithm also includes: calculating the time series similarity score after aligning the time series curve of the SAR complex phenological index with the time series curve of the standard phenological period; when the time series similarity score is higher than a preset threshold, the corrected phenological period determination result is output.

8. The method according to any one of claims 1 to 7, characterized in that, The crop is rice, and the synthetic aperture radar (SAR) image data comes from the Sentinel-1 satellite.

9. The method according to claim 5, characterized in that, In step 4, the value of the coupling coefficient λ is determined through the following iterative process: Set an initial value for λ and calculate the corresponding fusion feature Ffusion; Calculate the first mutual information between the fusion feature and the backscattering coefficient time series, and the second mutual information between the fusion feature and the interference coherence time series; Using the sum of the first mutual information and the second mutual information as the objective function, an optimization algorithm is used to iteratively adjust λ until the objective function value is maximized. The λ at this point is then used as the determined coupling coefficient.

10. The method according to claim 8, characterized in that, After extracting the interference coherence timing sequence in step 1 and before step 2, the following steps are also included: The interference coherence during the flooded empty field period before rice transplanting was used as the reference water surface coherence. Based on the backscattering coefficient time series, the temporal canopy coverage is calculated. Canopy coverage characterizes the proportion of water surface obstructed by the rice canopy. By utilizing the coherence of the reference water surface, the temporal canopy coverage, and the interference coherence time series, the coherence time series of the canopy structure is obtained through coherence decomposition. The canopy structure coherence timing sequence is used as an alternative interference coherence timing sequence to perform steps 2 to 5.

Images