Automatic Detection and Positioning of Ground Control Points Using TerraSAR-X Multi-Aspect Acquisitions

IEEE TRANSA CTIONS ON GEOSCIENCE AND REMO TE SENSING, IN PRESS 1 Automatic Detection and Positioning of Ground Control Points Using T erraSAR-X Multi-Aspect Acquisitions Sina Montazeri, Christoph Gisinger , Michael Eineder , F ellow , IEEE and Xiao Xiang Zhu, Senior Member , IEEE, Abstract — This is the pre-acceptance version, to read the ﬁnal version please go to IEEE Transactions on Geoscience and Remote Sensing on IEEE XPlore. (DOI: 10.1109/TGRS.2017.2769078) Geodetic stereo Synthetic A perture Radar (SAR) is capable of absolute three-dimensional localization of natural Persistent Scatterer (PS)s which allows for Ground Control Point (GCP) generation using only SAR data. The pr erequisite for the method to achieve high precision results is the correct detection of com- mon scatterers in SAR images acquired from different viewing geometries. In this contribution, we describe thr ee strategies f or automatic detection of identical targets in SAR images of urban areas taken from different orbit tracks. Moreo ver , a complete work-ﬂow for automatic generation of large number of GCPs using SAR data is presented and its applicability is shown by exploiting T erraSAR-X (TS-X) high resolution spotlight images over the city of Oulu, Finland and a test site in Berlin, Germany . Index T erms —Geodetic stereo SAR, Ground Contr ol Point, positioning, synthetic aperture radar , T erraSAR-X. I . I N T RO D U C T I O N Synthetic Aperture Radar (SAR) imaging geodesy and geodetic stereo SAR are relatively new techniques which aim at high precision absolute positioning of point targets in SAR images in two-dimensions (2-D) and three-dimensions (3-D), respectiv ely [1]–[3]. The accuracy of both methods, when coupled with data from T erraSAR-X (TS-X) and T anDEM- X, is in the centimeter regime for targets with accurately known phase centers such as corner reﬂectors [3]. This le vel of accuracy is achiev able due to the precise orbit determination [4] and instrument calibration of the aforementioned satellites followed by a thorough correction scheme which quantiﬁes and removes the most prominent error sources af fecting radar Manuscript receiv ed June 29, 2017; revised September 8, 2017; accepted October 12, 2017. This work was supported in part by the European Research Council through the European Union Horizon 2020 Research And Innovation Program under Grant ERC-2016-StG-714087, in part by the Helmholtz Association through the frame work of the Y oung In vestigators Group SiPEO under Grant VH-NG-1018, and in part by Munich Aerospace e.V . Fakultt fr Luft- und Raumfahrt. (Corresponding author: Xiao Xiang Zhu.) S. Montazeri is with the Remote Sensing T echnology Institute, German Aerospace Center , 82234 W essling, Germany (e-mail: sina.montazeri@dlr .de). C. Gisinger is with the Remote Sensing T echnology Institute, German Aerospace Center , 82234 W essling, Germany . M. Eineder is with the Remote Sensing T echnology Institute, German Aerospace Center, 82234 W essling, Germany , and also with the Chair of Re- mote Sensing T echnology , T echnische Universit ¨ at M ¨ unchen, 80333 Munich, Germany . X. X. Zhu is with the Remote Sensing T echnology Institute, German Aerospace Center , 82234 W essling, Germany , and also with the Signal Processing for Earth Observation (SiPEO), T echnische Uni versit ¨ at M ¨ unchen, 80333 Munich, Germany (e-mail: xiao.zhu@dlr .de). timing measurements. This pav es the way for remotely sensed generation of Ground Control Point (GCP)s using only SAR data. The essential prerequisite for applying the geodetic stereo SAR method is the correct detection of identical scatterers in SAR images acquired from different geometries. In this regard, a target can be visible only from same-heading orbits, i.e. exclusiv ely ascending or descending orbits, or also from cr oss-heading orbits, which include combinations of ascending and descending orbits. Conceptually , a target localized from the latter is fa vorable because of the more robust intersection geometry when compared to the former . This fact is demon- strated in Fig. 1a where the intersection angle occurs at almost 90 ◦ because of the large baseline between the satellites from cross-heading tracks. In Fig. 1b, the target is localized with satellites from same-heading tracks which force the baseline to be smaller and consequently the system of equations to solve for the 3-D coordinates to be less sensitiv e for the perpendicu- lar height direction. Howe ver , the rare occurrence of identical scatterers visible from cross-heading conﬁgurations as well as the challenging task of automatically detecting such targets, either from same- or cross-heading tracks, currently limit the applicability of geodetic stereo SAR for localization of large number of Persistent Scatterer (PS)s. T o ov ercome the limitation to some e xtent, this paper describes an automatic algorithm for detection and absolute positioning of large number of natural PSs in SAR images of urban areas. The candidates are selected from both same- heading and cross-heading geometries based on methods rely- ing on fusion of multitrack Persistent Scatterer Interferometry (PSI) point clouds, correspondence detection with optical data and utilizing vectorized road network data. The candidates are mainly chosen from the same-heading conﬁguration because of the fact that for many PSs the phase centers are assumed to remain unchanged in SAR images. On the other hand, additional candidates are chosen from cross-heading geome- tries, although in a small number , because conceptually they can be localized more precisely compared to the candidates from same-heading geometries. Coupled with the subsequent geodetic Stereo SAR, the proposed processing chain delivers sets of absolutely localized PSs in an in vestigated area. The remainder of the paper is organized as follo ws. Section II re views the theoretical background of the techniques utilized in this study and giv es an overvie w of the recent adv ances and the motiv ation for this work. Section III describes three methods for detecting identical PSs visible in SAR images IEEE TRANSA CTIONS ON GEOSCIENCE AND REMO TE SENSING, IN PRESS 2 (a) Cross-heading (b) Same-heading Fig. 1: Localization of a point tar get (red dot) from (a) cross-heading and (b) same-heading satellite tracks. The satellites are shown by black dots; their trajectories are presented by dashed lines and the baselines are depicted by solid lines between the satellite positions. The black circles are deﬁned by the range-Doppler equations and their intersection leads to the 3-D position of the target. from same- and cross-heading tracks. In Section IV, the complete work-ﬂo w for generating high precision absolute GCPs is explained. In section V, the applicability of the algorithm is demonstrated by exploiting TS-X high resolution spotlight images over the city of Oulu, Finland and a test site in Berlin, Germany , and ﬁnally the conclusions are drawn in section VI. I I . H I G H P R E C I S I O N A B S O L U T E 2 - D A N D 3 - D P O S I T I O N I N G W I T H T S - X At the core of high precision absolute positioning of can- didate GCPs using SAR data are the imaging geodesy and the stereo SAR methods. These methods are described in this section follo wed by the recent adv ances and applications which rely on absolute localization capability of TS-X. It is important to note that the complete explanation of the theory of the methods and their practical implementations are not in the scope of this paper . F or full treatment of these topics, the interested reader is referred to [1]–[7]. A. Backgr ound The SAR imaging geodesy technique aims at achie ving 2-D absolute pix el localization [1]. Based on the SAR measurement principle, a single pixel in a focused complex SAR image, processed to zero-Doppler coordinates, is characterized with two time tags: in the along-track direction, the time relative to the time of the closest approach deﬁnes the azimuth coordinate t az and in the across-track direction, the dif ference in the time trav el of the transmitted and the received chirp at t az describes the range coordinate τ rg [8]. If we measure the radar timing coordinates ( t az , τ rg ) for a point target located within the mentioned pixel, the following equations hold: τ rg = 2 R c + δ τ S D + δ τ O + δ τ F + δ τ I + δ τ T + δ τ G , (1) t az = t + δ t S D + δ t O + δ t F + δ t G , (2) where R is the geometric distance from the satellite to the center of the pixel in meters and c is the speed of light in vacuum in m/s while all the other terms are expressed in seconds; t is the raw acquisition time, δ τ S D and δ t S D are delays caused by satellite dynamics and electronics, δ τ O and δ t O are the orbit inaccuracies, δ τ F and δ t F are the feature localization error, δ τ G and δ t G include the geodynamic effects all on range and azimuth timings, respectiv ely while δ τ I and δ τ T are the ionospheric and the tropospheric delays considered only for range timings. The magnitude of the indi vidual effects can be scaled to units of length by multiplying the range error terms with c 2 and the azimuth error terms with the platform’ s velocity . The outcomes v ary from a couple of centimeters for the ionospheric effect, if the satellite operates in X-band, followed by decimeter regimes for satellite electronic delays and geodynamic ef fects for both components, to up to four meters for the tropospheric effect depending on the a verage incidence angle of the acquired TS-X images. Imaging geodesy corrects for all the error terms in (1) and (2) thus obtaining absolute range and azimuth timings. In this regard, the technique reduces the satellite dynamics effects by av oiding the stop-go approximation in the T erraSAR-X multimode SAR processor (TMSP) and by taking into account the non-zero duration of the pulses and the internal delay caused by the instrument cables [9]. The propagation errors are estimated based on the path delays deri ved from the near- by Global Navigation Satellite System (GNSS) stations or 3-D integration through weather models follo wed by appropriate mapping functions [2], [3], [10]. For the geodynamic effects such as solid earth tides, plate tectonics, ocean loading and atmospheric loading, which change the position of a target on the ground, the corrections are applied based on models issued by the International Earth Rotation and Reference Systems Service (IERS) [11]. T aking into account all the mentioned factors, SAR imaging geodesy is currently capable of providing range and azimuth measurements with 1.16 cm and 1.85 cm standard deviations, respectively [12]. If a target is visible in SAR images acquired from two or more dif ferent viewing geometries, then stereo SAR retrieves the 3-D position of the target by combining the extracted timing information of the tar get from each SAR image. Furthermore, if the timing coordinates hav e been a-priori corrected for the error sources expressed in (1) and (2), the method is called geodetic stereo SAR which allows for absolute 3-D localization [3]. The relation between the 2-D radar time coordinates of a speciﬁc target in the SAR image x T = ( t az , τ rg ) and its corresponding 3-D coordinates on the ground X T = ( X , Y , Z ) is deﬁned by the range-Doppler equation system [8]: | X S − X T | − c · τ rg = 0 , (3) ˙ X S ( X T − X S ) | ˙ X S || X S − X T | = 0 , (4) with X S and ˙ X S being the position and velocity v ector of the satellite relati ve to t az , and τ rg being the calibrated two-way traveled time from the satellite to the target. t az is implicitly included in (4) relating the state-vector of the satellite to the time of the acquisition via a polynomial model [3]. Equation (3) deﬁnes a sphere centered on X S which reduces to a circle perpendicular to the satellite trajectory when coupled with the zero-Doppler plane described in (4). Therefore X T can be retriev ed by including another set of IEEE TRANSA CTIONS ON GEOSCIENCE AND REMO TE SENSING, IN PRESS 3 timing observations from a different satellite position which ev aluates the intersection point of the two circles, see Fig. 1. The estimation of the coordinates is carried out by means of least squares plus stochastic modeling of range and azimuth using the v ariance component estimation (VCE) [3]. Precision of the estimated 3-D coordinates depends on the Signal- to-Clutter-Ratio (SCR) of the target, the precision of the external radar timing corrections, the separation in the viewing geometries and the number of acquisitions. Geodetic stereo SAR has been proven to be able to localize corner reﬂectors with 3-D precision better than 4 cm and an absolute accuracy of 2-3 cm when compared to independently surve yed reference positions [3]. B. Recent advances and motivation In our previous research, geodetic stereo SAR has been also applied to small number of natural PSs in urban areas where it could localize targets in 3-D with a precision better than one decimeter using TS-X high resolution spotlight products [3]. The PSs were manually extracted from SAR images and originated from building facades for candidates visible in same-heading tracks or from the base of street lights for candidates visible from cross-heading tracks. In [13], the ﬁrst attempt for automatic timing extraction and matching of limited number of PSs originated from a building facade visible in TS-X images from two same-heading tracks was reported. In the study , the geodetic stereo SAR method was extended to include the secular movement of the PSs as well as their 3-D absolute positions. The av eraged 3- D precision was reported to be below one decimeter with encouraging results for estimating the plate tectonics using SAR data. In [14], the concepts of imaging geodesy and stereo SAR were used to transform the relative estimates of SAR tomog- raphy into absolute 3-D point clouds by absolutely localizing the manually extracted reference point. The method, termed geodetic SAR tomography , allows for generation of dense point clouds with an absolute localization accuracy in the order of 20 cm and is the basis for geodetic fusion of multi- aspect Interferometric SAR (InSAR) point clouds. The latter enables the decomposition of deformation estimates from SAR tomography into highly detailed 3-D displacement maps [15]. Automatic extraction of GCPs from SAR products have been carried out for TS-X and COSMO-SkyMed in [16] and [17], respectiv ely . Both methods focus on detection of stereo candidates that presumably originate from street lights or trafﬁc signs and are visible in SAR images as bright isolated points. Therefore, the majority of PSs in urban areas which stem from building corners or facades are not considered as potential candidates in these methods due to complex radar reﬂection properties in such scenarios. Furthermore, in [16] the 3-D positioning is done only with same-heading geometry conﬁgurations and therefore the error ellipsoid of the scatterers’ coordinates is highly ske wed in the cross-range direction [3], [18]. Ne vertheless, the retrie ved 3-D coordinates of several candidates were compared to their true positions observed with GNSS which conﬁrmed an absolute accuracy better than 20 cm in each coordinate component [16]. Based on the above-mentioned studies, the motiv ation for carrying out this work is four -fold: • The manual detection, extraction and matching of PS candidates from SAR images acquired from different geometries is cumbersome and should be replaced by an automatic process. • The algorithm should be able to automatically detect and match identical PSs visible from cross-heading geome- tries in order to boost the precision of the retriev ed 3-D coordinates. • The total number of high quality PS candidates to be localized as GCPs should be as large as possible. This indicates that the majority of PSs in urban areas which stem from buildings should also be considered as can- didates for 3-D absolute localization from same-heading tracks. • The distribution of the GCPs should be as homogeneous as possible in the entire in vestigated area. This allows for generation of an absolute reference network to be integrated into relativ e InSAR techniques. I I I . D E T E C T I O N O F I D E N T I C A L P S S I N M U LT I - A S P E C T S A R I M A G ES Detection of identical PSs from SAR images acquired with different vie wing geometries is a challenging task. This is be- cause of the geometrical distortions of SAR images due to the oblique viewing geometry and less importantly the presence of speckle. Moreover , in urban areas captured by SAR sensors, which is the focus of this study , the backscattering mechanism is highly complex because of existence of se veral phase centers close to each other . Therefore, identical PS matching becomes ev en more difﬁcult for multi-aspect SAR images of urban areas. In recent years, there have been several studies which ex- plored the possibility to match features between SAR images. In [19], the capabilities of the con ventional scale in variant feature transform (SIFT) algorithm [20], which is commonly used for feature extraction and matching between optical images, were extended to be suitable for SAR images. In [21], the SAR-SIFT algorithm has been proposed which focuses on the ef ﬁcient extraction of local descriptors from SAR images by modifying the SIFT algorithm to take into account the statistical properties of speckle. Howe ver , both of the aforementioned methods are applicable only to SAR images taken from same-heading orbits with small difference in the respectiv e incidence angles. Speciﬁcally for the task of auto- matic 3-D positioning, in [17] the authors have proposed to identify identical scatterers based on detection of local features using the Harris corners. This is follo wed by constraining the search space by geocoding the local features, using an e xternal digital elev ation model and orbit information, and eventually using SIFT for the feature matching. Although the method is promising in terms of detection and positioning of tar gets ev en from cross-heading tracks, it only works on isolated PSs. In the follo wing we describe in detail the three strategies we apply for detection of identical PSs in SAR images acquired from same- and cross-heading orbits. The methods do not IEEE TRANSA CTIONS ON GEOSCIENCE AND REMO TE SENSING, IN PRESS 4 tackle the detection problem directly within the SAR images but instead rely on external geospatial data and on limiting the search space on geo-referenced positions of the PSs. A. Multitrack PSI point cloud fusion In [22], a method for geometrical fusion of multitrack PSI point clouds stemming from PSI has been proposed. The fusion algorithm, which is based on the geocoded PSI point cloud solutions of each geometry as well as information on heading and looking angle of the satellites, consists of three major parts, namely: 1) generation of initial point correspon- dences, 2) restricted least squares adjustment to minimize the distance between assumed identical points visible from different vie wing geometries, and 3) adding a range-dependent shift to all PSs using the result of the previous step for the ﬁnal registration. A summary of the method is described in the following. F or a detailed description of the algorithm the reader is referred to [22], [23]. Since we are interested in the detection of large number of point correspondences, only the ﬁrst part of the algorithm is relev ant. This coarse registration is performed based on cross-correlation of a subset of geocoded PS point clouds from different geometries, after projection on a regular grid, in the xy-, xz- and yz-planes. The subset is chosen based on precision of height update estimates av ailable for each PS after carrying out PSI [23], [24]. The resulting horizontal and vertical shifts from the mentioned cross-correlation procedure are applied to the PSs of one point cloud to align them with the PSs of the other point cloud. The coarse shifts are further reﬁned prior to the selection of corresponding PS pairs. The reﬁnement is carried out inside a small neighborhood around each PS which includes the PSs from the other point cloud and tends to accomplish it by performing a statistical search to ﬁnd the best ﬁt between both 3-D point clouds [23]. The reﬁned shift is applied to the PS point cloud of one acquisition geometry and a one-by-one PS correspondence is detected in the other point cloud. At the ﬁnal step, the 3-D coordinates of the geocoded PSs hav e to be projected on the SAR images of each orbit track, a process called radar-coding. Since the matching of the PSI results is performed on coordinates in the Universal T ransverse Mercator (UTM) map projection, the coordinates are ﬁrst con verted to the Cartesian geocentric system as ( X, Y , Z ) i for the ith PS. Subsequently , the range-Doppler equations described in (3) and (4) are in verted to obtain the azimuth and range timing coordinates ( t az , τ rg ) i which can be easily expressed in pixels in the radar coordinate system ( L, P ) i by knowledge of the range sampling frequency , pulse repetition frequency , ﬁrst sampled azimuth time and the ﬁrst sampled range time for each acquisition. The latter information is stated in the product annotation ﬁles accompanied by the TS-X image products [25]. For the same-heading tracks, this method typically generates 200 to 2000 of point correspondences per km 2 depending on how densely constructed is a city which directly affects the total number of PSs in each point cloud. B. T emplate matching on optical data Giv en the av ailability of suitable remotely sensed optical data, one can detect candidate objects from optical images which are probable to be observed in SAR images from different viewing geometries. In urban areas, scatterers which are good candidates to be visible from both same-heading and cross-heading tracks usually originate from lamp poles or other cylindrical objects that are vertically oriented tow ards the sensor . Therefore, the basic idea when using optical data for the aid of GCP identiﬁcation is to detect lamp poles and match the detected objects to the corresponding bright points in SAR images. The method identiﬁes lamp poles based on their distincti ve shadows in optical images using a template matching scheme [26]. Prior to extracting the template, common pre-processing steps such as noise ﬁltering and histogram equalization are carried out on the optical image. Additionally , in order to make the shadows of lamp poles more prominent, a simple sharpening procedure is carried out as follo ws: I = I o + a I m , (5) where I is the sharpened image, I o is the pre-processed original image, a is the scalar sharpening factor and I m is the un-sharp mask. I m is calculated as the difference between I o and its blurred version. Higher values of factor a , means higher level of sharpening. The process expressed in (5) is called high-boost ﬁltering [27]. After the sharpening, the template is extracted based on the shadow of an arbitrary lamp pole visible in the optical image. The template is then correlated with the reference image to calculate the following similarity measure for each pixel ( u, v ) in the reference image [28]: ρ ( u, v ) = Σ x,y  I ( x, y ) − ¯ I u,v  T ( x − u, y − v ) − ¯ T  q Σ x,y  I ( x, y ) − ¯ I u,v  2 Σ x,y  T ( x − u, y − v ) − ¯ T  2 , (6) where I ( x, y ) and T ( x, y ) denote pixel values of the refer- ence and the template image at ( x, y ) , respectiv ely , and Σ x,y stands for Σ N 1 x =1 Σ N 2 y =1 with N 1 × N 2 being the size of the template. Furthermore, ¯ I u,v and ¯ T denote the mean intensity values of the original image and the template, respecti vely . Equation (6) allows for calculation of the Normalized Cross- Correlation (NCC) value ρ ( u, v ) which leads to the detec- tion of the template location in the reference image after proper thresholding. It is important to note that due to the normalization carried out in the denominator of (6), ρ ( u, v ) is independent of changes in brightness or contrast of the image and therefore improves the result of template matching. After detection of pixels which belong to the shadows of lamp poles, the result is geo-referenced in the UTM coordinate system. Since for each lamp pole in the optical image more than one pixel exists which represent the object, a subsequent clustering is performed. The clustering is carried out non- parametrically using the mean shift concept [29]: IEEE TRANSA CTIONS ON GEOSCIENCE AND REMO TE SENSING, IN PRESS 5 M ( p i ) = Σ n j =1 p j g  k p i − p j h k 2  Σ n j =1 g  k p i − p j h k 2  − p i , (7) where p i denotes a 3-D point for which the shift vector M ( p i ) is calculated. p j represents the points in a neigh- borhood of p i , g is a kernel function with the bandwidth h and k·k is the Euclidean distance operator . The main idea of the algorithm is to shift each point in a small neighborhood tow ards its weighted mean value and thus representing each cluster by its centroid [29]. The process is carried out itera- tiv ely until the length of M ( p i ) becomes equal or close to zero. For our application, since in an y case there will be a mismatch between the detected points on optical data and the corresponding bright points in the SAR image, utilizing a ﬂat kernel in equation (7) suf ﬁces. This means the algorithm is simpliﬁed by calculating the sample mean in a speciﬁed radius of p i and shifting the desired point towards the estimated center . In the next step, the clustered points with UTM coordinates should be radar -coded to all the available SAR images. As it was mentioned earlier, the positions of the detected lamp poles on the optical data and the bright PSs in the SAR image will most probably not coincide after radar -coding. This can be explained by height uncertainties of the geo- referenced optical data and the fact that the data may not be perfectly orthorectiﬁed. Therefore in the ﬁnal step of the algorithm, the detected lamp poles are registered on the corresponding bright dots in the SAR image by employing the iterativ e closest point (ICP) algorithm [30]. T o this end, binary masks are generated based on thresholding on the bright points on the SAR image and keeping only the detected lamp poles from the optical data. The ICP algorithm then ﬁnds for each indi vidual point its closest point in the corresponding point set. It iterati vely estimates the transformation parameters (translation and rotation) to minimize the mean squared error between the two point sets and ﬁnally registers one point cloud onto the other point cloud with the reﬁned transformation parameters. It is note worthy that the positioning accuracy of the utilized optical image does not necessarily have to be high. A hori- zontal positioning accuracy in the order of couple of meters and an approximate knowledge of height based on freely av ailable sources usually sufﬁce for the procedure described in this subsection. If a mismatch occurs due to the low positioning accuracy , this will be compensated by the ﬁnal step of the algorithm with applying the ICP . Howev er, the spatial resolution of the optical data should be strictly high, 10 cm or better, in order to be able to accurately detect the shadows of the lamp poles. C. V ector r oad network data In urban areas, the cylindrical objects of our interest (lamp poles, road signs, traf ﬁc lights, etc.) are typically located along the roads. Therefore, with the av ailability of geospatial road data, either obtained from OpenStreetMap or country-speciﬁc geoportals, and the projection of such maps on SAR images, one can search for bright points in the neighborhood of the road data points. The method is applied to co-registered stacks of SAR images. If SAR stacks from multiple viewing geometries are av ailable, ﬁrst the road data, which is usually delivered in the UTM coordinate system, is radar-coded based on the master orbit information of each stack. It is important to note that if the road data do not have any information about the ellipsoidal height, then for the radar-coding a constant height value based on prior knowledge can be chosen for all the road data points. Furthermore, the data with horizontal positioning accuracy in the order of couple of meters will suf ﬁce for the PS matching procedure as the PS correspondences are detected on a neighborhood-analysis basis, as is depicted in the following, which does not require the exact position of the road data point. After radar-coding, a circular neighborhood is considered around each road data point. The radius of the circle depends on the typical width of streets and highways. Subsequently , for each pixel within the neighborhood the Amplitude Dispersion Index (ADI) is ev aluated [31]: D a ≈ σ a ¯ a , (8) where σ a and ¯ a are the temporal standard de viation and the temporal mean of calibrated amplitude values of the pixel, respectiv ely , and D a approximates the phase dispersion. The pixel with the lo west value of D a , i.e. the one with the highest phase stability is chosen as potential PS candidate. This process is carried out for all of the av ailable road data points. Since at this point, it is possible that many false pixels with relativ ely low D a values in the neighborhoods are categorized as potential GCP candidates, a further thresholding on D a is performed in SAR images from all av ailable viewing geometries. This operation, in addition to constraining the approximate elev ation of the PS candidates to be close to the ground, causes a dramatic decrease in the total number of detected candidates b ut improv es the accuracy of the detection. Finally , the presumable identical PSs in all av ailable geome- tries are geocoded using the respectiv e master orbit informa- tion. In the geocoded results the PSs which are close enough, in terms of coordinate differences, are selected as the ﬁnal GCP candidates. I V . A U T O M A T I C G C P G E N E R A T I O N : T H E P R O C E S S I N G C H A I N The processing chain for automatic detection and position- ing of GCPs includes a set of procedures which starts from single-look slant range comple x SAR images and their corre- sponding product annotation ﬁles to absolute 3-D coordinates of the chosen GCPs. The ﬂowchart of the algorithm is sho wn in Fig. 2. It consists of the following major processes which hav e to be carried out in the stated order: 1) identiﬁcation of identical scatterers visible in multi- aspect SAR images. 2) precise extraction of scatterers’ azimuth and range po- sitions from SAR images at sub-pixel level. IEEE TRANSA CTIONS ON GEOSCIENCE AND REMO TE SENSING, IN PRESS 6 S i n g l e - L o o k S l a n t r a n g e C o m p l e x ( S S C ) i m a g e s L 1 B X M L p r o d u c t a n n o t a t i o n f i l e S a m e - h e a d i n g G e o m e t r y ? Y e s P S I p r o c e s s i n g M a s t e r s e l e c t i o n a n d c o - r e g i s t r a t i o n I n t e r f e r o g r a m g e n e r a t i o n a n d P S d e t e c t i o n P S r e f e r e n c e n e t w o r k c o m p u t a t i o n A P S e s t i m a t i o n a n d r e m o v a l N o O p t i c a l d a t a a v a i l a b l e ? N o Y e s L a m p p o l e d e t e c t i o n C o n t r a s t s t r e t c h i n g a n d s h a r p e n i n g T e m p l a t e m a t c h i n g b y N C C e s t i m a t i o n G e o c o d i n g i n U T M M e a n s h i f t c l u s t e r i n g P S m a t c h i n g o n S A R i m a g e s R a d a r - c o d i n g B i n a r y m a s k s b a s e d o n o p t i c a l a n d S A R d a t a P o i n t s e t r e g i s t r a t i o n b a s e d o n I C P H i g h r e s o l u t i o n o p t i c a l d a t a R o a d n e t w o r k d a t a i n v e c t o r f o r m a t B r i g h t P S c a n d i d a t e d e t e c t i o n a n d m a t c h i n g M a s t e r s e l e c t i o n a n d c o - r e g i s t r a t i o n R a d a r - c o d i n g o f r o a d d a t a A D I e s t i m a t i o n i n r o a d p o i n t n e i g h b o r h o o d s T h r e s h o l d i n g o n A D I a n d e l e v a t i o n o f P S s G e o c o d i n g a n d o n e - b y - o n e P S m a t c h i n g F i n a l h e i g h t a n d d e f o r m a t i o n e s t i m a t i o n P S c o r r e s p o n d e n c e d e t e c t i o n G e o c o d i n g P S s o n 2 - D g r i d s C o a r s e c r o s s - c o r r e l a t i o n L o c a l r e f i n e m e n t o f t h e r o u g h s h i f t s O n e - b y - o n e P S m a t c h i n g R a d a r - c o d i n g P r e c i s e t i m i n g e x t r a c t i o n o f P S f r o m S S C i m a g e s C o r r e c t i o n o f P S t i m i n g s i n t h e s t a c k A b s o l u t e 3 - D l o c a l i z a t i o n o f P S O v e r s a m p l e b y a f a c t o r o f 3 2 i n r a n g e a n d a z i m u t h 2 - D p a r a b o l o i d k e r n e l i n t e r p o l a t i o n o f p e a k P S v i s i b i l i t y c h e c k a n d i n i t i a l o u t l i e r r e m o v a l S C R a n d p h a s e n o i s e e s t i m a t i o n R o b u s t o u t l i e r r e m o v a l f r o m p h a s e n o i s e t i m e s e r i e s R e m o v e d a t a t a k e s w i t h a b o v e 0 . 5 r a d p h a s e n o i s e C o m p u t a t i o n o f e x t e r n a l c o r r e c t i o n s S t e r e o S A R O b s e r v a t i o n r e s i d u a l a n a l y s i s a n d o u t l i e r r e m o v a l Fig. 2: The ﬂo wchart of the automatic GCP generation algorithm. The input data is the Single-Look Slant range Complex (SLC) SAR images and their accompanying L1B product annotation ﬁles. The identical PS detection part is drawn inside the gray dashed rectangle. In case of PS detection from cross-heading viewing geometries, the algorithm recei ves as input the optical image or road network data of the scene upon av ailability . The processes are depicted with big blue rectangles while the sub-processes are shown with small gray rectangles. The name of each process step is written in bold letters. The last blue rectangle includes the processes and sub-processes that all the detected candidates should go through independent of the detection methodology . 3) scatterer visibility check and initial remov al of outliers from time series of phase noise. 4) correction of radar timings for all the perturbing signals. 5) estimation of the 3-D absolute coordinates of the scat- terers. This section discusses each step. It is important to note that since some of the steps are well-established and well- documented techniques in the SAR community , all the details IEEE TRANSA CTIONS ON GEOSCIENCE AND REMO TE SENSING, IN PRESS 7 will not be repeated here. Instead, for these processes, the relev ant references are pro vided. A. Identiﬁcation of identical PS This part has been already cov ered in Section III and is shown as the processes inside the gray dashed rectangle in the ﬂowchart of Fig. 2. F or detection of PS candidates from same- heading orbits, the PSI processing is carried out following the guidelines from [24], [31]–[33] and the identical PS matching is done using the PSI multitrack fusion algorithm described in [22], [23]. In the case of localization from cross-heading tracks, the PS candidate selection is carried out using the methods based on the optical data or on the road network data. Regardless of the detection method, the output of this step is the approximate radar coordinates of identical PSs in terms of lines and samples in all non-coregistered images of different orbit tracks. B. Precise extraction of PS timings The rough radar coordinates of the PSs from the previous step should be reﬁned to sub-pixel le vel in order to extract the timings precisely . T o this end, a process called Point T arget Analysis (PT A) is carried out [8]. In each image of the scene, a 32 × 32 window centered on the PS is extracted. In both range and azimuth direction, an o versampling by a factor of 32 is performed and the integer peak position of the PS response is measured. Subsequently in a 3 × 3 windo w centered on the peak position, a paraboloid interpolation is performed to reﬁne the values around the maximum. This method allo ws the retriev al of the peak position of the PS with a sensitivity better than 1 1000 of a pixel in each dimension [5]. These v alues are then con verted to radar timings based on the product annotation ﬁles [25]. Based on the result of the PT A for each PS, the reﬁned peak power P peak and the clutter power P clutter can be computed. These values, if expressed in dB , are related to the SCR of the tar get as: SCR = 10  P peak − P clutter 10  (9) which is expressed as a digital number . C. PS visibility chec k and initial outlier r emoval The SCR of potential PS candidates should be high enough in the stack of SAR images so that the PS can be localized with high precision. Therefore, we analyze the time series of phase noise values σ φ of the PSs to exclude potential outliers and check if the scatterer is visible in one data take or not. The σ φ of a PS in acquisition i is related to the SCR of the target as [34]: σ φ i ≈ 1 √ 2 SCR i (10) which is expressed in radians. The values of σ φ , for a spe- ciﬁc PS in all data takes, are non-negati ve and follow a right- ske wed distribution. Removal of outliers based on statistical measures such as mean or median is not recommended since many regular v alues can be categorized as outliers. Therefore, we use a method called adjusted boxplot which allo ws for robust elimination of outliers in uni variate sk e wed distrib utions [35]. The main idea of the adjusted boxplot is to modify the original boxplot method, described in [36], to include infor- mation about the ske wness of the data. Therefore, instead of classifying an observation as outlier if it lies outside of the interval deﬁned by the boxplot method [36]:  Q 1 − 1 . 5 I QR ; Q 3 + 1 . 5 I QR  , (11) the adjusted boxplot method declares an observation as outlier if its value exceeds the follo wing interv al [35]:  Q 1 − 1 . 5 e ( − 4 M C ) I QR ; Q 3 + 1 . 5 e (3 M C ) I QR  . (12) In (11), Q 1 and Q 3 are the ﬁrst and the third quartiles of the data, respectiv ely and I QR = Q 3 − Q 1 denotes the interquartile range. In (12), M C is the medcouple, a robust measure of the skewness of a univ ariate sample which for right-ske wed distributions is always non-negati ve [37]. The exponential functions in (12) which depend on the M C as well as the included coefﬁcients are chosen experimentally based on some well-kno wn ske wed distrib utions. For more details on the theory and implementation of the adjusted boxplot method, the reader is referred to [35]. After the automatic identiﬁcation and e xclusion of σ φ values which do not lie within the interval of (12), the remaining time series is analyzed to remove the data takes in which the speciﬁc PS is not visible. This is done by removing all σ φ values which are abov e 0.5 radians ( ≈ 30 ◦ ) as is stated in [24]. D. Correction of PS timings in a stack The correction of the extracted PS timings is performed using the imaging geodesy technique [1] which was brieﬂy introduced in Subsection II-A. It is worth mentioning that the tropospheric and ionospheric effects are corrected based on global numerical weather models and global ionospheric maps, respectiv ely , if local GNSS receiv ers are not av ailable in the vicinity of the inv estigated area. Along with these corrections comes the corresponding geometrical calibration of the SAR sensor in range and azimuth which ensures centimeter localization accuracy . The calibration is based on corner reﬂectors with known reference coordinates [1]. The output of this part is the absolute 2-D radar timing coordinates. E. Absolute 3-D localization of PS At the ﬁnal step in the processing chain, the corrected range and azimuth timings from the entire multi-aspect set of SAR images are combined to retriev e the absolute 3-D position of the PS with the stereo SAR method described in Subsection II-A. Apart from the 3-D position of the target, stereo SAR reports on the standard deviation of each coordi- nate component ( S X , S Y , S Z ) as the by-product of the least squares adjustment. Furthermore, observation quality of each IEEE TRANSA CTIONS ON GEOSCIENCE AND REMO TE SENSING, IN PRESS 8 PS i.e. the azimuth and range standard deviations ( S az , S rg ), retriev ed from applying VCE to residuals, are deli vered. It is important to note that the VCE is carried out indi vidually for each geometry which allows to judge the consistency of the observed geometries with respect to the underlying assumption that the intersection occurs at a common PS. The residuals of the adjusted range and azimuth obser - vations are the basis for the elimination of outliers after stereo SAR processing. Therefore, the processing is carried out repeatedly , where ﬁrst the initial 3-D coordinates are estimated using the provided input timings. The range and azimuth residuals are analyzed to exclude observations with residual values larger than 0.6 m in range or larger than 1.1 m in azimuth. The upper bounds correspond to the nominal spatial resolution of TS-X high resolution spotlight products used in this study . Then stereo SAR is performed again with the cleaned observations. This time, observations which show residuals larger than two times the S az and the S rg are removed. Additionally , to remov e the PSs which are not considered ideal for stereo SAR due to wrong correspondence matching caused by several scatterers being too close, a third step of data cleaning is performed. PSs having an S az higher than 20 cm in any of their azimuth geometries are removed based on the assumption that the discrepancy should not exceed the typical size of the PS object, for instance, a lamp pole. The estimated variance-co v ariance matrix of the 3-D po- sition of the target is further used for error analysis. The matrix giv es important information about the stability of the coordinates results and is affected by the factors stated in Subsection II-A. V . E X P E R I M E N TA L R E S U L T S In this section, the work-ﬂo w described in Section IV is applied on real data to produce remotely sensed SAR-based GCPs. In Subsection V -A, the results are reported for a small test site in Berlin where the detection of GCP candidates are carried out using cross-heading geometries. In Subsection V -B, the processing results are shown for the entire city of Oulu, Finland, where the detection and positioning are performed on PS candidates detected from both same- and cross-heading orbit tracks using the methods described in Subsections III-A and III-C, respectively . A. Berlin The ﬁrst test site includes a small area close to the Berlin central railway station. The SAR data set, 214 images in total, consists of two stacks of TS-X very high resolution spotlight images acquired with a range bandwidth of 300 MHz. The images cov er a period of approximately six years from April 2010 to March 2015. In terms of viewing geometry , one stack was acquired from a descending orbit with images recorded at 05:20 coordinated universal time (UTC), and one stack was acquired from an ascending track with images recorded at 16:50 UTC. The acquisition parameters of each stack are summarized in T ab . I. T ABLE I: A C Q U I S IT I O N PAR A M E T E R S O F S TAC K S O F S A R I M AG E S I N B E R L I N Beam Nr . Incidence angle (degree) Heading angle (degree) T rack type Nr . of images 57 41.9 350.3 Asc 107 42 36.1 190.6 Dsc 107 For the selected test site, an aerial optical image with ground spacing of 7 cm is also av ailable. The optical image is orthorectiﬁed and was used in a stereo matching process to produce a digital surface model with decimeter accuracy [38]. The optical image of the test site and the corresponding SAR images are shown in Fig. 3. The GCP candidate selection is carried out based on the optical data which includes the detection of lamp poles and their projection onto cross-heading SAR images (see Subsec- tion III-B). The individual steps of this process, applied on the Berlin test site, are shown in Fig. 4. After the extraction of the template (see Fig. 4a), the NCC map is calculated and pixels with values higher than 0.6 are classiﬁed as parts of shadows of lamp poles as illustrated in Fig. 4b. It is important to note that the exact tuning of the threshold is not necessary as long as the value is chosen low enough. If the threshold is strictly chosen as a high value, although we are selecting the most similar pix els to the template, we may ignore all the other candidates which sho w less similarity to the template b ut might hav e been potential candidates for stereo SAR processing and 3-D localization. Therefore, in our processing chain the default value is set to 0.6, which is slightly higher than the half of the NCC range [0,1], to guarantee a certain degree of similarity while selecting a large number of pixels. In the Berlin case study , this leads to the selection of 2030 pixels which are further clustered to represent the 44 detected objects in the UTM coordinate system (see Fig. 4c). The objects include the lamp poles along the bridge and at the street perpendicular to the bridge as well as ﬂag poles at the top left of the optical image close to the Berlin railway station. After radar-coding of the results onto the entire ascending and descending SAR images, the mismatch between the projected points and the actual bright points in the SAR images, depicted in Fig. 4d and Fig. 4e, is resolved using ICP . The detection outcome is marked with green circles in Fig. 4f and Fig. 4g. It is seen from the results that not all of the av ailable street lights can be detected using the mentioned strategy as some can be occluded by cars or the object’ s shadow is not distinctiv e enough to match the extracted template. Howe ver , this is of low importance in our application since in such a small test site, with an area less than two km 2 , two or three GCPs will certainly sufﬁce. Moreover , if the method detects wrong candidates which do not fall in the category of PSs, the subsequent PT A and phase noise analysis will discard these points. The precise radar timings of the PSs are extracted using PT A (see Subsection IV -B). Subsequently for each PS candidate in each data take, SCR and σ φ values are ev aluated using (9) and (10), respectively . After excluding potential outliers IEEE TRANSA CTIONS ON GEOSCIENCE AND REMO TE SENSING, IN PRESS 9 Fig. 3: The optical image (left) and the SAR amplitude images (right) of the test area in Berlin. The contrast of the optical image has been adjusted to illustrate the shadows of the lamp poles prominently . with the adjusted boxplot method, the data takes in which the scatterer is not visible are discarded by thresholding on the σ φ values (see Subsection IV -C). At this point, one might argue that analyzing the time series of the remaining σ φ values can already gi ve a hint if the scatterer is a suitable candidate for positioning or not. This statement is partially true since the mentioned analysis is only useful to separate the time- coherent scatterers from the non-coherent ones. It is possible that sev eral scatterers, located close to each other , are mapped as one bright point which results in low σ φ value b ut of course not a suitable candidate for stereo SAR. These candidates are usually discarded in the stereo SAR processing due to large S az values which indicate that not the same object has been detected from multiple viewing geometries. Therefore, in our processing chain, the detected candidates are not entirely remov ed based on their average SCR and they are all passed to the ﬁnal stereo SAR processing. As for the corrections, the ionospheric delay is estimated using global ionospheric maps. The tropospheric delay is estimated using the zenith path delay information of the closest permanent GNSS station in Potsdam which is situated approximately 35 km away from the test site. The positioning of the 44 PSs is carried out using stereo SAR (see Subsection II-A) and is follo wed by outlier elim- ination according to criteria described in Subsection IV -E. Starting from the corrected input timings, the ﬁrst solution is analyzed for gross outliers in the observ ations exceeding the resolution of the underlying TS-X high-resolution spotlight product. Applying these thresholds to the residuals of the adjusted observations reduces the number of solvable PSs from 44 to 42, because in the case of an obvious mismatch, all the observ ations of one geometry are removed. Moreover , the total number of observ ations is reduced by 25 % but this strongly varies across the individual PSs. After re-computation and application of the 2 σ test using the estimated standard deviations from the VCE ( S az , S rg ), the number of solvable PSs remains 42 and the total amount of observ ations is reduced by another 8 % . At this stage, the data is fairly cleaned at the observation lev el regarding the individual range and azimuth geometries, but their consistency has not been considered so far . Looking at the observ ation residuals of the two PSs displayed in Fig. 5 rev eals that there are cases for which the azimuth of one geometry is clearly biased because we try to combine data from different phase centers. In the ideal case, the algorithm yields a coordinate solution for which all sets of observations (two sets of azimuth, two sets of range) can fulﬁll the range-Doppler positioning model of (3) and (4). For a mismatch or a spatial separation of the phase centers, the usually more precise range observations dominate the solution, and only one of the azimuth data sets may ﬁt the estimated coordinates without a bias, but not the second set of azimuth data. Such a situation is illustrated by PS 41 (see Fig. 5a), where the 36 ◦ azimuth displays a prominent bias of about -30 cm. T o a certain degree this must be accepted since we can not expect ideal multi-directional PSs, e.g. the lamp-poles hav e a certain diameter . In this study , as was mentioned in Subsection IV -E we deﬁne an empirical limit of 20 cm of what we consider acceptable which remov es candidates like the PS 41 during the ﬁnal processing step. Therefore, the residual results of scatterers like PS 43 (see Fig. 5b) may be seen as a best case scenario. The remaining difference in quality between range and azimuth is due to the non-square product resolution, i.e. the TS-X spotlight SLC data has a resolution of 0.6 m × 1.1 m in range and azimuth, respectiv ely [25]. For the aforementioned reasons, only 9 PSs remain which we consider as good GCP candidates. The bar graphs of Fig. IEEE TRANSA CTIONS ON GEOSCIENCE AND REMO TE SENSING, IN PRESS 10 (a) (b) Easting Northing (c) (d) (e) (f) (g) Fig. 4: Demonstration of PS correspondence detection in Berlin based on high resolution optical data. (a) shows the pre-processed optical image after negati ve intensity transformation and the extracted template. (b) is the calculated NCC map after correlating the extracted lamp pole template with the reference image in which the detected objects are marked by yellow rectangles. (c) shows the 44 detected objects in the UTM coordinate system after clustering. In (d) and (e), the radar -coded results are depicted by yellow circles which sho w offsets with respect to the bright points in the SAR images. (f) and (g) sho w the results of matching after using the ICP algorithm on the descending and the ascending image, respectiv ely . In the last two subﬁgures, it can be seen that the detected objects from the optical image (green circles) coincide with the visible bright points in the SAR images. 6 and Fig. 7 summarize the quality of these PSs. Fig. 6 shows the posterior estimated standard deviations of observations. The S az values vary from 3 cm to 19 cm with an av erage of 12 cm while the S rg values range from 1 cm to 24 cm with an av erage of 6 cm. This indicates that, for these natural PSs, the removal of error terms, as expressed in (1) and (2), and discarding the outliers allow sub-decimeter and decimeter precision in range and azimuth, respecti vely . In Fig. 7, the positioning quality is assessed by reporting the precision of the estimated coordinates. The standard deviations are deﬁned in the local coordinate system of Berlin in north, east, and height ( S N , S E , S H ) in the conﬁdence level of 95 % . The mean values of S N , S E and S H are 2.7 cm, 2.8 cm and 2.2 cm, respectiv ely . The higher precision in the height direction is merely the effect of the cross-heading geometry used for the positioning. Until no w the discussion was mostly focused on analyzing the relative accuracy of the coordinates based on the posterior precision estimates. Although, the most reasonable procedure to validate the absolute accuracy is a point-wise comparison of stereo SAR coordinates with respect to the corresponding GNSS-surve yed ones, this was not applicable at the time of the study . Instead, the stereo SAR estimated ellipsoidal heights of the 9 GCPs were compared with the LiD AR heights of IEEE TRANSA CTIONS ON GEOSCIENCE AND REMO TE SENSING, IN PRESS 11 (a) PS 41 (b) PS 43 Fig. 5: Range and azimuth residuals for two PS examples of the Berlin test case after the stereo SAR processing. Step 1 (gross outlier detection) and step 2 ( 2 σ test) hav e already been applied. the same area. W e assume the phase centers of the detected GCPs are at the base of the lamp poles on the ground. Therefore, the cross-comparison includes ﬁnding the nearest neighbors of the GCP candidate in the LiD AR point cloud within the radius of 1 m, excluding the LiD AR points with large height values which originate from the top of the lamp pole, estimating the mode of the LiD AR heights to represent the reference height and ev aluating the difference between the ellipsoidal height of stereo SAR results with respect to the reference height. The radius of the neighborhood is chosen in a way that the reference height calculation includes a reasonable number of samples and still be small enough to possibly prev ent the inclusion of different objects in the search window . It is also worth to note that the calculation of mode is carried out with the assumption that the majority of samples in the window stem from the ground. The results of the cross- comparison are reported in T ab . II. The estimated stereo SAR and approximated LiD AR reference heights are denoted by h S and h L while their offset is represented by h o . It is seen that for all except for one of the GCPs the height offset is below 20 cm. The results report a bias of 13 cm and a precision of 5 cm overall with respect to the LiD AR data which roughly implies the absolute accuracy of the height estimation using the stereo SAR method. B. Oulu The second test site covers the entire city of Oulu. The SAR data include four stacks of TS-X high resolution spotlight products with 177 images in total. The images were acquired from May 2014 to October 2016, from two ascending orbits IEEE TRANSA CTIONS ON GEOSCIENCE AND REMO TE SENSING, IN PRESS 12 PS 1 PS 7 PS 11 PS 13 PS 14 PS 18 PS 19 PS 20 PS 43 Candidate GCPs 0 0.05 0.1 0.15 0.2 0.25 Standard deviation [m] S A42 S A57 S R42 S R57 Fig. 6: The posterior azimuth and range standard deviations ( S az , S rg ) of the 9 best GCP candidates estimated by VCE in the geodetic stereo SAR processing. PS 1 PS 7 PS 11 PS 13 PS 14 PS 18 PS 19 PS 20 PS 43 Candidate GCPs 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 Standard deviation [m] S N S E S H Fig. 7: The posterior standard de viations scaled to 95 % conﬁdence lev el of the 9 best GCP candidates estimated posterior to the geodetic stereo SAR processing. The standard de viations are deﬁned in the local coordinate system of Berlin in north, east and vertical direction ( S N , S E , S H ). and two descending orbits. The acquisition parameters of the Oulu data set are reported in T ab . III while the mean scene cov erage and the acquisition time plot of the TS-X images are shown in Fig. 8. No images were ordered during the periods from November 2014 to March 2015 and November 2015 to March 2016 due to the accumulation of sno w expected for Oulu during the winter months. For Oulu, no optical images with sufﬁcient spatial resolution were a v ailable for the detection of PS candidates from cross- heading geometries. Therefore, the road network data of Oulu was used instead. The data was freely accessed from the Finnish Transport Agency [39]. It was delivered in vector format in the UTM coordinate system and includes the main streets and highways of Oulu. The detection of PS candidates from the same-heading tracks was carried out using the multitrack PSI fusion algo- rithm described in Subsection III-A. As the prerequisite of the algorithm, the PSI processing was performed by the PSI- GENESIS of the German Aerospace Center (DLR) [33]. For each detected PS, the ele vation and the deformation parameters (in this study only a linear trend) were estimated. As an example, the radar-coded PS elev ation map of the ascending stack of beam 30 is visualized in Fig. 9. After geocoding, the 3-D point clouds obtained from either ascending-ascending (AA) or descending-descending (DD) geometries form the initial input of the fusion algorithm as Fig. 10 demonstrates. Fig. 10a shows the geocoded PS point clouds from the DD geometries, visualized in white and gray . The yellow points represent the identiﬁed PS pairs from the fusion algorithm. The total number of the correspondences is approximately IEEE TRANSA CTIONS ON GEOSCIENCE AND REMO TE SENSING, IN PRESS 13 (a) May 2014 Aug 2014 Nov 2014 Feb 2015 May 2015 Aug 2015 Nov 2015 Feb 2016 May 2016 Aug 2016 Nov 2016 Beam30 Beam54 Beam69 Beam94 (b) Fig. 8: (a) The mean scene coverage of the TS-X images overlaid on the optical image of Oulu taken from Google Earth. (b) The acquisition time plot of the TS-X images of Oulu. Fig. 9: PS elev ation map obtained from PSI processing of an ascending stack of Oulu (beam 30). The total number of scatterers is approximately 540000 after selecting only the PSs with posterior coherence values equal or higher than 0.7. 32000 and the Euclidean distances between the matched PSs vary from 1.5 to 5 meters. In order to reduce the number of PS correspondences to the ones with higher quality and closer distance, as well as to preserv e the homogeneity of the distribution, a regular grid was imposed on the point clouds. Inside the 10 m × 10 m grid, the PS pairs which were closer together and had lo wer ADI values were selected to reduce the number of pairs from 32000 to 10000. The comparison between PS pairs before and after reduction can be seen in a zoomed-in area in Fig. 10b and Fig. 10c, respectively . The same procedure depicted in Fig. 10 was also carried out for the SAR images from the AA geometries and close to 9500 PS correspondences were detected. The results then were radar- coded for both geometry conﬁgurations. The PS candidates to be localized from the cross-heading geometries, either ascending-descending (AD) conﬁguration or quad geometry (AD AD) conﬁguration, were selected based on the detection of bright points along the roads using the road network data as was explained in Subsection III-C. The road network data was ﬁrst radar-coded on the master scenes of all the four geometries as is seen on the mean intensity images of each stack in Fig. 11. A circular neighborhood with a radius of 70 pixels was then considered around each road network node within which the ADI was ev aluated for all pixels for all the four stacks. The neighborhood is chosen based on a rough knowledge on the maximum width of highways in Oulu ( ≈ 35 m) and was adapted to the SAR data by taking into account the pixel spacing in range and azimuth direction IEEE TRANSA CTIONS ON GEOSCIENCE AND REMO TE SENSING, IN PRESS 14 (a) (b) (c) Fig. 10: Depiction of PS correspondence detection from SAR images of same-heading orbit tracks of Oulu using the ﬁrst step of the multitrack PS fusion algorithm proposed in [22]. (a) shows the geocoded PS point clouds of beams 54 and 94 in a DD conﬁguration as white and gray points as well as the detected PS correspondences in yellow . (b) and (c) show a zoom-in area of (a) marked by the red rectangle before and after imposing a 10 m grid in which the pairs with closest distance are chosen. and the ov ersampling factor used in the processing. After selecting the pixel with the lowest ADI in each neighborhood, a further threshold of 0.25 on ADI values, typically used in PSI processing [31], selects the stable bright points from each stack. After geocoding, the PSs from dif ferent stacks which are closer than a threshold of three meters are chosen as the ﬁnal stereo candidates and are subsequently radar-coded on all the SAR images. The distance threshold depends on the geometry conﬁguration from which the user is interested to localize the targets. If same-heading geometries are considered the value should be lower than three meters in order to ensure correct PS correspondence detection. An example of PS candidates visible from ADAD conﬁguration in Oulu is gi ven in Fig. 12. The candidates are all assumed to be bases of lamp poles and can be seen as bright points inside the green circles. The explained procedure produced 107 and 52 initial candidates from ADAD and AD geometry conﬁgurations, respectiv ely . The quantity is lo wer in the latter because of the strict distance threshold of 1.5 m imposed on coordinate differences of PSs visible in different stacks. The threshold value is chosen IEEE TRANSA CTIONS ON GEOSCIENCE AND REMO TE SENSING, IN PRESS 15 T ABLE II: T H E R E S U L T O F C R OS S - C O M P A R I S ON B E T W E E N T H E E S T I MAT ED H E I GH T S O F S T E R E O S A R h S A N D T H EI R C O R R E - S P O N DI N G L I DA R H E I G H T S h L . T H E O FF S E T h o I S A N I N D I C A T O R F O R T H E AB S O L U T E A C CU R AC Y O F h S GCP h L [m] h S [m] h o = h S − h L [m] PS 1 74.64 74.80 0.16 PS 7 74.45 74.54 0.09 PS 11 74.83 74.99 0.16 PS 13 75.40 75.62 0.22 PS 14 73.87 73.96 0.09 PS 18 73.87 73.95 0.08 PS 19 75.59 75.76 0.17 PS 20 75.02 75.18 0.16 PS 43 79.78 79.85 0.07 Mean 0 . 13 Standard de viation 0 . 05 T ABLE III: A C Q U I S I TI O N PA R A ME T E R S O F S TAC K S O F S A R I M - AG E S I N O U L U Beam Nr . Incidence angle (degree) Heading angle (degree) T rack type Nr . of images 30 30.9 346.1 Asc 44 54 41.1 191.4 Dsc 44 69 46.2 350.0 Asc 38 94 53.4 187.5 Dsc 51 empirically based on the histogram of minimum Euclidean distances e v aluated between the PS pairs. The PT A was performed on all the detected PS candidates in all the SAR images in which the candidate was visible. For each candidate a time-series of σ φ was ev aluated using (10). Fig. 13 demonstrates the initial outlier remo val on two PS candidates based on the method of adjusted boxplot explained in Subsection IV -C. The σ φ values of each scatterer are sorted in time in Fig.13a and Fig.13b . The distribution of σ φ is right- ske wed in both cases as can be seen in Fig.13c and Fig.13d. From the distribution plots, one can detect the samples which are not connected to the tail of the distribution and mark them as potential outliers. This process is done automatically and without the need for manual intervention using the interval deﬁned in (12). In Fig.13e and Fig.13f, the detected outliers are marked with red rectangles. For all the candidates the mentioned procedure was carried out. The atmospheric corrections were carried out using global ionospheric maps and the global tropospheric zenith path delays provided with the V ienna mapping function [40] since there was no access to the Oulu GNSS receiver at the time of the study . The geodynamic ef fects were fully considered according to the IERS con ventions and all the ef fects were remov ed from the PS timings. The ﬁnal positioning was carried out by stereo SAR with all the mentioned outlier removal steps in Subsection IV -E. The last criterion, which includes the removal of PSs based on S az values higher than 20 cm, reduces the amount of total PSs to only those points which can be considered ideal stereo candidates. The av eraged quality of the estimated 3-D coor - dinates for all the remaining high quality PSs are reported in T ab. IV. The scatterers are categorized based on the geometry conﬁguration that is used for their positioning. The results are all expressed in centimeter and are deﬁned in the local east, north and vertical coordinates. The standard deviations are all deﬁned within 95 % conﬁdence interval. From T ab . IV, it is seen that the averaged precisions are smaller than two decimeter for all the cases. As it was expected, the localization quality boosts as the difference in the viewing geometries becomes larger which is the case when changing from AA or DD to the AD and AD AD geometry conﬁgurations. It is also evident that for cross-heading geometries the retriev al of the height component is the most precise one as for the same-heading cases, the precision in the north component is the highest. Therefore, in general localization of targets from cross-heading tracks are desirable. The only remaining concern regarding localization using cross-heading tracks is the diameter of the lamp poles which may worsen the accuracy in the east coordinate component. This bias can be estimated and remov ed if the scatterer is also visible from same-heading tracks which is usually the case. The distribution of the total 2049 generated GCPs is visu- alized on the optical image of Oulu in Google Earth in Fig. 14. The scatterers are color-coded based on the underlying ge- ometry conﬁguration used for their localization. It is seen that almost the entire area of Oulu is covered with the generated GCPs. The ones from the same-heading geometries co ver the built areas while the ones from cross-heading conﬁguration include the base of lamp poles, street lights and traf ﬁc lights. V I . S U M M A RY , C O N C L U S I O N A N D O U T L O O K In this paper we described a processing chain for automatic detection and positioning of opportunistic PSs which are visible from multi-aspect TS-X images. This paves the way for generation of GCPs from SAR data only . Three algorithms ha ve been recommended for identical PS detection which are different in terms of number of generated GCPs and the subsequent positioning precision. The method based on the PSI fusion algorithm is able to provide point correspondences e ven on buildings and areas with complex scattering mechanisms in SAR images. Therefore, a large number of potential PS pairs can be obtained which normally cov er the entire scene. The downside of the method is that SAR image stacks are required for which a complete PSI processing has to be performed separately before being able to ﬁnd the PS correspondences. Furthermore, the method is usually applicable only for same-heading PS point clouds. Consequently in terms of localization precision with the subse- quent stereo SAR, the relati ve error in the cross-range is larger than the error for range and azimuth components. Another disadvantage is that many of the initial PS correspondences cannot be considered useful candidates for stereo SAR and hav e to be eliminated later in the processing, because the registration is only performed within the limits of the PSI 3-D localization quality for which PS pairs with distances of up to ﬁv e meters are detected. The detection algorithm based on external optical data is quite straightforward to implement and provides identical scatterers that are visible from cross- heading orbits leading to higher localization precision. The IEEE TRANSA CTIONS ON GEOSCIENCE AND REMO TE SENSING, IN PRESS 16 (a) (b) (c) (d) Fig. 11: Projection of the road network data of Oulu onto the master scene of (a) beam 30, (b) beam 54, (c) beam 69 and (d) beam 94. The road data is represented by red points and is the basis for detection of identical PS candidates from cross-heading orbit geometries. T ABLE IV: A V E R AG E D S TA T I S T I C S B A S ED O N T H E L E A S T S Q UA R E S E S T I M A T E D 3 - D C O O R D I NA T E S T A N DA R D D E V I A T I O N S I N O U LU . T H E L E T TE R S A A N D D S TA N D F O R A SC E N D I N G A N D D E SC E N D I N G G E O ME T R I E S , R E S PE C T I V E L Y . T H E S A M PL E M E A N A N D S TAN D A RD D E V I A T I O N A R E D E N OT E D B Y µ A N D σ , A N D S [ E N H ] R E P R ES E N T T H E L O C A L C O O R D I NA T E S S T A N DA R D D E V I A T I O N S W I T H I N 9 5 % C O N FI DE N C E L EV E L Geometry Nr . Scatterers µ s E [ cm ] µ s N [ cm ] µ s H [ cm ] σ s E [ cm ] σ s N [ cm ] σ s H [ cm ] AA 565 17.73 5.04 15.87 11.98 2.63 11.09 DD 1417 15.08 3.80 16.71 10.38 2.10 11.30 AD 24 2.26 2.50 1.75 0.99 1.11 0.83 AD AD 43 1.17 1.40 1.12 0.42 0.55 0.37 disadvantage of the method is that for reliable detection of lamp poles, the spatial resolution of the optical image should be in the sub-decimeter regime. Moreover , the method is highly prone to detecting other linear structures as shado ws of lamp poles and therefore more sophisticated object detection algorithms are recommendable. The method based on vector road network data, similar to the optical method, provides candidates which are suitable to be localized from cross- heading geometries. Also the external data is freely accessible for most locations. The disadvantage of the method is that a co-registration on one master has to be carried out for each stack and the amplitude data must be calibrated. It has been sho wn that the GCP generation processing chain is quite ﬂexible as it allows the user to constraint the number and the quality of the candidate GCPs either from the start of the procedure, by selecting different distance thresholds for detection or trimming the data based on estimated phase noise time series, or at the ﬁnal step of the procedure based on the posterior quality indicators obtained from stereo SAR. By applying the algorithm to two test sites, it has been demonstrated that it is capable of positioning natural PSs with precision values ranging from 2 cm and 4-5 cm, for cross-heading AD and AD AD conﬁgurations respectively , to approximately 20 cm for candidates extracted from same- IEEE TRANSA CTIONS ON GEOSCIENCE AND REMO TE SENSING, IN PRESS 17 (a) (b) (c) (d) Fig. 12: PS correspondence detection from ADAD geometry conﬁguration. For each image, the respecti ve a veraged heading angle α and the av eraged incidence angle θ inc are stated. heading geometries. As it was expected, the difference in the viewing geometries of the observed PS has the highest impact on the ﬁnal localization precision followed by number of acquisitions used in stereo SAR processing, the SCR of the target and the quality of external error corrections. Furthermore, as a preliminary cross-comparison, the estimated ellipsoidal height of the retrieved candidates in Berlin were compared to the corresponding height of a LiD AR data which reported an average bias of approximately 13 cm. The produced absolute GCPs have ample of applications in geodesy and absolute mapping. These points may substitute the conv entional GCPs that are required for geo-referencing of satellite imagery which are usually surveyed at the ﬁeld by GNSS. They can be further integrated as absolute reference points into multi-pass InSAR techniques. Furthermore, they can be used to detect long-term ground motions with small magnitudes and low spatial frequency which are in visible to phase-based InSAR methods. The future work will focus on smart pre-selection of the PS candidates by including the information obtained from PT A, using integrated side-lobe ratio (ISLR), to robustly remove the PS candidates which are located too close to each other . Furthermore, the stereo SAR processing could beneﬁt from weighting the initial timing observations of the PSs based on their respectiv e SCR or ADI values and also can be carried out with robust parameter estimation schemes like M-estimator . Moreov er , it is desirable to carry out GNSS measurements at selected test sites to be able to correctly validate the absolute accuracy of the generated GCPs. Finally , it is important to note that the proposed methodology is tailored to detection and absolute localization of GCPs in urban area where a large number of PSs are av ailable. In applications where the in vestigated scene includes mainly non-urban area, it is recommendable to employ artiﬁcial PSs such as corner reﬂectors or active transponders. A C K N O W L E D G M E N T From the Remote Sensing T echnology Institute (IMF) of the DLR, the authors would like to thank Dr . Ulrich Balss for pro- viding the routines for PT A, Mr . Fernando Rodriguez Gonzalez for his technical support on PSI processing of Oulu using the PSI-GENESIS and Mr . Nan Ge for re-ordering the T erraSAR- X images of Berlin with updated L1B product ﬁles. W e are also grateful to Dr . Heiko Hirschm ¨ uller of the DLR robotics institute for pro viding us with the optical data of Berlin. The LiD AR data of Berlin ha ve been provided by Land Berlin (EU EFRE project) and Landesamt f ¨ ur V ermessung und Geoin- formation Bayern. The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V . (www .gauss-centre.eu) for funding this project by providing computing time on the IEEE TRANSA CTIONS ON GEOSCIENCE AND REMO TE SENSING, IN PRESS 18 May 2014 Aug 2014 Nov 2014 Feb 2015 May 2015 Aug 2015 Nov 2015 Feb 2016 May 2016 Aug 2016 Nov 2016 0 0.02 0.04 0.06 0.08 0.1 [rad] (a) May 2014 Aug 2014 Nov 2014 Feb 2015 May 2015 Aug 2015 Nov 2015 Feb 2016 May 2016 Aug 2016 Nov 2016 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 [rad] (b) (c) (d) May 2014 Aug 2014 Nov 2014 Feb 2015 May 2015 Aug 2015 Nov 2015 Feb 2016 May 2016 Aug 2016 Nov 2016 0 0.02 0.04 0.06 0.08 0.1 [rad] (e) May 2014 Aug 2014 Nov 2014 Feb 2015 May 2015 Aug 2015 Nov 2015 Feb 2016 May 2016 Aug 2016 Nov 2016 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 [rad] (f) Fig. 13: Example of initial outlier removal based on phase noise time series of the detected candidates. The outliers are identiﬁed and remov ed automatically based on the distribution of σ φ . GCS Supercomputer SuperMUC at Leibniz Supercomputing Centre (LRZ, www .lrz.de). R E F E R E N C E S [1] M. Eineder, C. Minet, P . Steigenberger , X. Y . Cong, and T . Fritz, “Imaging Geodesy - T ow ard Centimeter-Le vel Ranging Accuracy With T erraSAR-X, ” IEEE T ransactions on Geoscience and Remote Sensing , vol. 49, no. 2, pp. 661–671, Feb. 2011. [2] X. Y . Cong, U. Balss, M. Eineder, and T . Fritz, “Imaging Geodesy - Centimeter-Lev el Ranging Accuracy W ith T erraSAR-X: An Update, ” IEEE Geoscience and Remote Sensing Letters , vol. 9, no. 5, pp. 948– 952, Sep. 2012. [3] C. Gisinger, U. Balss, R. Pail, X. X. Zhu, S. Montazeri, S. Gernhardt, and M. Eineder, “Precise Three-Dimensional Stereo Localization of Corner Reﬂectors and Persistent Scatterers With T erraSAR-X, ” IEEE T ransactions on Geoscience and Remote Sensing , vol. 53, no. 4, pp. 1782–1802, Apr . 2015. [4] S. Hackel, O. Montenbruck, P . Steigenberger , U. Balss, C. Gisinger, and M. Eineder, “Model improvements and validation of T erraSAR-X precise orbit determination, ” J ournal of Geodesy , Dec. 2016. [5] U. Balss, X. Y . Cong, R. Brcic, M. Rex er , C. Minet, H. Breit, M. Eineder , and T . Fritz, “High precision measurement on the absolute localization accuracy of T erraSAR-X, ” in Geoscience and Remote Sensing Sympo- sium (IGARSS) . IEEE, Jul. 2012, pp. 1625–1628. [6] U. Balss, C. Gisinger , X. Y . Cong, R. Brcic, S. Hackel, and M. Eineder , “Precise measurements on the absolute localization accurac y of T erraSAR-X on the base of far-distrib uted test sites, ” in 10th Eur opean Confer ence on Synthetic Aperture Radar EUSAR , 2014. [7] U. Balss, H. Breit, T . Fritz, U. Steinbrecher , C. Gisinger , and M. Eineder, “ Analysis of internal timings and clock rates of T erraSAR-X, ” in Geoscience and Remote Sensing Symposium (IGARSS) . IEEE, Jul. 2014, pp. 2671–2674. [8] I. G. Cumming and F . H.-c. W ong, Digital pr ocessing of synthetic apertur e radar data: algorithms and implementation , ser . Artech House remote sensing library . Boston: Artech House, 2005. [9] H. Breit, E. B ¨ orner , J. Mittermayer, J. Holzner , and M. Eineder, “The T erraSAR-X Multi-Mode SAR Processor - Algorithms and Design, ” in 5th Eur opean Confer ence on Synthetic Aperture Radar EUSAR , 2004. [10] C. Gisinger, “Atmospheric corrections for T erraSAR-X deriv ed from GNSS observations, ” Master Thesis, T echnische Uni versit ¨ at M ¨ unchen, M ¨ unchen, 2012. [11] G. Petit and B. Luzum, “IERS Conv entions, ” V erlag des Bundesamts f ¨ ur Kartographie und Geod ¨ asie, Frankfurt am Main, IERS T echnical Note No. 36, 2010. [12] M. Eineder , C. Gisinger, U. Balss, X. Y . Cong, S. Montazeri, S. Hackel, F . Rodriguez Gonzalez, and H. Runge, “SAR Imaging Geodesy– Recent Results for TerraSAR-X and for Sentinel-1, ” in FRINGE 2017, 10th International W orkshop on Advances in the Science and Applications of SAR Interfer ometry and Sentinel-1 InSAR , Jun. 2017. [13] C. Gisinger , S. Gernhardt, S. Auer , U. Balss, S. Hackel, R. Pail, and IEEE TRANSA CTIONS ON GEOSCIENCE AND REMO TE SENSING, IN PRESS 19 Fig. 14: T otal number of 2049 GCPs in Oulu color-coded based on the geometry conﬁguration used for their positioning. The underlying optical image is taken from Google Earth. M. Eineder , “ Absolute 4-D positioning of persistent scatterers with T erraSAR-X by applying geodetic stereo SAR, ” in Geoscience and Remote Sensing Symposium (IGARSS) . IEEE, Jul. 2015, pp. 2991– 2994. [14] X. X. Zhu, S. Montazeri, C. Gisinger , R. F . Hanssen, and R. Bamler , “Geodetic SAR T omography, ” IEEE T ransactions on Geoscience and Remote Sensing , vol. 54, no. 1, pp. 18–35, Jan. 2016. [15] S. Montazeri, X. X. Zhu, M. Eineder, and R. Bamler, “Three- Dimensional Deformation Monitoring of Urban Infrastructure by T o- mographic SAR Using Multitrack T erraSAR-X Data Stacks, ” IEEE T ransactions on Geoscience and Remote Sensing , vol. 54, no. 12, pp. 6868–6878, Dec. 2016. [16] U. Balss, H. Runge, S. Suchandt, and X. Y . Cong, “ Automated extraction of 3-D Ground Control Points from SAR images - an upcoming novel data product, ” in Geoscience and Remote Sensing Symposium (IGARSS) . IEEE, Jul. 2016, pp. 5023–5026. [17] D. O. Nitti, A. Morea, R. Nutricato, M. T . Chiaradia, C. La Mantia, L. Agrimano, and S. Samarelli, “ Automatic GCP extraction with high resolution COSMO-SkyMed products, ” in SAR Image Analysis, Mod- eling, and T echniques XVI, SPIE 10003 , C. Notarnicola, S. Paloscia, N. Pierdicca, and E. Mitchard, Eds., Oct. 2016, p. 1000302. [18] P . Dheenathayalan, D. Small, A. Schubert, and R. F . Hanssen, “High- precision positioning of radar scatterers, ” Journal of Geodesy , vol. 90, no. 5, pp. 403–422, May 2016. [19] S. Suri, P . Schwind, J. Uhl, and P . Reinartz, “Modiﬁcations in the SIFT operator for effecti ve SAR image matching, ” International Journal of Image and Data Fusion , v ol. 1, no. 3, pp. 243–256, Sep. 2010. [20] D. G. Lowe, “Distinctiv e Image Features from Scale-Inv ariant K ey- points, ” International Journal of Computer V ision , vol. 60, no. 2, pp. 91–110, Nov . 2004. [21] F . Dellinger, J. Delon, Y . Gousseau, J. Michel, and F . Tupin, “SAR- SIFT: A SIFT-Like Algorithm for SAR Images, ” IEEE T ransactions on Geoscience and Remote Sensing , vol. 53, no. 1, pp. 453–466, Jan. 2015. [22] S. Gernhardt, X. Y . Cong, M. Eineder, S. Hinz, and R. Bamler, “Geometrical Fusion of Multitrack PS Point Clouds, ” IEEE Geoscience and Remote Sensing Letters , vol. 9, no. 1, pp. 38–42, Jan. 2012. [23] S. Gernhardt, “High Precision 3D Localization and Motion Analysis of Persistent Scatterers using Meter-Resolution Radar Satellite Data, ” Doctoral Thesis, T echnische Uni versit ¨ at M ¨ unchen, M ¨ unchen, 2012. [24] B. M. Kampes, Radar Interfer ometry: P ersistent Scatter er T echnique . Dordrecht, the Netherlands: Springer, 2006. [25] T . Fritz and M. Eineder , “T erraSAR-X Ground Segment Basic Product Speciﬁcation Document, ” German Aerospace Center (DLR), Oberpfaf- fenhofen, T echnical report, 2008. [26] S. Montazeri, X. X. Zhu, U. Balss, C. Gisinger, Y . W ang, M. Eineder , and R. Bamler, “SAR ground control point identiﬁcation with the aid of high resolution optical data, ” in Geoscience and Remote Sensing Symposium (IGARSS) . IEEE, Jul. 2016, pp. 3205–3208. [27] R. C. Gonzalez and R. E. W oods, Digital imag e pr ocessing , 2nd ed. Upper Saddle River , N.J: Prentice Hall, 2002. [28] K. Briechle and U. D. Hanebeck, “T emplate matching using fast normalized cross correlation, ” in Optical P attern Recognition XII, SPIE 4387 , D. P . Casasent and T .-H. Chao, Eds., Mar. 2001, pp. 95–102. [29] D. Comaniciu and P . Meer , “Mean shift: a robust approach toward feature space analysis, ” IEEE T ransactions on P attern Analysis and Machine Intelligence , vol. 24, no. 5, pp. 603–619, May 2002. [30] P . Besl and N. D. McKay , “ A method for registration of 3-D shapes, ” IEEE T ransactions on P attern Analysis and Machine Intelligence , vol. 14, no. 2, pp. 239–256, Feb. 1992. [31] A. Ferretti, C. Prati, and F . Rocca, “Permanent scatterers in SAR interferometry , ” IEEE T ransactions on Geoscience and Remote Sensing , vol. 39, no. 1, pp. 8–20, Jan. 2001. [32] N. Adam, B. M. Kampes, M. Eineder, J. W orawattanamateekul, and M. Kircher, “The Development of a Scientiﬁc Permanent Scatterer System, ” in Proceedings of the Joint ISPRS/EARSeL W orkshop - High Resolution Mapping from Space 2003 , Hannover , 2003. [33] N. Adam, F . R. Gonzalez, A. Parizzi, and R. Brcic, “W ide area Per- sistent Scatterer Interferometry: Current developments, algorithms and examples, ” in Geoscience and Remote Sensing Symposium (IGARSS) . IEEE, Jul. 2013, pp. 1857–1860. [34] N. Adam, B. M. Kampes, and M. Eineder, “The development of a scien- tiﬁc persistent scatterer system: Modiﬁcations for mixed ERS/ENVISA T time series, ” in Proceedings of the Envisat & ERS Symposium , Salzburg, Austria, 2004, pp. 1–9. [35] M. Hubert and E. V andervieren, “An adjusted boxplot for skewed distributions, ” Computational Statistics & Data Analysis , vol. 52, no. 12, pp. 5186–5201, Aug. 2008. [36] J. W . Tuk ey , Exploratory data analysis , ser . Addison-W esley series in behavioral science. Reading, Mass: Addison-W esley Pub. Co, 1977. [37] G. Brys, M. Hubert, and A. Struyf, “A Robust Measure of Skewness, ” Journal of Computational and Graphical Statistics , vol. 13, no. 4, pp. 996–1017, Dec. 2004. [38] H. Hirschmuller, “Stereo Processing by Semiglobal Matching and Mu- IEEE TRANSA CTIONS ON GEOSCIENCE AND REMO TE SENSING, IN PRESS 20 tual Information, ” IEEE T ransactions on P attern Analysis and Mac hine Intelligence , v ol. 30, no. 2, pp. 328–341, Feb . 2008. [39] “Finnish T ransport Agency, ” Jan. 2017. [Online]. A vailable: http: //www .liikennevirasto.ﬁ/web/en [40] J. Boehm, B. W erl, and H. Schuh, “Troposphere mapping functions for GPS and very long baseline interferometry from European Centre for Medium-Range W eather Forecasts operational analysis data, ” Journal of Geophysical Resear ch: Solid Earth , vol. 111, no. B2, pp. 1–9, Feb . 2006.

Automatic Detection and Positioning of Ground Control Points Using TerraSAR-X Multi-Aspect Acquisitions

Original Paper

Comments & Academic Discussion

Leave a Comment

Original Paper

Related Papers

Comments & Academic Discussion

Leave a Comment