압축 차주파수 빔형성 기법을 이용한 도래각 추정

2.1 차주파수 빔형성기법(FD)

하다마드 곱(z∈ℂ^M)은 식 (3)과 같이 f_l과 f_u(f_l < f_u) 두 주파수 성분에 대한 요소별 곱(elementwise)으로 정의된다^(3,9).

기호 ‘◦’는 벡터 성분에 대한 요소별 곱을, ‘*’는 켤레 복소화(complex conjugation)을 의미한다. 두 주파수 f_u와 f_l은 신호 대역폭 내에 속해야 하며(f_L ≤ f_l < f_u ≤ f_U), 중심 주파수 f_c = (f_u + f_l)/2 및 차주파수 ∆f=f_u-f_l를 이용한 등가표현은 f_L + ∆f/2 ≤ f_c ≤ f_U - ∆f/2, 0 ≤ ∆f ≤ f_U-f_L과 같다.

하다마드 곱은 차주파수에서 음향 헬름홀츠 방정식의 해를 모방(mimic)한다⁽⁹⁾. 이는 진폭 x(f)를 갖는 평면파 y(f,r) = x(f)exp(j2πfe_sㆍr)를 이용하여 증명 가능하다(e_s는 파 진행방향 단위벡터, r은 위치벡터를 의미함). 하다마드 곱은 y(f_u,r)y*(f_l,r) = x(f_u)x*(f_l)exp(j2πfe_sㆍr)와 같은데, 진폭 x(f_u)x*(f_l), 주파수 ∆f를 가지는 평면파와 같다. 부연하면, 하다마드 곱은 그것이 의도하는 음장과 단위/진폭 측면에서 다르더라도, 여전히 저 주파수인 ∆f에서의 평면파를 묘사하고 있는 것이다(이에 상기에서 ‘모방’이라는 용어를 사용하였다).

따라서, CBF가 조향벡터 a(θ,f)와 측정벡터 y사이의 유사성을 평가하듯이, FD의 빔파워(B_FD)를 산출하기 위한 자연스러운 선택은 식 (4)처럼 y대신 z를, f 대신 ∆f에서의 조향벡터 a(θ,∆f)를 사용하는 것이다.

여기서, 상첨자 ‘H’는 허미시안(hermitian) 전치이다. B_FD는 조향각도(look angle) θ가 음원의 DOA와 일치할 때 최대화된다. ∆f가 신호의 주파수 대역보다 훨씬 낮기 때문에, 해당 파장 ∆λ(=c/∆f)이 공간 앨리에이싱 방지조건인 d < ∆λ/2를 충족하는 범위에서 ∆f를 선택한다. 저 주파수 영역 신호처리는 CBF가 취급하기 어려운 희박 배열에 대한 접근을 가능하게 할 뿐만이 아니라, 사용된 음원 모형에 대한 오류를 저감시키는 장점을 제공한다.

하다마드 곱의 특징 중 하나는 참 DOA 외에도 유령(ghost) DOA로 귀결되는 교차항(cross-term)이 발생한다는 것이다⁽⁹⁾. 이 유령 DOA θ_g는 공간 나이키스트(nyquist) 기준을 위배하여 발생하는 앨리에이싱 DOA와 구별된다. K개 음원에 대해 식 (3)의 하다마드 곱은 K²개의 항을 생성하는데, 이 중 K개의 자기항(auto-term)은 참 DOA에, 나머지 K(K-1)개의 교차항은 유령 DOA에 해당한다. 참 DOA와 달리 유령 DOA는 주파수 의존성을 가지므로⁽¹⁰⁾, 식 (4)에 대한 주파수 평균을 통해 유령 DOA의 영향을 억제할 수 있다. 즉, 식 (4)에 대한 광대역 표현은 식 (5)와 같다.

〈 〉_fc,∆f는 두 주파수 f_c,∆f에 대한 이중 평균을 의미하며, 식 (4)와 구분하기 위해 막대기호 ‘－’를 사용하였다. 상기 설명한 FD 기법은 단일 스냅샷(L=1)에 대한 것이다. 신호 대역폭이 넓은 경우(충분한 주파수 평균이 가능하므로), 단일 스냅샷 상황에서도 유령 DOA의 억제가 가능하다. 다중 스냅샷(L>1)의 경우, 스냅샷 평균화를 통해 빔형성 결과의 안정성을 더욱 향상시킬 수 있다.

한편, FD는 처리 주파수 감소에 따라 주엽폭이 넓어져, 서로 인접한 DOA의 식별이 어렵다. 다음절에서는 FD의 해상도 한계를 개선하기 위한 방법인 압축 차주파수 빔형성 기법(CFD)에 대해 논의한다.

2.2 압축 차주파수 빔형성기법(CFD)

압축센싱(compressive sensing, CS)은 정보통신 분야에서 제안된 이론으로서, 해가 희소한 경우 부정 문제(underdetermined problem)를 풀이할 수는 방법이다. 음원(해)은 공간 내에서 작은 범위를 차지하므로, CS를 DOA 추정에 적용할 수 있다. 압축 빔형성 기법으로^(6,8) 알려진 해당 방법은 l₁-norm 기반 최적화 문제⁽⁷⁾, 즉, LASSO를 풀이하여 DOA를 추정한다. 식 (2)의 신호모형에 대한 LASSO 문제는 식 (6)과 같이 정식화된다.

여기서, $$ {‖\mathit{x}‖}_p={\left(\sum _{n=1}^N {\left|x_n\right|}^p\right)}^{1/p}$$는 l₁-norm을 의미한다. 측정행렬 A(f)=[a(θ₁,f)a(θ₂,f)...a(θ_N,f)]∈ℂ^M×N는 주파수 f에서 모든 조향각도에 대한 조향벡터를 수집하여 구성된다. 정규화인자 μ(>0)는 모형오차항(식 (6)의 우변 첫째항)와 미지수 x = [x₁(f) x₂(f) … x_N(f)]^T∈ℂ^N의 희소성(식 (6)의 우변 둘째항) 사이의 상대적인 기여도를 조정한다. μ를 증가시키면 해 x_cs의 희소성을 촉진하고, 반대의 경우 모형적합성을 지향하여 x_cs에서 비영요소(nonzero element)의 수가 증가한다. 따라서, 적절히 선정된 μ에 대해, x_CS가 생성하는 DOA 스펙트럼은 초해상도(super-resolution) 특성을 가지게 된다.

이 연구에서 제안하는 CFD에서는, CS(식 (6))에서 배열데이터 y를 하다마드 곱 z로 대체하고, 차주파수에서 측정행렬을 산출한다, A(∆f) = [a(θ₁,∆f) a(θ₂,∆f) ... a(θ_N,∆f)]. 즉, 식 (7)과 같다.

CFD는 CS와 달리 2차 오더(order) 하다마드 곱을 사용하여 정의되므로 1차 오더 x와 구분하기 위해 2차 오더 미지수 기호 q를 사용한다. CS가 CBF의 희소화 표현인 것과 마찬가지로, CFD는 FD의 희소화 표현으로 간주할 수 있다. 식 (6)과 식 (7)을 비교하면, 측정벡터와 산출주파수만이 변경됨을 알 수 있다. 이로써 CFD는 CS의 장점을 계승하여 우수한 공간 해상도 뿐만 아니라 상관음원/단일 스냅샷 상황에도 적용가능하다는 장점을 갖는다. 그러나, 기저 불일치(basis mismatch) 및 상관기저(coherent basis)와 관련된 CS의 단점은⁽⁸⁾ CFD 역시 회피할 수 없다.

식 (7)은 컨벡스(convex) 최적화 문제로서, 공개 소프트웨어인 CVX⁽¹¹⁾ 또는 l₁-magic 등을 이용하여 풀이할 수 있다⁽¹²⁾. 충분히 수렴한 경우, 해 q_CFD는 K²개의 비영요소를 가지며(K²≪N), 나머지 N-K²개 요소들은 0 또는 0에 가까운 값을 가진다.

CFD의 빔형성 파워는 q_CFD의 제곱으로 표현되며(B_CFD(θ)=|q_CFD|²), 식 (4)의 B_FD와 견줄 수 있다. CFD 역시 교차항이 유발하는 유령 DOA를 억제하기 위해 주파수 평균을 수행해야하며, 이에 식 (5)에 상응하는 CFD 표현은 식 (8)과 같다.

FD와 동일하게, 식 (8)은 단일 스냅샷에 대한 결과이다. 당연히, 다중 스냅샷에 대해서는 추가적인 스냅샷 평균을 수행하여 추정결과의 안정성을 향상시킬 수 있다.

3.1 시뮬레이션

Fig. 1은 FD 및 CFD의 성능비교를 위한 환경설정을 보이고 있다. 자유 수중 공간내(c = 1510 m/s) M = 16개의 센서를 d = 3.75 m 간격으로 배치하였다. 이 배열에 대해, 공간 앨리에이싱을 회피하기 위한 공간 나이키스트 주파수는 f_N(=c/(2d)) = 200 Hz이다. 두 개 평면파 음원 θ₁, θ₂를 가정하되, θ₁은 0°에 고정되어 있으며 θ₂는 0° ~ 25° 범위에서 변화시키며 제안기법의 음원분리능력을 평가하고자 한다. 음원의 주파수 대역은 [f_L, f_U] = [10, 20] kHz이며, 해당 범위에서 d/λ는 25 ~ 50로써 기존방법(CBF, CS)을 DOA 추정에 적용하기는 어렵다. L = 20개 스냅샷을 식 (2)를 이용하여 생성하였으며, 신호대 잡음비(signal-to-noise ratio, SNR = $$ 20{log}_{10}\left(E\left[{∥\sum _{k=1}^K\mathbf{a}\left({\theta }_k\right)x_k ∥}_2\right]/E\left[{∥\mathbf{n}∥}_2\right]\right)$$)는 30 dB로 설정하였다.

조향각도 탐색범위 및 간격은 각각 θ∈[-90°, 90°], 0.1°이다(N = 1801). 공간 앨리에이싱을 회피함과 동시에 FD 주엽폭을 가급적 줄이기 위해 배열의 공간 나이키스트 주파수(f_N = 200 Hz)를 차주파수 ∆f로 사용한다. ∆f를 일정하게 유지함에 따라, 중심주파수 f_c에 대해서만 평균화를 수행하며(식 (5), 식 (8)), 식 (9)와 같은 겹침(overlapping) 규칙을 이용하여 f_c(=(f_l + f_u)/2)를 변화시켰다.

식 (9)에서, F=1+(f_U-f_L-∆f)/[(1-q)∆f]는 (f_l, f_u) 주파수 쌍 개수를 의미하며, 겹침비율 q를 75 %로 설정함에 따라 F = 197, 즉 197회의 주파수 평균을 취하게 된다. 이 이상의 평균은 결과에 큰 영향을 주지 않음을 확인하였다.

식 (7)의 컨벡스 최적화 문제 풀이를 위해 CVX를 이용하였으며⁽¹¹⁾, L-curve 방법을 이용하여 정규화인자 μ를 0.1로 선정하였다⁽¹³⁾. 엄밀하게는, 각 주파수 쌍 마다 μ를 세밀히 조율해야 하지만, 편의상 그 값을 일정하게 유지하였다. 선정치 기준 수십 % 범위에서 변화시키더라도, 아래 설명할 시뮬레이션 결과와 양상이 크게 다르지 않다.

첫 시뮬레이션에서는 비상관 음원을 다룬다. 즉, 두 음원간 상대위상은 균일분포 ∪[0, 2𝜋]을 따른다. Fig. 2는 세가지 대표적인 경우에 대한 FD와 CFD의 비교를 보이고 있다. 왼쪽 패널은 스냅샷에 대해 평균한 빔형성 파워 B를 중심주파수 f_c 및 조향각도 θ에 대해 표현한 2D 결과이다. FD의 경우 f_c에 대해 일정하게 나타나는 참 DOA의 주변에서, 동적으로 변화하는 유령 DOA를 볼 수 있다. 중간 및 오른쪽 패널에서 보듯, f_c축 주파수 평균을 통해 얻은 DOA 스펙트럼은 유령 DOA를 억제하여 참 DOA 성분만 남긴다. CFD에서는 2D 결과의 희소성이 향상됨에 따라, DOA 스펙트럼이 높은 해상도를 갖는 것을 확인할 수 있다.

Fig. 2

Comparison of FD and CFD for DOA differences between two sources. θ1 is at 0°, while another changes(left: beampower B using Eq. (4) for FD and BCFD (θ) = |qCFD|2 for CFD as a function of the look angle θ and center frequency fc, mid: frequency averaged-beampower B using Eq. (5) for FD and Eq. (8) for CFD, right: detailed view for the region [-10°, 20°] in the mid. True DOAs are represented by dotted lines).

FD에 대한 rayleigh 해상도 한계는 ∆θ_FD = sin^-1 (∆λ/Md) = 7.2°로서⁽²⁾, 두 음원이 이 각도 이상 분리되어야 음원식별이 가능함을 의미한다. Fig. 2(a)와 Fig. 2(b)는 분리각도가 rayleigh 한계보다 큰 경우, Fig. 2(c) ~ Fig. 2(f)는 반대 경우에 속한다. 넓은 분리간격에서는 FD, CFD 두 방법 모두 DOA를 정확하게 추정하는데, CFD가 FD 보다 고해상도 추정결과를 제공함을 알 수 있다. 분리각도가 ​∆θ_FD이하(Fig. 2(c))이면 FD의 두 주엽이 합쳐져 마치 하나의 음원만이 존재하는 것으로 해석될 수 있는 오해를 야기한다. 반면, Fig. 2(d)에서 볼 수 있듯이, CFD에서는 두 DOA를 정확하게 식별할 수 있음을 알 수 있는데, 이는 제안 기법의 우수성을 입증하고 있다.

그러나, 분리각도가 매우 작은 경우, FD 뿐만이 아니라 CFD도 DOA 분리가 어렵다(Fig. 2(e), Fig. 2(f)). Fig. 3은 시뮬레이션에서 고려된 모든 DOA에 대한 추정오차를 보이고 있다. FD는 분리각도가 ∆θ_FD = 7.2°보다 커야 두 음원을 구별하기 시작하지만, CFD는 구분 가능한 최소분리각도(분해능)가 2.2°로서 FD보다 우수한 성능을 보인다.

Fig. 3

Comparison of DOA estimation(left) and estimation error(right) for the range of test DOAs

물론, CFD의 분해능 ∆θ_CFD는 FD와는 달리 SNR에 대해 의존적인 특성을 보이며, 이에 추가 시뮬레이션을 통해 SNR에 따른 ∆θ_CFD를 조사하였다. Fig. 4에 도시한 바와 같이, SNR이 감소함에 따라 CFD의 분해능도 감소하는 경향을 볼 수 있는데, 이는 잡음이 커질수록 참 DOA에 해당하는 조향벡터 식별이 모호해지는 상관기저 특성과 관련이 깊다⁽⁸⁾. CFD의 분해능은 FD보다 SNR = 4 dB 이하에서 열등하지만, 그 이상의 SNR 조건에서는 CFD의 분해능이 우수함을 강조할 수 있다.

Fig. 4

Analysis of CFD resolution limit(∆θCFD) with increasing SNR. Dotted line represents the rayleigh resolution limit

다음 시뮬레이션은 다중경로(multi-path) 상황에서 발생할 수 있는 상관음원에 대한 것으로서, 음원간 상대 위상을 0°로 설정하였다. 음원 간 완벽한 상관관계는 분리각도에 따라 두 주엽의 주기적인 상쇄/보강 간섭을 유발한다. 보강간섭의 경우, FD의 두 주엽이 서로를 밀어내어 추정된 DOA가 참 DOA 외부에 위치하고, 상쇄간섭의 경우에는 그 반대가 발생한다. 따라서, FD의 추정오차는 주기적으로 진동하는 거동을 보인다(Fig. 5(a)). CFD의 해상도 한계는 음원의 상관관계와 관련이 없다. 따라서 CFD의 추정오차는 진동을 보이지 않으며, 분리각도 ∆θ_CFD 이상 영역에서 FD 대비 추정도가 향상된다(Fig. 5(b)).

Fig. 5

Comparison of DOA estimation(left) and estimation error(right) for the range of test DOAs (Coherent sources, SNR = 30 dB)

3.2 FAF06 실험

CFD의 성능을 실험적으로 검증하기 위해, FAF06 (focused acoustic field 06)로⁽¹⁴⁾ 지칭되는 수심 100 m 가량 실해역에서 측정한 데이터를 이용한다. 이 실험은 시뮬레이션과 동일한 형태의 수직 ULA를 사용한다. 총 길이 56.25 m인 배열의 중심은 수면으로부터 깊이 52.8 m에 있고, 수면에 가까운 몇 센서를 제외하고는 모두 음속 c = 1510 m/s로 간주할 수 있는 등음속(iso-soundspeed) 깊이 범위내에 있다. 신호 샘플링 주기는 1/50 000 s(= 50 kHz)이다.

FAF06 실험의 주 목적은 배열로부터 멀리 떨어진 음향 프로젝터에서 송출된 신호를 기록하는 것이다. 그러나, 의도된 송출신호가 아닌 수중생물의 먹이활동 또는 배열 케이블과 고정용 부이 간 마찰에서 기인한 충격성 잡음이 자주 포착된다. 이 연구에서는 그 충격 사건들을 합성하여 두 가지 경우를 시연하고자 한다.

Fig. 6(a)에는 첫 시연을 위해 합성된 신호를 보였다. 충격신호 형태를 부각시키기 위해 9 kHz의 고역통과필터(high-pass filter)를 적용하였다. 구형파면(spherical wavefront)을 닮은 선행신호는 겉보기 DOA +1.1° ~ +3.3°를, 후속 평면파 신호는 겉보기 DOA -17°를 가진다. 따라서, 첫 시연은 넓은 분리각도(약 20°)에 대한 것으로 볼 수 있다. 참고로, 별도의 빔형성기법을 거치지 않고 Fig. 6(a)와 같은 워터폴 플롯(waterfall plot)만으로도 DOA 추정이 가능하다. 따라서, 합성된 시계열 데이터(Fig. 6(a) 및 Fig. 7(a))는 FD와 CFD의 성능을 비교하기 위한 용도로만 사용되어야 함에 유의해야 한다.

Fig. 6

Experimental verification for two DOAs with a wide separation

Fig. 7

Experimental verification for two DOAs with a narrow separation

시계열 신호에 대한 2048개 FFT를 통해 주파수 대역 신호를 마련하였고, 조향각도 간격을 제외하고(θ = [-90° : 0.5 : 90°], N = 361) 나머지 해석 설정은 시뮬레이션과 동일하다.

Fig. 6(b)의 FD의 2D 결과에서, f_c를 따라 두 개 일정한 DOA 성분이 나타남을 볼 수 있다. 그 두 성분은 오른쪽의 주파수 평균 결과에 나타난 것처럼 겉보기 DOA와 일치한다. Fig. 6(b)의 CFD의 2D 결과 역시 희소화 처리에 의해 가느다란 두 개의 줄기를 보이며, Fig. 6(c)와 같이 해당 DOA 스펙트럼은 부엽수준을 억제함과 동시에 두 개의 피크만을 보인다. 참고로, 희소결과의 가시성 향상을 위해 Fig. 6(c) 왼쪽의 B_CFD와 3×3- ‘1’ 행렬(3 × 3 ones matrix)의 합성곱 적분(convolution integral) 결과를 Fig. 6(d)에 나타내었다.

두 번째 시연의 겉보기 DOA는 Fig. 7(a)와 같이 -12°, -17°로서 분리각도가 작은 상황을 다룬다. Fig. 7(b)에서 보듯, FD의 2D 결과는 Fig. 6(b) 대비 유의미한 특징이 없다. 작은 분리각도로 인해 주엽들이 서로 무작위로 보강/상쇄 간섭을 반복하고 있으며, 이에 FD의 DOA 스펙트럼은 -13° 주변 하나의 넓은 주엽을 형성하며 동시에 높은 부엽도 생성한다. 이에 비해, CFD는 개선된 추정결과를 보인다. Fig. 7(c)와 Fig. 7(d)와 같이 f_c를 따라 겉보기 DOA에 해당하는 두 개의 미세한 줄기를 생성하고 있으며, 결국 CFD의 DOA 스펙트럼은 서로 인접한 두 DOA 성분을 구분하게 된다.