CN108387899B - Ground control point automatic selection method in synthetic aperture radar interferometry - Google Patents

Ground control point automatic selection method in synthetic aperture radar interferometry Download PDF

Info

Publication number
CN108387899B
CN108387899B CN201810340738.0A CN201810340738A CN108387899B CN 108387899 B CN108387899 B CN 108387899B CN 201810340738 A CN201810340738 A CN 201810340738A CN 108387899 B CN108387899 B CN 108387899B
Authority
CN
China
Prior art keywords
gcp
interference
pixel
candidate
window
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201810340738.0A
Other languages
Chinese (zh)
Other versions
CN108387899A (en
Inventor
盛业华
周良辰
张棫
张卡
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Nanjing Normal University
Original Assignee
Nanjing Normal University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Nanjing Normal University filed Critical Nanjing Normal University
Priority to CN201810340738.0A priority Critical patent/CN108387899B/en
Publication of CN108387899A publication Critical patent/CN108387899A/en
Application granted granted Critical
Publication of CN108387899B publication Critical patent/CN108387899B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
    • G01S13/9021—SAR image post-processing techniques
    • G01S13/9023—SAR image post-processing techniques combined with interferometric techniques

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Physics & Mathematics (AREA)
  • Electromagnetism (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • General Physics & Mathematics (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

The invention discloses a ground control point automatic selection method in synthetic aperture radar interferometry, which comprises the following steps of randomly selecting two images from a multi-sequence SAR S L C image sequence to be subjected to interferometry to form an interference pair, obtaining an interference pattern through interference phase calculation and calculating a coherence coefficient pattern, partitioning the interference pattern according to the required GCP number, searching a region with stable or continuously and uniformly changed interference phases in each partition on the basis of solving the interference phase derivative of each pixel of the interference pair and counting the variance of the interference phase derivative in a local window of each pixel, and forming a candidate subset S of the GCP1(ii) a On the coherence coefficient map, pixels having a high coherence coefficient and a high coherence coefficient in a local window are selected to constitute a candidate subset S of GCP2(ii) a Solving a GCP candidate subset S1And S2Obtaining a preliminary GCP selection result by intersection; and optimizing the GCP initial selection result according to the quantity and distribution of the GCP to obtain the final GCP automatic selection result of the GCP. The method has important value for improving the automation degree and the precision of the radar interferometry of different modes.

Description

Ground control point automatic selection method in synthetic aperture radar interferometry
Technical Field
The invention belongs to the technology of mapping and remote sensing information application, and particularly relates to an automatic selection method of a ground control point in synthetic aperture radar interferometry.
Background
The method comprises the steps of utilizing an airborne or satellite-borne SAR antenna (or an antenna for multi-time sequence repeated observation) with interference imaging capability to transmit microwave signals to a target area, then receiving echoes reflected by the target, obtaining two or more Single view complex (S L C) images with certain view angle difference in the same target area and coherence, carrying out conjugate multiplication on the obtained complex images to obtain an interference diagram, or further carrying out differential processing (D-InSAR) on the interference diagram, comprehensively utilizing the height of a SAR sensor, the Radar wavelength, the beam view direction and the antenna base distance according to the phase value or the phase difference of the interference diagram, considering the geometrical relation between the earth and the topography, calculating the two or more path differences in the obtained interference diagram, and calculating the corresponding three-dimensional earth surface deformation of the SAR sensor, thereby calculating the three-dimensional earth surface deformation of the SAR sensor, and further calculating the corresponding micro ground surface deformation of the SAR and the micro ground surface deformation of the SAR.
The InSAR technology has been developed to date, and in addition to an interferometric processing (InSAR) mode for two complex images and a differential interferometric processing (D-InSAR) mode for three complex images, a differential interferometric measurement technology including a multiple time series SAR image (tens of images) such as a Permanent Scatterer (PS) radar interferometric measurement (PS-InSAR) mode, a Small baseline set (SBAS) radar interferometric measurement (SBAS-InSAR) mode, and the like, is developed, and aims to obtain long-time series surface deformation information.
For satellite-borne radar interferometry, whether InSAR, D-InSAR, PS-InSAR or SBAS-InSAR technology, in order to obtain high-precision terrain or surface deformation information and improve the precision of radar interferometry results, the influence caused by satellite orbit errors must be eliminated or weakened, the satellite orbit and phase offset must be corrected, orbit refinement and phase offset calculation are performed, and possible slope phases are eliminated.
The current mainstream commercial radar interferometry software, namely commercial GAMMA, SARScape, Earth View INSAR, open source Doris, ROI _ PAC and the like, adopts a manual interaction mode to select a plurality of GCPs from a single image or an interference pair in related processing links related to GCP selection. However, the manual interaction of the selection of the GCP in the SAR image or the interferogram not only has a large workload of human-computer interaction, but also depends on the experience of processing personnel in the point selection process, has strong subjectivity, and has extremely adverse effect on the processing result. Therefore, an automatic GCP selection method needs to be developed for each link related to GCP interactive selection in InSAR, D-InSAR and time sequence InSAR (including PS-InSAR and SBAS-InSAR) processing so as to improve the automation degree and precision of radar interferometry.
Disclosure of Invention
The invention aims to solve the defects in the prior art and provide an automatic selection method of ground control points in synthetic aperture radar interferometry.
In InSAR processing of different modes, the GCP is selected to meet the technical requirements that ① GCP is uniformly distributed in the overlapping area of multi-time sequence SAR S L C images (at least two images), ② GCP is selected at the position without residual terrain interference fringes or deformation interference fringes and is far away from a deformation area unless the deformation rate of the point is known, ③ GCP point is positioned at the position without Phase jump, if the GCP point is positioned on an isolated Phase and the quality of unwrapping is not ideal, the position can be a part of a Phase Ramp (Phase Ramp), and is not suitable to be selected as GCP, ④ in the process of multi-time sequence InSAR (such as SBAS-InSAR), the same GCP is difficult to be applied to all interference pairs (two single-view complex images with coherence), because different interference pairs have different coherence, therefore, the GCP is selected as much as possible in the time sequence InSAR to be applied to different interference pairs, and in general, at least 30 SAR repeated images in the same time sequence covering area are required to be selected in the same GCP.
The technical scheme is as follows: the invention discloses a method for automatically selecting ground control points in synthetic aperture radar interferometry, which is characterized by comprising the following steps of: the method sequentially comprises the following steps:
step one, two SAR S L C images C are selected from a multi-time sequence SAR S L C image sequence to be subjected to interferometry at will1And C2Forming interference pair, using window cross-correlation algorithm pair C1And C2Carrying out automatic registration and calculating to obtain an interference pattern; the specific process is as follows:
(2.1) two SAR S L C images C1And C2Using window cross-correlation algorithm pair C1And C2Carrying out automatic registration to form an interference pair;
(2.2) C in registration as per equations (2) and (3)1And C2Respectively calculating the interference phase of each pixel to obtain an interference pattern;
Figure GDA0001671732300000031
wherein A is1、A2Are respectively C1、C2The intensity value of (a) of (b),
Figure GDA0001671732300000032
are respectively C1、C2Q is an imaginary number, i.e.
Figure GDA0001671732300000033
Figure GDA0001671732300000034
Wherein the content of the first and second substances,
Figure GDA0001671732300000035
is C2The complex number of the conjugate of (a),
Figure GDA0001671732300000036
multiplying the complex conjugates of the two images;
step two: according to the number k of control points needed by different processing modes or links, relatively uniformly and automatically searching areas with smooth and continuous phases on the interferogram, using the areas as potential distribution areas of GCP, and generating a GCP candidate subset S1(ii) a This ensures that the required number of GCPs are approximately evenly distributed throughout the interferogram and fall into regions where the interference phase is stable or continuously evenly varying; the specific process is as follows:
(3.1) uniformly partitioning the interferogram obtained in the step one according to m × n rectangular ranges according to the required number k of control points, so that each GCP point only falls into one of the partitions, wherein m × n is approximately equal to 2k, and m is approximately equal to n;
(3.2) calculating the interference phase derivative of each pixel in the window (4) by using a window l × l (such as 31 × 31, 33 × 33, etc.) with a certain size one by one, namely the gradient of the interference phase of the adjacent pixels in the x and y directions, and then counting the interference phase derivative variance Z in the window (5)r,cThe variance image is given to a pixel p (r, c) to obtain an interference phase derivative variance image;
Figure GDA0001671732300000037
Figure GDA0001671732300000038
wherein
Figure GDA0001671732300000041
Respectively averaging the derivatives of the interference phase in the x direction and the y direction in the window; i. j is row and column control variables respectively when the interference phase derivative of each pixel is calculated one by one in a window, and the row and column control variables are integers which are sequentially increased by 1; l is an odd number and has a value range of 29-35;
obviously, Zr,cZero or a small value, indicating that the interference phase is continuous; when the variance value of the derivative of the interference phase is larger, the interference phase is discontinuous and changes violently.
(3.3) sequencing each pixel in the interference phase derivative variance image from small to large according to the numerical value of the interference phase derivative variance, recording by adopting a linear table, sequentially accessing the linear table until 30% of the total pixel number of the image, and then taking the interference phase derivative variance corresponding to the position as a segmentation threshold value tZBinarizing the interference phase derivative variance of each pixel in the interference phase derivative variance image according to equation (6), and assigning to each pixel:
Figure GDA0001671732300000042
here, selecting "30% of the total pixel number of the image" enables only about 30% of the pixels to be binarized to "1" during the subsequent binarization, thereby reducing the GCP candidate area by no more than 30% of the total pixels of the whole image or each block at most;
(3.4) counting the total number of pixels with the value of 1 according to the binarization result of each pixel for each block in the step (3.1), and if the total number is less than 30% of the total number of pixels of the block, not considering the selection of GCP in the block; otherwise, the block becomes a potential block of the GCP;
(3.5) marking a connected region by adopting an 8-neighborhood connected region marking method aiming at each GCP potential block in the step (3.4) to obtain a plurality of continuous regions with the value of 1 pixel; then, the area of each continuous region is counted, namely, the number of pixels of the continuous region is counted, then the continuous region with the largest area ratio is selected as a possible falling region of the GCP, and the GCP candidate subset S is formed by recording through a linear table1。
Step three: for registered interference pair C1And C2Calculating the coherence coefficient gamma of each pixel in the interference pair, marking the pixels with high coherence coefficient as possible points of GCP, and generating a GCP candidate subset S2(ii) a This enables the GCP point to fall on a high coherence pixel; the specific process is as follows:
(4.1) calculating the coherence coefficient gamma of each pixel in the interference pair in a certain window range (the window range is an odd number of square windows, and the width value range of the window range is 11-15, such as 11 × 11 and 13 × 13) according to the formula (6);
Figure GDA0001671732300000043
(4.2) making the coherent coefficient gamma of the interference pair larger than a set threshold t1(t1The value range is 0.8-1.0, such as 0.8), and the average value of the coherence coefficient gamma of each pixel in the local window (odd square window, width range is 7-9, such as 7 × 7, 9 × 9) is larger than the set threshold t2(t2Pixels with values ranging from 0.5 to 0.7, e.g., 0.6) are marked, and the high-coherence pixels are recorded by a linear table to form a GCP candidate subset S2，0.5≤t1≤1，0.5≤t2≤1；
Step four: using GCP candidate subset S1And GCP candidate subset S2The formula (7) is adopted to carry out the set intersection operation to obtain the pixels with uniform distribution, stable phase difference or continuous change and high coherence, and the pixels (c) are usedi,ri) As GCP candidate points:
S＝S1∩S2formula (7)
Step five: optimizing the GCP candidate point set S to make the number of GCPs equal to the number k of GCPs required or actually input, and finally obtaining the number of GCPs meeting the requirements and the distribution thereof, wherein the specific process comprises the following steps:
(5.1) if the number of GCP candidate points is less than k, subdividing each GCP potential block generated in step two (2.4) into subblocks of 2 × 2 (the division of 2 × 2 is simple and feasible for each subblock without considering the width of the subblock), searching the GCP potential subblocks in the same way, marking continuous regions in each GCP potential subblock in step two (2.5), and reconstructing a GCP candidate subset S1Then with the candidate subset S2Solving the intersection of the candidate points to obtain a new candidate point set S;
(5.2) if the number of GCP candidate points is more than k, the continuous region in each block or sub-block of the interference pattern with the smallest continuous region area is selected from the GCP candidate subset S1Directly delete and thenAnd candidate subset S2Solving the intersection of the candidate points to obtain a new candidate point set S;
(5.3) repeating the steps until the number of the candidate points in the GCP candidate point set S is equal to the required number k of the GCP, and obtaining the final GCP automatic point selection result.
The method has the advantages that on the basis of calculating the coherence coefficient and the interference phase, the method automatically searches the areas with high coherence coefficient and stable or uniform phase change from the interference pair only according to two SAR S L C images capable of forming the interference pair, uniformly configures the required number of GCPs in the areas, and meets the requirements of processing links related to the GCPs on the number of the GCPs and the point positions in different modes of InSAR in an automatic mode.
Drawings
FIG. 1 is a flow chart of the present invention;
FIG. 2 is a schematic diagram of two SAR S L C images of the same region at different times in the embodiment;
FIG. 3 is a diagram illustrating interference phase derivative variance uniformity and continuous region labeling results in an embodiment;
FIG. 4 is a diagram illustrating the GCP automatic selection result in the interference pair of the embodiment;
FIG. 2(a) is an SAR S L C image C in the embodiment1FIG. 2(b) is an example SAR S L C image C2。
Detailed Description
The technical solution of the present invention is described in detail below, but the scope of the present invention is not limited to the embodiments.
Example (b):
in this embodiment, the requirements of links such as land leveling effect, satellite orbit parameter correction and re-land leveling, slant-to-ground distance conversion, geocoding, and the like on the GCP are taken as examples when a multi-temporal SAR S L C image sequence (the present invention is also suitable for SAR S L C images acquired by other satellites) acquired by a COSMO SkyMed satellite in italy and covering a transverse mountain part area is subjected to multi-sequence SAR analysis (SBAS-SAR mode) to acquire the deformation of the ground surface of the area.
As shown in fig. 1, the automatic selection method of the ground control point in the synthetic aperture radar interferometry includes the following steps:
step one, image preprocessing and calculation of interference phase and coherence coefficient
Selecting two SAR S L C images C with close time or short time interval from a multi-time sequence COSMO SkyMed SAR S L C (or other SARS L C images) image sequence repeatedly covering the same region1And C2Forming an interference pair, as shown in fig. 2, under the support of commercial or open source InSAR processing software, preprocessing the original SAR image such as Multi-view (Multi L ooking) and filtering, accurately registering the two S L C images by using window cross-correlation algorithm and using approximate orbit parameters attached to the original image, and calculating the interference phase of each pixel according to the formula (2) and the formula (6) for the interference pair formed by the two SAR S L C images
Figure GDA0001671732300000061
And a coherence coefficient γ (the window size is set to 11 × 11), and an interference pattern and a coherence coefficient pattern are obtained.
Step two: screening the required number of GCP uniformly distributed areas by using an interference pattern, which comprises the following specific steps:
(1) according to the required (input parameter) number k of GCPs (for example, 17 GCPs are required when the InSAR processing of two SAR images is carried out to construct the region DEM in the example), the interference pair image is averagely divided into 42 blocks of 7 × 6 (the number of the blocks is about 2 times of the required number of the GCPs);
(2) for each pixel p (r, c) on the interference image full image, in the window range with the pixel as the center and the size of 33 × 33, the interference phase derivative of each pixel in the window is calculated one by one according to the formula (3), and the derivative variance Z of each pixel in the window range is counted according to the formula (4)r,cThe variance image is given to a pixel p (r, c) to obtain an interference phase derivative variance image;
(3) setting a linear table equal to the size of interference phase derivative variance image (line number × column number), sorting the numerical values of each pixel in the interference phase derivative variance image from small to large, filling the linear table in sequence, accessing the linear table in sequence until 30% of the linear table, and taking the corresponding interference phase derivative variance as image binary divisionCut threshold tZBinarizing the interference phase derivative variance image according to equation (5), taking the value of each pixel as '1' (small interference phase derivative variance) or '0' (large interference phase derivative variance), and setting a binary image T with the same size as the interference phase derivative variance imageB(as shown in fig. 3), recording the binarization result;
(4) counting the number of pixels with the median value of 1 in each block according to the binary image corresponding to each block range by using the blocking result of the interference image in the step two (1), if the number of the pixels with the median value of 1 in a certain block is less than 30 percent of the total number of the pixels of the block, ignoring the block, and otherwise, setting the block as a potential block of the GCP;
(5) aiming at each potential block of GCP, in the range of the block, an 8-neighborhood connected domain marking method is adopted to carry out interference phase derivative variance binary image TBMarking a connected region to obtain a plurality of continuous regions with the value of 1 pixel; counting the number of pixels in each continuous region, selecting the continuous region with the largest area ratio as a possible falling region of the GCP, wherein each pixel in the continuous region is a candidate pixel of the GCP, and using a linear table S1The row and column numbers (c, r) of these pixels are recorded to obtain the GCP candidate subset S1。
Step three: and (3) screening high coherence points according to the coherence coefficient graph obtained in the first step, wherein the process comprises the following steps:
(1) on the coherence coefficient map, a coherence coefficient threshold t is set for each pixel1(in this example t10.8) and sets the average coherence coefficient threshold t within the local window range of each pixel (window size 9 × 9 in this example)2(in this example t2＝0.6)；
(2) Traversing the coherent coefficient graph, and comparing the coherent coefficient with t1And the average coherence coefficient of each pixel in the pixel window is also greater than t2Is marked with a linear table S2Recording the row and column numbers (c, r) of all pixels satisfying both conditions to obtain the GCP candidate subset S2。
Step four: for the GCP candidate subset S generated in step two1And step three generated GCP candidate subset S2Solving the intersection set by using the formula (7) to obtain a candidate point set S of the GCP;
step five: optimizing the GCP candidate point set S generated in step four to make the number of GCPs equal to the required number k of GCPs (in this example, k is 17), and obtaining the number of GCPs satisfying the requirement and their distribution, the process is:
(1) if the number of the GCP candidate points is less than k, each GCP potential block generated in the second step (4) is divided into subblocks of 2 × 2, the GCP potential subblocks are searched in the same way, continuous region marking is carried out on each GCP potential subblock according to the second step (5), and the GCP candidate subset S is reconstructed1And then with the candidate subset S generated in step three2And solving the intersection to obtain a GCP candidate point set S.
(2) If the number of GCP candidate points is more than k, the continuous region in each block or sub-block of the interference pattern with the smallest area is selected from the GCP candidate subset S1Directly deleting the candidate subset S generated in the step three2And solving the intersection to obtain a GCP candidate point set S.
(3) Repeating the above steps until the number of candidate points in the GCP candidate point set S is equal to the required number k of GCPs, and obtaining the final automatic point selection result of GCP (as shown in fig. 4, the center of the red cross corresponds to the position in the figure).
It can be seen from the implementation that, in the application of the radar interferometry of different modes, the Ground Control Point (GCP) with uniform distribution, continuous interference phase, high coherence coefficient and meeting the quantity requirement can be automatically selected by only using the information of the interference image to the Ground Control Point.

Claims (4)

1. A ground control point automatic selection method in synthetic aperture radar interferometry is characterized in that: the method comprises the following steps:
step one, two SAR S L C images C are selected from a multi-time sequence SAR S L C image sequence to be subjected to interferometry at will1And C2Forming interference pair, using window cross-correlation algorithm pair C1And C2Carrying out automatic registration and calculating to obtain an interference pattern;
step two: according to the GCP number k required by different processing modes or links, relatively uniformly and automatically searching for areas with smooth and continuous phases on the interferogram, using the areas as potential distribution areas of GCP, and generating a GCP candidate subset S1；
Step three: for registered interference pair C1And C2Calculating the coherence coefficient gamma of each pixel in the interference pair, marking the pixels with high coherence coefficient as possible points of GCP, and generating a GCP candidate subset S2；
Step four: using GCP candidate subset S1And GCP candidate subset S2Performing set intersection operation by adopting the formula (1) to obtain pixels which are uniformly distributed, have stable phase difference or continuously change and have high coherence, and taking the pixels p (c, r) as candidate points of GCP:
S＝S1∩S2formula (1)
Step five: optimizing the GCP candidate point set S to ensure that the number of GCPs is equal to the number k of GCPs required or actually input, and finally obtaining the number of GCPs meeting the requirements and the distribution thereof;
the detailed method of the second step comprises the following steps:
(3.1) uniformly partitioning the interferogram obtained in the step one according to m × n rectangular ranges according to the required GCP number k, so that each GCP point only falls into one of the partitions, wherein m × n is approximately equal to 2k, and m is approximately equal to n;
(3.2) on the whole image of the interference image, for each pixel p (c, r), where c and r are the row number and column number of the pixel, respectively, calculating the interference phase derivative of each pixel in a window (4) by adopting a window l × l with a certain size one by one, namely the gradient of the interference phase of adjacent pixels in the x and y directions, and then counting the variance Z of the interference phase derivative in the window according to the formula (5)c,rThe variance image of the interference phase derivative is obtained by assigning the variance image to a pixel p (c, r);
Figure FDA0002452862690000011
Figure FDA0002452862690000012
wherein
Figure FDA0002452862690000021
Respectively averaging the derivatives of the interference phase in the x direction and the y direction in the window; i. j is row and column control variables respectively when the interference phase derivative of each pixel is calculated one by one in a window, and the row and column control variables are integers which are sequentially increased by 1; l is an odd number and has a value range of 29-35;
(3.3) sequencing each pixel in the interference phase derivative variance image from small to large according to the numerical value of the interference phase derivative variance, recording by adopting a linear table, sequentially accessing the linear table until 30% of the total pixel number of the image, and then taking the interference phase derivative variance corresponding to the position as a segmentation threshold value tZBinarizing the interference phase derivative variance of each pixel in the interference phase derivative variance image according to equation (6), and assigning to each pixel:
Figure FDA0002452862690000022
(3.4) counting the total number of pixels with the value of 1 according to the binarization result of each pixel for each block in the step (3.1), and if the total number is less than 30% of the total number of pixels of the block, not considering the selection of GCP in the block; otherwise, the block becomes a potential block of the GCP;
(3.5) marking a connected region by adopting an 8-neighborhood connected region marking method aiming at each GCP potential block in the step (3.4) to obtain a plurality of continuous regions with the value of 1 pixel; then, the area of each continuous region is counted, namely, the number of pixels of the continuous region is counted, then the continuous region with the largest area ratio is selected as a possible falling region of the GCP, and the GCP candidate subset S is formed by recording through a linear table1。
2. The method of claim 1, wherein the ground control point is selected automatically in the interferometry of synthetic aperture radar, the method comprising: the specific process of the step one is as follows:
(2.1) two SAR S L C images C1And C2Using window cross-correlation algorithm pair C1And C2Carrying out automatic registration to form an interference pair;
(2.2) C in registration as per equations (2) and (3)1And C2Respectively calculating the interference phase of each pixel to obtain an interference pattern;
Figure FDA0002452862690000023
wherein A is1、A2Are respectively C1、C2The intensity value of (a) of (b),
Figure FDA0002452862690000024
are respectively C1、C2Q is an imaginary number, i.e.
Figure FDA0002452862690000025
Figure FDA0002452862690000026
Wherein the content of the first and second substances,
Figure FDA0002452862690000027
is C2The complex number of the conjugate of (a),
Figure FDA0002452862690000028
is the complex conjugate multiplication of the two images.
3. The method of claim 1, wherein the ground control point is selected automatically in the interferometry of synthetic aperture radar, the method comprising: the specific process of the third step is as follows:
(4.1) calculating the coherence coefficient gamma of each pixel in the interference pair in a certain window range according to the formula (7);
Figure FDA0002452862690000031
the window range is odd number of square windows, and the width value range is 11-15;
(4.2) making the coherent coefficient gamma of the interference pair larger than a set threshold t1And the average value of the coherent coefficient gamma of each pixel in the local window is larger than the set threshold value t2Are marked, and these high-coherence pixels are recorded in a linear table to form a GCP candidate subset S2Here, t1The value range is 0.8-1.0, t2The value range is 0.5-0.7, the local window is a square window with odd number, and the width value range is 7-9.
4. The method of claim 1, wherein the ground control point is selected automatically in the interferometry of synthetic aperture radar, the method comprising: the concrete process of the step five is as follows:
(5.1) if the number of GCP candidate points is less than k, subdividing each GCP potential block into sub-blocks of 2 × 2, searching the GCP potential sub-blocks, performing continuous region marking on each GCP potential sub-block, and reconstructing a GCP candidate subset S1Then with the candidate subset S2Solving the intersection of the candidate points to obtain a new candidate point set S;
(5.2) if the number of GCP candidate points is more than k, the continuous region in each block or sub-block of the interference pattern with the smallest continuous region area is selected from the GCP candidate subset S1Direct deletion with the candidate subset S2Solving the intersection of the candidate points to obtain a new candidate point set S;
(5.3) repeating the steps until the number of the candidate points in the GCP candidate point set S is equal to the required number k of the GCP, and obtaining the final GCP automatic point selection result.
CN201810340738.0A 2018-04-17 2018-04-17 Ground control point automatic selection method in synthetic aperture radar interferometry Active CN108387899B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810340738.0A CN108387899B (en) 2018-04-17 2018-04-17 Ground control point automatic selection method in synthetic aperture radar interferometry

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810340738.0A CN108387899B (en) 2018-04-17 2018-04-17 Ground control point automatic selection method in synthetic aperture radar interferometry

Publications (2)

Publication Number Publication Date
CN108387899A CN108387899A (en) 2018-08-10
CN108387899B true CN108387899B (en) 2020-07-31

Family

ID=63064765

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810340738.0A Active CN108387899B (en) 2018-04-17 2018-04-17 Ground control point automatic selection method in synthetic aperture radar interferometry

Country Status (1)

Country Link
CN (1) CN108387899B (en)

Families Citing this family (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109444879A (en) * 2018-10-19 2019-03-08 西南交通大学 A kind of nearly tomography coseismic deformation extracting method of DInSAR
CN111059998B (en) * 2019-12-31 2020-11-13 中国地质大学（北京） High-resolution-based time sequence InSAR deformation monitoring method and system
CN111474544B (en) * 2020-03-04 2022-11-18 广东明源勘测设计有限公司 Landslide deformation monitoring and early warning method based on SAR data
CN111860158B (en) * 2020-06-15 2024-02-20 中国测绘科学研究院 Time sequence InSAR high coherence point extraction method fusing 1D-CNN and BiLSTM neural network
CN112130148B (en) * 2020-09-14 2021-12-28 北京大学 Land type-based DS self-adaptive selection method in InSAR time sequence analysis
CN112146622A (en) * 2020-10-23 2020-12-29 湖南航天智远科技有限公司 Method for monitoring geological settlement along power transmission line
CN112986994B (en) * 2021-02-06 2023-11-17 中国人民解放军国防科技大学 SAR chromatographic reference network rapid generation method
CN113311433B (en) * 2021-05-28 2022-08-02 北京航空航天大学 InSAR interferometric phase two-step unwrapping method combining quality map and minimum cost flow
CN114185048B (en) * 2022-02-15 2022-05-31 湖南吉赫信息科技有限公司 Method, system and storage medium for extracting landslide displacement vector by foundation InSAR
CN115143877B (en) * 2022-05-30 2023-04-25 中南大学 SAR offset tracking method, device, equipment and medium based on strong scattering amplitude suppression and outlier identification
CN115061136B (en) * 2022-06-08 2024-01-09 江苏省水利科学研究院 SAR image-based river and lake shoreline change point detection method and system
CN117092650A (en) * 2023-10-17 2023-11-21 月明星(北京)科技有限公司 Pixel interference method, device and equipment for synthetic aperture radar image
CN117609533B (en) * 2024-01-24 2024-04-05 航天宏图信息技术股份有限公司 Automatic screening method, device and equipment for synthetic aperture radar interference image pair

Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2001091650A (en) * 1999-09-20 2001-04-06 Mitsubishi Electric Corp Active ground control point device for synthetic aperture radar image precise geometrical correction
JP2008185375A (en) * 2007-01-29 2008-08-14 Mitsubishi Electric Corp 3d shape calculation device of sar image, and distortion correction device of sar image
CN102269813A (en) * 2011-06-23 2011-12-07 中国电子科技集团公司第三十八研究所 Interference processing technology of airborne non-vertical dual-antenna InSAR system
CN104361590A (en) * 2014-11-12 2015-02-18 河海大学 High-resolution remote sensing image registration method with control points distributed in adaptive manner
CN106203271A (en) * 2016-06-29 2016-12-07 南京大学 A kind of high ferro main line extracting method based on High Resolution SAR Images coherence
CN106338775A (en) * 2016-09-07 2017-01-18 民政部国家减灾中心（民政部卫星减灾应用中心） Building damage degree evaluation method based on interference synthetic aperture radar data
CN106443676A (en) * 2016-09-29 2017-02-22 中国科学院电子学研究所 Scarce control point space-borne synthetic aperture radar image ground positioning method
CN106780321A (en) * 2016-11-21 2017-05-31 中国测绘科学研究院 A kind of overall tight orientation of the satellite HR sensors images of CBERS 02 and correction joining method
CN107144823A (en) * 2017-06-16 2017-09-08 中国测绘科学研究院 A kind of interference calibrating method of airborne polarization interference synthetic aperture radar image

Patent Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2001091650A (en) * 1999-09-20 2001-04-06 Mitsubishi Electric Corp Active ground control point device for synthetic aperture radar image precise geometrical correction
JP2008185375A (en) * 2007-01-29 2008-08-14 Mitsubishi Electric Corp 3d shape calculation device of sar image, and distortion correction device of sar image
CN102269813A (en) * 2011-06-23 2011-12-07 中国电子科技集团公司第三十八研究所 Interference processing technology of airborne non-vertical dual-antenna InSAR system
CN104361590A (en) * 2014-11-12 2015-02-18 河海大学 High-resolution remote sensing image registration method with control points distributed in adaptive manner
CN106203271A (en) * 2016-06-29 2016-12-07 南京大学 A kind of high ferro main line extracting method based on High Resolution SAR Images coherence
CN106338775A (en) * 2016-09-07 2017-01-18 民政部国家减灾中心（民政部卫星减灾应用中心） Building damage degree evaluation method based on interference synthetic aperture radar data
CN106443676A (en) * 2016-09-29 2017-02-22 中国科学院电子学研究所 Scarce control point space-borne synthetic aperture radar image ground positioning method
CN106780321A (en) * 2016-11-21 2017-05-31 中国测绘科学研究院 A kind of overall tight orientation of the satellite HR sensors images of CBERS 02 and correction joining method
CN107144823A (en) * 2017-06-16 2017-09-08 中国测绘科学研究院 A kind of interference calibrating method of airborne polarization interference synthetic aperture radar image

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Automatic Detection and Positioning of Ground Control Points Using TerraSAR-X Multiaspect Acquisitions;Kamolratn Chureesampant,et al;《IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING》;20140131;p137-148 *
神东矿区开采沉陷干涉测量与时空变化规律分析;杨亚莉;《中国优秀硕士学位论文全文数据库 工程科技I辑》;20171215;全文 *

Also Published As

Publication number Publication date
CN108387899A (en) 2018-08-10

Similar Documents

Publication Publication Date Title
CN108387899B (en) Ground control point automatic selection method in synthetic aperture radar interferometry
Henriksen et al. Extracting accurate and precise topography from LROC narrow angle camera stereo observations
CN113624122B (en) Bridge deformation monitoring method fusing GNSS data and InSAR technology
CN104123464B (en) Method for inversion of ground feature high elevation and number of land subsidence through high resolution InSAR timing sequence analysis
CN106526590B (en) A kind of fusion multi-source SAR image industrial and mining area three-dimensional earth's surface deformation monitorings and calculation method
CN109388887B (en) Quantitative analysis method and system for ground settlement influence factors
Medjkane et al. High-resolution monitoring of complex coastal morphology changes: Cross-efficiency of SfM and TLS-based survey (Vaches-Noires cliffs, Normandy, France)
CN104122553B (en) Regional ground settlement monitoring method based on multiple track and long strip CTInSAR (coherent target synthetic aperture radar interferometry)
San et al. Digital elevation model (DEM) generation and accuracy assessment from ASTER stereo data
Zhang et al. Subsidence monitoring in coal area using time-series InSAR combining persistent scatterers and distributed scatterers
CN109100719B (en) Terrain map joint mapping method based on satellite-borne SAR (synthetic aperture radar) image and optical image
CN116047519B (en) Point selection method based on synthetic aperture radar interferometry technology
CN110703220B (en) Multi-baseline PolInSAR vegetation parameter inversion method considering time decoherence factors
CN112711021B (en) Multi-resolution InSAR (interferometric synthetic Aperture Radar) interactive interference time sequence analysis method
CN110109118B (en) Forest canopy biomass prediction method
CN108983231B (en) Interferometric video measuring method based on video synthetic aperture radar
CN113281749B (en) Timing sequence InSAR high coherence point selection method considering homogeneity
CN116299455A (en) Facility deformation analysis method based on PSInSAR and SqueseAR
CN108983239A (en) Spaceborne interferometric SAR digital elevation model reconstruction method
Nascetti et al. Fast terrain modelling for hydrogeological risk mapping and emergency management: the contribution of high-resolution satellite SAR imagery
Reinartz et al. Benchmarking and quality analysis of DEM generated from high and very high resolution optical stereo satellite data
CN103018728B (en) Laser radar real-time imaging and building characteristic extracting method
Re et al. CaSSIS-based stereo products for Mars after three years in orbit
Dong et al. Gapless-REMA100: A gapless 100-m reference elevation model of Antarctica with voids filled by multi-source DEMs
Momm et al. Methods for gully characterization in agricultural croplands using ground-based light detection and ranging

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant