Optimized GNSS Cal/Val Site Selection for Expanding InSAR Viability in Areas With Low Phase Coherence: A Case Study for Southern Louisiana

Interferometric synthetic aperture radar (InSAR) techniques can be used to derive spatially dense “relative” measurements of vertical land motion (VLM), whereas global navigation satellite system (GNSS) provides point-based “absolute” measurements of VLM. The combination of GNSS and InSAR observations can yield spatially dense VLM measurements in an absolute reference frame. In addition, GNSS observations can be used to correct atmospheric noise in InSAR deformation measurements and serve as a complementary measure to isolate deep and shallow subsidence components. Given the increasing spatial and temporal coverage available from InSAR satellites, there is a need to establish calibration/validation networks that enable the use of InSAR for measuring VLM in coherence-challenged areas such as many low-lying coastal lands. In this study, we provide a method for the selection of sites for new GNSS installations such that the resulting GNSS network can better serve as tie points and validation for InSAR in areas where low coherence prevents high-fidelity phase unwrapping. Our method is applied in a case study for expanding the existing GNSS network in southern Louisiana, using abandoned oil well sites as potential sites. Considering practical limitations, distribution among various land classes, and following National Geodetic Survey guidelines, our proposed GNSS network consists of 61 (45 existing + 16 new) stations spread over a 50 000 km2 area of southern Louisiana.

Abstract-Interferometric synthetic aperture radar (InSAR) techniques can be used to derive spatially dense "relative" measurements of vertical land motion (VLM), whereas global navigation satellite system (GNSS) provides point-based "absolute" measurements of VLM.The combination of GNSS and InSAR observations can yield spatially dense VLM measurements in an absolute reference frame.In addition, GNSS observations can be used to correct atmospheric noise in InSAR deformation measurements and serve as a complementary measure to isolate deep and shallow subsidence components.Given the increasing spatial and temporal coverage available from InSAR satellites, there is a need to establish calibration/validation networks that enable the use of InSAR for measuring VLM in coherence-challenged areas such as many low-lying coastal lands.In this study, we provide a method for the selection of sites for new GNSS installations such that the resulting GNSS network can better serve as tie points and validation for InSAR in areas where low coherence prevents high-fidelity phase unwrapping.Our method is applied in a case study for expanding the existing GNSS network in southern Louisiana, using abandoned oil well sites as potential sites.Considering practical limitations, distribution among various land classes, and following National Geodetic Survey guidelines, our proposed GNSS network consists of 61 (45 existing + 16 new) stations spread over a 50 000 km 2 area of southern Louisiana.Index Terms-Global navigation satellite system (GNSS), interferometric synthetic aperture radar (InSAR), NASA-ISRO synthetic aperture radar (NISAR), phase coherence, vertical land motion (VLM).

I. INTRODUCTION
L AND subsidence is a major geohazard that disproportion- ally affects low-lying coastal and deltaic areas [1], most of which are already at risk from decreased sediment input [2], [3].Given that land subsidence in combination with sea level rise has led to increased flood frequency and impact [4], [5], [6], more extensive, up-to-date, and accurate knowledge of surface elevation and subsidence trends is required to inform coastal sustainability and disaster management actions.Interferometric synthetic aperture radar (InSAR) is a powerful remote sensing technique for measuring relative surface elevation change [7], [8], which can be combined with reliable independent global navigation satellite system (GNSS) measurements to obtain absolute measurements in a geodetic reference frame [9], [10], [11], [12].Leveraging this technology for vertical land motion (VLM) measurements in at-risk areas is an attractive option given the availability of free and open InSAR data with near-global coverage from current and planned Earth-observing missions of the European Space Agency (ESA) [13], [14] and U.S. National Aeronautics and Space Administration / Indian Space Research Organisation (NASA/ISRO) [15].
Applying InSAR to measure VLM in undeveloped parts of deltas and coastal lowlands is often challenging [16] and references therein; [17].InSAR methods applied to nonurban areas with variable soil moisture and moderate-to-heavy precipitation and vegetation cover largely rely on pixels with stable returns in the SAR images known as persistent scatterers (PS) [18], [19].In general, well-distributed PS points are required throughout the SAR image for broad-scale measurements of VLM.This is not always possible in coastal areas.For example, the PS algorithms implemented in [19], [20] cannot identify a dense set of stable pixels in coastal Louisiana because the C-band Sentinel-1 SAR signal does not maintain InSAR coherence in the 12-day interval between acquisitions in many undeveloped areas [21].
There are two main issues limiting the use of InSAR for subsidence measurements in these areas.1) Land subsidence applications need VLM measurement over a large area, on the same scale or larger than that of the atmospheric noise, which must be corrected to obtain accurate results.However, the subsidence signals of interest (SOIs) are typically not uniform on that scale, because subsidence is driven by numerous natural and anthropogenic processes that act at different spatial and temporal scales [16], [22], [23], [24], [25].Even if InSAR can accurately resolve localized subsidence, it does not provide the information needed, which is accurate VLM across the entire area.Tie-points, generally provided by GNSS, are needed to separate the long-wavelength, spatially heterogeneous VLM from the long-wavelength, spatially heterogeneous noise (e.g., due to atmospheric delays).2) Vegetation, changing soil moisture, and tidal or seasonal inundation can cause temporal decorrelation (loss of interferometric phase coherence) even in the typical 12-to-24-day repeat orbit of InSAR satellites, which can lead to spatially disconnected InSAR phase measurements separated by areas where no VLM information is available.To avoid this issue, some studies have explored the possibility for attaching corner reflectors (CRs) or compact active transponders (CAT) to a GNSS station, effectively providing a reliable ground target for InSAR.Mahapatra et al. [26] established a theoretical framework for adding ground instrumentation in an area of known deformation.Their method finds the optimal locations for adding CRs to an area with low interferometric coherence such that subsidence parameter estimates (depth, source shape, etc.) can be retrieved to a desired accuracy.Such arrangements reasonably ensure that the return from the area covered by the InSAR pixel is mainly from the CR or CAT.This has mainly been applied in studies of small, isolated areas or localized features like landslides [27] or peatlands [28].In the absence of artificial targets, other instruments like extensometers, GNSS, rod surface elevation tables (RSET), and feldspar-marker horizons are used alone or in combination to measure elevation change.
In this article, we address the first issue, namely the need to correct long-wavelength noise in the presence of a heterogeneous long-wavelength VLM signal.We describe a method for selecting new sites for a GNSS network to serve as calibration and tie points for InSAR in areas challenged by a combination of low-interferometric phase coherence, sparse PS points, nonuniform long-wavelength deformation signals, and significant and variable atmospheric noise.The method uses the best-available information about the long-wavelength subsidence to optimally select GNSS station locations based on spatial gradients in VLM.We choose southern Louisiana as the case study (see Fig. 1) and identify specific locations where additional stations could be added to improve InSAR subsidence measurements in nonurban areas in the next decades.Harnessing the properties of SAR to determine locations of coherent InSAR measurements and installing a ground monitoring instrument (ex: GNSS) nearby is a path forward to improve knowledge of VLM in southern Louisiana.

II. INSAR VERSUS GNSS
InSAR measurements are relative to a local reference point and not tied to a geodetic reference frame.More specifically, In-SAR measures the cumulative line-of-sight (LOS) displacement at each ground pixel relative to a stable reference point over the period of observation, but only for those areas with coherent phase returns from the surface.A stable pixel or a set of stable pixels is chosen as the spatial phase reference, and any residual noise at the reference location propagates to all other pixels [11].Existing studies [10], [11] have chosen the InSAR reference point close to a GNSS station and used GNSS measurements to estimate and remove residual errors under the assumption that the GNSS and InSAR are measuring the same signal.This method works best if the extent of the area of interest is small and located near the chosen station, because InSAR tropospheric noise increases with distance from the reference point [29].
In areas where phase coherence is well-maintained between image acquisitions, the primary error sources for InSAR are long-wavelength errors due to changes in the troposphere [30] (particularly from turbulent water vapor, not correlated with topography), ionosphere [31], and errors in knowledge of satellite orbit or aircraft locations [32].Corrections are done using either geometric correlation [33], [34] or independent weather models for troposphere [35], [36], sensor frequency split-spectrum methods for ionosphere [37], and other empirical techniques for correcting orbital errors [38].GNSS zenith water vapor estimates within the SAR image frame can be used to correct tropospheric artifacts in interferograms [39].However, uncertainties in the corrections contribute to uncertainty in the VLM measurements [40], [41].
GNSS can be used to calibrate and validate InSAR measurements, as well as to convert relative VLM estimates to absolute VLM estimates.GNSS provides accurate point-location, 3-D displacement time series with high temporal resolution.GNSS measures the position of the station antenna, which is usually attached to a stable structure (rock, building) or anchored to a stable rock layer below the ground surface.Stations on the ground are sensitive to the displacements below the anchoring depth and stations on a building are sensitive to the displacements below the depth of the foundation.GNSS measurements do not capture movement above the footing depth, which can be significant depending upon the depth of footing and soil type and stratigraphy [42], [43].While most joint InSAR/GNSS studies for measuring long-wavelength deformation have focused on areas where there is a dense GNSS network and shallow subsidence does not dominate, one recent study has employed GNSS interferometric reflectometry (GNSS-IR) techniques to measure shallow subsidence rates at six sites in southern Louisiana [44].For example, Shen and Liu [45] were able to resolve horizontal tectonic motion in southern California to an uncertainty of 0.6-1.5 mm/yr by combining InSAR measurements with a GNSS network with spacing varying from 30-70 km.
InSAR measures the subsidence aggregated over the entire soil column, thus InSAR alone cannot separate the deep and the shallow subsidence components.It is possible to use a combination of InSAR plus GNSS to separate shallow and deep subsidence signals [46], [47], however, this approach requires careful selection of the tie points.Challenges for tie point networks come from non-uniform long-wavelength deformation SOIs whose spatial scale is comparable to that of the atmospheric noise signals.For many InSAR applications (land subsidence, tectonic, permafrost), the SOI covers a relatively large area (>100 km 2 ).Swath widths of typical SAR frames range from 30 km (Cosmo SkyMED) to 250 km [Sentinel-1, NASA-ISRO synthetic aperture radar (NISAR)], and the extent of long-wavelength noise sources can vary from ∼20 to 400 km [10], [37].Neely et al. [48] fitted a surface to the difference of In-SAR and GNSS measurements from sites uniformly distributed over the SAR footprint and added the surface back to the InSAR deformation field to filter out long-wavelength noise and obtain dense absolute land motion estimates, under the assumption that VLM above the GNSS footing depth was negligible.They determined that for their case, a small number of GNSS stations (at least 6 stations to estimate a 2nd order polynomial) spaced at ∼75 km and distributed across the SAR footprint can well-estimate the long-wavelength noise contributions.Thus, GNSS can be used to calibrate InSAR measurements to VLM by removing long-wavelength noise sources and validate InSAR measurements in locations with short-wavelength VLM.

III. SOUTHERN LOUISIANA GEOLOGY AND SUBSIDENCE
Southern Louisiana is a deltaic complex built mainly during the Holocene as a result of sea level fluctuations and riverine sediment deposition.The bulk of the sediments is terrigenous, deposited by the Mississippi River during the sea-level lows, and the river's meandering greatly influenced the formation of different deltaic complexes [49].Most of the shallow subsidence is concentrated in the Holocene layers deposited over a Pleistocene surface, and the thickness of the Holocene layers is directly related to the shallow subsidence [50], [51].Many studies have considered the top of the Pleistocene layer [∼0 to > 200 m deep, Fig. 1(b)] [52] as a sturdy base, and computed subsidence below it as deep subsidence and within the Holocene layer as shallow subsidence.Jankowski et al. [43] found that a significant portion of shallow subsidence occurs in the top 5-10 m of the Holocene layer.GNSS stations with their bases below 20 m are often considered to be measuring the deep subsidence component [42], [43], although this is not the top of the Pleistocene in much of eastern coastal Louisiana [see Fig. 1(b)].As a practical solution, the Louisiana Coastal Protection and Restoration Authority derives deep subsidence by a fit to VLM measured with GNSS at continuously operating reference stations (CORSs) and benchmarks located throughout southern Louisiana [53].Table B1 of [53] provides information on the depth of most of the GNSS sites they use for their deep subsidence measurements, which ranges between 10 and 25 m.
Present-day coastal Louisiana in its natural state is comprised of wetlands with different types of marshes (freshwater, brackish, and saline) and swamps [54].These play a major role in sediment retention and protect against marine erosion and inland flooding [55].In addition to consolidation of the soil column, the organic matter present in the soil is sensitive to water table changes and undergoes oxidation when dry. Physical properties such as thickness, compressibility, and lateral extent of each soil type influence the rate of sediment compaction and subsidence in the wetlands [56].In coastal Louisiana, marshes that cannot build and retain sediment sufficient to offset subsidence gradually convert into open water [57].
Southern Louisiana is also one of the hotspots for the increased threat of submergence and land loss due to modern-day sea level rise [58].The threat is further exacerbated by land subsidence on the coast due to factors such as soil consolidation, faulting, groundwater, and hydrocarbon extraction [59].Subsidence from these factors occurs at different depths and is distributed at different spatial scales.Leveling surveys from 1965 to 1982 showed that localized subsidence occurred on the hanging wall areas of the large tectonic faults in the Terrebonne basin and near subsurface salt domes [60].It has been suggested that a complex feedback-type relationship exists between fault and salt movement [61], causing increased subsidence near the salt domes.Jankowski et al. [43] reported subsidence of ∼5 mm/yr due to the consolidation of Holocene peat layers at a depth of 5-18 m.Using borehole pressure data, Mallman and Zoback [62] reported localized subsidence of 8 mm/yr due to sediment loading, with more subsidence closer to the mouth of the Mississippi River over oil and gas fields in southern Louisiana.Karegar et al. [63] reported increased GNSS horizontal motion of order >1 mm/yr due to downslope motion on listric normal faults south of Baton Rouge.Based on GNSS measurements over deep-rooted structures, Dokka et al. [64] attributed the tectonic subsidence in the area to downward movement along high-angle normal faults connected at depth to a detachment fault.Using airborne InSAR, localized subsidence due to anthropogenic activities was identified in urban New Orleans [59] and precursory subsidence before the collapse of the Bayou Carne sinkhole [65].Jankowski et al. [43] showed that subsidence differs between the Mississippi Delta and Chenier Plain with over 50% of the latter highly vulnerable to submergence in the future.
The existing subsidence measurements in southern Louisiana mostly rely on leveling surveys, GNSS, tide gauges, and rod surface elevation table-marker horizon (RSET-MH) stations.Leveling surveys measure land elevation changes from an established benchmark [60], but these surveys are labor-intensive and give only a few measurements.Tide gauges measure the sea level rise with respect to a land-based benchmark and differences in sea level rise estimates among the tide gauges isolate subsidence.Similar to GNSS, tide gauges do not measure the subsidence above the foundation depth of the benchmark (ranges 0.9-35.1 m; [66]).RSET instruments measure surface elevation change from an arm attached to the RSET-rod to the surface, which is more sensitive to sediment accretion than subsidence [67], [68].The RSET stations and the tide gauges are distributed mainly along the coast [see Fig. 1(a)] and are used to measure land loss/gain in the inundated wetlands.The subsidence measurements obtained from them can be skewed in the north-south direction due to their location distribution, but mainly fill in gaps where InSAR cannot be used.The Centre for Geoinformatics (C4G) at the Louisiana State University currently operates a state-wide GNSS network with 65 stations [see Fig. 1(a)] distributed over the entire state, in an area of ∼135 000 km 2 , which is sparse to monitor the ongoing subsidence given the 3-D-spatial heterogeneity of soil types and the multitude driving processes.Spaceborne InSAR surveys have been done before with Radarsat C-Band [69] and TerraSAR-X X-Band [70].Both studies detected localized subsidence in the urban New Orleans area but failed to provide any measurements in the vegetated delta plains and wetlands.More recent work based on a modified phase triangulation algorithm [71] has yielded subsidence estimates with ∼ 2 mm/year accuracy.Though at a coarse 1-km spatial resolution, these measurements are currently the best available region-wide estimates of subsidence in southern Louisiana.An illustration depicting the sensitivity of various monitoring instruments to the subsidence measurement is shown in Fig. 1(c).

IV. DATA AND PREPROCESSING
The method for identifying sites for new GNSS instrumentation sites that support correction of long wavelength InSAR noise in the presence of a nonuniform long wavelength SOI has three essential inputs.
1) Identification of locations where InSAR measurements of ground deformation are viable.This is necessary for the GNSS measurements to be used as tie-points for the InSAR.2) An estimation of the long wavelength SOI.This is necessary to identify locations where the spatial gradient is large so that new stations can better measure the change in those areas.
3) The candidate locations for new GNSS stations.The method is sufficiently general to be adapted for many other purposes where optimal site selection is not necessarily a uniform grid of stations.In the case study demonstrating the method, Sentinel-1 C-band data is used to identify areas where InSAR is viable based on interferometric coherence (see Section IV-A).Other bands or data sets could be used for this purpose, but since L-band SAR should maintain better coherence than C-band, the sites selected should work well for the upcoming NISAR mission that will image all coastal lands with 12-day repeats for sustained and comprehensive InSAR studies of low-lying lands.The background long-wavelength SOI for the case study used is that of Wang et al. [71] (see Section IV-B), but other estimates can be used, e.g., the deep subsidence estimate of [53] or for a different study area.The potential GNSS sites for our study are the abandoned oil and gas platforms in southern Louisiana without application of any other criteria (see Section IV-C), but sites could be selected from a different initial set or with additional criteria applied, e.g., adding that the site should be appropriate for GNSS-IR to measure the shallow subsidence component at the new GNSS sites [44], [72], [73].

A. Sentinel-1 Data and InSAR Processing
For this study, we analyze Sentinel-1 InSAR data from 2018-2021 covering southern Louisiana to determine the locations of C-band PS points.We limit our study to the state boundary and the four Sentinel-1/2 frames covering southern Louisiana [see Fig. 1(a)].We use single-look complex products (SLCs) from Sentinel-1A/B satellites in ascending pass configuration corresponding to Path165/Frame90 (P165F90), P165F95, P63F92, and P136F93, which in combination cover southern Louisiana including the city of New Orleans and Bird-foot delta [see Fig. 1(a)].More than 400 SAR scenes spanning four frames are used.Maximum spatial and temporal baselines in the stack are 1000 m and 60 days, respectively.Interferograms are processed using InSAR Scientific Computing Environment [74] and are unwrapped with SNAPHU [75].The interferogram network is connected in a hierarchical fashion such that each dataset forms interferograms with preceding and succeeding datasets, except the first and the last datasets which only form one interferogram.Fig. S1 shows the network for P63F92.The split-band spectrum technique is used to correct for the ionosphere errors [37].Unwrapping errors are identified and corrected using phase closure and island bridging [76], [77].Good coherence (>0.6) is maintained in urban areas like New Orleans, where there are Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
stable targets for InSAR processing, along the developed highways extending south towards the Gulf of Mexico, in isolated developed areas, and some PS points in wetlands.For some areas in southern Louisiana, especially in the Atchafalaya and Wax Lake Delta and nearby coastal wetlands and in Cameron Parish wetlands, InSAR coherence is low (<0.4), and the interferograms cannot be properly unwrapped.

B. Priori Estimate of Subsidence
To optimally locate the GNSS stations, we use a priori information to best estimate the spatial variation of subsidence.This allows us to identify the optimal locations for a minimum number of stations that can capture the spatial gradient of the subsidence to a specified accuracy.In our study area, the subsidence pattern cannot be defined as an analytical function, so we use the most recent, spatially extensive map, provided in [71].This is shown in Fig. 4(a), with the InSAR line of sight measurements converted to VLM assuming the horizontal deformation component is negligible.Wang et al. [71] combined adaptive Goldstein phase filtering and an eigenvalue decomposition of a weighted complex InSAR phase matrix to recover the spatially coherent phase measurements over densely vegetated terrain.Their subsidence map reveals significant subsidence in the Bird-foot delta area, from the eastern part of the Chenier plain to the west of Lafayette, and in the vicinity of Port Fourchon.

C. Potential GNSS Sites
Setting up a GNSS station network at new locations is a challenge as they need to be placed in sites with minimal disturbance, which are easily accessible, and at the same time maintain InSAR phase coherence for the 12-24 day orbit repeat period.We start by considering site accessibility.The Department of Natural Resources of the State of Louisiana maintains a database of well-site activity in the state.Hydrocarbon drilling was a major industry in southern Louisiana [78], and the state has many abandoned oil and gas well sites.These sites are viable for installing GNSS stations (locations can be accessed from https:// www.sonris.com/), because they present fewer hurdles to access and obtaining the owner's permission.Many of these well sites have a strong foundation to set up a GNSS monitoring station, and some of them have deeply drilled well casings onto which a GNSS instrument can be mounted.We incorporate the well-site information into our framework as possible locations for establishing new GNSS stations.Some of the sites can be eliminated later due to site conditions, redundancy, etc.

V. METHODS
We use the PS pixels and an existing InSAR surface deformation map to find viable sites to set up GNSS stations.Specifically, stations must be colocated with coherent SAR pixels or at least near an area where InSAR can measure coherent phase change given the time interval between repeated radar acquisitions.To derive the locations of stable InSAR returns (with Sentinel-1 C-band now and with NASA/ISRO L-band NISAR in future), we first identify stable InSAR PS pixels from the Sentinel-1 interferograms using interferometric properties (see Section V-A).From the orphaned well sites, we then identify sites that are closely located (e.g., within 5 km radius for our case study based on [53], [59], [71]) to the PS pixels and develop a methodology to identify the "best interpolation sites" among them (see Section V-B).The "best interpolation sites" are chosen to capture the spatial variation of ground displacement over the entire study area using simple interpolation when equipped with CORS GNSS.The final number of sites to instrument can be determined by the user considering available resources and the desired interpolation accuracy.

A. Identify InSAR PS
To achieve our goal of placing GNSS stations in close vicinity of stable InSAR measurements, we first need to determine pixels with stable returns in the interferogram stack.Each interferogram in the stack is a complex multiplication of two individual SAR SLCs calculated as where u i and u j are the complex representation of two different SAR SLCs forming the interferogram, A ij is the average amplitude of the two acquisitions, and Δϕ ij is the phase difference between the acquisitions.The amplitude represents the strength of the target return, and the phase contains the information about the LOS deformation that occurred between two radar acquisition times [30].The radar phase measurements can vary due to changes that occur during the interval between the acquisitions, and the quality of phase return can be measured by the spatial coherence (γ) defined as where E indicates the expected value computed as the moving average within a window of pixels (in this study 9 × 9).Rapidly changing targets such as water or heavily vegetated areas have lower γ, while stable targets with bright returns such as buildings have high γ.To identify pixels that maintain coherence over time, we compute for each pixel the mean of γ over the entire stack of interferograms, defined as the average spatial coherence, μ γ .
Interferometric phase, Δϕ ij , which is sensitive to the VLM between the acquisitions, is prone to decorrelation and can be biased by noise from various sources as described in Section II.Our objective is to identify pixels for which the decorrelation noise component is small, and the phase measurements are consistent over time.As γ is estimated on a window of pixels, some stable targets that are surrounded by incoherent areas can have low coherence.To overcome this, [79] proposed using another measure of radar return stability called Amplitude Dispersion (D A ), calculated as where σ A and μ A are the standard deviation and mean of the amplitude of all interferograms in the stack.The number of interferograms in a stack can vary for each SAR frame in the Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
study area.Low amplitude dispersion can indicate pixels less susceptible to temporal decorrelation and with high signal-tonoise ratio [79].To determine pixels with stable SAR returns, we identify pixels of low D A overlapping with pixels of high μ γ .Our methodology is similar to others used to identify stable In-SAR pixels as persistent scatters [18], [79], [80] and distributed scatterers [81] except that we do not use interferometric phase criteria (e.g., phase stability; [79]) as many of our interferograms are not properly unwrapped or have many isolated clusters of pixels that are well-unwrapped (connected components).Severe decorrelation results in such isolated coherent areas, which makes unwrapping challenging even with phase bridging techniques.From the InSAR stack of interferograms, we take the 1-yr dataset with the highest coherence, in our case June 2019 -June 2020 (Fig. S1 in supplement).We compute μ γ and D A for each pixel in the interferogram stack.Most of the vegetated areas, mainly swamps and coastal wetlands, have low coherence and high amplitude dispersion.For our analysis, we use μ γ threshold >0.6 and D A threshold <0.4 to choose our PS points.Pixels with μ γ and D A meeting the thresholds are shown in Fig. 2.Although some studies use a lower D A threshold, we use 0.4 to obtain some PS in the largely vegetated wetland parts of the study area such as Bird-foot and Atchafalaya deltas.In the future, L-band sensors such as NISAR and ROSE-L (ESA) should be able to perform better than C-band Sentinel-1 in terms of maintaining temporal coherence in these areas, so the PS pixels identified here should work for L-band sensors.We generate a map of low amplitude dispersion pixels overlapping areas of high coherence as PS pixels (127443 pixels) covering our study area [see Fig. 2(c)].For further analysis, we use the centers of these pixels as the PS points.Most PS points lie along the urban areas, small villages, roads, and concrete level structures.

B. Determining the Best GNSS Sites
A GNSS network can be established in an area with no prior GNSS stations or where new stations are needed to enhance the existing network.Here we define the methodology assuming no existing stations, and we then show how the method can be adjusted to accommodate the existing stations in the Louisiana case study in Section VI-A.
First, we extract the locations of orphaned well sites (categorized as not operational or plugged) in southern Louisiana [712 sites-Fig.3(a)].From these sites, we choose the sites that are located near the PS points, setting a threshold of 5 km radius around the PS for sufficiently close well sites [465 wells marked in green in Fig. 3(a)].These comprise the initial pool of potential sites for setting up a GNSS network.The threshold is assumed reasonable based on prior studies of subsidence in southern Louisiana [53], [59], [71], however, other conditions can be imposed based on the geologic conditions for different study areas.
Next, we apply the method described below to select the best site locations from the initial pool such that interpolating the GNSS measurements at those new locations can capture the long-wavelength VLM over the study area to a desired accuracy.If no prior information about subsidence rates is available, then a uniform zero-subsidence map can be used.This provides a tie-point network optimized for removing long-wavelength noise in areas with nonuniform long-wavelength VLM.These new GNSS locations thus determined are the "best interpolation sites" for a GNSS network.
Taking the existing surface deformation measurements [71] with m pixels as the best a priori estimate of ground displacement (Y = [y 1 , y 2 , y 3 , ….., y m ] T ) of the study area [see Fig. 4(a)], we define an initial GNSS network from a pool of n (465) available sites that best estimates Y using interpolation.We obtain the ground displacement measurements at the n sites ( To estimate Y using V sites, we design an inverse distance weighted (IDW) mapping function, G.The best interpolation sites among V are obtained by solving a linear regression system of the form Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.where Each row-element in G corresponds to the distance-weighted contribution of measurement at each site v j to the interpolated estimate for y i .X is the vector of multiplicative weights for each site to determine its relative importance to the interpolation, which is to be solved for.The sites with higher weight are then presumed to be the best for setting up a GNSS station.
In this study we use IDW interpolation to obtain Y from V sites, but there are many other spatial interpolation methods that can be used [82].Note that we weight the measurement at each site twice for interpolation, once based on the distance (IDW ijknown from the location of the pixel corresponding to y i and the location corresponding to v j ) and again with a site-specific weight x j .When the number of observations (m) exceeds the number of sites (n), the system is overdetermined, and we solve for X in a least-squares routine.We use a regularized least square solver (Ridge regression, [83]) to solve for X, which minimizes the cost function Here, . 2 denotes the Euclidean (L2) norm.λ is a penalty factor that controls the variability of individual site weights in X and results in an even distribution of weights.This prevents bias resulting from fewer sites having higher weight than the others [83].Some guidelines for choosing the optimal λ are given in supplement S4.X is then computed as where I is an identity matrix of order n.Uncertainties in Y can be included in a matrix, C, and propagated into the system, in which case the inversion becomes a regularized weighted least squares inversion and X is determined as Common examples of C are the quality of PS pixels, spatial coherence, or InSAR LOS phase errors.Alternatively, spatial correlation in the subsidence map can be used, calculated via a semi-variogram as shown in Supplement S6.In this case study, we applied PS coherence as weights in C so that well sites closer to pixels with high coherence are preferred.Here Y consists of vertical ground motion, but measurements in all three dimensions can be included if available.After solving for X, the sites that better contribute to the interpolation have a greater weight.We then order the station sites by their weights to choose a subset of the sites as "best interpolation sites."The order of the stations is called "site rank." Defining Ŷ to be the interpolated subsidence obtained using the first k sites from V, we calculate the root mean square error (RMSE) between Y and Ŷ to estimate the interpolation quality.
Redistribution of the stations: Though the method gives the GNSS station sites that can better define the gradients in the deformation pattern, some of the sites that are located very near to each other get equal weight.In such cases, we remove some of the closely spaced sites to increase the overall spatial distribution and accuracy.To detect the sites located closely, we run the DBSCAN algorithm ( [84]; available as Python package) to automatically detect clusters of sites that are within a threshold distance (10 km) from each other.We then iteratively identify the best site in each cluster by higher weight or local site conditions and remove the other sites in the cluster.At this point, we check whether the removal of a site would significantly affect the interpolation quality.Such a situation may arise if closely located sites are needed to infer the spatial gradient of the ground deformation.We calculate the RMSE before and after removing a site and retain a site if the RMSE increases over 0.05 mm after the removal.The process also provides flexibility to choose from a set of closely spaced alternate sites if other considerations (e.g., access) preclude installing instruments at the selected location.Fig. 3(b) shows the final sites after redistribution, which are called the "best interpolation sites" assuming there is not an existing GNSS network in the area.

A. Accounting for an Existing GNSS Network
The method described in Section V-B assumes that there are no existing GNSS instruments in the study area.Southern Louisiana has existing GNSS stations (∼45 stations in the study area) maintained by C4G, and the new stations should enhance the existing network and monitor areas experiencing significant subsidence.To make sure that the new sites are not closely spaced to the existing sites, we also removed the sites within 10 km of a C4G GNSS station.We are then left with an initial pool of 346 (n) potential sites.
The existing stations are already collecting data and some of them are very close to the locations of potential future sites.Adding or removing any site can change the weights for the remaining sites from which optimal sites are selected.To account for the contribution of the existing GNSS stations, we consider that the existing stations are already collecting measurements, subtract their IDW interpolation contribution from the observed value (Y ) i.e., subsidence obtained from [71], and perform linear regression for best interpolation sites as described in Section V-B.Alternatively, one can assign the existing GNSS network locations a higher weight (x j ) than the new sites and then invert for the regression weights.This can be done by Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
adjusting λ such that weights for all the new sites are less than the weight assigned for the existing GNSS station locations.
After running the regression algorithm and redistributing by removing stations within 10 km of each other, we are left with 58 "best interpolation sites" for installing GNSS instruments [squares in Fig. 4(b)].
Choosing the right number of stations: Having more stations can always improve the subsidence estimation, yet choosing the "right" number of new sites for expansion is in practice influenced by a combination of technical, programmatic, and economic considerations.For the purposes of this discussion and since the scientific objective is to capture the spatial trend of the long-wavelength subsidence, we ignore implementationspecific considerations here, which will limit the number of stations that can be added practically.Those can be downselected based on the ranking of locations that best meet science objectives.Therefore, we limit the optimal number of new stations using the accuracy to which the subsidence can be interpolated.The RMSE of the fit to Y is plotted against the number of stations used in the interpolation.The plot shows the decrease in RMSE as more stations are added to the network [see Fig. 4(c)].The RMSE is not the total achievable InSAR accuracy, which depends upon many factors not just the tie point network.Rather, this value determines the contribution of the GNSS to the InSAR-derived VLM uncertainty when used to remove long-wavelength noise and tie the InSAR to a geodetic frame.
We observe the RMSE does not reduce significantly beyond 61 sites (16 new + 45 existing).Hence, we propose 16 sites as "optimal sites" for setting up new GNSS stations and enhancing the network.Fig. 4(b) shows the interpolated subsidence map from the 61 sites.The locations of the sites are listed in the supplement (Section S2).Plots in Supplement Section S3 show the variation in the interpolated subsidence with the number of sites used for interpolation.

B. Station Distribution Among Subsidence-Causing Factors
While the new GNSS sites are selected considering only their ability to capture the long-wavelength subsidence based on the current best measurement, it is valuable to also assess how well these stations can capture different subsidence mechanisms.To assess the network for differentiating processes driving subsidence, we evaluate the distribution of GNSS sites (existing + proposed) among multiple tectonic, geologic, and marsh types.We examine the distribution of proposed GNSS station sites among various geologic units [Holocene, Pleistocene, and Miocene-Fig.5 Forty-three percent (43%) [see Fig. 5(e)] of the sites are in the Pleistocene layers and are expected to only capture the deep subsidence, and the balance are in the varying thickness of the Holocene units.The deep subsidence measured at the GNSS stations footed in the Pleistocene can be combined with InSAR, which measures both deep and shallow components, to isolate shallow subsidence components in the entire Holocene layer and that part of the shallow subsidence measured by GNSS station footed in the Holocene units.This could also be used to measure the correlation between the thickness of the Holocene layer and the amount of shallow subsidence [44].
Subsidence due to soil layer compaction varies widely over different marsh types [85].Studies have shown that soils with more organic content and higher oxidation rates are susceptible to higher subsidence [86], [87], [88].In marshlands of coastal Louisiana, subsidence is reported in depths of the top 1 m in freshwater swamps to 14 m in thick peat deposits [89].The sites span multiple marsh types with two stations in freshwater marshes, three stations in intermediate marshes, one station in brackish marsh, and three stations in saline marshes [see Fig. 5(f)].The GNSS stations' anchorage depth can be varied in the marshes to detect the subsidence occurring at various depths, e.g., by installing several stations with different footing depths at a site.
Clayey soils typically exhibit shrink-swell behavior when in contact with water, making them vulnerable to compaction.Much of the study area is overlain by clayey soil units, but some parts along the Mississippi River channels are rich in terrigenous silt deposited by the river.Those areas experiencing consistent sediment deposition are subject to high overburden pressure.We define clayey soil with high organic content as peat.Our network is able to separate subsidence among the multiple soil types [see Fig. 5(g)].
In addition, numerous faults run through the area, and it has been proposed that they accommodate a significant portion of the deep subsidence [60], [69].The two major fault systems in the region are the Tepatate-Baton Rouge fault system and the Lake Hatch fault system, and the GNSS network stations are distributed between the upthrown and downthrown spans [see Fig. 5(d)] to identify the differential movements along the hanging and foot walls [see Fig. 5(h)].However, the number of stations south of the Lake Hatch fault system is limited by a lack of PS pixels in the Bird-foot delta.
Other than providing a better spatial distribution of stations in relationship to the faults, and with the exception of marsh type [see Fig. 5(b) and (f)] and peat soil [in Fig. 5(c) and (g)], the new stations do not significantly add to the diversity of the existing network.However, the weights (x j ) can be adjusted to prioritize site selection by different criteria if desired.Additional details regarding instrument type, data collection strategies, and ideal site conditions are discussed in Appendix A.

C. Application to InSAR CR/CAT Site Selection
The methodology discussed in Section V-B to determine the best sites for installing GNSS instruments can be extended to install other ground instruments such as extensometers, CRs, and CATs.InSAR CR/CAT are used to obtain InSAR measurements in areas that would otherwise have too low coherence for accurate measurements, thereby increasing the coverage of InSAR-based ground deformation measurements [26].
The Bird-foot delta in southern Louisiana is underlain by thick Holocene sediments [see Fig. 1(b)] and has been identified as Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.having high rates of subsidence [53].However, our method has identified very few GNSS site locations (2 sites-Fig.4) in the area even though InSAR measurements show a strong subsidence signal [see Fig. 6(a)].This is mainly due to the fact that our initial pool of possible sites is limited by the requirement to have an InSAR PS point within at least 5 km radius.Dense measurements are needed to bring out the sharp gradient of subsidence in the lower Bird-foot delta area.Vegetation, soil moisture, periodic inundation, and the lack of built structures limit the number of PS in the area, necessitating the use of CR/CAT to obtain quality InSAR returns.Placement of this instrumentation can be chosen using a simplified version of the method described in Section V-B.
Since we do not have predefined sites for CR, we divide the area of interest into uniform 5 km × 5 km grids and determine the best grids to install CR.We implement the procedure described in Section V-B by subtracting the IDW contribution of measurement at PS locations [black "+"s in Fig. 6(b)] to the interpolated estimate (G ps ) from Y , and (4) can be rewritten as Here in G, we use the locations of grid centers and solve for their weights in X.The grids with higher weight are then presumed to be the best for setting up CR.To avoid redundancy, we did not consider grids within 5 km of a PS.Clearly, more measurements are needed in the lower regions of Bird-foot delta to better infer the subsidence [see Fig. 6(b)].
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

D. Influence of PS Estimation and Deformation Interpolation Methods
In our study, we used coherence and amplitude dispersion as criteria to identify stable InSAR measurements, and we then chose well sites closer to them.However, some of the properties inherent to the methods of coherence and amplitude dispersion estimation can influence the PS pixel selection.Recent advances in the selection of stable scatterers from InSAR data [81], [90], can result in more PS pixels than the ones used in the study.We also expect that this is the case when NISAR and ROSE-L satellites become operational and the abundant L-band InSAR data can offer better estimation of stable scatters.More stable scatterers may increase n.
The weights calculated for well sites in X are also dependent on the interpolation method used in (4b).We used a nongeostatistical method IDW, which solely relies on the distance between the well-site and the measurement point, whereas a geostatistical method such as kriging [82] can consider the relative distribution of the well-sites and the spatial autocorrelation in measurements between the well-sites.Hence, using a different interpolation technique can result in a different X vector.In the supplement S5, we illustrate the effect of the interpolation method on X and the resulting RMSE.Studies have shown that kriging can be a better interpolation technique, but its implementation involves many parameters estimating the semi-variogram function which can vary from study area to study area [91].So, when implementing the procedure given in this study, one can vary the interpolation method to determine the effect of the choice.

VI. SUMMARY
This study provides a methodology for identifying the best locations for installing GNSS stations in areas with heterogeneous long-wavelength subsidence that must be distinguished from long-wavelength atmospheric noise.In many areas where coastal subsidence measurements are needed, InSAR analysis is challenging not only because of severe tropospheric noise but also temporal decorrelation.In these areas, the GNSS network can tie together spatially disjoint InSAR ground displacement measurements to improve surface deformation estimates and avoid phase unwrapping issues associated with conventional InSAR processing in rapidly decorrelating areas.Our method uses InSAR coherence and amplitude dispersion measurements as a guide to identify locations of possible future measurements.
As a test case, our method is applied to design an expanded GNSS network suitable for southern Louisiana.The algorithm design is motivated by the anticipated L-band regular radar acquisitions from the NISAR mission (expected to launch in 2024) to enable better use of InSAR in low-coherence areas, specifically coastal areas facing permanent land loss through subsidence combined with sea level rise.The GNSS network can serve as ground truth for future InSAR measurements and can be used to connect spatially relative InSAR measurements to a geodetic reference frame.Studies have shown that installing many GNSS stations as a network is cost-effective.The study can also serve as a guide to incorporate site accessibility and reliability considerations into setting up a GNSS station and can be applied to similar locations worldwide.
Our study provides a methodology to optimally sample the ground deformation based on a priori knowledge as opposed to uniform sampling, which might not adequately sample highgradient areas.Even with long wavelength error corrections, GNSS and InSAR do not measure exactly the same ground displacement because InSAR also detects that part of the movement that occurs above the GNSS foundation depth.While a GNSS network with many stations spread across a broad area is insensitive to localized subsidence at a few stations, in cases where shallow subsidence (above the foundation) is spatially correlated over large areas and affects a large percentage of the sites in the network, measurements of shallow subsidence are needed.Section III describes several ways in which those shallow subsidence measurements can be made.
L-Band SARs are less subject to InSAR decorrelation than C-band SAR, the current standard given free and open access to Sentinel-1 data.Therefore, they can provide better subsidence measurements over vegetated areas because there will be more PS points, particularly in areas like Louisiana [92].Better subsidence measurements tied to a geodetic frame will yield land surface elevation measurements that can be updated more frequently, informing disaster management and coastal restoration projects.Furthermore, improved spatial and temporal resolution of subsidence measurements can help differentiate between the

Fig. 1 .
Fig. 1.(a) Sentinel-1 data location and the existing network of GNSS and RSET-MH stations in Louisiana.Blue dashed rectangles show the footprints of Sentinel-1 SAR images corresponding to tracks (P165F90-95; P63F92; P136F93).Dark-red dots indicate the GNSS CORS station network in Louisiana operated by C4G; red diamonds show the coastal RSET-MH station network operated by the State of Louisiana; grey dotted rectangle shows the extent for b.(b) Depth of Holocene-Pleistocene interface in the study area (from [52]) and the locations of the areas discussed in the text.(c) Illustration depicting the subsidence in the study area and various monitoring mechanisms.

Fig. 2 .
Fig. 2. InSAR pixels with (a) spatial coherence exceeding 0.6, and (b) amplitude dispersion below 0.4 for four Sentinel-1 SAR scenes in southern Louisiana (plotted at pixel centers).(c) PS points in the study area.A total of 127443 PS points over four InSAR data frames covering southern Louisiana and 465 viable well sites (within 5 km radius of PS pixels) are identified for solving for the optimal locations of future GNSS stations.The PS distribution only encompasses 0.3% of the total pixels (∼46 million), with most concentrated in the urban New Orleans area.

Fig. 3 .
Fig. 3. (a) Location of orphaned oil well sites (712) in southern Louisiana and well-sites with at least one PS point in a 5 km radius (465); PS points are from Fig. 2(c).(b) Best interpolation well sites (75) selected after redistribution assuming there is no prior GNSS network.

Fig. 4 .
Fig. 4. (a) VLM rates for ∼4 years of InSAR data (Jan.2017-Dec.2020) spanning four Sentinel-1 tracks (P165F90-95; P63F92; P136F93) derived using Eigen decomposition-based phase reconstruction in [71].Dark red dots indicate the existing GNSS CORS station network in southern Louisiana operated by C4G.(b) Interpolated subsidence map for the study area using the 61 best sites for GNSS CORS station installation (orange squares).(c) RMSE reduction with GNSS installation at best sites.The first 45 are the existing C4G stations and the red dotted line indicates 61 stations, beyond which the RMSE reduction is not significant.

Fig. 5 .
Fig. 5. Distribution of the GNSS station sites among (a) geologic units; (b) marsh types; (c) soil types; and (d) around major fault systems.Histograms showing the number of GNSS stations in each category of (e) geologic units; (f) marsh types and (g) soil types.(h) Histogram showing the number of sites located on both upthrown (North) and downthrown (South) sides of Tepetate-Baton Rouge fault system and Lake Hatch fault system.
Optimized GNSS Cal/Val Site Selection for Expanding InSAR Viability in Areas With Low Phase Coherence: A Case Study for Southern Louisiana Bhuvan K. Varugu , Cathleen E. Jones , Member, IEEE, Ke Wang , Student Member, IEEE, Jingyi Chen , Senior Member, IEEE, Randy L. Osborne , and George Z.Voyiadjis