Working system and method for performing vegetation coverage anomaly analysis based on remote sensing image

Abstract

The invention provides a technical method for performing vegetation coverage anomaly analysis based on a remote sensing image, which mainly comprises the following steps of: firstly, performing registration and preprocessing on a multi-temporal SAR (Synthetic Aperture Radar) image, and performing adaptive threshold coherence detection between adjacent time phases to quickly identify suspicious felling or severely sparse regions; secondly, dual-channel polarization InSAR height inversion is only carried out on the suspicious areas, and the height change of a forest canopy is accurately obtained so as to distinguish forest damage from other slight disturbances; and time sequence analysis is further performed on the multi-temporal image in combination with a small baseline subset InSAR technology, baseline optimization is realized in cooperation with vegetation coverage segmentation and coherence screening, and progressive or sudden vegetation coverage anomalies can be identified. And finally, according to comprehensive characteristics such as height change, coherence anomaly and phase jump, outputting the type and distribution of vegetation anomaly.

Description

Technical Field

[0001] The present invention relates to the technical field of signal processing, and in particular to a working system and method for analyzing vegetation coverage anomalies based on remote sensing images. Background Art

[0002] With the increasing demand for monitoring illegal logging, pest outbreaks, and forest health conditions, there is an urgent need for a technical means that can adapt to seasonal vegetation fluctuations, take into account multi-temporal analysis, and is easy to deploy in large-scale mountainous areas. In the research on monitoring forest vegetation cover changes using synthetic aperture radar (SAR), traditional methods usually face the following problems:

[0003] Coherence detection with a fixed threshold: If a constant threshold is used to judge the change of interference coherence, it is often impossible to distinguish the coherence fluctuations caused by seasonal leaf changes, and false alarms or missed alarms are likely to occur.

[0004] Lack of accurate height information of polarimetric InSAR: Single-intensity or single-polarization data cannot reliably distinguish different types of vegetation changes such as deforestation, pest damage, or natural attenuation; while using polarimetric interference technology can invert forest height, but the calculation amount is large, and if inversed point by point in the whole region, the efficiency is low.

[0005] Limited application of time-series InSAR in forests: Time-series InSAR technologies such as SBAS are mostly used for surface deformation monitoring. When directly applied to high-vegetation coverage areas, it is often difficult to obtain stable time-series information due to decorrelation and phase clutter. The lack of distinction of vegetation seasonal differences leads to improper baseline or coherence screening.

[0006] Realistic error factors: The terrain of mountain forests is complex, and the atmosphere, baseline error, and noise may all interfere with the interference coherence and height inversion results. The lack of targeted processing means often affects the detection accuracy. Summary of the Invention

[0007] The present invention aims to at least solve the technical problems existing in the prior art, and particularly innovatively proposes a working method for analyzing vegetation coverage anomalies based on remote sensing images, including:

[0008] S1, Coherence detection with an adaptive threshold: For the time phase of two adjacent SAR acquisitions or a specified period, calculate the interference coherence of the pixel region to obtain the observed coherence coefficient;

[0009] Set a weighted coherence threshold according to the vegetation coverage or historical statistical value. If the current coherence drops more than the weighted threshold relative to the reference coherence, it is determined as a significant change;

[0010] S2, Polarimetric InSAR Height Inversion: Perform polarimetric interference processing on two polarimetric SAR images in the same coordinate system, and calculate the complex interference coherence of each polarimetric channel;

[0011] For each local patch, use the Fourier–Legendre polynomial expansion of the random volume-ground model and the vertical scattering function for non-linear inversion to obtain the parameters of vegetation height and ground phase;

[0012] Based on the inversion results, generate a canopy height map or height change map of the target area;

[0013] S3, Temporal InSAR Analysis with Baseline Optimization: Divide multiple SAR acquisition time phases into two categories: high vegetation period and low vegetation period. Perform coherence filtering and baseline optimization on the formed small baseline interferometric pairs respectively, and eliminate the interferometric pairs with too low coherence or too large baseline;

[0014] If a phase mutation appears in the temporal analysis and is consistent with coherence or height change, it is considered that there is a significant vegetation anomaly.

[0015] In a preferred embodiment of the present invention, the coherence detection of the adaptive threshold includes:

[0016]

[0017] Among them, is the interference coherence corresponding to the polarimetric channel p observed from the SAR image;

[0018] μ p (x, y) is the pixel weighting factor, which is comprehensively determined by the signal-to-noise ratio, local incident angle, and filtering output confidence. Its calculation method is:

[0019]

[0020] Among them, κ p is the global calibration coefficient under the polarimetric channel p, which is used to adjust the overall magnitude of the weighting factor between different polarizations;

[0021] SNR p (x, y) represents the signal-to-noise ratio under the polarimetric channel p;

[0022] θ inc (x, y) is the local incident angle of the pixel, θ 0,p is the best incident angle of the pixel, χ p is used to control the attenuation intensity when deviating from the best angle;

[0023] Ξ p (x, y) is the filtering output confidence, and its value range is Ξ p(x, y) ∈ [0, 1], where 0 represents extreme suspicion of the pixel quality and 1 represents full trust in the pixel quality;

[0024] Ω is a local window that performs weighted superposition of data of multiple adjacent pixels within a local range, and expands or shrinks according to the computing power of the system;

[0025] s 1,p (x, y) and s 2,p (x, y) are the complex pixel values in the polarization channel p at the coordinate (x, y) for two different dates respectively;

[0026] represents the conjugate complex number of the complex pixel value;

[0027] represents performing phase correction within a local range, Δθ p (x, y) represents the phase correction value calculated for the pixel at (x, y) within a local range in the polarization channel p, and j represents the imaginary unit, satisfying j 2 = -1;

[0028] α p (x, y) and β p (x, y) represent the factors for weighting the two signal powers respectively;

[0029] |*| represents the complex modulus, such that will fall within the range of [0, 1];

[0030] For images of two adjacent times or a specified time period, calculate the difference in interferometric coherence:

[0031]

[0032] Δγ(x, y, t) represents the difference in reference coherence under stable conditions at the coordinate (x, y);

[0033] represents the interferometric coherence at the coordinate (x, y) during the latest SAR acquisition;

[0034] γ prev (x, y) represents the average coherence under stable conditions in the reference historical data;

[0035] Adopt an adaptive threshold τ c , and use prior knowledge of vegetation cover to weight the threshold:

[0036] τ c (t) = ηΦ(FVC(x, y, t), NDVI(x, y, t), Γ(x, y)) + δ;

[0037] where τc (t) is the coherence threshold at time t;

[0038] FVC(x, y, t) is the vegetation coverage at (x, y) at time t;

[0039] NDVI(x, y, t) is the optimized vegetation index, which is used to mitigate the influence of the boundary problem on NDVI. It performs piecewise mapping based on the size of NDVI. According to prior knowledge, parameters a1, b1, c1, d1, e1, f1, and g1 are adjusted so that their value ranges are [0, 1]:

[0040]

[0041] where v(x, y, t) is the original normalized vegetation index, and its value range is [0, 1], v th is the piecewise threshold of the NDVI value. In the region where v(x, y, t) ≤ v th , the vegetation information is weak, and fractional polynomial mapping is adopted. In the region where v(x, y, t) ≥ v th , the vegetation tends to be saturated, and an exponential decay-type nonlinear function is adopted, so that the optimized normalized vegetation index can more finely distinguish vegetation categories;

[0042] F(v(x, y, t)) is the fractional polynomial mapping, and G(v(x, y, t)) is the exponential decay-type nonlinear function. Its constraint conditions include: having continuity and derivative continuity at the segmentation point, that is, F(v th ) = G(v th ), F′(v th ) = G′(v th ), F′(v th ) is the derivative function of the fractional polynomial mapping, and G′(v th ) is the derivative function of the exponential decay-type nonlinear function;

[0043] a1, b1, and c1 are polynomial coefficients, which are used to fit the relationship between the low-value section of NDVI and its corrected value;

[0044] ln(*) is the logarithmic function, and d1 is the logarithmic decay coefficient, which is used to adjust the growth rate of the fractional part to ensure smooth output when the NDVI value is low and prevent violent fluctuations due to the input value approaching 0;

[0045] e1 is the coefficient that regulates the overall nonlinear mapping amplitude, and f1 is the factor that controls the exponential decay rate, so that when NDVI exceeds the threshold, the output corrected value does not saturate quickly but shows a gentle upward or downward trend;

[0046] g1 is a constant offset term used to ensure that the corrected value in the high NDVI region is consistent with the actual observation;

[0047] Γ(x, y) is the statistical value of historical coherence;

[0048] η represents the dependence of the threshold on the vegetation coverage, and δ represents the reference offset;

[0049] Φ(*) is a non-linear mapping function that takes into account FVC, NDVI, and historical coherence simultaneously. The higher its output, the lower the coherence level at this moment and this location under normal circumstances, so the threshold is set lower; otherwise, the threshold is increased. Its calculation method is as follows:

[0050]

[0051] where λ1, λ2, and λ3 are the adjustment factors for each term respectively;

[0052] If is lower than τ c γ(t) or Δγ(x, y, t) is less than -τ c (t), it is determined that the coherence has an abnormal decrease and further analysis is required.

[0053] In a preferred embodiment of the present invention, the polarimetric InSAR height inversion includes:

[0054] Performing polarimetric interference processing on each pair of time phases, and calculating the complex coherence coefficient of each polarization channel according to the physical model of polarimetric InSAR

[0055]

[0056] where is the interference coherence constructed according to the physical model of polarimetric InSAR;

[0057] ρ p represents the ground-volume scattering ratio coefficient under polarization channel p, which can comprehensively reflect the surface scattering intensity and the proportion of relative volume scattering;

[0058] φ g represents the phase of ground scattering, which is affected by the dielectric properties of the surface;

[0059] φ(z) represents the phase of volume scattering;

[0060] j represents the imaginary unit, satisfying j 2 = -1;

[0061] ω p(z) is the multiple scattering correction factor, used to simulate the attenuation or amplification of scattering with depth caused by multiple scattering and occlusion effects;

[0062] h v is the canopy height map;

[0063] ΔC p is the phase correction term of strong local scatterers in polarization channel p;

[0064] ΔD p represents an additional correction to the total energy in polarization channel p, used to compensate for the scattering amplitude;

[0065] μ p (z) is the vertical scattering distribution function, which varies with height z and represents the scattering power distribution from the ground surface to the top of the forest canopy;

[0066] This function is expanded using Fourier–Legendre polynomials to fit the non-uniform vertical scattering distribution, and the canopy height map h is solved by an iterative algorithm v , and its calculation method is:

[0067]

[0068] where z is the vertical height coordinate, and its value range is [0, h v ;

[0069] a k (t) is the coefficient that varies with time for the k-th term of the polynomial expansion, reflecting the change of forest canopy density or structure with time;

[0070] P k (*) is the k-th order Legendre polynomial, used to expand the vertical scattering distribution;

[0071] Q is the highest order of the polynomial expansion;

[0072] Perform non-linear least squares solution on the planned interferometric observations of two temporal phases to estimate the canopy height h v and the surface phase h g , and this process processes the images in blocks in parallel.

[0073] In a preferred embodiment of the present invention, the time-series InSAR analysis with baseline optimization includes:

[0074] Screen the numerous interferometric pairs formed by multi-temporal SAR images, only leaving the interferometric pairs with relatively high coherence and both small temporal and spatial baselines, perform SBAS solution, and solve the multi-temporal phases of each pixel in the SBAS framework to detect phase jumps or abnormal trends:

[0075] Δφ ab(x,y) = ψ ab (x,y)[φ SBAS (x,y,t b ) - φ SBAS (x,y,t a )] + σ veg (x,y,t a ,t b );

[0076] Δφ ab (x,y) is the phase difference at the coordinates (x,y) between two times;

[0077] φ SBAS (x,y,t b ) and φ SBAS (x,y,t a ) respectively represent the interferometric phase values after SBAS solution at the coordinates (x,y) at time t b and time t a ;

[0078] ψ ab (x,y) is the weight factor used to measure the credibility of the interferogram pair at this pixel;

[0079] σ veg (x,y,t a ,t b ) is the vegetation phase correction, and its calculation method is:

[0080] σ veg (x,y,t a ,t b ) = σ0Δh v (x,y,t a ,t b ) + τ(t a ,t b );

[0081] σ0 is used to reflect the sensitivity of the influence of height change on vegetation indicators in different regions;

[0082] τ(t a ,t b ) represents the coherence fluctuation caused by seasonal leaf growth or defoliation;

[0083] If Δφ ab (x,y) is too large and the coherence is abnormal at the same time, it is regarded as a strong abnormal signal.

[0084] The present invention also discloses a computer system, including:

[0085] A processor;

[0086] A memory for storing processor-executable instructions;

[0087] Wherein, when the processor is configured to execute the executable instructions, the disclosed working method for analyzing vegetation cover anomalies based on remote sensing images is implemented.

[0088] The present invention also discloses a computer-readable storage medium, including:

[0089] A memory having a computer program stored thereon;

[0090] A processor for executing the program in the memory to implement the disclosed working method for analyzing vegetation cover anomalies based on remote sensing images.

[0091] In summary, due to the adoption of the above technical solutions, the beneficial effects of the present invention are:

[0092] By using the fractional vegetation cover (FVC) or historical statistics to perform seasonal weighting on the coherence threshold, it is possible to effectively distinguish normal seasonal leaf changes from true abnormal disturbances, and reduce the high false alarm or missed detection phenomena caused by a fixed threshold.

[0093] Through prior coherence screening, a large area scene is reduced to suspected abnormal areas, and then polarimetric InSAR height inversion is performed in these areas. This not only reduces the computational burden but also can accurately obtain the forest canopy height or height change with the help of multi-polarimetric interference information to distinguish real deforestation from general disturbances.

[0094] By introducing high / low vegetation period baseline optimization and coherence screening, low-quality interferometric pairs are removed from the time series to achieve early detection of events such as progressive forest degradation (such as pest damage) and concealed logging.

[0095] If a phase mutation appears in the SBAS solution result and coincides with the coherence and tree height change, the abnormal type can be further confirmed.

[0096] For the decorrelation caused by baseline errors, residual atmospheric phase, and vegetation scattering in mountain forests, compensation or elimination is performed through short baselines, small time differences grouping, and multi-channel polarimetric information. Combining parallel or block processing modes improves the recognition reliability and timeliness under complex terrain conditions.

[0097] In summary, the present invention organically integrates adaptive coherence detection, polarimetric InSAR height inversion, and SBAS time series analysis through a multi-level and multi-stage method. It can not only timely capture sudden deforestation and other drastic changes in large forest areas but also identify progressive pest damage or vegetation degradation, and has high robustness and practical value in complex mountain environments.

[0098] Additional aspects and advantages of the present invention will be given in part in the following description, become apparent in part from the following description, or be learned by practice of the present invention. BRIEF DESCRIPTION OF THE DRAWINGS

[0099] The above and / or additional aspects and advantages of the present invention will become apparent and be readily understood from the description of the embodiments in conjunction with the following drawings, in which:

[0100] Figure 1 is a schematic diagram of the system of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0101] Embodiments of the present invention will be described in detail below. Examples of the embodiments are shown in the drawings, wherein the same or similar reference numerals denote the same or similar elements or elements having the same or similar functions throughout. The embodiments described below by referring to the drawings are exemplary only for explaining the present invention and should not be construed as limiting the present invention.

[0102] As Figure 1 shown, the present invention discloses a working system and method for analyzing vegetation cover anomalies based on remote sensing images, including:

[0103] Monitoring the vegetation cover of mountain forests is crucial for early detection of anomalies (such as deforestation or pest outbreaks) that threaten the health of the ecosystem. These anomalies may precede or exacerbate environmental disasters (such as landslides after deforestation, wildfires after tree death caused by pests). The goal of this system is to provide a near-real-time analysis method that uses a single remote sensing mode to detect abnormal changes in forest vegetation cover with high reliability. Innovative solutions are designed based on advanced SAR (Synthetic Aperture Radar) remote sensing technologies (especially polarimetric InSAR and coherence change detection).

[0104] Main objectives:

[0105] Timely detection: Detect sudden or gradual vegetation changes by analyzing each new satellite acquisition to achieve near-real-time monitoring.

[0106] Robustness for mountainous terrain: Use technologies that can handle the complex terrain and dense vegetation common in mountainous areas. Mountain forests often exhibit seasonal decorrelation and scattering effects that must be addressed to avoid false alarms.

[0107] Sensitivity to different anomalies: Utilize structural and radiometric indices to detect sudden changes (e.g., complete deforestation) and gradual changes (e.g., thinning of the tree canopy related to pests).

[0108] High-confidence alerts: Reduce noise and false alarms by combining vegetation structure modeling and adaptive thresholds that take into account the normal seasonal variations of vegetation cover.

[0109] After evaluating various possible methods (such as optical multispectral and SAR), synthetic aperture radar (SAR) was selected as the sole data source for the system. SAR is an active microwave remote sensing method that can provide all-weather, day-night imaging, which is very useful for mountainous areas that are often shrouded in clouds or shadows. More importantly, SAR data contains information about vegetation structure and changes through its backscatter intensity and phase coherence between repeated passes:

[0110] Penetration and structure: SAR signals (especially those with longer wavelengths, such as C-band or L-band) penetrate the forest canopy and interact with tree trunks and the ground. Using polarimetric SAR interferometry (PolInSAR), forest height and vertical structure can be estimated. This structural insight cannot be directly obtained from multispectral optical data.

[0111] Interferometric coherence: By comparing the phase information of two SAR acquisitions, interferometric coherence can be measured, which is sensitive to changes in scatterers over time. Forest areas usually exhibit lower coherence due to vegetation movement and growth; sudden deforestation stabilizes the scatterers (bare ground), resulting in a significant change in coherence. Therefore, SAR coherence is a powerful change indicator that can complement intensity-based detection.

[0112] Frequent revisit: With the Sentinel-1 (C-band SAR) satellite, the system can update anomaly detection almost in real time. In contrast, optical detection may be hindered by clouds, thus delaying detection.

[0113] By using SAR as a single mode, the system can operate continuously in mountain forest environments and capture gradually emerging declines in health (through changes in backscatter or height over time) and sudden disturbances (through coherence or height drops). The following technical methods utilize advanced SAR technologies (polarimetric interferometry, coherence analysis, and time series InSAR processing) to detect vegetation cover anomalies with high confidence.

[0114] 1. Coherence change detection with an adaptive threshold

[0115] Acquire two or more SAR images, and after registration and removal of the flat-earth phase preprocessing;

[0116] Perform a weighted calculation of the interferometric coherence between adjacent (or specified period) time phases, and use an adaptive threshold (WCTM) to determine whether it is abnormal;

[0117] Mark the pixels / regions with significant coherence mutations (suspected deforestation, severe sparseness, etc.).

[0118] Coherence Change Detection with Adaptive Threshold: To capture subtler or partial changes (and verify altitude changes), the system analyzes the InSAR coherence between consecutive SAR images.

[0119]

[0120] Among them, is the interferometric coherence corresponding to the polarization channel p observed from the SAR image;

[0121] μ p (x, y) is the pixel weighting factor, which is comprehensively determined by the signal-to-noise ratio, local incident angle, and filter output confidence. Its calculation method is:

[0122]

[0123] Among them, κ p is the global calibration coefficient under the polarization channel p, which is used to adjust the overall magnitude of the weighting factor between different polarizations;

[0124] SNR p (x, y) represents the signal-to-noise ratio under the polarization channel p;

[0125] θ inc (x, y) is the local incident angle of the pixel, and θ 0,p is the optimal incident angle of the pixel, and χ p is used to control the attenuation intensity when deviating from the optimal angle;

[0126] Ξ p (x, y) is the filter output confidence, and its value range is Ξ p (x, y) ∈ [0, 1], where 0 represents extreme suspicion of the pixel quality, and 1 represents full trust in the pixel quality;

[0127] Ω is a local window that weights and superimposes the data of multiple adjacent pixels within a local range, and expands or shrinks according to the computing power of the system;

[0128] s 1,p (x, y) and s 2,p (x, y) are the complex pixel values in the polarization channel p at the coordinate (x, y) for two dates respectively;

[0129] represents the conjugate complex number of the complex pixel value;

[0130] represents phase correction within a local range, and Δθ p (x, y) represents the phase correction value calculated for the pixel at (x, y) within a local range in the polarization channel p. j represents the imaginary unit, satisfying j 2=-1. The main function of this correction is to solve the problem that when averaging signals in a certain area, misaligned phases will cancel each other out, resulting in the final interference coherence being lower than the true value. Introducing phase correction can fine-tune the phase of each pixel to ensure that before superposition, the phases of all pixels are consistent, so that the modulus value of the weighted complex sum can more truly reflect the coherence level of this area.

[0131] α p (x,y) and β p (x,y) respectively represent the factors for weighting the power of two signals;

[0132] |*| represents the complex modulus, such that will fall within the range of [0,1];

[0133] The actual observed coherence coefficient is obtained through calculation, and this result will be used as the input for subsequent inversion.

[0134] Measure the scene consistency between two measurements. Healthy, undisturbed forests usually have medium to low coherence (due to vegetation movement and growth), while areas of sudden deforestation will cause a decrease in coherence between the images before and after deforestation (because the scatterers have changed), and subsequently, after the ground is exposed, the coherence in subsequent image pairs will be higher.

[0135] For adjacent two or specified time period images, calculate the interference coherence:

[0136]

[0137] Δγ(x,y,t) represents the difference in reference coherence under stable conditions at coordinates (x,y);

[0138] represents the interference coherence at coordinates (x,y) during the latest SAR acquisition;

[0139] γ prev (x,y) represents the average coherence under stable conditions in the reference historical data;

[0140] The weighted consistency threshold method (WCTM) is used to determine whether the consistency change is significant, taking into account the normal consistency changes caused by vegetation density and seasonal phenology. Essentially, instead of calculating a fixed consistency threshold, an adaptive threshold τ c is calculated. This threshold is lower in seasons or regions with dense vegetation (even if the undisturbed consistency is very low) and higher in periods with sparse vegetation. This means using prior knowledge of vegetation cover to "weight" the threshold. For example, if FVC is the vegetation cover fraction (ranging from 0 - 1), set:

[0141] τ cτ(t) = ηΦ(FVC(x,y,t), NDVI(x,y,t), Γ(x,y)) + δ;

[0142] where τ c (t) is the coherence threshold at time t;

[0143] FVC(x,y,t) is the fractional vegetation cover at (x,y) at time t;

[0144] NDVI(x,y,t) is the optimized vegetation index, which is used to mitigate the impact of the NDVI being vulnerable to boundary problems. It is segmented and mapped based on the magnitude of the NDVI. According to prior knowledge, the parameters a1, b1, c1, d1, e1, f1, and g1 are adjusted so that their value ranges are [0,1]:

[0145]

[0146] where v(x,y,t) is the original normalized difference vegetation index, with a value range of [0,1], and v th is the segmentation threshold of the NDVI value. In the region where v(x,y,t) ≤ v th , the vegetation information is weak, and a fractional polynomial mapping is used. In the region where v(x,y,t) ≥ v th , the vegetation tends to be saturated, and an exponential decay-type non-linear function is used, so that the optimized normalized difference vegetation index can more finely distinguish vegetation categories;

[0147] F(v(x,y,t)) is the fractional polynomial mapping, and G(v(x,y,t)) is the exponential decay-type non-linear function. Its constraint conditions include: having continuity and derivative continuity at the segmentation point, that is, F(v th ) = G(v th ), F′(v th ) = G′(v th ), F′(v th ) is the derivative function of the fractional polynomial mapping, and G′(v th ) is the derivative function of the exponential decay-type non-linear function;

[0148] a1, b1, and c1 are polynomial coefficients used to fit the relationship between the low-value section of the NDVI and its corrected value;

[0149] ln(*) is the logarithmic function, and d1 is the logarithmic decay coefficient used to adjust the growth rate of the fractional part to ensure smooth output when the NDVI value is low and prevent violent fluctuations due to the input value approaching 0;

[0150] e1 is the coefficient that regulates the amplitude of the overall nonlinear mapping, and f1 is the factor that controls the exponential decay rate, so that when the NDVI exceeds the threshold, the output correction value will not be saturated quickly, but will show a gentle upward or downward trend;

[0151] g1 is a constant offset term used to ensure that the correction value in the high NDVI area is consistent with the actual observation;

[0152] Γ(x,y) is the statistical value of historical coherence;

[0153] η represents the dependence of the threshold on vegetation coverage, and δ represents the benchmark offset;

[0154] Φ(*) is a nonlinear mapping function that takes into account FVC, NDVI, and historical coherence. The higher its output, the lower the coherence level at this location at this moment without abnormalities, and thus the threshold is set lower; otherwise, the threshold is increased. The calculation method is:

[0155]

[0156] Among them, λ1, λ2 and λ3 are the adjustment factors of each item;

[0157] In practice, FVC(x,y,t) can be estimated from the SAR itself or by a seasonal assumption (e.g., summer = high FVC, winter = low). Any region where the coherence falls below an adaptive threshold (compared to a reference coherence under stable conditions) is flagged as changing. Similarly, the appearance of persistently high coherence in a region that once had low coherence (indicating that it has lost volume scattering) can signal an anomaly. This adaptive coherence filtering helps filter out false positives, e.g., a forest that typically has low coherence in the summer will not be incorrectly flagged as uncorrelated due to decorrelation, because the threshold there will be lower.

[0158] like Lower than τ c (t) or Δγ(x,y,t) is less than -τ c (t), it is judged that the coherence has dropped abnormally and further analysis is required.

[0159] 2. Inversion of vegetation height using dual-channel PolInSAR

[0160] Vegetation height inversion via dual-channel PolInSAR: A method for estimating forest canopy height from only two SAR images using polarimetric interferometry is integrated. This method uses a combination of the random volume on the ground (RVoG) model and the extended Fourier-Legendre polynomial representation of vertical structure to retrieve tree heights. By applying it to duplicate channel data (e.g., two acquisitions in the same area), the vegetation height h is derived. vmap.

[0161] 2.1 Polarimetric InSAR and Forest Height Estimation

[0162] The polarimetric InSAR model characterizes the relationship between interferometric coherence and forest parameters. In fact, the solution is a non-linear inversion process that requires an iterative algorithm to minimize the difference between the observed γ p and the model, and each local patch is processed independently, and this method can be accelerated in parallel;

[0163] Mathematically, the PolInSAR observation is the complex interferometric coherence γ of different polarization channels p p , and solve the inversion model:

[0164] For the marked highly suspicious areas, dual-channel (or multi-channel) polarimetric InSAR inversion is performed to accurately obtain the tree height or tree height change;

[0165] Perform polarimetric interferometry on each pair of time phases (or key time phases) (if the data has polarization information), and calculate the complex coherence coefficients of each polarization channel

[0166]

[0167] where is the interferometric coherence constructed based on the physical model of polarimetric InSAR;

[0168] ρ p represents the ground-volume scattering ratio coefficient under polarization channel p, which can comprehensively reflect the surface scattering intensity and the proportion of relative volume scattering;

[0169] φ g represents the phase of ground scattering, which is affected by the surface dielectric properties;

[0170] φ(z) represents the phase of volume scattering;

[0171] j represents the imaginary unit, satisfying j 2 = -1;

[0172] ω p (z) is the multi-scattering correction factor, which is used to simulate the attenuation or amplification of scattering with depth caused by multiple scattering and occlusion effects;

[0173] h v is the canopy height map;

[0174] ΔC p is the phase correction term of strong local scatterers under polarization channel p;

[0175] ΔD pRepresents an additional correction to the total energy in polarization channel p for compensating the scattering amplitude;

[0176] μ p (z) is the vertical scattering distribution function, which varies with height z and represents the scattering power distribution from the ground surface to the top of the forest canopy;

[0177] This function is expanded using Fourier–Legendre polynomials to fit the non-uniform vertical scattering distribution, and the canopy height map h is solved through an iterative algorithm v , and its calculation method is:

[0178]

[0179] where z is the vertical height coordinate, and its value range is [0, h v ;

[0180] a k (t) is the time-varying coefficient of the k-th term of the polynomial expansion, reflecting the change of forest canopy density or structure over time;

[0181] P k (*) is the k-th order Legendre polynomial used to expand the vertical scattering distribution;

[0182] Q is the highest order of the polynomial expansion;

[0183] Nonlinear least squares solution is performed on the planned interferometric observations of two temporal phases to estimate h v and the surface phase h g , and this process can process the image in blocks (N×M pixels) in parallel;

[0184] The canopy height map h is obtained v ;

[0185] Assume that a small area of N×M pixels shares the same h v and the surface phase h g , and then a nonlinear iterative solver is applied to estimate h v . Even using a single interferometric baseline, the robustness can be improved.

[0186] 2.2 Relative height difference

[0187] The system regularly obtains the forest height estimate or at least the relative height change of each area:

[0188] Δh v (x,y,t1,t2) = h v (x,y,t2) - h v (x,y,t1);

[0189] where Δh v(x, y, t1, t2) represents the change value of the canopy height at the coordinate (x, y) between time t2 and t1. If the relative height is negative, it indicates canopy loss;

[0190] As time progresses, the estimated h v A significant decrease is a strong indicator of deforestation.

[0191] Height inversion is only performed in suspicious areas to reduce the computational load and can more accurately determine the type of anomaly.

[0192] 3. Time series with baseline optimization

[0193] To systematically accumulate evidence of change and handle long-term monitoring, the Small Baseline Subset (SBAS) InSAR framework is adopted, which is enhanced by considering baseline optimization for vegetation phases. All available SAR images in the area are used to form multiple interferometric pairs ("interferogram network"). Then, the time series is divided into segments with high and low vegetation cover (e.g., long leaf period and deciduous period) and processed slightly differently. According to the patent, the method is as follows:

[0194] Each SAR acquisition is classified as occurring in a high-vegetation period or a low-vegetation period based on the time of year (or an external indicator such as NDVI if available). For example, the FVC is higher in the leafy summer and lower in winter or the dry season.

[0195] All possible interferograms (within specific time baseline limits) are formed and their coherence is calculated. Then, interferograms with extremely low coherence within each group are discarded. This step removes pairs that may have excessive decorrelation (due to too long a time interval or large geometric differences during the growing season) - ensuring high-quality remaining interferometric data ("optimized baseline network"). The result is two sets of interferograms: one mainly from the low-vegetation period and one from the high-vegetation period, with only relatively coherent pairs retained in each set.

[0196] Using the retained "high-quality" interferometric pairs, the multi-temporal phase φ SBAS (x, y, t) of each pixel is solved under the SBAS framework to detect phase jumps or abnormal trends:

[0197] Δφ ab (x, y) = ψ ab (x, y)[φ SBAS (x, y, t b ) - φ SBAS (x, y, t a )] + σ veg (x, y, t a , t b );

[0198] Δφ ab (x, y) is the phase difference at the coordinate (x, y) between two times;

[0199] φ SBAS (x, y, t b ) and φ SBAS (x, y, t a ) respectively represent the interferometric phase values after SBAS solution at the coordinate (x, y) at time t b and time t a ;

[0200] ψ ab (x, y) is the weight factor used to measure the credibility of this interferometric pair on this pixel;

[0201] σ veg (x, y, t a , t b ) is the vegetation phase correction, and its calculation method is:

[0202] σ veg (x, y, t a , t b ) = σ0Δh v (x, y, t a , t b ) + τ(t a , t b );

[0203] σ0 is used to reflect the sensitivity of the influence of height change on vegetation indicators in different regions;

[0204] τ(t a , t b ) represents the coherence fluctuation caused by seasonal leaf growth or defoliation;

[0205] If Δφ ab (x, y) is too large and the coherence is abnormal at the same time, it is regarded as a strong abnormal signal.

[0206] These optimized interferogram sets are input into the SBAS InSAR solver. The SBAS algorithm will invert the phase information to estimate the change of each pixel over time. Traditionally, SBAS is used for ground deformation; here, the phase or coherence change is interpreted as a proxy for vegetation change. For example, significant unwrapped phase changes across the suspected event date in the interferogram may indicate ground movement or substantial vegetation loss (which changes the effective scattering phase center height). However, by looking at the coherence timeline and the heights derived from PolInSAR, ground movement can be distinguished from vegetation change. Vegetation phase separation mainly comes from the combination of using short baselines and segmentation by vegetation state - this minimizes the mixing of true ground displacement phase and random vegetation phase. Introducing optical data stratified by vegetation season and using coherence as a selection criterion can greatly reduce the impact of vegetation cover on the radar signal. In summary, the time series analysis component provides a method to detect progressive trends: for example, a slow, continuous decrease in canopy height or coherence over several months (pest infestation) or a sudden jump between two dates (deforestation). It also allows time averaging to reduce noise.

[0207] 4. Anomaly Decision

[0208] Combine height change, coherence anomaly, and SBAS phase jump:

[0209] If a significant sharp drop in height or a substantial change in coherence is found, it may be large-scale deforestation;

[0210] If there is a progressive change in height / scattering, it may point to pest infestation or other chronic stress;

[0211] Then consider geographical proximity (such as whether continuous patches are formed) to exclude random noise.

[0212] Output the vegetation anomaly distribution map, marking the warning levels (high-confidence suspected deforestation or possible pest infestation, etc.).

[0213] 5. Implementation Cases

[0214] 5.1 Implementation Background

[0215] Study area: A mountain forest (altitude between 800 and 2000 meters), with a coverage area of about 200 km 2 , and the main vegetation is coniferous forest mixed with some broad-leaved forests.

[0216] SAR data:

[0217] Use Sentinel-1 dual-polarization (VV+VH) data, with a resolution of about 10×10m;

[0218] Two scenes of SAR images were acquired in March 202X (spring, early leaf germination) and May (early summer, lush canopy), and used as adjacent temporal phases or specified time period pairs;

[0219] There are also multiple scenes of data (including 5 scenes in July, September of the current year and January of the following year, etc.) for time series analysis.

[0220] DEM and auxiliary data:

[0221] A set of DEM with 30m accuracy is used for terrain correction;

[0222] Historical vegetation data: including NDVI sequence, statistics of forestry department (vegetation coverage FVC is about 0.7 - 0.9).

[0223] 5.2 Coherence Detection with Adaptive Threshold

[0224] The two scenes of Sentinel-1 images in March and May are processed with the core gpt or self-developed InSAR software:

[0225] Registered to the same coordinate system;

[0226] Use DEM to remove the flat-earth phase and most of the terrain phase;

[0227] Calculate the interferometric coherence at both VV and VH polarizations

[0228] The threshold becomes smaller in the high FVC area (>0.8) to avoid false alarms due to generally low coherence caused by summer leaves; the threshold is slightly larger in the low FVC area.

[0229] Compared with the previous temporal phase (March) or historical mean: If the decrease amplitude Δγ(x, y, t) of the current coherence relative to the historical reference coherence is greater than the set threshold, mark "possibly abnormal".

[0230] Through statistics, it is found that the normal coherence is relatively low in most areas during this period, but only some patches have a decrease exceeding the threshold, suspected of being illegally logged or showing signs of large-area sparseness.

[0231] 5.3 Polarimetric InSAR Height Inversion

[0232] Within the range of the suspicious patches, the SLC data (VV and VH polarizations) of the images in March and May are extracted respectively, and high-precision registration is performed again;

[0233] Only perform PolInSAR processing on these patches (about 100×100 pixels per patch) to reduce the overall computational amount.

[0234] Calculate the complex interferometric coherence in the VV and VH channels, and multi-look within the range of window size (5×5) to improve the signal-to-noise ratio.

[0235] Parallel iterative processing is performed on each patch on a GPU or multi-core CPU, and finally a "canopy height" or "height change" map per pixel is obtained.

[0236] If the inversion result of a certain patch shows that the height suddenly drops from 15 meters to less than 5 meters, it is determined that the trees in this area have been severely deforested.

[0237] After comparing the height results within the suspicious patches with the surrounding reference areas, a high-resolution canopy height change map is formed, indicating the positions of pixels with significant drops.

[0238] 5.4 Temporal InSAR Analysis with Baseline Optimization

[0239] The research team subsequently obtained more Sentinel-1 images at different times, including those in July, September of the same year, and January of the following year, for a total of 5 scenes.

[0240] They were divided into two categories: "high vegetation period (May, July, September)" and "low vegetation period (March, January of the following year)".

[0241] In the high vegetation period, only interferometric pairs with short temporal differences or small spatial baselines are retained to ensure that the coherence is not too low; similar operations are also performed in the low vegetation period.

[0242] Pairs with coherence < 0.15 or baseline > 200 meters are excluded, resulting in two subnetworks, each containing 6 - 8 high-quality interferometric pairs.

[0243] The SBAS algorithm is applied to the high / low vegetation subnetworks respectively to obtain the multi-temporal phase sequences of each pixel.

[0244] If there is a sudden phase jump (combined with S1 / S2 observations), it indicates significant scatterer changes; if there is a slow and continuous phase shift, it may also be a progressive process (such as pest damage, etc.).

[0245] Result determination:

[0246] By cross-verifying the SBAS phase mutations with the coherence anomalies in S1 and the sudden tree height drops in S2, the truly significant deforestation or vegetation degradation areas are identified.

[0247] If a certain area shows stability (no mutations) in the time series but has low coherence in a certain polarimetric InSAR, it is only a seasonal effect or noise, and false alarms are excluded.

[0248] 5.6 Implementation Effect

[0249] First, "coherence threshold detection" is completed for the entire area using S1, significantly reducing the computational amount, and only 2% of the suspicious patches are subjected to S2 polarimetric InSAR height inversion.

[0250] Finally confirmed that there are 5 places (about 1.3 km in total) 2 ) were illegally cut down between March and May, and the tree height decreased by >50%.

[0251] Through the analysis of subsequent images in July and September by S3, it was confirmed that the coherence dropped suddenly again in 2 of them, accompanied by a continuous decrease in the tree crown. Suspecting the spread of pests, a warning was issued in a timely manner.

[0252] The SBAS baseline optimization also ensures that there are still enough high-quality interferometric pairs during the multi-vegetation period in the mountainous area.

[0253] 5.6.1 Accuracy and reliability

[0254] The comparison results with on-site sampling surveys show that the detection accuracy for the deforested area is over 92%;

[0255] Due to the adoption of adaptive thresholds and temporal baseline optimization, the false alarm rate caused by seasonal noise and topographic disturbances has decreased significantly.

[0256] Although the embodiments of the present invention have been shown and described, those of ordinary skill in the art can understand that various changes, modifications, substitutions, and variations can be made to these embodiments without departing from the principles and spirit of the present invention. The scope of the present invention is defined by the claims and their equivalents.

Claims

1. A working method for analyzing abnormal vegetation coverage based on remote sensing images, characterized in that, include: S1, adaptive threshold coherence detection: for two adjacent SAR acquisition phases or specified time periods, the interference coherence of the pixel area is calculated to obtain the observed coherence coefficient; A weighted coherence threshold is set according to vegetation coverage or historical statistical values, and if the current coherence decreases by more than the weighted threshold relative to the reference coherence, it is determined to be a significant change; S2, polarimetric InSAR height inversion: performing polarimetric interferometry processing on the two polarimetric SAR images in the same coordinate system, and calculating the complex interferometric coherence of each polarization channel; For each local patch, the parameters of vegetation height and ground phase are obtained by nonlinear inversion using the random volume-ground model and the Fourier-Legendre polynomial expansion of the vertical scattering function. Based on the inversion results, a canopy height map or height change map of the target area is generated; S3, time series InSAR analysis with baseline optimization: multiple SAR acquisition phases are divided into two categories: high vegetation period and low vegetation period. The small baseline interference pairs formed are subjected to coherence filtering and baseline optimization respectively, and the interference pairs with too low coherence or too large baseline are eliminated; If a phase jump occurs in the time series analysis and is consistent with a change in coherence or height, a significant vegetation anomaly is considered to exist.

2. The working method for analyzing vegetation cover anomalies based on remote sensing images according to claim 1, characterized in that, The coherence detection of the adaptive threshold comprises: Among them, is the interference coherence corresponding to the polarization channel p observed from the SAR image; μ p (x, y) is the pixel weighting factor, which is comprehensively determined by the signal-to-noise ratio, the local incident angle, and the confidence of the filter output. Its calculation method is as follows: where κ p is the global calibration coefficient under polarization channel p, which is used to adjust the overall magnitude of the weighting factors between different polarizations; SNR p (x, y) represents the signal-to-noise ratio in polarization channel p; θ inc (x, y) is the local incident angle of the pixel, and θ 0,p is the optimal incident angle of the pixel, and χ p is used to control the attenuation intensity when deviating from the optimal angle; Ξ p (x, y) is the confidence of the filtered output, and its value range is Ξ p (x, y) ∈ [0, 1]. 0 represents extreme suspicion of the pixel quality, and 1 represents full trust in the pixel quality; Ω is a local window that performs weighted superposition of data from multiple adjacent pixels in a local range, expanding or shrinking according to the computing power of the system; s 1,p (x, y) and s 2,p (x, y) are the complex pixel values of the two dates in the polarization channel p at the coordinates (x, y), respectively; denotes the conjugate complex number of the complex pixel value; Indicates performing phase correction within a local range, Δθ p (x, y) represents the phase correction value calculated for the pixel at (x, y) within the local range in polarization channel p. j represents the imaginary unit, satisfying j 2 = -1; α p (x, y) and β p (x, y) respectively represent the factors for weighting the two signal powers; |*| represents the complex modulus such that will fall within the range of [0, 1]; For two consecutive images or images of a specified period of time, calculate the difference in interference coherence: Δγ(x,y,t) represents the difference in reference coherence when the coordinates (x,y) are under stable conditions; Indicates the interferometric coherence at the coordinate (x, y) during the latest SAR acquisition; γ prev (x, y) represents the average coherence under stable conditions in the reference historical data; Adopt an adaptive threshold τ c , and use the prior knowledge of vegetation cover to weight the threshold: τ c (t) = ηΦ(FVC(x, y, t), NDVI(x, y, t), Γ(x, y)) + δ; Among them, τ c (t) is the coherence threshold at time t; FVC(x,y,t) is the vegetation cover at (x,y) at time t; NDVI(x,y,t) is an optimized vegetation index used to reduce the vulnerability of NDVI to boundary problems. It is segmented based on the NDVI size. Based on prior knowledge, the parameters a1, b1, c1, d1, e1, f1, and g1 are adjusted to have a value range of [0,1]: Among them, v(x, y, t) is the original normalized vegetation index, and its value range is [0, 1]. v th is the segmentation threshold of the NDVI value. In the area where v(x, y, t) ≤ v th , the vegetation information is weak, and the fractional polynomial mapping is adopted. In the area where v(x, y, t) ≥ v th , the vegetation tends to be saturated, and the exponential decay type nonlinear function is adopted, so that the optimized normalized vegetation index can more finely distinguish vegetation categories; F(v(x, y, t)) is a fractional polynomial mapping, and G(v(x, y, t)) is an exponential decay-type nonlinear function. Its constraint conditions include: having continuity and derivative continuity at the segmentation points, that is, F(v th ) = G(v th ), F′(v th ) = G′(v th ), F′(v th ) is the derivative function of the fractional polynomial mapping, and G′(v th ) is the derivative function of the exponential decay-type nonlinear function; a1, b1, and c1 are polynomial coefficients used to fit the relationship between the NDVI in the low-value segment and its corrected value; ln(*) is a logarithmic function, and d1 is a logarithmic attenuation coefficient, which is used to adjust the growth rate of the fractional part to ensure smooth output when the NDVI value is low and does not fluctuate violently due to input values approaching 0; e1 is the coefficient that regulates the amplitude of the overall nonlinear mapping, and f1 is the factor that controls the exponential decay rate, so that when the NDVI exceeds the threshold, the output correction value will not be saturated quickly, but will show a gentle upward or downward trend; g1 is a constant offset term used to ensure that the correction value in the high NDVI area is consistent with the actual observation; Γ(x,y) is the statistical value of historical coherence; η represents the dependence of the threshold on vegetation coverage, and δ represents the benchmark offset; Φ(*) is a nonlinear mapping function that takes into account FVC, NDVI, and historical coherence. The higher its output, the lower the coherence level at this location at this moment without abnormality, and thus the threshold is set lower; otherwise, the threshold is raised. The calculation method is: Among them, λ1, λ2, and λ3 are the adjustment factors for each item; If is lower than τ c (t) or Δγ(x, y, t) is less than -τ c (t), it is determined that the coherence has abnormally decreased and further analysis is required.

3. The working method for analyzing abnormal vegetation coverage based on remote sensing images according to claim 1, characterized in that, The polarimetric InSAR height inversion includes: Perform polarimetric interferometric processing on each pair of temporal phases, and calculate the complex coherence coefficients of each polarimetric channel according to the physical model of polarimetric InSAR. Among them, is the interference coherence constructed based on the physical model of polarimetric InSAR; ρ p represents the ground-volume scattering ratio coefficient under polarization channel p, which can comprehensively reflect the surface scattering intensity and the proportion of relative volume scattering; φ g represents the phase of ground scattering, which is affected by the dielectric properties of the earth's surface; φ(z) represents the phase of volume scattering; j represents the imaginary unit, satisfying j 2 = -1; ω p (z) is the multiple scattering correction factor, which is used to simulate the attenuation or amplification of scattering with depth caused by multiple scattering and occlusion effects; h v is the canopy height map; ΔC p is the phase correction term of the strong local scatterer under the polarization channel p; ΔD p represents an additional correction to the total energy in polarization channel p for compensating the scattering amplitude; μ p (z) is the vertical scattering distribution function, which varies with the height z and represents the scattering power distribution from the ground surface to the top of the forest canopy; The function is expanded with Fourier-Legendre polynomials to fit the non-uniform vertical scattering distribution, and the canopy height map h is solved by an iterative algorithm v , and its calculation method is as follows: where z is the vertical height coordinate with a range of [0, h v . a k (t) is the time-varying coefficient of the polynomial expansion term k, reflecting the change of forest canopy density or structure over time; P k (*) is the k-th order Legendre polynomial, which is used to expand the vertical scattering distribution; Q is the highest order of polynomial expansion; Perform a non - linear least - squares solution for the planned interferometric observations at two time phases to estimate the canopy height h v and the surface phase h g , and in this process, the images are processed in blocks in parallel.

4. A working method for analyzing abnormal vegetation coverage based on remote sensing images according to claim 1, characterized in that, The time-series InSAR analysis with baseline optimization includes: Screening a large number of interferometric pairs formed by multi-temporal SAR images, only leaving the interferometric pairs with relatively high coherence, small temporal and spatial baselines, performing SBAS solution, and solving the multi-temporal phases of each pixel in the SBAS framework to detect phase jumps or abnormal trends: Δφ ab (x,y) = ψ ab (x,y)[φ SBAS (x,y,t b ) - φ SBAS (x,y,t a )] + σ veg (x,y,t a ,t b ); Δφ ab (x, y) is the phase difference at the coordinates (x, y) between two times; φ SBAS (x, y, t b ) and φ SBAS (x, y, t a ) respectively represent the interferometric phase values after SBAS solution at the coordinate (x, y) at time t b and time t a ; ψ ab (x, y) is a weighting factor used to measure the credibility of the interference pair on the pixel; σ veg (x, y, t a , t b ) is the vegetation phase correction, and its calculation method is as follows: σ veg (x, y, t a , t b ) = σ0Δh v (x, y, t a , t b ) + τ(t a , t b ); σ0 is used to reflect the sensitivity of the influence of height changes on vegetation indices in different regions; τ(t a ,t b ) represents the coherence fluctuations brought about by seasonal leaf growth or leaf fall; If Δφ ab (x, y) is too large and the coherence is abnormal at the same time, it is regarded as a strong abnormal signal.

5. A computer system, characterized in that, including: a processor; a memory for storing processor-executable instructions; Among them, the processor is configured to implement a working method for analyzing vegetation cover anomalies based on remote sensing images according to any one of claims 1 to 4 when executing the executable instructions.

6. A computer-readable storage medium, characterized in that, including: a memory having a computer program stored thereon; a processor for executing the program in the memory to implement a working method for analyzing vegetation cover anomalies based on remote sensing images according to any one of claims 1 to 4.

Images