CN112698328B - Phase unwrapping method and system for monitoring dam and landslide deformation GB-SAR - Google Patents

Phase unwrapping method and system for monitoring dam and landslide deformation GB-SAR Download PDF

Info

Publication number
CN112698328B
CN112698328B CN202011372417.2A CN202011372417A CN112698328B CN 112698328 B CN112698328 B CN 112698328B CN 202011372417 A CN202011372417 A CN 202011372417A CN 112698328 B CN112698328 B CN 112698328B
Authority
CN
China
Prior art keywords
phase
unwrapping
value
phase unwrapping
interference
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
CN202011372417.2A
Other languages
Chinese (zh)
Other versions
CN112698328A (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.)
Sichuan University
Original Assignee
Sichuan 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 Sichuan University filed Critical Sichuan University
Priority to CN202011372417.2A priority Critical patent/CN112698328B/en
Publication of CN112698328A publication Critical patent/CN112698328A/en
Application granted granted Critical
Publication of CN112698328B publication Critical patent/CN112698328B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO 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/00Systems 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/88Radar or analogous systems specially adapted for specific applications
    • G01S13/89Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO 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/00Systems 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/88Radar or analogous systems specially adapted for specific applications
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO 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/00Systems 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/88Radar or analogous systems specially adapted for specific applications
    • G01S13/89Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
    • G01S13/9021SAR image post-processing techniques
    • G01S13/9023SAR 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)
  • Testing Or Calibration Of Command Recording Devices (AREA)

Abstract

The invention provides a phase for monitoring the deformation GB-SAR of a dam and a landslideA method and a system for bit unwrapping, the method comprising: obtaining GB-SAR interference images in a deformation monitoring area, selecting N PS points with stable phases, and performing one-dimensional time phase unwrapping to obtain a time phase unwrapped interference phase value
Figure DDA0002807126050000011
Constructing a Delaunay triangulation network by using the selected PS points, and calculating a double-difference phase observation value of M edges in the triangulation network
Figure DDA0002807126050000012
Error equation construction based on indirect least square adjustment method
Figure DDA0002807126050000013
Converting error equation of M edges into matrix expression
Figure DDA0002807126050000014
Establishing a normal equation ATPA‑ATPL is 0, then Φ is (a)TPA)‑1ATPL. The invention provides a phase unwrapping method aiming at a GB-SAR monitoring system of a dam and a landslide in a high mountain canyon, so that the influence of atmospheric disturbance on unwrapping precision is reduced; integration is carried out in a mode of avoiding 'passing through' residual points, errors cannot be transmitted, and the unwinding precision is greatly improved; and a global solution method is adopted, so that the convergence speed is accelerated, and the unwrapping process is quicker.

Description

Phase unwrapping method and system for monitoring dam and landslide deformation GB-SAR
Technical Field
The invention relates to the field of GB-SAR interference monitoring, in particular to a phase unwrapping method and a phase unwrapping system for monitoring dam and landslide deformation GB-SAR.
Background
The Ground-based synthetic aperture radar (GB-SAR) is a technology which appears in the last ten years and is widely applied to the field of surface micro-deformation monitoring, can provide remote monitoring of thousands of kilometers, has continuous spatial coverage characteristics of deformation data, can obtain high spatial and temporal resolution and high-precision results in the process of monitoring all day, and is favorable for accurately judging the property, range and process of disaster deformation. The method is an effective supplement to satellite-borne SAR and traditional monitoring means.
Phase unwrapping refers to a process of restoring a phase main value or an interference phase difference value to a real phase value, and is the same as In-SAR technology, and phase unwrapping of an obtained interference image is a key and the most difficult link for SAR data processing. In the prior art, a branch tangent method proposed by Goldstein et al in 1988 is generally adopted for phase unwrapping, and the basic idea is to generate a shortest path by identifying positive and negative singular points of each pixel, and perform phase unwrapping along a shortest path integral. The problem that the unwinding island can be generated, so that some points can not be unwound normally, and the problem that the branch tangent line is not the shortest can be generated, so that the efficiency is low, and the effect is not ideal. In 1994, Pritt et al proposed a path-independent fast Fourier transform-based unweighted least squares method for phase unwrapping, and thereafter, many scholars followed the development of improved algorithms based thereon. The algorithm carries out smoothing treatment on the whole image through global fitting, so that the unwrapping phases of all points are continuous and the error is small. However, the method also has a problem that while the global fitting is used to improve the unwrapping precision of the whole image, the unwrapping result of each point is contaminated by the error therein, i.e., the unwrapping result contains errors, is an approximate value in nature, and is not accurate. Pepe et al put forward a minimum cost flow method in 2006 to perform phase unwrapping on a space-time combination Delaunay triangulation network under a multi-baseline-distance condition in a satellite-borne SAR. This method is not applicable to GB-SAR because the spatial baseline of GB-SAR is zero. Other phase unwrapping methods such as minimum network flow, Kalman filtering, etc. are also possible.
The methods are all proposed based on the basis of In-SAR instead of GB-SAR, so although the principle has common parts, the unwrapping effect of observation equipment is possibly far different due to different environmental influences on the observation equipment. For GB-SAR monitoring phases of dams and landslides In high mountain canyons, the GB-SAR monitoring phases are greatly influenced by the environment, the influence of climate change on signals is strong, and the spatial correlation of the atmosphere is very strong, so that the unwrapping effect of the phase unwrapping method directly using the In-SAR is not ideal.
Therefore, the problem to be solved in the field is to provide a method suitable for specific unwrapping of GB-SAR observation phases of dams and landslides.
Disclosure of Invention
Aiming at the defects in the prior art, the invention provides a phase unwrapping method for GB-SAR monitoring of dams and landslides in high mountain canyons, which has the advantages of strict description, high unwrapping precision and reliable quality, and avoids the accumulation of errors in the unwrapping process.
In order to achieve the purpose, the invention provides a phase unwrapping method for GB-SAR monitoring of dam and landslide deformation, which comprises the following steps:
obtaining GB-SAR interference images in a deformation monitoring area, selecting N PS points with stable phases, and performing one-dimensional time phase unwrapping to obtain a time phase unwrapped interference phase value
Figure GDA0002914540670000021
N={1,2,3…i,j…N};
Constructing a Delaunay triangulation network by using the selected PS points, and calculating a double-difference phase observation value of M edges in the triangulation network
Figure GDA0002914540670000022
Error equation construction based on indirect least square adjustment method
Figure GDA0002914540670000023
Wherein
Figure GDA0002914540670000024
For each selected PS pointUnwrapping the phase value;
converting error equation of M edges into matrix expression
Figure GDA0002914540670000025
Wherein V is an error matrix; a is a coefficient matrix; phi is a parameter matrix and corresponds to the unwrapping phase values of the N selected PS points; l is a parameter matrix and corresponds to the time phase unwrapping value of the M edges;
establishing a normal equation ATPA-ATPL is 0, then Φ is (a)TPA)-1ATPL, where P is the weight matrix.
In some preferred embodiments, the selecting N PS points with stable phase includes:
calculating a coherence coefficient of the interference image, and selecting pixels meeting high coherence;
calculating the absolute deviation value of the coherent coefficient of each pixel, setting a threshold value T, counting the proportion of pixels with absolute deviation values larger than the threshold value T to the total number of pixels in the image, and rejecting the image with overlarge proportion value;
calculating the amplitude dispersion index of each pixel time sequence, and determining candidate PS points according to the relation between the threshold T and the amplitude dispersion index;
establishing a Delaunay triangulation network according to the candidate PS points, and calculating the absolute difference and the standard deviation of the spatial differential interference phase value of each side in the triangulation network;
and setting a standard deviation threshold value of the differential interference phase value, and selecting the PS point smaller than the threshold value as the selected PS point with stable phase.
In some preferred embodiments, the one-dimensional time-phase unwrapping comprises: and summing the interference phase values of corresponding pixels in the time-adjacent interference images.
In some preferred embodiments, the computing computes double-difference phase observations of M edges in a triangular mesh
Figure GDA0002914540670000026
The method comprises the following steps: double-difference phase observed value on each edge in triangular network
Figure GDA0002914540670000031
The time phase unwrapping interference phase value of the edge terminal point, i.e. the selected PS point
Figure GDA0002914540670000032
By taking the difference, i.e.
Figure GDA0002914540670000033
In some preferred embodiments, the coefficient matrix a is a [2N (N-1) (N × N) ] matrix that performs spatial phase unwrapping transformation in the x and y directions for the phase of each selected PS point.
In some preferred embodiments, the weight matrix P has M rows and M columns, where each element value PijIs the reciprocal of the corresponding edge length value in the triangular net and the element value p corresponding to the non-main diagonalijIs 0.
The invention also provides a phase unwrapping system for monitoring dam and landslide deformation GB-SAR, which comprises an interference image acquisition module, a one-bit time phase unwrapping module and a two-dimensional space phase unwrapping module;
the interference image acquisition module is connected with a dam and a landslide deformation GB-SAR monitoring system, acquires a GB-SAR interference image, selects N PS points with stable phases and transmits the PS points to the one-bit time phase unwrapping module;
the one-bit time phase unwrapping module is configured to perform one-dimensional time phase unwrapping on each selected PS point to obtain an interference phase value and transmit the interference phase value to the two-dimensional space phase unwrapping module;
and the two-dimensional space phase unwrapping module is configured to perform two-dimensional space phase unwrapping on each selected PS point according to a method equation established by an indirect least square adjustment method.
In some preferred embodiments, the two-dimensional spatial phase unwrapping module further includes a Delaunay triangulation network construction unit, an error matrix construction unit, and a normal equation construction unit;
the Delaunay triangulation network construction unit is set to construct a Delaunay triangulation network according to each selected PS point, calculate a time phase unwrapping value of each edge in the triangulation network, and transmit the time phase unwrapping value to the error matrix construction unit;
the error matrix construction unit is configured to construct an error equation according to an indirect least square adjustment method, convert the error equation into an error matrix and transmit the error matrix to the normal equation construction unit;
the normal equation building unit is configured to build a normal equation according to the error matrix and calculate the unwrapping phase value of each selected PS point.
The invention has the beneficial effects that:
1. a phase unwrapping method is provided for a GB-SAR system of a dam and a landslide in a high mountain canyon, so that the influence of atmospheric disturbance on unwrapping precision is reduced;
2. the method for unwrapping the phase based on the indirect least square adjustment is provided, integration is carried out in a mode of preventing the phase from passing through residual points, errors cannot be transmitted, and unwrapping precision is greatly improved;
3. the selection and the setting of branch tangent lines are avoided, and a global solution method is adopted, so that the convergence speed is accelerated while the island problem is better processed, and the unwrapping process is quicker;
drawings
FIG. 1 is a flow chart of a phase unwrapping method for GB-SAR monitoring of dam and landslide deformation in one embodiment of the present application;
FIG. 2 is a block diagram of a phase unwrapping system for GB-SAR monitoring of dam and landslide deformations in one embodiment of the present application;
fig. 3 is a structural diagram of a two-dimensional spatial phase unwrapping module in a phase unwrapping system for monitoring dam and landslide deformation GB-SAR in an embodiment of the present application;
FIG. 4 is a schematic diagram of a Delaunay triangulation network constructed from selected PS points in a preferred embodiment of the present application;
FIG. 5 is an interference image after temporal and spatial phase unwrapping in a preferred embodiment of the present application;
FIG. 6 is a statistical histogram of errors in phase unwrapping of each PS point in another 1 interference image according to a preferred embodiment of the present application;
FIG. 7 is a statistical histogram of errors in phase unwrapping of each PS point in another 1 interference image according to a preferred embodiment of the present application;
FIG. 8 is a schematic diagram of a Delaunay triangulation network constructed from selected PS points in another preferred embodiment of the present application;
FIG. 9 is an interference image after temporal and spatial phase unwrapping in another preferred embodiment of the present application;
FIG. 10 is a statistical histogram of errors in phase unwrapping of PS points of an interference image according to another preferred embodiment of the present application;
Detailed Description
In order to make the objects, technical solutions and advantages of the present invention more apparent, the present invention is further described with reference to the accompanying drawings. In the description of the present invention, it is to be understood that the terms "upper", "lower", "front", "rear", "left", "right", "top", "bottom", "inner", "outer", and the like, indicate orientations or positional relationships based on the orientations or positional relationships shown in the drawings, are merely for convenience in describing the present invention and simplifying the description, and do not indicate or imply that the device or element being referred to must have a particular orientation, be constructed and operated in a particular orientation, and thus, should not be construed as limiting the present invention.
As shown in fig. 1, a phase unwrapping method for monitoring dam and landslide deformation GB-SAR comprises the steps of:
s1, acquiring a GB-SAR interference image in a deformation monitoring area, selecting N PS points with stable phases, and performing one-dimensional time phase unwrapping to obtain a time phase unwrapping interference phase value
Figure GDA0002914540670000041
N={1,2,3…i,j…N};
It should be understood by those skilled in the art that the deformation monitoring area is specifically an area to be monitored in a dam or a landslide, and it is understood that the application of the technical solution of the present application to other fields of monitoring the micro deformation of the earth surface is also within the protection scope of the present application.
The PS points are permanent scatterer (PermanentScatterers) points in a GB-SAR deformation monitoring area, the backward scattering characteristics of the PS points hardly change along with time, and the PS points have the characteristic of high interference coherence. Based on the scattering characteristics of the PS points, the PS points with stable phases are selected, so that the phase components of all components can be accurately obtained by decomposition from the interference phases, the final deformation information can be obtained, the problems of atmospheric disturbance, time and space loss coherence can be indirectly solved, and the influence of system noise can be weakened to a certain extent.
There are many methods for selecting PS points, such as coherence coefficient thresholding, phase dispersion thresholding, and amplitude dispersion index thresholding. The application does not further limit the specific PS point selection method.
In a preferred embodiment, a comprehensive selection method of PS points is provided, which combines a coherence coefficient threshold method and an amplitude dispersion index threshold method, considers the stability of a target in time while considering the strong scattering characteristic, combines amplitude characteristics and phase characteristics, and proposes a differential interference phase standard threshold method to ensure the stability and continuity in space, thereby selecting the PS points with stable phases. The method specifically comprises the following steps:
calculating a coherence coefficient of the interference image, and selecting pixels meeting high coherence;
calculating the absolute deviation value of the coherent coefficient of each pixel, setting a threshold value T, counting the proportion of pixels with absolute deviation values larger than the threshold value T to the total number of pixels in the image, and rejecting the image with overlarge proportion value;
calculating the amplitude dispersion index of each pixel time sequence, and determining candidate PS points according to the relation between the threshold T and the amplitude dispersion index;
establishing a Delaunay triangulation network according to the candidate PS points, and calculating the absolute difference and the standard deviation of the spatial differential interference phase value of each side in the triangulation network;
and setting a standard deviation threshold value of the differential interference phase value, and selecting the PS point smaller than the threshold value as the selected PS point with stable phase.
Further, the time phase unwrapping is that according to the situation of obtaining interference images in dam and landslide deformation monitoring, one-dimensional time phase unwrapping is firstly carried out, and then two-dimensional phase unwrapping on the space is carried out, so that the three-dimensional unwrapping problem is converted into one-dimensional and two-dimensional problems, and the processing method is simplified. More recently, the algorithm construction mode is adopted, and mainly for solving the interference of weather in the high mountain canyon on the GB-SAR image, the problem of phase ambiguity of most interference phase observation values is solved by one-dimensional time phase unwrapping, and then the position where the phase unwrapping is incomplete or inaccurate is solved by a two-dimensional space phase unwrapping method.
In some preferred embodiments, the one-dimensional time-phase unwrapping comprises: and summing the interference phase values of corresponding pixels in the time-adjacent interference images. The specific method comprises the following steps:
based on the following one-dimensional time-phase unwrapping rationale: for m discrete winding behavior observation sequences
Figure GDA0002914540670000051
Unwrapping interference phase
Figure GDA0002914540670000052
The method can be obtained by performing interference coherence on the single-view complex images acquired at different times and then summing the images:
Figure GDA0002914540670000061
wherein k is 1,2, …, m;
as can be seen from the above formula, it is only necessary to ensure the phase integer ambiguity n between two consecutive monoscopic imagesk-1,kThe unwrapped interference phase value at any time can be calculated as zero
Figure GDA0002914540670000062
It will be appreciated that in radar signal acquisition n is guaranteed as long as the sampling rate of the radar signal is such that phase fringes can be detectedk-1,kZero, so that the correct unwrapping result is obtained. That is, the sampling rate can be used to detect the circumferential fringe as long as the Nyquist criterion is satisfied, so that the interference phase of the continuous image element
Figure GDA0002914540670000063
Thus, the one-dimensional time phase unwrapping of the interference phase is completed preliminarily. In this application, the time phase unwrapping interference phase value obtained by unwrapping the one-dimensional time phase is recorded as
Figure GDA0002914540670000064
Where N ═ {1,2,3 … i, j … N }.
S2, constructing a Delaunay triangulation network by using the selected PS points, and calculating a double-difference phase observation value of M edges in the triangulation network
Figure GDA0002914540670000065
The Delaunay triangulation network is a set of connected but non-overlapping triangles, and the circumcircle of these triangles does not contain any other point of this area, belonging to an irregular triangulation network. The construction of the triangular network meets the optimal triangular condition, even if three inner angles and three side lengths of each formed triangle are approximately equal, namely, the angles can not be large obtuse angles or small acute angles, the three angles are all acute angles as far as possible, and other additional conditions are selected according to specific application requirements to construct the triangular network meeting the requirements. The construction is very numerous and one of them is given in the present application in a preferred embodiment.
The method for establishing the irregular triangulation network (Delaunay triangulation network) of the PS points in the GB-SAR image by adopting an angle judgment method and combining a distance condition is as follows. In the method, GB-SAR image data and the specific requirements of actual dam or landslide detection are added, and parameters are set according to the GB-SAR image data and the specific requirements.
First, a theoretical analysis procedure of the construction is given. When two PS points are known, a base line side can be determined, the inner angle size of the triangle with the third adjacent alternative PS point as the corner vertex is calculated by using the cosine theorem, and the PS point corresponding to the maximum point is selected as the third vertex of the triangle. Further, for monitoring the dam and landslide GB-SAR, the spatial autocorrelation of the atmospheric states of adjacent PS points in the same image and the density and the distance of the point positions during deformation analysis need to be considered in the triangulation network of the triangulation network, so that the distance value between the adjacent PS points needs to be considered. I.e. the requirement of the side length also needs to be taken into account in case the angle requirement is met.
Secondly, a specific flow of construction is given:
(1) the obtained PS point data is blocked, so that only the points adjacent to the target triangle are searched in the subsequent steps, and all data does not need to be searched. And (3) searching all adjacent points, and setting different intervals aiming at different areas according to the change condition of the atmosphere in the environment where the dam and the landslide are positioned and the requirement of object deformation characteristic analysis.
(2) And determining a network construction starting point. A point a is randomly selected from the plurality of discrete PS points as a first point of the detection triangle. In order to ensure the consistency of the operation results and avoid the inconsistency of the triangulation networks formed at different times under the same data source and influence on subsequent data processing and computational analysis, the position of the point a may be fixed, for example, by taking the first PS point or the first PS point in the first search grid. Further, according to the minimum distance principle and the actual situation of the dam or the landslide site, a point with the distance A being less than 30 meters is selected as a second point B of the triangle. Then PS point C near the triangle side formed by A, B pointsiCalculating the angle C by using the cosine theoremi
Figure GDA0002914540670000071
Wherein, ai=BCi;bi=ACi;ci=ABi. If < C > max { < C { (C) }iAnd selecting the point C as the other point of the triangle.
(3) The triangle is expanded. And (3) expanding outwards from the triangle constructed in the step (2) to form a triangular net. The expansion process follows the principle of maximum angle, the requirement of limiting the length of the side of the triangle is considered, the correctness of the spatial topological relation of the triangular network is ensured, and the problems of repetition, capping, intersection and the like are avoided. The specific expanding method comprises the following steps: and expanding the newly added two sides of each generated triangle outwards according to the angle maximum principle, and checking the topological relation such as repetition, intersection and the like.
The method for outward expansion comprises the following steps:
if the vertex is P1(x1,y1),P2(x2,y2),P3(x3,y3) Triangle P1P2The edge is extended outwards and should be located on a straight line P1P2And P3Points on opposite sides. P1P2The equation of a straight line is
F(x,y)=(y2-y1)(x-x1)-(x2-x1)(y-y1)=0
If the coordinate of the alternative point P is (x, y), then F (x, y) is F (x)3,y3) When < 0, P and P3On a straight line P1P2On condition that the distance between P and P1 and P2 is less than 30 meters, the point can be used as an alternative expansion vertex.
The method for checking the topological relation comprises the following steps:
because any side can only be a common side of two triangles at most, only the expansion times of each side need to be recorded, and when the expansion times of the side exceeds 2 times, the expansion is invalid; otherwise, the sub-extension is valid. When all the new edges of the generated triangle are subjected to the expansion processing, all the discrete data points are connected into an irregular triangular network.
Thus, the description of the construction of the Delaunay triangulation network is completed.
In a preferred embodiment, the computing of double-difference phase observations of M edges in a triangular mesh is performed
Figure GDA0002914540670000081
The method comprises the following steps:
double-difference phase observed value on each edge in triangular network
Figure GDA0002914540670000082
The time phase unwrapping interference phase value of the edge terminal point, i.e. the selected PS point
Figure GDA0002914540670000083
By taking the difference, i.e.
Figure GDA0002914540670000084
S3, constructing an error equation based on an indirect least square adjustment method
Figure GDA0002914540670000085
Wherein
Figure GDA0002914540670000086
The unwrapped phase value for each selected PS point.
The indirect adjustment method is an adjustment calculation method widely used in the measurement field, and when determining the most probable value of a plurality of unknowns, an independent quantity without any conditional relation between them is selected as an unknowns composition, a function relation of the unknowns expression measurement is formed, an error equation is listed, and the most probable value of the unknowns is obtained according to the principle of least square method. The theoretical derivation process of the method is as follows: the t independent unknowns having a certain relation with the observed values are selected as parameters, each observed value is respectively expressed as a function of the t parameters, a function model is established, and the most probable value of the parameters is solved by a method of solving a free extreme value according to the least square principle, so that the adjustment value of each observed value is obtained.
Specifically, in the present application, an error equation is first constructed by using an indirect least squares adjustment method according to observation conditions, and unwrapping interference phase values of selected PS points are selected
Figure GDA0002914540670000087
For parameters, the mathematical model of phase unwrapping can be uniquely determined by these N independent parameters:
Figure GDA0002914540670000088
it should be understood that each edge of the triangle corresponds to an error equation. As can be seen from the preceding steps,
Figure GDA0002914540670000089
therefore, the temperature of the molten metal is controlled,
Figure GDA00029145406700000810
further, when the indirect adjustment method is applied to phase unwrapping of GB-SAR monitoring of dam and landslide deformation, an indirect adjustment mathematical model needs to be constructed by combining the characteristics of a monitoring environment and GB-SAR monitoring equipment, and relevant descriptions are provided in the application specification.
S4, converting the error equation of the M edges into a matrix expression
Figure GDA00029145406700000811
Wherein V is an error matrix; a is a coefficient matrix; phi is a parameter matrix and corresponds to the unwrapping phase values of the N selected PS points; and L is a parameter matrix corresponding to the time phase unwrapping value of the M edges.
So far, the unwrapping phase value of the selected PS point can be obtained by solving the elements in the parameter matrix phi, and the phase unwrapping of the interference image is completed.
The elements in the matrix expression are further described below.
The coefficient matrix a includes M rows and N columns, and here, the values of the elements of the a matrix are determined by the pattern composition condition of the triangular network. In a preferred embodiment, the coefficient matrix a is a [2N (N-1) (N × N) ] matrix for performing spatial phase unwrapping transformation on the phase of each selected PS point in x and y directions, and the specific expression is as follows:
Figure GDA0002914540670000091
further, the parameter matrix Φ corresponding to the N selected PS point unwrapped phase values specifically includes: the N selected PS are obtained by double-differencing with other points in sequence, so that NN elements are shared in total, and the specific expression is as follows:
Figure GDA0002914540670000092
further, the parameter matrix L corresponding to the M-edge double-difference phase observed value specifically includes: n selected PS are subjected to double difference with other points in sequence, and the double difference between the selected PS and the other points is removed, so that 2N (N-1) elements are shared in total, and the specific expression is as follows:
Figure GDA0002914540670000101
s5, establishing a normal equation ATPA-ATPL is 0, then Φ is (a)TPA)-1ATPL, where P is the weight matrix.
It should be understood that the parameter matrix Φ needs to satisfy V according to the basic principle of the least squares methodTThe corresponding equation A can be established according to the requirement of PV ═ minTPA-ATAnd PL is 0, wherein P is a weight matrix of M rows and M columns obtained according to the relation between points of the triangular net.
In some preferred embodiments, each element P in the weight matrix P is determined according to the size of each side length in the triangular meshijIn which the element p on the main diagonalijIs the corresponding edge length value DijIs inverse of (i.e.
Figure GDA0002914540670000102
Further, in other preferred embodiments, the phase values of the selected PS points at the endpoints of the triangular network are independent from each other, so that the element P on the non-principal diagonal in the weight matrix P isijIs 0.
Thus, a weight matrix P can be determined, and further, the true interference phase value of each selected PS point can be determined according to the equation: phi ═ aTPA)-1ATPL。
The invention provides a phase unwrapping method aiming at a GB-SAR system of a dam and a landslide in a high mountain canyon, so that the influence of atmospheric disturbance on unwrapping precision is reduced; integration is carried out in a mode of avoiding 'passing through' residual points, errors cannot be transmitted, and the unwinding precision is greatly improved; the method has better processing effect on the island problem, and simultaneously accelerates the convergence speed, so that the unwrapping process is faster.
In some preferred embodiments, as shown in fig. 2, the present invention further discloses a phase unwrapping system for GB-SAR monitoring of dam and landslide deformations, comprising: the device comprises an interference image acquisition module, a one-bit time phase unwrapping module and a two-dimensional space phase unwrapping module;
the interference image acquisition module is connected with a dam and a landslide deformation GB-SAR monitoring system, acquires a GB-SAR interference image, selects N PS points with stable phases and transmits the PS points to the one-bit time phase unwrapping module;
the one-bit time phase unwrapping module is configured to perform one-dimensional time phase unwrapping on each selected PS point to obtain an interference phase value and transmit the interference phase value to the two-dimensional space phase unwrapping module;
and the two-dimensional space phase unwrapping module is configured to perform two-dimensional space phase unwrapping on each selected PS point according to a method equation established by an indirect least square adjustment method.
In other preferred embodiments, as shown in fig. 3, the two-dimensional spatial phase unwrapping module further includes a Delaunay triangulation network constructing unit, an error matrix constructing unit, and a normal equation constructing unit;
the Delaunay triangulation network construction unit is set to construct a Delaunay triangulation network according to each selected PS point, calculate a time phase unwrapping value of each edge in the triangulation network, and transmit the time phase unwrapping value to the error matrix construction unit;
the error matrix construction unit is configured to construct an error equation according to an indirect least square adjustment method, convert the error equation into an error matrix and transmit the error matrix to the normal equation construction unit;
the normal equation building unit is configured to build a normal equation according to the error matrix and calculate the unwrapping phase value of each selected PS point.
In the following, the invention provides a theoretical process for studying the precision of phase unwrapping by the least square method.
According to the error theory, in order to evaluate the observation accuracy, the calculation unit is needed firstlyEstimation of weight variance
Figure GDA0002914540670000111
And then calculating a covariance matrix of each observed value, and finally calculating the precision of the parameter matrix phi according to the covariance propagation law and the error propagation law.
Error value v of differential phase observations from the edge of a triangulation networki,jThe unit weight variance can be calculated:
Figure GDA0002914540670000112
the error in unit weight for the differential phase observations thus obtained is:
Figure GDA0002914540670000113
the variance of each differential phase observation can be obtained according to the weighting method of the differential phase observations:
Figure GDA0002914540670000114
the above formula can also be expressed by the side length value in the triangular net:
Figure GDA0002914540670000115
from this, we can get the double difference phase observation on each side
Figure GDA0002914540670000116
The medium error of (2):
Figure GDA0002914540670000117
in order to determine a covariance matrix of the parameter matrix phi, the mean error of each PS-point interference phase unwrapping value is obtained. The co-factor matrix Q of the parameter phi can be obtained according to the co-factor propagation lawΦΦ=(ATPL)-1. Further, the covariance of the parameter Φ is obtained from the relationship between the covariance matrix and the covariance matrix
Figure GDA0002914540670000121
Wherein, covariance matrixDΦΦThe elements on the middle main diagonal line correspond to the variance of the spatial unwrapping phase of each PS point, and the accuracy of the unwrapping value can be obtained.
In order to verify the effectiveness of the phase unwrapping method provided by the invention, the following is given of a GB-SAR image phase unwrapping embodiment applied to deformation monitoring of a GY arch dam and a HS large-scale landslide body.
Example 1
In this embodiment, the phase unwrapping is performed according to the technical scheme provided by the invention on the GB-SAR monitoring image obtained by detecting the GY hydropower station engineering gravity arch dam of the north and the lake of china for a total of 147 hours for 7 days, and the unwrapping precision is analyzed.
As shown in fig. 4, the generated interference image is obtained, 4289 PS points are extracted, and one-dimensional time phase unwrapping is performed on the PS points to obtain time phase unwrapped interference phase values, so as to construct a Delaunay triangulation network, and 14433 baselines and 9617 triangles are obtained in total.
For the time phase unwrapping interference phase value of 37 GB-SAR monitoring interference images in the dam body part, according to the phase unwrapping method provided by the invention, an error equation is established and an unwrapping value is solved, so that the average mean error of spatial phase unwrapping of 37 images is +/-0.004 mrad, the maximum is +/-0.027 mrad, and the minimum is +/-0.001 mrad.
As shown in fig. 5, the 1 st and 9 th images are taken as examples before and after the temporal and spatial phase unwrapping. Fig. 5a shows the interference image with the 1 st phase winding, which has a serious phase winding, because the data acquisition is still in a debugging state at the beginning of the experiment, the phase continuity between two adjacent time instants is poor, the interference phase value of some pixels exceeds 2 pi radian, and the fuzzy number of the whole cycle is not 0. As can be seen from fig. 5b, before and after the time phase unwrapping, there are still large discontinuities in the phase values at many points, and there are large jumps, indicating that the unwrapping of these points is not complete and further unwrapping is required. After the phase unwrapping method provided by the invention is used for carrying out time phase unwrapping and then carrying out space phase unwrapping based on indirect least square adjustment, the unwrapped image shown in figure 5c is obtained, and as can be seen from the figure, each PS point has very good spatial continuity and is in an interval of [ -pi, pi). FIG. 5d shows the 9 th phase-wrapped interference image during the steady monitoring period, with better phase temporal continuity. After the Itoh time phase is adopted for unwrapping, as can be seen from FIG. 5e, the time phase unwrapping result at this time is better, the fuzzy number of the whole cycle contained in the interference phase of each PS point is basically 0, only the phase at the individual point position does not meet the requirement of consistency in space, the space phase unwrapping method of the text is continuously adopted for unwrapping, the final unwrapping result shown in FIG. 5f is obtained, and at this time, all the PS points are better consistent in space, and a better unwrapping result is obtained.
According to the error analysis method provided by the present invention, the phase unwrapping errors of the PS points in the 1 st and 9 th interference images are calculated respectively, and the distribution thereof is counted, as shown in fig. 6 and 7. As can be seen from the statistical histogram of the errors, most of the PS point unwrapping errors in the 1 st interference image are concentrated between 0.002mrad and 0.004mrad, and most of the PS point unwrapping errors in the 9 th interference image are concentrated between 0.001mrad and 0.003 mrad. As can be seen from fig. 6 and 7, both the unwrapping results have no gross error, the error distribution conforms to the distribution characteristics of the accidental errors, the unwrapping results have strong robustness, and the unwrapping precision is high. Because the 1 st interference image has larger fluctuation of the phase observation value during data acquisition and unstable data, the error value is slightly larger than the unwrapping precision of the 9 th interference image, which also shows that the phase unwrapping method adopted by the text can fully reflect the quality of the data, thereby effectively unwrapping the data.
Example 2
In this embodiment, a GB-SAR monitoring image obtained by monitoring a large landslide body in the three gorges reservoir region for 2 months from 20 days 4 and 2019 and 20 days 6 and 2019 for 20 days is subjected to phase unwrapping according to the technical solution provided by the present invention, and its unwrapping accuracy is analyzed. The basic situation of the monitored subject is as follows: the elevation of the front edge is about 65.0m, the elevation of the rear edge is about 840.0m, and the total volume is about 1800 ten thousand cubic meters. The terrain surface is basically parallel to the bedrock surface, and the gradient is steeper.
As shown in fig. 8, PS point selection is performed on the interference image, and a PS triangulation network is constructed. According to the phase unwrapping method provided by the invention, an error equation is established and an unwrapping value is calculated, so that the average mean error of phase unwrapping in the whole monitoring period is +/-0.031 mrad, maximally +/-0.301 mrad and minimally +/-0.001 mrad.
As shown in fig. 9, a GB-SAR monitoring image acquired at 18 pm at 30 pm in 4/2019 is taken as an example. As shown in FIG. 9a, the phase of each PS point in the original interference image is discontinuous and is in the range of [ - π, π). As shown in fig. 9b, the result of the one-dimensional time phase unwrapping performed by the phase unwrapping method according to the present invention is that the monitoring time is long, and weather changes such as rainfall occur many times in the process, so that the PS points vary in space, the spatial consistency condition cannot be strictly met, and the phases of some PS points are not continuous in space after the time phase unwrapping. After two-space phase unwrapping, the result shown in fig. 9c is obtained, at which point the PS dots are perfectly continuous in space.
According to the error analysis method provided by the present invention, the phase unwrapping errors of each PS point in the image are calculated respectively, and the distribution of the phase unwrapping errors is counted as shown in FIG. 10, wherein most of the unwrapping errors of the PS points are concentrated between + -0.06 mrad and + -0.18 mrad. Compared with GY arch dam deformation monitoring, the slip mass has long monitoring time, the phase observation value is more obviously disturbed by the atmosphere, and the unwrapping precision is greatly influenced by the atmospheric disturbance. However, the overall unwrapping precision analysis shows that the unwrapping result has no gross error, conforms to the distribution characteristics of accidental errors and still has high unwrapping precision. The conclusion shows that the phase unwrapping method provided by the invention has better applicability to spatial phase unwrapping in continuous long-time landslide mass monitoring data.
The foregoing shows and describes the general principles, essential features, and advantages of the invention. It will be understood by those skilled in the art that the present invention is not limited to the embodiments described above, which are described in the specification and illustrated only to illustrate the principle of the present invention, but that various changes and modifications may be made therein without departing from the spirit and scope of the present invention, which fall within the scope of the invention as claimed. The scope of the invention is defined by the appended claims and equivalents thereof.

Claims (8)

1. A phase unwrapping method for GB-SAR monitoring of dam and landslide deformation is characterized by comprising the following steps:
obtaining GB-SAR interference images in a deformation monitoring area, selecting N PS points with stable phases, and performing one-dimensional time phase unwrapping to obtain a time phase unwrapped interference phase value
Figure FDA0002807126020000011
Constructing a Delaunay triangulation network by using the selected PS points, and calculating a double-difference phase observation value of M edges in the triangulation network
Figure FDA0002807126020000012
Error equation construction based on indirect least square adjustment method
Figure FDA0002807126020000013
Wherein
Figure FDA0002807126020000014
Unwrapped phase values for each selected PS point;
converting error equation of M edges into matrix expression
Figure FDA0002807126020000015
Wherein V is an error matrix; a is a coefficient matrix; phi is a parameter matrix and corresponds to the unwrapping phase values of the N selected PS points; l is a parameter matrix and corresponds to the time phase unwrapping value of the M edges;
establishing a normal equation ATPA-ATPL is 0, then Φ is (a)TPA)-1ATPL, where P is the weight matrix.
2. The phase unwrapping method according to claim 1, wherein the selecting of the N PS points with stable phase comprises the steps of:
calculating a coherence coefficient of the interference image, and selecting pixels meeting high coherence;
calculating the absolute deviation value of the coherent coefficient of each pixel, setting a threshold value T, counting the proportion of pixels with absolute deviation values larger than the threshold value T to the total number of pixels in the image, and rejecting the image with overlarge proportion value;
calculating the amplitude dispersion index of each pixel time sequence, and determining candidate PS points according to the relation between the threshold T and the amplitude dispersion index;
establishing a Delaunay triangulation network according to the candidate PS points, and calculating the absolute difference and the standard deviation of the spatial differential interference phase value of each side in the triangulation network;
and setting a standard deviation threshold value of the differential interference phase value, and selecting the PS point smaller than the threshold value as the selected PS point with stable phase.
3. The phase unwrapping method of claim 1, wherein the one-dimensional time phase unwrapping includes: and summing the interference phase values of corresponding pixels in the time-adjacent interference images.
4. The phase unwrapping method of claim 1, wherein the computing of double-difference phase observations of M edges in a triangular mesh
Figure FDA0002807126020000016
The method comprises the following steps: double-difference phase observed value on each edge in triangular network
Figure FDA0002807126020000017
The time phase unwrapping interference phase value of the edge terminal point, i.e. the selected PS point
Figure FDA0002807126020000018
By taking the difference, i.e.
Figure FDA0002807126020000019
5. The phase unwrapping method as recited in claim 1, wherein: the coefficient matrix A is a [2N (N-1) (N × N) ] matrix which performs spatial phase unwrapping transformation in the x and y directions according to the phase of each selected PS point.
6. The phase unwrapping method as recited in claim 1, wherein: the weight matrix P has M rows and M columns, wherein each element value PijIs the reciprocal of the corresponding edge length value in the triangular net and the element value p corresponding to the non-main diagonalijIs 0.
7. A phase unwrapping system that is used for dam and landslide to warp GB-SAR monitoring which characterized in that: the device comprises an interference image acquisition module, a one-bit time phase unwrapping module and a two-dimensional space phase unwrapping module;
the interference image acquisition module is connected with a dam and a landslide deformation GB-SAR monitoring system, acquires a GB-SAR interference image, selects N PS points with stable phases and transmits the PS points to the one-bit time phase unwrapping module;
the one-bit time phase unwrapping module is configured to perform one-dimensional time phase unwrapping on each selected PS point to obtain an interference phase value and transmit the interference phase value to the two-dimensional space phase unwrapping module;
and the two-dimensional space phase unwrapping module is configured to perform two-dimensional space phase unwrapping on each selected PS point according to a method equation established by an indirect least square adjustment method.
8. The bit unwrapping system as recited in claim 7, wherein: the two-dimensional space phase unwrapping module further comprises a Delaunay triangulation network construction unit, an error matrix construction unit and a normal equation construction unit;
the Delaunay triangulation network construction unit is set to construct a Delaunay triangulation network according to each selected PS point, calculate a time phase unwrapping value of each edge in the triangulation network, and transmit the time phase unwrapping value to the error matrix construction unit;
the error matrix construction unit is configured to construct an error equation according to an indirect least square adjustment method, convert the error equation into an error matrix and transmit the error matrix to the normal equation construction unit;
the normal equation building unit is configured to build a normal equation according to the error matrix and calculate the unwrapping phase value of each selected PS point.
CN202011372417.2A 2020-11-30 2020-11-30 Phase unwrapping method and system for monitoring dam and landslide deformation GB-SAR Active CN112698328B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011372417.2A CN112698328B (en) 2020-11-30 2020-11-30 Phase unwrapping method and system for monitoring dam and landslide deformation GB-SAR

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011372417.2A CN112698328B (en) 2020-11-30 2020-11-30 Phase unwrapping method and system for monitoring dam and landslide deformation GB-SAR

Publications (2)

Publication Number Publication Date
CN112698328A CN112698328A (en) 2021-04-23
CN112698328B true CN112698328B (en) 2021-08-10

Family

ID=75506004

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011372417.2A Active CN112698328B (en) 2020-11-30 2020-11-30 Phase unwrapping method and system for monitoring dam and landslide deformation GB-SAR

Country Status (1)

Country Link
CN (1) CN112698328B (en)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114595192B (en) * 2022-03-10 2023-02-28 青海省地质调查院 Intelligent data real-time gathering method and system suitable for regional geological survey
CN115308656B (en) * 2022-08-09 2024-05-28 上海电气控股集团有限公司智惠医疗装备分公司 Method and device for determining B0 field diagram, electronic equipment and storage medium
CN117435902B (en) * 2023-12-20 2024-04-02 武汉华威科智能技术有限公司 Method and device for determining RFID tag movement behavior based on machine learning

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108957456A (en) * 2018-08-13 2018-12-07 伟志股份公司 Landslide monitoring and EARLY RECOGNITION method based on multi-data source SBAS technology
CN111474544A (en) * 2020-03-04 2020-07-31 广东明源勘测设计有限公司 Landslide deformation monitoring and early warning method based on SAR data

Family Cites Families (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6150973A (en) * 1999-07-27 2000-11-21 Lockheed Martin Corporation Automatic correction of phase unwrapping errors
CN104091064B (en) * 2014-07-02 2017-02-22 北京航空航天大学 PS-DInSAR ground surface deformation measurement parameter estimation method based on optimal solution space search method
CN104765026A (en) * 2015-04-29 2015-07-08 天津市测绘院 Method for extracting ground attribute data in interferometry synthetic aperture radar data
CN105158761A (en) * 2015-08-31 2015-12-16 西安电子科技大学 Radar synthetic phase unwrapping method based on branch-cut method and surface fitting
CN105738896B (en) * 2016-02-25 2018-06-19 内蒙古工业大学 A kind of ground SAR Slope with multi-step interferometric phase unwrapping methods and device
CN108627833B (en) * 2018-05-15 2021-08-24 电子科技大学 Atmospheric phase compensation method based on GB-InSAR
CN109031301A (en) * 2018-09-26 2018-12-18 云南电网有限责任公司电力科学研究院 Alpine terrain deformation extracting method based on PSInSAR technology
CN109884635B (en) * 2019-03-20 2020-08-07 中南大学 Large-range high-precision InSAR deformation monitoring data processing method
CN110109106A (en) * 2019-04-23 2019-08-09 中国电力科学研究院有限公司 A kind of InSAR interferometric phase unwrapping method in region with a varied topography
CN110425995B (en) * 2019-08-30 2020-04-14 四川大学 Landslide deformation monitoring method based on point cloud average domain vector algorithm

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108957456A (en) * 2018-08-13 2018-12-07 伟志股份公司 Landslide monitoring and EARLY RECOGNITION method based on multi-data source SBAS technology
CN111474544A (en) * 2020-03-04 2020-07-31 广东明源勘测设计有限公司 Landslide deformation monitoring and early warning method based on SAR data

Also Published As

Publication number Publication date
CN112698328A (en) 2021-04-23

Similar Documents

Publication Publication Date Title
CN112698328B (en) Phase unwrapping method and system for monitoring dam and landslide deformation GB-SAR
CN110058237B (en) InSAR point cloud fusion and three-dimensional deformation monitoring method for high-resolution SAR image
CN109917356B (en) Airborne laser scanning system error calibration method
Wang et al. Improved SAR amplitude image offset measurements for deriving three-dimensional coseismic displacements
Catalão et al. Merging GPS and atmospherically corrected InSAR data to map 3-D terrain displacement velocity
Liu et al. Estimating Spatiotemporal Ground Deformation With Improved Persistent-Scatterer Radar Interferometry $^\ast$
CN112986993B (en) InSAR deformation monitoring method based on space constraint
CN113340191B (en) Time series interference SAR deformation quantity measuring method and SAR system
CN103454636B (en) Differential interferometric phase estimation method based on multi-pixel covariance matrixes
CN115201825B (en) Atmospheric delay correction method in InSAR (interferometric synthetic aperture radar) inter-seismic deformation monitoring
CN109239710B (en) Method and device for acquiring radar elevation information and computer-readable storage medium
CN112444188B (en) Multi-view InSAR sea wall high-precision three-dimensional deformation measurement method
CN114689015B (en) Method for improving elevation precision of optical satellite stereoscopic image DSM
Yin et al. Multi-dimensional and long-term time series monitoring and early warning of landslide hazard with improved cross-platform SAR offset tracking method
CN115877421A (en) Deformation detection method and device for geological sensitive area of power transmission channel
Kampes et al. Velocity field retrieval from long term coherent points in radar interferometric stacks
Liu et al. Dynamic estimation of multi-dimensional deformation time series from Insar based on Kalman filter and strain model
Li et al. A precise underwater positioning method by considering the location difference of transmitting and receiving sound waves
Seube et al. Multibeam echo sounders-IMU automatic boresight calibration on natural surfaces
Wang et al. Accurate persistent scatterer identification based on phase similarity of radar pixels
Duan et al. Adaptively selecting interferograms for SBAS-InSAR based on graph theory and turbulence atmosphere
CN103760582B (en) A kind of optimization method blocking satellite double-difference observation structure under environment
CN112363166A (en) InSAR phase unwrapping method and system based on reliable redundant network
Estahbanati et al. A phase unwrapping approach based on extended Kalman filter for subsidence monitoring using persistent scatterer time series interferometry
CN116184463A (en) Wide swath sea surface wind field inversion method, device, equipment and medium

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
CB03 Change of inventor or designer information
CB03 Change of inventor or designer information

Inventor after: Xiang Xia

Inventor after: Chen Chen

Inventor after: Liu Chao

Inventor after: Nie Ruihua

Inventor after: Wang Hui

Inventor after: Yang Zhengli

Inventor after: Chen Jiankang

Inventor before: Xiang Xia

Inventor before: Chen Chen

Inventor before: Liu Chao

Inventor before: Nie Ruihua

Inventor before: Wang Hui

Inventor before: Yang Zhengli

Inventor before: Chen Jiankang

GR01 Patent grant
GR01 Patent grant