InSAR baseline optimization method, device and computer-readable storage medium

Abstract

The present invention relates to the field of InSAR radar technology, and in particular to an InSAR baseline optimization method, device and computer-readable storage medium. An InSAR baseline optimization method, device and computer-readable storage medium are provided, wherein the method comprises: interferometrically superimposing first image data to obtain a plurality of interference image pairs; dividing a time period corresponding to second image data into a high vegetation coverage time period and a low vegetation coverage time period, and dividing the interference image pairs into a first interference image pair and a second interference image pair; calculating the coherence coefficient of the interference image pair to obtain a first optimized interference image pair and a second optimized interference image pair; performing SBAS-InSAR technology processing on the first optimized interference image pair and the second optimized interference image pair, dividing the high vegetation coverage time period and the low vegetation coverage time period by introducing optical remote sensing images, and introducing the coherence coefficient to correct the processing process of the first image data, so as to reduce the influence of vegetation coverage on radar signals.

Description

Technical Field

[0001] The present invention relates to the field of InSAR radar technology, and in particular to an InSAR baseline optimization method, device and computer-readable storage medium. Background Art

[0002] InSAR technology is an active remote sensing technology that uses a radar system to transmit microwave signals and record echo signals. Compared with traditional field measurements and optical remote sensing technologies, InSAR technology has the advantage of conducting large-scale, long-term sequence monitoring in harsh environments and weather conditions.

[0003] By comparing the phase differences of two or more radar beams in repeated orbits in the same area to infer surface elevation or surface deformation information, SAR images at multiple times can be processed to form multiple interferograms, obtain time series deformation information, and capture the temporal evolution of the surface. However, interferometric superposition of SAR images from different periods will result in a baseline network composed of interferometric image pairs of different qualities, which may contain some image pairs with poor coherence, directly affecting the accuracy of InSAR results, so the interferometric baseline network needs to be optimized.

[0004] At present, there are three main types of InSAR baseline optimization methods: the first type is based on expert experience, manually eliminating interferometric image pairs with poor coherence. This type of method is time-consuming and labor-intensive, and is not suitable for long-term InSAR processing with large amounts of data. In addition, this type of method relies on the prior knowledge of experts and has a certain degree of subjectivity. The second type of method is to impose threshold constraints on the temporal and spatial baselines. The accuracy of this type of method depends on the setting of the temporal and spatial baseline thresholds. When the threshold is set too small, the number of connected interferometric image pairs will be reduced, resulting in the loss of deformation information; when the threshold is set too large, temporal and spatial decoherence will occur, which will reduce the quality of the interferometric image pairs. Therefore, this type of method has great instability, and the setting of temporal and spatial baseline thresholds still relies more on past experience. The third type is a method that uses the average coherence coefficient as a threshold to screen interferometric image pairs, namely the average coherence coefficient threshold method. When the coherence of the interferometric image pairs shows obvious seasonal changes and the range is large, this method may cause the interferometric baseline to be disconnected, and subsequent inversion of surface deformation information cannot be performed.

[0005] At the same time, dense vegetation will cause radar signals to be reflected or scattered multiple times between leaves and branches, making it difficult for them to reach the surface. In addition, the moisture and other components in the vegetation will also absorb the radar signals to a certain extent, causing radar signals to attenuate. Therefore, the impact of vegetation coverage on radar signals also needs to be considered. Summary of the invention

[0006] (I) Purpose of the invention

[0007] The purpose of the present invention is to provide an InSAR baseline optimization method, device and computer-readable storage medium, which can reduce the impact of dense vegetation on radar signals.

[0008] (II) Technical solution

[0009] To solve the above problems, the first aspect of the present invention provides an InSAR baseline optimization method, comprising:

[0010] Step 100: interference stacking processing step, after fully connecting the first image data of the study area, interference stacking processing is performed to obtain a plurality of interference image pairs, wherein the first image data is a single-view complex image;

[0011] Step 200: a vegetation division step, according to the second image data of the study area in the same period as the first image data, the time period corresponding to the second image data is divided into a high vegetation coverage period and a low vegetation coverage period, wherein the interference image pair corresponding to the high vegetation coverage period is the first interference image pair, and the interference image pair corresponding to the low vegetation coverage period is the second interference image pair, and the second image data is an optical remote sensing image;

[0012] Step 300: a step of calculating a coherence coefficient, calculating the coherence coefficient of the interference image pair, and respectively calculating a first average coherence coefficient of the first interference image pair and a second average coherence coefficient of the second interference image pair;

[0013] Step 400: a step of eliminating low-coherence interference image pairs, eliminating the first interference image pairs whose coherence coefficient is lower than the first average coherence coefficient to obtain a first optimized interference image pair, and similarly, eliminating the second interference image pairs whose coherence coefficient is lower than the second average coherence coefficient to obtain a second optimized interference image pair;

[0014] Step 500: SBAS-InSAR technology processing step, performing SBAS-InSAR technology processing on the first optimized interferometric image pair and the second optimized interferometric image pair.

[0015] Furthermore, the step 100 in the above InSAR baseline optimization method further includes:

[0016] Step 110: a redundant data processing step, before the interference superposition processing, firstly setting the spatial baseline of the first image data to a first preset value, and setting the time baseline to a second preset value;

[0017] Step 120: Atmospheric correction step, in which an atmospheric model is used to correct the atmospheric delay phase during the image acquisition period during the interference superposition processing, and a corrected image is obtained.

[0018] Furthermore, step 200 in the above InSAR baseline optimization method is specifically as follows:

[0019] The vegetation coverage and average vegetation coverage were calculated using pixel dichotomy, where the vegetation coverage was calculated as follows:

[0020]

[0021] Among them, FVC is the vegetation coverage, S is the single pixel information; S soil is soil information; that is, remote sensing pixel information without vegetation coverage; S veg It is the remote sensing pixel information of the part covered by vegetation;

[0022] The step 200 further includes:

[0023] Step 210: Normalization processing step. After normalization processing, formula (3) is improved to:

[0024]

[0025] Among them, NDVI is the normalized difference vegetation index; NDVI soil There is no vegetation cover information; NDVI veg To include information about green vegetation cover;

[0026] The vegetation coverage and the average vegetation coverage are calculated monthly. The calculation formula of the average vegetation coverage is as follows:

[0027]

[0028] Among them, FVC avg is the average vegetation coverage; FVC i is the vegetation coverage in the i-th month.

[0029] So the vegetation coverage is greater than the average vegetation coverage FVC avg The time period is the high vegetation coverage period, which is lower than the average vegetation coverage FVC avg The time period is the period of low vegetation cover.

[0030] Furthermore, in step 300 of the above InSAR baseline optimization method,

[0031] The calculation formula of the first average coherence coefficient is as follows:

[0032]

[0033] Among them, C avg1 is the first average coherence coefficient; C1 i is the coherence coefficient of the first interferometric image pair in the i-th month;

[0034] Similarly, the calculation formula of the second average coherence coefficient is as follows:

[0035]

[0036] Among them, C avg2 is the second average coherence coefficient; C2 i is the coherence coefficient of the second interferometric image pair in the i-th month.

[0037] Furthermore, step 500 in the above InSAR baseline optimization method specifically includes:

[0038] Step 510: Obtain the deformation amount in each time period, obtain (N+1) SAR images of the same orbit but different time phases of the study area, and calculate the deformation amount of each time period. A and the second moment T B The SAR images are interferometrically superimposed to generate the i-th (i=1,2…,N+1) interferogram. The interferometric phase of the pixel with the azimuth coordinate x and the range coordinate y is:

[0039]

[0040] in, is the interference phase of the pixel with azimuth coordinate x and range coordinate y; and T B and T A The interference phase of the pixel with the azimuth coordinate x and the distance coordinate y at the moment; is the accumulated deformation in the radar line of sight direction, is the residual topographic phase in the differential interferogram, is the atmospheric delay phase, is the decoherent noise.

[0041] The average phase velocity is then calculated as:

[0042]

[0043] Among them, v i For the T i and T i-1 The average phase velocity between moments; is the phase at the i-th moment; is the phase at the i-1th moment; T i -T i-1 is the time interval between the i-th and i-1-th moments;

[0044] The phase value of the i-th interference pattern can be expressed as:

[0045]

[0046] Among them, Tp -T p-1 is the time interval between the pth and p-1th moments; v p is the average phase velocity at the pth moment; is the total phase change of the i-th interference pattern;

[0047] Therefore, the integral of each time period over the time interval between the master and slave images is represented by an M×N matrix:

[0048]

[0049] Where B is an M×N matrix; v is a vector containing the speed in each time period; is a vector containing the phase values ​​at each time point;

[0050] Then, the singular value decomposition method is used to calculate the generalized inverse matrix of the matrix B, and then the minimum norm solution of the velocity vector is solved. Finally, the deformation of each time period can be obtained by integrating the velocity between the time periods.

[0051] Furthermore, step 500 in the above InSAR baseline optimization method further includes:

[0052] Step 520: a shadow part and an overlapping part identification step, obtaining a local incident angle θ of the radar beam and an incident angle α of the radar line of sight, and determining a shadow part and an overlapping part in the SAR image according to the local incident angle θ of the radar beam and the incident angle α of the radar line of sight;

[0053] Step 530: Geometric distortion area identification step, according to the radar satellite parameters, based on the LSM algorithm, radar visibility analysis is performed on the study area, the geometric distortion area is identified, and mask processing is performed on the geometric distortion area.

[0054] Furthermore, the judgment rules of the shadow part and the overlap part in the above InSAR baseline optimization method are as follows:

[0055]

[0056] Where θ is the local incident angle of the radar beam; α is the incident angle of the radar line of sight.

[0057] A second aspect of the present invention provides an InSAR baseline optimization device, comprising:

[0058] An interference superposition processing module, which performs interference superposition processing on the first image data of the research area after full connection, to obtain a plurality of interference image pairs, wherein the first image data is a single-view complex image;

[0059] The vegetation division module is used to divide the time period corresponding to the second image data into a high vegetation coverage time period and a low vegetation coverage time period according to the second image data of the study area in the same time period as the first image data, wherein the interference image pair corresponding to the high vegetation coverage time period is the first interference image pair, and the interference image pair corresponding to the low vegetation coverage time period is the second interference image pair;

[0060] a coherence coefficient calculation module, used to calculate the coherence coefficient of the interference image pair, and respectively calculate a first average coherence coefficient of the first interference image pair and a second average coherence coefficient of the second interference image pair,

[0061] A low-coherence interference image pair elimination module is used to eliminate the first interference image pair whose coherence coefficient is lower than the first average coherence coefficient to obtain a first optimized interference image pair, and similarly, to eliminate the second interference image pair whose coherence coefficient is lower than the second average coherence coefficient to obtain a second optimized interference image pair;

[0062] The SBAS-InSAR technology processing module is used to perform SBAS-InSAR technology processing on the first optimized interferometric image pair and the second optimized interferometric image pair.

[0063] A third aspect of the present invention provides a computer-readable storage medium having computer-readable instructions stored thereon, and when the computer-readable instructions are executed by a processor, the steps of the InSAR baseline optimization method according to any one of claims 1 to 7 are implemented.

[0064] The above technical solution of the present invention has the following beneficial technical effects: the present invention takes into account the change of surface vegetation coverage, divides the interference image pairs into two categories for screening respectively, and its basic idea is to divide all interference image pairs into interference image pairs of high vegetation coverage period and interference image pairs of low vegetation coverage period according to vegetation coverage and combined imaging time of interference image pairs, and then use the average coherence coefficient of the two categories of image pairs as a threshold to screen the two categories of image pairs, and basically eliminate all interference images of high vegetation coverage period. Although the method without baseline optimization has more interference image pairs, it contains a large number of low coherence image pairs, which reduces the accuracy of SBAS-InSAR surface deformation inversion. BRIEF DESCRIPTION OF THE DRAWINGS

[0065] Figure 1 is a schematic diagram of the location of a study area according to an embodiment of the present invention;

[0066] Figure 2 is a flow chart according to an embodiment of the present invention;

[0067] Figure 3 is a baseline connection diagram according to an embodiment of the present invention;

[0068] Figure 4 is a baseline connection diagram optimized according to an embodiment of the present invention;

[0069] Figure 5 is a schematic diagram of radar visibility analysis according to an embodiment of the present invention;

[0070] Figure 6 is a schematic diagram of landslide hazard identification results according to an embodiment of the present invention;

[0071] Figure 7 is a schematic diagram of detailed information of a landslide creep body according to an embodiment of the present invention;

[0072] Figure 8 is a block diagram of a device according to an embodiment of the present invention. DETAILED DESCRIPTION

[0073] In order to make the purpose, technical solutions and advantages of the present invention clearer, the present invention is further described in detail below in conjunction with specific embodiments and with reference to the accompanying drawings. It should be understood that these descriptions are exemplary only and are not intended to limit the scope of the present invention. In addition, in the following description, the description of known structures and technologies is omitted to avoid unnecessary confusion of the concept of the present invention. In the description of the present invention, it should be noted that the terms "first", "second" and "third" are only used for descriptive purposes and cannot be understood as indicating or implying relative importance.

[0074] The embodiments of the present invention will be described below in conjunction with the accompanying drawings.

[0075] refer to Figure 1 In this example, the Yuanmou Dry-Hot Valley (25°21′~26°6′N, 101°35′~102°6′E) in Chuxiong Yi Autonomous Prefecture, Yunnan Province, China is taken as the study area. Figure 1As shown. The Yuanmou Dry-Hot Valley is located in the northern part of the Yunnan Plateau, the eastern part of the Hengduan Mountains, and the lower reaches of the Longchuan River, a first-level tributary of the Jinsha River. It is a typical high mountain canyon terrain, with a general terrain feature of high surroundings and low in the middle. The overall altitude of the area is 868m~2838m, with a relative height difference of 1970m. The Yuanmou Fault Zone runs through Yuanmou County from north to south, with a large number of fault planes and broken zones, forming basins, hills, terraces and other landforms. Because it is located in the hinterland, the basin landform with a large relative height difference closes the source of water and steam. Under the combined effect of the southwest monsoon and continental airflow, the climate of the Yuanmou Dry-Hot Valley is distinct in dry and wet seasons, with an average annual precipitation of about 624mm, mainly concentrated in June-October, accounting for about 90% of the annual precipitation, and the dry season precipitation is only 10% of the annual precipitation. The annual evaporation is about 3848mm, which is 6 times the precipitation. The water evaporation ratio is extremely unbalanced, and the contradiction between water and heat is prominent, forming a typical dry and hot climate. Affected by the dry and hot climate, the vegetation in Yuanmou dry and hot valley has a unique community appearance and flora composition. The appearance of the vegetation changes significantly with the alternation of dry and wet seasons. The vegetation is well developed in the rainy season and sparse in the dry season. The main trees are Yunnan pine (Pinusyunnanensis), Phyllanthus emblica, and Ziziphus mauritiana. The main shrubs are Acacia farnesiana, Vitex negundo, and Osteomelesanthyllidifolia. The main herbaceous species is Heteropogon contortus, forming a plant community that has adapted to drought conditions for a long time.

[0076] At the same time, the method shown in this embodiment is verified. The data used in this embodiment include: ① 74 scenes of C-band Sentinel-1A orbit-raising data (hereinafter referred to as Sentinel-1A data) released by the Copernicus program of the European Space Agency (ESA), single look complex (SLC) images in the interferometric wide scanning mode (IW), polarization mode is VV, spatial resolution is 5m×20m, time span is from January 2021 to May 2023, used for interferometric stacking processing to obtain surface deformation information of Yuanmou Dry Hot Valley, which can be downloaded for free at https: / / search.asf.alaska.edu / # / . ② Sentinel-2 optical remote sensing images (hereinafter referred to as Sentinel-2 data) released by the Copernicus program of the European Space Agency (ESA), with a spatial resolution of 10m, used to calculate and obtain FVC of the study area, which can be downloaded for free at https: / / scihub.copernicus.eu. ③ Copernicus Sentinel POD precise orbit determination ephemeris, used to improve satellite orbit accuracy, can be obtained online at https: / / scihub.copernicus.eu / gnss / # / home; ④ 30-meter resolution digital elevation model (DEM) provided by (Japan Aerospace Exploration Agency, JAXA) ALOS WORLD 3D, used to correct terrain phase effects in Sentinel-1A data processing, can be obtained online at https: / / www.eorc.jaxa.jp / ALOS / en / aw3d30 / data / index.htm. ⑤ GACOS data, used to correct the interferometric image pairs generated by Sentinel-1A data interferometric stacking processing, eliminate atmospheric delay errors, and improve InSAR accuracy, can be obtained free of charge at http: / / www.gacos.net / ⑥ Google satellite images, used to assist InSAR results in landslide identification, can be accessed at http: / / www.google.cn / intl / zh-CN / earth / . Detailed data introduction is shown in Table 1.

[0077] Table 1 Data introduction

[0078]

[0079] Reference below Figure 2 The InSAR baseline optimization method shown in this embodiment includes the following steps:

[0080] Step 100: Interference superposition processing step, after fully connecting the first image data of the study area, interference superposition processing is performed to obtain several groups of interference image pairs, the first image data is; in this embodiment, the first image data is Sentinel-1A data, and the time period is: January 2021 to May 2023.

[0081] Theoretically, all first image data are fully connected for interference superposition to generate as many interference image pairs as possible, which can provide more options for baseline optimization. However, as the interval between two SAR images in the first image data increases, the temporal and spatial decoherence of the interference image pairs formed by their combination will become more and more significant, resulting in many unusable low-coherence image pairs. To avoid generating redundant data, step 100 also includes:

[0082] Step 110: Redundant data processing step. Before interferometric stacking, the spatial baseline of the first image data is set to a first preset value, and the temporal baseline is set to a second preset value. In this embodiment, the spatial baseline of the Sentinel-1A data is set to 50%, and the temporal baseline is 108 days. A total of 423 interferometric image pairs are obtained. On this basis, the InSAR baseline optimization is performed. The baseline connection is as follows: Figure 3 shown.

[0083] The influence of the atmosphere on the electromagnetic wave signal causes the signal to be refracted during the propagation process, thereby causing the signal's propagation path and propagation direction to change. In the single-track InSAR mode, the atmospheric influences on the two images are very similar, and can basically cancel each other out after interference processing. However, in the repeated track mode, due to the change in the atmospheric state during the image acquisition time period, the atmosphere has a very large impact on the generated interference phase, resulting in an atmospheric delay phase. Therefore, the GACOS atmospheric model is used in the interference superposition processing to correct the atmospheric delay phase caused by the special natural climate conditions of the Yuanmou Dry Hot Valley. So step 100 also includes:

[0084] Step 120: Atmospheric correction step, in which an atmospheric model is used to correct the atmospheric delay phase during the image acquisition period during the interference superposition processing, and a corrected image is obtained.

[0085] In this embodiment, the GACOS atmospheric model (InSAR atmospheric correction system developed by the team of Professor Li Zhenhong of Newcastle University) is used. Based on the 0.1°×0.1° grid of the European Center for Medium-Term Weather Forecasts with high resolution, the iterative tropospheric decomposition model is used to separate the stratification and turbulence components from the tropospheric delay. The first image data is processed to obtain a 90m resolution zenith troposphere delay ZTD (Zenith Troposphere Delay) corrected image. The ZTD corrected image is defined as follows:

[0086]

[0087] Where k is the coordinate position; ZTD k is the total zenith delay of the vertical stratified component and the horizontal turbulent component; T represents the turbulence signal, x k is the site coordinate vector in the local geocentric coordinate system; L0 represents the layered component delay at sea level; β is the attenuation factor, which represents the rate of attenuation with increasing altitude, and ε k represents the remaining unmodeled residual, including the unmodeled stratification and turbulence signals; is the height scale, and the calculation formula is as follows:

[0088]

[0089] Among them, h k h is the height scale; min is the minimum elevation in the area; h max is the maximum elevation in the area.

[0090] The GACOS products of the same period as the SAR images were downloaded separately, and atmospheric delay phase correction was performed on the generated interferometric image pairs during the interferometric stacking process.

[0091] Step 200: Vegetation division step, according to the second image data of the study area in the same period as the first image data, the time period corresponding to the second image data is divided into a high vegetation coverage period and a low vegetation coverage period, wherein the interference image pair corresponding to the high vegetation coverage period is the first interference image pair, and the interference image pair corresponding to the low vegetation coverage period is the second interference image pair, and the second image data is an optical remote sensing image.

[0092] In this embodiment, the pixel dichotomy method is used to calculate the vegetation coverage, that is, the FVC (Fractional VegetaionCover) value, and the average vegetation coverage, that is, the FVC average value. Assuming that the surface consists of a pixel including a vegetation-covered part and a non-vegetation-covered part, the proportion of the vegetation-covered part in the pixel is calculated, that is, the calculation formula of the vegetation coverage is:

[0093]

[0094] Among them, S is the single pixel information; S soil is soil information; that is, remote sensing pixel information without vegetation coverage; S veg It is the remote sensing pixel information of the part covered by vegetation.

[0095] The Normalized Difference Vegetation Index (NDVI) is also a remote sensing index that reflects the surface vegetation conditions and is calculated by remote sensing sensors based on the spectral information of ground objects. The NDVI value of a pixel can be expressed as the information containing green vegetation coverage. veg and NDVI without vegetation cover soil The NDVI index can be used to analyze the two parameters S in the pixel binary model. soil and S veg Improvement is made, so step 200 also includes:

[0096] Step 210: Normalization processing step. After normalization processing, formula (3) is improved to:

[0097]

[0098] Among them, NDVI is the normalized difference vegetation index; NDVI soil There is no vegetation cover information; NDVI veg Contains information about green vegetation cover.

[0099] The FVC value and the average FVC value were calculated monthly. The formula for calculating the average FVC value is as follows:

[0100]

[0101] Among them, FVC avg is the mean value of FVC; FVC i is the FVC value at month i.

[0102] So the FVC value is greater than the mean FVC avg The time period is a period of high vegetation coverage, which is lower than the average FVC avg The time period is the period of low vegetation cover.

[0103] In this embodiment, July, August, September, October, November, and December of 2021 and July, August, October, and November of 2022 are classified as months with high vegetation coverage, and the remaining months are classified as months with low vegetation coverage.

[0104] Step 300: Calculating the coherence coefficient step, calculating the correlation coefficient of the interference image pair, and respectively calculating the first average coherence coefficient of the first interference image pair and the second average coherence coefficient of the second interference image pair.

[0105] The calculation formula of the first average coherence coefficient is as follows:

[0106]

[0107] Among them, C avg1 is the first average coherence coefficient; c1i is the coherence coefficient of the first interferometric image pair in the i-th month.

[0108] Similarly, the calculation formula of the second average coherence coefficient is as follows:

[0109]

[0110] Among them, C avg2 is the second average coherence coefficient; C2 i is the coherence coefficient of the second interferometric image pair in the i-th month.

[0111] The first average coherence coefficient C avg1 and the second average coherence coefficient C avg2 The value range of is [0,1], 0 means complete loss of coherence, 1 means complete coherence. The larger the coherence coefficient is, the higher the coherence of the interference image pair is, the higher the similarity between the two echoes is, and the interference phase difference can accurately reflect the distance difference between the echoes, thus obtaining more accurate surface deformation information.

[0112] Step 400: a step of eliminating low coherence interference image pairs, eliminating the first interference image pairs with a coherence coefficient lower than the first average coherence coefficient to obtain the first optimized interference image pairs, and similarly, eliminating the second interference image pairs with a coherence coefficient lower than the second average coherence coefficient to obtain the second optimized interference image pairs.

[0113] Therefore, in this embodiment, 64 interference image pairs of the first interference image pairs in the high vegetation coverage period are eliminated, and 74 interference image pairs of the second interference image pairs in the low vegetation coverage period are eliminated. The optimized baseline network is as follows: Figure 4 shown.

[0114] Step 500: SBAS-InSAR technology processing step, performing SBAS-InSAR technology processing on the first optimized interferometric image pair and the second optimized interferometric image pair.

[0115] The specific steps of step 500 are as follows:

[0116] Step 510: Obtain the deformation variables for each time period. Obtain (N+1) SAR images of the same orbit but different time phases in the study area. A and the second moment T B The SAR images are interferometrically superimposed to generate the i-th (i=1,2…,N+1) interferogram. The interferometric phase of the pixel with the azimuth coordinate x and the range coordinate y is:

[0117]

[0118] in, is the interference phase of the pixel with azimuth coordinate x and range coordinate y; and T B and T A The interference phase of the pixel with the azimuth coordinate x and the distance coordinate y at the moment; is the accumulated deformation in the radar line of sight direction, is the residual topographic phase in the differential interferogram, is the atmospheric delay phase, is the decoherent noise.

[0119] The average phase velocity is then calculated as:

[0120]

[0121] Among them, v i For the T i and T i-1 The average phase velocity between moments; is the phase at the i-th moment; is the phase at the i-1th moment; T i -T i-1 is the time interval between the i-th and i-1-th moments.

[0122] The phase value of the i-th interference pattern can be expressed as:

[0123]

[0124] Among them, T p -T p-1 is the time interval between the pth and p-1th moments, v p is the average phase velocity at the pth moment; is the total phase change of the i-th interference pattern;

[0125] Therefore, the integral of each time period over the time interval between the master and slave images is represented by an M×N matrix:

[0126]

[0127] Where B is an M×N matrix; v is a vector containing the speed in each time period; is a vector containing the phase values ​​at each time point;

[0128] Then, the singular value decomposition method is used to calculate the generalized inverse matrix of matrix B, and then the minimum norm solution of the velocity vector is solved. Finally, the deformation variable of each time period can be obtained by integrating the velocity between each time period.

[0129] Since SAR is a side-view imaging sensor, the radar beam is directed at the target object at a certain depression angle, and the image is recorded in the order of the reflection information of the target object in the distance upward. Even a slight change in elevation may cause a large-scale distortion in the image, which will affect the accuracy of the final result. Therefore, step 500 also includes:

[0130] Step 520: Shadow and overlapped part identification step, obtaining the local incident angle θ of the radar beam and the incident angle α of the radar line of sight, and judging the shadow and overlapped parts in the SAR image according to the local incident angle θ of the radar beam and the incident angle α of the radar line of sight, and the judging rule is as follows:

[0131]

[0132] Step 530: Geometric distortion area identification step, according to the radar satellite parameters, based on the LSM (Layover and Shadow Map) algorithm, radar visibility analysis is performed on the study area to identify the geometric distortion area and perform mask processing on the geometric distortion area.

[0133] refer to Figure 5 In this embodiment, the geometric distortion area is mainly concentrated in the northeast part of the Yuanmou dry hot valley, where the overlapping part accounts for 4.4% of the entire study area and the shadow part accounts for 0.2% of the entire study area.

[0134] After atmospheric correction, baseline optimization, geometric distortion area identification, and mask processing, the deformation rate of the radar line of sight (LOS) in the Yuanmou dry-hot valley from January 2021 to May 2023 was obtained, as shown in Figure 2. Figure 6 As shown in the figure, the overall deformation rate range is: -73.87mm / yr~127.35mm / yr. A positive deformation rate indicates that the deformation is in the direction of the satellite, that is, the surface is uplifted; a negative deformation rate indicates that the deformation is in the opposite direction of the satellite, that is, the surface is subsided. A total of 15 landslides and potential landslide areas were identified, named 1-15 respectively. The detailed deformation information of each landslide creep body is as follows Figure 7 shown.

[0135] The present invention also provides an InSAR baseline optimization device, comprising:

[0136] An interference superposition processing module, used for performing interference superposition processing on the first image data of the research area after full connection, to obtain a plurality of interference image pairs, wherein the first image data is a single-view complex image;

[0137] The vegetation division module is used to divide the time period corresponding to the second image data into a high vegetation coverage time period and a low vegetation coverage time period according to the second image data of the study area in the same time period as the first image data, wherein the interference image pair corresponding to the high vegetation coverage time period is the first interference image pair, and the interference image pair corresponding to the low vegetation coverage time period is the second interference image pair;

[0138] A coherence coefficient calculation module, used for calculating the coherence coefficient of the interference image pair, and respectively calculating a first average coherence coefficient of the first interference image pair and a second average coherence coefficient of the second interference image pair;

[0139] A low-coherence interference image pair elimination module is used to eliminate the first interference image pair whose coherence coefficient is lower than the first average coherence coefficient to obtain a first optimized interference image pair, and similarly, to eliminate the second interference image pair whose coherence coefficient is lower than the second average coherence coefficient to obtain a second optimized interference image pair;

[0140] The SBAS-InSAR technology processing module is used to perform SBAS-InSAR technology processing on the first optimized interferometric image pair and the second optimized interferometric image pair.

[0141] Preferably, the interference superposition processing module further includes:

[0142] A redundant data processing module, used for setting the spatial baseline of the first image data to a first preset value and setting the temporal baseline to a second preset value before interference superposition processing;

[0143] The atmospheric correction module is used to correct the atmospheric delay phase during the image acquisition period using the atmospheric model during interference superposition processing and obtain a corrected image.

[0144] Preferably, the vegetation division module is specifically used to calculate the vegetation coverage and the average vegetation coverage by using the pixel dichotomy method, wherein the calculation formula of the vegetation coverage is:

[0145]

[0146] Among them, FVC is the vegetation coverage, S is the single pixel information; S soil is soil information; that is, remote sensing pixel information without vegetation coverage; S veg It is the remote sensing pixel information of the part covered by vegetation;

[0147] The vegetation division module also includes:

[0148] The normalization processing module is used to improve formula (3) to:

[0149]

[0150] Among them, NDVI is the normalized difference vegetation index; NDVI soil There is no vegetation cover information; NDVI veg To include information about green vegetation cover;

[0151] The vegetation coverage and the average vegetation coverage are calculated monthly. The calculation formula of the average vegetation coverage is as follows:

[0152]

[0153] Among them, FVC avg is the average vegetation coverage; FVC i is the vegetation coverage in the i-th month.

[0154] So the vegetation coverage is greater than the average vegetation coverage FVC avg The time period is the high vegetation coverage period, which is lower than the average vegetation coverage FVC avg The time period is the period of low vegetation cover.

[0155] In the module for calculating the coherence coefficient,

[0156] The calculation formula of the first average coherence coefficient is as follows:

[0157]

[0158] Among them, C avg1 is the first average coherence coefficient; C1 i is the coherence coefficient of the first interferometric image pair in the i-th month;

[0159] Similarly, the calculation formula of the second average coherence coefficient is as follows:

[0160]

[0161] Among them, C avg2 is the second average coherence coefficient; C2 i is the coherence coefficient of the second interferometric image pair in the i-th month.

[0162] SBAS-InSAR technology processing modules include:

[0163] The module for obtaining the deformation variables in each time period is used to obtain (N+1) SAR images of the same orbit but different phases in the study area. A and the second moment T B The SAR images are interferometrically superimposed to generate the i-th (i=1,2…,N+1) interferogram. The interferometric phase of the pixel with the azimuth coordinate x and the range coordinate y is:

[0164]

[0165] in, is the interference phase of the pixel with azimuth coordinate x and range coordinate y; and T B and T A The interference phase of the pixel with the azimuth coordinate x and the distance coordinate y at the moment; is the accumulated deformation in the radar line of sight direction, is the residual topographic phase in the differential interferogram, is the atmospheric delay phase, is the decoherent noise.

[0166] The average phase velocity is then calculated as:

[0167]

[0168] Among them, v i For the T i and T i-1 The average phase velocity between moments; is the phase at the i-th moment; is the phase at the i-1th moment; T i -T i-1 is the time interval between the i-th and i-1-th moments.

[0169] The phase value of the i-th interference pattern can be expressed as:

[0170]

[0171] Among them, T p -T p-1 is the time interval between the pth and p-1th moments; v p is the average phase velocity at the pth moment; is the total phase change of the i-th interference pattern;

[0172] Therefore, the integral of each time period over the time interval between the master and slave images is represented by an M×N matrix:

[0173]

[0174] Where B is an M×N matrix; v is a vector containing the speed in each time period; is a vector containing the phase values ​​at each time point;

[0175] Then, the singular value decomposition method is used to calculate the generalized inverse matrix of the matrix B, and then the minimum norm solution of the velocity vector is solved. Finally, the deformation of each time period can be obtained by integrating the velocity between the time periods.

[0176] Preferably, the SBAS-InSAR technology processing module further includes:

[0177] A shadow part and overlapping part recognition module is used to obtain the local incident angle θ of the radar beam and the incident angle α of the radar line of sight, and judge the shadow part and overlapping part in the SAR image according to the local incident angle θ of the radar beam and the incident angle α of the radar line of sight;

[0178] The geometric distortion region identification module is used to perform radar visibility analysis on the study area based on radar satellite parameters and LSM algorithm, identify geometric distortion regions, and perform mask processing on the geometric distortion regions.

[0179] Preferably, the judgment rules of the shadow part and the overlapped part in the SBAS-InSAR technology processing module are as follows:

[0180]

[0181] Where θ is the local incident angle of the radar beam; α is the incident angle of the radar line of sight.

[0182] The present invention also provides a computer-readable storage medium, characterized in that:

[0183] The computer-readable storage medium stores computer-readable instructions, and when the computer-readable instructions are executed by a processor, the steps of any of the above-mentioned InSAR baseline optimization methods are implemented.

[0184] It should be understood that the above specific embodiments of the present invention are only used to illustrate or explain the principles of the present invention, and do not constitute a limitation of the present invention. Therefore, any modifications, equivalent substitutions, improvements, etc. made without departing from the spirit and scope of the present invention should be included in the protection scope of the present invention. In addition, the appended claims of the present invention are intended to cover all changes and modifications that fall within the scope and boundaries of the appended claims, or the equivalent forms of such scope and boundaries.

Claims

1. An InSAR baseline optimization method, characterized in that: include: Step 100: interference stacking processing step, after fully connecting the first image data of the study area, interference stacking processing is performed to obtain a plurality of interference image pairs, wherein the first image data is a single-view complex image; Step 200: a vegetation division step, according to the second image data of the study area in the same period as the first image data, the time period corresponding to the second image data is divided into a high vegetation coverage period and a low vegetation coverage period, wherein the interference image pair corresponding to the high vegetation coverage period is the first interference image pair, and the interference image pair corresponding to the low vegetation coverage period is the second interference image pair, and the second image data is an optical remote sensing image; The step 200 is specifically as follows: The vegetation coverage and average vegetation coverage were calculated using pixel dichotomy, where the vegetation coverage was calculated as follows: Among them, FVC is the vegetation coverage, S is the single pixel information; S soil is soil information; that is, remote sensing pixel information without vegetation cover; S veg It is the remote sensing pixel information of the part covered by vegetation; The step 200 further includes: Step 210: Normalization processing step. After normalization processing, formula (3) is improved to: Among them, NDVI is the normalized difference vegetation index; NDVI soil There is no vegetation cover information; NDVI veg To include information about green vegetation cover; The vegetation coverage and the average vegetation coverage are calculated monthly. The calculation formula of the average vegetation coverage is as follows: Among them, FVC avg is the average vegetation coverage; FVC i is the vegetation coverage in the i-th month; So the vegetation coverage is greater than the average vegetation coverage FVC avg The time period is the high vegetation coverage period, which is lower than the average vegetation coverage FVC avg The time period is the period of low vegetation cover; Step 300: a step of calculating a coherence coefficient, calculating the coherence coefficient of the interference image pair, and respectively calculating a first average coherence coefficient of the first interference image pair and a second average coherence coefficient of the second interference image pair; Step 400: a step of eliminating low-coherence interference image pairs, eliminating the first interference image pairs whose coherence coefficient is lower than the first average coherence coefficient to obtain a first optimized interference image pair, and similarly, eliminating the second interference image pairs whose coherence coefficient is lower than the second average coherence coefficient to obtain a second optimized interference image pair; Step 500: SBAS-InSAR technology processing step, performing SBAS-InSAR technology processing on the first optimized interferometric image pair and the second optimized interferometric image pair.

2. The InSAR baseline optimization method according to claim 1, characterized in that: The step 100 further includes: Step 110: a redundant data processing step, before the interference superposition processing, firstly setting the spatial baseline of the first image data to a first preset value, and setting the time baseline to a second preset value; Step 120: Atmospheric correction step, in which an atmospheric model is used to correct the atmospheric delay phase during the image acquisition period during the interference superposition processing, and a corrected image is obtained.

3. The InSAR baseline optimization method according to claim 1, characterized in that: In step 300, The calculation formula of the first average coherence coefficient is as follows: Among them, C avg1 is the first average coherence coefficient; C1 i is the coherence coefficient of the first interferometric image pair in the i-th month; Similarly, the calculation formula of the second average coherence coefficient is as follows: Among them, C avg2 is the second average coherence coefficient; C2 i is the coherence coefficient of the second interferometric image pair in the i-th month.

4. The InSAR baseline optimization method according to claim 1, characterized in that: The step 500 specifically includes: Step 510: Obtain the deformation amount in each time period, obtain (N+1) SAR images of the same orbit but different time phases of the study area, and calculate the deformation amount of each time period. A and the second moment T B The SAR images are interferometrically superimposed to generate the i-th (i=1,2…,N+1) interferogram. The interferometric phase of the pixel with the azimuth coordinate x and the range coordinate y is: in, is the interference phase of the pixel with azimuth coordinate x and range coordinate y; and T B and T A The interference phase of the pixel with the azimuth coordinate x and the distance coordinate y at the moment; is the accumulated deformation in the radar line of sight direction, is the residual topographic phase in the differential interferogram, is the atmospheric delay phase, is the decorrelation noise; The average phase velocity is then calculated as: Among them, v i For the T i and T i-1 The average phase velocity between moments; is the phase at the i-th moment; is the phase at the i-1th moment; T i -T i-1 is the time interval between the i-th and i-1-th moments; The phase value of the i-th interference pattern can be expressed as: Among them, T p -T p-1 is the time interval between the pth and p-1th moments; v p is the average phase velocity at the pth moment; is the total phase change of the i-th interference pattern; Therefore, the integral of each time period over the time interval between the master and slave images is represented by an M×N matrix: Where B is an M×N matrix; v is a vector containing the speed in each time period; is a vector containing the phase values ​​at each time point; Then, the singular value decomposition method is used to calculate the generalized inverse matrix of the matrix B, and then the minimum norm solution of the velocity vector is solved. Finally, the deformation of each time period can be obtained by integrating the velocity between the time periods.

5. The InSAR baseline optimization method according to claim 4, characterized in that: The step 500 further includes: Step 520: a shadow part and an overlapping part identification step, obtaining a local incident angle θ of the radar beam and an incident angle α of the radar line of sight, and determining a shadow part and an overlapping part in the SAR image according to the local incident angle θ of the radar beam and the incident angle α of the radar line of sight; Step 530: Geometric distortion area identification step, according to the radar satellite parameters, based on the LSM algorithm, radar visibility analysis is performed on the study area, the geometric distortion area is identified, and mask processing is performed on the geometric distortion area.

6. The InSAR baseline optimization method according to claim 5, characterized in that: The judgment rules of the shadow part and the overlapping part are as follows: Where θ is the local incident angle of the radar beam; α is the incident angle of the radar line of sight.

7. An InSAR baseline optimization device, characterized in that: include: An interference superposition processing module, which performs interference superposition processing on the first image data of the research area after full connection, to obtain a plurality of interference image pairs, wherein the first image data is a single-view complex image; The vegetation division module is used to divide the time period corresponding to the second image data into a high vegetation coverage time period and a low vegetation coverage time period according to the second image data of the study area in the same time period as the first image data, wherein the interference image pair corresponding to the high vegetation coverage time period is the first interference image pair, and the interference image pair corresponding to the low vegetation coverage time period is the second interference image pair; The vegetation division module is specifically as follows: The vegetation coverage and average vegetation coverage were calculated using pixel dichotomy, where the vegetation coverage was calculated as follows: Among them, FVC is the vegetation coverage, S is the single pixel information; S soil is soil information; that is, remote sensing pixel information without vegetation coverage; S veg It is the remote sensing pixel information of the part covered by vegetation; The vegetation division module also includes: Normalization processing module, after normalization processing, formula (3) is improved to: Among them, NDVI is the normalized difference vegetation index; NDVI soil There is no vegetation cover information; NDVI veg To include information about green vegetation cover; The vegetation coverage and the average vegetation coverage are calculated monthly. The calculation formula of the average vegetation coverage is as follows: Among them, FVC avg is the average vegetation coverage; FVC i is the vegetation coverage in the i-th month; So the vegetation coverage is greater than the average vegetation coverage FVC avg The time period is the high vegetation coverage period, which is lower than the average vegetation coverage FVC avg The time period is the period of low vegetation cover; A coherence coefficient calculation module, used for calculating the coherence coefficient of the interference image pair, and respectively calculating a first average coherence coefficient of the first interference image pair and a second average coherence coefficient of the second interference image pair; A low-coherence interference image pair elimination module is used to eliminate the first interference image pair whose coherence coefficient is lower than the first average coherence coefficient to obtain a first optimized interference image pair, and similarly, to eliminate the second interference image pair whose coherence coefficient is lower than the second average coherence coefficient to obtain a second optimized interference image pair; The SBAS-InSAR technology processing module is used to perform SBAS-InSAR technology processing on the first optimized interferometric image pair and the second optimized interferometric image pair.

8. A computer-readable storage medium, characterized in that: The computer-readable storage medium stores computer-readable instructions, and when the computer-readable instructions are executed by a processor, the steps of the InSAR baseline optimization method according to any one of claims 1 to 7 are implemented.

Images