Determination of the Earthquake Epicenter from line-of-sight displacement images obtained by Sentinel-1A/B radar data using the GMTSAR software package

. This article described the technology of determining earthquake epicenter with radar remote sensing on the example of Sentinel-1A/B. To determine the epicenter of the earthquake, the Earth's crust displacements were analyzed using radar remote sensing data obtained for the ascending and descending flight orbits. Coordinates of Earthquake epicenters were found according to line-of-sight displacement images via its maximum value. Displacement of the Earth's crust was obtained by processing in the GMTSAR package in the VirtualBox virtual machine of the Linux Ubuntu 16.04 operation system. Two earthquakes that occurred in 2020 were studied to determine the accuracy of finding epicenters using the ascending and descending orbits Sentinel-1A/B. These earthquakes occurred in Western Xizang, China, and Doganyol, Turkey. The maximum deviation from the officially registered epicenter coordinates was 15.38 km for Doganyol and 3.2 km for the Western Xizang earthquake. The negative displacement was 90 mm for Doganyol and 50 mm for Western Xizang.


Introduction
Based on the difference in phase images obtained from remote sensing radar data, it is possible to detect an earthquake's epicenter visually. The results of radar image processing show that the location of the epicenter is visible in the form of repeated bar spectra in the color palette of the phase image [1,2]. The pattern that has a rainbow pattern characteristic of an earthquake is concentrated in the epicenter zone in the form of concentric phase transitions, in the center of which the earthquake center can be assumed, and it is proportional to the magnitude of the earthquake. Sentinel-1A/B radar remote sensing satellites launched in 2014 can be used to determine these earthquake epicenters [3,4]. The main feature of these radar remote sensing data; first, they can be used to detect an earthquake's epicenter with greater coverage of the territory than others, and second, they are free for end-users. Monitoring of a single territory can be repeated every 12 days with the same view or angle geometry. The epicenter of an earthquake is easier to detect than the phenomena caused by landslides and subsidence of the earth's surface [5][6][7]. The received radar images can be processed by a special software, such as ROI_PAC, GMTSAR, GAMMA, ENVI Sarscape, ESA SNAP [8][9][10][11]. Some commercial software is licensed, has for more than $ 10,000, and may not be available to the end-user. In contrast, the GMTSAR package, a set of additional GMT scripts developed by David Sandwell and Xiaohua Xu of the University of California, San Diego, can also be successfully used to determine the epicenter of an earthquake. In the future, we will use it to determine the location of the earthquake epicenter and show an algorithm for processing remote sensing radar data.

Study area and data
Two earthquakes are being investigated. The first of them occurred in Turkey near Doganyol on January 24, 2020. The second was in China, western Xizang, in the same year on July 22. They were with magnitudes M6. 7 and 6.3. Sentinel-1 radar remote sensing data before and after the earthquake and DEMfile were used to create the phase difference of VV polarization channel (see Figure 3). Sentinel-1A/B radar remote sensing satellites operate at a wavelength of approximately 5.6 cm. Sentinel-1A and Sentinel-1B radar remote sensing data were used with ascending and descending orbits for two earthquakes (see Table 1). All this data was downloaded from the Copernicus Open Access Hub (https://scihub.copernicus.eu/). Data were obtained for the territory of Turkey and Western Xizang, as shown in Figure 1 03.02.2020 S1B_IW_SLC__1SDV_20200203T032513_20200203T032540_020 099_0260A8_412D Figure 1. The footprint of selected Sentinel-1A master scenes for ascending and descending orbits.

Radar image processing
Radar images are processed using the GMTSAR package. Radar image processing involves the sequential execution of scripts in Fortran and C, one of the most important is the procedure of joint registration of each sub-swath IW1, IW2, IW3 and the formation of an interferogram from Sentinel-1 data [13]. The output products of interferogram formation are the phase difference (see Figure 3). The phase difference is calculated as follows for each pixel for a pair of jointly co-registered radar images: where 1 , 2the real part of a complex number from a pair of radar images, and 1 , 2the imaginary part of the complex number of each pixel of these radar images [14]. The phase difference is the sum of several components associated with the displacement ∆ , related by topographic phase ∆ , by atmospheric phase delay ∆ and noise effect ∆ [15]: ∆ = ∆ + ∆ + ∆ + ∆ An earthquake causes a displacement of the earth's surface, where the phase associated with the displacement is significantly larger than the phase associated with the atmosphere and noise ∆ ≫ ∆ + ∆ . Knowing this, we can conclude that sum of atmospheric and noise phase effects is negligible compared to what is caused by displacement from an earthquake. Therefore, the phase associated with the offset is calculated as subtracting the topographic phase from the phase difference: ∆ = ∆ − ∆ Subtraction of the topographic phase is used to calculate the phase associated with displacement in radar image processing. The topographic phase depends on the spatial baseline ⊥ [16] (See Figure 4). The topographic phase is calculated as follows: The procedure for subtracting the topographic phase is performed using a digital elevation model (DEM) and phase difference image. The DEM image contains the height value ℎ. The DEM image is used to calculate the topographic phase component ∆ [17][18][19]. The phase associated with displacement can be distorted by noise. Noise can occur due to temporal decorrelation, geometric decorrelation, volume scattering, and processing error. The Goldstein filter is used to eliminate or partially reduce the influence of the phase noise effect [20][21][22][23][24]. To generate an interferogram in the GTMSAR package, we need to run the script "p2p_S1_TOPS_Frame.csh", the command line looks like this for the first pair of the VV polarization interferogram for July 18 and July 30, 2020.
p2p_S1_TOPS_Frame.csh S1A_IW_SLC__1SDV_20200718T121528_20200718T121555_033509_03E207_57A0.SAFE S1A_OPER_AUX_POEORB_OPOD_20200807T120800_V20200717T225942_20200719T005942.EOF S1A_IW_SLC__1SDV_20200730T121529_20200730T121556_033684_03E765_141D.SAFE S1A_OPER_AUX_POEORB_OPOD_20200819T121202_V20200729T225942_20200731T005942.EOF config.s1a.txt vv 0 Before running the script, two folders called "raw" and "topo" should be prepared. The raw folder must contain Sentinel-1 satellite data in the directory with the name*. SAFE at the end. The topo folder must contain a DEM file generated from https://topex.ucsd.edu/gmtsar/demgen/.
The execution time of the "p2p_s1_tops_frame.csh" script depends on the computer's performance and takes the longest time, unlike other scripts. In general, the duration is 2 hours in VirtualBox Ubuntu 16.04 on an Intel® Core™ i7 -8565U 1.8 GHz processor with RAM 8 GB.
As a result of executing the script "p2p_s1_tops_frame. csh", three folders F1, F2, and F3, are created for each sub-swath IW1, IW2, and IW3. There are corresponding interferograms in those directories. In addition to these directories, the "merge" directory is formed, where the result of gluing interferograms of all sub-swaths are located (See Figure 3).
The epicenter of an earthquake exists in one of the sub-bands, and it is further processed to calculate the absolute unwrapped phase. The location of the earthquake epicenter is visible in the filtered phase difference images associated with the displacement, as a rainbow pattern periodically repeated in the color palette, their concentration consists of a large number of phase transitions from -π to +π (See Figure 5). Using the script " snaphu. csh .12 " in the GMTSAR package, the phase difference image is converted to an absolute unwrapped phase, the program code is developed by Chen and Zebker [25][26][27]. This script is written in C and must be executed separately.

snaphu.csh .12
The absolute unwrapped phase difference image has a negative part -8 and a positive part up to 12 radians in Figure 6. The script ".12" value means that the part of the image that has coherence less than 0.12 will be masked and will not be subjected to absolute phase calculation. Figure 6. The absolute unwrapped phase difference of IW2 sub-swath Sentinel-1A with an ascending orbit for the Xizang earthquake on July 18 and July 30, 2020 The following script is " geocode. csh .12" calculates the line-of-sight displacement image (See Figure  7) and georeferences all data from radar coordinates (range, azimuth) to latitude and longitude geographic coordinates, along with kml-files for Google Earth.

geocode.csh .12
The line-of-sight displacement radar is the difference between the distances before and after an earthquake, the shortest way a radio wave propagates from the antenna to the earth's surface and back. Displacement on line-of-sight may differ depending on ascending (asc) and descending (desc) orbits of Sentinel-1A/B satellite , . These displacements at each point of the earth's surface depend on the radio waves' angle of incidence for the ascending and descending orbits , , respectively. The essence of line-of-sight displacement caused by the Earth's crust movement is illustrated in Figure 8 when deforming up.   The accuracy of determining the epicenters of earthquakes was carried out according to the US Geological Survey's official data using a priori information about the coordinates of the epicenters. Western Xizang: 33.144°N 86.864°E; Doganyol: 38.431°N 39.061°E. The distance between the official epicenter of the earthquake and the data obtained after processing the Sentinel-1A/B radar remote sensing data is presented in Table 3 and Figure 12. The distance between the near epicenters S1A and S1B with ascending and descending orbits is 0.32 km for the Doganyol earthquake. 3 км Note: S1A -Sentinel-1A, S1B -Sentinel-1B radar remote sensing satellites Changes in displacement values from negative to positive were considered for the Western Xizang and Doganyol earthquakes (see Table 4). Negative values were used to determine the maximum changes that characterize the uplift of the Earth's crust. The absolute value displacement along the direct line-of-sight, which has negative values for Doganyol, is 90 mm, for the Western Xizang Earthquake is 50 mm. These values are completely consistent with the magnitude of the earthquake, which shows that the larger the magnitude, the greater the displacement we observe. The Doganyol earthquake magnitude of 6.7 M is greater than Western Xizang of 6.3 M.

Conclusion
Sentinel-1A/B radar remote sensing data is applicable for determining earthquake epicenter. The maximum error of deviation from the official one is 15.38 km for the Doganyol earthquake. The distance between the ascending orbits for the Doganyol earthquake is up to 0.32 km. A different geometric view with different orbits (ascending and descending) leads to a greater deviation of the distance error. In the case of the Doganyol earthquake, we have volumetric scattering from vegetation; therefore, the phase is with noise in some places than in Xizang. Western Xizang has a distance error of 3.2 km, which is better than the Doganyol case.
Determining the earthquake epicenter location associated with displacement, we can find the maximum value (in millimeters) of pixels to displacement in line-of-sight on the resultant image. This method will determine earthquake epicenter using modulo the maximum displacement along the direct line of sight of the radar. As the results show, the displacement values along the line of sight (see Table 4) are directly proportional to earthquake magnitude. Doganyol Earthquake with magnitude 6.7 M has a line-of-sight displacement -90 mm, Western Xinjiang with a magnitude 6.3 M -50 mm.
Despite the GMTSAR package's availability, a drawback is the lack of a user-friendly interface in the ESA SNAP software. There are inconveniences, and they are associated with the manual launch of GMTSAR scripts and pre-preparation of data. This GMTSAR package is more intended for experts than ordinary users who do not have enough experience with commands in the Linux terminal. Radar image processing in GMTSAR to create displacement images and determine the maximum displacement coordinates takes less elapsed time than ESA SNAP. The GMTSAR package is undemanding to use more powerful computers than other software on Windows operation system. The original Sentinel-1A/B data takes up 14 gigabytes of data, and it can increase several times due to data processing, so more capacious memory drives are needed. Therefore, more productive computing systems are required to work with such data arrays in some cases. We need a supercomputer with more RAM and high-speed SSD-based hard drives for fast, automatic radar image processing, a program code that implements the entire algorithm as in Figure 9, and a friendly interface for interacting with radar remote sensing data.