1. 서 론
일반적으로 철도 궤도는 자갈도상 궤도와 콘크리트도상 궤도로 구분되지만, 초기 부설 비용이 저렴하고 적절한 탄성 유지 및 관리보수가 쉬운 자갈도상 궤도 구조가 전세계적으로 오랫동안 사용되고 있다(Shin et al., 2017; Lee et al., 2020). 국내 철도 궤도의 경우, 약 80% 이상이 자갈궤도로 이루어져 있으며, 갈수록 증가하는 철도 사용으로 인해 자갈도상재료의 효율적인 유지보수기술에 대한 중요성 및 관심도가 증대되고 있다. 궤도의 자갈도상재료는 열차의 반복적 진동하중과 보수작업으로 인한 파쇄 및 마모, 외부물질의 침투, 노반재료의 상승 등의 원인으로 Fig. 1과 같이 지반 내부의 파울링(Fouling)층이 형성된다(Selig and Waters, 1994). 이러한 파울링층은 자갈입자 간 마찰효과(Interlocking effect)를 저감시켜 궤도의 안전 성능 악화를 유발한다(Leng and Al-Qadi, 2010). 현장에서는 철도지반 내부에 형성된 파울링층을 검측하기 위해 굴착 작업을 주로 진행해왔지만, 최근에는 지반투과레이더(Ground Penetrating Radar, GPR) 장비를 이용하여 도상 내부의 파울링 양상을 검측하는 시도가 이루어지고 있으며, 이와 관련된 실내 및 현장시험 연구들이 수행되었다(Shin et al., 2017; Lee et al., 2020; Benedetto et al., 2017; Robert et al., 2009; Sadeghi et al., 2018). GPR 장비는 전자기펄스를 지반 내부에 송출한 후 지반 내부의 변이양상에 대한 반응 신호를 고분해능으로 영상화할 수 있는 장비로서 현장에서 지반 변이 양상을 확인하는데 효과적이며, 단시간에 1차적인 이상체 검측이 가능하여 다양한 지반탐사분야에 활용되고 있다(Choi et al., 2001; Hong et al., 2015). 하지만, GPR 장비의 안테나 종류, 지반구조, 재료물성, 함수비 등에 따라 신호의 전파거동이 매우 복잡한 양상을 나타내기 때문에 관측된 신호 자료에 대한 면밀한 신뢰성 검토 및 검증작업이 요구되며, 이를 위해서는 복잡한 지반조건이 혼재된 상태에서도 내부 상태를 비교적 정확히 평가·검증할 수 있는 GPR 모델링 해석기술이 필요한 실정이다.
국외에서는 파울링층 식별을 위한 GPR 수치 모델링 연구가 시도된 바 있다. 일례로, Huang and Tutumluer(2011)는 이미지 지원 이산요소모델(image-aided based DEM) 기법을 통해 실험실에서 관측된 이미지를 바탕으로 오염된 골재 모델을 제안한 후 실내시험 및 DEM 시뮬레이션 결과를 비교하여 모델 매개변수를 보정하고 검증하였다. 또한, Ngo et al.(2016)은 DEM을 사용하여 오염된 파울링층의 미세한 거동을 모델링하였고, 파울링 수준이 증가함에 따라 자갈입자의 최대 접촉력이 감소한다는 사실을 확인하였다. 국내에서는 GPR 수치모델링 기술 개발 및 응용 연구가 점차적으로 활발해지고 있지만 아직까지는 지하공동(Cavity) 탐지 분야쪽으로 적용된 사례가 주를 이루고 있다. Cheon(2016)은 도심지 도로 천부의 지하공동 탐지를 위한 GPR 해석연구를 수행한 후 그 신호반응을 실대형 실험결과와 비교분석한 바 있으며, Jang et al.(2016)은 GPR 장비의 전자기펄스 방사를 모사하는 3차원 모델을 개발하여 단순화된 공동모델에 대한 GPR 반응을 검토하였고, 개발된 모델이 도심지 도로하부 지하공동의 부존 정보를 얻는데 효과적으로 활용될 수 있음을 확인하였다. 상기 언급된 바와 같이 국내에서 GPR 수치모델링 기술을 지하공동 탐지에 적용한 사례는 있지만 아직까지 철도궤도 자갈도상 내부의 파울링층을 검측하기 위한 목적으로 GPR 수치모델이 활용된 사례는 전무한 실정이다. 따라서 본 연구에서는 python 언어 기반의 오픈소스코드인 gprMax 소프트웨어와 RSA(Random sequential Absorption) 알고리즘을 이용하여 자갈도상 재료의 파울링 정도에 따른 GPR 신호특성을 정량적으로 분석하고 평가할 수 있는 수치해석 기법을 제안하였다. 또한, 해석 결과의 예측 신뢰성을 검증하고자 모형시험체에 대한 GPR 신호 특성을 모형실험 데이터와 비교·검토하였고, 다양한 지반조건 별 시험과 해석에서의 파울링층의 식별 여부를 확인하였다.
2. 수치모델링
본 연구에서 사용된 gprMax는 GPR 장비에서 송출, 수신되는 신호의 전파 거동을 모사할 수 있는 해석 소프트웨어로서(http://www.gprmax.com) Maxwell 방정식을 구성방정식으로 사용하며, 방정식의 해를 근사화하기 위해 공간 및 시간의 연속성을 이산화(discretization)한 시간영역유한차분법(Finite-Difference Time-Domain) 기법을 도입하고 있다(Yee, 1996).
여기서, E는 전기장 강도(V/m), H는 자기장 강도(A/m), D는 전속밀도(C/m2), B는 자속밀도(T), JC는 안테나 요소에서의 전류 밀도(A/m2), JS는 저항 소자에 흐르는 전류 밀도(A/m2), t는 시(sec), qv는 체적전하밀도(C/m3)를 의미한다.
gprMax는 기본적으로 시간영역유한차분법(Finite Difference Time Domain, FDTD)을 통해 3차원 공간 상의 전가지장 해를 구하며 이를 TEM(transverse electromagnetic) 모드라고 칭한다. 하지만 본 연구에서는 복잡 다양한 입자형상 및 입도분포를 가지는 자갈도상에 대한 해석의 효율화를 위해 2차원 TM(transverse magnetic) z-direction 모드를 이용한 cell을 사용하였다(Yee, 1996). 이 경우 Maxwell 방정식은 다음과 같이 3개의 상호 연계된 편미분방정식 형태로 표현될 수 있으며, 시간과 공간에 대해 이산화되어 Fig. 2에 표시된 각각의 FDTD 셀에 적용된다.
Fig. 2
Single FDTD Yee cell showing electric (red), magnetic (green), and zeroed out (grey) field components for 2D transverse magnetic (TM) z-direction mode (Yee, 1966)
여기서, ɛ는 유전율(farads/meter), μ는 투자율(henries/meter), σ는 전기전도도(Siemens/meter), σ*는 자기손실(Ohms/meter)를 의미한다.
한편, 모델에서 이산화된 요소의 크기인 Δx, Δy, Δz 및 미소시간 Δt의 크기를 적절하게 정의하는 것이 매우 중요하다. 요소의 크기 및 미소시간간격이 작을수록 실제 신호의 전파거동과 더 유사해지지만 해석을 수행하는 전산장비의 제한된 메모리 공간과 유한한 처리속도를 감안하기 위해서 적정 수준의 크기 및 간격을 정의하는 것이 요구된다. 해석에서 요구되는 요소(cell) 크기는 파장 길이(λ)의 1/10보다 작아야하며 파장길이 λ는 식 (3)을 통해 산정된다.
여기서, c는 빛의 속도 2.99×108m/s이며, f는 신호의 최대주파수(Hz), εr은 지반재료의 상대유전율을 의미한다. 요소크기가 정의되면 최소 미소시간 간격인 Δt는 다음의 Courant 안정조건을 만족해야한다(Taflove, 1995; Sullivan, 2000). 본 연구에서 적용된 TMz 모드의 경우 ∆z→∞로 간주한다.
2.2 입자 모델링 및 재료 물성
본 연구에서는 Fig. 3과 같이 파울링층의 오염 수준을 3단계로 구분하여 자갈도상 지반을 모델링하였고, 자갈층과 노반의 경계에는 ‘*Ballast Pocket’를 모사하여 실제 지반과 유사한 형태로 해석도메인을 구축하였다(*Ballast Pocket은 자갈도상이 노반을 침투하여 형성된 굴곡을 말하며, 지표에서 침투된 수분이 모이게 되는 경향이 있음). Fig. 4는 파울링층의 오염수준을 보여주는데, Clean Ballast(이하 CB), Fully Fouled Ballast(이하 FFB), Floated Ballast(이하 FB)로 자갈도상 재료의 파울링 단계를 구분하고 있다. CB는 세립분이 없이 자갈입자의 공극이 공기로 채워져 있는 깨끗한 상태를, FFB는 공극안에 세립분이 차있는 상태를, FB는 세립분이 자갈사이의 공극 부피 이상으로 포함되어 자갈입자가 세립분 사이에 떠있는 것과 같은 매우 오염도가 심각한 상태를 의미한다.
또한, 각 지반재료별 입도분포를 산정하기 위해 사전에 체분석시험(Sieve analysis)을 수행하였고, 각 지반재료별 입도분포곡선(Particle size distribution)은 Fig. 5와 같다. 본 연구에서는 각 재료별 입도분포가 합리적으로 모사되도록 자갈 입자들을 생성하여 모델의 입력인자로 적용하였다. 여기에는 자체적으로 개선된 RSA(Random sequential Absorption) 알고리즘이 적용되었는데 평균적으로 반환할 입자의 개수 및 입자 직경의 범위를 고려하여 입자를 무작위로 생성하되, 입자들이 서로 중첩될 경우 전체 도메인이 아닌 중첩된 입자 인근 영역에서 입자를 재생성하도록 하여 전체 입자생성 시간을 상당부분 단축시킬 수 있는 효과를 구현하였다(Fig. 6). 한편, 본 연구에서 사용된 지반재료의 전자기 물성은 Table 1과 같다. 이 중 각 재료의 상대 유전율은 실대형 모형 시험에서 계측된 등가의 전자기 물성을 VMF(Volumetric Mixing Formular)에 적용한 후, 개별 요소의 물성 값을 역산(Back-calculation)하여 산정하였다(Fig. 6 참조).
Table 1.
Material properties used in the study
2.3 GPR 안테나 제원
본 해석연구에서는 총 2가지 종류의 GPR 장비(1,600MHz, 400MHz)를 고려하였고, 파형의 종류는 해당 GPR 안테나에서 일반적으로 사용되는 ‘Gaussian waveform’을 사용하였다. 각 안테나의 사양은 Table 2에 제시되었다. 소스와 관련된 중요한 입력 정보 중 하나는 바로 중심주파수(centre frequency, Hz)이며, 식 (5)와 같이 장비의 High pass filter 주파수(f1)와 Low pass filter 주파수(f2) 사이의 관계를 기반으로 산정하거나, Taguchi method(Taguchi et al., 2005; Warren and Giannopoulos, 2011)과 같은 최적화 알고리즘을 통해서도 산정될 수도 있다. 본 연구에서는 Taguchi method를 통해 중심주파수를 산정하였다.
Table 2.
Specification of GPR Antenna
| Type | GSSI 1.6 GHz | GSSI 400 MHz | Unit |
| Range | 12 | 50 | ns |
| High pass filter, f1 | 250 | 100 | MHz |
| Low pass filter, f2 | 3000 | 800 | MHz |
| Samples per scan | 512 | 512 | - |
| Scans per second | 100 | 120 | - |
이러한 소스 정보는 모델의 입력자료(input parameter)로 사용된다. TMz 모드에서 GPR 장비의 안테나의 물리적인 구조를 다 형상화하지 않고 송신부와 수신부가 각각 하나의 점(point)으로 정의된다. 본 연구에서 사용된 실대형 모의시험 장비는 Ground-coupled 방식이기 때문에 해석 모델에서 안테나 위치는 지표면과 완전히 맞닿은 것으로 간주하였고 송신부와 수신부 사이의 이격거리는 0.15m로 설정하였다.
3. 자갈도상 및 파울링층 모형시험체
본 연구에서는 제안된 해석모델의 예측 신뢰성을 검증하고자 자갈도상에서 발생할 수 있는 파울링층의 오염 수준을 3단계로 구분하여 모형시험체를 제작하였다. Fig. 7은 모형시험제 조성 순서를 보여준다. 모형시험체에서 CB층은 실제 철도현장에 투입되는 자갈도상재료를 공급처에서 구입하여 세립분을 전혀 추가하지 않은 상태로 조성하였고, FFB층은 자갈도상재료와 세립분의 부피피를 1:1로, FB층은 자갈도상재료와 세립분의 부피피를 1:2로 혼합하여 모형시험체를 조성하였다. 또한, Fig. 3의 ‘Ballast Pocket’을 모사하기 위해 노반(Subgrade)면에 깊이 300mm, 폭 900mm의 홈을 만들어 CB, FFB, FB 재료로 다짐을 하고 각각의 섹션에 CB, FFB, FB 시료로 약 200mm 높이까지 다짐작업을 하였다. 침목 두께에 해당되는 200mm 지층은 CB로 조성하여 현장의 자갈궤도 자갈도상층을 모사하였다. GPR 장비의 탐사 측선은 모형체의 횡방향으로 선정하였고, 안테나의 주파수 별로 동일 위치를 관측하여 주파수 변화에 따른 신호 변화를 비교 분석할 수 있도록 하였다. 모형시험체 시험에서 사용된 메인 컨트롤러는 미국 GSSI(Geophysical Survey Systems, Inc.)사에서 제작된 4채널 SIR-30 System이며, GPR 안테나는 1,600MHz 및 400MHz 장비를 사용하였다(Lee et al., 2020).
4. 결과 및 분석
본 연구에서 제안된 수치해석모델을 이용하여 1회당 총 22ns의 해석시간을 고려하여 전자기 해석을 수행하였다. 신호분석에서, 측선의 깊이방향의 전자기파 진폭을 1차원적으로 나타낸 것을 A-scan이라 칭하며, 탐사 측선을 따라 깊이 방향의 전자기파 신호를 2차원으로 나타낸 것을 GPR 신호의 B-scan이라고 표현한다. B-scan 산정 시 안테나의 이동 간격은 0.005m, 반복 연산회수는 1,200회로 정의하였다. Fig. 8은 시간 경과에 따른 GPR 신호의 전파 양상을 보여주는 snapshot 이미지이다. gprMax로부터 얻어진 출력자료는 Paraview(https://www.paraview.org)와 같은 후처리 소프트웨어를 통해 시각화(visualization)할 수 있다. 이러한 snapshot 이미지는 GPR 전파 신호가 서로 다른 매질(지층)을 만날 때의 반응 신호 특성을 시각적으로 확인하는데에 도움이 되지만 이것을 정량적으로 분석하고 평가하는 데에는 한계가 있다. 따라서, 실무에서는 특정위치에서 시간에 따른 반응 신호 변화를 관측할 수 있는 A-scan이 더 유용한 분석 자료가 될 수 있다.
Fig. 9는 본 연구에서 고려된 모형 시험체 조건(Fig. 3)에 따른 도메인 중심부에서 산정된 A-scan 결과를 나타내며, 이를 모형체 시험 결과와 비교한 결과를 보여준다. 각 시험 조건에서 해석과 실험의 A-scan 결과에 약간의 오차가 있음을 알 수 있는데 이는 본 연구에서 고려된 2차원 수치해석모델이 3차원 형상을 지닌 안테나의 송수신부를 각각 하나의 점(point)으로 간주하였기 때문에 발생한 것으로 사료된다. 보다 신뢰성 있는 A-scan 자료를 확보를 위해서는 향후 안테나의 형상을 온전히 고려하는 3차원 해석모델로의 확장이 필요할 것으로 판단된다.
Fig. 10～12은 400MHz 저주파 장비 사용 시의 각 시험조건 별 도출된 B-scan 이미지를 보여준다. B-scan 도출 시, 신호의 세기가 깊이방향으로 감쇠한다는 가정아래, 관측데이터와 해석데이터 모두 동일한 비율의 감쇠특성보정(Gain compensation)이 적용되었다. Fig. 10～12에서 수치해석으로부터 얻어진 B-scan은 고려된 모든 파울링 상태의 자갈층과 노반의 경계면을 비교적 명확하게 식별하였지만 CB층과 FFB층 사이 경계면은 해석결과에서 분명한 반사파 신호로 식별되지 않음을 볼 수 있다(Fig. 11). 이는 FFB층이 공극 내 세립분 함유량이 작아 상부 CB층과의 재료적 특이성이 명확히 구분되기 어려웠기 때문으로 분석된다. 그럼에도 수치해석모델은 CB층과 FB층 사이 경계면은 비교적 선명하게 식별하였는데(Fig. 12), 400MHz 장비에서 송출되는 전자기파가 자갈도상의 입도분포, 공극 내에 존재하는 세립분의 전자기적 특성을 식별하기에 충분히 높은 진폭(Amplitude)을 가지고 있기 때문인 것으로 사료된다. Fig. 13(a)는 CB-FB-Subgrade 조건과 CB-FFB-Subgrade 조건에서의 시간 도메인 상 A-scan(x=4m)을 비교해서 보여주고 있다. 구간 t1은 안테나 내 직접파(direct wave)에 의한 반사파 신호를 의미하며 t2 이후부터가 지반매질에 의한 반사파 신호이다. Fig. 13(b)는 Fig. 13(a)의 t2 이후 영역에서의 반사파 신호를 시간 축이 아닌 깊이 축으로 변환하여 그래프를 도출한 결과이며, 자갈층에서 파울링층으로 매질이 바뀌는 경계에서 FB층이 FFB층보다 상대적으로 높은 진폭 에너지 값을 나타내는 것을 확인하였다. 즉, 자갈층과 FB층 경계에서 반사파 강도 및 변동성이 크게 나타나면서 B-scan 상에서 파울링층 경계가 두드러지게 나타난 것으로 판단된다. 특히, 모형체 시험의 B-scan 자료에서 다량의 노이즈로 인해 다소 불분명하게 나타났던 파울링층의 경계가 해석데이터에서는 보다 명확하고 선명하게 식별된다는 점에서 향후 제안된 수치해석결과를 참고자료로 활용하여 계측자료를 분석, 검토한다면 복잡한 지반조건 및 노이즈가 혼재된 GPR 신호 자료 분석 시 경계층 판독의 어려움이 상당부분 해소될 수 있을 것이라 판단된다.
한편 Fig. 14는 1,600MHz 고주파 장비를 사용하였을 때의 시험 및 해석으로부터 도출된 B-scan 이미지를 보여주고 있다. 고주파 신호의 경우, 실험과 해석 모두 스캔 B-scan 이미지의 해상도는 높아지는 반면 재료 경계면의 식별은 어려웠다. 이는 경계면에서 반사되는 신호의 반사파의 상대적 차이가 크지 않기 때문이며, 이를 보완하기 위해 향후 필터링, 증폭 등 최적의 신호처리 기법의 개발 등이 필요할 것이다. 그러나 고주파 신호의 경우 높은 해상도로 자갈도상 재료의 입자경계, 공극 등의 특성 변화를 반영한 산란 특성을 활용하여 분석 기법 개발의 데이터로 유용하게 활용될 수 있을 것으로 사료된다.
5. 결 론
본 연구에서는 gprMax와 RSA 알고리즘을 이용하여 철도자갈도상 재료의 파울링 정도에 따른 GPR 신호특성을 정량적으로 분석하고 평가할 수 있는 해석 기법을 제안하였다. 해석 결과의 예측 신뢰성을 검증하고자 실대형 모형 시험체에 대한 GPR 신호 특성을 실험자료와 비교·검토하였고, 각 시험 조건 별 파울링층의 식별 여부를 확인하였다. 이를 통해 도출된 결론은 다음과 같다.
(1) 각 시험 조건에서의 A-scan은 해석과 모형시험 간 약간의 오차가 존재하였는데 이는 2차원 수치해석모델 지닌 근본적인 한계(3차원 형상을 지닌 안테나를 두 개의 점으로 가정) 때문인 것으로 사료된다. 해석결과의 신뢰성을 높이기 위해서는 안테나의 실제 형상을 반영할 수 있는 3차원 모델로의 확장이 요구된다.
(2) 자갈도상층과 노반 경계층 또는 파울링층과의 경계면에 대한 반사파 이미지(B-scan)를 이용한 분석을 위해서는 상대적으로 저주파인 400MHz의 신호가 용이한 특성을 보였다. 이는 자갈도상의 상태가 깨끗할수록, 파울링이 심하게 진행된 FB층에서 반사파 강도가 크게 발생되는 것을 확인할 수 있었다.
(3) 1,600MHz 고주파 장비 사용 시 스캔자료의 해상도는 높아지는 반면 파울링층 및 암반층 경계의 식별은 불가능하였다. 이는 경계면에서 반사되는 신호의 반사파의 상대적 차이가 크지 않기 때문이며, 이를 보완하기 위해 향후 필터링, 증폭 등 최적의 신호처리 기법의 개발 등이 필요하다.
(4) 1,600MHz 고주파 신호의 경우, 실험과 해석 모두 스캔 B-scan 이미지의 해상도는 높아지는 반면 재료 경계면의 식별은 어려웠다. 이는 경계면에서 반사되는 신호의 반사파의 상대적 차이가 크지 않기 때문으로, 이를 위해서는 신호처리기법의 추가적인 연구가 필요할 것으로 판단된다. 반면 고주파 신호의 경우 높은 해상도로 자갈도상 재료의 입자경계, 공극 등의 특성 변화를 반영한 산란 특성을 활용하여 분석 기법 개발의 데이터로 유용하게 활용될 수 있을 것으로 기대된다.
