CN108983233B - PS point combination selection method in GB-InSAR data processing - Google Patents

Abstract

The invention discloses a PS point combination selection method in GB-InSAR data processing, comprising the following steps: step S1: coherence coefficient calculation; step S2: high coherence pixel extraction; step S3: calculation of the absolute deviation value of the coherence coefficient of each pixel Step S4: Evaluation and rejection of low-quality images; Step S5: PSC detection and selection; Step S6: Triangulation network establishment; Step S7: Calculation of spatial differential interference phase value; Step S8: absolute difference of differential interference phase value and Step S9: select PS points by threshold method; solve the problem that the real PS points cannot be accurately obtained in the past.

Description

Method of selecting PS point combination in GB-InSAR data processing

technical field

The invention relates to the field of GB-InSAR data processing, in particular to a method for selecting a combination of PS points in GB-InSAR data processing.

Background technique

The selection of PS points by the combination method combines the coherence coefficient threshold method, the amplitude dispersion threshold method and the differential phase standard deviation threshold method to identify the PS points. In the combination method PS selection, not only the acquisition of high signal-to-noise ratio pixels or even high-quality images is considered, but also the amplitude and phase information of each pixel, and the temporal stability and spatial continuity are considered. Choose the real PS point.

Firstly, with the coherence coefficient value calculated in the process of interference processing, the coherence coefficient threshold method is used to eliminate the seriously decoherent targets in the area, such as waters, and even the entire low-quality image is eliminated, and the image preprocessing of the PS selection work is completed. Then, the threshold method of amplitude dispersion index is used to select PS candidate points (ie PSC) to select stable points in time. Finally, a delaunay triangulation network is established for the PSC points. According to the standard deviation threshold of the differential phase value of each triangular network edge Select spatially continuous PS points that are less affected by atmospheric disturbance. .

SUMMARY OF THE INVENTION

In order to solve the problems existing in the prior art, the present invention provides a method for selecting a combination of PS points in GB-InSAR data processing, which solves the problem that the real PS points cannot be accurately obtained in the past.

The technical scheme adopted in the present invention is that the PS point combination selection method in the GB-InSAR data processing comprises the following steps:

Step S1: Obtain N-1 interference image pairs, and for any one of the available N-1 interference image pairs, obtain the coherence coefficient γ of any resolution unit in the image overlapping area according to the coherence coefficient calculation formula;

将大于平均相干系数阈值γ_T的平均相干系数

提取出来Step S2: According to the calculated coherence coefficient and the calculation formula of the average coherence coefficient, calculate the average coherence coefficient of the pixel (i, j) in the time series

will be greater than the average coherence coefficient threshold γ _T of the average coherence coefficient

extracted

The rest of the cells are deleted;

Step S3: Calculate the absolute deviation value Δγ _i,j according to the coherence coefficient and the average coherence coefficient of each pixel;

Step S4: Analyze the distribution of the absolute deviation value of the pixel coherence coefficient in each image, determine the threshold T according to the actual image quality, and calculate the number of pixels with △γ _i,j ≥T in the image accounting for the total image of the entire image. The ratio of the element is used as the standard for image quality evaluation, and the images with too large ratio value are eliminated to complete the preprocessing of PS detection and selection;

Step S5: extract the amplitudes of the pixels in the overlapping area of the processed M SAR images one by one to form an amplitude time series, and calculate the dispersion index D _A from the amplitude time series pixel by pixel;

Set a dispersion index threshold, when the calculated dispersion index D _A is less than the given threshold, the pixel is a PSC point, otherwise it is a non-PSC point;

Step S6: connect the adjacent PSC points according to the obtained PSC point according to "the delaunay triangulation TIN establishment method", and establish the triangulation;

Step S7: Calculate the spatial differential interference phase value of the triangular mesh edge PQ formed by two PSC points such as P point and Q point

计算时间序列上的差分干涉相位值的绝对差值及中误差m；Step S8: Interferometric phase value according to the spatial difference of PQ

Calculate the absolute difference and the median error m of the differential interference phase values on the time series;

Step S9: According to the calculated middle error result, a middle error threshold is set, and when the middle error of each side is less than the middle error threshold, the PSC point at the end point is used as the last PS point.

Preferably, step S6 includes the following steps:

Step S61: first construct a triangle according to the obtained PSC points to obtain a triangle;

Step S62: Establish a delaunay triangulation according to the acquired triangles.

Preferably, step S61 includes the following steps:

Step S611: according to the retrieval requirements of adjacent points, namely dividing the PSC points whose distance between two points is less than 30m into blocks;

Step S612: select a point A from several discrete points, and select the nearest point with a distance less than 30m near it as point B;

Step S613: Calculate ∠C _i according to the formula of the cosine theorem, and the formula of the cosine theorem is

Among them, a _i =BC _i ; b _i =AC _i ; c=AB; when ∠C=max{∠C _i }, then C is the third vertex of the triangle.

Preferably, step S62 includes the following steps:

Step S621: Extend the obtained triangle P ₁ P ₂ P ₃ outward, and expand the triangles of vertices P ₁ (x ₁ , y ₁ ), P ₂ (x ₂ , y ₂ ), P ₃ (x ₃ , y ₃ ) The P ₁ P ₂ side expands outward, and the equation of the P ₁ P ₂ line is:

F(x, y)=(y ₂ -y ₁ )(xx ₁ )-(x ₂ -x ₁ )(yy ₁ )=0

If the coordinates of the selected point P are (x, y), then when

When F(x, y) F(x ₃ , y ₃ )<0, P and P ₃ are on the opposite side of the straight line P ₁ P ₂ , and this point can be used as the vertex of the selected extension;

Step S622: According to any side of the triangle can only be the common side of two triangles at most, write down the number of times of expansion of any side, if the number of expansion exceeds 2, the expansion is invalid, otherwise the expansion is valid;

Step S623: After expanding the new edges of all the generated triangles, all discrete data points obtained are connected into an irregular triangular network.

Preferably, the calculation formula of the coherence coefficient in step S1 is:

In the formula, M(i, j) represents the main image complex number set, S(i, j) represents the secondary image complex number set, "*" represents the complex conjugate operator, which represents the coherence coefficient, m represents the number of interference image pairs, n represents the number of interference image pairs, and the value range of the coherence coefficient γ is [0, 1]; the coherence coefficient γ=0 indicates that the two images are completely incoherent; the coherence coefficient γ=1 indicates that the two images are completely coherent.

的计算公式为：Preferably, the average coherence coefficient of step S2

The calculation formula is:

In the formula, γ _i,j represents the coherence coefficient, and K represents the number of interference image pairs.

Preferably, the calculation formula of the absolute deviation value in step S3 is:

表示平均相干系数。where γ _{i, j} represent the coherence coefficient,

represents the average coherence coefficient.

Preferably, the calculation formula of the dispersion index D _A in step S5 is:

where D _A is the dispersion index, m _A is the mean amplitude, and σ _A is the standard deviation of the amplitude.

的空间差分干涉相位值

计算公式为：Preferably, in step S7

The spatial differential interference phase value of

The calculation formula is:

为P点相位，

为q点相位。In the formula,

is the phase at point P,

is the phase at point q.

Preferably, the absolute difference of the differential interference phase value in step S8 is the difference in time between the sides at the same position in space, and the middle error calculation formula is:

差分干涉相位值，

为差分干涉相位值的算术平均值。In the formula, m is the median error, ν is the correction number of the absolute difference of the differential interference phase value, n is the number of correction numbers, [νν] represents the sum of the squares of the correction numbers,

Differential interference phase value,

is the arithmetic mean of the differential interference phase values.

The beneficial effects of the PS point combination selection method in the GB-InSAR data processing of the present invention are as follows:

The invention adopts the PS point combination selection method in GB-InSAR data processing, which has the optimal detection ability and positioning method, and belongs to the non-contact measurement method, and does not need to manually set observation points.

Description of drawings

1 is a flowchart of a method for selecting a combination of PS points in GB-InSAR data processing of the present invention.

FIG. 2 is a flow chart of the triangulation network construction of the PS point combination selection method in the GB-InSAR data processing of the present invention.

FIG. 3 is a flow chart of triangle construction of the method for selecting a combination of PS points in the GB-InSAR data processing of the present invention.

FIG. 4 is a partial flow of the triangulation network construction of the PS point combination selection method in the GB-InSAR data processing of the present invention.

Detailed ways

The embodiments of the present invention will be described in detail below with reference to the accompanying drawings.

The specific embodiments of the present invention are described below to facilitate those skilled in the art to understand the present invention, but it should be clear that the present invention is not limited to the scope of the specific embodiments. For those skilled in the art, as long as various changes Such changes are obvious within the spirit and scope of the present invention as defined and determined by the appended claims, and all inventions and creations utilizing the inventive concept are within the scope of protection.

As shown in Figure 1, the method for selecting PS point combination in GB-InSAR data processing, step S1: acquiring N-1 interference image pairs, for any image pair in the available N-1 interference image pairs, according to the coherence The coefficient calculation formula is used to obtain the coherence coefficient γ of any resolution unit in the overlapping area of the image;

将大于平均相干系数的阈值像元γ_T提取出来，Step S2: According to the calculated coherence coefficient and the calculation formula of the average coherence coefficient, calculate the average coherence coefficient of the pixel (i, j) in the time series

Extract the threshold pixel _γT that is greater than the average coherence coefficient,

The rest of the cells are deleted;

Step S3: Calculate the absolute deviation value Δγ _i,j according to the coherence coefficient and the average coherence coefficient of each pixel;

Step S4: Analyze the distribution of the absolute deviation value of the pixel coherence coefficient in each image, determine the threshold T according to the actual image quality, and calculate the number of pixels with △γ _i,j ≥T in the image accounting for the total image of the entire image. The ratio of the element is used as the standard for image quality evaluation, and the images with too large ratio value are eliminated to complete the preprocessing of PS detection and selection;

Step S5: extract the amplitudes of the pixels in the overlapping area of the processed M SAR images one by one to form an amplitude time series, and calculate the dispersion index D _A from the amplitude time series pixel by pixel;

Set a dispersion index threshold, when the calculated dispersion index D _A is less than the given threshold, the pixel is a PSC point, otherwise it is a non-PSC point;

Step S6: connect the adjacent PSC points according to the obtained PSC point according to "the delaunay triangulation TIN establishment method", and establish the triangulation;

Step S7: Calculate the spatial differential interference phase value of the triangular mesh edge PQ formed by two PSC points such as P point and Q point

计算时间序列上的差分干涉相位值的绝对差值及中误差m；Step S8: Interferometric phase value according to the spatial difference of PQ

Calculate the absolute difference and the median error m of the differential interference phase values on the time series;

Step S9: According to the calculated middle error result, a middle error threshold is set, and when the middle error of each side is less than the middle error threshold, the PSC point at the end point is used as the last PS point.

When this embodiment is implemented, the method combines the coherence coefficient threshold method, the amplitude dispersion threshold method and the differential phase standard deviation threshold method to identify PS points. In the combination method PS selection, not only the acquisition of high signal-to-noise ratio pixels or even high-quality images is considered, but also the amplitude and phase information of each pixel, and the temporal stability and spatial continuity are considered. Choose the real PS point.

Firstly, with the coherence coefficient value calculated in the process of interference processing, the coherence coefficient threshold method is used to eliminate the seriously decoherent targets in the area, such as waters, and even the entire low-quality image is eliminated, and the image preprocessing of the PS selection work is completed. Then, the threshold method of amplitude dispersion index is used to select PS candidate points (ie PSC) to select stable points in time. Finally, a delaunay triangulation network is established for the PSC points, and the spatially continuous PS points with less atmospheric disturbance are selected according to the standard deviation threshold of the differential phase values of each triangular network edge.

As shown in Figure 2, the "network construction criteria" in step S6 includes the following steps:

Step S61: first construct a triangle according to the obtained PSC points to obtain a triangle;

Step S62: According to the acquired triangles, construct a triangular net to obtain a triangular net.

As shown in Figure 3, step S61 includes the following steps:

Step S611: According to the retrieval requirements of adjacent points, the PSC points whose distance between the two points is less than 30m are divided into blocks; the original data (obtained PSC data) is divided into blocks, so as to retrieve the points adjacent to the processed triangle, without having to retrieve All data and retrieval of all adjacent points require the distance between two points to be less than 30m.

Step S612: select a point A from several discrete points, and select the nearest point with a distance less than 30m near it as point B;

Step S613: Calculate ∠C _i according to the formula of the cosine theorem, and the formula of the cosine theorem is

Among them, a _i =BC _i ; b _i =AC _i ; c=AB; when ∠C=max{∠C _i }, then C is the third vertex of the triangle.

When this step is implemented, the establishment of the delaunay triangulation should be based on the conditions of the optimal triangle, that is, it should be ensured that each triangle is an acute triangle or the lengths of the three sides are approximately equal to avoid excessively large obtuse angles and too small acute angles. The following is the use of angle judgment method to establish TIN.

In this method, when two vertices (ie, one side) of a triangle are known, the law of cosines is used to calculate the size of the interior angle of the triangle whose third vertex is an angular vertex, and the point corresponding to the largest one is selected as the third vertex of the triangle. However, for the construction of PSC points, the spatial self-coherence of the atmospheric state of adjacent PS targets in the same image needs to be considered, that is, the distance value between adjacent PSC points needs to be considered. That is, in the case of satisfying the angle requirements, the side length requirements also need to be considered.

As shown in Figure 4, step S62 includes the following steps:

Step S621: Extend the obtained triangle P ₁ P ₂ P ₃ outward, and expand the triangles of vertices P ₁ (x ₁ , y ₁ ), P ₂ (x ₂ , y ₂ ), P ₃ (x ₃ , y ₃ ) The P ₁ P ₂ side expands outward, and the equation of the P ₁ P ₂ line is:

F(x, y)=(y ₂ -y ₁ )(xx ₁ )-(x ₂ -x ₁ )(yy ₁ )=0

If the coordinates of the selected point P are (x, y), then when

When F(x, y) F(x ₃ , y ₃ )<0, P and P ₃ are on the opposite side of the straight line P ₁ P ₂ , and this point can be used as the vertex of the selected extension;

Step S622: According to any side of the triangle can only be the common side of two triangles at most, write down the number of times of expansion of any side, if the number of expansion exceeds 2, the expansion is invalid, otherwise the expansion is valid;

Step S623: After expanding the new edges of all the generated triangles, all discrete data points obtained are connected into an irregular triangular network.

When this step is implemented, the first triangle is expanded outward to form a triangular net with all discrete points, and it is necessary to ensure that there are no repeated and intersecting triangles in the triangular net. The method is to expand the newly added two sides of each generated triangle in turn according to the principle of the largest angle, and to check whether it is repeated.

The calculation formula of the coherence coefficient in step S1 of this scheme is:

In the formula, M(i, j) represents the main image complex number set, S(i, j) represents the secondary image complex number set, "*" represents the complex conjugate operator, which represents the coherence coefficient, m represents the number of interference image pairs, n represents the number of interference image pairs, and the value range of the coherence coefficient γ is [0, 1]; the coherence coefficient γ=0 indicates that the two images are completely incoherent; the coherence coefficient γ=1 indicates that the two images are completely coherent.

的计算公式为：The average coherence coefficient of step S2 of this scheme

The calculation formula is:

In the formula, γ _i,j represents the coherence coefficient, and K represents the number of interference image pairs.

The calculation formula of the dispersion index D _A in step S5 of this scheme is:

where D _A is the dispersion index, m _A is the mean amplitude, and σ _A is the standard deviation of the amplitude.

的空间差分干涉相位值

计算公式为：Step S7 of this scheme

The spatial differential interference phase value of

The calculation formula is:

为P点相位，

为q点相位。In the formula,

is the phase at point P,

is the phase at point q.

的空间差分干涉相位值

计算公式为：Step S7 of this scheme

The spatial differential interference phase value of

The calculation formula is:

为P点相位，

为q点相位。In the formula,

is the phase at point P,

is the phase at point q.

The absolute difference of the differential interference phase value in step S8 of this scheme is the difference in time between the sides at the same position in space, and the middle error calculation formula is:

差分干涉相位值，

为差分干涉相位值的算术平均值。In the formula, m is the median error, ν is the correction number of the absolute difference of the differential interference phase value, n is the number of correction numbers, [νν] represents the sum of the squares of each correction number,

Differential interference phase value,

is the arithmetic mean of the differential interference phase values.

When this embodiment is implemented, the method combines the coherence coefficient threshold method, the amplitude dispersion threshold method and the differential phase standard deviation threshold method to identify PS points. In the combination method PS selection, not only the acquisition of high signal-to-noise ratio pixels or even high-quality images is considered, but also the amplitude and phase information of each pixel, and the temporal stability and spatial continuity are considered. Choose the real PS point.

Firstly, with the coherence coefficient value calculated in the process of interference processing, the coherence coefficient threshold method is used to eliminate the seriously decoherent targets in the area, such as waters, and even the entire low-quality image is eliminated, and the image preprocessing of the PS selection work is completed. Then, the threshold method of amplitude dispersion index is used to select PS candidate points (ie PSC) to select stable points in time. Finally, a delaunay triangulation network is established for the PSC points, and the spatially continuous PS points with less atmospheric disturbance are selected according to the standard deviation threshold of the differential phase values of each triangular network edge.

Claims (6)

1. PS point combination selection method in GB-InSAR data processing, is characterized in that, comprises the steps: Step S1: Obtain N-1 interference image pairs, and for any one of the available N-1 interference image pairs, obtain the coherence coefficient γ of any resolution unit in the image overlapping area according to the coherence coefficient calculation formula; Step S2: According to the calculated coherence coefficient and the calculation formula of the average coherence coefficient, calculate the average coherence coefficient of the pixel (i, j) in the time series

will be greater than the average coherence coefficient threshold γ _T of the average coherence coefficient

extracted

The rest of the cells are deleted; Step S3: Calculate the absolute deviation value Δγ _i,j according to the coherence coefficient and the average coherence coefficient of each pixel; Step S4: Analyze the distribution of the absolute deviation value of the pixel coherence coefficient in each image, determine the threshold T according to the actual image quality, and calculate the number of pixels with △γ _i,j ≥T in the image accounting for the total image of the entire image. The ratio of the element is used as the standard for image quality evaluation, and the images with too large ratio value are eliminated to complete the preprocessing of PS detection and selection; Step S5: Extract the amplitudes of the pixels in the overlapping area of the processed M SAR images one by one to form an amplitude time series, calculate the dispersion index D _A from the amplitude time series pixel by pixel, and set a dispersion index threshold, When the calculated dispersion index D _A is less than the given threshold, the pixel is a PSC point, otherwise it is a non-PSC point; Step S6: connect the adjacent PSC points according to the obtained PSC point according to "the delaunay triangulation TIN establishment method", and establish the triangulation; Step S7: Calculate the spatial differential interference phase value of the triangular mesh edge PQ formed by two PSC points such as P point and Q point

Step S8: Interferometric phase value according to the spatial difference of PQ

Calculate the absolute difference and the median error m of the differential interference phase values on the time series; Step S9: According to the calculated middle error result, a middle error threshold is set, and when the middle error of each side is less than the middle error threshold, the PSC point at its endpoint is used as the last PS point; The step S6 includes the following steps: Step S61: first construct a triangle according to the obtained PSC points to obtain a triangle; Step S62: establish a delaunay triangulation according to the acquired triangle; The step S61 includes the following steps: Step S611: According to the retrieval requirements of adjacent points, that is, dividing PSC points whose distance between two points is less than 30m into blocks; Step S612: select a point A from several discrete points, and select the nearest point with a distance less than 30m near it as point B; Step S613: Calculate ∠C _i according to the formula of the cosine theorem, and the formula of the cosine theorem is

Wherein, a _i =BC _i ; b _i =AC _i ; c=AB; when ∠C=max{∠C _i }, then C is the third vertex of the triangle; The step S62 includes the following steps: Step S621: Extend the obtained triangle P ₁ P ₂ P ₃ outward, and expand the triangles of vertices P ₁ (x ₁ , y ₁ ), P ₂ (x ₂ , y ₂ ), P ₃ (x ₃ , y ₃ ) The P ₁ P ₂ side expands outward, and the equation of the P ₁ P ₂ line is: F(x, y)=(y ₂ -y ₁ )(xx ₁ )-(x ₂ -x ₁ )(yy ₁ )=0 If the coordinates of the selected point P are (x, y), then when When F(x, y) F(x ₃ , y ₃ )<0, P and P ₃ are on opposite sides of the straight line P ₁ P ₂ , and this point is used as the vertex of the selected extension; Step S622: According to any side of the triangle can only be the common side of two triangles at most, write down the number of times of expansion of any side, if the number of expansion exceeds 2, the expansion is invalid, otherwise the expansion is valid; Step S623: After expanding the new edges of all the generated triangles, all discrete data points are obtained and connected into an irregular triangular network; The absolute difference of the differential interference phase value in the step S8 is the time difference between the sides at the same position in space, and the middle error calculation formula is:

In the formula, m is the median error, ν is the correction number of the absolute difference of the differential interference phase value, n is the number of correction numbers, [νν] represents the sum of the squares of each correction number,

Differential interference phase value,

is the arithmetic mean of the differential interference phase values. 2. the PS point combination selection method in GB-InSAR data processing according to claim 1, is characterized in that, the calculation formula of the coherence coefficient of described step S1 is:

In the formula, M(i, j) represents the main image complex number set, S(i, j) represents the secondary image complex number set, "*" represents the complex conjugate operator, which represents the coherence coefficient, m represents the number of interference image pairs, n represents the number of interference image pairs, and the value range of the coherence coefficient γ is [0, 1]; the coherence coefficient γ=0 indicates that the two images are completely incoherent; the coherence coefficient γ=1 indicates that the two images are completely coherent. 3. the PS point combination selection method in GB-InSAR data processing according to claim 1, is characterized in that, the average coherence coefficient of described step S2

The calculation formula is:

In the formula, γ _i,j represents the coherence coefficient, and K represents the number of interference image pairs. 4. the PS point combination selection method in GB-InSAR data processing according to claim 1, is characterized in that, the calculation formula of the absolute deviation value of described step S3 is:

where γ _{i, j} represent the coherence coefficient,

represents the average coherence coefficient. 5. the PS point combination selection method in the GB-InSAR data processing according to claim 1, is characterized in that, the calculation formula of the dispersion index D _A of described step S5 is:

where D _A is the dispersion index, m _A is the mean amplitude, and σ _A is the standard deviation of the amplitude. 6. the PS point combination selection method in the GB-InSAR data processing according to claim 1, is characterized in that, the step of described step S7

The spatial differential interference phase value of

The calculation formula is:

In the formula,

is the phase at point P,

is the phase at point q.

Images