Mountain adaptive distributed scatterometer interferometry method, device and medium

Abstract

This disclosure relates to the field of remote sensing and mapping technology, and provides an adaptive distributed scatterer interferometry method, apparatus, and medium for mountainous areas. The method includes: acquiring radar image data and digital elevation model (DEM) data of the target monitoring mountainous area; calculating the radar vegetation index of each pixel in the dual-polarization synthetic aperture radar (DAP) image sequence and classifying it into multiple vegetation coverage levels to obtain a vegetation coverage grading atlas; identifying statistically homogeneous pixels in each image scene under the spatial constraints of vegetation coverage levels according to the grading atlas to determine a candidate set of distributed scatterers; further, using the grading atlas, performing adaptive coherence threshold filtering on each candidate point to obtain multiple distributed scatterers; and determining the long-term surface deformation results of the mountainous area based on the multiple distributed scatterers, the DAP image sequence, and the DEM data. This embodiment effectively improves the accuracy and reliability of distributed scatterer extraction and is suitable for deformation monitoring in complex mountainous areas.

Description

Technical Field

[0001] This disclosure relates to the field of remote sensing mapping technology, and more specifically, to a method, apparatus, and medium for adaptive distributed scatterer interferometry in mountainous areas. Background Technology

[0002] Interferometric Synthetic Aperture Radar (InSAR), as a non-contact surface deformation monitoring technology, has been widely applied in fields such as geological disaster early warning and engineering safety assessment due to its advantages of all-weather, all-time, wide-area, and high precision. Temporal InSAR technology further enhances deformation monitoring capabilities, with Persistent Scatterer InSAR (PS-InSAR) and Small Baseline Set InSAR (SBAS-InSAR) being the two most representative methods. PS-InSAR models deformation based on point targets with strong reflection and stable phase, while SBAS-InSAR obtains continuous planar deformation fields by constructing a small time-space baseline interferometric atlas.

[0003] The aforementioned methods have achieved good results in typical scenarios such as cities and seismic zones, but they have inherent limitations in complex mountainous applications: PS-InSAR has sparse monitoring point distribution in low-coherence areas such as vegetation and bare soil, resulting in insufficient spatial coverage; SBAS-InSAR is sensitive to phase jumps and noise interference, easily introducing estimation bias. The vegetation cover characteristics of complex mountainous areas further amplify these problems: dense vegetation in mountainous areas leads to severe spatiotemporal decoherence of radar signals, and related methods often use a global single threshold to select measurement points, which easily results in sparse monitoring points, uneven distribution, or even large-scale "coherence holes"; the selection of monitoring points faces a contradiction between quantity and quality, lowering the threshold can increase the number but introduce low-quality pixels, while maintaining a high threshold makes it difficult to ensure the monitoring point density in densely vegetated areas; the diverse types of surface cover in mountainous areas, with different vegetation cover having significantly different effects on coherence, and related methods have not fully considered this difference, leading to poor quality of measurement point selection. Summary of the Invention

[0004] This disclosure provides at least one method, apparatus, and medium for adaptive distributed scatterer interferometry in mountainous areas. By introducing radar vegetation index to classify vegetation coverage in dual-polarization SAR images of mountainous areas, and on this basis, implementing statistical homogeneous pixel identification and adaptive coherence threshold screening under spatial constraints, the accuracy and reliability of distributed scatterer extraction are effectively improved, making it suitable for deformation monitoring in complex mountainous areas.

[0005] This disclosure provides an adaptive distributed scatterer interferometry method for mountainous areas, including:

[0006] The radar image data and digital elevation model data of the target monitoring mountain area are acquired; wherein, the radar image data includes a single-polarization synthetic aperture radar image sequence and a dual-polarization synthetic aperture radar image sequence, the dual-polarization synthetic aperture radar image sequence includes multiple dual-polarization SAR images, and each dual-polarization SAR image includes multiple pixels.

[0007] For the dual-polarization synthetic aperture radar image sequence, the radar vegetation index of each pixel in each dual-polarization SAR image is calculated, and based on the calculation results, the pixels in each dual-polarization SAR image are divided into multiple vegetation coverage levels to obtain a vegetation coverage classification atlas corresponding to the dual-polarization synthetic aperture radar image sequence.

[0008] Based on the vegetation coverage grading atlas, statistical homogeneous pixel identification is performed on each pixel in each dual-polarization SAR image under the spatial constraint of vegetation coverage level, and a set of distributed scatterer candidate points corresponding to the dual-polarization synthetic aperture radar image sequence is determined based on the identification results; and, based on the vegetation coverage grading atlas, an adaptive coherence threshold is used to filter each distributed scatterer candidate point in the set of distributed scatterer candidate points to obtain multiple distributed scatterers.

[0009] Based on the multiple distributed scatterers, the single-polarization synthetic aperture radar image sequence, and digital elevation model data, the long-term surface deformation results of the target monitoring mountain area are determined.

[0010] This disclosure provides an adaptive distributed scatterer interferometry device for mountainous areas, comprising:

[0011] The data acquisition module is used to acquire radar image data and digital elevation model data of the mountainous area for target monitoring; wherein, the radar image data includes a single-polarization synthetic aperture radar image sequence and a dual-polarization synthetic aperture radar image sequence, the dual-polarization synthetic aperture radar image sequence includes multiple dual-polarization SAR images, and each dual-polarization SAR image includes multiple pixels.

[0012] The pixel classification module is used to calculate the radar vegetation index of each pixel in each dual-polarization synthetic aperture radar image sequence, and to classify the pixels in each dual-polarization synthetic aperture radar image sequence into multiple vegetation coverage levels based on the calculation results, so as to obtain a vegetation coverage classification map set corresponding to the dual-polarization synthetic aperture radar image sequence.

[0013] The scatterer identification module is used to identify statistically homogeneous pixels in each pixel of each dual-polarization SAR image under the spatial constraint of vegetation coverage level, based on the vegetation coverage grading atlas, and to determine a set of distributed scatterer candidate points corresponding to the dual-polarization synthetic aperture radar image sequence based on the identification results; and to perform adaptive coherence threshold filtering on each distributed scatterer candidate point in the set of distributed scatterer candidate points based on the vegetation coverage grading atlas to obtain multiple distributed scatterers.

[0014] The result determination module is used to determine the long-term surface deformation results of the target monitoring mountain area based on the multiple distributed scatterers, the single-polarization synthetic aperture radar image sequence, and digital elevation model data.

[0015] This disclosure provides a computer device, including a processor, a memory, and a bus. The memory stores machine-readable instructions executable by the processor. When the computer device is running, the processor communicates with the memory via the bus. When the machine-readable instructions are executed by the processor, the mountain adaptive distributed scatterer interferometry method as described in any of the above possible embodiments is executed.

[0016] This disclosure provides a computer-readable storage medium storing a computer program that, when executed by a processor, implements the mountain adaptive distributed scatterer interferometry method as described in any of the above possible embodiments.

[0017] The adaptive distributed scatterer interferometry method, apparatus, and medium for mountainous areas provided in this disclosure specifically utilize radar vegetation index-based hierarchical processing to finely characterize the spatial distribution differences of vegetation on complex underlying surfaces in mountainous areas, providing effective prior constraint information for subsequent pixel identification. By performing statistical homogeneous pixel identification within the vegetation coverage level space, the interference of heterogeneous materials on pixel selection is reduced, making the candidate point set of distributed scatterers more accurately reflect the scattering characteristics under different vegetation coverage conditions. Furthermore, by combining the vegetation coverage hierarchical results with adaptive coherence threshold screening, the screening criteria can be dynamically adjusted according to the coherence attenuation characteristics of different regions, thereby retaining more effective measurement points and avoiding the applicability limitations of a single threshold in complex mountainous environments.

[0018] In this way, by introducing radar vegetation index to classify vegetation coverage in dual-polarization SAR images of mountainous areas, and on this basis implementing statistical homogeneous pixel identification and adaptive coherence threshold screening under spatial constraints, more dense and reliable distributed scatterers can be extracted, thereby effectively improving the spatial continuity and solution accuracy of long-term surface deformation monitoring in mountainous areas.

[0019] To make the above-mentioned objects, features and advantages of this disclosure more apparent and understandable, preferred embodiments are described below in detail with reference to the accompanying drawings. Attached Figure Description

[0020] To more clearly illustrate the technical solutions of the embodiments of this disclosure, the accompanying drawings referenced in the embodiments will be briefly described below. These drawings are incorporated in and constitute a part of this specification. They illustrate embodiments conforming to this disclosure and, together with the specification, serve to explain the technical solutions of this disclosure. It should be understood that the following drawings only show some embodiments of this disclosure and should not be considered as limiting the scope. Those skilled in the art can obtain other related drawings based on these drawings without creative effort.

[0021] Figure 1 A flowchart of an adaptive distributed scatterer interferometry method for mountainous areas provided by an embodiment of this disclosure is shown;

[0022] Figure 2 A flowchart of a vegetation cover grading method provided in an embodiment of this disclosure is shown;

[0023] Figure 3 A flowchart of a method for identifying statistical homogeneous pixels and determining candidate point sets of distributed scatterers provided in an embodiment of this disclosure is shown;

[0024] Figure 4 A flowchart of a statistical homogeneity test method provided by an embodiment of this disclosure is shown;

[0025] Figure 5 A flowchart of an adaptive coherence threshold filtering method provided by an embodiment of this disclosure is shown;

[0026] Figure 6 A flowchart of a method for determining long-term surface deformation results provided by an embodiment of this disclosure is shown;

[0027] Figure 7 A schematic diagram of an experimental monitoring grading map of vegetation cover in mountainous areas provided by an embodiment of this disclosure is shown;

[0028] Figure 8 A schematic diagram showing a comparison of homogeneous pixel recognition results using different methods provided in an embodiment of this disclosure is illustrated.

[0029] Figure 9 This illustration shows a schematic diagram comparing the spatial distribution of homogeneous samples using different methods, as provided in an embodiment of this disclosure.

[0030] Figure 10A schematic diagram showing a comparison of the selection results of a fixed threshold and an adaptive threshold distributed scatterer provided in an embodiment of this disclosure is illustrated.

[0031] Figure 11 A schematic diagram of a target monitoring surface deformation rate map of a mountainous area provided by an embodiment of this disclosure is shown;

[0032] Figure 12 A schematic diagram of optical image interpretation of a localized deformation region provided by an embodiment of this disclosure is shown;

[0033] Figure 13 A schematic diagram illustrating the interpretation of local features in a key deformation region provided by an embodiment of this disclosure is shown;

[0034] Figure 14 A schematic diagram of the structure of an adaptive distributed scatterer interferometry device for mountainous areas provided in an embodiment of this disclosure is shown.

[0035] Figure 15 A schematic diagram of the structure of a computer device provided in an embodiment of this disclosure is shown. Detailed Implementation

[0036] To make the objectives, technical solutions, and advantages of the embodiments of this disclosure clearer, the technical solutions of the embodiments of this disclosure will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of this disclosure, and not all of them. The components of the embodiments of this disclosure described and shown in the accompanying drawings can generally be arranged and designed in various different configurations. Therefore, the following detailed description of the embodiments of this disclosure provided in the accompanying drawings is not intended to limit the scope of the claimed disclosure, but merely represents selected embodiments of this disclosure. All other embodiments obtained by those skilled in the art based on the embodiments of this disclosure without inventive effort are within the scope of protection of this disclosure.

[0037] It should be noted that similar labels and letters in the following figures indicate similar items. Therefore, once an item is defined in one figure, it does not need to be further defined and explained in subsequent figures.

[0038] In this document, the term "and / or" merely describes a relationship, indicating that three relationships can exist. For example, A and / or B can represent three cases: A alone, A and B simultaneously, and B alone. Furthermore, the term "at least one" in this document means any combination of at least two of any one or more elements. For example, including at least one of A, B, and C can mean including any one or more elements selected from the set consisting of A, B, and C.

[0039] Synthetic Aperture Radar Interferometry (InSAR), as a non-contact, all-day, all-weather surface deformation monitoring technology, can acquire surface displacement information with millimeter-level accuracy over a large area through multi-temporal SAR image interferometry. It has been widely used in fields such as geological disaster early warning and engineering safety assessment.

[0040] Research has shown that temporal InSAR technology, with its advantages of high precision, wide coverage, and high spatiotemporal resolution, has become a core means of monitoring surface deformation in mountainous areas. Persistent scatterer InSAR (PS-InSAR) and small baseline set InSAR (SBAS-InSAR) are the two most representative technologies. PS-InSAR models deformation based on point targets with strong reflection and stable phase, while SBAS-InSAR obtains continuous planar deformation fields by constructing interferometric atlases of small spatiotemporal baselines and using least squares estimation. Although both have achieved good results in typical scenarios such as urban areas and seismic zones, they both have significant limitations. The vegetation cover characteristics of complex mountainous areas further amplify these limitations, with the core constraints concentrated in three aspects: prominent vegetation incoherence problems, unresolved contradictions in the selection of monitoring points, and insufficient adaptability of surface cover.

[0041] Specifically, PS-InSAR has a sparse point distribution in low-coherence areas such as vegetation and bare soil, resulting in limited spatial coverage. SBAS-InSAR is sensitive to phase jumps and noise interference, easily introducing estimation bias. Traditional time-series InSAR uses a global single threshold to select measurement points, which easily leads to sparse and uneven distribution of monitoring points, or even large-scale "coherence holes," making it impossible to fully capture deformation details. Lowering the point selection threshold can increase the number of monitoring points, but it may introduce a large number of low-quality pixels, affecting the solution accuracy. If a high threshold is maintained, it is difficult to ensure the monitoring point density in densely vegetated areas, resulting in insufficient spatial sampling and an inability to balance the quantity and quality of monitoring points. Mountainous areas have diverse surface cover types, with significant differences in scattering characteristics between different cover types, and vegetation cover exhibits obvious hierarchical differences, each with different impacts on coherence. The aforementioned methods do not fully consider these differences, resulting in poor quality of measurement point selection and difficulty in adapting to complex surface cover scenarios.

[0042] To overcome the aforementioned limitations, some technologies have proposed InSAR phase optimization methods based on distributed scatterers (DS). Early methods relied on amplitude information statistical verification to identify statistically homogeneous pixels, but amplitude and phase have no direct mathematical correlation, and similar amplitude does not necessarily indicate phase stability, easily leading to misselection in areas with severe deformation such as landslide boundaries. Some improved methods introduced land cover type optimization for point selection, but did not specifically consider the impact of vegetation cover grade differences on coherence, resulting in insufficient adaptability. In recent years, while pixel recognition methods based on phase behavior are closer to the physical meaning of phase, they still face key challenges: First, phase data is sensitive to noise, easily fluctuating drastically in low coherence areas or under strong noise interference, leading to unstable point selection and inaccurate covariance estimation; second, the thresholds for indicators such as correlation coefficient and phase consistency, on which point selection depends, lack clear physical or statistical interpretations, are highly empirical and uncertain, and are difficult to uniformly set in different scenarios, resulting in insufficient repeatability and interpretability of point selection results.

[0043] Based on the above research, this disclosure provides an adaptive distributed scatterer interferometry method, device, and medium for mountainous areas, comprising: firstly acquiring radar image data and digital elevation model data of the target monitoring mountainous area, wherein the radar image data includes a single-polarization synthetic aperture radar (SAR) image sequence and a dual-polarization SAR image sequence, the dual-polarization SAR image sequence containing multiple dual-polarization SAR images and each image consisting of multiple pixels; then calculating the radar vegetation index of each pixel in each image sequence, and dividing the pixels in each image into multiple [unclear] based on the calculation results. A vegetation cover grading atlas is generated for each image sequence, based on the vegetation cover level. Then, statistical homogeneous pixel identification is performed on each pixel in each image scene under the spatial constraint of vegetation cover level, according to the grading atlas. Based on the identification results, a candidate set of distributed scatterers corresponding to the image sequence is determined. Adaptive coherence thresholding is then applied to each candidate point according to the grading atlas to obtain multiple distributed scatterers. Finally, based on the selected distributed scatterers, single-polarization synthetic aperture radar image sequences, and digital elevation model data, the long-term surface deformation results of the target monitoring mountain area are determined.

[0044] In this embodiment, by introducing radar vegetation index to classify vegetation coverage of dual-polarization SAR images in mountainous areas, and on this basis implementing statistical homogeneous pixel identification and adaptive coherence threshold screening under spatial constraints, more dense and reliable distributed scatterers can be extracted, thereby effectively improving the spatial continuity and solution accuracy of long-term surface deformation monitoring in mountainous areas.

[0045] To facilitate understanding of this embodiment, the executing entity of the adaptive distributed scatterer interferometry method for mountainous areas provided in this disclosure will first be described in detail. The executing entity of the adaptive distributed scatterer interferometry method for mountainous areas provided in this disclosure is a computer device. This computer device can be a terminal device or a server. The terminal device can also be a mobile device, a user terminal, a terminal, a handheld device, a computing device, an in-vehicle device, a wearable device, etc. The server can be an independent physical server, a server cluster or distributed system composed of multiple physical servers, or a cloud server providing basic cloud computing services such as cloud services, cloud databases, cloud computing, cloud storage, big data, and artificial intelligence platforms. Optionally, this method can also be applied to an implementation environment composed of computer devices and servers.

[0046] The adaptive distributed scatterer interferometry method for mountainous areas provided in this application will be described in detail below with reference to the accompanying drawings. See also: Figure 1 The diagram shows a flowchart of an adaptive distributed scatterer interferometry method for mountainous areas provided in this embodiment of the present disclosure. The method includes the following steps S101 to S104:

[0047] S101 acquires radar imagery data and digital elevation model data of the mountainous area for target monitoring.

[0048] Understandably, "target monitoring mountainous area" refers to complex terrain areas requiring surface deformation monitoring. These could be densely vegetated volcanic areas, mining subsidence areas, or landslide-prone areas, such as the Tianchi volcanic area of ​​Changshan Mountain, the Panzhihua mining area, or the landslide zone along the Jinjiang River. Radar image data refers to the dataset composed of multi-temporal synthetic aperture radar (SAR) images covering the target monitoring mountainous area. This can include single-polarization SAR image sequences and dual-polarization SAR image sequences. Specifically, a dual-polarization SAR image sequence represents a dataset composed of multiple images acquired at different times covering the same target monitoring mountainous area. It contains information on both vertical and horizontal polarization, and features all-weather, all-day imaging and sensitivity to ground object scattering characteristics. It can be acquired through spaceborne or airborne SAR systems and focused to obtain single-view multiple images. Single-polarization synthetic aperture radar (SAR) image sequences are single-polarization image sequences used for subsequent interferometric processing. They are radar data containing only a single polarization mode (such as vertical or horizontal polarization). These can be obtained by extracting polarization channels from dual-polarization data or acquired individually by other SAR systems. Digital elevation model (DEM) data refers to a digital model describing the surface elevation information of a target monitoring mountainous area. It stores ground elevation data in a grid format, with each grid cell recording the altitude of its corresponding location. This data can be acquired through methods such as space shuttle radar topographic mapping missions, high-resolution optical stereo image pairs, or lidar measurements, for example, using the TanDEM-X or AW3D30 digital elevation models.

[0049] Here, the dual-polarization synthetic aperture radar (SAR) image sequence includes multiple dual-polarization SAR images (Synthetic Aperture Radar), which are multiple images acquired at certain time intervals (e.g., 12 or 24 days) within the same monitoring period. To accurately reflect the impact of vegetation growth in different seasons on radar backscattering characteristics, the time span of the image sequence can be set to cover the complete vegetation growth cycle, thus providing a data foundation that characterizes the annual changes in vegetation for subsequent radar vegetation index calculations. Simultaneously, to ensure the effectiveness and reliability of time-series statistical analysis, the number of images can also meet certain scale requirements, typically no less than 20 images, to obtain a sufficient sample size for statistical homogeneous pixel identification and coherence estimation, reducing the impact of random errors. Each dual-polarization SAR image includes multiple pixels. A pixel is the smallest unit constituting a radar image, and each pixel corresponds to the echo information of a ground resolution unit, including the amplitude and phase values ​​of that unit.

[0050] Similarly, like the dual-polarization synthetic aperture radar image sequence mentioned above, the single-polarization synthetic aperture radar image sequence also includes multiple SAR images. Each SAR image is also composed of multiple pixels, and each pixel corresponds to the echo information of the same resolution unit on the ground and includes amplitude and phase values.

[0051] In some possible embodiments, the above-mentioned dual-polarization synthetic aperture radar image sequence can be acquired by a C-band synthetic aperture radar system using a dual-polarization imaging mode. For example, based on the interferometric wide-swath mode of the Sentinel-1 satellite of a certain regional air traffic control bureau, repeated orbit observations are performed on the mountainous area for target monitoring. After single-view complex image focusing processing, dual-polarization SAR data containing vertical and horizontal polarization information is obtained, so as to capture the scattering characteristics of different surface cover types in the mountainous area.

[0052] In some possible embodiments, to obtain a spatially consistent image sequence for subsequent temporal analysis, after acquiring initial image data via radar, registration processing can be performed separately on all dual-polarization synthetic aperture radar (SAP) images and single-polarization SAP images. For single-polarization SAP images, one image can be selected as the master image, and the remaining images as slave images. The slave images are then registered to the radar coordinate system of the master image through geometric transformation and resampling, ensuring that the pixels of all images are spatially aligned. For dual-polarization SAP image sequences, the registration process needs to be performed independently for the vertical and horizontal polarization channels to ensure that the images of both polarization modes achieve sub-pixel level registration accuracy. The registered dual-polarization image sequence is used for subsequent radar vegetation index calculation. Since registration ensures that the ground target correspondence of the same pixel is consistent across different time phases and different polarization channels, the radar vegetation index of each pixel in each time phase can be accurately calculated. At the same time, the registered single-polarization synthetic aperture radar image sequence is the basis for subsequent generation of differential interferograms, providing spatially unified input data for homogeneous pixel identification and deformation inversion.

[0053] In some possible embodiments, initial single-view complex image data acquired by the same single-polarization synthetic aperture radar (SAP) system can be processed separately for polarization channel extraction and pre-interferometric processing to obtain the aforementioned single-polarization SAP image sequence and dual-polarization SAP image sequence. Specifically, this may include: extracting images of both vertical and horizontal polarization channels from the initial single-view complex image data to form a dual-polarization SAP image sequence for subsequent radar vegetation index calculation and vegetation cover classification; simultaneously, selecting an image of one polarization channel (such as the vertical polarization channel) from the initial single-view complex image data to form a single-polarization SAP image sequence for interferometric processing. Both processing paths share the same set of raw data, ensuring consistency between the two image sequences in terms of time and space, while avoiding resource waste caused by repeated data acquisition. The specific data acquisition method and processing path can be adjusted according to the characteristics of the actual SAP system used and the monitoring task requirements, and are not specifically limited here.

[0054] S102, for the dual-polarization synthetic aperture radar image sequence, calculate the radar vegetation index of each pixel in each dual-polarization SAR image, and based on the calculation results, divide the pixels in each dual-polarization SAR image into multiple vegetation coverage levels to obtain a vegetation coverage classification map set corresponding to the dual-polarization synthetic aperture radar image sequence.

[0055] Here, the radar vegetation index refers to an index calculated based on dual-polarization synthetic aperture radar data to quantitatively represent the degree of surface vegetation cover. It uses the ratio between the vertical polarization and horizontal polarization backscattering coefficients to reflect the intensity of vegetation scattering. The vegetation cover level is represented by the result of classifying the density of surface vegetation according to the radar vegetation index value range, which is used to distinguish areas with different vegetation cover conditions. It can be divided according to the preset radar vegetation index threshold range.

[0056] Specifically, for each pixel in each dual-polarization SAR image in a dual-polarization synthetic aperture radar image sequence, its vertical and horizontal polarization backscattering coefficients can be extracted and substituted into the radar vegetation index calculation formula to calculate its corresponding radar vegetation index. Then, according to the preset threshold range in which the radar vegetation index value falls, the pixel is divided into the corresponding vegetation coverage level to distinguish areas with different vegetation density, thereby obtaining a vegetation coverage classification map set corresponding to the image sequence. This map set includes the vegetation coverage classification map corresponding to each dual-polarization SAR image.

[0057] The formula for calculating the radar vegetation index can be expressed as:

[0058] ;

[0059] In the formula, Represented as the radar vegetation index, it is a dimensionless numerical value used to quantitatively describe the degree of vegetation cover on the ground. The value range is usually between 0 and 1, with a larger value indicating a higher degree of vegetation cover. It is expressed as the vertical polarization backscattering coefficient, which is the dimensionless value or decibel value obtained after radiometric calibration and normalization of the echo intensity of the ground target when the radar wave is transmitted and received in a vertical polarization manner. It is expressed as the horizontal polarization backscattering coefficient, which is the echo intensity of a ground target when a radar wave is transmitted in a vertical polarization mode and received in a horizontal polarization mode. It has also undergone radiometric calibration and normalization.

[0060] For example, since different vegetation cover levels result in different scattering characteristics of radar signals, and this difference affects the accuracy of subsequent homogeneous pixel identification, in order to enable spatial constraints based on vegetation cover in subsequent steps, when classifying vegetation cover, reference is made to... Figure 2 As shown, the steps S201~S204 may be included:

[0061] S201, for each dual-polarization SAR image, extract the vertical polarization backscattering coefficient and the horizontal polarization backscattering coefficient of each pixel in the dual-polarization SAR image.

[0062] Here, the vertical polarization backscattering coefficient refers to the intensity of radar wave reflection from a ground target when the radar wave is transmitted and received in a vertical polarization manner, usually expressed in decibels; the horizontal polarization backscattering coefficient is expressed as the intensity of radar wave reflection from a ground target when the radar wave is transmitted and received in a horizontal polarization manner, also expressed in decibels.

[0063] S202, for each pixel, calculate the radar vegetation index of the pixel based on the vertical polarization backscattering coefficient and the horizontal polarization backscattering coefficient corresponding to the pixel.

[0064] Understandably, based on the vertical and horizontal polarization backscattering coefficients, the radar vegetation index value for each pixel can be calculated using the aforementioned radar vegetation index calculation formula. This index can effectively distinguish land surface types with different vegetation cover levels. Here, before calculation, to ensure the correct physical meaning of the calculation results, the decibel values ​​of the two polarization coefficients can be converted to linear values.

[0065] S203, based on multiple preset radar vegetation index threshold ranges, classify pixels whose radar vegetation index values ​​fall within the same threshold range into the same vegetation coverage level.

[0066] Furthermore, after obtaining the radar vegetation index for each pixel, the vegetation coverage level of each pixel can be classified according to preset threshold intervals. These preset radar vegetation index threshold intervals serve as boundary values ​​for mapping continuously distributed radar vegetation index values ​​to discrete vegetation coverage levels, and can be set according to different regional vegetation types and monitoring needs. By classifying pixels whose radar vegetation index values ​​fall within the same threshold interval into the same vegetation coverage level, the land surface can be classified according to vegetation density, thereby determining the vegetation coverage type of each pixel.

[0067] For example, in this disclosure, vegetation coverage is divided into six levels based on radar vegetation index: bare land, low coverage, low-medium coverage, medium coverage, medium-high coverage, and high coverage. The corresponding radar vegetation index threshold ranges are set as follows: 0 ≤ radar vegetation index < 0.45 is classified as bare land, 0.45 ≤ radar vegetation index < 0.55 is classified as low coverage, 0.55 ≤ radar vegetation index < 0.65 is classified as low-medium coverage, 0.65 ≤ radar vegetation index < 0.75 is classified as medium coverage, 0.75 ≤ radar vegetation index < 0.85 is classified as medium-high coverage, and 0.85 ≤ radar vegetation index ≤ 1 is classified as high coverage.

[0068] In some other embodiments, the preset multiple radar vegetation index threshold ranges can also be adjusted according to the vegetation type distribution characteristics of the specific monitoring area. For example, different threshold division standards can be used for arid or high-altitude areas, which are not specifically limited here.

[0069] S204, the vegetation coverage level of each pixel in each dual-polarization SAR image is marked in the radar coordinate system to generate a vegetation coverage classification map corresponding to each dual-polarization SAR image; and, the vegetation coverage classification maps corresponding to all dual-polarization SAR images are combined to obtain the vegetation coverage classification map set.

[0070] Specifically, a radar coordinate system refers to a spatial reference frame established using the range and azimuth axes of a radar image as coordinate axes, used to describe the positional relationship of image pixels in imaging geometry. The radar coordinate system referred to here is the same spatial reference after registration of all dual-polarization SAR images, ensuring a one-to-one spatial correspondence of pixels in all images. Therefore, to obtain vegetation cover grading information with a unified spatial reference, after obtaining the vegetation cover level of each pixel, the vegetation cover level can be labeled under the registered radar coordinate system, thereby generating a vegetation cover grading map corresponding to each image scene.

[0071] It is understandable that by combining the vegetation coverage grading maps corresponding to all dual-polarization SAR images according to the image time sequence, a vegetation coverage grading map set covering the entire monitoring period and with a consistent spatial reference can be obtained. This map set includes the vegetation coverage grading map corresponding to each dual-polarization SAR image, which is the vegetation level labeling result of each image in the radar coordinate system, and records the vegetation coverage level information of each pixel at different time phases.

[0072] S103, based on the vegetation coverage grading atlas, statistical homogeneous pixel identification is performed on each pixel in each dual-polarization SAR image under the spatial constraint of vegetation coverage level, and a set of distributed scatterer candidate points corresponding to the dual-polarization synthetic aperture radar image sequence is determined based on the identification results; and, based on the vegetation coverage grading atlas, adaptive coherence threshold screening is performed on each distributed scatterer candidate point in the set of distributed scatterer candidate points to obtain multiple distributed scatterers.

[0073] Here, statistical homogeneous pixel identification is a method for screening pixels with similar scattering characteristics based on temporal amplitude information. By comparing the backscattering coefficient distribution of pixels over time, it can identify neighboring pixels that are statistically homogeneous with non-water body pixels, thus determining a stable sample set that can be used for subsequent phase estimation. Statistical homogeneity means that the backscattering coefficient amplitude sequences of two pixels are statistically insignificantly different across all time phases; that is, their temporal amplitude distributions come from the same continuous distribution, indicating that the two pixels have similar scattering characteristics and stability features. The vegetation cover level spatial constraint means that in the homogeneous pixel identification process, only neighboring pixels belonging to the same vegetation cover level as non-water body pixels are considered for statistical testing, thereby avoiding mismatches across vegetation types.

[0074] Specifically, based on the vegetation coverage level information of each image in the vegetation coverage grading atlas, statistical homogeneous pixel identification under spatial constraints is performed on each pixel in each dual-polarization SAR image. This allows the determination of the homogeneous pixel set corresponding to each pixel. Pixels with a sufficient number of homogeneous pixels are then identified as candidate points for distributed scatterers. This yields a set of candidate points for distributed scatterers corresponding to the dual-polarization synthetic aperture radar image sequence. This set of candidate points includes all pixels identified as candidate points at different time phases and their spatial location information.

[0075] For example, pixel scattering characteristics vary significantly across regions with different vegetation cover levels. Directly identifying homogeneous pixels can easily lead to misjudgments. Therefore, to achieve more accurate homogeneous pixel identification based on vegetation cover grading results, a method can be referenced... Figure 3 As shown, the following steps S301~S305 may be included when performing statistical homogeneous pixel recognition:

[0076] S301, construct a water body mask according to a preset backscattering coefficient threshold, and use the water body mask to remove water body pixels in each scene of dual-polarization SAR image.

[0077] Understandably, water pixels appear as low backscattering coefficient regions in radar images and exhibit poor temporal stability, making them unsuitable as candidate points for distributed scatterers. Therefore, they can be removed before processing. A water mask is a binary image used to identify water pixels. By setting a preset backscattering coefficient threshold, pixels with backscattering coefficients below that threshold can be marked as water, thus obtaining the water mask. Furthermore, by applying the water mask to each dual-polarization SAR image scene, water pixels can be removed, resulting in image data containing only non-water pixels for subsequent homogeneous pixel identification.

[0078] The preset backscattering coefficient threshold can be set according to different regions and different radar bands. For example, in C-band SAR images, pixels with a vertical polarization backscattering coefficient below -17 dB can usually be identified as water bodies; or in high-altitude and cold regions, the backscattering coefficient may increase due to water freezing, so the threshold needs to be adjusted according to the actual situation. Specifically, it can be calibrated and set according to the backscattering characteristics of typical ground features in the monitoring area, and no specific limitation is made here.

[0079] S302, for each non-water body pixel, a neighborhood window of a preset size is set with the non-water body pixel as the center, and neighborhood pixels belonging to the same vegetation coverage level as the non-water body pixel are selected within the neighborhood window.

[0080] Specifically, after removing water body pixels from each image, the remaining non-water body pixels are the candidate objects to be processed. For each non-water body pixel, in order to find pixels with similar scattering characteristics in the spatial neighborhood, a preset-sized neighborhood window can be set around the non-water body pixel. All neighboring pixels are traversed within this window, and based on the vegetation cover level information recorded in the vegetation cover grading atlas, neighboring pixels belonging to the same vegetation cover level as the central non-water body pixel are selected to narrow down the candidate range for subsequent statistical tests. Here, neighboring pixels belonging to the same vegetation cover level as the central non-water body pixel mean that they have similar land cover types, which can be used for subsequent statistical homogeneity tests.

[0081] The preset size of the neighborhood window can be set according to the image resolution and the spatial distribution characteristics of ground features. For example, for a Sentinel-1 image with a resolution of 20 meters, the window size can be set to 9×9 or 11×11, without any specific limitation.

[0082] S303, for each non-water body pixel, perform a statistical homogeneity test on the selected neighboring pixels corresponding to the non-water body pixel, and determine the neighboring pixels that pass the test and are spatially connected to the non-water body pixel as homogeneous pixels of the non-water body pixel.

[0083] Specifically, after each non-water body pixel has completed the screening of neighboring pixels with the same vegetation coverage level, a statistical homogeneity test can be performed on the screened neighboring pixels. That is, a statistical method can be used to determine whether the temporal amplitude distribution of the neighboring pixels and the central non-water body pixels comes from the same continuous distribution.

[0084] Here, refer to Figure 4 As shown, the statistical homogeneity test among pixels may include the following steps S401~S404:

[0085] S401, for each selected neighboring pixel corresponding to the non-water body pixel, extract the backscattering coefficients of the neighboring pixel and the non-water body pixel in all phases of the dual-polarization synthetic aperture radar image sequence, and construct the neighborhood temporal amplitude sequence of the neighboring pixel and the center temporal amplitude sequence of the non-water body pixel respectively.

[0086] Understandably, dual-polarization synthetic aperture radar (DAP) image sequences contain image data from multiple temporal phases. Each pixel has a corresponding backscattering coefficient record at different temporal phases. When performing statistical homogeneity tests, the amplitude information from all temporal phases can be used for judgment. Here, to construct a temporal amplitude sequence for statistical testing, the backscattering coefficients of each selected neighboring pixel at all temporal phases can be extracted. Simultaneously, the backscattering coefficients of the central non-water body pixel at all temporal phases can be extracted. By converting the decibel values ​​to linear values ​​and arranging them in temporal order, the neighborhood temporal amplitude sequence of the neighboring pixel and the central temporal amplitude sequence of the central non-water body pixel are constructed, respectively. The neighborhood temporal amplitude sequence of the neighboring pixel is represented as a vector of the amplitude values ​​of the neighboring pixel at all temporal phases arranged in chronological order, which can be used to characterize the temporal scattering characteristics of the pixel. The central temporal amplitude sequence of the non-water body pixel is a vector of the amplitude values ​​of the non-water body pixel at all temporal phases arranged in chronological order, which can be used as a benchmark for statistical comparison.

[0087] S402, perform a two-sample position test on the neighborhood temporal amplitude sequence of the neighboring pixels and the center temporal amplitude sequence of the non-water body pixels to determine whether the two sequences come from the same continuous distribution.

[0088] Specifically, the two-sample location test is a nonparametric statistical test that can be used to determine whether two independent samples come from the same population distribution without making assumptions about the distribution pattern. Here, by comparing the neighborhood temporal amplitude sequence with the center temporal amplitude sequence using the two-sample location test, it is possible to determine whether there is a significant difference in the distribution location of the two sequences, thereby determining whether the scattering characteristics of neighborhood pixels and non-water body pixels are statistically homogeneous.

[0089] S403, if the test statistic is less than the preset critical value, then the neighboring pixels are determined to have passed the statistical homogeneity test; otherwise, they are determined to have failed.

[0090] Here, if the test statistic calculated by the two-sample location test is less than the preset critical value, it indicates that there is no significant difference in the distribution of the two time-series amplitude sequences. In this case, it can be determined that the neighboring pixels and the non-water body pixels are statistically homogeneous. Otherwise, it indicates that there is a significant difference in the distribution of the two sequences. In this case, the scattering characteristics of the neighboring pixels and the non-water body pixels are different. Therefore, in this case, the neighboring pixels do not pass the statistical homogeneity test.

[0091] S404, among the neighboring pixels that have passed the test, the neighboring pixels that are directly adjacent to the non-water body pixel, and the neighboring pixels that are indirectly connected to the non-water body pixel through other neighboring pixels that have passed the test, are determined as homogeneous pixels that are spatially connected to the non-water body pixel.

[0092] Specifically, spatial connectivity refers to the formation of continuous regions between pixels through direct or indirect adjacency in spatial location. By determining whether there is a direct adjacency or indirect connection through other tested neighboring pixels between the tested neighboring pixels and the central non-water body pixel, it is possible to further filter out pixels from the tested neighboring pixels that form spatially connected regions with the non-water body pixel, thereby determining the set of homogeneous pixels that form spatially connected regions with the non-water body pixel.

[0093] Here, region growing can be performed starting from non-water body pixels using four-neighbor or eight-neighbor connectivity rules. The neighboring pixels that pass the test and are directly adjacent to the non-water body pixels, as well as other neighboring pixels that pass the test and are indirectly connected through these neighboring pixels, are collectively identified as homogeneous pixels that are spatially connected to the non-water body pixels.

[0094] In some other embodiments, statistical homogeneity testing of the selected neighborhood pixels can also be performed using other methods, such as the Kolmogorov-Smirnov test or the generalized likelihood ratio test. Taking the Kolmogorov-Smirnov test as an example, its specific steps may include: calculating the empirical cumulative distribution function of two time-series amplitude sequences, finding the maximum vertical difference between the two cumulative distribution functions as the test statistic, determining a critical value based on the sample size and significance level, and comparing the test statistic with the critical value to determine whether the two sequences originate from the same distribution. The specific test method used can be selected based on data characteristics and computational efficiency requirements, and is not specifically limited here.

[0095] Furthermore, after completing the statistical homogeneity test and spatial connectivity screening, neighboring pixels that pass the test and are spatially connected to the central non-water body pixel can be identified as homogeneous pixels of the non-water body pixel. Homogeneous pixels have similar temporal scattering characteristics and spatial continuous distribution features as non-water body pixels.

[0096] S304, count the number of homogeneous pixels corresponding to each non-water body pixel, and determine the non-water body pixels with a number of homogeneous pixels greater than the preset homogeneous pixel threshold as candidate points of distributed scatterers.

[0097] Here, after identifying homogeneous pixels for each non-water body pixel, to select stable points with sufficient sample support from a large number of non-water body pixels for subsequent processing, the number of homogeneous pixels corresponding to each non-water body pixel can be counted. Non-water body pixels with a number of homogeneous pixels greater than a preset homogeneous pixel threshold are identified as candidate distributed scatterers. These candidate points have a sufficient number of statistically homogeneous samples for subsequent phase optimization and coherence estimation. The preset homogeneous pixel threshold is used to select stable pixels with sufficient sample support and can be set according to image resolution and monitoring requirements, for example, it can be set to 15 or 20.

[0098] S305, the candidate points of distributed scatterers determined in all dual-polarization SAR images are combined according to the image time sequence to obtain the set of candidate points of distributed scatterers.

[0099] Furthermore, by organizing the candidate distributed scatterers identified in each image according to the order of image acquisition time, a candidate point set covering the entire monitoring period can be constructed, thereby obtaining a candidate distributed scatterer set corresponding to the dual-polarization synthetic aperture radar image sequence. Here, combining images according to time sequence means arranging the candidate point information from different time phases according to the time dimension. The purpose is to provide a data foundation for subsequently calculating the temporal coherence of each candidate point, so as to realize adaptive threshold screening based on time series analysis.

[0100] In this embodiment, by introducing a spatial constraint of vegetation cover level during the statistical homogeneous pixel identification process, only neighboring pixels belonging to the same vegetation cover level as the central pixel are selected for statistical testing within the neighborhood window. This effectively avoids the problem of pixel misjudgment caused by differences in scattering characteristics across vegetation type regions. Furthermore, by employing a two-sample location test to statistically compare the temporal amplitude sequences and combining it with spatial connectivity filtering, the set of homogeneous pixels with similar temporal scattering characteristics to the central pixel and forming a continuous spatial region can be identified more accurately, significantly improving the accuracy and reliability of homogeneous pixel identification. Based on this, by statistically analyzing the number of homogeneous pixels for each non-water body pixel and setting a reasonable threshold for homogeneous pixels, candidate points for distributed scatterers with sufficient statistical sample support can be selected, thereby improving the overall accuracy and stability of distributed scatterer extraction.

[0101] Understandably, after obtaining the candidate set of distributed scatterers corresponding to the dual-polarization synthetic aperture radar image sequence, due to differences in land cover type, vegetation density, and scattering stability among different candidate points, some candidate points may have low coherence and poor phase quality, which may not meet the quality requirements of subsequent deformation inversion for monitoring points. Therefore, adaptive coherence threshold screening can be performed on each candidate point in the candidate set of distributed scatterers according to the vegetation cover classification atlas to further screen out reliable distributed scatterers from the candidate points, so as to obtain multiple distributed scatterers for subsequent phase unwrapping and deformation inversion processing.

[0102] Among them, adaptive coherence threshold screening is a method that sets coherence screening thresholds according to different vegetation coverage levels. It can determine a reasonable screening threshold for each vegetation coverage level, thereby achieving a dynamic balance between the number and quality of monitoring points under different vegetation conditions.

[0103] In some possible embodiments, considering that coherence calculation needs to be performed based on differential interferograms, a differential interferogram sequence for calculating coherence can be generated before adaptive coherence threshold screening. The specific process may include the following steps (1) to (2):

[0104] (1) Based on the single-polarization synthetic aperture radar image sequence and the digital elevation model data, generate an initial differential interferogram sequence;

[0105] (2) Perform the first adaptive filtering on the initial differential interferogram sequence to obtain the differential interferogram sequence for candidate point screening, which is used to calculate the coherence coefficient of each candidate point of the distributed scatterer.

[0106] Specifically, the initial differential interferogram sequence refers to the data set consisting of interferograms generated from different image pairs in a single-polarization synthetic aperture radar (SAP) image sequence after interferometry processing, with terrain phase removed. This data is used for subsequent coherence calculations and deformation information extraction. When generating the initial differential interferogram sequence, image pairs that meet the preset spatiotemporal baseline thresholds are first selected from the SAP image sequence to form interferometric pairs. These interferometric pairs need to meet certain requirements on both the temporal and spatial baselines to ensure interferogram quality. Simultaneously, each interferometric pair undergoes interferometry processing to obtain an interferogram. This involves generating an interferometric fringe pattern containing phase difference information through conjugate multiplication, and then simulating the terrain phase using digital elevation model (DEM) data and removing it from the interferogram to obtain an initial differential interferogram containing only the deformation phase.

[0107] Here, using preset spatiotemporal baseline thresholds to filter interferometric pairs ensures high coherence in the generated interferograms, effectively reducing the impact of spatiotemporal decoherence on subsequent processing. The preset spatiotemporal baseline thresholds limit the range of image pairs participating in interferometric processing in terms of time interval and spatial vertical baseline. These thresholds can be set according to radar system parameters and monitoring area characteristics; for example, the temporal baseline threshold can be set to 36 or 48 days, and the spatial baseline threshold can be set to 150 meters or 200 meters. Interferometric processing refers to the process of generating an interferogram by conjugate multiplication of two registered single-polarization synthetic aperture radar images. This process extracts phase difference information between the two images. The interferogram is represented as a complex image containing amplitude and phase information, where the phase information records the combined effects of surface deformation, topographic relief, and atmospheric delay. Removing the topographic phase from the interferogram using digital elevation model data eliminates the contribution of topographic relief to the interferometric phase, ensuring that the remaining phase primarily reflects surface deformation information and better serves subsequent deformation inversion.

[0108] Furthermore, after obtaining the initial differential interferogram sequence, a first adaptive filtering process can be performed on it. Different filtering intensities are applied based on the local coherence characteristics of each pixel to suppress noise while preserving as much detail as possible in the interference fringes, resulting in a differential interferogram sequence for candidate point selection. This sequence is used to calculate the coherence coefficients of each distributed scatterer candidate point, enabling candidate point quality assessment and selection based on coherence indices.

[0109] For example, the first adaptive filtering process described above can use methods such as the Goldstein filtering method, the improved Goldstein filtering method, or the adaptive window filtering method, without specific limitations. The Goldstein filtering method divides the interferogram into blocks in the frequency domain and adaptively adjusts the filtering parameters according to the coherence of each block, thus suppressing noise while preserving the edge features of the interference fringes. The improved Goldstein filtering method further optimizes the determination of filtering parameters, employing more accurate coherence estimation and more flexible filtering window settings, improving the adaptability of the filtering effect. The adaptive window filtering method can dynamically adjust the size and weight of the filtering window in the spatial domain according to the local coherence level, using a small window in high-coherence regions to preserve details and a large window in low-coherence regions to enhance noise suppression.

[0110] Specifically, after obtaining the differential interferogram sequence and vegetation cover grading atlas for candidate point selection, adaptive coherence threshold screening of distributed scatterer candidate points can be performed based on these data, referring to... Figure 5 As shown, the steps S501 to S505 may be included:

[0111] S501, for each candidate distributed scatterer, based on the homogeneous pixel set of the candidate distributed scatterer, calculate the coherence coefficient of the candidate distributed scatterer in each differential interferogram in the differential interferogram sequence for candidate point screening, and perform an arithmetic mean of all coherence coefficients to obtain the temporal average coherence of the candidate distributed scatterer.

[0112] As is understandable, the coherence coefficient is an indicator of interferogram quality, reflecting the degree of phase consistency between two images, and can be used to evaluate the phase stability of candidate points for distributed scatterers (DSS). The two images are a pair of single-polarization synthetic aperture radar (SAP) images used to generate the k-th differential interferogram, including a primary image as a reference and a secondary image registered with it. These two images are selected from the SAP image sequence and form the interferometric pair. The coherence coefficient of each candidate point in each differential interferogram can be calculated using the maximum likelihood estimation method based on a homogeneous pixel set. To obtain the overall coherence level of DSS candidate points in the time dimension, the coherence coefficients in all differential interferograms can be arithmetically averaged, thus determining the temporal average coherence of each DSS candidate point. This value represents the average coherence level of the DSS candidate point over the entire monitoring period.

[0113] For example, for each candidate point of the distributed scatterer, the coherence coefficient of the candidate point in the k-th differential interferogram can be expressed as:

[0114] ;

[0115] in, It is represented as the coherence coefficient of candidate point P of the distributed scatterer in the k-th differential interferogram. It is a dimensionless value between 0 and 1, used to measure the stability of the phase of the candidate point in the k-th interferogram. It represents the total number of pixels in the homogeneous pixel set of the candidate point P of the distributed scatterer, that is, the number of statistical samples participating in the coherence estimation; This represents the complex observation value of the p-th homogeneous pixel in the first single-polarization synthetic aperture radar image used to generate the k-th differential interferogram, including amplitude and phase information; It represents the complex observation value of the p-th homogeneous pixel in the second single-polarization synthetic aperture radar image used to generate the k-th differential interferogram; Represented as The conjugate of complex numbers.

[0116] The temporal average coherence of candidate points of the distributed scatterer is obtained by taking the arithmetic mean of the coherence coefficients of all differential interferograms, and its expression can be expressed as:

[0117] ;

[0118] in, The time-series average coherence of candidate point P of the distributed scatterer reflects the average coherence level of the candidate point throughout the entire monitoring period; M represents the total number of differential interferograms used for coherence calculation, i.e., the number of interferograms participating in the time-series averaging.

[0119] S502, based on the vegetation coverage classification atlas, the temporal average coherence of all candidate points of distributed scatterers within each vegetation coverage level is statistically analyzed, and the regional spatial mean of each vegetation coverage level is calculated respectively.

[0120] Here, the regional spatial mean is the arithmetic mean of the temporal average coherence of all candidate points of distributed scatterers within the same vegetation cover level, which can be expressed as the average coherence level of candidate points of distributed scatterers within that level. By calculating the regional spatial mean of each vegetation cover level, the overall coherence distribution characteristics of candidate points of distributed scatterers under different vegetation conditions can be quantified, so as to determine the benchmark level of coherence for each level.

[0121] For example, the regional spatial mean for each vegetation cover level can be expressed as:

[0122] ;

[0123] in, It is represented as the regional spatial mean of the I-th vegetation cover level, used to characterize the average coherence level of all candidate points of distributed scatterers within that vegetation cover level; It represents the total number of candidate points for distributed scatterers within the I-th vegetation cover level.

[0124] S503: Calculate the temporal average coherence of all candidate points of distributed scatterers and the global spatial mean.

[0125] Here, the global spatial mean refers to the arithmetic mean of the temporal average coherence of all distributed scatterer candidate points within the entire target monitoring mountain area. It is used to characterize the overall coherence level of distributed scatterer candidate points in the entire monitoring area and serves as a reference benchmark for subsequent threshold setting. The global spatial mean can be calculated by statistically analyzing the temporal average coherence of all distributed scatterer candidate points and performing an arithmetic mean.

[0126] For example, the global spatial mean of all candidate points of the distributed scatterer can be expressed as:

[0127] ;

[0128] in, It is represented as the global spatial mean, used to characterize the overall coherence level of all distributed scatterer candidate points within the entire target monitoring mountain area; This represents the total number of candidate points for the distributed scatterer.

[0129] S504, for each vegetation coverage level, the smaller value between the regional spatial mean and the global spatial mean of the vegetation coverage level is determined as the coherence screening threshold of the vegetation coverage level.

[0130] Specifically, since the surface scattering characteristics of different vegetation cover areas are different, the coherence levels of candidate distributed scatterers will also vary significantly. Therefore, in order to achieve reasonable monitoring point selection under different vegetation conditions, an appropriate coherence screening threshold can be determined for each vegetation cover level. This can be achieved by comparing the regional spatial mean of each level with the global spatial mean and taking the smaller value as the screening threshold for that level, thereby obtaining an adaptive threshold for each level for subsequent distributed scatterer screening.

[0131] In this way, by taking the smaller value between the regional spatial mean and the global spatial mean as the screening threshold, a relatively lenient global mean standard can be used in high coherence areas with low vegetation coverage to retain a sufficient number of monitoring points, while the mean standard of the area itself can be used in low coherence areas with high vegetation coverage to avoid excessive removal of monitoring points due to an excessively high threshold. This enables a dynamic balance between the number and quality of monitoring points under different vegetation conditions.

[0132] For example, the coherence screening threshold for each vegetation cover level can be determined by comparing the regional spatial mean of that level with the global spatial mean and taking the smaller value. The process can be expressed as follows:

[0133] ;

[0134] in, This represents the coherence screening threshold for the I-th vegetation cover level, used to screen candidate points of distributed scatterers within that level.

[0135] S505, candidate points of distributed scatterers whose time-series average coherence is greater than or equal to the coherence screening threshold corresponding to the vegetation coverage level are identified as distributed scatterers.

[0136] Furthermore, after obtaining the temporal average coherence of each candidate point and the coherence screening threshold for each level, the temporal average coherence of each candidate point can be compared with the coherence screening threshold of its corresponding vegetation cover level. By judging whether the temporal average coherence is greater than or equal to the threshold, it serves as the basis for screening distributed scatterers. Distributed scatterers have high phase stability and reliability and can be used for subsequent high-precision deformation inversion.

[0137] In the above embodiments, by introducing vegetation coverage classification and adopting an adaptive coherence threshold screening strategy, the screening criteria can be dynamically adjusted according to the coherence distribution characteristics of different vegetation coverage areas. While ensuring the quality of monitoring points, the density of monitoring points in dense vegetation areas can be significantly improved, thereby obtaining distributed scatterers with more complete spatial coverage and more uniform distribution.

[0138] S104, based on the multiple distributed scatterers, the single-polarization synthetic aperture radar image sequence, and the digital elevation model data, determine the long-term surface deformation results of the target monitoring mountain area.

[0139] Understandably, the long-term surface deformation results for target monitoring mountainous areas refer to the displacement information of the target area's surface along the radar line of sight during the entire monitoring period. This information can reveal the spatiotemporal evolution characteristics of geological hazards such as volcanic activity, mining subsidence, or landslide deformation, and can include surface deformation distribution maps and time-series deformation curves. The surface deformation distribution map is a graphical representation of the spatial distribution of the average deformation rate or cumulative deformation across the entire monitoring area, visually reflecting the geographical extent and intensity of deformation. The time-series deformation curve is a line graph showing the cumulative deformation of a single distributed scatterer over time during the monitoring period, illustrating the deformation evolution process at that point. By performing phase unwrapping and deformation inversion based on selected high-quality distributed scatterers, time-series analysis of the phase information of each distributed scatterer can be conducted, thereby obtaining deformation information covering the entire monitoring area and achieving accurate monitoring of surface deformation in complex mountainous areas.

[0140] In some possible embodiments, in order to obtain more accurate deformation results and reduce the interference of noise and errors on deformation inversion, when performing deformation inversion, a reference is made to... Figure 6 As shown, the steps S601 to S605 may be included:

[0141] S601, Based on the single-polarization synthetic aperture radar image sequence and the digital elevation model data, an initial differential interferogram sequence is generated.

[0142] It is understandable that deformation inversion also needs to be based on differential interferograms. The specific process of generating the initial differential interferogram based on the single-polarization synthetic aperture radar image sequence and digital elevation model data can be referred to the description of the initial differential interferogram generation in the above steps. The specific process has been clearly explained in the above steps and will not be repeated here.

[0143] S602, based on the plurality of distributed scatterers and the homogeneous pixel set of the plurality of distributed scatterers, the initial differential interferogram sequence is subjected to a second adaptive filtering process to obtain a differential interferogram sequence for deformation inversion.

[0144] Here, since the distributed scatterers have already undergone vegetation cover grading constraints and adaptive coherence threshold screening, they are reliable monitoring points. In order to further improve the quality of the interferograms used for deformation inversion, a second adaptive filtering is performed on the initial differential interferogram based on the selected distributed scatterers and their homogeneous pixel set. This allows the statistical information of the distributed scatterers and their homogeneous pixels to guide the determination of the filtering parameters, making the filtering process more accurate and resulting in a higher quality differential interferogram sequence that is more suitable for deformation inversion.

[0145] The second adaptive filtering can select an appropriate filtering method according to actual needs. It can use the same filtering method as the first adaptive filtering (such as Goldstein filtering, adaptive window filtering, or improved Goldstein filtering), or it can use different filtering methods according to the data characteristics. No specific restrictions are made here.

[0146] S603, perform phase unwrapping on the differential interferogram sequence used for deformation inversion to obtain an unwrapped phase map; and perform time-series deformation inversion processing on the unwrapped phase map to obtain the deformation rate and cumulative deformation of each distributed scatterer.

[0147] It is understandable that phase unwrapping refers to the process of recovering the true phase value from the interference phase wrapped between -π and π. By performing spatial phase recovery on each unwrapped map in the differential interferogram sequence for deformation inversion, the phase ambiguity problem can be eliminated. The unwrapped phase map represents the spatial distribution image of the true phase value of each pixel, which can reflect the actual phase change caused by surface deformation.

[0148] Furthermore, after obtaining the unwrapped phase map, it can be subjected to time-series deformation inversion processing. Using methods such as least squares or singular value decomposition, multiple unwrapped phase maps can be modeled and solved in the time dimension. By constructing a linear relationship between deformation rate and phase change and estimating parameters, the deformation rate and cumulative deformation of each distributed scatterer can be obtained. The deformation rate represents the average displacement of the distributed scatterer per unit time, which can be used to assess the severity of deformation; the cumulative deformation refers to the total displacement from the start of monitoring to a certain point in time, reflecting the cumulative effect of deformation over time.

[0149] S604 performs atmospheric phase correction and orbital error removal on the deformation rate and cumulative deformation of each distributed scatterer obtained by inversion to obtain the corrected deformation information.

[0150] Here, the initial inversion results may contain non-deformation components such as atmospheric delay phase and orbital residual errors. Therefore, in order to obtain purer deformation information, after obtaining the deformation rate and cumulative deformation of each distributed scatterer, atmospheric phase correction and orbital error removal can be performed to optimize the accuracy of the deformation information and obtain the corrected deformation information, which mainly includes the corrected deformation rate and cumulative deformation.

[0151] Specifically, atmospheric phase correction is a commonly used deformation refinement method that analyzes the spatial and temporal distribution characteristics of atmospheric delay and uses filtering or modeling methods to separate and remove atmospheric contributions from deformation information. Orbit error removal can eliminate residual phase slope caused by orbit inaccuracies through polynomial fitting or ground control point-based methods. Here, the method for correcting deformation information can be selected based on the size of the monitoring area and the data characteristics, choosing an appropriate atmospheric correction model and orbit error removal method; no specific limitations are imposed here.

[0152] S605, the corrected deformation information is converted from the radar coordinate system to the geographic coordinate system to generate the long-term surface deformation results of the target monitoring mountain area; wherein, the long-term surface deformation results include a surface deformation distribution map and a time-series deformation curve.

[0153] Furthermore, since the corrected deformation information is still in the radar coordinate system, geocoding can be performed after the correction is completed to transform the deformation information from the radar coordinate system to the geographic coordinate system, so as to generate long-term surface deformation results with geographic reference significance, which can be used for geological disaster risk assessment and monitoring report output, thereby completing the interferometric measurement task of target monitoring mountainous areas based on distributed scatterers, and realizing high-precision monitoring of surface deformation in complex mountainous areas.

[0154] For example, in the deformation inversion process, the minimum cost flow method can be used to unwrap the differential interferogram sequence used for deformation inversion. The minimum cost flow method is a phase unwrapping algorithm based on network optimization theory. By constructing a network flow model for phase unwrapping, it minimizes the accumulation of unwrapping errors while satisfying the phase gradient consistency constraint. When applying this method, a coherence threshold can be set to mask low-coherence areas, removing pixels with low coherence due to snow cover and remaining vegetation cover, thus avoiding interference from these areas with the unwrapping results. Unwrapped interferograms of ascending and descending orbit data are then obtained to form an unwrapped phase map covering the monitoring area.

[0155] Furthermore, residual orbital errors and trend errors caused by long-wave atmospheric delay can be removed from the unwrapped results using a quadratic polynomial fitting method, eliminating the contamination of deformation information by systematic errors. Then, the unwrapped interferograms after error removal are quality-screened, and interferograms of poor quality are discarded based on indicators such as residual magnitude and phase continuity. Finally, an annual average deformation rate map can be obtained using a phase stacking algorithm. This algorithm suppresses atmospheric noise by weighted averaging of multiple unwrapped interferograms to obtain the regional average deformation rate; alternatively, the regional cumulative deformation can be obtained using a small baseline subset interferometry method. This method uses a combination of small baseline interferometrics from multiple master images to construct a system of equations, and solves the deformation time series through singular value decomposition. Geocoding the inversion results transforms the deformation information from the radar coordinate system to the geographic coordinate system, yielding long-term surface deformation results in the WGS-84 coordinate system, including surface deformation distribution maps and time-series deformation curves.

[0156] The following example, using the Changshan Mountain area, illustrates the adaptive distributed scatterer interferometry method for mountainous regions provided in this disclosure. This embodiment selects the Changshan Mountain volcanic region in northeastern my country as the target monitoring area. This region has dense vegetation cover, with a forest coverage exceeding 80%, exhibiting significant vertical vegetation zonation, transitioning sequentially from low-altitude mixed coniferous and broad-leaved forests to high-altitude alpine meadows and bare rocks. This presents a severe spatiotemporal incoherence problem for synthetic aperture radar (SAR) interferometry deformation monitoring. Traditional temporal SAR interferometry methods use a global single threshold to select monitoring points, resulting in sparse monitoring points and frequent coherence voids in mid-to-high vegetation cover areas, making it difficult to capture minute deformation signals.

[0157] To conduct research on time-series InSAR methods applicable to complex mountainous areas, this publication collects two sets of C-band Sentinel SAR datasets: 211 Sentinel-1A rising orbit images from January 2020 to December 2023, and 58 Sentinel-1B falling orbit images from January 2020 to December 2021. The main parameters of these data include: both Sentinel-1A and Sentinel-1B satellites have an orbital altitude of 693 km and an orbital inclination of 98.16°, and both use the C-band (wavelength approximately 5.55 nm). The imaging modes are all 5.405 GHz (cm), with a regression period of 12 days. All imaging modes are wide-swath interferometric modes, with a range resolution of 5 m and an azimuth resolution of 20 m. Sentinel-1A is an ascending orbit with a heading angle of -13.72° and an incident angle of 39.0984°, containing 211 images covering the period from January 2020 to December 2023. Sentinel-1B is a descending orbit with a heading angle of 193.73° and an incident angle of 42.3223°, containing 58 images covering the period from January 2020 to December 2021. Furthermore, the digital elevation model data used in this disclosure is 90 m resolution TanDEM-XDEM, which, compared to traditional InSAR systems, avoids phase noise caused by various factors such as decorrelation, atmospheric delay, and orbital errors.

[0158] In the specific implementation process, this disclosure first performs radiometric calibration and polarization calibration on the acquired sentinel images, converting the vertical and horizontal polarization backscattering coefficients from decibel values ​​to linear values; based on the threshold of vertical polarization backscattering coefficients being below -17 dB, a water body mask is constructed to remove water body pixels such as crater lakes; a 30 m × 30 m median filter is used to denoise the images, reducing noise interference while preserving vegetation scattering characteristics.

[0159] Furthermore, for the dual-polarization synthetic aperture radar image sequence, the radar vegetation index of each pixel in each image is calculated. Then, according to the preset radar vegetation index threshold range, the vegetation coverage is divided into six levels: bare land, low coverage, medium-low coverage, medium coverage, medium-high coverage, and high coverage. Specifically, a radar vegetation index less than 0.45 is classified as bare land, greater than or equal to 0.45 and less than 0.55 as low coverage, greater than or equal to 0.55 and less than 0.65 as medium-low coverage, greater than or equal to 0.65 and less than 0.75 as medium coverage, greater than or equal to 0.75 and less than 0.85 as medium-high coverage, and greater than or equal to 0.85 and less than or equal to 1 as high coverage. Then, the vegetation coverage level of each pixel in each image is labeled in the registered radar coordinate system to generate a vegetation coverage grading map corresponding to each image. The vegetation coverage grading maps of all images are then combined according to the image time sequence to obtain a vegetation coverage grading map set corresponding to the dual-polarization synthetic aperture radar image sequence.

[0160] Reference Figure 7 As shown, the vegetation cover data for the monitoring area is as follows: 1 represents bare land, 2 represents low vegetation cover, 3 represents medium-low vegetation cover, 4 represents medium vegetation cover, 5 represents medium-high vegetation cover, and 6 represents high vegetation cover. The data shows that the average vegetation cover in this monitoring area ranges from 0 to 0.962, exhibiting a radial distribution pattern with higher cover around the perimeter and lower cover in the center. Areas with medium-high vegetation cover account for over 94%, while bare land and low vegetation cover areas account for only about 5%, mainly concentrated in the crater and surrounding mountains. There are significant differences in vegetation on the east and west sides of the crater. The east side has concentrated volcanic debris accumulation and poor soil development, dominated by dwarf alpine mosses and meadows; the west side has gentle slopes, stable surfaces, and a continuous transition in vegetation cover. The southern and western outer mountain areas have gentle terrain and fertile soil, dominated by mixed coniferous and broad-leaved forests, with an average vegetation cover of 0.76–0.88, locally approaching 0.94. The patchy bare land on the eastern edge of the study area corresponds to sparsely vegetated mountains, villages, and farmland, consistent with optical imagery verification, confirming the good applicability of the dual-polarization radar vegetation index (RVI). The southern mountainous area, due to its steep terrain and intense erosion, is mostly continuous bare land; some mountains exhibit significant differences in vegetation cover between sunny and shady slopes due to slope aspect. The piedmont transition zone, with its gentle terrain and well-developed drainage system, forms strip-shaped areas of high vegetation cover. Therefore, the vegetation cover in the Changshan volcanic area is dominated by topography and vegetation type, exhibiting a clear vertical zonation. High-value areas are concentrated in gently sloping areas with abundant water and heat, while low-value areas are distributed in densely populated areas and steep mountains.

[0161] Here, after labeling vegetation cover levels, the vegetation cover classification results can be used as constraints for homogeneous pixel identification. During homogeneous pixel identification, a 9×9 (azimuth × distance) detection window is set, with a significance level of 0.05, and the homogeneous pixel number threshold is calculated to be 16. To evaluate the effectiveness of the proposed homogeneous pixel identification and distributed scatterer selection method with vegetation cover as a constraint, Sentinel-1A ascending orbit image covering a volcanic area of ​​Changshan Mountain is used as experimental data. The homogeneous pixel selection method combined with vegetation cover is compared with the traditional fast statistical homogeneous pixel identification method and the Baumgartner-Weiβ-Schindler test method. The spatial distribution of homogeneous pixels obtained by the three methods is as follows: Figure 8 As shown, the three homogeneous pixel recognition methods exhibit significant differences in performance: Figure 8 (a) is the result obtained by the Baumgartner-Weiβ-Schindler test method. The method extracts a small number of homogeneous pixels with discontinuous spatial distribution. There are obvious discontinuities in the terrain transition area, which makes it difficult to support the subsequent deformation details. Figure 8 (b) is the result image extracted using a rapid statistical homogeneous pixel identification method. While this method increases the number of pixels, excessive smoothing leads to a large Type II error, resulting in numerous misidentifications in areas such as bare rock in volcanic craters. The result image extracted using the vegetation cover constraint method proposed in this disclosure is... Figure 8 (c) Although the number of pixels is slightly less than that of the fast statistical homogeneous pixel identification method, it closely matches the characteristics of real land cover, showing a continuous distribution in areas with the same vegetation level. The number of pixels in the transition zone between bare land and vegetation is significantly reduced, effectively avoiding cross-type mismatches. This method can also remove water body pixels from crater lakes, improving the rationality of the set, and has the best computational efficiency (only 1.53546 seconds), which is much faster than the Baumgartner-Weiβ-Schindler test method (720 seconds) and the fast statistical homogeneous pixel identification method (22.61637 seconds).

[0162] Specifically, the homogeneous pixel identification results under different vegetation coverage levels show that the method provided in this disclosure can select homogeneous pixels in a targeted manner at each vegetation coverage level, taking into account both quantity and quality, and providing more reliable basic data for subsequent high-precision deformation processing. The overall stability and effectiveness are significantly better than traditional methods. The data includes: 221,813 homogeneous pixels in bare land areas with a regional coherence of 0.28419; 574,382 homogeneous pixels in low vegetation cover areas with a regional coherence of 0.27576; 2,678,386 homogeneous pixels in low to medium vegetation cover areas with a regional coherence of 0.21691; 4,572,872 homogeneous pixels in medium vegetation cover areas with a regional coherence of 0.18307; 269,179 homogeneous pixels in medium to high vegetation cover areas with a regional coherence of 0.16338; and 173,310 homogeneous pixels in high vegetation cover areas with a regional coherence of 0.15842. The above data shows that this method can select a sufficient number of homogeneous pixels in areas with high coherence, such as bare land, to ensure quality, and can also maintain a certain number of homogeneous pixels in densely vegetated areas to fill monitoring gaps, thus achieving a dynamic balance between the number and quality of monitoring points under different vegetation cover conditions.

[0163] To more intuitively compare the differences among the three homogeneous pixel identification methods, pixel coordinates [104, 469] were selected based on the temporal average intensity map and optical imagery of the study area. A reference point was chosen at the boundary of different vegetation cover levels. Within an 800m × 800m window, centered on this pixel, the extraction results of the Baumgartner-Weiβ-Schindler test, the fast statistical homogeneous pixel identification method, and the homogeneous pixel identification method constrained by vegetation cover were compared. Figure 9 As shown. The reference pixel is mainly bare land, while the neighborhood includes woodland, sparse vegetation, and bare land on ridges, indicating significant heterogeneity in scattering mechanisms. The Baumgartner-Weiβ-Schindler test method (i.e., Figure 9 (a) and fast statistical homogeneous pixel recognition method (i.e. Figure 9 (b) The recognition effect is similar, but only sporadic homogeneous pixels are extracted on the slope side where the reference pixel is located, resulting in poor continuity and failure to capture the linear features of local bare land and surface roads. In contrast, the method proposed in this disclosure (such as...) Figure 9 (c) This method extracts the most homogeneous samples, resulting in more complete spatial coverage. Pixels are distributed only within the same bare land level as the reference pixels, closely matching the actual land cover boundaries and avoiding misselection across vegetation types. The method also successfully extracts continuous homogeneous pixels along road directions and identifies stable pixel sets consistent with the true distribution in valleys and slopes. Experiments further validate the method's excellent applicability and stability in volcanic environments with complex land cover types, accurately capturing homogeneous pixels on heterogeneous surfaces.

[0164] After selecting candidate distributed scatterers (DSARs), the coherence of each candidate DSAR is calculated. The study area is divided into zones based on the surface vegetation cover level, and the coherence threshold for each zone is calculated and compared with the global coherence. A secondary screening is then performed on all candidate targets to obtain all DSARs selected based on the adaptive threshold. To accurately measure the effectiveness of the adaptive threshold-based DSAR selection, candidate DSARs with a time-averaged coherence greater than the corresponding level threshold are selected as the final monitoring points, ensuring a balance in the quality and density of monitoring points across different vegetation areas. The coherence threshold is set based on the global average coherence (approximately 0.21392). The experimental area is expanded to more clearly illustrate the spatial distribution differences of the DSARs selected by the two methods. The results of the DSARs selected using the fixed coherence threshold-based synthetic aperture radar interferometry method and the adaptive threshold method are shown below. Figure 10 As shown in the figure. Referring to the content in the figure, it can be seen that the spatial distribution obtained by the fixed threshold method... Figure 10 The distributed scatterer point density of B is unevenly distributed, concentrated in high coherence areas, with no points observed in areas with high vegetation cover or complex scattering. The volcanic core area has insufficient cover, making it difficult to characterize deformation details. A total of approximately 1,026,000 points were selected. The spatial distribution obtained by the adaptive thresholding method proposed in this disclosure... Figure 10 Method A dynamically adjusts the threshold according to vegetation cover, resulting in a regular, uniform, and continuous distribution of distributed scatterer points that can fit complex surface features. A total of 3.927 million points were selected, nearly four times that of the fixed threshold method. This method sets thresholds by zoning vegetation cover areas, avoiding the rejection of low-coherence but stable scatterers, and significantly improving the retention rate of distributed scatterer points while ensuring quality.

[0165] Finally, the deformation results of the study area were extracted using an adaptive distributed scatterer synthetic aperture radar interferometry method based on vegetation cover. Figure 11 This study presents the surface deformation rate of the Changshan volcanic area from December 2020 to July 2023. The results show that the crater and its surrounding area are generally stable (deformation rate ±8 mm / year). Localized patchy deformation exists in the transition zone on the southwest and east slopes of the volcano, mostly caused by small-scale landslides or residual atmospheric errors. To quantitatively analyze the long-term deformation trend of the Changshan volcanic area, targets at six typical monitoring points around the volcano were also mapped. Figure 11The time-series deformation curves of P1~P6 show fluctuating characteristics at all six monitoring points, with overall deformation ranging from -10 mm to 10 mm. Monitoring points near the crater exhibit smaller deformation fluctuations, indicating overall stability in the crater area. Monitoring points located in the transition zone between the southwest and east slopes of the volcano show more significant deformation fluctuations, with some periods exhibiting cumulative deformation trends. These trends correspond to the higher-than-background seismic activity during the second low-level disturbance in 2020. To clarify the deformation characteristics and causes and verify the reliability of the method, this experiment selected four severely locally deformed areas for analysis. Figure 12 ). Figure 12 a1、 Figure 12 a2 is located in a transitional zone with significant topographic changes. The soil and rock are loose and eroded, and the deformation of the rising and falling tracks is consistent, indicating that it was caused by a landslide. Figure 12 a3、 Figure 12 Located on the crater cone and ridge slope, area a4 is prone to weathering and landslides due to volcanic clastic rocks. The deformation trend matches the characteristics, further confirming the reliability of the method. Furthermore, this experiment conducted statistical analysis on the monitoring points in the four areas mentioned above, finding that the adaptive distributed scatterer synthetic aperture radar interferometry method can still provide relatively rich small-area deformation information even when conducting large-scale monitoring, increasing the number of monitoring points by approximately 30% compared to the original method. In areas with no or small deformation, the deformation levels of the original and adaptive methods are similar, almost identical; however, in areas with large deformation, there are significant differences in the deformation levels. Figure 13 It can be seen that, regarding the above Figure 12 a1、 Figure 12 a2、 Figure 12 a3、 Figure 12 In the extraction results of a4, the deformation distribution features extracted by the original method are not significantly different from those of the adaptive method proposed in this disclosure, mainly manifested in the lack of prominent differences in color distribution in the image. Specifically, Figure 13 (a) indicates as to Figure 12 The surface deformation results for region a1 were extracted using the method proposed in this disclosure based on the elevated orbit data. Figure 13 (b) indicates as to Figure 12 The surface deformation results for region a1 were extracted using the method proposed in this disclosure based on the down-track data; Figure 13 (c) indicates that in Figure 12 The surface deformation results for region a1 were extracted using the original method and elevated orbit data. Figure 13 (d) indicates that... Figure 12 The surface deformation results for region a2 were extracted using the method proposed in this disclosure based on the elevated orbit data. Figure 13 (e) indicates that... Figure 12 The surface deformation results for region a2 were extracted using the method proposed in this disclosure based on the down-track data. Figure 13 (f) indicates that in Figure 12 The surface deformation results for region a2 were extracted using the original method and elevated orbit data. Figure 13 (g) indicates that... Figure 12 The surface deformation results for region a3 were extracted using the method proposed in this disclosure based on the elevated orbit data. Figure 13 (h) indicates that... Figure 12 The surface deformation results for region a3 were extracted using the method proposed in this disclosure based on the down-track data. Figure 13 (i) indicates that in Figure 12 The surface deformation results for region a3 were extracted using the original method and elevated orbit data. Figure 13 (j) represents the pair Figure 12 The surface deformation results for region a4 were extracted using the method proposed in this disclosure based on the elevated orbit data. Figure 13 (k) represents the pair Figure 12 The surface deformation results for region a4 were extracted using the method proposed in this disclosure based on the down-track data. Figure 13 (l) indicates that in Figure 12 The surface deformation results for region a4 were extracted using the original method and elevated orbit data.

[0166] This experiment selected 940,763 identical geographically located measurement points of the same name to compare the deformation rates obtained by the original method and the method of this application. The standard deviation and root mean square error of the deformation rate error were calculated to be 12.4296 mm and 14.2161 mm, respectively. These results indirectly verify the reliability of the adaptive distributed scatterer synthetic aperture radar interferometry method combined with vegetation cover in monitoring surface deformation in the Changshan volcanic area. Therefore, based on the embodiment of deformation monitoring in the Changshan area, compared with the original method, the adaptive method proposed in this disclosure can provide richer spatiotemporal deformation details, better characterize small-area deformation features, and identify and locate local deformation areas over a large scale.

[0167] The adaptive distributed scatterer interferometry method, device, and medium for mountainous areas provided in this embodiment introduce radar vegetation index to classify vegetation coverage of dual-polarization SAR images in mountainous areas, and on this basis, implement statistical homogeneous pixel identification and adaptive coherence threshold screening under spatial constraints. This can extract denser and more reliable distributed scatterers, thereby effectively improving the spatial continuity and solution accuracy of long-term surface deformation monitoring in mountainous areas.

[0168] Those skilled in the art will understand that, in the above-described method of the specific implementation, the order in which each step is written does not imply a strict execution order and does not constitute any limitation on the implementation process. The specific execution order of each step should be determined by its function and possible internal logic.

[0169] Based on the same inventive concept, this disclosure also provides a mountain adaptive distributed scatterer interferometry device corresponding to the mountain adaptive distributed scatterer interferometry method. Since the principle of the device in this disclosure is similar to the mountain adaptive distributed scatterer interferometry method described above, the implementation of the device can refer to the implementation of the method, and the repeated parts will not be described again.

[0170] Reference Figure 14 The diagram shown is a schematic of an adaptive distributed scatterer interferometry device 1400 for mountainous areas provided in an embodiment of this disclosure. The device includes:

[0171] The data acquisition module 1401 is used to acquire radar image data and digital elevation model data of the mountainous area for target monitoring; wherein, the radar image data includes a single-polarization synthetic aperture radar image sequence and a dual-polarization synthetic aperture radar image sequence, the dual-polarization synthetic aperture radar image sequence includes multiple dual-polarization SAR images, and each dual-polarization SAR image includes multiple pixels.

[0172] The pixel classification module 1402 is used to calculate the radar vegetation index of each pixel in each dual-polarization synthetic aperture radar image sequence, and divide the pixels in each dual-polarization synthetic aperture radar image into multiple vegetation coverage levels based on the calculation results, so as to obtain a vegetation coverage classification map set corresponding to the dual-polarization synthetic aperture radar image sequence.

[0173] The scatterer identification module 1403 is used to identify statistically homogeneous pixels in each pixel of each dual-polarization SAR image under the spatial constraint of vegetation coverage level according to the vegetation coverage classification atlas, and to determine the distributed scatterer candidate point set corresponding to the dual-polarization synthetic aperture radar image sequence based on the identification results; and to perform adaptive coherence threshold screening on each distributed scatterer candidate point in the distributed scatterer candidate point set according to the vegetation coverage classification atlas to obtain multiple distributed scatterers.

[0174] The result determination module 1404 is used to determine the long-term surface deformation results of the target monitoring mountain area based on the multiple distributed scatterers, the single-polarization synthetic aperture radar image sequence, and digital elevation model data.

[0175] In some possible embodiments, the pixel classification module 1402 is specifically used for:

[0176] For each dual-polarization SAR image, the vertical polarization backscattering coefficient and the horizontal polarization backscattering coefficient of each pixel in the dual-polarization SAR image are extracted respectively.

[0177] For each pixel, the radar vegetation index of the pixel is calculated based on the vertical polarization backscattering coefficient and the horizontal polarization backscattering coefficient corresponding to the pixel.

[0178] Based on multiple preset radar vegetation index threshold ranges, pixels whose radar vegetation index values ​​fall within the same threshold range are classified into the same vegetation coverage level.

[0179] The vegetation coverage level of each pixel in each dual-polarization SAR image is labeled in the radar coordinate system to generate a vegetation coverage classification map corresponding to each dual-polarization SAR image; and the vegetation coverage classification maps corresponding to all dual-polarization SAR images are combined to obtain the vegetation coverage classification map set.

[0180] In some possible embodiments, the scatterer identification module 1403 is specifically used for:

[0181] A water body mask is constructed based on a preset backscattering coefficient threshold, and water body pixels in each dual-polarization SAR image are removed using the water body mask.

[0182] For each non-water body pixel, a neighborhood window of a preset size is set with the non-water body pixel as the center, and neighborhood pixels belonging to the same vegetation coverage level as the non-water body pixel are filtered within the neighborhood window.

[0183] For each non-water body pixel, a statistical homogeneity test is performed on the selected neighboring pixels corresponding to the non-water body pixel, and the neighboring pixels that pass the test and are spatially connected to the non-water body pixel are determined as homogeneous pixels of the non-water body pixel.

[0184] The number of homogeneous pixels corresponding to each non-water body pixel is counted, and non-water body pixels with a number of homogeneous pixels greater than a preset homogeneous pixel threshold are identified as candidate points for distributed scatterers.

[0185] The candidate distributed scatterers identified in all dual-polarization SAR images are combined according to the image time sequence to obtain the candidate distributed scatterer set.

[0186] In some possible embodiments, the scatterer identification module 1403 is further configured to:

[0187] For each selected neighboring pixel corresponding to the non-water body pixel, the backscattering coefficients of the neighboring pixel and the non-water body pixel are extracted in all phases of the dual-polarization synthetic aperture radar image sequence, and the neighborhood temporal amplitude sequence of the neighboring pixel and the center temporal amplitude sequence of the non-water body pixel are constructed respectively.

[0188] A two-sample position test is performed on the neighborhood temporal amplitude sequence of the neighboring pixels and the center temporal amplitude sequence of the non-water body pixels to determine whether the two sequences come from the same continuous distribution.

[0189] If the test statistic is less than the preset threshold, the neighboring pixel is determined to have passed the statistical homogeneity test; otherwise, it is determined to have failed.

[0190] Among the neighboring pixels that pass the test, the neighboring pixels that are directly adjacent to the non-water body pixel, and the neighboring pixels that are indirectly connected to the non-water body pixel through other neighboring pixels that pass the test, are identified as homogeneous pixels that are spatially connected to the non-water body pixel.

[0191] In some possible embodiments, the scatterer identification module 1403 is specifically used for:

[0192] Based on the single-polarization synthetic aperture radar image sequence and the digital elevation model data, an initial differential interferogram sequence is generated;

[0193] The initial differential interferogram sequence is subjected to a first adaptive filtering process to obtain a differential interferogram sequence for candidate point screening, which is used to calculate the coherence coefficient of each candidate point of the distributed scatterer;

[0194] The generation of the initial differential interferogram sequence includes:

[0195] Based on a preset spatiotemporal baseline threshold, image pairs that satisfy the preset spatiotemporal baseline threshold are selected from the single-polarization synthetic aperture radar image sequence to form interferometric pairs. Interferometric processing is performed on each interferometric pair to obtain an interferogram. The terrain phase in the interferogram is removed using the digital elevation model data to obtain the initial differential interferogram sequence.

[0196] In some possible embodiments, the scatterer identification module 1403 is specifically used for:

[0197] For each candidate distributed scatterer, based on the homogeneous pixel set of the candidate distributed scatterer, the coherence coefficient of the candidate distributed scatterer in each differential interferogram in the differential interferogram sequence for candidate point screening is calculated, and the arithmetic mean of all coherence coefficients is performed to obtain the temporal average coherence of the candidate distributed scatterer.

[0198] Based on the vegetation cover classification atlas, the temporal average coherence of all candidate points of distributed scatterers within each vegetation cover level is statistically analyzed, and the regional spatial mean of each vegetation cover level is calculated respectively.

[0199] Statistically calculate the temporal average coherence of all candidate points of distributed scatterers and then calculate the global spatial mean.

[0200] For each vegetation cover level, the smaller of the regional spatial mean and the global spatial mean of the vegetation cover level is determined as the coherence screening threshold of the vegetation cover level.

[0201] Distributed scatterer candidate points whose temporal average coherence is greater than or equal to the coherence screening threshold corresponding to the vegetation coverage level are identified as distributed scatterers.

[0202] In some possible embodiments, the result determination module 1404 is specifically used for:

[0203] Based on the single-polarization synthetic aperture radar image sequence and the digital elevation model data, an initial differential interferogram sequence is generated;

[0204] Based on the plurality of distributed scatterers and the homogeneous pixel set of the plurality of distributed scatterers, the initial differential interferogram sequence is subjected to a second adaptive filtering process to obtain a differential interferogram sequence for deformation inversion;

[0205] The deformation inversion sequence is phase-unwrapped to obtain an unwrapped phase map; and the unwrapped phase map is subjected to time-series deformation inversion processing to obtain the deformation rate and cumulative deformation of each distributed scatterer.

[0206] Atmospheric phase correction and orbital error removal are performed on the deformation rate and cumulative deformation of each distributed scatterer obtained by inversion to obtain the corrected deformation information;

[0207] The corrected deformation information is converted from the radar coordinate system to the geographic coordinate system to generate long-term surface deformation results for the target monitoring mountain area; wherein, the long-term surface deformation results include a surface deformation distribution map and a time-series deformation curve.

[0208] Based on the same technical concept, this disclosure also provides a computer device. (See also...) Figure 15 The diagram shown is a structural schematic of a computer device 1500 provided in an embodiment of this disclosure, including a processor 1501, a memory 1502, and a bus 1503. The memory 1502 stores execution instructions and includes a main memory 15021 and an external memory 15022. The main memory 15021, also called internal memory, is used to temporarily store computational data in the processor 1501 and data exchanged with external memory 15022 such as a hard disk. The processor 1501 exchanges data with the external memory 15022 through the main memory 15021.

[0209] In this embodiment, the memory 1502 is specifically used to store application code that executes the solution of this application, and its execution is controlled by the processor 1501. That is, when the computer device 1500 is running, the processor 1501 communicates with the memory 1502 through the bus 1503, so that the processor 1501 executes the application code stored in the memory 1502, and then executes the method described in any of the foregoing embodiments.

[0210] The memory 1502 may be, but is not limited to, random access memory (RAM), read-only memory (ROM), programmable read-only memory (PROM), erasable programmable read-only memory (EPROM), electrically erasable programmable read-only memory (EEPROM), etc.

[0211] Processor 1501 may be an integrated circuit chip with signal processing capabilities. The aforementioned processor can be a general-purpose processor, including a Central Processing Unit (CPU), a Network Processor (NP), etc.; it can also be a Digital Signal Processor (DSP), an Application Specific Integrated Circuit (ASIC), a Field Programmable Gate Array (FPGA), or other programmable logic devices, discrete gate or transistor logic devices, or discrete hardware components. It can implement or execute the methods, steps, and logic block diagrams disclosed in the embodiments of this invention. The general-purpose processor can be a microprocessor or any conventional processor.

[0212] It is understood that the structures illustrated in the embodiments of this application do not constitute a specific limitation on the computer device 1500. In other embodiments of this application, the computer device 1500 may include more or fewer components than illustrated, or combine some components, or split some components, or have different component arrangements. The illustrated components may be implemented in hardware, software, or a combination of software and hardware.

[0213] This disclosure also provides a computer-readable storage medium storing a computer program. When executed by a processor, the computer program performs the steps of the mountainous adaptive distributed scatterer interferometry method described in the above-described method embodiments. The storage medium can be a volatile or non-volatile computer-readable storage medium.

[0214] This disclosure also provides a computer program product carrying program code. The program code includes instructions that can be used to execute the steps of the mountain adaptive distributed scatterer interferometry method described in the above method embodiments. For details, please refer to the above method embodiments, which will not be repeated here.

[0215] The aforementioned computer program product can be implemented through hardware, software, or a combination thereof. In one optional embodiment, the computer program product is specifically embodied in a computer storage medium; in another optional embodiment, the computer program product is specifically embodied in a software product, such as a software development kit (SDK), etc.

[0216] Those skilled in the art will clearly understand that, for the sake of convenience and brevity, the specific working processes of the systems and devices described above can be referred to the corresponding processes in the foregoing method embodiments, and will not be repeated here. In the several embodiments provided in this disclosure, it should be understood that the disclosed systems and methods can be implemented in other ways. The device embodiments described above are merely illustrative. For example, the division of units is only a logical functional division; in actual implementation, there may be other division methods. Furthermore, multiple units or components may be combined or integrated into another system, or some features may be ignored or not executed. Another point is that the displayed or discussed mutual coupling or direct coupling or communication connection may be through some communication interfaces; the indirect coupling or communication connection of devices or units may be electrical, mechanical, or other forms.

[0217] The units described as separate components may or may not be physically separate. The components shown as units may or may not be physical units; that is, they may be located in one place or distributed across multiple network units. Some or all of the units can be selected to achieve the purpose of this embodiment according to actual needs.

[0218] In addition, the functional units in the various embodiments of this disclosure can be integrated into one processing unit, or each unit can exist physically separately, or two or more units can be integrated into one unit.

[0219] If the aforementioned functions are implemented as software functional units and sold or used as independent products, they can be stored in a processor-executable, non-volatile, computer-readable storage medium. Based on this understanding, the technical solution of this disclosure, in essence, or the part that contributes to the prior art, or a portion of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute all or part of the steps of the methods described in the various embodiments of this disclosure. The aforementioned storage medium includes various media capable of storing program code, such as USB flash drives, portable hard drives, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical disks.

[0220] Finally, it should be noted that the above-described embodiments are merely specific implementations of this disclosure, used to illustrate the technical solutions of this disclosure, and not to limit them. The protection scope of this disclosure is not limited thereto. Although this disclosure has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that any person skilled in the art can still modify or easily conceive of changes to the technical solutions described in the foregoing embodiments, or make equivalent substitutions for some of the technical features within the scope of the technology disclosed in this disclosure. Such modifications, changes, or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of this disclosure, and should all be covered within the protection scope of this disclosure.

Claims

1. A mountainous terrain adaptive distributed scatterer interferometric method, characterized in that, include: The radar image data and digital elevation model data of the target monitoring mountain area are acquired; wherein, the radar image data includes a single-polarization synthetic aperture radar image sequence and a dual-polarization synthetic aperture radar image sequence, the dual-polarization synthetic aperture radar image sequence includes multiple dual-polarization SAR images, and each dual-polarization SAR image includes multiple pixels. For the dual-polarization synthetic aperture radar image sequence, the radar vegetation index of each pixel in each dual-polarization SAR image is calculated, and based on the calculation results, the pixels in each dual-polarization SAR image are divided into multiple vegetation coverage levels to obtain a vegetation coverage classification atlas corresponding to the dual-polarization synthetic aperture radar image sequence. Based on the vegetation coverage grading atlas, statistical homogeneous pixel identification is performed on each pixel in each dual-polarization SAR image under the spatial constraint of vegetation coverage level, and a set of distributed scatterer candidate points corresponding to the dual-polarization synthetic aperture radar image sequence is determined based on the identification results; and, based on the vegetation coverage grading atlas, an adaptive coherence threshold is used to filter each distributed scatterer candidate point in the set of distributed scatterer candidate points to obtain multiple distributed scatterers. Based on the multiple distributed scatterers, the single-polarization synthetic aperture radar image sequence, and digital elevation model data, the long-term surface deformation results of the target monitoring mountain area are determined.

2. The method of claim 1, wherein, The calculation of radar vegetation index for each pixel in each dual-polarization SAR image and the classification of pixels in each dual-polarization SAR image into multiple vegetation coverage levels based on the calculation results include: For each dual-polarization SAR image, the vertical polarization backscattering coefficient and the horizontal polarization backscattering coefficient of each pixel in the dual-polarization SAR image are extracted respectively. For each pixel, the radar vegetation index of the pixel is calculated based on the vertical polarization backscattering coefficient and the horizontal polarization backscattering coefficient corresponding to the pixel. Based on multiple preset radar vegetation index threshold ranges, pixels whose radar vegetation index values ​​fall within the same threshold range are classified into the same vegetation coverage level. The vegetation coverage level of each pixel in each dual-polarization SAR image is labeled in the radar coordinate system to generate a vegetation coverage classification map corresponding to each dual-polarization SAR image; and the vegetation coverage classification maps corresponding to all dual-polarization SAR images are combined to obtain the vegetation coverage classification map set.

3. The method of claim 1, wherein, The step of identifying statistically homogeneous pixels in each pixel of each dual-polarization SAR image under spatial constraints of vegetation cover level, and determining the distributed scatterer candidate point set corresponding to the dual-polarization synthetic aperture radar image sequence based on the identification results, includes: A water body mask is constructed based on a preset backscattering coefficient threshold, and water body pixels in each dual-polarization SAR image are removed using the water body mask. For each non-water body pixel, a neighborhood window of a preset size is set with the non-water body pixel as the center, and neighborhood pixels belonging to the same vegetation coverage level as the non-water body pixel are filtered within the neighborhood window. For each non-water body pixel, a statistical homogeneity test is performed on the selected neighboring pixels corresponding to the non-water body pixel, and the neighboring pixels that pass the test and are spatially connected to the non-water body pixel are determined as homogeneous pixels of the non-water body pixel. The number of homogeneous pixels corresponding to each non-water body pixel is counted, and non-water body pixels with a number of homogeneous pixels greater than a preset homogeneous pixel threshold are identified as candidate points for distributed scatterers. The candidate distributed scatterers identified in all dual-polarization SAR images are combined according to the image time sequence to obtain the candidate distributed scatterer set.

4. The method of claim 3, wherein, The step of performing a statistical homogeneity test on the selected neighboring pixels corresponding to the non-water body pixels, and determining the neighboring pixels that pass the test and are spatially connected to the non-water body pixels as homogeneous pixels of the non-water body pixels, includes: For each selected neighboring pixel corresponding to the non-water body pixel, the backscattering coefficients of the neighboring pixel and the non-water body pixel are extracted in all phases of the dual-polarization synthetic aperture radar image sequence, and the neighborhood temporal amplitude sequence of the neighboring pixel and the center temporal amplitude sequence of the non-water body pixel are constructed respectively. A two-sample position test is performed on the neighborhood temporal amplitude sequence of the neighboring pixels and the center temporal amplitude sequence of the non-water body pixels to determine whether the two sequences come from the same continuous distribution. If the test statistic is less than the preset threshold, the neighboring pixel is determined to have passed the statistical homogeneity test; otherwise, it is determined to have failed. Among the neighboring pixels that pass the test, the neighboring pixels that are directly adjacent to the non-water body pixel, and the neighboring pixels that are indirectly connected to the non-water body pixel through other neighboring pixels that pass the test, are identified as homogeneous pixels that are spatially connected to the non-water body pixel.

5. The method of claim 1, wherein, Before performing adaptive coherence threshold screening on each candidate point of the distributed scatterer in the candidate point set based on the vegetation cover grading atlas, the process includes: Based on the single-polarization synthetic aperture radar image sequence and the digital elevation model data, an initial differential interferogram sequence is generated; The initial differential interferogram sequence is subjected to a first adaptive filtering process to obtain a differential interferogram sequence for candidate point screening, which is used to calculate the coherence coefficient of each candidate point of the distributed scatterer; The generation of the initial differential interferogram sequence includes: Based on a preset spatiotemporal baseline threshold, image pairs that satisfy the preset spatiotemporal baseline threshold are selected from the single-polarization synthetic aperture radar image sequence to form interferometric pairs. Interferometric processing is performed on each interferometric pair to obtain an interferogram. The terrain phase in the interferogram is removed using the digital elevation model data to obtain the initial differential interferogram sequence.

6. The method according to claim 5, characterized in that, The step of adaptively filtering candidate points of distributed scatterers based on the vegetation cover grading atlas and the candidate points of each distributed scatterer in the candidate point set includes: For each candidate distributed scatterer, based on the homogeneous pixel set of the candidate distributed scatterer, the coherence coefficient of the candidate distributed scatterer in each differential interferogram in the differential interferogram sequence for candidate point screening is calculated, and the arithmetic mean of all coherence coefficients is performed to obtain the temporal average coherence of the candidate distributed scatterer. Based on the vegetation cover classification atlas, the temporal average coherence of all candidate points of distributed scatterers within each vegetation cover level is statistically analyzed, and the regional spatial mean of each vegetation cover level is calculated respectively. Statistically calculate the temporal average coherence of all candidate points of distributed scatterers and then calculate the global spatial mean. For each vegetation cover level, the smaller of the regional spatial mean and the global spatial mean of the vegetation cover level is determined as the coherence screening threshold of the vegetation cover level. Distributed scatterer candidate points whose temporal average coherence is greater than or equal to the coherence screening threshold corresponding to the vegetation coverage level are identified as distributed scatterers.

7. The method of claim 1, wherein, The determination of long-term surface deformation results of the target monitoring mountain area based on the multiple distributed scatterers, the single-polarization synthetic aperture radar image sequence, and digital elevation model data includes: Based on the single-polarization synthetic aperture radar image sequence and the digital elevation model data, an initial differential interferogram sequence is generated; Based on the plurality of distributed scatterers and the homogeneous pixel set of the plurality of distributed scatterers, the initial differential interferogram sequence is subjected to a second adaptive filtering process to obtain a differential interferogram sequence for deformation inversion; The deformation inversion sequence is phase-unwrapped to obtain an unwrapped phase map; and the unwrapped phase map is subjected to time-series deformation inversion processing to obtain the deformation rate and cumulative deformation of each distributed scatterer. Atmospheric phase correction and orbital error removal are performed on the deformation rate and cumulative deformation of each distributed scatterer obtained by inversion to obtain the corrected deformation information; The corrected deformation information is converted from the radar coordinate system to the geographic coordinate system to generate long-term surface deformation results for the target monitoring mountain area; wherein, the long-term surface deformation results include a surface deformation distribution map and a time-series deformation curve.

8. A mountain-adaptive distributed scatterer interferometry device, comprising: include: The data acquisition module is used to acquire radar image data and digital elevation model data of the mountainous area for target monitoring; wherein, the radar image data includes a single-polarization synthetic aperture radar image sequence and a dual-polarization synthetic aperture radar image sequence, the dual-polarization synthetic aperture radar image sequence includes multiple dual-polarization SAR images, and each dual-polarization SAR image includes multiple pixels. The pixel classification module is used to calculate the radar vegetation index of each pixel in each dual-polarization synthetic aperture radar image sequence, and to classify the pixels in each dual-polarization synthetic aperture radar image sequence into multiple vegetation coverage levels based on the calculation results, so as to obtain a vegetation coverage classification map set corresponding to the dual-polarization synthetic aperture radar image sequence. The scatterer identification module is used to identify statistically homogeneous pixels in each pixel of each dual-polarization SAR image under the spatial constraint of vegetation coverage level, based on the vegetation coverage grading atlas, and to determine a set of distributed scatterer candidate points corresponding to the dual-polarization synthetic aperture radar image sequence based on the identification results; and to perform adaptive coherence threshold filtering on each distributed scatterer candidate point in the set of distributed scatterer candidate points based on the vegetation coverage grading atlas to obtain multiple distributed scatterers. The result determination module is used to determine the long-term surface deformation results of the target monitoring mountain area based on the multiple distributed scatterers, the single-polarization synthetic aperture radar image sequence, and digital elevation model data.

9. A storage medium having stored thereon a computer program, characterized in that When the computer program is executed by a processor, it implements the method of any one of claims 1 to 7.

10. A computer device comprising a storage medium, a processor, and a computer program stored on the storage medium and executable on the processor, characterized in that, When the processor executes the computer program, it implements the method of any one of claims 1 to 7.

Images