I. 서 론

소나는 수중에서 음파를 이용하여 수중 표적을 탐지하는 장비이며, 따라서 수중음향 전달 채널은 소나의 성능을 좌우하는 중요한 요인이다. 수중음향 전달 채널은 음원에서 송신한 신호가 소나의 수신기에 도달하였을 때, 그 크기와 위상이 달라지는 것을 표현하며, 해수면 및 해저면의 경계조건, 해저면의 구성 성분, 수온 분포 등의 다양한 환경 요인들의 영향을 받는다. 이러한 수중환경의 복잡성으로 인해 다양한 수중음향 채널에 따른 소나의 계측 신호를 모의하는 것은 쉽지 않다.

일반적으로 해양 환경에 대한 매질 특성과 경계 조건이 주어진 경우, 물리적인 파동방정식의 해를 구함으로써 수중음향 전달 채널을 예측할 수 있다.^[1] 다만 이러한 파동방정식은 종종 주파수 영역에서 표현되어, 광대역 시계열 신호 모의시 바로 활용이 어렵다.^[2,3] 반면 음선 모델에서는 고유 음선별 도달 시간 및 진폭을 계산할 수 있으므로, 다중경로 신호의 합으로부터 광대역 신호를 손쉽게 모의할 수 있다.^[3]

다만 음선모델을 활용하는 경우, 물리적 환경 파라미터(음속 프로파일, 경계 특성 등)를 채널 모델에 손쉽게 반영할 수 있으며, 선형 시불변 가정 하에서 수중음향 전달 채널의 임펄스 응답만을 가지고 모든 정보를 간단하게 표현할 수 있는 장점이 있다.^[4]

본 논문에서는 음선 모델로부터 얻어진 고유 음선 분석 결과를 소나 수신기에서의 유한 임펄스 응답으로 근사하여 모의하는 방법을 제안하였다. 특히 고유 음선의 위상과 손실, 그리고 연속시간에서 산출된 음선의 도달 시간지연을 이산 시간으로 근사할 때의 오차를 최소화 할 수 있는 모의 방법을 제안하였다.

II. 고유 음선 분석

음선 이론은 파동 방정식의 고주파수 근사 해로서 얻어지지만, 그 결과가 직관적으로 이해가 쉬우며 연산 부담이 적다는 점에서, 실제로 많이 사용된다. 다만 고주파 근사를 적용함으로써, 저주파 영역에서는 적용이 제한될 수 있다.^[1]

음선 해는 파동 방정식을 아래와 같이 두 개의 방정식으로 분리함으로써 얻어진다. 공간의 좌표를 $\mathbf{x}=(\mathrm{x},\mathrm{y},\mathrm{z})$, 각주파수를 $\omega$, 매질의 음속을 $c\left(\mathbf{x}\right)$라 하면, 음선 계산을 통해 얻는 음압 $p\left(\mathbf{x}\right)=A\left(\mathbf{x}\right)e^{j\omega \tau \left(\mathbf{x}\right)}$는, $\left|{\nabla }^{2}A/A\right|\ll {\omega }^2/c^2\left(\mathbf{x}\right)$라는 가정 하에서 다음 방정식을 만족한다.^[1]

Eq. (1)은 아이코날 방정식으로 불리며, 파면에 수직하여 진행하는 음선의 경로 및 그에 따른 도달 시간 지연 $\tau \left(\mathbf{x}\right)$를 결정한다. Eq. (2)는 수송 방정식이며, 음선의 진행에 따른 에너지의 확산 및 수축을 모델링하여 음선 진폭 $A$를 결정한다. 한편 위 방정식으로 주어진 음선 근사는 $\left|{\nabla }^{2}A/A\right|\ll {\omega }^2/c^2\left(\mathbf{x}\right)$의 가정에서 성립하므로, 일반적으로 저주파수 영역에서는 근사 정확도가 떨어지게 된다. 특히 음속 $c\left(\mathbf{x}\right)$의 구배가 한 파장 안에서 급격히 변하는 경우, 근사가 맞지 않는 것으로 알려져 있으며, 음속 구배의 변동이 작은 경우에도, 진폭의 공간적인 변동이 큰 표면 도파관이나 음선 집중(caustic) 영역에서 정확도가 떨어짐이 알려져 있다.^[5] 일반적으로 음선을 이용한 계산 결과는 저주파 · 원거리 영역에서 정확도가 저하되므로, 이 경우 정상 모드 또는 포물선 방정식 기반의 모델의 적용을 고려하여야 한다.^[1]

그러나 음선을 이용한 계산 결과는 이상적인 경계 조건에서 등음속 프로파일이 주어진 경우 가상 음원을 통해 얻은 결과와 동등하므로, 저주파를 포함한 전 대역에서 정확한 결과를 얻을 수 있다. 또한 경계면 조건이 달라지는 경우에도, 균질한 평면 경계면을 가정하는 경우, 주파수 대역에 따른 영향을 크게 받지는 않는다. 이를 확인하기 위해, Fig. 1의 Pekeris 도파관에서, 정상 모드에서 계산된 음압과 음선으로부터 도출한 음압을 100 Hz ~ 5000 Hz의 주파수에서 비교하였다.^[1] 이 때 음원과 수신기 사이의 거리를 5 km의 원거리로 설정하고, 정상 모드 중 진행 모드만을 음압 계산에 활용하였다.^[1] 음선 계산의 경우 수평면 기준 ± 30°이내의 401개의 경로를 이용하였다.

Fig. 1.

Pekeris waveguide with the water depth D.

Fig. 2는 수심 $D$가 1000 m인 심해의 상황이며, Fig. 3은 수심이 50 m인 천해 환경이다. 다중경로가 많지 않은 심해 환경에서는 음선 근사를 적용한 결과가 정상모드의 해로부터 계산된 정확한 결과와 큰 차이를 보이지 않는다. 다만 복잡한 다중경로가 발생하는 천해 환경에서는, Fig. 3에서 볼 수 있듯이 음선 근사에 따른 오차가 상대적으로 크게 발생한다. 따라서 음선을 이용한 근사는 다중경로 환경이 복잡한 경우에는 적용이 일부 제한될 수 있다. 그렇지만 이상적인 수평면 경계면 환경에서, 음속 구배의 변화가 심하지 않은 경우에는, 제한적으로 음선 이론을 저주파 대역까지 확장하여 적용할 수 있다.

Fig. 2.

(Color available online) Frequency response of acoustic field at the water depth of 1000 m. (a) amplitude, (b) phase response.

Fig. 3.

(Color available online) Frequency response of acoustic field at the water depth of 50 m. (a) amplitude, (b) phase response.

본 논문에서는 $\left|{\nabla }^{2}A/A\right|\ll {\omega }^2/c^2\left(\mathbf{x}\right)$ 조건이 모든 주파수 대역에서 충족된다고 가정하고, 음선 계산 결과로부터 광대역 신호를 모의하는 방법을 제안하였다. 이 경우 매질의 속력 $c\left(\mathbf{x}\right)$이 주파수에 따라 변하지 않는다고 가정하면, Eqs. (1)과 (2)는 주파수와 무관하므로, 음선 경로로부터 얻어진 해는 광대역 소음의 모든 주파수에서 동일하게 적용할 수 있다.

일반적으로 수중환경에서 음속은 주파수에 영향을 받는 인자는 아니지만,^[1] 음속의 허수부로 표현되는 체적 손실은 주파수의 영향을 받는다.^[1] 다만 음선 경로 $\tau \left(\mathbf{x}\right)$의 섭동 급수해 $\tau ={\tau }_0+ϵ{\tau }_1+{ϵ}^2{\tau }_2\cdots$를 가정하여 2차항 이하를 무시하는 경우, 체적 손실은 음선 해의 진폭의 절대값 $\left|A\left(\mathbf{x}\right)\right|$에만 영향을 주며, 음선의 경로 $\tau \left(\mathbf{x}\right)$ 및 위상 $\measuredangle \left|A\left(\mathbf{x}\right)\right|$은 영향을 받지 않는다.^[1] 이에 따라 체적 손실의 크기만을 반영하는 경우 음선의 음압 $p\left(\mathbf{x}\right)$는 아래와 같이 표현된다.^[1]

단, $\alpha \simeq 0.0033+\frac{0.11f^2}{1+f^2}+\frac{44f^2}{4100+f^2}+0.0003f^2$

$\alpha :$감쇄율 $[dB/km]$, $f$: 주파수 $\left[kHz\right]$

$L={\int }_{}^{}ds$: 음선 경로길이 $\left[km\right]$, $\frac{d\mathbf{x}}{ds}=c\nabla \tau$

한편 적층된 이상적인 경계면에서 평면파 입사를 가정하는 경우, 경계면 반사 계수의 위상 및 크기는 주파수에 무관하게 입사각에 따라 주어지며, 음선을 평면파로 가정하는 경우 이를 음선의 진폭 $A$에 반영할 수 있다.^[1] 다만 입사파가 평면파가 아니라는 점은 경계면 작용시 음선의 정확도를 저하시키는 원인이 되며, 일반적으로 적층된 경계면에서 발생하는 음선의 복반사 등을 음선 분석에 반영하는 것도 쉽지 않다.^[1,3]

그 외에 음선 이론의 한계점으로는 앞서 언급한 음선 집중 및 음영 구역을 들 수 있다. 음선 집중 시에는 진폭이 무한대가 되며 음영구역에서는 진폭이 0이 되므로, 고유 음선이 계산되지 않는다. 이를 보완하기 위해서 가우시안 빔 경로 기법이 제안되었다.^[6] 한편 기존의 기하학적 음선 기법을 적용하는 경우, 음선 집중이 발생한 이후 위상 $\pi /2$를 보상해주어야 한다.^[1]

한편 위에서 언급된 모든 음선 분석은 인터넷에 공개된 Bellhop 프로그램을 이용하여 손쉽게 수행할 수 있다.^[3] Bellhop 프로그램에 음속 프로파일, 해저 지형, 경계면에서의 반사 특성 및 송수신기 위치를 입력하면, 고유 음선의 도달 시간지연, 진폭의 크기 와 위상, 고각, 반사 횟수 등의 정보를 얻는다.^[3]

III. 음선의 위상 및 손실 근사

이상적인 경계면을 가정하여 주파수 별로 동일하게 음선의 위상이 주어진 경우, 반사된 음파는 원래 신호와 위상이 $\pi /2$만큼 천이된 힐버트 변환 신호의 가중 합으로 표현된다.^[1,2,3] 이산시간에서 힐버트 변환의 주파수 응답 $H_h\left(e^{j\omega }\right)$는 다음과 같이 주어진다.^[4]

그런데 힐버트 변환의 위상과 크기를 정확하게 만족하는 유한 차수 필터는 존재하지 않으므로, 적절한 근사가 요구된다. 한편, 4형 선형 위상 필터(linear phase filter, Type 4)의 경우, 힐버트 변환과 동일한 위상이 보장되며, $\omega =0,\pi$에서 주파수 응답이 영이 되는 3형 필터와는 달리 $\omega =0$에서만 영이 되므로, 힐버트 변환 근사에 유리하다.^[4] 이 경우 필터 차수는 홀수 임이 요구된다.

한편 차수 $N$인 임의의 필터를 (N + 1)⨯1 크기의 벡터 $\mathbf{h}=\left[h\right(0\left)h\right(1)...h(N){]}^T$로 나타내면, 이는 $N$이 홀수인 경우 다음과 같이 2형 및 4형 선형 위상 필터의 합으로 분해된다.

단, $h_2\left(n\right)=h_2(N-n)$ (even symmetric),

$h_4\left(n\right)=-h_4(N-n)$ (odd symmetric).

4형 선형 위상 필터는 2형 필터에 비해 $\pm \pi /2$만큼 위상 천이가 발생하므로, 임의 위상의 복소 진폭은 $A\left(e^{j\omega }\right)\simeq H_2\left(e^{j\omega }\right)\pm j\left|H_4\left(e^{j\omega }\right)\right|$으로 근사된다.^[4] 단, $N$차의 선형 위상 필터의 적용에 따른 시간지연 $N/2$은 보상된 것으로 가정하였다.

그런데 2형 필터의 경우, $\omega =\tau$에서 주파수 응답이 0이 된다. 따라서 $\omega =0,\pi$에서는 임의의 크기와 위상을 갖도록 필터를 근사하기 어려우므로, 이를 제외한 통과 대역에서 근사가 요구된다. 또한 필터 계수가 실수(real)이면, 주파수 응답의 대칭성을 활용하여 양의 주파수($0<\omega <\pi$) 대역만을 최적화 대상으로 할 수 있다.

한편 주파수 영역에서 필터 응답을 최적 근사하기 위해 최소제곱 방법을 적용할 수 있다. 이 방법은 근사가 요구되는 주파수 대역에서만 최적화를 실시하므로, 근사 대역 내에서 정확도가 증대되는 장점이 있다.^[7] 최소 제곱 방법의 적용을 위해 근사 목표인 주파수 영역에서 이상적인 필터응답을 $H_{id}\left(e^{j\omega }\right)$라고 하면, 가중 최소제곱 오차 $E$는 다음과 같이 주어진다.^[7]

단, 오차가 정의된 근사 대역은 $\left[\alpha \pi ,\beta \pi \right]$이며, $W\left(\omega \right)$는 값이 0 이상인 주파수영역 가중함수이다.

위 오차를 최소화하는 차수가 $N$인 최적 필터 계수를 구하기 위해, (N + 1)⨯1 크기의 벡터로 주어지는 푸리에 변환 계수 $\mathrm{e}$를 아래와 같이 정의한다.

이를 대입하여 정리하면 Eq. (6)의 가중최소제곱 오차 $E$는 아래와 같이 차수가 N인 필터를 나타내는 (N + 1)⨯1 크기의 벡터 $\mathbf{h}$의 이차식으로 주어진다.^[7] 이차식의 계수는 (N + 1)⨯(N + 1) 크기의 행렬 $\mathbf{P}$와 (N + 1)⨯1 크기의 벡터 $\mathbf{q}$ 및 스칼라 $r$로 표기하였다.

단, $\mathbf{P}=\frac{1}{\pi }{\int }_{\alpha \pi }^{\beta \pi }W\left(\omega \right)\mathrm{Re}\left\{{\mathbf{ee}}^H\right\}d\omega$,

$\mathbf{q}=\frac{1}{\pi }{\int }_{\alpha \pi }^{\beta \pi }W\left(\omega \right)\left[\mathrm{Re}\left\{H_{id}^{\mathit{*}}\mathit{\left(}{\mathit{e}}^{\mathit{j}\mathit{w}}\mathit{\right)}\mathbf{e}\right\}\right]d\omega$,

$r=\frac{1}{\pi }{\int }_{\alpha \pi }^{\beta \pi }W\left(\omega \right){\left|H_{id}\left(e^{jw}\right)\right|}^{2}d\omega$.

위 식에서 어깨 글자인 ${}^H$는 복소공액 전치행렬, ${}^T$는 전치행렬, $*$는 복소공액을 의미한다. Eq. (8)는 최적화 대상인 필터 계수 $\mathbf{h}$의 이차식이므로, 그 구배를 0으로 두면 필터 계수를 다음과 같이 얻는다.^[7]

만약 가중함수를 단순하게 $W\left(\omega \right)=1$로 두면, 이차식의 계수 행렬 $\mathbf{P}$의 원소 $P_{ij}$는 이상적인 필터 응답 $H_{id}\left(e^{jw}\right)$와 무관하게 다음과 같이 주어진다.

단, $i,j=0,1,2,...,N$.

한편, 이상적인 필터 응답인 $H_{id}\left(e^{jw}\right)$에는 고유 음선으로부터 계산된 체적 손실, 해수면 및 해저면 반사에 따른 손실 및 위상, 그리고 음선 집중에 의한 위상 천이가 반영된다. 이는 일반적으로 주파수의 함수이므로, 최적 필터 설계를 위해서는 Eq. (8)에 따른 벡터 $\mathbf{q}$의 적분을 수치적으로 계산하여야 한다.

만약 음선의 손실과 위상이 주파수에 무관하게 주어져, 진폭이 $A=R{e}^{j\varphi }$으로 고정된 경우, 이상적인 응답은 $H_{id}\left(e^{jw}\right)=R{e}^{j\varphi }{e}^{-jN\omega /2}$로 표현할 수 있다. 그러면 Eq. (8)에 따른 벡터 $\mathbf{q}$의 원소 $q_n$은 다음과 같은 수식으로 주어지며, 최적 필터 응답은 Eq. (9)을 사용하여 계산할 수 있다.

단, $n=-N/2,-N/2+1,...,N/2$, $N$은 홀수.

설계된 필터의 특성을 알아보기 위해, 근사 대역을 $\alpha =0.1,\beta =0.9$로 하고 필터 차수를 $N=25$으로 하여 분석을 실시하였다. 위상만이 천이되는 경우를 가정하면 $H_{id}\left(e^{j\omega }\right)=e^{j\varphi }{e}^{-jN\omega /2}$이며, 주파수 영역 최소제곱 기법으로 설계된 필터의 주파수 대역 응답을 Fig. 4에 나타내었다. 단, 시간지연 $N/2$에 해당하는 위상 $e^{-jN\omega /2}$은 보상하였다.

Fig. 4로부터 위상이 $\varphi =\pi /2$ 또는 $\varphi =0,\pi$에서 설계된 필터의 위상 오차가 없음을 알 수 있다. 이는 4형 또는 2형 선형위상 필터가 Eq. (6)으로 주어진 조건 하에서도 최적이기 때문이다. 한편 위상이 $\pi /2$의 배수인 특수한 경우를 제외하면, 유한 차수 필터의 특성상 근사 대역 밖에서는 오차가 발생함을 확인할 수 있다.

Fig. 4.

(Color available online) Frequency response of the phase shift filter. (a) magnitude, (b) phase response.

추가로 표적과의 거리가 10 km, 샘플링 주파수가 30 kHz인 경우에 Eq. (3)으로 주어진 체적 손실을 반영하여 설계된 필터 응답을 Fig. 5에 나타내었다. Eq. (3)에서 체적 손실은 진폭의 크기에 반영되며, 위상은 영향을 받지 않는다. 즉, $H_{id}\left(e^{j\omega }\right)={10}^{-\alpha L/20}{e}^{j\varphi }{e}^{-jN\omega /2}$으로 주어진다. Fig. 5에서 Eq. (3)에 따른 체적 손실 크기는 검은색 점선으로 나타내었으며, 필터 설계 결과, 근사 대역 내에서 목표로 한 체적 손실이 정확하게 반영되는 것을 확인할 수 있다.

Fig. 5.

(Color available online) Frequency response of the volume attenuation and phase shift filter. (a) magnitude, (b) phase response.

IV. 음선 도달 시간지연 근사

음선의 도달 시간지연은 Eq. (1)의 미분방정식의 해로부터 계산되므로, 도달 시간지연이 이산시간의 샘플링 주기로 나누어 떨어지지 않는 경우, 나머지 부분 샘플의 시간지연에 대한 반영이 필요하다. 특히 배열내의 센서간 시간지연은 상대적으로 작은 값을 가지므로, 샘플 주기내의 작은 시간지연 근사 오차가 빔형성 등의 성능 평가에 영향을 줄 수 있다.

Shannon의 샘플링 이론에 따르면, 대역 제한된 이산 시간 신호를 sinc 함수를 이용하여 보간하면 원래 신호를 얻는다.^[7] 따라서 이상적인 도달 시간지연 필터는 Fig. 6과 같이 샘플링된 sinc 함수로 주어진다. 이 때 도달 시간지연이 정수이면 필터는 Fig. 6(a)와 같이 하나의 임펄스로 주어지나, 도달 시간지연이 샘플링 주기로 나누어지지 않는 경우에는 Fig. 6(b)와 같이 무한 차수가 된다. 따라서 도달 시간지연을 이산 시간으로 모의하려면, 우선 이를 유한 임펄스 응답으로 근사하여야 한다.^[7]

Fig. 6.

(Color available online) Impulse response of an ideal fractional delay filter with the (a) delay = 0 and (b) delay = 0.5.

한편 III장에서 설계된 필터의 차수가 홀수이므로 $N/2$의 지연이 발생한다. 따라서 도달 시간지연 필터의 차수 역시 홀수이면, 두 필터를 합하여 발생한 정수개의 샘플 지연을 쉽게 보상할 수 있다. 이에 시간지연을 근사하는 필터 $\mathbf{g}=\left[g\right(0\left)g\right(1)...g(M){]}^T$의 차수 $M$은 홀수로 가정하였다.

시간지연 필터는 Fig. 6에서 볼 수 있듯이, 기본적으로 보간 필터이다. 따라서 저주파수 대역 통과 특성을 가지므로, 근사 목표 주파수 영역에서 이상적인 필터응답을 $G_{id}\left(e^{j\omega }\right)$라고 하면, 가중 최소제곱 오차 $E$는 Eq. (6)과 마찬가지로 다음과 같이 주어진다.

단, 오차가 정의된 저주파수 대역은 $\left[0,\beta \pi \right]$이며, $G_{id}\left(e^{j\omega }\right)$는 이상적인 필터응답, $W\left(\omega \right)$는 값이 0 이상인 주파수영역 가중함수이다. 이 때 푸리에 변환 계수는 Eq. (7)에서 정의한 바와 같다. 이를 대입하여 정리하면, 마찬가지 과정을 거쳐 최적 필터 계수는 Eq. (9)와 같이 얻는다.^[7]

이 때 만약 가중함수를 단순하게 $W\left(\omega \right)=1$로 두면, 행렬 $\mathbf{P}$의 원소 $P_{ij}$는 이상적인 필터 응답 $G_{id}\left(e^{j\omega }\right)$와 무관하게 다음과 같이 주어진다.^[7]

단, $i,j=0,1,...,N$.

한편 샘플 주기 내의 나머지 시간지연을 $D$ 샘플($\left|D\right|<1$)이라고 하면, 이상적인 필터의 응답으로 $G_{id}\left(e^{j\omega }\right)=e^{-jN\omega /2}{e}^{-jD\omega }$를 얻는다.^[7] 마찬가지로 가중함수를 $W\left(\omega \right)=1$로 두면, Eq. (8)에 따른 적분 결과인 벡터 $\mathbf{q}$의 원소 $q_n$은 다음과 같이 주어지며, 최적 필터 응답은 Eq. (9)을 사용하여 계산할 수 있다.^[7]

단, $n=0,1,...,N$.

설계된 도달 시간지연 필터의 특성을 알아보기 위해, 근사 대역을 $\beta =0.9$로 하고, 필터 차수를 $M=25$로 두어 분석을 실시하고 그 결과를 Fig. 7에 나타내었다. 단, 시간지연 $N/2$에 해당하는 위상 $e^{-jN\omega /2}$은 보상하였다.

Fig. 7의 결과로부터 근사된 대역 내에서는 설계된 필터가 정확한 위상 시간지연을 가짐을 확인할 수 있다. 특히 도달 시간지연이 $D=0,0.5$인 경우, 전대역에서 정확한 위상지연을 가진다. 한편 근사 대역 이외의 영역($\omega \ge 0.9$)에서는 오차가 발생하나, 주파수 중첩(aliasing) 문제로 인해 신호처리 시 해당 대역을 보통 활용하지 않으므로, 오차를 용인할 수 있다.

Fig. 7.

(Color available online) Frequency response of fractional delay filter. (a) magnitude, (b) phase delay.

한편 Fig. 7에서는 설계된 필터의 위상 시간지연을 도시하였는데, 이는 근사 대역내의 협대역 신호 성분별 시간지연에 관심이 있기 때문이다. 만약 Detection of Envelop Modulation On Noise(DEMON) 신호와 같이 포락선의 시간지연이 중요한 경우에는 위상 시간지연 대신 군 시간지연의 설계 결과가 검토되어야 한다.^[7]

V. 시뮬레이션

앞장에서 도출한 수중음향 전달 채널의 이산 시간 근사 방법을 블록도로 표현하면 Fig. 8와 같다. 제안 방법은 우선 Bellhop 프로그램을 이용하여 고유음선을 찾은 뒤, 각각의 고유음선을 필터로 근사하는 과정을 거친다. 이후 각각의 고유음선 별로 근사된 필터 응답을 합하면, 전달 채널의 모든 정보가 포함된 유한 임펄스 응답을 얻으며, 임의의 음원의 신호와 유한 임펄스 응답의 콘볼루션을 통해 소나에서 수신한 신호를 모의할 수 있다. 참고로 필터 차수가 긴 경우의 콘볼루션 연산은 신호의 분할과 FFT를 활용한 주파수 영역의 연산으로 효율적으로 계산할 수 있다.^[8]

Fig. 8.

(Color available online) Block diagram for the discrete time modeling of underwater propagation channel.

Fig. 8의 블록도에서 음선의 진폭 $A=R{e}^{j\varphi }$은 주파수의 함수이나, 시간지연 $\tau$는 앞서 II장에서 언급한 바와 같이, 무손실의 아이코날 방정식의 해로부터 주파수와 무관하게 주어진다. 단, Bellhop의 경우 수치분석 단위를 자동으로 지정한 경우, 주파수에 따라 단위가 달라져 시간지연 $\tau$의 연산오차가 발생할 수 있다.^[3] 또한 Bellhop에서는 음선의 진폭 $A=R{e}^{j\varphi }$ 역시, 사용자가 반사계수 테이블을 주파수 별로 제공하지 않는 한 평면파 반사를 가정하므로, 체적 손실을 제외하면 주파수의 함수가 아니다. 따라서 Fig. 8과 같이 고유음선 분석을 주파수별로 실시하지 않고 1회만 실시한 뒤, 필터 설계 과정에서 별도로 주파수별 손실에 따른 크기 및 위상을 보상하는 방법도 가능하다.

제안한 방법으로 음선 추정 결과를 활용하여 모델링한 광대역 수중음향 전달 채널을 Figs. 9와 10에 도시하였다. 수심은 1000 m이고, 음원 깊이는 10 m, 수신기 깊이는 200 m로 두었다. 음원과 수신기는 5 km 이격하였다. 해수면 경계는 진공으로 두었으며, 해저면은 음속이 1600 m/s, 물의 밀도의 1.8배인 음향학적 반공간으로 설정하였다. 가우시안 빔 경로 기법을 적용하였고, 수직각 ±12° 이내에 101개의 음선을 사용하였으며, Eq. (3)의 체적손실을 반영하였다. Fig. 9에서는 가상의 등음속 프로파일을, Fig. 10에서는 동해에서 겨울철에 측정한 실제 음속 프로파일을 적용하였다.

Fig. 9.

(Color available online) Impulse response of underwater acoustic channel with iso-speed profile.

Fig. 10.

(Color available online) Impulse response of underwater acoustic channel with measured sound speed profile.

Figs. 9과 10의 시뮬레이션 결과 그림은 좌측의 음속 프로파일, 중간의 고유 음선도, 우측 상단의 Bellhop 고유음선의 진폭 및 우측 하단의 임펄스 응답 근사 결과로 구성되어 있다. 우측 상단 및 하단 그림의 x축은 음원에서 출발하여 음파가 수신기에 도달하는 도달 시간 지연(단위 : 초)을 나타내며, 우측 하단의 임펄스 응답의 샘플링 주파수는 30 kHz이다. 우측 상단에서는 고유 음선의 진폭에서 위상을 제외한 절대값 만을 도시하였다. 중간의 고유 음선도에서는, 직접 경로는 빨강색, 해수면 반사 경로는 녹색, 해저면 반사 경로는 파랑색, 해수면/해저면에 모두 반사되는 경로는 검은색으로 표시하였다.

고유 음선 분석 결과 Fig. 9의 가상의 등음속 프로파일 환경의 경우, 음원과 수신기 사이에 직접 경로와 해수면 반사 경로 두 가지가 존재하며, 시간이 앞선 음선이 직접 경로가 된다. 두 경로의 크기는 유사하나, 해수면 반사 경로는 위상이 $\pi$만큼 반전된다. 따라서 직접 경로의 임펄스 응답과 그 뒤의 해수면 반사 경로의 임펄스 응답은 부호를 제외하면 유사한 형태를 나타낸다.

Fig. 10은 측정된 음속 프로파일을 적용한 결과이다. 음속 프로파일의 특성상 해수면 부근의 덕트가 형성되어 대부분의 에너지가 층심도 이하로 전파되지 않으나, 일부 해저면 반사에 따른 경로가 존재한다.

한편, Fig. 10의 해저면 반사 경로의 임펄스 응답은 Fig. 9의 직접경로 또는 해수면 반사 경로와 상이한 임펄스 응답을 보이며, 따라서 해저면 반사 경로를 통해 전달되는 신호는, 직접 경로와 다른 왜곡을 거쳐서 수신됨을 알 수 있다. 이는 Fig. 9의 등음속 프로파일의 결과에서처럼 직접 경로와 해수면 반사가 동일한 신호의 시간지연 합으로 이해될 수 있는 상황과 달리, 수중 환경에서 해저면 반사에 따른 신호는 직접 경로와 상관도가 떨어질 수 있음을 의미한다. 실제로도 해상시험을 통해, 여섯 가지 조합의 수중 다중경로의 시간차를 이용하여 표적의 위치추정을 실시한 사례에서, 해저면 반사 신호와 해저면-해수면 반사 신호의 시간차를 이용한 표적의 위치 추정 결과가 가장 좋은 성능을 보였다.^[9]

VI. 결 론

본 논문에서는 음선 경로법에 기반한 광대역 수중음향 전달 채널을 모델링하는 방법을 다루었다. 수중음향 전달 채널은 종종 주파수 영역에서 취급되어, 시간 영역 신호 모의 시 활용이 어렵다. 한편 음선 모델은 고유 음선별 도달 시간 및 진폭을 계산할 수 있어, 다중경로 신호의 합으로부터 광대역 신호를 손쉽게 모의할 수 있으므로, 광대역 신호의 시계열 모의에 적합하다. 다만 고유 음선은 연속 시간 방정식의 해로써 산출되므로, 이를 이산시간 시스템에 적용하기 위해 적절한 근사가 요구된다. 이에 본 논문에서는 고유 음선 분석을 통해 산출된 주파수별 손실 및 위상, 그리고 음선의 도달 시간지연을 이산시간 유한 임펄스 응답으로 근사하는 방법을 제안하였으며, 목표 주파수 대역 내에서는 정확한 근사가 가능함을 확인하였다. 다만, 본 연구 결과는 시불변 음파 전달 채널만 모의가 가능하여, 실제 환경에서 발생하는 도플러 효과를 모의할 수 없으므로, 추후 시변 채널에서 모의가 가능한 방법에 대한연구가 필요하다.