RESEARCH ON COORDINATE TRANSFORMATION METHOD OF GB-SAR IMAGE SUPPORTED BY 3D LASER SCANNING TECHNOLOGY

In the image plane of GB-SAR, identification of deformation distribution is usually carried out by artificial interpretation. This method requires analysts to have adequate experience of radar imaging and target recognition, otherwise it can easily cause false recognition of deformation target or region. Therefore, it is very meaningful to connect two-dimensional (2D) plane coordinate system with the common three-dimensional (3D) terrain coordinate system. To improve the global accuracy and reliability of the transformation from 2D coordinates of GB-SAR images to local 3D coordinates, and overcome the limitation of traditional similarity transformation parameter estimation method, 3D laser scanning data is used to assist the transformation of GB-SAR image coordinates. A straight line fitting method for calculating horizontal angle was proposed in this paper. After projection into a consistent imaging plane, we can calculate horizontal rotation angle by using the linear characteristics of the structure in radar image and the 3D coordinate system. Aided by external elevation information by 3D laser scanning technology, we completed the matching of point clouds and pixels on the projection plane according to the geometric projection principle of GB-SAR imaging realizing the transformation calculation of GB-SAR image coordinates to local 3D coordinates. Finally, the effectiveness of the method is verified by the GB-SAR deformation monitoring experiment on the high slope of Geheyan dam. * Corresponding author


INTRODUCTION
Ground-based Synthetic Aperture Radar (GB-SAR) interferometry is a ground active microwave remote sensing technology which has been developing for more than 20 years.Compared with the spaceborne InSAR technology, the structure of GB-SAR sensor is much simpler.There is no spatial baseline in continuous monitoring mode of GB-SAR.Therefore, when processing GB-SAR image data sequences there is generally no need to consider the influence of spatial baseline error, flat earth phase error and terrain phase error (Pieraccini et al., 2004;Rödelsperger, 2011).In the application of deformation monitoring in small ground surface or engineering structures, such as dams, high slopes, dangerous rocks and landslides, GB-SAR can set up specific geometric scenes for monitoring targets and determine time interval of image acquisition flexibly according to the deformation characteristics of different targets (Margottini et al., 2015;Atzeni et al., 2015;Tarchi et al., 2005;Traglia et al., 2014).GB-SAR is an effective complement to spaceborne SAR interferometry in practice of displacement monitoring.It can be applied well when spaceborne SAR is difficult to meet the requirements of spatial and temporal resolution and accuracy.
At present, the synthetic aperture of GB-SAR system is usually realized by a radar sensor moving on a fixed length guide at a constant speed.The radar sensor transmits and receives modulated electromagnetic signals during orbit movement, collects echo signals at various sampling positions and frequencies, and achieves azimuth resolution through digital signal processing in later stage (Tarchi et al., 2003).The location of the guideway directly determines the axis of the imaging plane, and then affects the projection calculation from 2D pixel coordinates of the GB-SAR image to the actual local 3D coordinates.Unlike the orthographic projection of traditional topographic maps and the central projection of photogrammetry, GB-SAR system can distinguish different targets according to two main parameters.One is the oblique distance of the target to the radar centre, the other is the deviation angle from the centre line of the radar antenna beam.Moreover, the maximum detection distance of GB-SAR is generally small, for example, 10 kilometres.As a result, the deformation target is generally located in the near field area of radar antenna radiation, which does not meet the far-field approximation condition, forming a special sector grid coordinate system of GB-SAR image finally.
In the 2D plane image of GB-SAR, identification of deformation distribution is usually carried out by artificial interpretation.This method requires analysts to have adequate experience of radar imaging and target recognition, otherwise it can easily cause false recognition of deformation target or region.Therefore, it is very meaningful to connect 2D plane coordinate system with the common 3D terrain system, not only helps to improve the accuracy of deformation identification, but also contribute to the analysis of long-term deformation trend and integration of various monitoring results from different means.
In fact, the special projection mode of GB-SAR imaging leads to image coordinate migration relative to the actual plane coordinates.The offset is different for the topographic points with different plane positions and heights.To achieve high accuracy conversion from GB-SAR image coordinates to local 3D terrain coordinates, we need to use external elevation data (such as DTM/DEM, terrain point cloud, etc.) (Martínez, 2008).
On the basis of analysing the geometric projection principle of GB-SAR imaging, this paper expounds the conversion steps of image pixel coordinate to local 3D coordinate.First, collect the point cloud data of the monitoring area, transform the point cloud coordinate into local 3D coordinate system by using common points in both coordinate, and then establish the terrain surface model.According to the same projection method as GB-SAR imaging, the projection calculation of the 3D laser point cloud data into such imaging plane is completed.In this process, we propose a new method to calculate the horizontal angle between GB-SAR and local coordinate by using characteristic lines of monitoring region, which could improve the reliability of the calculation of the horizontal rotation angle.In the GB-SAR imaging plane, a minimum distance criterion is used to match the pixel and the projected point cloud, and the 3D coordinate of the projection point cloud are directly given to the corresponding pixel.An experimental campaign carried out with a GB-SAR to monitor a high slope in dam area, 3D point cloud data were collected at the same time.We complete the 3D coordinate transformation processing of the GB-SAR image verifying the reliability of the methods and steps in this paper.Obviously the results are beneficial to the 3D interpretation and analysis of the deformation, so as to better understand the temporal and spatial regularity of deformation in actual GB-SAR monitoring applications.

GENERAL METHOD
Point cloud data collected directly by 3D laser scanning technology is fixed in 3D coordinate system which takes the scanner centre as the original point.For the uniformity of coordinate system between radar image and point cloud, the point cloud should be first converted to 3D coordinate system of the terrain.And then point cloud data in 3D terrain coordinate system can be converted into plane coordinates in radar projection plane according to projection principle of GB-SAR imaging.Therefore by using matching relation between this plane coordinates and GB-SAR image pixels, 3D spatial coordinate of any pixel can be calculated.The main processing steps are as follows.After determining 3D terrain coordinates of radar centre, 3D movement can make point cloud have the same coordinate origin as the GB-SAR image.But at this time, the direction of plane coordinate axes of such two planes is not consistent, so it is necessary to calculate the horizontal rotation angle between 2D GB-SAR image coordinate system and 3D terrain coordinate system firstly.Figure 1 shows the relationship between the GB-SAR image plane coordinate system and the plane component of terrain coordinate system.Plane x-y in Figure 1 is 2D plane coordinate system of radar pixel, while plane E-N is the plane component of 3D terrain coordinate system.θ is the angle that target point deviates from the direction of radar centre line in GB-SAR image plane.αx-y is the polar angle of target point in the x-y coordinate system and αE-N is the polar angle of target point in the E-N coordinate system.βE-N is the horizontal rotation angle of the two different coordinate system we need to calculate.

Horizontal Angle Calculation
Generally, the calculation of the horizontal rotation angle is carried out by the common point (also known as control point) (Martínez, 2008).Echo signal of such a common point has a certain stability, meanwhile, it has a relative high thermal signal to noise ratio (TSNR), easily distinguished from other surrounding surface pixels.The scope of the common point should be as small as possible to ensure that its central position is determined with high accuracy.However, in the actual monitoring work, it is difficult to have enough ideal natural surface as a common point.There is a need to set up some corner reflector and measure its 3D terrain coordinates, at the same time, identify the radar plane coordinates of the corresponding pixels in the GB-SAR image.

Polar Coordinate Conversion for 3D Point Cloud
It is assumed that there is a pixel point P located in monitoring region of radar image, its radar plane coordinates are (xP, yP), and the corresponding 3D terrain coordinates are (EP, NP, ZP).
The radar plane coordinates of the point P (xP, yP) can also be expressed in polar form.
where rP is the distance between the target pixel P and the radar centre, and θP is the angle of the target pixel P deviating from the radar centre line.
Calculation of 3D coordinates of image pixel P is actually to obtain its corresponding point in 3D terrain points from radar image.In order to achieve the above purpose, all 3D terrain point cloud data is usually transformed into polar coordinates on a 2D plane according to the formula (2).The polar coordinates of the 3D terrain data can be compared directly with the polar coordinates of radar image pixel in the same plane coordinate system.This process is called polar coordinate conversion for 3D point cloud coordinates.
(2) where P (E, N, Z) is 3D coordinate of the terrain point.(Eradar, Nradar, Zradar) is coordinate of the radar centre in the 3D terrain coordinate system, and the βE-N is the horizontal rotation angle of the axis system.rterrain is horizontal distance between terrain point P and radar centre in projection plane of radar imaging, θterrain is the horizontal angle of the terrain point P deviating from radar centre line.Through this calculation, 3D coordinate system of the terrain point cloud relative to the radar centre and plane coordinate system of radar image has been unified.After transforming 3D terrain data into polar coordinates, the corresponding terrain point of P (rP, θP) are found in the polar coordinate system according to minimum distance criterion, and 3D coordinates of this terrain point are taken as the coordinates of pixel P (xP, yP), as shown in figure 2. In order to ensure the reliability of the results, when extracting 3D information of pixels, the first few points with minimum distance can be selected in the target terrain point cloud, and the average value of 3D coordinates is calculated as the final coordinates of pixels.

Azimuth difference calculation of common point
The horizontal angle is generally calculated by common point between different coordinate systems.The selected common points should have coordinates in the different coordinate systems at the same time, and all of them can be clearly and accurately identified in both.It is usually possible to choose the rock mass with higher reflection strength and the edge point of low reflection intensity as a common point for the calculation of horizontal angle.But in actual monitoring, sometimes it is difficult to have a small enough natural surface as a control point, a corner reflector is generally used as a common point.It is necessary to measure the 3D spatial coordinates of the corner reflector, identify the corresponding pixel in radar image and calculate its radar plane coordinates.
Gauss plane coordinates are usually used in the plane coordinates of the 3D terrain coordinate system.Suppose that the plane coordinates of common point in the local Gauss coordinate system are (Xref, Yref), Xref is North coordinate, and Yref East.Plane coordinates of radar centre in this terrain coordinate system are (Xradar, Yradar).Coordinates of common point in the radar plane coordinate system are (xref, yref), and in this radar imaging plane, coordinates of radar centre are (0, 0).Considering the different axis definition of the two coordinate system, it is necessary to reverse the transverse ordinate of the Gauss coordinates in actual calculation.At this point coordinates common point in the local Gauss plane coordinate system are (Yref, Xref), and coordinates of the centre of the radar are (Yradar, Xradar).The angle of the common point in the two coordinate system can be calculated separately.

arctan( ) arctan( )
where αx-y is the polar angle of common point in radar coordinate system and αE-N is the polar angle of common point in Gauss plane coordinate system.
Then the horizontal axis angle βE-N of the two coordinate axis can be calculated.This rotation determines the rotation relationship between Gauss plane coordinate system and radar plane coordinate system, and it will not change in the whole continuous monitoring process of GB-SAR. = At this stage, distance resolution of GB-SAR system is very small (usually around 0.5m).The corner reflector has strong reflection signal and the window processing in imaging process is difficult to completely eliminate the sidelobe effect.Therefore, the corner reflector in radar image will occupy many resolution units, as shown in Figure 3.The azimuth resolution of GB-SAR increases with the increase of the distance from the centre point of the radar.For example, the azimuth resolution of the IBIS-L system is 4.4mrad, the azimuth resolution at 100m from the radar centre point is 0.44m, and for 1km it is 4.40m.If the corner reflector is closer to the radar system, horizontal angle is more difficult to determine because of the more azimuth resolution units.Usually, only the geometric centres of multiple resolution units can be calculated as the common point target coordinates.However, it is difficult to achieve ideal results using the relationship between multiple resolution units TSNR, and generally only achieve half pixel accuracy.It needs to be adjusted many times in the follow-up coordinate transformation work according to the actual situation.
On the contrary, when the angle reflector is far away from the radar system, the actual resolution of the pixel orientation will be very large, and the side lobe effect will make the location of the corner reflector difficult to determine accurately.The calculation error caused by the horizontal rotation is also inevitable.To sum up, selecting a common point is hard to get a precise horizontal rotation angle, so it needs to select more common points that can be correctly identified in two coordinate systems, but in practice, it is often difficult to get enough ideal public points directly and accurately.

Characteristic straight line method
The common point which is used to calculate horizontal angle can also be regarded as a special line, that is, the straight line connecting common point and radar centre.The direction of this special straight line is entirely determined by the accuracy of the common point selected.However, it is difficult to accurately select the common point in the field measurement, so the use of common point is prone to cause errors.If the imaging of linear structure target in monitoring area can be used, the horizontal rotation angle can be obtained accurately.Actually, in the practice of dam area monitoring, there are a lot of linear structures in radar field, such as Power plant building and stilling pool in water retaining wall.
where αx-y|fitting is azimuth of linear structure in radar coordinate system and αE-N|fitting is azimuth of linear structure in Gauss plane coordinate system.
In the method adopted in this paper, calculation of the horizontal angle and the polar coordinate from 3D terrain coordinate are carried out simultaneously, and the 3D coordinates are converted to the horizontal plane with the radar centre as original point.The feature line is made up of a point set with 3D coordinates, so the azimuth of point set after projection can be accurately determined.The linear building structure in radar image also occupies multiple pixels, and its azimuth can also be accurately calculated.

General introduction
In this experiment an advanced GB-SAR system was used for deformation monitoring of Geheyan dam body and near dam area.Geheyan water conservancy project is located in the main stream of the Yangtze River branch in Changyang County, Hubei, China.It is 207km away from Enshi city and 50km from Gaobazhou hydropower station.Development of this project is mainly for power generation, flood control, shipping guidance and other comprehensive benefits.The project is composed of concrete gravity arch dam, discharge building, water diversion station on the right bank and vertical ship lift on the left bank.The hydropower station is located on the right bank, and the 4 tunnels are connected to 4 water turbine units respectively.Prestressed concrete lining is used in the tunnel.The penstock are excavated to form a high slope of 170m.The experiment mainly monitored the high side slope on the right bank of the dam and the dam area.In order to better interpret image data, the radar coordinates should be converted to the 3D coordinate system of the dam area.Both the high slopes and right bank are placed in the optimal field range of radar, as shown in Figure 5.

Horizontal rotation angle calculation
The transverse and longitudinal axes of the dam coordinate system are opposite to the radar plane coordinate system, terrain coordinates (x, y) need to exchange firstly, obtaining coordinates in local 3D coordinate system.The radar imaging reference plane can be considered to be completely parallel to the horizontal plane because the angle of the radar guide transverse axis and the pitch angle of the radar centre line are within the allowable range.At this time, there is only an acute angle difference between longitudinal axis of local coordinate system and radar plane coordinate system.The edge of the roof of a workshop can be regarded as a space line which becomes a curve with a very small change in curvature after projection calculation into the radar imaging plane.Calculate the curvature changes at both ends of the curve, as in Figure 8.The calculation results show that the curvature value is very small, and the maximum angle change of the angle of the two ends is only 0.0052 degrees, which can be regarded as a straight line.In fact, for the edge of a regular linear structure, we only need to measure the 3D coordinates of two ends using total station, so that the horizontal angle can be calculated with high precision.In figure 9, after extracting some point data from the top of the workshop building, the angle between the fitting line and the xpolar axis is calculated by regression analysis method, which is 0.262230rad.At the same time, extract the edge of the house in GB-SAR image and make fitting calculation.For improving the reliability of the selection, we can calculate the average TSNR of multi GB-SAR images, and use coherence threshold and phase stability threshold method to remove the unstable pixels and false signals.The angle between the fitting line and the X axis is 0.085088 rad, and the final horizontal rotation angle is 0.177142rad.Finally, the coordinate rotation formula is used to calculate the coordinates of the topographic points after clockwise rotation of 0.177142rad.

3D transformation of GB-SAR image coordinates
After completing coordinates projection and horizontal angle rotation, we can find the corresponding points of each GB-SAR image pixel in the radar plane coordinate system as shown in figure 10 (b) and (c).In the experiment, we adopted the minimum distance criterion for 3D point target selection.
Finding the closest point to pixel in radar imaging coordinate system after the projection calculation, the 3D coordinates of this selected point in the polar coordinate system can be used as pixel coordinate.In this case, the elevation extraction of the GB-SAR image pixel is completed.All the related images such as diagram of TSNR, coherent and deformation can be mapped to the local 3D surface models based on the above work.We selected the radar image of the high slope area, and the TSNR threshold was set to 10dB to eliminate the weak part of the signal in the image.As the angle of the high slope in the experiment deviates from radar centre line, the reflected signal in right part of side slope is relatively ideal.GB-SAR image pixels of edge beside the third horse road is clear, which is directly connected to the outer wall of the 4 water diversion tunnels, with relatively high TSNR value, as shown in figure 9.These pixels in figure 9 (a) and the terrain point cloud in (b) are already in the same radar imaging plane coordinate system.By using minimum distance criterion, the matching calculation can be used to determine the 3D terrain coordinates corresponding to each pixel.The terrain surface model of the slope area was established, and the high TSNR pixel was drawn to the surface of the model.In order to remove the pixel that is not in the model surface, the elimination distance was set to 3m before the calculation, and the pixel beyond this distance will not appear on 3D terrain surface model, as shown in figure 11.

CONCLUSION
The corresponding relationship between the GB-SAR image coordinate system and the 3D terrain coordinate system is helpful to the understanding of the image.We used 3D laser scanning technology to realize the coordinate conversion from GB-SAR image plane to 3D terrain space.A straight line fitting method for calculating horizontal angle was proposed in this paper.After projection into a consistent imaging plane, we can calculate horizontal rotation angle by using the characteristics of linear structure in radar image and 3D coordinate system.Experimental results on 3D transformation of radar image coordinates in Geheyan high slope demonstrate that precision of this method is easy to reach resolution level.And it does not need to identify the reflector, or to install corner reflector.But there are also some defects in this method.It can be used only when there are some linear structures within GB-SAR monitoring area, and it has quite a lot restrictions on the spatial attitude of linear structure.It is worth further research how to use any complex structure to reduce the limit conditions and improve the precision of horizontal rotation angle calculation.

Figure 1 .
Figure 1.Relationship between coordinate axes and horizontal angle rotation.

Figure 2 .
Figure 2. Coordinate polarization of point cloud and matching with GB-SAR pixel.

Figure 3 .
Figure 3. Corner reflector and its GB-SAR image pixels.

Figure 4 .
Figure 4. Linear structural targets in the dam area: (a) power plant; (b) stilling pool in water retaining wall.

Figure 5 .
Figure 5. Relation of dam and radar plane coordinate system.

Figure 6 .
Figure 6.Projection calculation of 3D terrain coordinate of slope point cloud to radar imaging plane.

Figure 7 .
Figure 7. 3D coordinates on the edge of power plant and its corresponding GB-SAR image.

Figure 8 .
Figure 8. Projection calculation of point set on workshop roof edge into radar imaging plane: (a) workshop edge point set and projection of workshop edge; (b) angle change caused by curvature change after linear structure projection; (c) curvature change after projection calculation.

Figure 9 .
Figure 9. Straight line fitting calculation: (a) fitting using point set projected to radar imaging plane of linear structure; (b) fitting using pixels of linear structure in radar image.

Figure 10 .
Figure 10.Matching of terrain point and pixel based Minimum distance measurement: (a) TSNR image in high slope area; (b) projection of 3D terrain coordinate of high slope to radar plane.

Figure 11 .
Figure 11.Mapping of TSNR map to digital elevation model.