서 론

자연재료인 지반은 그 특성이 공간에 따라 변하는 불균질성(heterogeneity)을 내포하고 있어 설계 및 시공 시 불확실성(uncertainty)을 피할 수 없으며, 이로 인해 굴착 안정성 분석을 위한 지반물성의 산정에 어려움이 따른다. 지반물성 산정을 위한 실내실험의 경우, 시료의 대표성 문제와 작은 크기의 시료로부터 획득된 물성을 현장의 불연속면을 포함한 물성으로 변환 시 생기는 오차 등으로 인해 설계 시 예측된 지반거동과 시공 중 관측된 거동간의 차이가 발생할 수 있다. 실제 현장의 지반거동이 굴착 전 예측된 거동과 크게 상이한 경우에는 굴착설계의 적정성을 재검토할 필요가 있다. 설계 적정성 검토를 위해서는 현장 지반거동을 모사할 수 있는 재료물성이 재산정되어야 하며, 추가적인 지반조사 및 현장실험에 대한 대안으로서 역해석(inverse analysis) 기법이 활용될 수 있다.

지하굴착에 대한 초기의 역해석 연구는 굴착문제를 단순한 형태의 수학적 모델로 표현할 수 있는 재료의 선형탄성 거동(linear elastic behavior)을 가정한 탄성 역해석에 기반하였다(Jang and Kim, 1998; Jun et al., 1994; Kaiser et al., 1990; Kim and Jang, 1995; Kim et al., 2000; Kim, 2005, Lee and Kim, 1991a, 1991b; Park and Park, 2016; Sakurai and Takeuchi, 1983; Yang and Jeon, 2002). 이러한 탄성 재료거동 가정으로부터 추정된 암반의 변형계수는 탄성한계, 즉 항복강도(yield strength) 이후의 소성거동(plastic behavior)을 고려한 경우에 비해 그 값이 보수적으로 낮게 산정되므로 변형계수 기반의 굴착설계에 활용될 수 있다(Park and Park, 2017). 그러나 재료의 탄성거동 가정은 재료에 작용하는 응력이 작거나 응력변화가 크지 않은 경우에 타당한 조건으로서, 암반상태와 지하굴착 규모에 따라 굴착 후 소성변형이 일어날 수 있고, 이 경우에는 재료의 강도특성을 고려하여 굴착 안정성을 분석해야 하므로 암반의 강도정수 추정이 가능한 탄소성 역해석의 적용이 필요하다. 특히, 지반굴착 후 붕괴된 구간의 사고원인 조사와 붕괴 시 강도특성을 파악하기 위해서는 비선형 변형거동을 고려한 탄소성 역해석의 적용이 적합하다고 할 수 있다.

본 연구에서는 지하굴착 주변암반의 강도를 추정하기 위한 지반변위 기반의 역해석 기법을 제안하였다. 제안된 역해석 기법은 응답면 기법(response surface method)을 기반으로 한다. 유한차분 수치코드인 FLAC2D(Itasca, 2018)를 굴착해석 솔버로 사용하였으며, 응답면 기법과 수치해석코드를 연계하기 위한 시뮬레이션 툴을 개발하였다. Mohr-Coulomb 탄소성 재료모델에 기반한 지하굴착 예제를 통해 굴착주변 암반의 강도 매개변수(점착력과 내부마찰각)를 추정하여 계측변위 기반 역해석에 의한 암반강도 추정의 적용성 및 한계점을 분석하였다.

직접탐색법에 기반한 역해석

역해석 방법은 정해석 문제의 지배방정식을 역으로 수학적 정식화(mathematical formulation)하여 분석하는 역산법(inverse method)과 역해석 매개변수 추정을 최적화 문제로 정의하고 해당 목적함수가 최적값에 도달하는 매개변수 조합을 찾는 직접탐색법(direct search method)으로 분류된다. 역산법은 역해석 문제를 수학적으로 정식화하여 미지변수 추정을 위한 반복연산이 없어 계산시간이 상대적으로 짧은 장점이 있지만, 지반거동의 비선형성, 시간의존성, 점탄성거동 등을 고려하는 경우, 역 지배방정식의 정식화에 어려움이 있어 주로 수학적으로 단순한 문제(예, 탄성지반을 가정한 경우)에 한정되어 사용된다(Park and Park, 2016). 직접탐색법은 특정 물리적 특성에 대한 관측값과 수치해석적 계산값간의 오차를 목적함수로 정의한 후, 목적함수 값이 최소가 되도록 매개변수 값을 변경하면서 반복적으로 정해석을 수행하여 최적해를 찾는 방법이다. 직접탐색법은 역산법과는 달리 지반거동 지배방정식에 대한 수학적 정식화가 필요하지 않아 비선형 역해석 문제에도 적용할 수 있는 이점이 있다. 최근 컴퓨터 계산성능의 향상에 따라 반복연산에 대한 계산부담이 줄어들어 비선형 역해석 문제와 복잡한 해석모델에 적용이 쉽지 않은 역산법에 비해 직접탐색법의 활용빈도가 높은 편이다.

최적화 목적함수의 일차 또는 이차미분을 이용하는 도함수 기반 최적화(derivative-based optimization) 기법이 역해석 문제에서 많이 활용되어 왔으나, 이 기법은 최적해 추정을 위한 매개변수의 탐색방향(search direction)과 이동거리(step size) 결정조건에 따라 최적해 탐색이 실패할 수 있는 한계점이 있다. 즉, 탐색방향과 이동거리 결정이 적절하지 않은 경우, 최적점을 훨씬 벗어난 공간으로 매개변수를 이동시켜 최적화 과정이 수렴하지 않고 발산해 버릴 수 있다. 본 연구에서는 이러한 도함수 기반 최적화 기법의 한계점을 보완할 수 있는 응답면 기법을 최적해 탐색 알고리즘으로 적용하였다. 응답면 기법은 시스템의 응답이 특정 매개변수들에 의해 영향을 받는 경우, 이들 매개변수에 대한 응답면을 근사 모델링(approximate modeling)하여 시스템 거동을 해석하는 수학적 기법이다. 이 기법은 화학공정에서 반응이 최적(최대 또는 최소)이 되는 조건을 결정하기 위해 제안된 방법이다(Box and Wilson, 1951). 응답면 기법에서는 매개변수의 영역 내에서 추출된 몇 개의 점들에 대해 시스템 응답값들을 계산하고, 식 (1)과 같은 다항식 형태로 시스템 응답면을 추정한다. 추정된 응답면의 수학적 관계를 이용하여 매개변수의 탐색방향과 이동거리를 결정하므로 앞서 기술한 도함수 기반 최적화 기법의 문제점을 피할 수 있다.

     (1)

여기서 f(x)는 시스템의 응답, β는 회귀계수(regression coefficient), x는 설계 매개변수, k는 설계 매개변수의 개수, 𝜀은 오차항이다.

응답면 기법을 이용한 최적해 탐색은 식 (2)와 같이 시스템 응답의 관측값과 계산값 간의 오차를 목적함수로 정의하고, 이 함수의 값이 최소로 수렴하도록 역해석 매개변수 값을 조정하면서 반복적으로 정해석을 수행하여 이루어진다.

     (2)

여기서 F(x)는 최소자승잔차(least squares residual)를 나타내는 목적함수, x는 역해석으로부터 추정되어야 할 매개변수, N은 총 계측지점의 수, f_i(x)와 F_i는 각각 i번째 계측지점에 대해 해석적으로 예측된 값과 실제 관측된 값, ΔF_i는 계측오차와 관련된 인자이다.

시스템 응답면의 근사 정확도는 설계 매개변수의 관심영역의 크기에 의존하며, 관심영역 크기가 작을수록 정확도는 높아지는 경향을 보인다. 일반적으로 Fig. 1과 같이 초기에 작은 관심영역을 설정하여 응답면을 근사한 후 임시 최적해를 추정하고, 시스템 응답의 관측값과 계산값 간의 차이를 토대로 최적해 접근정도를 분석하여 점진적으로 관심영역을 변경시키면서 최적해를 탐색한다. 각 반복연산단계에서 관심영역은 식 (3)과 (4)에 의해 설정되며, 관심영역의 크기는 식 (4)와 같이 관측값과 계산값의 오차에 의해 결정된다.

     (3)

     (4)

여기서 x는 역해석 매개변수, µ는 관심영역의 보정계수, A는 수축모델 상수, 𝜀는 수축계수로서 해석적으로 예측된 값(f_i)과 실제 관측된 값(F_i)과의 평균 오차를 나타낸다.

Fig. 1.

Schematic of parameter identification process using the response surface method.

본 연구에서는 응답면 기법을 적용하여 지하굴착 주변암반의 강도특성을 파악하고자 하였다. 연속체 해석모델에 기반한 지하굴착 문제를 고려하였으며, 굴착해석을 위한 솔버로서 유한차분 수치코드인 FLAC2D를 사용하였다. 수치코드에 의한 굴착해석과 응답면 기법을 이용한 최적해 탐색 과정을 처리하기 위한 시뮬레이션 툴을 Fig. 2와 같이 작성하였다. 시뮬레이션 툴은 수치해석을 위한 소스파일과 최적화와 관련된 데이터파일의 입출력 처리, 굴착 해석 및 결과 출력, 응답면 생성 및 최적해 탐색 등 일련의 반복적인 매개변수 추정 과정을 자동적으로 처리한다. 본 연구에서는 Mohr-Coulomb 탄소성 모델의 강도 매개변수, 즉 점착력과 내부마찰각을 역해석 매개변수로 정의하였으며, 선형탄성 변형거동을 가정한 탄성 역해석의 매개변수인 재료의 변형계수와 초기응력(측압계수)은 알려진 값이라고 가정하여 강도 매개변수 추정을 위한 탄소성 역해석 시뮬레이션 툴을 작성하였다.

Fig. 2.

Simulation tool for inverse analysis.

암반의 강도정수 추정을 위한 지반변위 기반 역해석의 적용성 분석

지하굴착 예제

계측 지반변위 기반의 역해석을 이용한 암반강도 추정의 적용성을 분석하기 위해 Fig. 3과 같이 지하굴착 문제를 설정하였다. 암반의 거동은 Mohr-Coulomb 탄소성 재료모델을 따르는 것으로 가정하였으며, 적용한 암반물성은 Table 1과 같다. 수치모델의 해석영역은 가로 80 m, 세로 50 m이었고, 굴착 직경은 10 m, 지표에서 굴착 천단까지의 심도는 15 m로 설정하였다. 전단면 지반굴착을 고려하였고, 역해석을 위한 입력자료는 Fig. 3에 표시된 천단부 P1 지점과 측벽부 P2 지점의 굴착 후 지반변위를 이용하였다. 입력 지반변위를 산정하기 위해 Table 1의 물성을 사용하여 수치계산을 수행하였으며, 암반 측압계수는 1.0으로 가정하였다. 이 해석조건에 대한 계측지점 P1과 P2에서 굴착 후 지반변위는 각각 3.962 mm와 8.691 mm로 계산되었다. 본 예제는 암반강도 추정 탄소성 역해석의 유용성을 분석하기 위한 가상의 굴착 문제이므로 굴착 직후 변위계측이 수행되는 것으로 가정하였으며, 굴착 전후 지반 보강 및 지보는 고려하지 않았다. 또한 계측 과정에서 오차가 발생하지 않는 것으로 가정하여 위에서 계산된 변위값들(3.962 mm와 8.691 mm)을 역해석 입력자료로 그대로 사용하였다. 본 연구에서는 Mohr-Coulomb 탄소성 재료모델의 강도 매개변수인 점착력과 내부마찰각을 역해석에 의해 추정해야 할 매개변수로 설정하였으며, 따라서 입력 지반변위 산정을 위해 사용된 점착력 100 kPa과 내부마찰각 30°가 역해석에 의해 추정되어야 할 매개변수 값들이다.

Fig. 3.

Numerical model used for inverse analysis.

강도 매개변수 추정 결과

Fig. 4는 역해석에 의해 추정된 점착력(c, 단위 kPa)과 내부마찰각(φ, 단위 °) 값들을 나타낸다. 이 그림으로부터 입력 지반변위 산정을 위해 사용되고 역해석으로부터 추정되어야 할 강도 매개변수 조합 Si = (c, φ) = (100, 30) 외에도 다양한 조합의 매개변수들이 추정된 것을 알 수 있다. 추정된 매개변수 조합에 대한 계측지점 P1과 P2(Fig. 3)에서의 굴착 후 지반변위는 역해석을 위해 입력된 계측변위와 유사한 수준을 나타냈다. 예를 들어 Fig. 4에 표시된 매개변수 조합 S1과 S2는 각각 (c, φ) = (144.0, 20)과 (57.6, 40)이고, 이 조합들에 대해 굴착으로 인해 발생하는 지반변위는 Fig. 5와 같으며 최종 수렴변위가 입력 계측변위에 근접한 것을 알 수 있다. 한편, Fig. 4로부터 역해석에 의해 추정된 점착력과 내부마찰각이 반비례, 즉 음의 상관관계(negative correlation)를 보이는 것을 알 수 있다. 즉, 내부마찰각의 증가에 따라 점착력이 감소하는 경향을 보이는데, 이는 내부마찰각이 커지면 점착력이 작아지는 조합, 반대로 내부마찰각이 작아지면 점착력이 커지는 조합으로 소성 포텐셜(plastic potential)이 유사한 수준으로 유지되어 계측 지반변위량을 만족시키는 다양한 조합의 강도 매개변수 값들이 존재할 수 있음을 나타낸다.

Fig. 4.

Strength parameters identified by inverse analysis.

Fig. 5.

Ground displacements at points P1 and P2 for the strength parameter sets S1 and S2.

이러한 강도 매개변수들간의 상관관계는 Gioda and Maier(1980)와 Arai(1993)의 연구에서도 확인할 수 있다. Fig. 6은 Gioda and Maier(1980)의 연구결과를 나타내는 그림으로, Mohr-Coulomb 재료모델을 토대로 압력터널(pressure tunnel)에 대한 탄소성 역해석을 수행하여 추정된 강도 매개변수(점착력과 내부마찰각)를 나타낸다. 이 그림에서 경계 II를 벗어나 강도 매개변수 값이 커지는 영역에서는 탄성 지반변위만 발생하여 입력 굴착변위에 근접하지 못하고, 반대로 경계 III를 벗어나 강도 매개변수 값이 작은 영역에서는 작은 강도정수 값으로 인해 굴착변위가 크게 발생하여 계측변위에 접근하지 못하는 매개변수 공간들을 나타낸다. 경계 II와 III은 강도 매개변수가 존재할 수 있는 매개변수 공간을 제한하는 상한과 하한의 경계로서 점착력과 내부마찰각이 음의 상관관계를 보임을 알 수 있다. 참고로, 경계 III의 첨자 p와 u는 각각 압력터널 내부에 압력이 작용하는 조건(loading)과 작용하지 않는 조건(unloading)을 나타내고, 경계 I의 내부는 공학적 또는 경험적으로 의미가 있는 강도 매개변수의 사전정보(prior information)에 해당하는 영역을 나타낸다. Arai(1993)는 Fig. 7(a)와 같이 수치모델을 설정하여 지반 지지력(bearing capacity) 분석을 위한 평판재하시험(plate loading test)을 대상으로 역해석을 수행하였다. 역해석으로부터 추정된 Mohr-Coulomb 재료모델의 강도 매개변수는 Fig. 7(b)와 같으며, 본 연구의 해석결과와 유사하게 점착력과 내부마찰각이 반비례적 관계를 나타냄을 알 수 있다.

Fig. 6.

Inverse analysis for the pressure tunnel (Gioda and Maier, 1980).

Fig. 7.

Inverse analysis for the conventional plate loading test (Arai, 1993).

암반 변형계수가 강도 매개변수 추정에 미치는 영향

본 연구에서는 암반 변형계수가 강도 매개변수의 추정결과에 미치는 영향을 분석하기 위해 입력 지반변위 산정 시 사용된 변형계수(800 MPa) 이외의 값들을 적용하여 강도 매개변수를 추정하였다. 추가로 고려한 암반의 변형계수 값들은 700, 1000, 1200 MPa이었다. Fig. 8은 역해석에 의해 추정된 변형계수별 강도 매개변수 값들을 나타낸다. 앞서의 변형계수 800 MPa 경우와 유사하게 다양한 조합의 매개변수들이 추정된 것을 알 수 있다. Fig. 8의 변형계수 1200 MPa에 대한 매개변수 값들 중 점선으로 표시된 부분은 입력 지반변위에 수렴하지 않은 최적해가 아닌 매개변수 조건을 나타낸다. Fig. 8에 표시된 매개변수 조합 S3, S4, S5, S6에 대한 계측지점 P1과 P2(Fig. 3)에서의 굴착 후 지반변위는 Fig. 9와 같다. Fig. 9로부터 최적해에 해당하는 매개변수 조합 S3, S4, S5의 경우는 역해석을 위해 입력된 계측변위에 수렴하는 경향을 보인 반면, 최적해가 아닌 S6의 경우는 계측지점 P2에서 입력 계측변위에 근접하지 못하는 것을 알 수 있다. 변형계수의 변화에 따라 추정된 점착력과 내부마찰각은 앞절의 역해석과 유사하게 반비례적 관계를 나타냈다. 또한 변형계수가 커짐에 따라 입력 계측변위를 만족시키는 강도 매개변수의 최적값이 작아지는 경향을 나타냈다. 이는 계측변위를 만족시키기 위해서는 변형계수가 커짐에 따라 재료의 소성변형량이 증가해야 하기 때문이다(Fig. 10). Fig. 10은 변형계수 조건은 다르지만 최종 발생변위는 동일한 응력-변형률 관계들을 나타낸다. 이 그림에서 볼 수 있듯이 서로 다른 변형계수 조건(변형계수 크기: E₃ > E₂ > E₁)에서 측정된 변형수준 S_m만큼 최종 재료변형이 발생하기 위해서는 변형계수의 증가에 따라 소성변형량이 늘어나야 한다. 소성변형의 증가는 강도 매개변수의 감소에 의해 발생하므로 Fig. 8과 같이 변형계수의 증가에 따라 추정된 최적 강도 매개변수의 값은 작아지는 경향을 보이게 된다.

Fig. 8.

Strength parameters identified under different elasticity conditions.

Fig. 9.

Ground displacements at points P1 and P2 for different strength parameter sets.

Fig. 10.

Stress-strain relationships under different elasticity conditions, which show the same total deformation level at a given stress.

역해석 결과 분석 및 토의

본 연구의 분석결과로부터 지반변위 기반의 탄소성 역해석에 의한 강도 매개변수 추정 시 암반의 변형계수가 알려진 값으로 설정되더라도 다양한 조합의 매개변수값들이 추정될 수 있고, 설정한 암반 변형계수에 오차(또는 시간에 따른 변화)가 있는 경우 매개변수 추정의 다양성이 확대될 수 있는 것으로 파악되었다. 즉, 암반 변형계수의 적용조건에 따라 재료 강도특성의 유일해(unique solution)를 찾지 못하고 다양한 조합의 강도 매개변수가 추정될 수 있는 것으로 조사되었다. 유일해 추정을 위해서는 재료의 변형특성에 대한 사전정보와 강도 매개변수간 상관관계에 대한 자료가 추가적으로 필요할 것으로 판단된다.

본 연구의 역해석에 사용한 지반변위는 굴착면(excavation surface)에서의 계측자료를 기반으로 한다. Fig. 4와 Fig. 8에 표시된 강도 매개변수 조합 Si와 S1-S5는 역해석에 입력된 지반변위와 유사한 수준의 지반변형량을 제공하는 매개변수 조건들이다. 입력 계측변위를 만족시키는 이 결과는 암반 굴착면(excavation surface)에서의 계측자료에 기반한 것으로서, 굴착면 이외의 지중변형에 대해서도 유사한 결과가 나타나는지 검토할 필요가 있다. 강도 매개변수 조건에 따른 지중변형 거동을 비교 분석하기 위해 굴착 천단부 상부측(계측선 A-A' 길이 15 m)과 측벽부 우측(계측선 B-B' 길이 20 m)에 Fig. 11과 같이 계측지점을 설정하였다. Fig. 12는 강도 매개변수 조합 Si와 S1-S5에 대해 측정된 지중변형 자료를 나타낸다. 천단부 상부측 계측선에서는 수직 지중변위를 측정하였고, 측벽부 우측 계측선에서는 수평 지중변위를 측정하였다. Fig. 12의 좌측에 음영으로 표시된 부분을 보면, 측정위치의 거리가 0인 굴착면에서는 앞서의 역해석 결과에서 알 수 있듯이 6가지 매개변수 조합에서 유사한 수준의 지반변형이 발생한 것을 알 수 있다. 그러나 굴착면으로부터 멀어짐에 따라 역해석 입력자료로 활용된 Si 조합의 변형거동과는 다른 형태를 보이는 것을 알 수 있다. 즉, 강도 매개변수 조건에 따라 굴착면에서의 변형거동은 유사하지만 굴착면 이외의 구역에서는 변형거동이 상이할 수 있는 것으로 파악되었으며, 이로부터 지중변형거동을 역해석 입력자료로 활용하면 다양한 매개변수 추정의 가능성을 피하거나 줄일 수 있는 매개변수 정제(parameter refinement)가 가능할 것으로 판단된다.

Fig. 11.

Locations for measuring ground displacements.

Fig. 12.

Ground displacements measured for different strength parameter sets.

한편, 굴착 후 계측기 설치 이전의 미계측자료 및 계측 오차 등을 고려한 정확한 입력 계측변위 산정이 어렵고, 탄소성 역해석에 의한 강도정수 추정 시 암반의 변형계수 및 측압계수가 영향을 미칠 수 있기 때문에 입력변위, 변형계수, 측압계수의 불확실성을 고려한 강도정수 추정이 합리적일 것으로 판단된다. 즉, Park and Park (2107)이 제안한 것처럼, 불확실성 개념을 적용하여 암반물성 매개변수의 단일해를 추정하지 않고 해당 매개변수가 존재할 수 있는 범위를 추정하는 것이다. 기존의 신뢰도 기반 안정성 평가는 설계단계에서 산정된 입력물성의 불확실성에 기반하므로 불확실성 정보가 충분하지 않고 대표성이 결여될 수 있는 한계가 있는데, 역해석을 통해 파악된 현장의 암반물성 변동성 자료와 조합하여 불확실성 정보를 업데이트함으로써 개선된 신뢰도 분석결과를 도출할 수 있을 것으로 판단된다(Park and Park, 2017).

결 론

본 연구에서는 최적해 탐색 알고리즘인 응답면 기법과 연속체 수치해석코드인 FLAC2D를 연계하여 지하굴착 주변암반의 강도 매개변수를 추정할 수 있는 계측변위 기반의 역해석 시뮬레이션 툴을 개발하였다. 시뮬레이션 툴을 Mohr-Coulomb 탄소성 지반거동 모델에 기반한 지하굴착 예제에 적용하여 주변암반의 강도 매개변수를 추정하였다. 분석결과에 따르면, 계측변위 기반의 탄소성 역해석은 기본적으로 재료 강도특성의 유일해를 제공하지 못하고, 다양한 조합의 강도 매개변수가 추정될 수 있는 것으로 파악되었다. 유일해 추정을 위해서는 변형특성에 대한 사전정보, 강도 매개변수간 상관관계와 같은 추가적인 정보가 필요하거나, 굴착면 이외의 지중계측자료를 활용한 매개변수 정제 기술의 적용이 필요한 것으로 파악되어 향후 이와 관련된 연구가 수행되어야 할 것으로 판단된다.

일반적으로 계측기 설치작업으로 인해 계측개시 시점과 지반굴착 시점간의 시간 차이가 발생하므로 계측실시 이전의 암반거동을 파악하는 것은 용이하지 않다. 이러한 미계측 자료는 역해석에 의한 지반물성 추정시 불확실성 요인으로 작용하게 된다. 따라서 계측개시 이전의 지반거동(미계측 자료)을 합리적으로 추정할 수 있는 기법에 대한 연구가 수행되어야 할 것이다. 또한 지반거동 측정을 위한 계측기 자체의 절대 측정오차와 서로 다른 계측기 사용으로 인한 상대 측정오차가 발생할 수 있으므로 실제 현장사례에 대해 역해석을 수행하는 경우 이를 고려할 필요가 있다. 한편, 실제 현장의 굴착변위는 3차원적으로 발생하므로 굴착 종방향 지반거동을 가정한 2차원 수치해석보다는 3차원 수치해석에 기반한 역해석 연구가 향후 수행될 필요가 있을 것으로 판단된다.