1 Introduction

The most recent nuclear tests in the world were announced by the Democratic People Republic of Korea (DPRK) having been carried out on October 9th, 2006 and May 25th, 2009. In the following, these tests will be called DPRK06 and DPRK09, respectively. The clear seismological signals detected by many stations of the globe clearly characterized these tests as underground explosions. The analysis of the seismic waveforms allowed the location and magnitude estimation for both events. For the 2006 test a body wave magnitude m b = 4.0 ± 0.1 was obtained, while the relatively high long-period seismic noise covered the surface waves, so that only an upper limit of 3.5 was estimated for M s . The higher amplitudes of the recorded seismic signals for the 2009 event lead to the estimation of 4.5 ± 0.1 and 3.2 ± 0.2 for m b and M s, respectively.

The magnitudes estimated for the two events by the international agencies are fairly well consistent with ours. The m b magnitude ranged from 4.1 (IDC/REBFootnote 1) to 4.2 (USGS/NEICFootnote 2) for DPRK06 and from 4.5 (IDC/REB) to 4.7 (USGS/NEIC) for DPRK09. The surface-wave magnitude Ms for the explosion, determined from the IDC REB based on 15 stations, was 3.6 (significantly larger than expected for an explosion with that small m b). Also Murphy et al., (2010 ) noted somewhat unusually large amplitudes for surface waves of both events, compared to the historical explosion surface wave measurements, and they justified this peculiarity with the circumstance that the explosions were small and embedded in high-velocity hard rock.

The epicenters of both tests were located in a mountainous area of scarce natural seismicity, geologically belonging to a tectonic region constituted by granitic massifs and not far from volcanic structures, some of which still active at the present time (Fig. 1). For a conversion of the m b magnitude into yield, numerous relations developed in the past are available. The proper relation to be used in specific cases, like this, depends on an accurate knowledge of the geological framework. Unfortunately, this is not our case, due to the vicinity of granitic intrusions to lava flows and pyroclastic material coming from the volcanoes present in the area. As an example, adopting the m b values reported above and m b yield relations obtained for the former Soviet Republic (Ringdal et al., 1992) and the Nevada test site (Murphy 1981), the possible values of yield range between 0.25 and 1.0 kt for the 2006 test and between 1 and 5 kt for the 2009 test.

Fig. 1
figure 1

Geological map of the DPRK region where the 2006 and 2009 tests were carried out. Red mesozoic granite, pink precambrian rock, purple quaternary lava. The white star shows the location of the test sites

Full size image

In the present study, the DPRK09 test has been investigated also by applying DInSAR (Synthetic Aperture Radar Differential Interferometry) technique. Since its first use in 1992 (Massonnet et al., 1993) DInSAR is nowadays considered an effective tool in Earth Sciences studies able to detect centimetric and/or millimetric surface movements over large areas. Recent studies report that although earthquakes of magnitude <4.8 are unlikely to be observable, coseismic surface deformations induced by very shallow events are detectable (Dawson 2007; Dawson et al., 2008). In recent years the capabilities of DInSAR were improved in terms of spatial resolution, accuracies and revisiting time. Indeed, the new generation of high-resolution (HR) and very HR (VHR) SAR sensors is available since 2007, with the Japanese ALOS PALSAR, the Canadian Radarsat-2, and up to the most recent sensors, like the VHR German TerraSAR-X and the Italian COSMO-SkyMed constellation, both achieving 1 m of spatial resolution. The DPRK09 is the first ever nuclear test where VHR DInSAR has been applied.

Previous studies concerning the co-seismic surface deformation signal due to underground nuclear tests have been conducted at the Nevada Test Site in 1992, where SAR data were collected over a 14-month time span to cover three tests (Vincent et al., 2003; Vincent 2005). These studies reported a subsidence of a limited area above the detonation point as common co-seismic feature of the underground nuclear explosions.

The present study has exploited the capabilities of DInSAR applied to a dataset consisting of three TerraSAR-X images acquired in Stripmap Mode (at 3 m spatial resolution). The DPRK09 case study has encountered some difficulties which were not present in the former Nevada cases. First, as written above, the location of the nuclear test is known approximately only within the typical uncertainty of seismic location. Second, the topography of the investigated area is far from the flat plate of Nevada test sites. This implies a more complex scenario needing a digital elevation model (DEM) to fully remove the topographic phase. Moreover, the use of three SAR data is expected at least to make some cross-analysis concerning the presence of the atmospheric phase (in particular, the wet troposphere) and possible topographic residual phase.

2 Single Locations and Systematic Travel Time Residuals

In the context of verifying compliance of any nuclear test ban treaty, the step immediately following the detection of a waveform event is its location. We collected P wave recordings from as many as possible stations belonging to international monitoring networks, and picked the arrival times.

In the framework of the International Scientific Studies (ISS09), an initiative launched by the Comprehensive Nuclear Test-Ban-Treaty Organization (CTBTO) to stimulate the scientific community to improve the Treaty verification, we evaluated some data collected by the International Monitoring System (IMS) established by the Preparatory Commission of CTBTO. Moreover, to evaluate the potential of the Geotool software package developed within the Organization, we picked various phases from the waveforms using Geotool as well. Table 1 reports the arrival times picked up for both the 2006 and 2009 tests. As reported in Fig. 2, it clearly appears that these two events were recorded at mostly, or even only, teleseismic distances.

Table 1 Arrival times picked up for the 2006 and 2009 DPRK tests
Full size table
Fig. 2
figure 2

Stations used for event location

Full size image

We located the events using a least-squares location algorithm developed at the INGV. In this algorithm epicentral distances are computed from geographical coordinates through the WGS84 ellipsoid model of the Earth, and travel-times are based on the IASPEI91 tables (Kenneth and Engdahl 1991). Moreover, travel-times are corrected both for the ellipticity of the Earth by the formulation of Dziewonski and Gilbert (1976), and for the station elevation. Our tests showed that these two corrections, even if they are of the order of only few tenths of second, reduce the RMS of the time residuals at the recording stations. We didn’t apply any other static station correction nor specific source site correction because such information is unavailable at the present status.

In developing our location algorithm, particular care has been devoted to the computation of the parameters (semiaxis sizes and azimuth) of the epicenter error ellipses. We applied the theoretical framework available in literature for this problem (e.g. Stein and Wysession 2002). Moreover, we carried out numerous tests with a Monte Carlo method, by simulating the location of up to 1,000 events. These locations were obtained starting from synthetic arrival times from a given epicenter, and changing these arrival times in random way with the same standard deviation of the real observation data. The error ellipses statistically obtained from the synthetic epicenters match quite well the theoretical one obtained in connection with the original location.

Our algorithm has been implemented in a FORTRAN 90 code. The locations of the two events using all the available arrival times, and having assigned a zero depth to both of them, are listed in the first two rows of Table 2 and shown in the map of Fig. 3.

Table 2 Epicentral coordinates and error ellipse parameters obtained for the 2006 and 2009 DPRK tests using different data sets
Full size table
Fig. 3
figure 3

Epicentral locations of the 2006 (a) and 2009 (b) DPRK tests using the different data sets reported in Table 2. The epicenters are tagged with the same ID numbers as in Table 2. The error ellipses are at a 95 % confidence level. The green marks point to the preferred locations reported by NEDB (see Sect. 5)

Full size image

Table 2 (lines 1 and 2) and Fig. 3 (orange symbols) show a modest difference between the epicenters of the two events, obtained with the arrival times picked on all the available waveforms. In fact the 2006 event (eight stations used) is located nearly 7 km to the south west of the 2009 event (17 stations used). However, the error ellipses at the 95 % confidence level of the two events partly overlap. We can easily show that this distance is mostly due to the use of different data sets, associated with systematic station-depending discrepancies existing between the real travel-times in the Earth and the theoretical travel-times used in our location algorithm.

By repeating the locations using only seven of the eight common stations for the two events (bold phase in Table 1) the new epicenters, listed in the third and fourth row of Table 2 and showed in Fig. 3 (red symbols), are both located near the epicenter of the 2006 event obtained using all the eight available arrival times, and much closer to each other. In fact, in this case the 2006 event appears located only 3 km to the east of the 2009 event. This distance is contained within the error ellipses at the 95 % confidence level. In this exercise, station MKAR was excluded due to its uncertain time pickings.

We repeated the locations by means of a set of P wave arrival times reported in the ISC bulletins. In this case the number of stations used was 27 and 69 for DPRK06 and DPRK09, respectively. The results reported in Table 2 (lines 5 and 6) and Fig. 3 (blue symbols) show epicenters shifted a few kilometers to the north-east of the previous locations and smaller error ellipses. For sake of completeness, in Table 2 (lines 7 and 8) and in Fig. 3 (purple symbols) we report also the locations obtained by the ISC. These epicenters are very close and their error ellipses are comparable to those reported at lines 5 and 6 of Table 2, which were obtained with a smaller set of arrival times.

Based on the above results, we proceeded to study in more detail the systematic residuals affecting the arrival times at the seven common stations listed in Table 1.

We used a computer code (in the following named wave-shifter) specifically developed at the INGV for comparing the relative locations of two or more seismic events. The input for this code consists in the location coordinates, depth and origin time of the events to be compared, besides the digital waveforms recorded at a number of stations from these events, and outputs a plot of such waveforms referred to a time scale the origin of which is, for each single station, the theoretical arrival time computed through the IASPEIFootnote 391 travel-times tables.

Figure 4 shows the output of wave-shifter for the seven common stations used for obtaining the first two epicentral coordinates reported in Table 2. The plots of Fig. 4 show that:

Fig. 4
figure 4

P arrival waveforms plotted by program wave-shifter for K06 (red) and K09 (blue) events. The zero-value on the x-axis coincides with the expected arrival time computed for both events at each of the seven common stations, from their respective locations and origin times obtained through all the available data

Full size image
  1. (a)

    The first onsets of the P wavelets at some stations are variably shifted by several tenths of second with respect to the computed arrival times, pointing out the existence of errors in time pickings and/or systematic differences between the IASPEI91 model and the real travel-times in the Earth;

  2. (b)

    There are also significant shifts between the first onsets of DPRK06 (red line) and those of DPRK09, meaning that at least one or both events are mislocated.

The difference between the two epicenters is greatly reduced, if we use the arrival times of the same stations for both events. As a further test, we applied the wave-shifter program by using the epicentral location obtained by the seven common stations as reported in the third and fourth rows of Table 2. The results are plotted in Fig. 5. Now, the first onsets for DPRK09 are much closer to those of DPRK06. Fairly surprisingly, we note a systematic shift between the first onsets for the two events. More precisely, the wavelets from DPRK06 seem to arrive about 0.2–0.3 s earlier than those from DPRK09 at all the stations, the same that we would note delaying the origin time of DPRK06 by that time interval. A possible explanation for this circumstance is that the arrival times picked up for DPRK06 and used for its location (Table 1), due to the lower signal-to-noise level, have been read by the analyst with an average delay of 0.2–0.3 s relatively to those of the other event.

Fig. 5
figure 5

As in Fig. 4, but using the locations and origin times obtained through the data of the seven common stations only

Full size image

If we applied a correction of this size to the computed origin time of DPRK06, the overlapping of the waveforms of the two events shown in Fig. 5 would be almost perfect. It confirms that the relative location of these epicenters, as reported in the third and fourth rows of Table 2, is very close to the true one.

We may confidently assume that, among our results, the epicentral locations obtained by the largest number of arrival times, reported in the fifth and sixth lines of Table 2, are the most reliable. Nevertheless, based on these data only, and without ground truth information, we ignore the actual size of the systematic mislocation that could still affect the results.

The satellite-based approach introduced earlier and described in more detail later in the following sections is aimed at establishing the missing ground truth.

3 Relative Locations

In the previous section we showed that the plots obtained by our wave-shifter program allow considerations about the relative location of DPRK06 and DPRK09. In the following, we shall show that it is possible to obtain a reliable relative location by a trial-and-error procedure, changing the epicentral coordinates and origin time one by one, until a satisfying overlapping of the wavelets at all the stations is obtained. However, this method requires considerable time and skill. At first, we simply tested the hypothesis that the two tests were carried out exactly on the same point; therefore, we applied the wave-shifter program, giving the same epicentral coordinates to both of them. For this exercise, we adopted the coordinates and the origin time of the DPRK09 event reported in the fourth line of Table 2. Because the origin time to be associated to the simulated DPRK06 epicenter was unknown, we proceeded tentatively changing the origin time until a very good coincidence of the two wavelets at station PDAR was obtained. The output of the wave-shifter program is shown in Fig. 6.

Fig. 6
figure 6

As in Fig. 4, but assuming identical locations for DPRK06 and DPRK09, and assigning the origin time to DPRK06 in order to obtain a perfect coincidence of the two wavelets at station PDAR

Full size image

Figure 6 shows that not only stations NVAR and PDAR (in Northern America), but also ASAR and WRA (in Australia) exhibit good correlations. This indicates that the distance between the real epicenters of the two considered events and each of these four stations is respectively about the same. This is not the case for the other three stations AKASG, FINES and GERES, all of them located in the Euro-Asian continent, to the northwest of Northern Korea. The time shift observed for these stations is of the order of 0.2 s, with the signals from DPRK06 (red line) arriving later than those from DPRK09. From this exercise we may conclude that the two tests were not carried out on exactly the same point, and the epicenter of DPRK06 is probably some kilometers to the east-southeast of that of DPRK09.

For relative location of two or more events, the double-difference joint hypocenter determination (DDJHD) method (Waldhauser and Ellsworth 2000) has become very popular in the last decade for high-resolution imaging of clustered seismicity in active areas. At the INGV, we developed our own DDJHD algorithm (Console and Giuntini 2006). We applied it in a global environment, again making use of the IASPEI91 travel-time tables. For both the local and global scales, the DDJHD method is based on the principle that the hypocenters to be located relative to each other are closely spaced in comparison with the distance between the hypocenters themselves and the recording stations. The advantage of the DDJHD method consists in removing the influence of the systematic travel time residuals from the location process, which, as we have seen before, may be relevant when the different seismic events are located using different sets of stations.

The DDJHD method can be applied to the P wave arrival times picked up by the analyst on the seismograms, but its accuracy can be significantly improved by applying a cross-correlation technique to a suitable segment of the P waveforms. In such a way, the accuracy of the results is not limited by the skill of the analyst in subjectively reading the first onsets of the P arrivals, but is objectively assured by the automatic procedure carried out by the computer in finding out the time by which one segment must be shifted with respect to the other in order to yield the maximum correlation coefficient between these waveform segments. Moreover, the information used by this method comes from a waveform segment of suitable length, and not only from a single sample where the first onset is detected. Nevertheless, according to our experience, the use of the method requires some expertise to assess, for example, the most appropriate filtering band and the length of the waveform segments that must be correlated. Multiple maxima of the correlation function of similar size are a common problem faced when processing the signals recorded from the two considered events. In our algorithm, implemented in a MATLAB code, the solution is obtained through a simple least-squares best fit.

The application of the DDJHD algorithm to the waveforms of the seven common stations has provided the relative location shown in Fig. 7, which is fairly well consistent with what was inferred from the application of the wave-shifter program in the test described above.

Fig. 7
figure 7

Relative location of the two Northern Korea tests obtained from the DDJHD algorithm. The origin of the coordinates has been arbitrarily put on the epicenter of DPRK09

Full size image

The RMS of the time shifts residuals obtained by the best fit algorithm is equal to 0.167 s, a small value in comparison with the RMS obtained from all the locations reported in Table 2. This confirms the advantage of using the DDJHD algorithm for relative location of seismic events.

For testing the quality of our DDJHD solution, we applied once again the wave-shifter program. In this case we assume that the coordinates of DPRK09 are those reported in the fourth row of Table 2, and the coordinates of DPRK06 are obtained from the former by shifting them by 2.55 km to the south and 1.45 km to the east (see Fig. 7). The origin time of DPRK06 is obtained by the trial-and-error procedure described above in order to obtain a perfect coincidence of the two wavelets at station PDAR. The output of the wave-shifter program is shown in Fig. 8.

Fig. 8
figure 8

As in Fig. 4, but assuming for DPRK06 the location and origin time obtained through the DDJHD algorithm relatively to DPRK09

Full size image

Figure 8 shows that at all the stations the P waves arrive about 0.1 s earlier than expected from the theoretical travel-times, except for station ASAR and WRA, where the actual arrivals are about 0.1 s later. This means that our DDJHD locations are still slightly different from the true ones. As the time differences obtained from the cross-correlation technique are certainly accurate enough, a possible explanation of the misfit observed in Fig. 8 is that the geological structure in the test area is inhomogeneous, representing a violation of the requirements for the application of the DDJHD method.

In order to correct the misfit still existing in Fig. 8 with the DDJHD location, we used the wave-shifter program with a procedure of trial-and-error, as explained above. By visual inspection of the plot in Fig. 8, we found out that to improve the relative location of the two events, the epicenter of DPRK06 has to be shifted to the North. Of course, in addition to the shift of the epicenter, it is also necessary to adjust the origin time of DPRK06 to maintain a perfect overlapping of the two wavelets at station PDAR. After a number of trials, we considered acceptable the solution shown in the map of Fig. 9. The related wave-shifter plot is shown in Fig. 10.

Fig. 9
figure 9

Relative location of the two DPRK tests obtained adjusting the epicentral coordinated of DPRK06 for the maximum correlation of its waveforms with those of DPRK09. The origin of the coordinates has been arbitrarily put on the epicenter of DPRK09

Full size image
Fig. 10
figure 10

As in Fig. 4, but assuming for DPRK06 the location and origin time obtained adjusting the epicentral coordinated of DPRK06 for the maximum correlation of its waveforms with those of DPRK09

Full size image

4 Analysis of Satellite Data

Remote sensing methods have demonstrated their ability to support the localization of underground nuclear test sites. Both optical and synthetic aperture radar (SAR) data can be used for this purpose (Canty et al., 2009). In particular, in the last decade, optical satellite sensors have increased the resolution of the data imagery, reaching sub-meter per pixel resolution, thus giving a new opportunity for better identify changes directly induced by nuclear test (e.g., damage or strong surface modifications), or indirectly by human activities during site preparation before the nuclear test (e.g., new buildings or roads and spreading of material in the surrounding). Mainly, the change detection approach and visual comparison methods (Schlittenhardt et al., 2010 and references therein) have been applied on optical data to locate area affected by these types of changes in the proximity of nuclear test sites.

Beside optical imagery, SAR data can successfully exploited. SAR are active radar imaging sensors, working on the microwave region of the electromagnetic spectrum and thus they can operate in almost all weather conditions and during day and night time.

4.1 InSAR Rationale

The SAR is a radar imaging sensor exploiting satellite orbit paths to achieve a spatial resolution much better (tens of meters to meters) than standard radar systems. SAR processing significantly improves the resolution of point targets in both the cross-track (range) and along-track (azimuth) directions by focusing the raw radar echoes (Elachi 1988; Curlander and McDonough 1991). In order to exploit SAR data, since 1991 the SAR signal processing technique referred to as SAR interferometry (InSAR) has been developed, and it is now not far from the truth to say that InSAR revolutionized a relevant part of Earth sciences. InSAR is widely used in seismology, volcanology, hydrogeology and landslide studies (Stramondo 2008). Further frameworks wherein InSAR provides relevant contributions are the monitoring of mining regions, urban areas, strategic infrastructures (bridges, dams, nuclear power plants), etc.

In the last 20 years, InSAR was demonstrated to have unique capabilities for mapping the topography and the deformation of the Earth's surface. The InSAR technique is based on extracting the phase component of the complex SAR data (a two-dimensional record of both the amplitudes and the phases of the returns from targets) to compute the pixel-by-pixel difference of SAR signal relative to a specific area and imaged from two nearby geometric conditions. The interferogram, i.e., the result of the interferometric processing, contains the measurement of the sensor-to-target distance and of any possible change in distance. The amplitude stands for the reflectivity, while the phase is a term proportional to the sensor-to-target distance and records possible surface movements.

The interferogram corresponds to the phase difference of two images having comparable viewpoints, and it can accurately measure any shifts of the returned phase, thus computing the Earth’s surface movement towards or away from the satellite. In order to reliably measure the effects of natural disasters generating surface displacements, such as earthquakes, volcanic eruptions, landslides, or man-made activities (e.g., nuclear tests), two images acquired in two different times (one before and one after the event) are needed. This approach is the so called repeat-pass interferometry, characterized by the temporal baseline parameter which corresponds to the time separation between the two SAR scenes.

Describing the technique in deeper detail,, we can state that the interferogram is the combination of the signals S 1 and S 2 received at SAR Sensor 1 and 2, respectively [see Fig. 1 in Stramondo (2008)]. More precisely,

$$ S_{1} = A_{1} e^{{ - j\frac{{4\pi r_{1} }}{\lambda }}} \,S_{2} = A_{2} e^{{ - j\frac{{4\pi r_{2} }}{\lambda }}}, $$

where \( A_{1} \) and \( A_{2} \) are the amplitudes, and λ is the wavelength. The interferogram is the difference of the phase component and is the product of S 1 versus the complex conjugate of S 2,

$$ S_{1} \cdot S_{2}^{ * } = A_{1} A_{2} e^{{ - j\frac{{4\pi (r_{1} - r_{2} )}}{\lambda }}}. $$

Therefore, the interferometric phase \( \varphi_{\text{int}} \) can be schematically split into five terms:

$$ \varphi_{\text{int}} = \varphi_{f} + \varphi_{\text{topo}} + \varphi_{\text{displ}} + \varphi_{\text{atm}} + \varphi_{\text{err}}, $$

where \( \varphi_{f} \) is the flat Earth component (the orbital phase), the topographic phase is \( \varphi_{\text{topo}} \), the displacement phase is\( \varphi_{\text{displ}} \), the atmospheric term \( \varphi_{\text{atm}} \) and the error phase \( \varphi_{\text{err}} \). Except for this latter and the \( \varphi_{f} \), each term contains information relevant to specific issues. The \( \varphi_{\text{displ}} = \frac{4\pi }{\lambda }\Updelta R \) is the phase component accounting for the satellite-to-target distance change \( \Updelta R \).

Conversely from optical data, DInSAR can retrieve terrain modification induced by underground explosions, even if no surface changes are visible by visual analysis. Cong et al. (2007) demonstrated the suitability of DInSAR to measure the surface deformation caused by the nuclear test that took place in Nevada Test Site (NTS). Co-seismic and post-seismic displacement fields have been highlighted in this work, with a typical Gaussian shape pattern of subsidence. Similarly, Riechmann et al. (http://www.ctbto.org/specials/the-international-scientific-studies-project-iss/scientific-contributions/on-site-inspectionposters/) applied DinSAR for two case study, NST and Lop Nor (China) and, in particular for the second case study, they inferred the position of the nuclear test by observing the pattern of the interferometric coherence loss (i.e., noisy interferogram) and few surrounding fringes.

4.2 InSAR Data and Results

In order to investigate possible surface deformation caused by the registered event, three SAR images were elaborated by means of a differential SAR interferometry technique (DInSAR). In particular, the German TerraSAR-X satellite imaged the investigated area. The available data were acquired in strip-map mode, which is characterized by a swath of about 30 km cross-track and with a ground resolution of 3 × 3 m per pixel (see Table 3). The frames of the SAR images cover an area of about 30 × 60 km (see Fig. 11) which guarantees an optimal coverage of the entire investigation area.

Table 3 Characteristics of the three TerrSAR-X images used for this work
Full size table
Fig. 11
figure 11

Map of the DPRK region where the 2006 and 2009 tests were carried out. The red rectangle shows the area investigated by means of SAR interferometry technique (InSAR)

Full size image

The pre-event acquisition is dated May 18, 2009, 1 week prior to the nuclear test, while the post-event roughly span two (July 23) and three (August 14) months after the test. DInSAR has been applied to the overall possible combinations. Both the co-event pairs have very short spatial baselines. Notwithstanding the fact that the looking geometries are very close to each other, and this strongly limits the spatial decorrelation, the high and steep relief, together with the vegetation coverage, affects the interferometric coherence. In order to reduce such effects, a multilook of 13 × 15 pixels has been applied which does not have any impact, ensuring the minimum spatial resolution (24 m) able to achieve the deformation signal. The topographic phase has been removed using the SRTM (Shuttle Radar Topographic Mission) digital elevation model at 90-m spacing.

In Fig. 12, the entire TSX interferograms of the co-event (upper panel) and of the post-event pairs (bottom panel) are reported. Both interferograms are very noisy, except for some areas: on the central-south section some coherent signals are visible, mainly located inside valleys. The signals of these regions are clearly related to atmospheric artifacts. As far as the co-event interferogram (20090518–20090723) is concerned, we found a significant fringe pattern only in the northern portion of the interferogram (detail in Fig. 13). The result from InSAR shows three fringes of deformation corresponding to 40-45 mm displacement along the satellite LOS (Line of Sight), i.e., 35° from nadir. The fringe pattern is well positioned in the range/azimuth directions, reducing distortion-related phenomena and ensuring the exploitation of the full resolution of SAR sensor.

Fig. 12
figure 12

Interferograms of the whole elaborated frames. Upper panel refers to the co-event SAR pair, and bottom panel refers to the post-event data pair. Both interferograms are noisy except for some coherent areas

Full size image
Fig. 13
figure 13

Interferogram showing the detected displacement over the investigated area. A three-fringes pattern is present (red ellipse), roughly corresponding to 45-mm LOS movement. The upper right sub-panel reports the relative position of the showed interferogram and the epicenters locations (rows 5 and 6 in table 2, magenta and green points, respectively)

Full size image

The analysis of the detected signal after applying a phase unwrapping procedure can be interpreted either as a movement towards the satellite or, in other words, an inflation of the top of the mountain with respect to its southern bottom, or as a lowering of the valley to the south of the mountain with respect to the top (Fig. 14). The very short spatial baselines (64 and 55 m for co- and post-event interferograms, respectively) ensure that the TerraSAR-X pairs are poorly dependant on the topography. However, in order to verify possible residual topography, the post-event pair (20090723–20090814) has been processed too, but no fringe is present. This second interferogram is more coherent than the co-event one, because of the shorter temporal and spatial baselines. Looking closely in the same region where we found some fringes in the co-event interferogram, we detect a small coherent area, but no fringes are present. In any case, the complexity of the area in terms of topography and the dense vegetation coverage have prevented us from obtaining better results. Indeed, the capabilities of X-Band SAR, like TerraSAR-X are hampered by scene properties like those mentioned above ,and the longer temporal baseline, mainly for the co-event data, fosters the loss of interferometric coherence.

Fig. 14
figure 14

Map showing the results of the unwrapping procedure applied to the fringes of Fig. 12. This map shows a downward movement of the bottom of the valley with respect to the top of the mountain. The upper right sub-panel reports the relative position of the unwrapped interferogram and the epicenters locations (rows 5 and 6 in Table 2, magenta and green points, respectively)

Full size image

5 Discussion

The nuclear explosion database (NEDB, Bennett 2010) maintained by the Research and Development Support Services (RDSS, http://www.rdss.info/) reports, when available, ground truth information for the past nuclear tests. Events DPRK06 and DPRK09 are classified in the NEDB as GT1 (uncertainty of ±1 km on the epicentral coordinates). For these events, the ground truth location estimates are mainly based on the satellite imagery of the entrance to the tunnel that is supposed to have been used for both tests. The evidence that this tunnel was actually used for the nuclear tests is supported by considerations that can not be verified through the methodologies applied in this study (http://www.globalsecurity.org/wmd/world/dprk/kilju.htm). The identification of the tunnel entrance has been combined with analysis of the relative seismic locations of the two events, along with careful assessment of the topography in the vicinity of the tunnel complex and conditions necessary to achieve event containment (assuming normal nuclear testing practice). This analysis produced the so called “preferred” locations for the two events, which are, respectively, 41.287°N, 129.090°E for DPRK06, and 41.293°N, 129.066°E for DPRK09. They are indicated by green marks in Figs. 3 and 15. Both locations are consistent with the epicenters obtained by our seismological analysis, as reported in lines 5 and 6 of Table 2.

Our visual comparison of the pre- and post-event satellite images available from the web for the 2006 test (http://www.globalsecurity.org/wmd/world/dprk/html/kilchu-punggye-yok_comp01.htm) did not provide us with an objective recognition of the NEDB preferred locations as the real test sites. We also applied an accurate analysis on a restricted area around the preferred locations by the interferometric satellite techniques described in the previous section. As a result, we found the following: (1) a majority of the area within the frame does not show any signal coherence, which implies lack of information relevant to potential deformation patterns; (2) the areas characterized by coherence do not show evidence of deformation except for the location focused in Fig. 12.

In the previous section we showed that a few fringes clearly visible on a limited area, covering about 3 km in the east-west direction, denote a displacement of the order of 45 mm occurred between May 18 and July 23, 2009. This time period a little longer than 3 months encompasses the date of the second and larger nuclear test carried out in DPRK on 25 May 2009. The observed displacement is located near the eastern edge of a narrow mountain ridge between two nearly parallel fluvial valleys, at the confluence of two rivers (Fig. 13).

The displacements observed through the DInSAR technique are relative. Therefore, these measurements are usually based on the assumption of absolute stability of a given area with respect to which the displacements of other places are estimated. These displacements are obtained by applying the so-called unwrapping procedure to the detected fringe pattern. In our case, a stable reference zone is not available. Then we assume, as a working hypothesis, the stability of the top of the mountain, which has as a consequence the subsidence of the small red spot shown in Fig. 14 at the bottom of the valley.

As an alternative hypothesis to the preferred solutions reported by the NEDB, we may assume that the observed subsidence was induced by the May 2009 nuclear test and it occurred not far from the place where the nuclear device was detonated. There is no way to prove our assumption objectively, but it appears reasonable at least for the following considerations:

  • Subsidence is usually observed as post-seismic effect in the epicentral area of underground nuclear explosions (Houser 1969; Vincent et al., 2003);

  • It is reasonable to assume that the epicenter of the explosion is closer to the base of the mountain than to its top;

  • The fringes of the interferometric image could close more to the south, but they can not be observed due to the scarce coherence of the InSAR signals in the southern area.

Another possible cause of the measured movement can be the re-activation of a landslide located on the southern flank of the mountain (Moro et al., 2011), triggered by the nuclear test.

If our assumption is true, we must explain the circumstance that the ground truth of the May 2009 test happens to be located about 10 km to the north of the seismological location of such event, out of the 95 % confidence-level error ellipse (Fig. 15). In this respect, it must be recalled that the size of the error ellipse is just a statistical assessment of the consistence among arrival time readings, and it does not represent the absolute reliability of the location. The results of seismological location procedures strongly depend on the azimuthal variability of the travel-times due to lateral dishomogeneity in the Earth's structure, and on local time corrections for single stations. These discrepancies between theoretical and real travel-times are not taken into account in our location procedures, and could be responsible of systematic mislocations even larger than the size of the statistical error ellipse, as often noted in the seismological praxis.

Fig. 15
figure 15

Map of the DPRK region where the 2006 and 2009 tests were carried out. The two white tags marked with 5 and 6 point to the epicenters of the DPRK06 and DPRK09 events reported at lines 5 and 6 of Table 2. The respective 95 % confidence level error ellipses are shown in blue. The green tags point to the preferred locations reported by NEDB. The red square shows the area where a significant displacement has been observed through the SAR interferometry technique

Full size image

To assess the size of the effect of lateral variations of wave velocity on the epicenter locations, we performed tests on seismic events observed in Japan, an area which belongs to the same subduction region as North Korea. We compared the epicenters of two sets of densely clustered earthquakes reported by the ISC with those obtained by the Japan Meteorological Agency (JMA). In consideration of the high density of the JMA national seismological network and the good quality of its locations, these locations were taken as ground truth in our tests. The magnitude of most of the selected events ranged from 4.0 to 5.0 (similar to the m b magnitudes of DPRK06 and DPRK09), so as to assure a coverage of international seismological stations similar both in azimuth and distance distribution to that used for the nuclear tests. The results of these tests showed that systematic shifts of the epicenters computed by a network composed mostly by teleseismic and regional stations can be as large as 10 km from the ground truth (Giuntini et al., 2012).

The systematic effects due to the dishomogeneity of the Earth structure can be reduced for relative locations, as shown in Sect. 3, by applying such methods as the DDJHD. The application of a relative location method to the waveforms recorded for DPRK06 and DPRK09 showed a shift of the order of 2 km to the east-southeast of the former with respect to the latter. Assuming the ground truth location obtained through the InSAR technique, such shift would locate the DPRK06 detonation point on the same valley and nearly at the same altitude as that of DPRK09. In fact, the trend of the valley bottom in that place is west-northwest to east-southeast.

A distance of 2 km between the detonation points is not very large. However, the strong horizontal variability of the geological conditions in the area (Fig. 1) and a possibly different burial depth of the explosive device could explain the difference between the magnitude of the two tests, even if the yield of the explosions are similar. As we have seen in the introduction, a difference of 0.3 units for the m b magnitude can be easily obtained for the same yield, assuming regression laws suitable for different lithological environments.

6 Conclusions

In Sect. 2 we obtained single locations of the DPRK06 and DPRK09 events by a standard seismological method. In Sect. 3 we compared the locations of the two events by means of the DDJHD and waveform alignment methods, concluding that the two tests were carried out at a distance of 2 km from each other, the DPRK06 epicenter being shifted to the ESE with respect to the DPRK09 epicenter. However, the relative location seen in Sect. 3 doesn’t improve the absolute accuracy of any of the two locations.

The NEDB available on the web site (http://www.rdss.info/) reports for DPRK06 and DPRK09 ground truth locations based on the identification of the entrance to the tunnel where the nuclear tests are supposed to be carried out (http://www.globalsecurity.org/wmd/world/dprk/kilju.htm). This identification was supported by activities of intelligence carried out during the years preceding the 2006 test and cannot be verified by means of our purely scientific methods. These ground truth locations are consistent with the epicentral locations obtained by standard seismological methods, within their error ellipses.

The analysis of SAR satellite images does not show any displacement on or close to the NEDB ground truth locations, whereas it puts in evidence an area of subsidence centered about 10 km to the north of the seismological locations. The observed displacement took place during 3 months encompassing the date of the 2009 test. We considered the hypothesis that this subsidence area is strictly related to the nuclear test, possibly as a landslide triggered at a certain distance from the detonation point. This hypothesis can not be rejected just by the separation of 10 km between the epicenter of the event and the subsidence area, by taking into account the systematic shift of the seismological location due to non-homogeneity of wave velocity affecting the travel times of seismic waves at a global scale.

The present study, together with other recent studies about the use of satellite imageries in underground explosion monitoring, supports the idea of the improvement of the verification regime by additional monitoring technologies such as satellite monitoring, in accordance with Art. IV, Section A, Paragraph 11 of the CTBT.