• Moho depth distribution in Gangwon Province consistent with Airy’s isostacy

• Shallow low-velocity layers observed beneath sedimentary basins

• Velocity inversion within 4-12 km depths with mid-crustal discontinuities

2016년 경주지진, 2017년 포항지진, 2022년 괴산지진 등 한반도 내에서 모멘트 규모(Mw) 4.0 이상의 지진이 최근 잇따라 발생하면서 한반도 내 지체구조에 대한 관심이 증가하고 있다. 한반도 남부에는 최소 450개 이상의 단층이 분포하고 있으며(Bae, 2016) 특히 강원도 지역은 북쪽으로 인제단층, 금왕단층이 지나고 있고, 삼척시를 중심으로 영월군과 평창군을 따라 단층선의 밀도가 매우 높은 복잡한 지체구조를 형성하고 있기에(Fig. 1) 잠재적인 지진 발생 위험이 늘 상존하고 있다(Park, 2007; Lee et al., 2008; Kim and Lee, 2017).

Figure 1. Locations of broadband seismic stations (blue and red triangles) and faults (black thin lines) in and around Gangwon Province. Two red thick lines indicate the locations of (a) Inje fault and (b) Geumwang fault, respectively.

일반적으로 강원도는 영남권에 비해 지진 발생 빈도수는 적으나 2007년에 규모 4.8의 오대산 지진이 발생하였고 2019년에는 동해 바다에서 규모 4.3의 지진이 발생하는 등 중규모의 지진이 꾸준히 발생하고 있으며, 조선왕조실록의 기록에 의하면 1681년 양양에서 추정 규모 6.5의 대지진이 발생한 적이 있다고 보고되었다(Lee et al., 2004). 따라서 강원도에서 앞으로 발생할 수 있는 중규모 이상의 지진의 정밀한 분석을 위해 이전보다 더 정밀한 지각속도구조에 대한 연구가 필요하다.

지각속도구조는 지진이 발생한 진원의 위치와 지진의 규모를 정확히 알아내는데 필수적인 요소(Kim et al., 2006)로서 한반도와 같이 지진 활동이 비교적 적은 지역은 인공파원을 이용한 탄성파 굴절법 및 반사법 탐사를 적용하여 지각속도구조에 대한 2차원 영상화가 활발하게 진행되었으며(Cho et al., 2006; Cho et al., 2013; Kim et al., 2007), 지진 토모그래피를 이용한 3차원 영상화도 진행된 바 있다(Park et al., 2002). 하지만 굴절법 및 반사법 탐사의 경우 많은 시간적, 금전적 비용이 소모되고 지진파 토모그래피는 수직적인 해상도가 낮다는 문제점이 존재한다. 이에 비해 원거리 지진파를 이용하는 수신함수는 별도의 파원을 필요로 하지 않으면서 속도변화에 대한 해상도가 높아 관측소 하부에 대한 국소적인 지각속도구조를 저비용 고효율로 분석하는데 널리 사용되어왔다(Kim and Lee, 2001; Chang and Baag, 2005; Yoo et al., 2007; Jeon et al., 2013).

이전의 몇몇 연구들은 한반도 남부에 설치된 관측소 하부에서의 모호면 깊이를 계산하고 이를 보간한 모호면 등고선을 제시하였으나, 당시 강원도 지역의 경우 한국지질자원연구원(KIGAM: Korea Institute of Geoscience and Mineral Resources)과 기상청(KMA: Korean Meteorological Administration)이 설치한 광대역 관측소가 10개소도 채 되지 않아 지역적 해상도가 떨어지거나(Chang and Baag, 2005), 단주기 관측소의 주파수 한계로 인해 일부 부정확한 결과가 나오는 경우도 발생하였다(Yoo et al., 2007; Park et al., 2009). 하지만 2010년대 들어 기상청에서 강원도 지역에 다수의 광대역 관측소를 설치하면서 보다 높은 해상도의 강원도 지역 지각속도구조를 도출할 수 있게 되었다.

수신함수를 사용하여 지각속도구조를 분석하는 방법은 원거리 지진파의 상대적인 주시를 이용하기 때문에 모호면의 깊이와 지각속도 사이의 상충효과(trade-off)로 인해 속도 모델의 비유일성(non-uniqueness) 문제가 발생한다(Ammon et al., 1990). Chang et al. (2004)은 이 문제를 보완하기 위하여 전단파 속도 대비에 민감한 수신함수와 전단파 평균 속도에 민감한 표면파 분산을 함께 역산하였다. 또한, 선험적인 속도구조 정보가 거의 없는 지역에서 지역 최저값(local minima)에 빠질 우려가 있는 선형화 역산을 대신하여 전역 최적화(global optimization) 기법인 유전 알고리즘을 이용하였고, 수치 테스트를 통해 좋은 해상도를 얻을 수 있음을 보여주었다. 또한, 이 방법에서는 V_p/V_s 비(ratio)를 이용하여 수신함수와 표면파분산에 더 많은 영향을 주는 S파 속도구조를 P파 속도구조로 변환한다. P파 속도구조가 진원의 위치 측정 등 기본적인 지진 분석에 더 많이 쓰이고 있기 때문이다.

따라서 본 연구에서는 강원도 지역과 그 주변에 설치된 21개의 관측소(Fig. 1; Table 1)에 대하여 먼저 Zhu and Kanamori (2000)가 제안한 H-κ 중합법을 통해 1차적으로 모호면의 깊이와 V_p/V_s 비를 추정하고, 이를 Chang et al. (2004)이 제시한 수신함수와 표면파 분산의 유전 알고리즘 연합 역산에 적용하여 각 관측소 하부에 대한 1차원 P파 속도구조를 도출하였다. 역산 결과를 통해 강원도 지역에 설치된 관측소의 천부에 나타나는 저속도층과 중간지각에 나타나는 속도역전층을 확인할 수 있었다.

Table 1 . Locations of broadband seismic stations in and around Gangwon Province.

Station name	Network	Location
		Latitude (°N)	Longitude (°E)	Elevation (m)
BURB	KMA	37.2317	127.7486	88.0
CGUB	KMA	37.5583	127.7136	160.0
CHC2	KMA	37.7776	127.8145	269.0
CHNB	KIGAM	38.2710	127.1211	184.9
CHYB	KMA	36.9439	128.9144	366.0
GAPB	KMA	37.8442	127.4949	123.8
GEJB	KMA	38.4165	128.3858	135.0
GMHB	KMA	38.2500	127.4209	251.0
IMWB	KMA	37.2380	129.3419	55.0
JECB	KMA	37.1595	128.1945	317.0
JSB	KIGAM	37.3175	128.6854	322.9
KOSB	KMA	38.6028	128.3596	62.1
KSA	KIGAM	38.5953	128.3516	105.2
NAMB	KMA	37.7711	128.3863	656.0
PYCB	KMA	37.3775	128.3947	344.0
SH2B	KMA	38.2686	128.2525	423.0
SND	KIGAM	37.1579	128.7962	779.1
SNGB	KMA	37.3465	129.0862	124.0
YAYB	KMA	38.0255	128.7160	45.0
YKB	KIGAM	38.2254	127.9903	303.2
YOJB	KMA	36.8720	128.5166	252.0

2.1. 수신함수(Receiver Function)

진앙거리 30° 이상에서 발생하는 원거리 지진파가 관측소에 거의 수직으로 입사하면서 지진계에 진원의 시계열 함수와 관측소 하부의 구조, 기기응답함수에 의한 정보들을 포함하는 파형으로 기록된다(Cassidy, 1992). 각 관측소에 대한 지진파의 후방위각(backazimuth)을 통해 수직성분(vertical component), 법선성분(radial component), 접선성분(transverse component)의 3성분으로 변환시키고 나면 시간영역에서의 지진 파형은 다음과 같이 나타낼 수 있다(Langston, 1979).

(1)$$\begin{array}{'}{D}_V\left(t\right)=I\left(t\right)*S\left(t\right)*E_V\left(t\right)\\ D_R\left(t\right)=I\left(t\right)*S\left(t\right)*E_R\left(t\right)\\ D_T\left(t\right)=I\left(t\right)*S\left(t\right)*E_T\left(t\right)\end{array}$$

D_V,R,T는 수직성분, 법선성분, 접선성분의 지진파형이며 I(t)는 기기응답함수, S(t)는 진원시계열함수, E_V,R,T는 수직성분, 법선성분, 접선성분에서의 구조응답함수 즉, 수신함수를 나타낸다. 지진계에 거의 수직으로 입사하는 원거리 P파가 모호면과 같은 속도 대비층을 만나게 되면 파의 변환이 일어나게 되고, 이를 Ps변환파라고 일컫는다. 이때 변환되지 않은 직접 P파만이 수직성분에 강하게 반응하고, S파로 변환된 위상들은 수직성분의 구조응답함수 E_V(t)에는 거의 나타나지 않게 된다. 이는 E_V(t)가 하나의 디락델타함수로 근사될 수 있음을 의미하며 이에 따라 수직성분의 파형은 다음과 같이 나타낼 수 있다(Langston, 1979).

(2)$$D_V\left(t\right)=I\left(t\right)*S\left(t\right)*\delta \left(t\right)$$

(3)$$D_V\left(t\right)\simeq I\left(t\right)*S\left(t\right)$$

이를 이용하여 나머지 두 성분에 대한 구조응답함수를 다음과 같이 나타낼 수 있다.

(4)$$E_R\left(t\right)=\frac{D_R\left(t\right)}{I\left(t\right)*S\left(t\right)}-\frac{D_R\left(t\right)}{D_V\left(t\right)}$$

(5)$$E_T\left(t\right)=\frac{D_T\left(t\right)}{I\left(t\right)*S\left(t\right)}-\frac{D_T\left(t\right)}{D_V\left(t\right)}$$

즉, 수신함수는 법선성분 혹은 접선성분을 수직성분으로 곱풀기(deconvolution)하여 형성된다.

2.2. 표면파 분산(Surface Wave Dispersion)

수신함수만을 사용하여 속도구조를 추정할 경우 수신 함수에는 직접 P파와 Ps변환파의 상대적인 주시정보가 포함되어 있기 때문에 속도와 모호면 깊이 사이에 상충효과가 발생한다(Ammon et al., 1990). 즉, 모호면이 깊고 지각의 지진파 속도가 빠른 경우와 모호면이 얕고 지각의 지진파 속도가 느린 경우를 분간할 수 없게 된다. 따라서, 지각의 평균 지진파 속도에 민감한 표면파 분산곡선 자료를 추가로 사용하면 이러한 수신함수의 비유일성(non-uniqueness) 문제를 해결할 수 있다. 수신함수와 표면파 분산곡선 자료를 함께 사용하여 역산하면 비슷한 수신함수를 가지는 여러 속도 모델 중 주어진 평균 S파 속도와 일치하는 하나의 모델을 선택함으로써 비유일성 문제를 보완할 수 있다.(Julià et al., 2000).

2.3. 유전 알고리즘(Genetic Algorithm)

유전 알고리즘은 자연계의 진화 핵심 원리인 자연 선택과 유전자의 개념에 기초한 전역 최적화 기법으로, 다양한 초기 모델을 생성하고 선택(selection), 교차(crossover), 변이(mutation) 과정을 통해 최적화된 해를 찾는다(Holland, 1975; Stoffa and Sen, 1991). 이 방법은 초기 모델에 대한 의존성이 낮기 때문에 지각속도구조에 대한 사전 정보가 부족한 한반도에서 선형 역산보다 더 유용하게 사용될 수 있다(Chang et al., 2004).

본 연구에서는 적은 인구(population) 설정으로 기존 유전알고리즘에 비하여 빠르게 수렴하는 마이크로 유전 알고리즘을 채택하였다(Krishnakumar, 1990; Chang et al., 2004). 이 방법은 더 이상 수렴하지 않을 경우 최적의 모델을 제외하고 다른 모델들을 무작위로 생성된 모델로 교체함으로써 기존의 유전 알고리즘의 수많은 인구와 변이 단계를 대체한다. 또한, 각 유전자(gene)마다 교차 확률을 동일하게 설정하는 균일 교차(uniform crossover)를 사용할 경우 좋은 결과를 보여준다고 보고되었다(Carroll, 1996). 본 연구의 유전 알고리즘으로 생성되는 속도 모델은 지표부터 6 km 깊이까지 1 km 두께, 6 km부터 40 km 깊이까지 2 km 두께로 총 23개의 층으로 구성되어 있다.

수신함수 역산을 진행할 경우 종종 발생하는 깊이에 따른 급격한 속도 변화를 억제하기 위하여 거칠기 기준(roughness norm)을 도입하였다(Ammon et al., 1990).

(6)$$roughness=\sum _{i=1}^{n-2}\frac{|{\alpha }_i-2{\alpha }_{i+1}+{\alpha }_{i+2}|}{\left(n-2\right)}$$

α_i는 i번째 층의 P파 속도이며, n은 속도 모델 층의 총 개수이다. 이때 거칠기 기준이 커질수록 파형 일치도 (waveform fitness)가 증가하고 거칠기 기준이 작아질수록 파형 일치도가 감소하는 경향이 있기 때문에 두 개의 매개변수 사이에 상충효과가 발생한다. 따라서, 본 연구에서는 파형 일치도의 평균 제곱근(root mean square; rms) 오차가 최소가 되게 하는 평활도 상충효과 매개변수(smoothness trade-off parameter)를 선택하는 Ammon et al. (1990)의 기준을 따라 최종적인 평활도 매개변수를 결정하였다. 따라서 본 연구에 대한 유전 알고리즘의 적합도 함수(fitness function)은 다음과 같다.

(7)$$fitness=C-\left[\sqrt{\sum _j{\left(O_j-C_j\right)}^2}+\omega \times roughness\right]$$

O_j와 C_j는 각각 i번째 샘플링 시간에서의 실제 수신함수와 유전 알고리즘에 의해 생성된 합성 수신함수를 의미하고 w는 평활도 상충효과 매개변수를 의미한다. 또한, 모델이 수렴할 때 값이 양수로 나오도록 하기 위하여 기준상수인 C를 채택하였다(Chang et al., 2004).

2.4. H-κ 중합법(H-κ Stacking)

Zhu and Kanamori (2000)는 수신함수의 직접 P파와 Ps 변환파, PpPs, PpSs+PsPs와 같은 다중반사파 사이의 전파 시간 차이와 미리 결정된 지각의 평균 P파 속도 Vp를 이용하여 지각의 두께 H와 V_p/V_s 비인 κ를 추정하는 방법을 제시하였다.

(8)$$H=\frac{t_{p_s}}{\sqrt{\frac{1}{V_s^2}-p^2}-\sqrt{\frac{1}{V_p^2}-p^2}}$$

(9)$$H=\frac{t_{p_p{p}_s}}{\sqrt{\frac{1}{V_s^2}-p^2}-\sqrt{\frac{1}{V_p^2}-p^2}}$$

(10)$$H=\frac{t_{p_p{p}_s+p_s{p}_s}}{2\sqrt{\frac{1}{V_s^2}-p^2}}$$

여기서 p는 입사한 파의 파선 변수(ray parameter)이며, t는 각 위상의 도달 시간을 의미한다. 식 (8) – (10)에서 V_p와 p는 이미 알고 있는 값이므로 H와 V_s 그리고 각 위상의 도달 시간 t 사이의 관계를 파악하기 위하여 H-κ 영역에서 격자 탐색법(grid search method)을 통해 아래의 식을 이용한 수신함수 중합을 진행한다.

(11)$$s\left(H,k\right)=w_{1}r\left(t_{P_s}\right)+w_{2}r\left(t_{P_p{P}_s}\right)-w_{3}r\left(t_{P_p{S}_s+P_S{P}_S}\right)$$

w는 가중치 인자, r(t)는 각 위상 도달 시간에서 법선성분 수신함수의 진폭을 의미한다. s(H, κ)는 올바른 H와 κ값으로 정확하게 세 가지 위상이 중합되었을 때 최대값을 가진다. 이후 결정된 H와 κ 값을 식 (8) – (10)에 대입하면 각 위상의 도달 시간 역시 추정할 수 있다. 그러나 H-κ 중합법은 수평적인 속도가 일정하다고 가정하기 때문에 수평적인 속도 변화가 존재할 경우 부정확한 결과를 보일 수 있으므로(Zandt and Ammon, 1995) 다양한 진앙거리와 후방위각을 가진 수신함수를 중합하여 이를 억제하고자 한다(Zhu and Kanamori, 2000).

본 연구에서는 강원도 지역과 그 주변 설치된 16개의 기상청(KMA) 관측소와 5개의 한국지질자원연구원(KIGAM) 관측소 총 21개의 광대역 관측소(Fig. 1)를 사용하였으며, 일부 관측소가 설치된 날짜를 고려하여 2019년 3월 18일부터 2022년 12월 31일까지 IRIS(Incorporated Research Institution for Seismology)가 제공하는 지진목록 중에서 모멘트 규모(Mw) 5.8 이상, 진앙거리 30° – 90°에서 발생한 원거리 지진자료 139개를 선별하였다(Fig. 2).

Figure 2. Locations of 139 teleseismic events and 6 groups for receiver function stacking. The concentric circles centered on the average of the locations of 21 stations represent 30-degree intervals of angular distances.

획득한 지진자료에 관측소와 지진 발생 위치 정보 등을 입력하고 SAC(Seismic Analysis Code)에서 제공하는 SSS(Signal Stacking Subprocess) 명령어를 이용하여 P파 도착시간을 계산하였다. 이때 기준 속도 모델은 iasp91(Kennett and Engdahl, 1991)이 사용되었다. 계산된 P파 도착시간을 기준으로 –10초부터 100초까지 잘라준 뒤 수신함수 계산을 위하여 3성분 자료를 후방위각 정보를 사용하여 각각 수직성분, 법선성분, 접선성분으로 변환하였다.

이때, 시추공 관측소의 경우 수평성분 센서가 진북방향으로부터 틀어져 설치된 경우가 존재하기 때문에 방위각을 보정해 줄 필요가 있다. 본 연구에서는 방위각 보정에 사용되고 있는 방법 중 하나로 레일리파 극성을 분석하는 Stachnik et al. (2012)의 보정 방법을 사용하였다. 레일리파는 수직 혹은 법선성분에만 기록되기 때문에 수평성분 자료의 방향이 정확하다면 이론적으로 레일리파의 수직성분 자료와 9 0° 위상 편이된 법선성분 자료는 선형 관계를 갖는다. 이러한 특성을 고려한다면 발생 위치를 알고 있는 지진에 대한 레일리파의 수직성분 자료와 위상 편이된 법선성분 자료 사이의 상호상관계수(crosscorrelation coefficient)를 이용하여 방위각 보정값을 측정할 수 있다. 지진자료의 레일리파의 신호만을 추출하기 위하여 평균값과 경향성, 기기응답함수를 제거하는 전처리 과정을 거친 후 Stachnik et al. (2012)가 제시한 레일리파의 이론적 평균 속도인 4 km/s로 계산된 도착시간을 기준으로 –20초부터 600초까지 잘라준 자료에 대하여 0.04 – 0.08 Hz의 버터워스 필터(butterworth filter)를 적용하였다. 이후 자료처리가 완료된 수직성분 자료와 법선성분 자료의 상호상관계수를 이용하여 각 관측소에 대한 방위각 보정값을 측정하였다.

이후 앞서 잘라준 110초 분량의 3성분 자료에 평균값 제거, 경향성 제거와 같은 전처리를 수행한 후 방위각 보정값을 적용하여 법선성분을 계산하고 이 법선성분으로부터 수직성분을 다음 식을 이용하여 곱풀기 하였다.

(12)$$R\left(t\right)={\displaystyle \int \frac{D_R\left(\omega \right){\overline{D}}_V\left(\omega \right)}{\phi \left(\omega \right)}}\times G\left(\omega \right)e^{jwt}d\omega$$

이때, $\overline{D_V}\left(w\right)$는 D_V(w)의 켤레복소수이고 수위값(water level)에 해당하는 ϕ(w)는 매개변수 c를 최소값으로 하는 D_V(w)의 자기상관(autocorrelation)이다. G(w)는 가우시안 필터(Gaussian filter)를 의미한다.

(13)$$\phi \left(\omega \right)=\max \left\{D_V\left(\omega \right){\overline{D}}_V\left(\omega \right),C\times \max \left\{D_V\left(\omega \right){\overline{D}}_V\left(\omega \right)\right\}\right\}$$

(14)$$G\left(\omega \right)=\exp \left(-\frac{{\omega }^2}{4{\alpha }^2}\right)$$

여기서 매개변수 c를 0.01로 설정하여 곱풀기 과정에서 값이 무한대가 되는 것을 방지하였고, α는 주파수 한계를 0.5Hz로 제한하는 2.5로 설정하여 고주파 잡음을 억제하였다(Langston, 1979; Chang et al., 2004).

이를 각 관측소마다 139개의 모든 지진자료에 대해 수행하여 법선성분과 접선성분 수신함수를 추출하였다. 추출된 수신함수 중에서 구조의 이방성으로 인한 산란 효과나 낮은 신호대잡음비로 인해 나타나는 접선성분의 수신함수(Chang et al., 2004)가 거의 나타나지 않고, 직접 P파와 모호면으로 인한 Ps파가 뚜렷하게 나타나는 법선 성분 수신함수를 양질의 수신함수로 정의하였고 이후 양질의 수신함수만을 사용하여 수신함수의 연합 역산과 H-κ 중합법 계산을 진행하였다.

4.1. H-κ 중합법 결과 및 분석

각 관측소마다 추출된 139개의 수신함수를 확인하여 직접 P파와 모호면으로 인한 Ps파가 뚜렷하게 나타나는 양질의 수신함수를 걸러내었고, 수평적인 속도 변화 효과를 억제하기 위하여 후방위각과 진앙거리에 관계없이 모두 중합하였다. H-κ 중합법에 대한 식 (8) – (10)의 V_p 값은 한반도 남부의 지각 평균 P파 속도인 6.3 km/s를 사용하였으며(Chang and Baag, 2006; Chang and Baag, 2007; Lee, 2010; Jeon et al, 2012; Cho et al., 2013), 식 (11)의 가중치 인자는 진폭에 비례하여 각각 w₁ = 0.7, w₂ = 0.2, w₃ = 0.1 을 사용하였다(Zhu and Kanamori, 2000; Chang and Baag, 2007; Jeon et al., 2012).

21개 관측소에 대한 H-κ 분석 결과 얻은 모호면 깊이에서 각 관측소의 고도를 빼준 뒤 위도 37.6° 기준으로 북부와 남부로 나누어 각각 Fig. 3와 Fig. 4에 도시하였다. 강원도 지역의 모호면 깊이는 24.9 – 33.2 km로 나타났으며 V_p/V_s 비는 1.695 – 1.760으로 나타났다. 그중 CHC2와 CHNB 관측소는 각각 30.7 ± 1.2 km, 1.715 ± 0.015와 29.8 ± 2.2 km, 1.720 ± 0.021로, Chang and Baag (2007)이 기존에 제시한 30.8 ± 3.7 km, 1.71 ± 0.06과 28.8 ± 2.8 km, 1.72 ± 0.06과 비교하여 9 5% 신뢰구간에서 일치하였다. H-κ 분석 결과를 보간하여 모호면 깊이와 V_p/V_s 비에 대한 분포도를 그렸을 때(Fig. 5), 동해안에 위치한 GEJB, YAYB 관측소의 V_p/V_s 비가 가장 높게 관찰되며 모호면 깊이도 동해안에서 얕고 태백산맥을 지나면서 깊어지는 특성을 보인다.

Figure 3. Results of H-κ stacking in and around northern Gangwon Province.

Figure 4. Results of H-κ stacking in and around southern Gangwon Province.

Figure 5. Distribution of (a) Moho depths (b) V_p/V_s ratios from H-κ stacking.

4.2. 수신함수와 표면파 분산의 연합 역산

강원도 지역 관측소 하부에 대한 1차원 지각속도구조를 획득하기 위하여 수신함수와 표면파 분산곡선 자료를 결합하고 유전 알고리즘을 이용하여 연합 역산을 진행하였다. 이때 각 관측소에 대한 표면파 분산곡선 자료는 Chang and Baag (2005)이 제시한 한반도 남부의 레일리파 위상속도 자료 중 거리상으로 가장 가까운 관측소에서의 자료를 사용하였다(Table 2). 그리고 H-κ 중합법을 통해 구한 V_p/V_s 비를 사용하였다.

Table 2 . Rayleigh-wave phase velocity data (Chang and Baag, 2005) used in this study.

Stations where Rayleigh-wave phase velocities were measured	Stations applied	Rayleigh wave phase velocity (km/s)
		25 s	20 s	15 s	10 s
CHC	CGUB
	CHC2	3.75	3.60	3.42	3.29
	GAPB
CHC	BURB
	JECB	3.74	3.60	3.43	3.28
	YOJB
CHNB	CHNB	3.76	3.62	3.33	3.17
	GMHB
DGY	JSB	3.75	3.59	3.42	3.29
	NAMB
	PYCB
	YAYB
KSA	GEJB	-	-	3.46	-
	KOSB
	KSA
	SH2B
	YKB
ULC	CHYB	3.76	3.67	3.52	3.26
	IMWB
	SND
	SNGB

방위각에 따른 지각 구조의 이방성에 대한 영향 고려와, 수신함수 중합을 통한 신호대잡음비 향상을 위하여 139개의 수신함수를 후방위각에 따라 6개의 그룹으로 나누었다(Fig. 2). 이때 중합 범위는 후방위각 20° 이내를 기준으로 하였다. 이 중에서 그룹 5는 후방위각 범위가 30° 로 넓게 나타났고 진앙거리도 30° 이상으로 넓게 퍼져 있으며, 그룹 6은 신호대잡음비가 낮아서 최대 4개의 양질의 수신함수만을 추출할 수 있기 때문에, 수신함수 중합을 진행하기에 부적합하다고 판단하였고 연합 역산대상 그룹에서 제외하였다. 따라서 각 관측소마다 그룹 1부터 그룹 4중에서 최소 6개의 양질의 수신함수를 추출하는 3개의 그룹을 선정하고 각 그룹에 대한 그룹 속도 모델을 획득한 뒤 3개의 그룹 속도 모델을 중합함으로써 최종 속도 모델을 결정하였다. 이때 여러 그룹이 같은 개수의 양질의 수신함수를 추출하는 경우 수신함수 역산을 진행하였을 때 파형 일치도의 rms 오차가 더 작은 그룹을 우선적으로 선택하였다(Table 3).

Table 3 . Selected groups and number of high-quality receiver functions for each group.

Station name	Group used / Number of high-quality receiver functions
	Group name
	Group 1	Group 2	Group 3	Group 4
BURB	O / 7	X / 6	O / 10	O / 8
CGUB	O / 10	X / 7	O / 11	O / 10
CHC2	O / 9	X / 8	O / 10	O / 8
CHNB	O / 6	O / 6	O / 8	X / 4
CHYB	O / 9	X / 7	O / 9	O / 7
GAPB	O / 7	X / 5	O / 9	O / 6
GEJB	X / 6	O / 6	O / 10	O / 6
GMHB	O / 7	X / 4	O / 9	O / 6
IMWB	O / 7	X / 6	O / 8	O / 8
JECB	X / 6	O / 7	O / 9	O / 7
JSB	X / 6	O / 7	O / 8	O / 7
KOSB	X / 7	O / 7	O / 10	O / 7
KSA	X / 6	O / 7	O / 7	O / 8
NAMB	X / 8	O / 6	O / 10	O / 7
PYCB	X / 8	O / 7	O / 8	O / 10
SH2B	O / 8	X / 5	O / 10	O / 6
SND	O / 6	O / 6	O / 8	X / 6
SNGB	O / 8	X / 4	O / 10	O / 8
YAYB	O / 8	O / 7	O / 7	X / 4
YKB	O / 10	X / 5	O / 10	O / 6
YOJB	O / 10	O / 6	O / 8	X / 6

Chang et al. (2004)의 방법론에 따라 유전 알고리즘의 초기 인구를 50으로 설정하고 이를 총 300세대 진화시켰다. 이후 평활도 상충효과 매개변수를 0.01부터 0.5까지 로그 스케일(logarithmic scale)로 대입하여 파형 일치도의 rms 오차가 최소가 되는 값을 선택한 뒤 이를 고정시킨 후, 초기 모델을 변경하며 유전 알고리즘을 5회 반복하여 총 5개의 최적 속도 모델을 획득하였다. 5개의 최적 속도 모델을 평균하여 그 그룹에 대한 최적 속도 모델로 결정하였다. 위 과정을 각각의 그룹에 대해서 실행하였으며, 관측소에 대한 최종적인 속도 모델을 결정할 때는 3개의 그룹 최적 속도 모델을 결정할 때 사용되었던 총 15개의 최적 속도 모델을 평균하여 구하였다(Fig. 6).

Figure 6. (a) An example of how the final velocity model is constructed for the station PYCB. One-dimensional final velocity models are obtained for each group (Fig. 2) by averaging the 5 best velocity models from 5 runs of the genetic algorithm, and the final model for the station is constructed by averaging the 15 best velocity models from the 3 best groups (in this case, groups 2, 3, and 4). The solid and dashed lines in the velocity model panels show the average velocity model and the standard deviations of the velocities in each layer, respectively. The solid and dashed lines in the receiver function panels indicate the observed and synthetic receiver functions, respectively. (b) Comparison of observed and synthetic Rayleigh-wave phase velocities from the final model.

4.3. 연합 역산 결과 및 분석

강원도 지역과 그 주변에 설치되어 있는 21개의 광대역 관측소에 대하여 수신함수와 표면파 분산곡선의 연합역산을 시행하였다. 이때 모호면의 깊이는 처음으로 P파 속도가 7.5 km/s를 초과하는 층의 윗부분으로 결정하였다. 수신함수를 이용하여 P파 속도구조를 도출한 여러 연구가 7.5 km/s를 기준으로 모호면 깊이를 결정한 바 있고(Doloei and Roberts, 2003; Chang et al., 2004; Chang and Baag, 2005; Park et al., 2009 ), S파 속도구조의 경우도 일반적으로 사용되는 $V_p=\sqrt{3}{V}_s$으로 계산되는 4.3 km/s를 기준으로 하는 많은 연구들이 존재하는 것을 고려하여 모호면 깊이를 결정하였으며(Özalaybey et al., 1997; Jeon et al., 2012; Jeon et al., 2013; Kgaswane et al., 2009), 본 연구의 경우에도 7.5 km/s를 기준으로 사용하였을 때 H-κ 중합법과 유사한 모호면 깊이 결과값을 보여주었다.

연합 역산의 결과도 H-κ 분석과 같은 방법으로 모호면 깊이에서 각 관측소의 고도를 빼준 뒤 위도 37.6° 기준으로 북부와 남부로 나누어 각각 Fig. 7과 Fig. 8에 도시 V 하였다. 이에 따른 강원도 지역의 모호면 깊이는 25.9 – 33.7 km로 전반적으로 북부보다 남부에서의 모호면 깊이가 더 깊게 나타났다. 특히 CHNB, KSA, CHC2 관측소는 각각 29.8 km, 25.9 km, 31.7 km로 나타났으며 이는 기존에 Chang and Baag (2005)이 제시한 27.8 km, 25.9 km, 31.8 km와도 비교적 잘 일치한다.

Figure 7. The final velocity models in and around northern Gangwon Province. The results are sorted in the ascending order of Moho depth. Black, red and blue dashed lines indicate Moho depths, velocity inversion layers, and low-velocity layers, respectively.

Figure 8. The final velocity models in and around southern Gangwon Province. The results are sorted in the ascending order of Moho depth. Black, red and blue dashed lines indicate Moho depths, velocity inversion layers, and low-velocity layers, respectively.

역산 결과를 검증하기 위해 각 관측소마다 H-κ 중합법으로 추정한 모호면의 깊이와 연합 역산으로 추정한 모호면의 깊이를 비교했을 때의 표준편차는 ± 0.9 km로 나타났으며(Table 4), 이는 연합 역산에 사용된 속도 모델의 층의 두께가 모호면 근처에서 2 km임을 고려할 때, 그리고 모호면의 기준을 임의로 7.5 km/s로 정한 것을 고려할 때, 대체적으로 잘 일치하고 있음을 보여준다. 이후 연합 역산 결과를 보간하여 모호면 분포도를 그렸다(Fig. 9). 이를 H-κ의 분포도와 비교한 결과, 대체로 일치하는 양상을 보인다.

Table 4 . Differences in Moho depths from H-κ stacking and joint inversion.

Station name	Vp/Vs ratio	Moho depth from H-κ stacking (km)	Moho depth from Joint inversion (km)	Difference (km) (H-κ - Joint inversion)
BURB	1.745	30.9	31.9	-1.0
CGUB	1.740	29.8	29.8	-
CHC2	1.715	30.7	31.7	-1.0
CHNB	1.720	29.8	29.8	-
CHYB	1.705	31.1	31.6	-0.5
GAPB	1.730	30.9	31.9	-1.0
GEJB	1.760	27.9	27.9	-
GMHB	1.695	30.8	29.8	+1.0
IMWB	1.710	27.0	28.0	-1.0
JECB	1.745	32.2	33.7	-1.5
JSB	1.695	33.2	33.7	-0.5
KOSB	1.720	24.9	25.9	-1.0
KSA	1.720	26.9	25.9	+1.0
NAMB	1.745	29.8	31.3	-1.5
PYCB	1.725	32.2	33.7	-1.5
SH2B	1.695	30.1	31.6	-1.5
SND	1.715	30.7	31.2	-0.5
SNGB	1.725	31.9	31.9	-
YAYB	1.760	26.5	26.0	+0.5
YKB	1.735	30.2	29.7	+0.5
YOJB	1.735	32.7	31.7	+1.0
			Standard deviation	±0.9 km

Figure 9. Distribution of Moho depths from the joint inversion of receiver functions and surface wave dispersion.

본 연구는 강원도 지역의 지각속도구조를 분석하기 위해 현재 강원도와 그 주변에서 운영되고 있는 21개의 광대역 관측소에서 획득한 지진자료에 H-κ 중합법과 수신 함수와 표면파 분산의 연합 역산을 적용하여 다음과 같은 결론을 얻었다.

H-κ 중합법 적용 결과 강원도 지역의 모호면 깊이는 24.9 – 33.2 km로 나타났고, V_p/V_s 비는 1.695 – 1.760으로 나타났다. 일부 관측소는 기존에 제시된 결과와 9 5% 신뢰구간에서 일치하였다. 그중 동해안에 위치하고 있는 GEJB와 YAYB 관측소의 V_p/V_s 비가 1.760으로 내륙 지역에 비해 상대적으로 높게 나타났다. 지각을 이루는 암석의 종류가 고철질(mafic)이거나(Musacchio et al., 1997) 지각 내부에 유체가 존재할 경우(Wang et al., 2012) V_p/V_s비가 높게 나타나는 요인이 된다. GEJB와 YAYB 관측소가 해안가에 위치한 것을 고려하면 해수의 영향으로 인해 S파 속도가 감소하여 상대적으로 높은 V_p/V_s 비가 계산된 것으로 보일 수 있으나 해수는 매우 얕은 부분에서만 영향을 미칠 것으로 보이며 이는 지각 평균의 V_p/V_s 비가 높게 나타나는 요인이라고 보기 어렵다고 생각된다. 게다가 같은 해안가에 위치한 KSA, KOSB 관측소와 IMWB 관측소는 각각 1.720, 1.710로 비교적 낮은 V_p/V_s비가 나타난 것으로 볼 때 해수로 인한 영향보다는 암석의 종류에 의한 영향일 가능성이 더 높을 것으로 생각된다.

H-κ 중합법으로 얻은 V_p/V_s 비를 이용하여 수신함수와 표면파 분산의 연합 역산을 진행한 결과 모호면의 깊이는 25.9 – 33.7 km로 나타났으며, 이를 검증하기 위하여 각 관측소마다 H-κ 중합법과 연합역산으로 각각 추정한 모호면 깊이를 비교했을 때의 표준편차는 ± 0.9 km로 나타났다. 이는 연합 역산에 사용된 속도 모델의 모호면 근처에 대한 층의 두께가 2 km임을 고려하였을 때, 그리고 모호면의 기준을 임의로 7.5 km/s로 정한 것을 고려할 때, 두 독립적인 결과가 대체적으로 공간적인 일치를 잘 보여주고 있다. 이는 V_p/V_s 비의 영향이 그다지 크지 않다는 것을 의미하는데, H-κ 중합법으로 구한 V_p/V_s 비의 범위가 0.1 이내로 분포하는 것으로 보아(Fig. 4b) 예상되는 결과이다. 또한, 연합 역산 결과를 보간하여 모호면 분포도를 그린 후 H-κ의 분포도와 비교하였을 때 대체로 동해안에서 얕고 태백산맥을 지나면서 깊어지는 일치된 양상을 보인다.

중력과 지형 자료를 바탕으로 유효 탄성 두께를 이용하여 모호면 깊이를 계산한 Shin et al. (2006)에서도 태백산맥을 중심으로 동쪽은 얕고 서쪽은 깊은 구조를 보였으며, 기존에 수신함수를 사용하여 한반도 남부의 모호면 등고선을 제시한 Chang and Baag (2005)에서도 동해안의 모호면이 가장 얕은 것으로 나타났다. 이를 통해 강원도 지역의 모호면이 에어리(Airy)의 지각평형설(isostasy)을 대체적으로 따르고 있음을 보여준다.

강원도 남부에 위치한 CHYB, JECB, JSB, PYCB, YOJB 관측소는 산악지역에 설치되어 있음에도 불구하고, 직하부 2 km 깊이까지 P파 속도 5 km/s 이하의 저속도층이 관찰되었다. 이는 관측소의 위치를 고려하였을 때 태백산 분지의 영향으로 보인다. 태백산 분지는 고생대에 형성된 퇴적분지로서 이에 속해 있는 조선누층군은 퇴적당시 기반암의 블록 지구조 운동(block tectonics)으로 인해 1500 m에 달하는 석회암층이 쌓여 있는 것으로 보고된 바 있다(Yun, 1978). 강원도 북부의 CHNB, GAPB, SH2B 관측소 직하부 2 km 깊이에서도 이러한 저속도층이 관찰되었으며, CHNB 관측소가 위치한 철원군과 GAPB 관측소가 위치한 가평군은 한국지질자원연구원이 제공하는 지질도에 의하면 신생대 제 4기의 충적층이 존재하고 있으며 특히 Hwang(2010)은 철원군의 충적층이 백악기의 화산암류와 관입암류를 덮으면서 신생대 제 4기의 현무암으로 이루어진 용암대지 위에 쌓여 있다고 보고하였다. 하지만 SH2B 관측소가 위치한 서화면은 화강암과 호상편마암을 기반암으로 하는 소양강의 상류부 유역에 위치하여 있음에도 불구하고 5 km/s 이하의 저속도층이 나타났는데 이는 기반암의 화학적 풍화나 균열로 인하여 지진파 속도가 감소하였을 가능성이 있다(Clair et al., 2015). 연합 역산으로 도출된 지각속도구조는 강원도 지역의 4 – 12 km 깊이 사이에 깊이가 깊어질수록 P파 속도가 감소하는 속도역전층이 존재하고 있음을 암시한다. 이는 Kim et al. (2010)이 경기도 강화와 경상북도 영덕을 잇는 2 9 9 km 길이의 KCRT(Korean Crust Research Team)-2008 측선을 따라 인공지진파 초동주시 역산을 수행하여 옥천대와 경기육괴 하부에 깊이 7.2 km에 중심을 둔 저속도층이 해수면 하부 11.7 km 깊이까지 존재한다는 것을 발견한 결과와 유사하며, Kim et al. (2007)이 이전에 KCRT-2002, KCRT-2004 측선을 따라 실시한 굴절파 탐사에서 발견한 저속도층의 깊이와도 유사한 결과이다. 또한 CGUB, GAPB, JSB, PYCB, SND, SNGB, YKB 관측소의 하부 10 km 깊이에서 속도가 가파르게 상승하는 불연속면을 관찰하였다. 이러한 중간 지각에서의 불연속면은 대부분 암석의 특성이 변하는 경계로 볼 수 있고(Chang and Baag, 2005), 불연속면의 직전 깊이에 속도역전층이 4 – 5 km 두께로 존재하고 있음으로 보아 해당 관측소 하부 10 km 부근에 밀도차를 보이는 암석의 경계가 존재하고 있을 것이라고 생각된다. 이러한 한반도 남부의 중간지각 불연속면은 Kim and Kim (1983), Cho et al. (2006), Yoo et al. (2007)등의 연구에서 본 연구보다는 약간 더 깊은 15 km 깊이 부근에 존재함이 지속적으로 보고되어 왔다.

본 연구는 연합 역산 과정에서 Chang and Baag(2005)의 표면파 분산곡선 자료를 사용하였다. 당시에는 강원도 지역에 설치된 광대역 관측소의 수가 적어 본 연구에 필요한 모든 관측소에 대한 레일리파 위상속도를 얻을 수 없었고 따라서 가장 가까운 관측소에서의 측정값을 사용하였기 때문에 다소 부정확한 결과가 나타났을 가능성이 존재할 수 있다. 그러나 표면파 분산곡선 자료가 수평적으로 변동성이 크지 않은 특성을 고려하면 그 영향은 크지 않을 것으로 생각된다. 하지만 이후 연구에서는 본 연구에 사용된 모든 관측소에 대한 표면파 분산곡선 자료를 취득하여 속도 모델에 대한 정확도를 높일 필요성은 있다.

본 연구는 원거리 지진파를 사용하여 인공파원을 사용하는 반사법 혹은 굴절법 탐사보다 경제적이면서도 속도 변화에 민감하여 지진파 토모그래피보다 높은 수직 해상도를 가질 수 있는 수신함수를 사용하여 강원도 지역 하부의 지각속도구조를 도출하였다. 특히, Chang et al. (2004)의 방법을 사용하여 수신함수가 가지는 가장 큰 문제점인 모호면 깊이와 지각속도에 대한 상충효과를 표면파 분산 자료와의 연합 역산으로 억제하였고, 초기 모델에 영향을 받지 않는 유전 알고리즘을 사용함으로써, 초기 속도구조에 대한 정보가 부족한 강원도 일대에서 21개의 광대역 지진관측소 하부에 대한 높은 해상도의 1차원 지각속도구조를 도출하였다. 모호면 깊이는 25.9 – 33.7 km로 대체적으로 에어리의 지각평형설을 따르는 양상을 보인다. 1차원 지각속도 모델 해석 결과 태백산 분지와 강원도 북부 지표면 부근에서 저속도층이 관측되었고, 모든 관측소의 4 – 12 km 깊이 사이에서 속도역전층이 발견되었으며, 일부 관측소 하부에서는 중간지각 불연속면이 발견되었다.

이 논문은 행정안전부의 방재안전분야 전문인력 양성사업의 지원을 받았으며, 2023년도 정부(교육부)의 재원으로 한국연구재단의 지원을 받아 수행된 기초연구사업(No.2019R1A6A1A03033167), 2023년도 정부(과학기술정보통신부)의 재원으로 한국연구재단의 지원을 받아 수행된 기초연구사업(No.2019R1A2C208506111)입니다.