I. 서 론
음원 위치 추정 분야는 레이더, 소나, 통신, 지질학등 많은 분야에서 중요하게 다루어지는 분야이다.[1] 이중 많은 음원 위치 추정 문제들이 음원을 공간상의 점으로 모델링할 수 있다고 가정한다. 하지만 음원 신호가 다중경로와 같은 분산 환경을 거쳐 센서에 도달할 경우 반사와 산란으로 인해 점음원보다는 분산음원으로 모델링 되는 것이 적절하며 이러한 모델링 방법에 기반한 위치 추정 기법들이 다수의 논문에서 다루어 진 바 있다.[2-6] 분산음원 모델은 크게 음원의 크기와 위상만 다른 동일 신호의 합으로 보는 코히어런트 분산 모델[2,3]과 서로 상관관계가 없는 신호의 합으로 보는 인코히어런트 분산 모델[3,4]로 분류된다. 본 논문에서는 음원주위로 산란된 신호들이 채널의 변화주기보다 짧기 때문에 분산 신호를 크기와 위상만 다른 동일한 신호의 합으로 보는 코히어런트 분산 환경을 가정한다.
음원 위치 추정 문제에 있어 센서 배열의 종류를 보면 구조가 간단하고 해석이 쉬운 균일 선배열 센서가 주로 다루어졌으나 수평방향의 방위각만 추정가능하며 좌우 구분이 불가능하다. 반면 환형 배열 센서의 경우 방위와 고각을 추정할 수 있고 360° 전방위를 동일한 성능으로 탐색할 수 있는 장점이 있다. 환형 배열을 이용한 분산음원의 모델링과 위치 추정 문제는 다수의 논문에서 다루어졌다.[3-6] 그중 Valaee의 DSPE(Distributed Source Parameter Estimator)기법은 MUSIC(Multiple Signal Classification)기법을 확장한 대표적인 기법으로 다수의 음원에 대해 대표 방위와 분포를 동시에 추정할 수 있으나 분포 모델이 사전에 주어져야 하며 환형배열에서는 방위와 고각, 방위와 고각의 분산에 대해 4차원 탐색이 필요하여 매우 많은 연산량을 요구한다.[3] Lee의 SOS(sequential one-dimensional searching) 기법의 경우 대표 방위와 고각을 두 번의 1차원 탐색만으로 추정하여 상대적으로 적은 연산량을 요구하나 균일한 환형 배열센서 한 쌍을 필요로 하는 제한점이 있다.[5]
최근 균일 환형 배열 점음원 모델에서 상관함수와 최소 자승기법을 이용하여 음원의 위치를 추정하는 기법이 제안되었다.[7,8] Liao의 기법은 음원이 한 개라는 가정하에서 탐색과정 없이 음원의 방위와 고각을 닫힌 형태(closed form)으로 추정하여 적은 연산량을 요구한다.[7] 본 논문에서는 이를 코히어런트 분산 음원에 적용하여 탐색과정 없이 음원의 대표 방위와 대표 고각을 추정하는 기법을 제안한다. 특히 음원의 분산 모델이 가우시안 모델인 경우 비슷한 과정을 통해 방위와 고각의 분산 또한 탐색과정 없이 추정하여 매우 빠르게 표적의 모든 파라미터 추정이 가능하다.
II. 신호 모델
Fig. 1은 균일 환형 배열에서 M개의 전 방향성(omnidirectional) 센서가 단일 협대역 코히어런트 분산 음원의 신호를 수신하는 것을 나타낸다. M개의 센서는 xy평면에 반시계 반향으로 균등하게 분포하고 있으며 센서의 중심을 기준으로 r의 반지름을 가진다. 이때 센서배열에 수신된 신호 x(t)는 다음과 같다.[5]
| $$\chi(t)=\int_{-\pi/2}^{\pi/2}\int_0^{2\pi}a(\vartheta,\varphi)\zeta(\vartheta,\varphi,t)d\vartheta d\varphi+n(t).$$ | (1) |
는 조향 벡터, 는 음원의 각밀도함수(angular spread density function), n(t)는 백색 가우시안 잡음을 의미한다. 여기서 는 파장이며 이다. 는 아래 수식으로 표현된다.
| $$\zeta(\vartheta,\varphi,t)=s(t)\rho(\vartheta,\varphi;\mu),$$ | (3) |
여기서 s(t)는 음원 신호, 는 결정(deterministic) 각밀도함수이다. 는 대표 방위, 대표 방위의 분포, 대표 고각, 대표 고각의 분포를 의미한다. 조향 방위에 대해 수신 신호를 표현 할 경우 코히어런트 분산 음원은 아래와 같이 표현된다.
| $$x(t)=s(t)b(\mu)+n(t).$$ | (4) |
| $$b(\mu)=\int_{-\pi/2}^{\pi/2}\int_0^{2\pi}a(\vartheta,\varphi;\mu)d\vartheta d\varphi.$$ | (5) |
Reference [5]에 따르면 방위분포()와 고각분포()의 분포가 충분히 작은 값을 가질 경우 b(µ)는 다음과 같이 표현된다.
| $$b(\mu)=a(\theta,\phi)\otimes g(\mu),$$ | (6) |
여기서 \otimes는 Hadamard product이다. g(µ)는 결정 각분포함수로 가장 많이 이용되는 가우시안(Gaussian) 분포일 경우 아래와 같다.
III. 제안 기법
환형배열에서 임의의 센서 p,q에 대해 다음과 같은 상관함수를 정의한다.
여기서 는 켤레 복소수, 과 은 음원과 배경 잡음의 파워이며 배경 잡음은 백색잡음을 가정한다. 위 수식에서 상관함수의 위상은 대표방위와 대표고각의 함수로 나타나며 크기 성분은 신호의 파워와 결정 각밀도 함수로 나타낸다.
3.1 대표방위, 대표 고각 추정 기법
수신된 신호에서 잡음성분을 배제하고 신호 성분만 고려하면 Rp,q의 위상은 아래와 같다.
| $$\begin{array}{l}u_{p,q}=arg(R_{p,q})\\\;\;\;\;\;\;\;=\eta(\sin\phi\cos(\theta-\gamma_p)-\sin\phi\cos(\theta-\gamma_q))\end{array}.$$ | (9) |
Eq. (9)는 다음과 같이 벡터의 곱으로 정리된다.
| $$b(\mu)=\int_{-\pi/2}^{\pi/2}\int_0^{2\pi}a(\vartheta,\varphi;\mu)d\vartheta d\varphi.$$ | (10) |
p–q = l로 두고 다수의 센서 쌍에 적용하여 다시 정리하면 다음과 같이 행렬 연산으로 표현가능하다.
| $$u=Ad.$$ | (11) |
| $$u=\lbrack u_{1,1+l}\;u_{2,2+l}\;\dots\;u_{M-l,M}\rbrack^T.$$ | (12) |
| $$d=2\eta\begin{bmatrix}\cos\theta\sin\phi\\\sin\theta\sin\phi\end{bmatrix}.$$ | (14) |
𝜃와 𝜙를 추정하기 위해 u를 추정치 로 대체하면 는 다음의 최소자승법으로 표현된다.
| $$\widehat\theta=arg({\widehat d}_1+j{\widehat d}_2).$$ | (15) |
따라서 다음과 같이 음원의 방위와 고각을 구한다.
| $$\widehat\theta=arg({\widehat d}_1+j{\widehat d}_2).$$ | (16) |
| $$b(\mu)=\int_{-\pi/2}^{\pi/2}\int_0^{2\pi}a(\vartheta,\varphi;\mu)d\vartheta d\varphi.$$ | (17) |
Eqs. (16)과 (17)을 통해 대표방위와 대표 고각을 추정할 수 있으며 분산음원의 분포 종류와 상관없이 사용 가능하다.
3.2 방위 분포와 고각 분포 추정
분산 음원의 분포함수가 주어진 경우 다음 식과 같이 MUSIC 기법을 일반화한 DSPE 기법을 통해 음원의 분포값을 추정할 수 있다.
| $$\{\widehat\theta,\widehat{\sigma}_{\theta},\widehat\phi,\widehat{\sigma}_{\phi}\}=argmax\frac1{\parallel\alpha(\theta,\sigma_\theta,\phi,\sigma_\phi)V\parallel^2},$$ | (18) |
여기서 V는 잡음 공간을 의미한다. Eq. (18)에서 모든 파라미터를 한번에 추정하기 위해서는 4차원 탐색 과정이 요구되며 3.1절에서와 같이 대표방위와 고각이 주어진 경우에는 방위 분포와 고각 분포에 대해 2차원 탐색이 요구된다. 그 중 결정 각함수가 Eq. (7)의 가우시안 분포를 따를 경우 다음 과정을 통해 닫힌 형태로 구할 수 있다. Eq. (8)에 절대값과 자연로그를 적용하고 신호성분만 고려 할 경우 다음과 같다.
| $$\ln\left|R_{p,q}\right|=2\ln\left|\sigma_S\right|+\ln\left|g_p(\mu)\right|+\ln\left|g_q(\mu)\right|.$$ | (19) |
여기서 과 는 다음과 같다.
| $$\begin{array}{l}f_1(\theta,\phi,\gamma_p,\gamma_q)\\=\frac{-1}{2(\eta\cos\phi\cos(\theta-\gamma_p))^2}+\frac{-1}{2(\eta\cos\phi\cos(\theta-\gamma_q))^2}.\end{array}$$ | (21) |
| $$\begin{array}{l}f_2(\theta,\phi,\gamma_p,\gamma_q)\\=\frac{-1}{2(\eta\sin\phi\sin(\theta-\gamma_p))^2}+\frac{-1}{2(\eta\sin\phi\sin(\theta-\gamma_q))^2}.\end{array}$$ | (22) |
다수의 센서 쌍에 대해 다시 정리하면 다음과 같이 행렬 연산으로 표현된다.
| $$\upsilon=\left[\upsilon_{1,1+l}\;\;\upsilon_{2,2+l}\;\;\dots\;\upsilon_{M-l,M}\right].$$ | (23) |
| $$f=\begin{bmatrix}\ln\left|\sigma_s\right|\\\sigma_\phi^2\\\sigma_\theta^2\end{bmatrix}.$$ | (25) |
Eq. (24) 행렬 C의 변수 중 𝜃, 𝜙는 3.1절에서 추정한 로 대체하고 를 구하기 위해 다음의 최소자승기법을 적용한다.
| $$\widehat f=\begin{bmatrix}f_1\\f_2\\f_3\end{bmatrix}=(C^TC)^{-1}C^T\widehat\upsilon.$$ | (26) |
따라서 방위의 표준편차 와 고각의 표준편차 는 다음식과 같이 f2와 f3의 제곱근이다.
| $$\sigma_\theta=\sqrt{f_2}.$$ | (27) |
| $$\sigma_\phi=\sqrt{f_3}.$$ | (28) |
연산량의 경우 제안기법은 Reference [8]과 동일하며 DSPE기법의 경우 MUSIC 기법과 같은 연산량을 요구한다.
IV. 모의 실험
제안 기법의 타당성을 비교검증하기 위해 가장 대표적인 기법인 DSPE기법과 함께 모의 실험을 수행하였다. 음원의 파라미터는 =[30°, 50°, 2°, 3°]이다. 배열의 센서 수 M=10이며 배열의 반경 이다. 그리고 제안 기법은 l = 5를 이용하였다.
Fig. 2는 신호대 잡음비(Signal to Noise Ratio, SNR)에 따른 추정성능을 평균 제곱근 오차(Root Mean Square Error, RMSE)로 나타내었다. 시험은 100회의 Monte Carlo 시험을 수행하였으며 이때 사용한 스냅샷의 개수는 N = 100이다. 제안 기법은 적은 연산량에도 불구하고 모든 SNR에서 DSPE기법과 성능의 차이가 없다. 모의 실험에서 두 기법 모두 대표방위/고각 추정 성능에 비해 방위/고각 표준편차의 추정 성능이 떨어짐을 보였는데 이를 DSPE 탐색을 통해 Fig. 3에서 분석하였다. Fig. 3(a)는 대표 방위/고각에 대한 2차원 탐색 결과를 RMSE 등고선으로 표현 하였으며 Fig. 3(b)는 방위/고각 표준편차에 대한 결과를 RMSE 등고선으로 표현하였다. 이때 SNR = 10 dB이며 그 외 변수는 정확한 값을 대입하였다. 그림에서 대표 방위/고각 추정의 등고선보다 표준편차 추정의 등고선이 완만하게 분포해 있으며, 이로 인해 표준편차 추정 값이 잡음에 더 민감함을 확인하였다.
마지막으로 시뮬레이션을 이용하여 연산량 수치를 비교하였다. DSPE 기법의 경우 탐색 범위에 따라 연산량이 달라지므로 4가지 변수에 대해 모두 51번의 탐색과정을 고려했다. 모의실험은 i5 3.4 GHz PC에서 50번 시험을 평균하였으며 Matlab의 'tic', 'toc' 명령어로 측정되었다. Table 1에서 제안기법은 DSPE 기법과 비교하였을 때 매우 큰 연산시간의 차이를 확인하였다.
Table 1. Elapsed CPU times required for calculation.
| Proposed method | DSPE | DSPE/Proposed method | |
| Time (s) | 9.9578e-05 | 83.1497 | 8.3502e+05 |
