서 론

자력 탐사는 지하 매질의 대자율 대비에 의해 발생하는 자기 이상을 이용하여 광체 및 지질 구조를 탐지하는 효과적인 물리탐사 방법이다. 자기 이상은 관측점과 이상체 간 거리 증가에 따라 급격히 감소하므로, 지표 자력 자료만으로는 심부 자기 이상체의 정확한 위치와 형태를 결정하기 어렵다. 이를 보완하기 위해 시추공에서 관측한 자력 자료가 활용되고 있으며, 이를 지표 자력 자료와 함께 해석하기 위해 다중 공간 자료에 대한 복합 역산이 활용되고 있다(Li and Oldenburg, 2000). 다만, 실제 자력 탐사 현장에서는 비용 및 작업 환경의 제약으로 다수의 시추공을 확보하기 어려운 경우가 많다. 이에 따라 이 연구에서는 활용 가능한 시추공의 개수가 1개일 경우, 해당 시추공의 최적 위치를 결정하기 위한 방법을 제안하고자 한다.

지표 자력 자료와 시추공 자력 자료를 결합한 다중 공간 자력 자료의 복합 역산은 이상체의 수직 정보를 보완함으로써 역산 성능을 향상시키는 효과적인 방법으로 사용되었다(Li and Oldenburg, 2000; Zhang et al., 2025; Shi et al., 2025). 그러나 기존 연구는 대부분 다수의 시추공을 가정하고 있으며, 단일 시추공 조건에서의 연구는 수행되지 않았다. 실제 자원탐사 현장에서는 예산의 제약이나 환경적 요인으로 인해 다수의 시추공을 확보하는 것은 한계가 있다. 이 연구에서는 단일 시추공만을 확보할 수 있는 상황을 가정하고 시추공 위치 변화에 따른 복합 역산 결과의 차이를 정량적으로 분석하여 효율적인 탐사 설계를 위한 최적의 시추공 위치를 제안하고자 한다.

이론적 배경

자력 탐사에서 관측되는 자기 이상은 지하 매질의 대자율 분포와 지자기장, 그리고 관측점과 이상체 간의 거리에 의해 결정된다. 대자율 대비를 갖는 이상체는 자기장에 의해 유도 자화가 발생하며 이로 인해 주변 자기장 분포에 변화가 나타난다. 이러한 자기 이상은 관측점과 이상체 사이의 거리가 증가할수록 급격히 감쇠하는 특성에 의해 민감도가 크게 달라진다.

지표 자력 자료는 관측점이 지표 또는 상공에 분포하여 지표면에 대해 수평적으로 측정하여, 수평 방향의 이상체 분포에 대한 정보를 효과적으로 제공하지만 이상체의 심도나 두께와 같은 수직 방향 정보에는 상대적으로 한계를 가진다. 반면 시추공 자력 자료는 관측점이 지하에 설정되어 지표면에 수직한 방향으로 측정한다. 이로 인해 시추공 자력 자료는 이상체의 수직 방향 변화에 대한 민감도가 높으며 지표 자력 자료에서 부족한 정보를 보완할 수 있다(Li and Oldenburg, 2000). 따라서 지표 자력 자료와 시추공 자력 자료는 서로 상호보완적인 특성을 가지며 두 자료를 결합한 다중 공간 자력 자료는 지하 대자율 구조에 대한 해석의 신뢰도를 향상시킬 수 있다.

자력 역산을 위해서는 지하 매질을 다수의 요소로 분할하고, 각 요소 내의 대자율이 일정하다고 가정한다. 또한 각 요소는 잔류자화(remanent magnetization)를 가지지 않는다고 가정한다. 이 경우에 자력이상은 지하 매질의 대자율과 선형 관계를 가지며, 다음과 같은 관계를 갖는다.

여기서 d = (d₁, ..., d_N)^T는 N개의 측점에서 관측된 총 자력이상 값으로 구성된 열 벡터이고, m = (m₁, ..., m_M)^T는 M개로 이루어진 요소의 대자율 열 벡터이다. 행렬 G는 크기 N × M인 감도 행렬로서 각 요소의 대자율이 자력 이상 자료에 미치는 영향을 의미한다.

관측된 자력 이상에 부합하는 대자율 모형을 찾기 위해 본 연구에서는 역산에 대한 목적함수(φ)를 다음과 같이 설정하고, 이를 최소화하는 모델 변수를 추정하였다.

여기서 φd는 데이터 오차항(misfit)에 해당하고, φm은 모델 목적 함수이다. β는 데이터 오차항과 모델 목적함수 항 사이의 상대적인 중요도를 조절하는 가중치이다.

식 (3)은 L2 norm으로 표현되며, 여기서 G민감도 행렬, d_obs는 관측 데이터, W_d는 가중치를 결정하는 대각행렬이다.

모델 목적함수는 다음과 같이 Tikhonov 정규화를 통해 계산하였다(Li and Oldenburg, 1996; Li and Oldenburg, 2000).

여기서 w_s, w_x, w_y, w_z는 공간적으로 변화하는 가중함수로 첨자 s는 모델 파라미터의 크기를 조절하는 항이고, x, y, z는 공간 축방향 항이다. 이때, x는 북쪽을, y는 동쪽을, z는 연직 하향 방향으로 설정하는 직교좌표계를 이용하였다. α_s, α_x, α_y, α_z는 목적함수에서 각 항의 상대적 중요도를 조절하는 가중치 계수이다. w(r)은 깊이 가중함수(depth weighting function) 혹은 거리 기반 가중함수(distance weighting function)을 사용하였고, 이는 민감도가 거리에 따라 감쇠하는 효과를 보정하기 위해서 사용하였다.

모델 설정

역산 공간의 총 크기는 4 km × 4 km × 1 km에 영역으로 설정하였으며, 각 요소는 x 축과 y 축 방향으로 100 m크기로 분할하였으며, z축 방향은 50 m 크기로 설정하였다. 따라서 역산 공간에 해당하는 모델변수(직육면체)의 총 수는 32,000 (40 × 40 × 20)개이다. 합성 모델 이상체는 반지름 1000 m, 윗면의 깊이는 300 m, 두께 300 m의 원판 형태이고, 대자율은 0.05 (SI)로 설정하였다(Fig. 1). 시추공 위치에 따른 자력 반응의 차이만을 확인하기 위해 원판 형태의 이상체를 설정하여 이상체의 형태에 따른 영향을 최소화하였다.

Fig. 1

Disk model. Its radius and thickness are 1000 m and 300 m, respectively. The magnetic susceptibility is 0.05 (SI).

자력 반응은 균질한 대자율을 갖는 직육면체 프리즘(prism)으로 계산하였고, 직육면체 셀에 의한 자력 반응은 포아송(Poisson) 관계식(Blakely, 1996)을 포함한 식 (5)를 이용하였다.

여기서 γ는 중력 상수, ρ는 밀도, κ는 대자율, m은 자력 모멘트 벡터이다. U는 직육면체 프리즘에 의한 중력 퍼텐셜로서 Nagy et al. (2000)과 Fukushima (2020)가 제안한 식을 이용하여 계산하였다.

원판 이상체가 유발하는 합성 자력 이상은 동서 100 m, 남북 600 m 간격으로 구성된 격자 측점에 대해 계산하였다(Fig. 2(a)). 시추공은 원판 모델을 중심으로 1,150 m 떨어진 위치에서 방위각 15° 간격으로 총 24개로 분포시켰다(Fig. 2(b)). 모든 합성 자력 값은 공통적으로 표준 편차가 10 nT인 무작위 잡음을 더하였다. 이때 유도 자화를 발생시키는 지구자기장은 한반도 지구자기장 세기 및 복각인 49,550 nT, 복각 +51°이고 편각은 해석의 편의를 위하여 0°로 가정하였다.

Fig. 2

Horizontal view of the disk model. The black triangles in panel (a) represent the surface observation points, and those in panel (b) indicate the borehole locations.

이 연구에서는 복합 역산에서 시추공의 위치가 역산 해상도에 미치는 영향을 정량적으로 평가하고자 하였다. 이를 위해 지표 자력 자료와 단일 시추공 자력 자료를 결합하여 복합 역산을 수행하였다. 이 연구에서 실험한 모든 역산 결과들은 역산으로 추정한 자력 이상이 자력 이상 참값에 포함된 잡음 수준 내에서 자력 이상을 재현하고 있음을 확인하였다(Fig. 3과 4).

Fig. 3

Comparison of the total magnetic intensity (TMI) for the disk model. (a) Observed TMI; (b) predicted TMI obtained from the joint inversion using a borehole located at 0°; and (c) normalized misfit. The observed TMI are contaminated with Gaussian noise of 10 nT.

Fig. 4

Observed borehole total magnetic intensity (TMI) (blue point) and predicted borehole TMI (red line) obtained through the joint inversion using a borehole located at 0°. The observed TMI are contaminated with Gaussian noise of 10 nT.

역산 결과의 정량적 평가 방법

시추공 위치에 따른 복합 역산 결과를 비교하기 위해, 복원된 이상체 영역과 실제 이상체 영역 간 공간적 유사도를 정량적으로 평가하였다. 일반적으로 역산 결과는 단면도 또는 3차원 분포도를 통해 시각적으로 비교할 수 있으나, 이러한 비교는 해석자의 주관이 개입될 가능성이 있으며, 특히 시추공의 위치 변화에 따른 여러 역산 결과를 통해 최적의 해를 결정하는데 한계가 있다. 이에 이 연구에서는 정량적 평가 방법을 도입하여 역산 결과를 비교하였다. 역산 결과에 기초해 이상체의 경계를 특정하기 위해 임계값(threshold) 설정이 필요하다. 임계값을 설정하기 위해 Otsu (1979)의 방법을 사용하였다. Otsu 방법은 대자율의 히스토그램을 이용하여 배경과 이상체를 두 개의 클래스로 구분하고, 두 클래스 간의 분리도가 최대가 되는 지점을 임계값으로 선정하는 방식이다. 역산 결과, 대자율 대비가 임계값보다 높은 영역을 자기 이상체로 정의하였다(Fig. 7).

역산 결과로 복원된 대자율 모델과 원판 이상체 모델을 각각 이상체 영역과 배경 영역으로 분류한 후, 혼동행렬(confusion matrix)을 이용하여 정량적인 비교를 수행하였다(Table 1). 원판 이상체 모델과 역산 결과 모두에서 이상체로 분류된 영역은 True Positive (TP)로 정의하였으며, 원판 이상체 모델에서는 이상체 영역이지만 역산 결과에서는 배경으로 분류된 영역은 False Negative (FN)로 정의하였다. 반대로 원판 이상체 모델에서는 배경이지만 역산 결과에서 이상체로 분류된 영역은 False Positive (FP)로 정의하였다. 혼동행렬의 값을 이용하여 IoU (Intersection of Union), Dice 계수, 정밀도(Precision), 재현율(Recall)을 식 (6),(7),(8),(9)를 통해 계산하였다.

IoU는 원판 이상체 모델과 역산 결과의 이상체로 분류된 영역이 서로 겹치는 정도를 전체 이상체 영역 대비 비율로 나타내는 지표로, 이상체의 공간적 일치도를 종합적으로 평가하는 데 사용된다. 따라서, 시추공 위치 변화에 따른 복합 역산 결과가 실제 이상체의 공간적 범위를 얼마나 정확하게 복원하는지 판단하는 지표로 사용하였다. Dice 계수는 이상체 영역 간의 중첩 정도를 강조하는 지표로, 두 영역의 일치성이 증가할수록 높은 값을 나타낸다. 정밀도는 역산 결과에서 이상체로 분류된 영역 중 실제 이상체가 차지하는 비율을 의미하며, 역산 결과의 과대 추정 경향을 평가하는 데 유용하다. 반면 재현율은 실제 이상체 영역 중 역산 결과에서 재현된 비율을 의미하며, 이상체를 얼마나 놓치지 않고 복원했는지를 나타낸다. 이 연구에서는 이상체의 공간적 일치성과 과대 추정 여부를 우선적으로 평가하기 위해 IoU, Dice 계수 및 정밀도를 주요 판단 지표로 사용하여 역산 결과의 우수성을 평가하였다. 재현율은 이상체의 누락 여부를 보조적으로 확인하는 지표로 활용하여, 주요 지표들과 함께 종합적으로 해석하였다.

시추공 위치에 따른 역산 결과

지자기장 편각이 0°인 경우, 시추공의 위치가 0°와 180°인 수평 단면(Fig. 5(a), 5(c))과 수직 단면(Fig. 6(a), 6(c))에서는 실제 이상체의 형태와 역산 결과가 일치하는 것을 확인할 수 있다. 이를 통해 시각적으로 평가할 때 이상체의 형태와 경계가 비교적 명확하게 복원된 것으로 보인다. 반면, 시추공의 위치가 90° (Fig. 5(b), 5(d))와 270°(Fig. 6(b), 6(d))에서는 실제 이상체의 형태와 복원된 모델의 이상체의 경계가 일치하지 않아, 역산 결과가 비교적 좋지 않다는 것을 확인할 수 있다.

Fig. 5

Horizontal slices of the recovered susceptibility models at the middle of the disk (z = −475 m) for borehole locations of (a) 0°, (b) 90°, (c) 180°, and (d) 270° when the declination is 0°. The white lines indicate the boundaries of the true model.

Fig. 6

Vertical slices of the recovered susceptibility models along y = 50 m for borehole locations of (a) 0°, (b) 90°, (c) 180°, and (d) 270° when the declination is 0°. The black rectangles indicate represents the boundaries of the true model.

역산 결과의 정량적 평가를 위해 원판 이상체 모델과 역산 결과의 이상체와 배경 분리 기준을 설정하였다. 이를 위해 역산된 대자율 값 분포를 분석하고, 이상체와 배경을 구분하는 임계값을 Fig. 7과 같이 계산하였다. Fig. 7에서 파란색 곡선은 역산을 통해 계산된 대자율 수치의 확률 분포를 로그 스케일로 나타낸 것이며, 주황색 곡선은 각 대자율 값에 대한 클래스 간 분산을 나타낸다. 역산을 통해 계산된 대자율 수치는 배경 영역의 저대자율 값과 이상체에 해당하는 고대자율 값이 혼재된 분포를 가지며, 이에 따라 대자율 값의 히스토그램은 두 클래스(이상체와 배경)를 포함하는 형태를 보인다. 이 연구에서는 이러한 분포 특성을 기반으로, 클래스 간 분산이 최대가 되는 지점을 분리 기준으로 설정하였다.

Fig. 7

Histogram of the recovered susceptibility model recovered from the joint inversion using a borehole located at a declination of 0°. The threshold for distinguishing the target body from the background was determined using the Otsu method, which maximizes between-class variance.

시추공 위치 변화에 따른 이상체 영역의 유사도를 IoU, Dice 계수, 정밀도, 재현율의 네 가지 지표로 정량적으로 평가하였다. Fig. 8은 지자기장의 편각이 0°인 조건에서 시추공 위치의 방위각 변화에 따른 각 유사도 지표의 변화를 나타낸다. IoU와 Dice 계수는 원판 이상체 모델과 역산 결과로 복원된 이상체 영역 간의 공간적 중첩 정도를 평가하는 지표로, 시추공이 이상체를 기준으로 방위각 0° 인 곳에 위치할 때 가장 높은 값을 나타내며(Fig. 8(a), 8(b)), 이는 역산의 결과가 참값과 가장 근접하다는 것을 의미한다. 정밀도 또한 최대값을 보이며(Fig. 8(c)), 이는 복원된 모델의 이상체 영역 중 실제 이상체가 차지하는 비율이 가장 높음을 나타낸다. 반면, 재현율은 시추공이 이상체를 기준으로 방위각 180° 인 곳에 위치할 때 상대적으로 높은 값을 보이는데(Fig. 8(d)), 이는 실제 이상체의 상당 부분이 복원되었음을 의미한다. IoU와 정밀도가 낮게 나타나는 것은 역산이 과대 추정되었을 가능성을 시사한다. 종합적으로 볼 때, IoU, Dice 계수 및 정밀도가 동시에 최대가 되는 시추공 0° 방향이 역산 결과가 가장 우수하다고 해석할 수 있다.

Fig. 8

Quantitative evaluation at a declination of 0° and an inclination of 51°.

토의 및 결론

이러한 역산 결과는 시추공의 최적 위치를 결정함에 있어 외부 자기장의 편각이 반드시 고려되어야 함을 의미한다. 이를 확인하기 위해 지자기장의 편각을 15°로 바꿔 복합 역산을 수행하였다. Fig. 9에서 나타난 바와 같이, 지자기장의 편각이 15°인 조건에서는 시추공이 이상체를 기준으로 방위각 15°인 곳에 위치했을 때 IoU, Dice계수 및 정밀도가 가장 높은 값을 나타난다. 반면, 재현율은 다른 각도에 비해 다소 낮은 값을 보이지만, 이런 결과는 이상체가 누락된 것이 아닌 역산 과정에서 이상체의 경계를 보다 정밀하게 복원되었음을 시사한다. 즉, 지자기장의 편각과 시추공 위치의 방위각이 일치하는 조건에서 이상체의 분포가 정확하게 복원된다. 반면, 불확실성이 높은 경계 영역이 과도하게 포함되는 것을 방지하는 방향으로 역산 결과가 수렴한 것으로 해석할 수 있다. 시추공 위치를15° 간격으로 변화시킨 모든 결과에서 지자기장 편각 간의 동일한 연관성을 확인할 수 있었다.

Fig. 9

Quantitative evaluation of the inversion results at a declination 15° and an inclination of 51°.

앞선 결과를 통해, 시추공이 지자기장의 편각 방향에 위치할 때 역산 결과의 유사도가 가장 높게 나타나는 경향을 확인하였다. 이러한 경향이 단순히 지자기장 편각에 의한 효과인지 확인하기 위해 복각에 의한 영향을 추가적으로 분석하였다. 이를 위해 지자기장의 복각을 –51°로 편각은 0°로 설정하여 복합 역산한 결과, 시추공이 편각의 반대 방향인 180°에 위치할 때 유사도가 가장 높게 나타났다(Fig. 10). 즉, 시추공이 지자기 편각에 평행한 방향에 위치할 때 역산 결과의 유사도가 가장 높게 나타났으며, 북반구에서는 편각 방향, 남반구에서는 편각의 반대 방향으로 나타나고 둘 다 지하 이상체 방향을 가리킨다.

Fig. 10

Quantitative evaluation of the inversion results at a declination of 0° and an inclination of -51°, corresponding to the Southern Hemisphere.

정량적 평가 결과, 역산 결과의 유사도가 가장 높게 나타나는 시추공의 위치는 확인할 수 있었지만 원인에 대한 직접적인 분석에는 한계가 있었다. 이에 시추공의 위치에 따른 자력 반응의 차이를 통해 관찰된 결과를 해석하고자 하였다. 시추공 자력 이상은 이상체와의 기하학적 위치에 따라 자력 이상의 진폭이 크게 달라진다. 시추공 자력 이상만을 이용해 역산한 결과와 시추공 자력 이상의 절대적인 크기를 나타내는 평균제곱근(Root Mean Square, RMS) 값을 비교하였다(Fig. 11과 12). 이를 통해 자력 이상의RMS 값이 클수록 역산으로 복원된 모델이 원판 이상체 모델과 유사하다는 것을 확인하였다. 즉 시추공 자력 이상값이 클수록 원판 이상체 모델에 더 근사하게 복원하는 것을 알 수 있었다. 총자력값은 이상체에서 유도된 유도자화를 지자기 방향에 정사영하여 계산하므로 시추공 위치의 방위각이 편각과 일치할 때 지자기의 반응이 가장 크고 그에 따라서 복원된 역산 모델의 유사도도 가장 크게 나타난다.

Fig. 11

Quantitative evaluation of borehole-only inversion results at a declination of 0° and an inclination of 51°.

이 연구에서는 지자기장의 복각은 ±51°로 가정하여 분석을 수행하였지만 대부분의 자력자료는 자극화 변환을 통해 복각에 의한 영향을 제거한 후 분석을 시도한다. 따라서, 자극화 변환된 시추공 자료의 역산 성능을 평가하였고, 자극화 변환한 자력 이상의 크기가 자극화 변환하지 않은 자력 반응보다 상대적으로 작은 것을 확인하였다(Fig. 12). Fig. 13에 나타난 바와 같이, 지자기장 편각이 0°인 조건에서 시추공 방위각이 이상체 기준으로 0°인 지점의 자극화 변환한 자료와 자극화 변환하지 않은 자료를 비교했을 때, 자극화 변환하지 않은 자료가 IoU, Dice 계수 및 정밀도가 높게 나타나고 재현율은 비교적 낮은 값을 나타냈다. 이는 자극화 변환하지 않은 자료의 역산 결과가 유사도가 더 크다는 것을 의미한다. 따라서 자극화 변환하지 않은 자력 자료가 자극화 변환 자료에 비해 자력 반응이 더 크고 이에 따라 역산 결과가 더 우수하게 나타나는 것을 확인하였다.

Fig. 12

Root mean square (RMS) of magnetic anomalies as a function of borehole location at a declination of 0° and an inclination of 51°. The last data point corresponds to the reduction-to-the-pole (RTP) case.

Fig. 13

Quantitative evaluation of the inversion performance for inclination angles of 51° and 90° (RTP) .

실제 자력 탐사 현장에서는 비용 및 작업 여건의 제약으로 인해 다수의 시추공을 확보하기 어려운 경우가 많다. 이러한 조건을 고려할 때, 시추공 위치에 따른 복합 역산 해상도 차이는 효율적인 탐사 설계 측면에서 중요한 시사점을 제공한다. 이 연구에서는 원판 이상체 모델에 의한 지표, 시추공 자력 반응을 계산하고, 복합 역산을 수행하여 결과를 정량적으로 분석하였다. 결과적으로 시추공의 위치와 지자기장의 편각 사이에 나타나는 관계를 파악할 수 있었고, 복합 역산의 해상도가 가장 우수한 최적의 시추공 위치는 북반구에서는 편각 방향으로 남반구에서는 편각의 반대 방향이다. 실제 자력 탐사에서는 시추를 진행하기에 앞서 육상, 항공 및 드론 자력 탐사를 통해 이상체의 수평적 위치를 추정할 수 있다. 이후 이 연구에서 제시한 결과를 바탕으로, 지자기장 편각과 이상체 기준 시추공의 방위각을 일치시킴으로써 제한된 시추공 조건에서도 보다 효율적인 탐사가 가능할 것으로 기대된다.