In this paper, we propose an threshold-based Schur algorithm for estimating the media characteristics of sub-bottom multi-layers by using the signal generated by a parametric array transducer. We use the KZK model to generate a parametric array signal, and use the proposed threshold-based Schur algorithm for estimating the reflection coefficients of multiple sea bottom layers. Using computer simulation, we verify that the difference frequency component generated by the KZK model prevails over the signals of primary frequencies at long range. For the simulation, we use the transmit signal generated by the KZK and the reflected signal obtained from a lattice filter model for the seawater and sub-bottom of multi-level non-homogeneous layers. Through the simulation, we verify that the proposed threshold-based Schur algorithm can give much more accurate and efficient estimates of the reflection coefficients than methods using received signal, matched filter output signal, and normal Schur algorithm output. 지층 탐사기술은 지하자원 탐사, 지층 매질 분석, 은닉 목표물 검출 등에 많이 응용되는 기술이다. 해저지층탐사 기술의 핵심은 알려진 송신된 음파 신호와 해수층과 해저지층의 퇴적물의 다양한 배치 및 형태에 따라 반사되어온 수신된 음파 신호를 이용하여 해저지층의 특징을 추출하는 것이다( Lurton, 2002 ). 해저지층탐사를 위하여 일반적으로 사용되는 송신 신호는 Burst, Ricker pulse, Chirp 신호등이 사용된다. 이중 고해상도의 신호를 얻기 위한 Chirp sonar는 LeBlanc과 Mayer 연구진에 의해 Chirp sonar 시스템이 개발된 후 해저 지형의 표층 및 매설물 탐지, 퇴적물의 특성 분석 등 다양한 분야에서 응용 연구되고 있다( Tolsma, 2004 ). 1990년 후반부터 연구되어 온 파라메트릭 배열 sonar 기술은 해양 통신 분야 및 탐사분야에 적용되어 왔다. 그러나 국내에서는 파라메트릭 배열 sonar를 이용한 해저 지층 탐사 기술에 관한 전반적인 소개 및 성능 분석 연구가 많이 진행되어 있지 않을 뿐만 아니라 해저 지층 분해능 향상에 관한 연구는 거의 없는 실정이다. 해저면은 다양한 밀도, 속도, 감쇠 계수, 반사 계수 등을 갖는 다양한 형태의 퇴적물과 암반 등으로 이루어져 있으며 퇴적층에 따른 물성차이는 해저 지층 탐사장비에 각기 다른 산란 신호로 측정된다. 수신된 산란 신호로부터 퇴적물의 물성에 따른 임피던스, 감쇠 계수, 반사 계수 등을 추정할 수 있다. 이러한 지질음향 역산란 기법은 음향학, 지질음향학 그리고 지구물리 탐사 분야에서 확률적 접근, 최적화, 패턴인식, 특이값 분해 (Singular value decomposition, SVD), Layer peeling 등 다양한 연구방법으로 연구되고 있다 ( Lurton, 2002 ; Tolsma, 2004 ; Rakotonarivo et al., 2011 ; Schwetlick, 1983 ). 공진 주파수로 변조된 수신 신호에 포락선 검파를 하는 일반적인 지층의 분석 방법은 반사된 sonar 신호의 펄스폭에 따른 거리 해상도의 한계가 존재한다( Lurton, 2002 ; Tolsma, 2004 ). 이러한 거리해상도를 극복하기 위해 layer peeling 방법인 Schur 알고리즘을 사용하여 샘플링 시간의 해상도로 지층의 변화 특성을 계산할 수 있는 방법 등이 있다 ( Bruckstein and Kailath, 1987 ; Kailath, 1987 ). Schur 알고리즘을 이용하여 해수층과 해저지층의 이론적인 모델을 통해 반사 계수를 분석하는 논문이 최근 제안되었다( Bae et al., 2013 ). 1917년 I. Schur에 의해 제안된 Schur 알고리즘은( Kailath, 1987 ) Multirate 신호처리( Vetterli, 1987 ; Vaidyanathan, 1993 ), 역산란 문제( Bruckstein and Kailath, 1987 ), 음성 신호처리( Markel and Gray, 1976 ), 광 산란 필터( Dowling and MacFarlane, 1994 ), 지층 Profiling( Markel and Gray, 1976 ), 그리고 해수의 음향속도 Profiling ( Schwetlick, 1983 )과 같은 다양한 문제들에 적용되었다. 특히 Schur 알고리즘은 해수층과 해저지층의 다층 매질을 모델화 할 수 있는 Lattice 필터 구조의 반사계수를 빠르고 쉽게 구할 수 있고 고정된 Wordlength의 수행에서 낮은 Roundoff 잡음 특성과 양자화 잡음에 민감하지 않다는 특성을 갖는 장점이 있다 ( Alexander, 1986 ). 그러나 오차 및 잡음 등에 의해 행렬이 Positive definte 하지 않은 경우에 추정오차가 발생하여 알고리즘이 발산할 수 있는 약점이 있다. 이러한 단점의 극복을 위해 본 논문에서는 Schur 알고리즘에 임계치 설정 방법을 적용한 고속 해저지층탐사 알고리즘을 제안한다. 알고리즘 검증을 위해 비선형 음향 특성을 수치 해석적으로 분석한 KZK(Khokhlov-Zabolotskaya-Kuznetsov) 모델을 사용하여 비선형 송신 신호를 생성하였고, 생성된 송신 신호를 Lattice 필터 구조로 모델링된 해저지층을 투과 및 반사되게 하여 수신 신호를 생성하였다. 끝으로 제안된 임계치가 고려된 Schur 알고리즘과 생성된 송신 및 수신된 데이터 행렬을 이용하여 해저지층의 반사계수의 변화를 추정하게 된다. 본 논문의 2장에 KZK를 이용한 송신 신호 생성 및Llattice 필터 구조로 해저지층을 모델링하여 송수신 신호를 합성하는 방법을 제시하였고, 3장에 합성된 데이터로 부터 해저지층의 매질 특성인 반사계수의 변화를 추정하는 임계치가 적용된 Schur 알고리즘을 설명하였다. 4장에서 제시된 신호의 전파 특성 및 모델링된 해저지층의 전달특성으로부터 얻은 송수신 데이터를 이용하여 제안된 임계치가 적용된 Schur 알고리즘을 검증하였고 5장에서 결론을 맺는다. 탐사를 위해 사용되는 송신 신호는 파라메트릭 배열 현상에 의해 비선형적으로 생성된다. 따라서 파라메트릭 배열 현상을 해석하기 위해서는 일반적인 선형 미분방정식이 아닌 비선형 모델 미분방정식이 요구된다. 해석을 위한 비선형 모델에는 유체의 국소적 효과는 무시하고 거리를 진행하면서 축적되는 효과만을 고려한 Westervelt 방정식, 음파의 회절 효과를 제외한 비선형과 감쇠효과를 고려한 Burgers 방정식 그리고 1차 주파수의 고지향성과 거리에 따른 느린 변화를 가정한 KZK방정식 등이 존재한다. 지향성 빔의 회절, 감쇠, 그리고 비선형 현상을 고려한 KZK 방정식은 직교 좌표계에서 음원이 xy 평면에 위치한 후 z 방향으로 방사될 때 식 (1)과 같다( Hamilton and Blackstock, 1998 ; Lee, 1993 ). 여기서, τ = t - z/c ₀ 는 retarded time, δ 는 음향확산도(sound diffusivity), c ₀ 는 음속, 그리고 ρ ₀ 는 매질의 밀도를 나타낸다. 그리고 β 는 비선형 계수(Nonlinearity coefficient)로 해수는 3.8, 담수는 3.6, 그리고 공기 중은 1.2의 값을 갖는다( Hamilton and Blackstock, 1998 ). 1차 주파수의 비선형 현상에 의해 간접적으로 발생되는 차 주파수의 가상음원은 감쇠거리(Absor- ption range)와 회절거리 (Rayleigh distance)에 영향을 받는다. 감쇠거리 R_A 는 간접적으로 발생되는 차 주파수의 가상음원이 선 배열(Line array) 형태로 누적되는 구간으로 식 (2)와 같이 전 흡음계수( α_T : Total absorption coefficient)의 역수로 계산된다( Moffett and Mellen, 1977 ). 여기서, α _f1,f2 는 1차 주파수의 흡음계수(absorption coefficient) 그리고 α_s 는 차 주파수의 흡음계수이다. 본 논문에서 흡음계수의 계산은 Fisher와 Simmons의 모델을 사용하였다( Fisher and Simmons, 1977 ). 그리고 근거리 음장과 원거리 음장의 경계인 회절거리 R_R 는 식 (3)과 같다( Novikov et al., 1987 ). 여기서, S 는 트랜스듀서의 방사면적 그리고 λ_p 는 평균 1차 주파수( f_p = ( f ₁ + f ₂ )/2)의 파장이다. 회절거리를 벗어난 차 주파수의 음원은 더 이상 누적되지 않고 확산(Spreading)현상에 의해 감쇠된다. 본 논문에서는 위 식 (1)~(3)을 이용하여 송신 신호를 생성하였다( Lee et al., 2012 ). 해저지층 채널을 모델링하기 위해서는 신호의 투과 및 반사 를 적절히 고려해야만 한다. 이러한 신호 전달 특성을 잘 모델링 할 수 있는 것은 Lattice 모델이며 해저 지층탐사에 매우 유용한 모델임이 알려져 있다( Bae et al., 2013 ). 다층 채널에 의해 생기는 역산란 현상은 Fig. 1(a) 과 같이 손실이 없는 서로 다른 매질의 특성을 갖는 두층을 고려할 때, ”1”을 수직으로 입사하 였을 때 반사계수에 따라 일부에너지는 R = r ₁ 만큼 반사되고 나머지 에너지는 만큼 투과 되는 현상을 말한다 ( Bruckstein and Kailath, 1987 ). 이러한 현상은 Fig. 1(b) 와 같이 Lattice 필터 구조로 모델링 된다( Orfanidis, 1988 ). 여기서 ρ는 매질의 밀도를 나타내고, c 는 매질의 음속을 말하며, 이를 이용하여 아래 식 (4)와 같이 반사계수를 계산할 수 있다. Fig. 1 과 같은 산란현상의 수학적인 모델은 아래 식 (5)와 같이 행렬식으로 나타낼 수 있다( Orfanidis, 1988 ). 여기서 은 투과계수라고 하고, 이는 에너지 보존의 법칙을 만족한다. 만약 매질이 Fig. 2(a) 와 같이 일정한 시간지연을 가진 N +1층으로 변한다고 가정하면, 각 매질은 Fig. 2(b) 와 같이 lattice 구조로 모델링을 할 수 있다. 그리고 식 (5)를 이용해 식 (6)과 같은 행렬식으로 확장할 수 있다( Orfanidis, 1988 ; Bae et al., 2005 ; Bae, 2007 ; Choi et al., 2011 ). 여기서 는 디지털신호처리의 개념으로 시간 지연 τ 와 각주파수 w 에 의해 계산되는 단위시간지연(Unit delay)이다( Kailath, 1987 ; Dowling and MacFarlane, 1994 ; Bae et al., 2005 ). 식 (6)에서 U_i 와 V_i 는 i 층에서의 송신 및 수신된 시계열 신호이며, 각 층에서의 송신 및 수신의 시간지연은 τ /2이다. 이 시간지연 τ /2은 수신 시계열 신호의 샘플링 시간의 반이며, 이것은 신호가 송신되고 반사되는 시간을 모두 포함하는 의미이다. 또한 샘플링 시계열 신호의 개수가 N일 때 매질의 층은 N +1로 가정하여 Lattice 구조로 구성된다. 결과적으로, 다층 모델의 송수신 신호는 Fig. 2 의 모델에 따라 식 (6)을 이용하여 다층 매질에 입력 시계열 신호 U ₀ = { u ₀ , u ₁ , ⋯} _T 가 입사 할 때 반사 시계열 신호 V ₀ = { v ₀ , v ₁ , ⋯} _T 를 얻어 구한다. Fig. 2 에서 보인 다층 구조의 해수층 및 해저지층에 따른 매질의 변화는 수신 데이터 취득의 샘플링에 따라 변화된다( Bae et al., 2013 ). 본 논문에서 사용된 해저 지층 모델은 해수층과 해저지층을 고려한 Lattice 모델이다( Bae et al., 2013 ). 모델링된 해저 지층은 해수층과 다섯 층으로 구성된 퇴적층으로 개별 퇴적층의 특성은 표 Table 1 과 같다( Lurton, 2002 ). 각 해저 지층의 수심에 대하여 0m에서 402.5m를 Layer1, 402.5m에서 404m를 Layer2, 404m에서 405.5m를 Layer3, 405.5m에서 406.7m를 Layer4, 406.7m에서 411.9m를 Layer5 그리고 411.9m이상을 Layer6으로 가정하였다. 또한 수심에 따른 음속과 밀도의 분포를 Fig. 3 과 같이 가정하였는데 이때 수심은 음속 1500m/s를 가정하여 정규화 하였다. 반사계수는 식 (4)와 Fig. 3 의 특성을 적용하면 계산할 수 있다. 여기서 해수층의 반사계수 변화는 매질의 변화가 매우 적어 해저 퇴적층의 반사계수 변화에 비해 무시할만한 변화를 갖는 다고 가정할 수 있다( Bae et al., 2013 ). Fig. 4 는 퇴적층이 존재하는 400m에서 415m 깊이의 반사계수를 나타낸다. 반사계수는 샘플링 시간을 2μs 단위로 나누어 계산하였으며, 총 반사계수의 수는 5,000개가 된다. 계산된 반사계수를 이용하여 Fig. 2 의 Lattice 모델에 송신 신호를 입사 시키면 수신 신호를 얻을 수 있다. 따라서 반사되어 돌아오는 시간을 고려하면 수신 시계열 신호의 샘플 개수는 10,000개가 된다. 본 논문에서 사용된 해저 지층의 반사계수 추정 알고리즘은 임계치가 고려된 Shur 알고리즘을 기반으로 한다( Bae et al., 2013 ). 이 알고리즘은 송신과 수신의 시계열 신호로부터 에너지 보존의 법칙을 만족하는 Autocorrelation 행렬을 식 (7)과 같이 구성하여 원하는 행렬 T_A 를 계산한다. T_A 행렬은 Toepliz 행렬로 Positive definite를 만족해야한다( Kailath, 1987 ). Schur 알고리즘은 일반적으로 주어진 식 (7)의 T_A 행렬을 고속으로 Cholesky 분해하여 T_A = LL^T 로 계산하는 알고리즘이다. 여기서 L 은 하부삼각행렬(lower triangular matrix)이고, T 는 행렬의 전치(transpose)라고 한다. 행렬 T_A 는 Toepliz 행렬의 형태로 아래의 식 (8)과 같이 생성 행렬인 G ₀ = [ U ₀ V ₀ ]로 나타낼 수 있다. 식 (8)에서 S 와 J 는 아래와 같이 정의된 천이행렬과 J-unitary 행렬이다. 이렇게 만들어진 행렬 G ₀ 를 이용하여 아래와 같이 Schur 알고리즘을 이용하여 반사계수를 계산한다. 아래에 정리된 Schur 알고리즘은 식 (6)을 이용하여 바로 반사되는 신호로부터 반사계수를 찾아 하나씩 경계면을 제거하고, 다시 반사계수를 구하는 방식이다. 우선 r ₀ 는 첫 번째 송신신호와 첫 번째로 취득한 수신신호의 비로, 이것을 이용하여 역행렬을 구하여 다음 식 (10)과 같이 첫번째 경계면을 제거한다. 이때 구해진 r ₀ 값이 오차 및 잡음 등에 의해 정해진 임계치 이하가 되면 0으로 세팅한다. 그리고 시간지연에 따른 신호가 이동한 매질을 제거하고, 다음 경계면을 제거하는 방법으로 반사계수를 계산하는 방법이 아래에 임계치가 고려된 Schur 알고리즘으로 샘플 시간에 의해 계산되는 거리 단위로 매질을 분석할 수 있다. Threshold based Schur Algorithm ( Koltracht and Lancaster, 1988 ) 위에 설명한 것과 같이 임계치가 고려된 Schur 알고리즘을 이용하여 주어진 송신신호 U ₀ 와 수신신호 V ₀ 로부터 반사계수 { r _o , r ₁ , ⋯} 를 안정적으로 계산 할 수 있다. 즉, 위 임계치 조건이 적용이 되면 각 G_i , i = 0, 1, 2, ⋯ 의 첫 번째 열(Column) 벡터를 순차적으로 나열하여 만든 L 행렬의 대각선(Diagonal) 성분은 이 허수가 되는 상황이 발생이 되지 않아 Cholesky 분해 가능하여 송신 및 수신신호로 만들어진 행렬 T_A 는 Positive definite 하게 된다. 본 논문에서 파라메트릭 배열 송신 신호 생성을 위해 사용된 파라미터는 아래 Table 2 와 같다. 주어진 파라미터를 사용하여 송신 신호는 잘 알려진 KZK 모델을 이용하여 생성되었다. Table 2 의 파라미터를 사용할 경우, 두 개의 1차 주파수 265kHz와 235kHz의 흡음 계수는 각각 0.0104Np/m와 0.0096Np/m로 계산되고, 전 흡음계수는 0.0195Np/m로 계산되며, 식 (2)로부터 얻어진 감쇠거리는 51.28m이다. 또한 식 (3)에 의해 회절거리는 11.64m가 되는데, 이 경우 차 주파수 생성을 위한 가상음원 거리는 회절거리가 감쇠거리보다 짧기 때문에 회절거리까지 음압이 누적되어 증가된 후 회절거리 이후부터 확산현상에 의해 감쇠하게 된다( Lee et al., 2012 ). Fig. 5 는 이와 같은 파라메트릭 배열 현상으로 생성된 송신 신호의 거리에 따른 감쇠현상을 보인다. Fig. 5 에서 보인 바와 같이 송신 신호는 해수층을 전파하면서 주파수의 자승에 비례하는 전파손실로 감쇠되고, 매질의 비선형 현상에 의해 간접적으로 생성된 파라메트릭 배열 송신 신호는 회절거리인 약 11m까지 음압이 누적되어 증가된 후 전파거리가 증가함에 따라 서서히 감쇠되어 약 387m 이상의 지점부터는 1차 주파수의 음압준위보다 높은 음압준위를 보임을 알 수 있다. 또한, 약 400m 지점의 송신 신호가 100m 전파될 때 1차 주파수 235kHz의 음압준위는 184.2dB에서 176.5dB로 7.7dB 감소된 반면, 차 주파수 30kHz의 음압준위는 184.3dB에서 184dB로 0.3dB 감소되었음을 알 수 있다. Fig. 6 은 스펙트로그램을 통해 거리에 따른 1차 주파수 및 차 주파수의 감쇠를 보인 것이다. Fig. 7 은 KZK 모델를 이용하여 생성된 송신 신호가 거리가 0.1m와 400m 지점에서의 파형과 주파수 특성을 나타낸 것이다. 400m 지점에서 신호 파형은 비선형 현상에 의해 왜곡되었음을 확인 할 수 있으며, 차 주파수의 음압 준위가 1차 주파수 265kHz의 음압 준위보다 8.4dB 크고 235kHz의 음압 준위보다 0.1dB 큼을 주파수 특성으로부터 확인할 수 있다. 본 논문에서는 Fig. 7 의 (b) 신호를 입사신호로 사용하였고 Fig. 8 과 같은 구조로 형성된 퇴적층을 가정하였다. Fig. 8 에서 보는 봐와 같이 퇴적층은‘Trace 0’~‘Trace 22’번째 까지는 6개의 층으로, ‘ Trace 23 ’ ~ ‘ Trace 39 ’ 번째 까지는 5개의 해저 지층으로 구성되어 있다. Fig. 8 의 지층 구조에 따른 반사계수 특성은 표1 의 해저 지층 특성 및 식 (4)를 이용하면 얻을 수 있으며 그 결과는 Fig. 9 와 같다. 끝으로, 반사된 수신 신호는 Fig. 7 (b) 의 신호와 Fig. 9 의 반사계수 정보가 적용된 Lattice 채널 모델을 이용하면 얻을 수 있다. 일례로 ‘ Trace 1 ’ 을 통과하여 반사된 신호는 Fig. 10 과 같다. 이때 송신 신호가 퇴적층이 분포하는 지점 400m에서 415m지점까지 투과 및 반사되어 돌아오는 시간은 20ms이다. 다음 Fig. 11 는 본 논문에서 제안한 임계치가 고려된 알고리즘의 성능을 입증하기 위해 Fig. 4 에 사용된 퇴적층의 반사계수를 확대한 것이다. 그림에서 보는 바와 같이 크기가 큰 반사계수의 값 사이에 매우 작은 반사계수 오차가 존재함을 알 수 있다. 추정된 반사계수 오차의 표준편차는 0.0001이다. 따라서 본 논문에서는 이러한 오차를 고려하여 임계치를 0.001로 선정하였으며 이 임계치 값 이내의 값은 0으로 설정하였다( Koltracht and Lancaster, 1988 ). Fig. 12 은 본 논문에서 제안된 임계치가 고려된 Schur 알고리즘의 반사계수 추정결과이다. Fig.12 에서 실선은 설계된 반사계수, 점선은 추정된 반사계수를 나타내는데 정확하게 반사계수 및 퇴적층의 수심을 추정함을 알 수 있다. 본 논문에서 정규화된 퇴적층의 수심은 해수층의 음속 1,500m/s을 기준으로 정규 화된 수심을 의미한다. Fig. 13 는 본 임계치가 고려되지 않은 Schur 알고리즘의 반사계수 추정결과이다. Fig.13 에서 보는 바와 같이 반사계수 오차에 의한 영향으로 반사계수가 추정되지 않음을 알 수 있다. Fig. 14 과 15 는 수신된 신호를 이용하여 모델링된 해저 지층의 반사계수를 추정한 결과이다. Fig. 14 은 후처리를 하지 않은 수신 시계열 신호이며 Fig. 15 는 수신 시계열 신호에서 차단주파수 40kHz의 FIR 저역통과필터를 이용하여 여과한 후 차 주파수에 해당되는 30kHz 대역의 신호를 정합하여 추정된 해저 지층의 반사계수 분포이다. 이 두 그림을 Fig. 9 와 비교하여 보면, 대체적인 지층의 반사계수를 파악할 수 있으나 정확한 지층의 반사계수 추정은 불가능함을 알 수 있다. 특히 Layer 3과 Layer 4가 겹치는 ‘ Traces 23 ’ 의 미세한 반사계수의 변화는 구분하기 매우 어려움을 알 수 있다. Fig. 16 와 17 은 본 논문에서 제안된 임계치를 고려한 Schur 알고리즘을 이용하여 반사계수를 추정 결과와 기존의 논문( Bae et al., 2013 )에서 제안된 결과를 보인 것이다. 제안된 Schur 알고리즘을 이용한 반사계수 추정 기법은 정확한 지층 분포 변화 및 매질의 특성 변화를 파악할 수 있다. 그러나 반사계수 오차의 영향이 고려되지 않은 알고리즘을 사용할 경우에는 알고리즘이 발산하여 전혀 결과를 얻을 수 없음을 알 수 있다. 본 논문에서는 파라메트릭 배열 송신 신호 및 Lattice 모델을 이용하여 수신된 신호를 이용하여 다층으로 구성된 해저지층의 반사계수를 안정적으로 추정할 수 있는 알고리즘을 제안하였다. 제안된 알고리즘은 모델링 오차에 기인한 잡음에 둔감한 문턱치가 적용된 Schur 알고리즘으로, 송신과 수신신호로 만든 행렬 T_A 가 Positive definite 조건이 성립되지 않는 상황, 즉 Schur 알고리즘이 발산하여 추정된 반사계수의 절대값이 1보다 큰 경우에도 적용이 가능한 알고리즘이다. 모의실험을 위해 KZK 모델을 이용하여 송신 신호 생성하였고 모의실험 결과 문턱치가 고려되지 않은 알고리즘보다 안정적인 반사계수 추정이 가능함을 검증하였다. 또한, 제안된 알고리즘이 일반 신호 및 정합필터를 사용하여 반사계수를 추정하는 것보다 매우 정밀하고 효율적으로 해저 지층 분해하고 탐사할 수 있음을 확인하였다.