Open Access
Add to BibSonomy
Abstract
Widely distributed optical fibers, together with phase-sensitive laser interferometry, can expand seismic detection methods and have great potential for epicenter localization. In this paper, we propose an integral response method based on a forward transmission scheme. It uses spectrum analysis and parameter fitting to localize the epicenter. With the given shape of the fiber ring, the integral phase changes of light propagating in the forward and reverse directions can be used to determine the direction, depth, distance of the epicenter, and seismic wave speed. For the noisy case with SNR = 20 dB, the simulation results show ultrahigh precision when epicenter distance is 200 km: the error of the orientation angle is ∼0.003°±0.002°, the error of the P-wave speed is ∼0.9 ± 1.2 m/s, the error of the epicenter depth is ∼9.5 ± 400 m, and the error of the epicenter distance is ∼200 ± 760 m.
© 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
Forward transmission [1] and backscattering [2–11] are two main schemes of fiber vibration sensing technology. The backscattering scheme can use the discrete Direction of Arrival (DOA) technique to find the direction of incoming waves [8,12]. However, it has disadvantages of large optical loss and limited detection bandwidth, which makes it just suitable for dark fibers [4,5]. Furthermore, the low-loss fiber is a never-ending goal of mankind [13–15]. Virtually every Optical Fiber Communication conference features reports of hard-fought battles to eke out even a fraction of a dB/km lower loss each year, because this significantly impacts system reach and cost over long-haul distances [16]. However, these hard-fought battles are not good news for backscattering laser interferometry, which relies on Rayleigh scattering. The forward transmission scheme overcomes the above problems of backscattering. However, it has drawbacks in vibration localization. Since the forward transmission scheme detects the phase change of whole fiber, how to analyze vibrations occurring in different locations is the main problem.
Recently, the preliminary wave method was proposed to localize epicenters [1]. As shown in Fig. 1(a), on each fiber link, the phase difference of two laser interferometers can be obtained, and then points A1 and A2 where the seismic wave first arrives at link 1 and link 2 can be determined. The normal line of points A1 and A2 will point to the epicenter. When two fiber links are implemented, the intersection of two normal lines is the epicenter.
Fig. 1. Schematic diagrams of the preliminary wave method and its defects. (a) Schematic diagram of the preliminary wave method and the detection blind area of two fiber links. (b) The schematic diagram when the seismic wave first arrives. (c) The schematic diagram when the derived preliminary point A′ deviates from the real point A. (d) Schematic diagram of seismic envelope deformation.
There are several potential problems with this method. First, since seismic vibration is a nonstationary process, the intensity envelope of seismic waves gradually becomes stronger. This means that the strength of the preliminary wave is weak, and the affected fiber length is limited (Fig. 1(b)). This will lead to a low signal-to-noise ratio (SNR), which dramatically affects the localization result [17]. Second, with the propagation of the seismic wave (Fig. 1(c)), a larger section of fiber will be disturbed by the earthquake. The calculated point A′ determined by the preliminary wave method will deviate from the real point A, and the detected direction of the epicenter will deviate from the real direction. Third, for a certain fiber layout, vibration localization will have a large blind area using the preliminary wave method. For example, when the end point of a fiber link is the nearest point of the link to the epicenter, this method will fail to localize it. Figure 1(a) shows the blind area of the fiber link layout in Ref. [1]. The yellow part is the locatable area for both links. The green and blue parts are the locatable areas for just one link. The gray part indicates the detection blind area.
Furthermore, the seismic envelope consists of several seismic waves caused by different propagation paths, and the envelope may be formed by the wavefront in different directions. The prerequisite of the preliminary wave method is that the propagation speeds of the seismic waves are the same along different directions. Actually, this is usually the case when seismic waves propagate in the upper crust (epicenter distance less than 200 km and epicenter depth less than 15 km). In this case, the preliminary wave propagates directly from the epicenter to the ground, referred to as direct wave. With the increase of epicentral distance and depth, other seismic waves may overtake direct wave. In this case, the fastest wave is head wave, which travels through the mantle with faster speed than that in upper crust. Therefore, the seismic envelope is complex due to different crust thickness and seismic wave velocity in the mantle on a large scale (Fig. 1(d)). This will lead to a distorted seismic envelope and finally cause epicenter positioning errors. Additionally, the preliminary wave method cannot obtain the depth, wave speed and distance of the earthquake event, which limits its broad application.
In this paper, we propose an integral response method based on forward-transmission laser interferometry. The orientation, depth, distance of the epicenter, and seismic wave speed can be determined by integrating the phase changes throughout the fiber link. The simulation result shows that the proposed method can overcome the above 4 problems of the preliminary wave method. For the noisy case with SNR = 20 dB, when the epicenter distance is 200 km, the error of the orientation angle is ∼0.003°±0.002°, the error of the P-wave speed is ∼0.9 ± 1.2 m/s, the error of the epicenter depth is ∼9.5 ± 400 m, and the error of the epicenter distance is ∼200 ± 760 m.
2. Schematic and model
Figure 2 shows the schematic diagram of epicenter relative to the fiber optic link. TX1, TX2 and RX1, RX2 are the transmitting and receiving ports of two laser interferometers, respectively [18–20]. The phase change of the laser interferometer can be extracted from the interference of the sensing signal with reference signal using a digital IQ demodulator. Considering that the depth of epicenter E is H, the projection of epicenter E on the ground is point O, as shown in Fig. 2 (a). To analyze the earthquake effect on the ground, a cylindrical coordinate system is established with O as the origin and the vertical direction as the cylindrical axis z. Any horizontal direction can be set as the polar axis x. Considering the microelement dl on the fiber link at point , its position can be expressed using polar diameter r, angular position $\theta$ and height z, where height z = 0. The direction vectors are the radial direction vector $\widehat r$, the circumferential direction vector $\widehat \theta $, and the vertical direction vector $\widehat z$. The fiber length from the link’s end point to microelement dl is l (clockwise direction). The angle between the tangent direction and the radial direction $\widehat r$ on the microelement is β. The incident angle between the wave vector $\overrightarrow k $ and the vertical direction $\widehat z$ is ${i_p}$. The distance from the centroid C of the fiber link to the epicenter E is D. Figure 2(b) shows the top view of Fig. 2(a). The radial direction vector at point C is $\mathrm{\hat{r}^{\prime}}$. To describe the shape of the fiber optic link and its location relative to the epicenter E, we establish a polar coordinate system with C as the pole, and the direction from the observation station to the pole is the reference direction of the polar axis $x^{\prime}$. The polar angle of O in the polar coordinate system is $\alpha $.
Fig. 2. (a) Schematic diagram of the epicenter relative to the fiber optic link. (b) Top view of the schematic diagram.
To calculate the phase change of the laser interferometer, a model should be built to describe the earthquake-induced strain throughout the fiber link. Seismic waves can be classified into P-wave, SH-wave, and SV-wave. Normally, P-wave comes first. Taking the P-wave as an example, the displacement induced by the P-wave is only in the $\mathrm{\hat{r}\ -\ }\hat{z}$ plane. First, we analyze the radial displacement of centroid C in the frequency domain, which can be expressed as ${u_{C - r}}(\omega )$, where $\omega $ is the frequency of the seismic wave. Then, we analyze the earthquake effect on the micro element dl. The radial displacement on the micro element dl can be expressed as ${u_{dl - r}}(\omega )$ [21]:
The fiber in the horizontal plane is insensitive to the shear strain and the linear strain perpendicular to the tangential direction of the fiber [6,23]. Therefore, we can deduce the strain parallel to the tangential direction of the fiber ${\varepsilon _{dl}}(\omega )$ from applying Eqs. (2)–(4) to strain transformation model [24]:
3. Method and results
In the proposed integral response scheme, the phase change of two laser interferometers can be calculated theoretically as follows:
To normalize the amplitudes, we divide ${\varphi ^ + }(\omega )$ by ${\varphi ^ - }(\omega )$ and obtain the Link Response Function (LRF):
In Eq. (9), the terms $\xi$, D and ${u_{C - r}}(\omega )$ have been canceled because they have no relation with the variable dl. LRF is a function of the seismic frequency $\omega$, the epicentral direction $\alpha $, the P-wave speed vp, the distance from the centroid to epicenter D, and the depth of epicenter H. Actually, the terms H and $\alpha $ are hidden in variables $\eta$, $\beta$ and d.
Using the forward-transmission laser interferometry, the phase change ${\varphi ^ + }(t)$ and ${\varphi ^ - }(t)$ can be detected. Through Fourier transform, we can obtain their spectra ${\varphi ^ + }(\omega )$ and ${\varphi ^ - }(\omega )$, respectively. $DRF(\omega ) = {\varphi ^ + }(\omega )/{\varphi ^ - }(\omega )$ is the Detected Response Function (DRF). According to the minimum mean square error principle (MMSE), we define the loss function $Loss({\alpha ,{v_p},D,H} )$ as:
Here $angle[LRF(\alpha ,{v_p},D,H,\omega )]$ and $angle[DRF(\omega )]$ mean the phase spectrum part of LRF and DRF, respectively. We fit $DRF(\omega )$ with $LRF(\alpha ,{v_p},D,H,\omega )$ by the genetic algorithm and so on. In other words, we change the value of $\alpha ,{v_p},D,H$ to minimize the deviation of the LRF from the DRF. Theoretically, the LRF fits the DRF best when $\alpha ,{v_p},D,H$ is equal to their true value. We can use neural networks to increase the confidence and noise resistance. The detailed methodology of the neural network is shown in Supplement 1. As a summary, we use a flow chart to clearly show steps of the proposed method in Fig. 3.
To demonstrate the novelty of the proposed integral response method, we carry out a simulation using a circle fiber ring with a radius of 50 km. To improve the credibility of our simulation, we use a real earthquake waveform (event ID: 39338407) from Southern California Earthquake Data Center (SCEDC), to mimic the actual seismic wave [26]. The earthquake event is detected 200 km away from the epicenter, thus the simulation epicentral distance D is set as 200 km. Other parameters are set as follows: the epicenter depth H is 20 km, the P-wave speed is 5500 m/s, and the orientation of epicenter $\alpha$ is ${70^ \circ }$. We also simulate a 1/f noise to mimic the real noise of detection system [27]. From the detected seismic events of Ref. [1], we can roughly obtain the power ratio between the detected seismic signals and system noise. Thus, we define SNR as the ratio between the power of phase change signal and the power of added 1/f phase noise. Their SNRs are within the range of 10 ∼ 30 dB, approximately. Therefore, we carry out simulation with SNR = 10 dB, 20 dB and 30 dB to verify the method. The simulation results are shown in Fig. 4 and Table 1. Bias is the deviation from the true value. The halfwidth is the full width at half maximum of the inverse of Loss function in a single measurement. First, we analyze the result of epicenter orientation $\alpha$ and P-wave speed ${v_p}$ under different noise levels. In Fig. 4(a), plots I, II, and III are the localization cases when SNR = 10 dB, 20 dB and 30 dB, respectively. They are plotted as the reverse of the loss function. A larger value means a better fitting of the parameters at that point. The polar diameter represents the seismic wave speed ${v_p}$, and the polar angle is the epicenter orientation $\alpha$. The epicenter can be clearly extracted from the polar plot. The errors are given in Table 1. Figures 4(b), (c), (d), and (e) show the simulation errors of the epicenter depth H, the epicenter distance D, the epicenter orientation $\alpha $ and the seismic wave speed ${v_p}$ respectively, under different epicenter orientations and P-wave speeds when SNR = 20 dB. Because the circle fiber ring used for simulation is symmetric, the error of ${v_p}$, H and D is larger when the angle between the symmetry axis is within 0.5°, and this part is not included in the figure. This can be avoided by choosing an asymmetric fiber link. When the SNR is 20 dB, the error of the orientation angle is ∼0.003°±0.002°, the error of the P-wave speed is ∼0.9 ± 1.2 m/s, the error of the epicenter depth is ∼9.5 ± 400 m, and the error of the epicenter distance is ∼200 ± 760 m.
Fig. 4. Seismic parameter estimation results of the circle fiber ring under different conditions. (a) Errors of epicenter direction and the wave speed with different SNRs. (b) Errors of epicentral distance H with different epicentral directions and wave speeds. (c) Errors of the epicenter depth D with different epicentral directions and wave speeds. (d) Errors of the epicenter orientation $\alpha $ with different epicentral directions and wave speeds. (e) Errors of the seismic wave speed ${v_p}$ with different epicentral directions and wave speeds. The polar diameter represents the seismic wave speed ${v_p}$, and the polar angle is the epicenter orientation $\alpha $.
Table 1. The localization errors of $\alpha ,{v_p},D,H$ with different SNRs
Actually, the shape of fiber ring doesn’t affect the parameter fitting results. As illustrated in Eqs. (7) and (8), the only part related to fiber ring is to integrate the phase change along the fiber. As long as the fiber ring can be described using a coordinate system, it can be taken for arbitrary shape, like the normally existing optical fiber cable. We also carry out a simulation using a square-shaped fiber ring with a side length of 50 km, to demonstrate the practicability of the proposed integral response method. The earthquake parameters are the same as above (plot II of Fig. 4(b)). The simulation results are shown in Fig. 5(a) and Table 1.
Fig. 5. Seismic parameter estimation results under different conditions. (a) Errors of epicenter direction and the wave speed with square-shaped fiber ring. (b) Errors of epicenter direction and the wave speed when affected by both direct waves and head waves. (c) Errors of epicenter direction and the wave speed when affected by head waves only. The polar diameter represents the seismic wave speed ${v_p}$, and the polar angle is the epicenter orientation $\alpha $.
In conclusion, the integral response method can effectively overcome the disadvantages of the preliminary localization method. First, the proposed method is not affected by the low SNR problem of Fig. 1(b). Because it is not limited to detect the preliminary seismic wave, but can select a segment with relatively stronger seismic wave. Second, this method will take the phase changes caused by the seismic wave on the whole fiber into account, rather than only using the detected signal on a certain point. Thus, the deviation between the calculated preliminary point A′ from the actual point A can be avoided. This is demonstrated by the case of the square-shaped fiber ring in Fig. 5(a). Third, the preliminary method is influenced by the envelope distortion caused by the superposition of direct waves and head waves. On the contrary, the proposed integral response method is insensitive to head waves.
To demonstrate that, we carry out simulation under two situations: fiber affected by both direct waves and head waves simultaneously, and fiber only affected by head waves. Considering that the attenuation of head waves is usually much greater than that of direct waves, we set the amplitude ratio between head waves and direct waves to be 1:7. Other earthquake parameters are the same as above and the circle fiber ring is used. The simulation results are shown in Figs. 5(b) and (c). When the fiber is affected by both waves, the epicenter still can be correctly localized, as shown in Fig. 5(b). This is because that the proposed method is insensitive to the head waves. The simulation result of a noisier case can be found in Supplement 1. Figure 5(c) shows the localization result when the fiber is only affected by the head waves, we won’t get a valid localization peak.
Similar conclusions can be made for the SV-wave and SH-wave. The only difference between P-wave and S-wave is: the vibration directions of SV-waves and SH-waves are perpendicular to the direction of wave propagation [21]. For SV-wave, the vibration direction is in the $\mathrm{\hat{r}\ -\ }\hat{z}$ plane of Fig. 2(a). We just need to change the parameter $\eta$ and replace the wave number of P-wave with the wave number of SV-wave ${\textrm{k}_{SV}}$, to adjust the model from P-wave case to SV-wave case. LRFSV can be calculated as below:
Besides, we can use different LRF functions to fit the DRF in different seismic cases. For example, when P-wave is the main component of seismic waves, the correct localization result can be obtained by using the P-wave LRF to fit DRF, and no effective result can be obtained by using LRFSV to fit DRF. When detecting shallow focus earthquakes, the incidence angle is large. ${\eta _P}$ and ${\eta _{SV}}$ are much smaller than ${\eta _{SH}}$ so the effects of P-wave and SV-wave are negligible compared to SH-wave. The phase changes are mainly caused by SH-wave. In this case, we can fit the DRF with LRFSH without being affected by the P-wave and SV-wave.
4. Conclusion
In this paper, we analyze the problems of current epicenter localization method. To solve the problems of low SNR, blind zone and localization error of the current scheme, we propose an integral response method. With two counterpropagating laser interferometers, the integral phase changes can be analyzed. From their phase spectra and parameter fitting, the optimal parameters of the epicenter orientation $\alpha $, the seismic wave speed, the depth of epicenter H and the distance from centroid to epicenter D can be obtained without relying on a dense seismic network. This method can enrich the earthquake monitoring method and broaden the application fields of existing fiber network.
Funding
National Natural Science Foundation of China (62171249, 61971259); Tsinghua Initiative Scientific Research Program.
Disclosures
The authors declare no conflicts of interest.
Data availability
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.
Supplemental document
See Supplement 1 for supporting content.
References
1. G. Marra, C. Clivati, R. Luckett, A. Tampellini, J. Kronjager, L. Wright, A. Mura, F. Levi, S. Robinson, A. Xuereb, B. Baptie, and D. Calonico, “Ultrastable laser interferometry for earthquake detection with terrestrial and submarine cables,” Science 361(6401), 486–490 (2018). [CrossRef]
2. A. Sladen, D. Rivet, J. P. Ampuero, L. De Barros, Y. Hello, G. Calbris, and P. Lamare, “Distributed sensing of earthquakes and ocean-solid Earth interactions on seafloor telecom cables,” Nat. Commun. 10(1), 5777 (2019). [CrossRef]
3. N. J. Lindsey, T. C. Dawe, and J. B. Ajo-Franklin, “Illuminating seafloor faults and ocean dynamics with dark fiber distributed acoustic sensing,” Science 366(6469), 1103–1107 (2019). [CrossRef]
4. E. F. Williams, M. R. Fernandez-Ruiz, R. Magalhaes, R. Vanthillo, Z. Zhan, M. Gonzalez-Herraez, and H. F. Martins, “Distributed sensing of microseisms and teleseisms with submarine dark fibers,” Nat. Commun. 10(1), 5778 (2019). [CrossRef]
5. J. B. Ajo-Franklin, S. Dou, N. J. Lindsey, I. Monga, C. Tracy, M. Robertson, V. Rodriguez Tribaldos, C. Ulrich, B. Freifeld, T. Daley, and X. Li, “Distributed Acoustic Sensing Using Dark Fiber for Near-Surface Characterization and Broadband Seismic Event Detection,” Sci. Rep. 9(1), 1328 (2019). [CrossRef]
6. M. Karrenbach, S. Cole, L. LaFlame, E. Bozdağ, W. Trainor-Guitton, and B. Luo, “Horizontally orthogonal distributed acoustic sensing array for earthquake- and ambient-noise-based multichannel analysis of surface waves,” Geophys J Int. 222(3), 2147–2161 (2020). [CrossRef]
7. F. Walter, D. Graff, F. Lindner, P. Paitz, M. Kopfli, M. Chmiel, and A. Fichtner, “Distributed acoustic sensing of microseismic sources and wave propagation in glaciated terrain,” Nat. Commun. 11(1), 2436 (2020). [CrossRef]
8. N. J. Lindsey, E. R. Martin, D. S. Dreger, B. Freifeld, S. Cole, S. R. James, B. L. Biondi, and J. B. Ajo-Franklin, “Fiber-Optic Network Observations of Earthquake Wavefields,” Geophys. Res. Lett. 44(23), 11–792 (2017). [CrossRef]
9. A. Mateeva, J. Lopez, H. Potters, J. Mestayer, B. Cox, D. Kiyashchenko, P. Wills, S. Grandi, K. Hornman, B. Kuvshinov, W. Berlang, Z. Yang, and R. Detomo, “Distributed acoustic sensing for reservoir monitoring with vertical seismic profiling,” Geophysical Prospecting 62(4), 679–692 (2014). [CrossRef]
10. C. Yu, Z. Zhan, N. J. Lindsey, J. B. Ajo-Franklin, and M. Robertson, “The Potential of DAS in Teleseismic Studies: Insights From the Goldstone Experiment,” Geophys. Res. Lett. 46(3), 1320–1328 (2019). [CrossRef]
11. P. Jousset, T. Reinsch, T. Ryberg, H. Blanck, A. Clarke, R. Aghayev, G. P. Hersir, J. Henninges, M. Weber, and C. M. Krawczyk, “Dynamic strain determination using fibre-optic cables allows imaging of seismological and structural features,” Nat Commun 9(1), 2509 (2018). [CrossRef]
12. S. Rost, “ARRAY SEISMOLOGY: METHODS AND APPLICATIONS,” Rev. Geophys. 40(3), 2-1–2-27 (2002). [CrossRef]
13. X. Zhao, X. Wu, X. Lan, J. Luo, L. Zhang, P. Li, J. Xiang, W. Ma, and S. Wang, “5-tube hollow-core anti-resonant fiber with ultralow loss and single mode,” Opt. Commun. 501, 127347 (2021). [CrossRef]
14. Q. Liu, Z. Ren, W. Liu, Y. Sun, T. Lv, C. Liu, W. Lu, B. Li, Y. Jiang, T. Sun, and P. K. Chu, “A photonic quasi-crystal fiber composed of circular air holes with high birefringence and low confinement loss,” Optik (Munich, Ger.) 231, 166497 (2021). [CrossRef]
15. T. Akand, M. J. Islam, and M. R. Kaysir, “Low loss hollow-core optical fibers conjoining tube lattice and revolver structures,” Results in Optics 1, 100008 (2020). [CrossRef]
16. J. Ballato, “Optical Fiber: Through the Looking Glass,” Opt. Photonics News 33(3), 30–40 (2022). [CrossRef]
17. R. Yu, M. Yuan, and Y. Yu, “Developed empirical model for simulation of time-varying frequency in earthquake ground motion,” Earthquakes and Structures 8(6), 1463–1480 (2015). [CrossRef]
18. G. Wang, H. Si, Z. Pang, B. Zhang, H. Hao, and B. Wang, “Noise analysis of the fiber-based vibration detection system,” Opt. Express 29(4), 5588–5597 (2021). [CrossRef]
19. G. Wang, Z. Pang, B. Zhang, F. Wang, Y. Chen, H. Dai, B. Wang, and L. Wang, “Time shifting deviation method enhanced laser interferometry: ultrahigh precision localizing of traffic vibration using a urban fiber link,” Photonics Res. 10(2), 433–443 (2022). [CrossRef]
20. Z. Pang, G. Wang, B. Wang, and L. Wang, “Comparison Between Time Shifting Deviation and Cross-correlation Methods,” J. Lightwave Technol. 40(9), 3003–3009 (2022). [CrossRef]
21. S. Seth and W. Michael, An Introduction to Seismology, Earthquakes, and Earth Structure, 1st ed. (Blackwell, 2003), pp. 29–116.
22. M. Sadd, Elasticity: Theory, Applications, and Numerics, 4th ed. (Elsevier, 2021), pp. 31–55.
23. X. Wu, M. E. Willis, W. Palacios, A. Ellmauthaler, O. Barrios, S. Shaw, and D. Quinn, “Compressional- and shear-wave studies of distributed acoustic sensing acquired vertical seismic profile data,” The Leading Edge 36(12), 987–993 (2017). [CrossRef]
24. M. Sadd, Elasticity: Theory, Applications, and Numerics, 4th ed. (Elsevier, 2021), pp. 40–41.
25. B. N. Kuvshinov, “Interaction of helically wound fibre-optic cables with plane seismic waves,” Geophysical Prospecting 64(3), 671–688 (2016). [CrossRef]
26. Southern California Earthquake Data Center (SCEDC), Caltech. Dataset, ID=39338407, Caltech/USGS Southern California Seismic Network (SCSN) (2013), https://dx.doi.org/10.7909/C3WD3xH1.
27. E. Ip, Y. K. Huang, G. Wellbrock, T. J. Xia, M. F. Huang, T. Wang, and Y. Aono, “Vibration detection and localization using modified digital coherent telecom transponders,” J. Lightwave Technol. 40(5), 1472–1482 (2022). [CrossRef]
