Improvement of non-negative matrix factorization-based active sonar reverberation suppression method using L1-norm and majorization-minimization

Seokjin Lee; Wonnyoung Lee; Yena You; Seungheon Lee; Daekyung Kim; Junsub Nam

doi:10.7776/ASK.2025.44.2.094

Preview

Research Article

The Journal of the Acoustical Society of Korea. 31 March 2025. 94-108
https://doi.org/10.7776/ASK.2025.44.2.094

Improvement of non-negative matrix factorization-based active sonar reverberation suppression method using L1-norm and majorization-minimization

L1-노름과 주요화최소화 기법을 활용한 비음수 행렬 분해 기반 능동소나 잔향제거 개선 방법

Seokjin Lee¹^*

Wonnyoung Lee¹

Yena You¹

Seungheon Lee¹

Daekyung Kim²

Junsub Nam²

이 석진¹^*

이 원녕¹

유 예나¹

이 승헌¹

김 대경²

남 준섭²

¹경북대학교 전자전기공학부

²한화시스템 해양연구소

^{*Corresponding Author}

License (open-access, http://creativecommons.org/licenses/by-nc/4.0/):

This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.

ABSTRACT

In this paper, we conduct research on an algorithm to suppress reverberation using a non-negative matrix factorization technique for target detection using an active sonar system. In this paper, we focus on the disadvantage that conventional algorithm estimates the target echo and reverberation bases through iterative estimation, which takes a long time to calculate. To improve this, we use the L1-norm to quickly converge to the desired solution at the beginning of the iteration. In addition, considering the fact that the heuristic multiplicative update method used in the conventional algorithm does not mathematically guarantee the convergence of additional penalty functions, the update equation is developed through a majorization-minimization technique that can guarantee convergence. Simulations were performed to verify the algorithm, and the results showed that the signal-to-noise ratio performance of the proposed algorithm was improved by 3 dB compared to the conventional algorithm.

Keywords

Active sonar

Reverberation suppression

Non-negative matrix factorization

Majorization-minimization

본 논문에서는 능동 소나 시스템을 이용한 표적 탐지를 위해 비음수 행렬 분해 기법을 활용하여 잔향을 제거하는 기법에 대한 연구를 진행하였다. 본 논문에서는 기존의 기법이 반복 추정을 통해 표적 및 잔향 기저를 추정하여 연산 시간이 오래 걸린다는 단점에 주목하였다. 이를 개선하기 위하여 L1-노름을 활용하여 반복 초기에 원하는 해로 빠르게 수렴할 수 있도록 하였다. 또한, 기존의 기법이 활용하고 있는 휴리스틱 곱셈 갱신 기법은 추가적인 비용 함수를 활용할 경우 수학적으로 수렴을 보장하지 못하기 때문에, 수렴을 보장할 수 있는 주요화최소화 기법을 통해 갱신식을 새롭게 도출하였다. 본 알고리즘을 검증하기 위하여 시뮬레이션을 수행하였으며, 제안하는 알고리즘이 기존 대비 약 3 dB 향상된 신호대잔향비 성능을 보였다.

키워드

능동 소나

잔향 제거

비음수 행렬 분해

주요화최소화

MAIN

I. 서 론
II. 비음수 행렬 분해 기반의 기존 잔향 제거 알고리즘
2.1 비음수 행렬 분해와 휴리스틱 곱셈 갱신
2.2 비음수 행렬 분해 기반의 능동 소나 잔향 제거 기법
III. 적은 반복 횟수에 대한 성능 개선 방안
3.1 L1-노름을 이용한 추가 제약 조건
3.2 주요화최소화 기법을 활용한 갱신 기법
IV. 시뮬레이션
4.1 시뮬레이션 환경
4.2 결과
V. 결 론
부 록

I. 서 론

수중에서 표적을 탐지하기 위해서는 음향 신호를 활용한 소나 시스템을 활용한다. 특히 표적의 위치 정보 등을 효과적으로 탐지하기 위해서는 소나 시스템이 특정한 음향 신호를 방사하여 그 반향을 탐지하는 능동 소나 시스템을 활용하는 것이 효과적이다.

표적이 발생시키는 소음을 수신하여 이를 분석하는 수동 소나 시스템과 달리, 능동 소나 시스템은 반향하는 신호의 세기를 조절하여 탐지 거리를 확보할 수 있다. 송신 음파를 낮은 세기로 방사하는 경우, 탐지 확률은 주로 신호대잡음비에 비례한다. 따라서 음향 전파 거리에 따라 크기가 감쇠된 표적 반향음이 충분한 신호대잡음비를 확보하도록 송신 파형의 세기를 크게 하면 먼 거리의 표적도 탐지가 가능하다.^[1]

불행히도, 이와 같은 방법이 항상 적용되는 것은 아니다. 방사되는 음향 신호의 세기가 커짐에 따라 주변 산란체로부터의 잔향도 함께 커지게 되며, 잔향음의 세기가 주변 잡음의 세기보다 커지게 되면 탐지 성능은 신호대잡음비가 아닌 신호대잔향비에 의해 결정이 된다.^[1] 잔향의 세기는 방사 음향 신호의 세기와 비례하기 때문에, 이 경우에는 단순히 방사 음향의 신호를 강하게 하는 것으로는 탐지 성능을 확보할 수 없다.

잔향이 강한 환경에서 표적 반향을 탐지하기 위해서 여러 연구가 진행되었으며, 그 중 한 가지는 잔향에 강인한 송신 파형을 설계하여 잔향의 영향을 줄이는 방법이다. 기존의 능동소나 시스템에서는 단일 주파수 성분으로 구성되어 있어 도플러 탐지에 유리한 지속파 신호를 송신 파형으로 활용하고 있었는데, 이는 잔향에 취약한 단점이 있다. 이를 대신하여 잔향에 강인한 선형 주파수 변조 파형^[2] 및 정현파 주파수 변조 파형^[3] 등의 신호가 고안되었고, 최근에는 이를 발전시켜 일반화된 정현파 주파수 변조 파형 등의 신호가 개발되었다.^[4] 선형 주파수 변조 파형의 경우 능동 소나 시스템에서 널리 활용되고 있지만, 표적의 도플러 탐지가 어려운 등의 단점을 가지고 있어 기존의 지속파 파형을 완전히 대체하지 못하였다. 현재 활용되고 있는 능동 소나 시스템은 주로 지속파 신호와 선형 주파수 변조 파형을 함께 활용하는 형태로 운용된다.

이와 같이 능동 소나에서 지속파 신호가 유용하게 사용되기 때문에, 이를 효과적으로 활용하기 위해서 지속파 신호의 잔향을 제거하는 기법에 대해서도 연구가 지속적으로 수행되고 있으며, 주로 신호의 특성을 분석하는 신호처리 기법을 활용하는 방법으로 연구되었다. 확률 모델을 기반으로한 신호처리 기법이 발전함에 따라 Kay와 Salisbury^[5]는 잔향을 제거하기 위해 자기 회귀 모델에 기초한 백색화 기법을 제안하였으며, 주성분 분석 등의 신호 분석이 발전된 이후 Ginolhac과 Jourdain^[6]은 주성분 역산을 활용한 잔향 제거 기법을 개발하기도 하였다. 이는 추후 신호 부공간 추출의 기법으로 발전하였다.^[7]

최근에는 신호의 의미론적 분석 기법 중 하나로 비음수 행렬 분해(Non-negative Matrix Factorization, NMF) 기법이 제안되었다.^[8,9] 이는 하나의 비음수 행렬을 두 비음수 행렬의 곱으로 분해하는 것으로, 음향 신호의 스펙트로그램을 적은 개수의 음향 이벤트에 대한 주파수 기저와 시간 기저로 분석할 수 있기 때문에 음악 자동 채보^[10] 혹은 음성 잡음 제거^[11] 등의 기법으로 응용되었다.

앞서 언급한 바와 같이 비음수 행렬 분해 기법은 여러 음향 이벤트를 분리할 수 있고, 능동 소나 시스템에서 수신된 신호 또한 표적으로부터의 반향음과 산란체들로부터의 여러 반사음과 같은 음향 이벤트의 합으로 모델링될 수 있다. 이러한 점에 착안하여, 비음수 행렬 분해 기법에 기반한 능동 소나 잔향 제거 기법이 제안되었다.^[12] 비음수 행렬 분해 기반의 잔향 제거 기법은 지속파 신호에 대해 먼저 개발되었고, 이후 선형주파수변조 신호에 대해 적용될 수 있도록 확장되었다.^[13] 비음수 행렬 분해 기반의 잔향 제거 기법은 기존의 주성분 분석이나 신호 부공간 추출 등의 기법에 비해 좋은 성능을 보여주었다. 다만, 비음수 행렬 분해 기법은 하나의 정답을 추정하기 위하여 추정 알고리즘을 반복 수행해야하기 때문에 연산 시간이 오래 걸리는 문제가 있으며, 연산 시간을 줄이기 위하여 반복 시행 횟수를 줄이는 경우 성능이 저하되는 문제가 있다.

본 논문에서는, 비음수 행렬 분해 기반의 잔향 제거 알고리즘에서 기저 행렬 추정의 반복 시행 횟수가 적은 환경에서 성능을 개선하는 방법에 대해 연구를 진행하였다. 본 연구에서 제안하는 성능 개선 시법은, 반복 시행 횟수가 작아도 정답으로 빠르게 수렴할 수 있도록 추가적인 L1-노름 제약 조건을 활용하고, 이를 최적화하는 과정에서 기존의 휴리스틱 곱셈 갱신 대신 주요화최소화 기법을 활용하였다. 일반적으로 표적 반향음의 시간 기저는 희박성을 가지며, L1-노름 제약조건을 활용하면 이러한 특성을 활용하여 추정과정 초기부터 정답으로 빠르게 수렴하는 것을 도울 수 있다. 또한, 기존의 비음수 행렬 분해 기법에서 활용한 휴리스틱 곱셈 갱신 기법은 추가적인 제약조건이 있는 경우 수학적으로 수렴을 보장할 수 없기 때문에 수렴이 느릴 수 있다는 단점이 있으나, 주요화최소화 기법을 활용하면 매 갱신마다 비용함수가 증가하지 않음을 보장하기 때문에 정답으로 더욱 빠르게 수렴될 수 있다.

II. 비음수 행렬 분해 기반의 기존 잔향 제거 알고리즘

2.1 비음수 행렬 분해와 휴리스틱 곱셈 갱신

Lee와 Seung^[8,9]은 다음과 같이 하나의 비음수 행렬 $V \in ℝ_{+}^{K \times N}$ 을 두 비음수 행렬 $W \in ℝ_{+}^{K \times R}$ 및 $H \in ℝ_{+}^{K \times R}$ 의 곱으로 분해하는 비음수 행렬 분해 기법을 고안하였다.

(1)

V \approx WH .

위와 같은 모델에서 두 행렬 W와 H를 추정하기 위해서, 먼저 행렬 V와 행렬 WH 사이의 거리 함수 $D (V | W H)$ 를 선정한다. Lee와 Seung^[9]은 이와 같은 거리 함수로 먼저 유클리드 거리와 쿨백-라이블러 발산을 활용하였으며, 쿨백-라이블러 발산을 활용하는 경우 거리 함수는 다음과 같이 정의된다.

(2)

D (V ∣ WH) = \sum_{i, j} \{(V)_{i, j} \log \frac{(V)_{i, j}}{(WH)_{i, j}} - (V)_{i, j} + (WH)_{i, j}\},

여기서 $(A)_{i, j}$ 는 행렬 A의 $i, j$ 번째 원소를 의미한다. Lee와 Seung^[9]에 의해 고안된 휴리스틱 곱셈 갱신 기법은 경사하강법의 일종으로, Eq. (2)와 같은 거리 함수의 미분을 이용해 두 행렬 W와 H를 추정한다. 예를 들어, Eq. (2)에 대해 일반적인 경사하강법을 활용하면 다음과 같이 $(H)_{r, n}$ 를 추정할 수 있다.

(3)

(H)_{r, n} \leftarrow (H)_{r, n} - η_{r, n} \nabla_{(H)_{r, n}} D (V ∣ WH),

여기서 $\nabla_{(H)_{r, n}} D (V ∣ WH)$ 는 $(H)_{r, n}$ 에 대한 Eq. (2)의 미분을 의미하며, 이를 계산하면 다음과 같은 경사하강법 기반의 갱신식을 얻을 수 있다.^[9]

(4)

(H)_{r, n} \leftarrow (H)_{r, n} + η_{r, n} \{\sum_{k} (W)_{k, r} \frac{(V)_{k, n}}{(WH)_{k, n}} - \sum_{k} (W)_{k, r}\} .

스텝사이즈를 $η_{r, n} = (H)_{r, n} / \sum_{k} (W)_{k, r}$ 와 같이 설정하면, Eq. (4)의 갱신식은 다음과 같이 나타낼 수 있다.^[9]

(5)

(H)_{r, n} \leftarrow (H)_{r, n} \frac{\sum_{k} (W)_{k, r} \frac{(V)_{k, n}}{(WH)_{k, n}}}{\sum_{k} (W)_{k, n}} .

동일한 방법으로 $(W)_{k, r}$ 의 갱신식도 다음과 같이 얻을 수 있다.^[9]

(6)

(W)_{k, r} \leftarrow (W)_{k, r} \frac{\sum_{n} (H)_{r, n} \frac{(V)_{k, n}}{(WH)_{k, n}}}{\sum_{n} (W)_{r, n}} .

한편, 행렬 V, W, H 의 원소가 모두 비음수이기 때문에, Eq. (4)의 식은 다음과 같이 양의 기울기 $\nabla_{(H)_{r, n}}^{+} D (V ∣ WH)$ 와 음의 기울기 $\nabla_{(H)_{r, n}}^{-} D (V ∣ WH)$ 의 조합으로 일반화할 수 있다.^[9,14]

(7)

(H)_{r, n} \leftarrow (H)_{r, n} + η_{r, n} \{\nabla_{(H)_{r, n}}^{-} D (V ∣ WH) - \nabla_{(H)_{r, n}}^{+} D (V ∣ WH)\} .

Eqs. (4)와 (7)을 비교해 보면, Eq. (4)의 중괄호 안쪽 두 번째 항 $\sum_{k} (W)_{k, r}$ 은 Eq. (7)의 중괄호 안쪽 두 번째 항 $\nabla_{(H)_{r, n}}^{+} D (V ∣ WH)$ 에 해당한다. 따라서, Eq. (4) 로부터 Eq. (5)를 얻을 때 사용한 스텝 사이즈 $η_{r, n} = (H)_{r, n} / \sum_{k} (W)_{k, r}$ 는 Eq. (7)에 대해 다음과 같이 일반화 될 수 있다.

(8)

η_{r, n} = \frac{(H)_{r, n}}{\nabla_{(H)_{r, n}}^{+} D (V ∣ WH)} .

이 경우, Eq. (7)의 식은 다음과 같이 변형된다.

(9)

(H)_{r, n} \leftarrow (H)_{r, n} \frac{\nabla_{(H)_{r, n}}^{+} D (V ∣ WH)}{\nabla_{(H)_{r, n}}^{-} D (V ∣ WH)} .

마찬가지로, $(W)_{k, r}$ 에 대한 갱신식도 다음과 같이 얻어진다.

(10)

(W)_{k, r} \leftarrow (W)_{k, r} \frac{\nabla_{(W)_{k, r}}^{+} D (V ∣ WH)}{\nabla_{(W)_{k, r}}^{-} D (V ∣ WH)} .

이와 같이 임의의 비용 함수 $D (V ∣ WH)$ 에 대해 Eqs. (9)와 (10)을 반복하여 행렬 W와 H를 추정하는 방법을 휴리스틱 곱셈 갱신 기반의 비음수 행렬 분해 기법이라고 한다.

2.2 비음수 행렬 분해 기반의 능동 소나 잔향 제거 기법

Fig. 1에서 보는 바와 같이, 비음수 행렬 분해 기반의 능동 소나 잔향 제거 기법은 입력 음향 신호의 크기 스펙트로그램 $V \in ℝ_{+}^{K \times N}$ 를 행렬 $W \in ℝ_{+}^{K \times R}$ 및 $H \in ℝ_{+}^{R \times N}$ 의 곱으로 분해한다. 이 때 $K$ 는 주파수 빈의 개수, $N$ 은 시간 프레임의 개수가 된다. 이 경우 행렬 W와 H는 각각 음향 이벤트들의 주파수 특성과 시간 특성을 나타내는 주파수 기저 행렬과 시간 기저 행렬이 되며, $R$ 은 비음수 행렬 분해 기법의 파라미터 중 하나로 분석하고자 하는 음향 이벤트의 개수가 된다.

https://cdn.apub.kr/journalsite/sites/ask/2025-044-02/N0660440202/images/ASK_44_02_02_F1.jpg

Fig. 1.

An illustrative diagram for NMF-based active sonar reverberation method.^[15]

능동 소나 잔향 제거 기법에서는 음향 이벤트를 표적 반향과 주변 산란체로부터의 잔향으로 나누어 분석을 하고자 하며, 이를 위하여 다음과 같이 각 기저 행렬을 두 그룹으로 나눈다.^[12]

(11)

W = [W_{P} ∣ W_{R}],

(12)

H = {[H_{P} {}^{T}∣ {H_{R}}^{T}]}^{T},

여기서 W_P와 H_P는 각각 표적 반향 신호의 주파수 및 시간 기저 행렬, W_R와 H_R은 각각 산란체 잔향 신호의 주파수 및 시간 기저 행렬을 의미한다. 이 경우, Eq. (1) 의 비음수 행렬 분해 모델은 다음과 같은 의미를 가진다.

(13)

V \approx WH = \underset{V_{P}}{\underset{⏟}{W_{P} H_{P}}} + \underset{V_{R}}{\underset{⏟}{W_{R} H_{R}}} .

즉, 입력 신호의 크기 스펙트로그램을 표적 반향 성분과 산란체 잔향 성분으로 나눌 수 있으며, 적절한 제약조건의 활용을 통해 W_P와 H_P를 추정하면 그 곱을 통해 표적 반향 신호의 크기 스펙트로그램을 구할 수 있다.

Reference [12]는 이와 같이 W_P와 H_P를 추정하기 위하여 다음과 같이 가정한다. 1) 지속파 송신 신호의 반향에 대한 주파수 특성을 나타내는 W_P는 송신 신호가 도플러 이동된 레플리카의 주파수 특성을 가진다. 2) 잔향의 크기 스펙트로그램인 V_R은 시간축에서 무작위로 요동치는 특성을 가질 수 있는 반면, 표적 반향의 크기 스펙트로그램인 V_P는 시간축에서 상대적으로 더 연속적인 경향을 보인다. 3) 표적 반향의 크기 스펙트로그램인 V_P는 시간축에서 길이가 한정적이며, 이는 송신 핑 길이와 유사한 정도의 길이를 가진다.

위의 가정을 활용하여, Reference [12]의 알고리즘은 다음과 같이 도출되었다. 먼저, 2)의 가정에 의해 H_P는 시간 연속성을 가지게 되며, 3)의 가정에 의해 H_P는 시간 축에서 길이가 한정적인 특징을 가지게 된다. 따라서, H_P 를 추정하기 위한 비용함수는 다음과 같이 구성한다.^[12]

(14)

C (H_{P}) = C_{E} (V ∣ WH) + α C_{T} (H_{P}) + γ C_{L} (H_{P}),

여기서 𝛼와 𝛾는 각각 시간 연속성 비용함수와 시간 길이 제한 비용함수의 가중치를 나타낸다. 시간 연속성 비용함수는 시간 축의 변화량이 클수록 더 큰 값을 가지도록 다음과 같이 설계되었다.^[12]

(15)

C_{T} (H_{P}) = \sum_{r = 1}^{R} \frac{\sum_{n = 2}^{N} {\{{(H_{P})}_{r, n} - {(H_{P})}_{r, n - 1}\}}^{2}}{\frac{1}{N} \sum_{n = 1}^{N} {\{{(H_{P})}_{r, n}\}}^{2}} .

또한, 시간 길이 제한 비용함수는 다음과 같이 특정 길이 $l_{n}$ 보다 긴 신호에 대해서 비용을 크게 하도록 설계되었다.^[12]

(16)

C_{L} (H_{P}) = \sum_{r = 1}^{R} (1 - \sum_{m = n}^{n + l_{n} - 1} {(\overset{\lor}{H_{P}})}_{r, m}),

여기서 $\overset{\lor}{H_{P}}$ 는 $l_{n}$ 구간 단위 기준으로 최대의 에너지를 가지는 구간에서 1에 가까운 값을 가지고 나머지 구간에서 0에 가까운 값을 가지는 구간 최대값 지시 행렬로, 다음과 같이 softmax 함수를 활용하여 정의된 값이다.

(17)

{(\overset{\lor}{H_{P}})}_{r, n} = \frac{\exp \{{({\bar{H}}_{P})}_{r, n}\}}{\sum_{m = 1}^{N} \exp \{{({\bar{H}}_{P})}_{r, m}\}},

여기서 ${\bar{H}}_{P}$ 는 다음과 같이 $l_{n}$ 길이로 이동-평균을 취한 값이다.

(18)

{({\bar{H}}_{P})}_{r, n} = \sum_{m = n - l_{n} + 1}^{n} {(H_{P})}_{r, m} .

시간 길이 제한 비용함수는 Reference [12]에서 독창적으로 설계된 비용함수로, 구체적인 의미는 Reference [12]의 Appendix 그래프에서 확인할 수 있다.

위의 관계를 활용하여, H_P는 Eq. (14)를 미분한 후 이에 대해 Eq. (9)와 같은 휴리스틱 곱셈 갱신을 적용하면 다음과 같이 갱신식을 얻을 수 있다.^[12]

(19)

{(H_{P})}_{r, n} \leftarrow {(H_{P})}_{r, n} \frac{\{\nabla_{{(H_{P})}_{r, n}}^{-} C_{E} (V ∣ WH) + α \nabla_{{(H_{P})}_{r, n}}^{-} C_{T} (H_{P}) + γ \nabla_{{(H_{P})}_{r, n}}^{-} C_{L} (H_{P})\}}{\{\nabla_{{(H_{P})}_{r, n}}^{+} C_{E} (V ∣ WH) + α \nabla_{{(H_{P})}_{r, n}}^{+} C_{T} (H_{P}) + γ \nabla_{{(H_{P})}_{r, n}}^{+} C_{L} (H_{P})\}} .

1)의 가정에 의해, W_P는 Fig. 1에서 보는 바와 같이 지속파 송신 신호의 도플러 이동된 레플리카의 주파수 특성을 가지는 것으로 설정되며, 추정 과정에서 값이 변하지 않고 고정된다. W_R과 H_R은 2)와 3)의 가정을 따르지 않으므로 추가 제약조건을 적용하지 않고 추정 에러 $C_{E} (V ∣ WH)$ 만을 비용함수로 사용하여 최적화하며, 휴리스틱 곱셈 갱신을 적용하면 다음과 같은 갱신식을 얻는다.^[12]

(20)

{(H_{R})}_{r, n} \leftarrow {(H_{R})}_{r, n} \frac{\nabla_{{(H_{R})}_{r, n}}^{+} C_{E} (V ∣ WH)}{\nabla_{{(H_{R})}_{r, n}}^{-} C_{E} (V ∣ WH)} .

(21)

{(W_{R})}_{k, r} \leftarrow {(W_{R})}_{k, r} \frac{\nabla_{{(W_{R})}_{k, r}}^{+} C_{E} (V ∣ WH)}{\nabla_{{(W_{R})}_{k, r}}^{-} C_{E} (V ∣ WH)} .

위 식들의 각 미분 값들은 Reference [12]에서 확인할 수 있다. Eqs. (19), (20), (21) 의 갱신식을 수렴할 때까지 반복하여 계산하면 H_P, W_R, H_R을 추정할 수 있으며, 추정이 끝난 후 W_P와 H_P의 값을 곱하면 표적 반향의 크기 스펙트로그램 V_P를 추정할 수 있다.

III. 적은 반복 횟수에 대한 성능 개선 방안

3.1 L1-노름을 이용한 추가 제약 조건

위와 같은 비음수 행렬 분해 기반의 잔향 제거 기법은 기존의 기법에 비해 향상된 성능을 보여주지만, 실제 적용하기에는 Eqs. (19), (20), (21)의 추정 과정을 반복하여 연산하여야 한다는 문제가 있다. 따라서, 반복 연산의 횟수를 적게 수행하는 상황에서의 성능을 개선할 필요가 있다.

본 논문에서는 이와 같이 적은 반복 횟수, 즉 반복 추정 초기의 성능을 향상시키기 위해서 두 가지 측면을 주목하였다. 첫 번째는, 수신 신호의 크기 스펙트로그램에서 표적 반향음이 활성화된 시간-주파수 빈의 개수가 적다는 점이며, 이는 표적 반향음이 희박함을 의미한다. 따라서, 다음과 같이 L1-노름을 이용한 희박 제약조건을 추가하여 식을 도출한다.

(22)

C (H_{P}) = C_{E} (V ∣ WH) + α C_{T} (H_{P}) + δ C_{S} (H_{P}) + γ C_{L} (H_{P}),

(23)

C_{S} (H_{P}) = \sum_{k, r} {||h_{P, r}||}_{1} {(W_{P})}_{k, r} = \sum_{r, n} {||w_{P, r}||}_{1} {(H_{P})}_{k, r},

여기서 $h_{P, r}$ 은 행렬 H_P의 $r$ 번째 행벡터를, $w_{P, r}$ 은 행렬 W_P의 $r$ 번째 열벡터를 의미한다. 위의 L1-노름 제약조건은 Reference [14]를 참고하여 설계하였다.

3.2 주요화최소화 기법을 활용한 갱신 기법

적은 반복 횟수에서의 성능을 향상시키기 위해서 본 논문에서 고려한 두 번째 측면은, 기존에 활용하던 휴리스틱 곱셈 갱신 기법이 수렴을 보장하지 않는다는 점이다. Eqs. (9)와 (10)과 같이 도출된 휴리스틱 곱셈 갱신 기법은 경사하강법에 기반하여 도출되었으며, 경사하강법은 스텝사이즈가 충분히 작은 경우 일반적으로 국소점에 잘 수렴한다. 그러나 Eq. (8)과 같이 정의된 스텝사이즈가 충분히 작은 값인가에 대해서는 검증된 바가 없다. 물론, 추가적인 제약조건이 없이 쿨백-라이블러 함수만을 최적화 비용 함수로 활용한 Eqs. (5)와 (6)의 경우에는 반복 갱신을 통해 발산하지 않고 수렴함을 Reference [9]에서 보였지만, 제약조건이 존재하는 경우에 대해서는 수렴을 보장할 수 없다.

Eq. (14)와 같이 구성된 비용함수 중 제약 조건에 대한 가중치 𝛼와 𝛾를 잘 조절하면 발산하지 않도록 조절할 수 있어, 기존의 알고리즘에서 수렴을 수학적으로 보장하지 않는 것이 큰 문제를 일으키지는 않는다. 그럼에도 불구하고, 후술할 주요화최소화 기법은 갱신 과정에서 $i$ +1번째 반복과정의 비용함수가 $i$ 번째 반복과정의 비용함수보다 작거나 같도록 하기 때문에, 이를 보장하지 않는 휴리스틱 갱신 기법에 비해 더 빠르게 수렴하기를 기대할 수 있다. 또한, 이는 우리가 다루고 있는 문제와 같이 적은 반복 횟수의 갱신 과정을 수행하였을 때 원하는 정답에 더 가깝기를 기대할 수 있게 한다.

3.2.1 주요화최소화 기법

주요화최소화 기법의 도출 과정은 크게 두 가지의 과정으로 이루어진다. 첫 번째, 즉 주요화 단계는 비용 함수에 대한 보조 함수를 설정하는 단계이다. 만약 임의의 정의역 $E$ 에서 정의되는 임의의 행렬 X 에 대한 비용함수를 $C (X)$ 라 하고, 현재 반복 단계에서 갱신되기 전 행렬을 $\tilde{X}$ 라 하자. 이 때, $C (X)$ 에 대한 보조 함수 $G (X ∣ \tilde{X})$ 는 다음과 같은 관계를 가지는 함수로 정의된다.^[14]

(24)

\forall X \in E, G (X ∣ \tilde{X}) \geq C (X) G (\tilde{X} ∣ \tilde{X}) = C (\tilde{X}) .

최소화 단계에서는 위의 보조 함수 $G (X ∣ \tilde{X})$ 를 최소화하는 행렬 $\tilde{X}$ 를 구하며, 일반적으로는 보조함수의 미분이 0이 되는 점 $\hat{X}$ 을 찾아서 구하게 된다. 이와 같은 점에서는 $G (\hat{X} ∣ \tilde{X}) \leq G (\tilde{X} ∣ \tilde{X})$ 의 관계를 만족하기 때문에, 다음과 같이 비용 함수를 최소화할 수 있게 된다.

(25)

C (\hat{X}) \leq G (\hat{X} ∣ \tilde{X}) \leq G (\tilde{X} ∣ \tilde{X}) = C (\tilde{X}) .

3.2.2 표적 반향 기저 추정을 위한 보조 함수

Eq. (13)에서 보는 바와 같이, 표적 반향 신호의 스펙트로그램은 표적 반향 기저 W_P와 H_P의 곱으로 이루어진다. 이 중 W_P는 송신 파형의 도플러 레플리카의 주파수 특성으로 사전에 결정 및 고정되어 있는 값이므로, 표적 반향 신호를 구하기 위해서는 H_P의 추정이 필요하다.

본 논문에서 설정한 최적화 비용 함수는 Eq. (22)에서 보는 바와 같이 추정 비용 함수 $C_{E} (V ∣ WH)$ , 시간 연속성 비용 함수 $C_{T} (H_{P})$ , 희박 비용 함수 $C_{S} (H_{P})$ , 그리고 시간 길이 제한 함수 $C_{L} (H_{P})$ 의 가중치 합으로 이루어져 있다. 이 중 $C_{L} (H_{P})$ 은 주요화최소화를 위한 적절한 보조 함수를 찾기 어렵기 때문에, 다음과 같이 이를 제외한 나머지 항에 대한 보조 함수 $G (H_{P} ∣ {\tilde{H}}_{P})$ 를 도출한다.

(26)

G (H_{P} ∣ {\tilde{H}}_{P}) = G_{K L} (H_{P} ∣ {\tilde{H}}_{P}) + α C_{T} (H_{P}) + δ C_{S} (H_{P}) .

이 중 $G_{K L} (H_{P} ∣ {\tilde{H}}_{P})$ 에 대한 보조 함수는 이미 많은 학자들에 의해 연구가 되었으며, 본 논문에서는 Reference [14]를 참고하여 다음과 같이 정의하였다.

(27)

G_{K L} (H_{P} ∣ {\tilde{H}}_{P}) = \sum_{r, n} \{{(W_{P} {}^{T}1_{K \times N})}_{r, n} {(H_{P})}_{r, n} - {(W_{P} {}^{T}{\frac{V}{W \tilde{H}}})}_{r, n} {({\tilde{H}}_{P})}_{r, n} \log \frac{{(H_{P})}_{r, n}}{{({\tilde{H}}_{P})}_{r, n}}\},

여기서 행렬에 대한 분수식 $\frac{A}{B}$ 는 두 행렬의 원소에 대한 나눗셈을 의미한다.

Eq. (26)의 $C_{T} (H_{P})$ 와 $C_{S} (H_{P})$ 의 함수도 필요한 경우 적절한 보조함수를 찾아야 하지만, $C_{S} (H_{P})$ 의 함수는 보조함수를 정의하지 않아도 갱신식의 도출이 가능하다. 시간 연속성 비용 함수 $C_{T} (H_{P})$ 의 경우 기존의 Eq. (15)에서 정의된 형태로는 갱신식의 도출이 어려워서 변형이 필요하며, 본 논문에서는 쿨백-라이블러 발산 함수를 활용하여 다음과 같이 변형된 시간 연속성 비용 함수를 정의하였다.

(28)

C_{T}^{'} (H_{P}) = \sum_{r, n} \{D ({(H_{P})}_{r, n - 1} ∣ {(H_{P})}_{r, n}) + D ({(H_{P})}_{r, n + 1} ∣ {(H_{P})}_{r, n})\} .

Eq. (15)는 인접한 두 값의 유클리드 거리를 이용하여 정의된 식이며 Eq. (28)은 인접한 두 값의 쿨백-라이블러 발산으로 정의된 식으로, 두 식 모두 시간 축에서 인접한 값의 차이가 클수록 함수의 값이 커지는 비용 함수라는 점에서 동일하다.

3.2.3 보조 함수를 최소화하는 갱신식

주요화최소화 기법의 두 번째 단계인 최소화 단계에서는, 다음과 같이 주요화 단계에서 얻어진 보조 함수를 최소화하는 ${\hat{H}}_{P}$ 를 찾는다.

(29)

{\hat{H}}_{P} = \underset{H_{P}}{\arg \min G (H_{P} ∣ {\tilde{H}}_{P})} .

Eq. (29)를 만족하는 값을 찾기 위해서, Eq. (26)의 식을 H_P로 미분한 뒤 이를 0으로 만드는 ${\hat{H}}_{P}$ 를 구하면 다음과 같이 얻어진다(상세한 과정은 부록을 참고할 수 있다).

(30)

{(H_{P})}_{r, n} = \frac{[{(W_{P} {}^{T}{\frac{V}{W \tilde{H}}})}_{r, n} {({\tilde{H}}_{P})}_{r, n} + α \{{(H_{P})}_{r, n - 1} + {(H_{P})}_{r, n + 1}\}]}{[{(W_{P} {}^{T}1_{K \times N})}_{r, n} + 2 α + δ {||w_{P, r}||}_{1}]} .

표적 반향 주파수 기저 행렬 W_P의 경우 사전에 설계된 값으로 고정되어 있기 때문에, 추정 과정에서 그 값이 변하지 않는다. 따라서, 사전 설계 과정에서 행렬 W_P의 열벡터의 합이 1이 되도록 설계한 경우 ${||w_{P_{, r}}||}_{1} = 1$ 이 된다.

Eq. (30)은 주요화최소화 기법으로 Eq. (26)의 비용함수를 최적화할 수 있는 갱신식이다. 다만, 앞서 언급한 바와 같이 해당 비용함수에는 시간 길이 제약 조건 $C_{L} (H_{P})$ 가 누락되어 있으며, 알고리즘이 제대로 동작하기 위해서는 이를 고려한 갱신식이 필요하다. 불행히도 이는 주요화최소화 기법으로 얻을 수 없으므로, 본 논문에서는 $C_{L} (H_{P})$ 에 한하여 휴리스틱 곱셈 갱신법을 활용하여 갱신식을 적용하였다. Eq. (9)와 Reference [12]에서 확인할 수 있는 바와 같이, 휴리스틱 곱셈 갱신법을 활용하면 Eq. (30)의 갱신식의 분자에는 ${(\nabla_{{(H_{P})}_{r, n}}^{+} C_{L} (H_{P}))}_{r, n} {({\tilde{H}}_{P})}_{r, n}$ 를 더해주고, 분모에는 ${(\nabla_{{(H_{P})}_{r, n}} C_{L} (H_{P}))}_{r, n}$ 를 더해주면 된다. 이를 정리하면 $C_{L} (H_{P})$ 까지 고려한 갱신식을 다음과 같이 쓸 수 있다.

(31)

P_{1} = {(W_{P}^{T} \frac{V}{WH})}_{r, n} {(H_{P})}_{r, n},

(32)

P_{2} = {(H_{P})}_{r, n - 1} + {(H_{P})}_{r, n + 1},

(33)

P_{3} = \sum_{m = n}^{n + l_{n} - 1} \frac{\exp ({({\bar{H}}_{P})}_{r, m})}{\sum_{i = 1}^{N} \exp ({({\bar{H}}_{P})}_{r, i})} {(H_{P})}_{r, n},

(34)

M_{1} = {(W_{P} {}^{T}1_{K \times N})}_{r, n},

(35)

M_{2} = \sum_{m = n}^{n + l_{n} - 1} {\{\frac{\exp ({({\bar{H}}_{P})}_{r, m})}{\sum_{i = 1}^{N} \exp ({({\bar{H}}_{P})}_{r, i})}\}}^{2},

(36)

{(H_{P})}_{r, n} \leftarrow \frac{P_{1} + α P_{2} + γ P_{3}}{M_{1} + γ M_{2} + 2 α + δ},

여기서 ${\bar{H}}_{P}$ 는 다음과 같이 $H_{P}$ 행렬의 값을 $l_{n}$ 길이로 시간 축에서 이동합산을 취한 값이다.^[12]

(37)

{({\bar{H}}_{P})}_{r, n} = \sum_{m = n - l_{n} + 1}^{n} {(H_{P})}_{r, m} .

Eqs. (33)과 (35)의 도출과정은 Reference [12]에서 살펴볼 수 있으며, 본 논문에서 도출된 알고리즘 또한 Reference [12]와 동일한 비용 함수에 대해 동일한 알고리즘, 즉 휴리스틱 곱셈 갱신법을 이용하여 도출하였기 때문에 그 과정 또한 동일하다. 따라서 편의를 위하여 본 논문에서는 도출과정을 생략하여 기술하였다.

3.2.4 잔향 기저에 대한 갱신식

Reference [12]에서 살펴볼 수 있는 바와 같이, 잔향 기저를 갱신할 때에는 추가적인 제약조건 없이 V와 WH 사이의 쿨백-라이블러 발산함수 $D (V ∣ WH)$ 를 최소화한다. 이는 Eq. (30)의 도출 결과에서 𝛼=0 및 𝛿=0인 경우와 동일하며, 다음과 같이 기술할 수 있다.

(38)

(H_{R}) \leftarrow \frac{{(W_{R}^{T} \frac{V}{WH})}_{r, n} {(H_{R})}_{r, n}}{{(W_{P}^{T} 1_{K \times N})}_{r, n}},

(39)

(W_{R}) \leftarrow \frac{{(\frac{V}{WH} H_{R}^{T})}_{r, n} {(W_{R})}_{r, n}}{{(1_{K \times N} H_{R}^{T})}_{r, n}} .

Eqs. (38)과 (39)는 휴리스틱 곱셈 갱신 기법으로 도출된 수식과 수학적으로 동일하며, 이와 같이 추가적인 제약 조건이 존재하지 않는 상황에서 주요화최소화 기법과 휴리스틱 곱셈 갱신 기법이 동일한 결과를 보이는 것은 Reference [14] 등에서 연구된 바와 일맥상통하는 결과이다.

본 논문에서 제안된 알고리즘이 Table 1에 요약되어 있다. Table 1에서 보는 바와 같이, 본 알고리즘은 수신된 소나 신호의 크기 스펙트로그램 V에 대해서 반복적인 갱신 과정을 통하여 W_R, H_R, H_P를 추정하며, 갱신 과정이 수렴하고 난 후 사전에 정의된 W_P와 추정된 H_P와의 곱을 통하여 잔향이 제거된 표적 반향 신호의 크기 스펙트로그램을 추정한다.

Table 1.

Summary of the proposed algorithm.

Majorization-minimization-based NMF algorithm for active sonar reverberation suppression

Initialization: W_P is initialized with frequency structures of transmitted ping and its Doppler-shifted replicas. W_R, H_P, and H_R are initialized by absolute values of Gaussian random numbers.

Input: Magnitude spectrogram

V \in ℝ_{+}^{K \times N}

of received signal.

Iterations:
1) W_R is updated using Eq. (39);
2) W is updated by Eq. (11);
3) H_P is updated using Eqs. (31), (32), (33), (34), (35), (36);
4) H_R is updated using Eq. (38);
5) H is updated by Eq. (12).

After convergence: Magnitude spectrogram V_P of the output signal is reconstructed by Eq. (13).

IV. 시뮬레이션

4.1 시뮬레이션 환경

본 논문에서 고안된 비음수 행렬 분해 기반의 잔향 제거 기법의 성능을 확인하기 위하여 PC 환경에서 MATLAB을 활용한 시뮬레이션을 진행하였다. 본 시뮬레이션에서 가정한 시뮬레이션 환경은 Reference [12] 및 Reference [15]와 유사하도록 설정하였으며, 소나 송수신기는 100 ms 길이의 지속파 송신파형을 송신하면서 $2 V / c = 0.033$ 이 되는 속력 $V$ 로 이동한다. 잔향 신호는 Reference [16]의 비-레일레이 잔향 환경 모델을 기반으로 합성되었으며, 표적 반향 신호의 도플러 주파수는 $f_{d} = 0.005 f_{0}$ 로 가정하였으며, 이는 송수신기의 속력에 의한 잔향 퍼짐에 의해 표적 반향의 탐지가 방해를 받는 상황이 되도록 가정한 것이다. 잔향을 포함한 수신 신호의 길이는 5 s로 설정되었다.

생성된 입력 신호는 75 % 중첩된 16 ms의 해밍 윈도우가 적용된 신호에 대해 128개의 주파수 빈을 가지는 단시간 푸리에 변환을 활용하여 크기 스펙트로그램으로 변환된다. 제안하는 알고리즘의 파라미터인 표적 반향 기저의 개수는 17개, 잔향 기저의 개수는 60개로, 전체 기저의 개수 $R$ 의 값이 77이 되도록 설정하였고, 시간 연속성 제약 조건의 가중치 𝛼는 0.001, 희박 제약 조건의 가중치 𝛿는 1, 시간 길이 제약 조건의 가중치 𝛾는 0.1로 설정하였다. 시간 길이 제한 요소 $l_{n}$ 은 100 ms에 해당하는 샘플 수로 설정하였다. 특히, 비음수 행렬 분해 과정의 반복 추정 횟수를 10회로 설정하여, 적은 연산량 환경에서 성능을 확보할 수 있는지 확인하고자 하였다.

알고리즘의 성능 비교 대상으로는 기존의 비음수 행렬 분해 기반 잔향 제거 기법^[12]과 주성분 역산^[6] 기법을 설정하였으며, 비음수 행렬 분해 기반 기법의 파라미터인 $R$ , 𝛼, 𝛾, 그리고 반복 추정 횟수는 제안하는 알고리즘과 동일하게 설정하였다. 주성분 역산 기법의 고유값 문턱값은 잔향 에너지를 알고 있다는 가정을 통해 이상적인 값으로 설정하였다.

신호 잔향 환경은 최대 –8 dB부터 최소 –23 dB의 신호대잔향비 환경에 대해 3 dB 간격으로 성능을 평가하였으며, 각 신호 환경에서 100회의 몬테-카를로 시뮬레이션을 수행한 후 성능지표의 평균값을 획득하는 방법으로 진행하였다. 출력 신호의 신호대잡음비를 성능 지표로 활용하였으며, 이는 다음과 같이 정의된다.

(40)

S N R = \frac{\sum_{n} {|s_{c} (n)|}^{2}}{\sum_{n} {|s_{c} (n) - {\hat{s}}_{c} (n)|}^{2}},

여기서 $s_{e} (n)$ 는 잔향이 없는 이상적인 핑 신호를, ${\hat{s}}_{e} (n)$ 는 알고리즘의 출력 신호를 나타낸다.

4.2 결과

Fig. 2는 출력 신호대잡음비(Signal-to-Noise Ratio, SNR)에 대한 성능 평가 결과를 도시하고 있다. Fig. 2의 가로축은 입력 신호대잔향비(Signal-to-Reverberation Ratio, SRR) 환경을, 세로축은 성능지표인 출력 신호대잡음비를 나타낸다. 적색(원형마크) 그래프는 제안하는 알고리즘의 성능을, 청색(삼각형마크) 및 흑색(X마크) 그래프는 각각 비교 대상인 기존 비음수 행렬 분해 알고리즘^[12] 및 주성분 역산^[6] 알고리즘의 결과를 나타낸다. Fig. 2의 결과를 통해 제안하는 알고리즘이 기존의 알고리즘에 비해 우수한 성능을 보이는 것을 확인할 수 있으며, 특히, 낮은 입력 신호대잔향비 환경에서 성능이 개선되는 것을 확인할 수 있다.

https://cdn.apub.kr/journalsite/sites/ask/2025-044-02/N0660440202/images/ASK_44_02_02_F2.jpg

Fig. 2.

(Color available online) Performance comparisons in various input SRR conditions.

Figs. 3, 4, 5는 –17 dB의 입력 신호대잔향비 환경을 기준으로 파라미터 변화에 따른 제안하는 알고리즘의 성능을 도시하고 있다. Fig. 3은 시간 연속성 제약조건의 가중치인 𝛼의 변화에 따른 성능 변화를 보여주고 있으며, 𝛼의 값을 10^-5에서 1까지 10배 간격으로 값을 변화시키면서 그 결과를 확인하였다. 결과 그래프에 따르면 𝛼의 값이 10^-4와 10^-2 사이일 때 안정적인 성능을 보이는 것을 확인할 수 있으며, 이보다 큰 값에서는 성능이 저하되는 것을 확인할 수 있다.

https://cdn.apub.kr/journalsite/sites/ask/2025-044-02/N0660440202/images/ASK_44_02_02_F3.jpg

Fig. 3.

(Color available online) SNR performance comparisons in various parameter values for temporal continuity constraints with input SRR condition of –17 dB,

https://cdn.apub.kr/journalsite/sites/ask/2025-044-02/N0660440202/images/ASK_44_02_02_F4.jpg

Fig. 4.

(Color available online) SNR performance comparisons in various parameter values for temporal length limitation constraint with input SRR condition of –17 dB.

https://cdn.apub.kr/journalsite/sites/ask/2025-044-02/N0660440202/images/ASK_44_02_02_F5.jpg

Fig. 5.

(Color available online) SNR performance comparisons for various parameter values for sparsity constraint with input SRR condition of –17 dB.

Fig. 4는 시간 길이 제한 제약조건의 가중치인 𝛾 에 대한 성능 변화를 보여주고 있으며, 𝛾의 값을 10^-5에서 10²까지 10배 간격으로 값을 변화시키면서 그 결과를 확인하였다. 결과 그래프에 따르면 𝛾의 값이 10^-3과 1 사이일 때 안정적인 성능을 보여주는 것을 확인할 수 있으며, 이보다 작거나 큰 값에서 성능이 저하되는 것을 확인할 수 있다.

Fig. 5는 희박 제약조건의 가중치인 𝛿에 따른 성능 차이를 보여준다. 희박 제약 조건의 경우 𝛿 가 0.1보다 큰 경우에서 성능이 향상되는 것을 확인할 수 있으며, 𝛿\ 값이 10인 경우의 성능이 가장 좋은 결과를 보인다. 다만, Fig. 6(a)에서 보는 바와 같이 𝛿=1인 경우 표적 반향에 해당하는 크기 스펙트로그램의 시간-주파수 빈 에너지가 약 12.6 dB이었으며, 이는 𝛿=10 인 경우를 제외한 대부분의 경우 유사한 경향을 보였다. 반면 𝛿=10의 경우 Fig. 6(b)에서 보는 바와 같이 표적 반향의 에너지가 –1.1 dB로 크게 감소된 결과를 보였으며, 이는 표적 반향의 에너지가 크게 제거되는 것으로 해석할 수 있다. 따라서, 알고리즘을 안정적으로 활용하기 위해서는 𝛿=1의 값을 활용하는 것이 권장된다.

https://cdn.apub.kr/journalsite/sites/ask/2025-044-02/N0660440202/images/ASK_44_02_02_F6.jpg

Fig. 6.

(Color available online) Magnitude spectrograms of output signals for (a) 𝛿=1 and (b) 𝛿=10 cases with input SRR condition of –17 dB.

Fig. 7은 –17 dB의 입력 SRR 조건에서 비음수 행렬 분해의 반복 추정 횟수에 따른 성능을 보여주고 있다. Fig. 7의 그래프에 의하면 제안하는 알고리즘이 기존의 비음수 행렬 분해 기반의 알고리즘에 비해 적은 횟수의 반복을 활용할 때 안정적인 성능을 보여주고 있으며, 특히 20 이하의 적은 반복 횟수에서 성능 이득이 큰 것을 확인할 수 있다. 반면, 일정 반복 횟수(Fig. 7 그래프에서는 약 80회) 이상에서는 제안하는 알고리즘이 기존 알고리즘에 비해 성능이 감소하는 것을 확인할 수 있는데, 이는 학습 초기의 성능을 향상시키기 위한 희박 제약 조건이 오히려 성능을 저하시키는 방향으로 작용한 것으로 판단된다. 즉, 반복 추정 초기의 경우 넓은 시간-주파수 영역에 걸쳐 존재하는 잔향 신호가 희박 제약 조건에 의해 제거되는 이득이 크게 나타나지만, 반복 추정 횟수가 커지게 되면 희박 제약 조건이 없어도 잔향을 효과적으로 제거하게 되기 때문에, 희박 제약 조건에 의한 부작용이 더 크게 나타나는 것으로 해석할 수 있다. 희박 제약 조건에 의한 부작용은, Fig. 6에서 살펴보는 바와 같이 실제로 여러 시간-주파수 빈에 걸쳐 있는 표적 반향 신호를 과도하게 희박한 분포로 간주하여 일부를 제거하게 되는 것으로 나타난다. 따라서, 반복 추정 횟수를 증가시키는 경우 희박 제약 조건의 가중치를 순차적으로 감소시키는 등의 전략이 유효할 것으로 판단된다.

https://cdn.apub.kr/journalsite/sites/ask/2025-044-02/N0660440202/images/ASK_44_02_02_F7.jpg

Fig. 7.

(Color available online) SNR performance comparisons according to the number of iterations for –17 dB input SRR condition.

Fig. 8은 기존의 알고리즘과 제안하는 알고리즘의 비음수 행렬 분해 연산에 걸리는 시간을 측정하여 도시한 결과이다. 제안하는 알고리즘의 비음수 행렬 분해 기법이 기존의 기법에 비해 추가적인 연산을 필요로 하기 때문에 각 반복 횟수 당 연산 시간이 기존의 기법에 비해 조금 더 긴 것을 확인할 수 있으나, 반복 횟수에 따른 연산 시간의 증가 추세에 비하면 그 차이가 확연히 작은 것을 확인할 수 있다. 따라서, 제안하는 기법을 통하여 반복 연산의 횟수를 줄이는 것이 연산량의 측면에서 유의미하다는 결론을 얻을 수 있다.

https://cdn.apub.kr/journalsite/sites/ask/2025-044-02/N0660440202/images/ASK_44_02_02_F8.jpg

Fig. 8.

(Color available online) Calculation time comparisons for NMF iterations between the conventional and proposed methods.

V. 결 론

본 논문에서는 비음수 행렬 분해 기법 기반의 능동 소나 잔향 제거 기법을 개선하는 기법에 대하여 연구를 진행하였다. 기존의 비음수 행렬 분해 기반 잔향 제거 기법은 반복 연산을 통해 잔향 및 표적 반향 기저를 추정하기 때문에, 추정 성능을 확보하기 위해서는 반복 연산을 충분히 수행해야 한다는 단점이 존재한다. 본 연구에서는 반복 연산 초기에 원하는 해로 빠르게 수렴하여 적은 반복 연산 횟수로도 충분한 성능을 얻을 수 있는 알고리즘을 개발하였다. 본 논문에서 제안하는 알고리즘은 표적 반향이 희박성을 가진다는 점에 착안, L1-노름을 활용하여 수렴 초기에 원하는 해로 빠르게 수렴하도록 하였다. 또한, 기존의 알고리즘이 활용하고 있는 휴리스틱 곱셈 갱신 기법은 추가적인 비용 함수를 활용하는 경우 수렴을 보장하지 못한다는 점에 착안하여, 수학적으로 매 추정 단계마다 수렴을 보장할 수 있는 주요화최소화 기법을 활용한 갱신식을 도출하였다.

본 논문에서 제안한 잔향 제거 기법을 검증하기 위하여 시뮬레이션을 통한 반복실험을 진행하였다. 본 시뮬레이션에서는 동일한 잔향 환경에 대해 100회의 몬테카를로 시뮬레이션을 진행하여 출력 SNR 성능의 평균을 산출하였으며, 특히 반복 추정 횟수를 10회로 제한하여 적은 반복 횟수에서의 잔향 제거 성능을 관찰하였다. 그 결과 제안하는 알고리즘이 기존 알고리즘 대비 약 3 dB의 성능 향상 결과를 얻었으며, 특히 낮은 신호대잔향비 환경에서 성능 향상의 정도가 더 큰 것을 확인할 수 있었다. 또한, MATLAB 기준으로 0.2 s가량의 낮은 연산 시간을 필요로 하는 것 또한 확인할 수 있었다.

제안하는 알고리즘이 기존 대비 향상된 성능을 얻었지만 그 향상 정도에 한계가 있다는 단점이 있다. 특히, 도출 과정에서 제약조건 중 하나인 시간 길이 제한 함수에 대한 보조함수를 도출하지 못하여, 이에 대하여 기존과 같은 휴리스틱 곱셈 갱신 기법을 적용할 수밖에 없었던 아쉬운 점이 있다. 향후 연구에서는 이와 같은 부분을 개선하여 더욱 향상된 성능을 확보하고자 한다. 또한, 실시간 시스템의 적용 가능 여부 또한 추가적인 연구가 필요하다.

Acknowledgements

본 연구는 2024년도 ㈜한화시스템의 재원을 지원받아 수행된 연구이다(계약번호 : U-24-020).

부 록

A.1 표적 반향 기저 갱신식의 도출 상세

Eqs. (22) 및 (28)로부터 다음과 같은 비용 함수를 얻을 수 있다.

(A1)

C (H_{P}) = C_{E} (V ∣ WH) + α C_{T}^{'} (H_{P}) + δ C_{S} (H_{P}) + γ C_{L} (H_{P}) .

Eq. (A1)에 대해 Eq. (9)와 같은 휴리스틱 곱셈 갱신 기법을 적용하면 다음과 같이 얻어진다.

(A2)

{(H_{P})}_{r, n} \leftarrow \frac{\{\nabla_{{(H_{P})}_{r, n}}^{-} C_{E} (V ∣ WH) {(H_{P})}_{r, n} + α \nabla_{{(H_{P})}_{r, n}}^{-} C_{T}^{'} (H_{P}) {(H_{P})}_{r, n} + δ \nabla_{{(H_{P})}_{r, n}}^{-} C_{S} (H_{P}) {(H_{P})}_{r, n} + γ \nabla_{{(H_{P})}_{r, n}}^{-} C_{L} (H_{P}) {(H_{P})}_{r, n}\}}{\{\nabla_{{(H_{P})}_{r, n}} C_{E} (V ∣ WH) + α \nabla_{{(H_{P})}_{r, n}} C_{T} (H_{P}) + δ \nabla_{{(H_{P})}_{r, n}}^{+} C_{S} (H_{P}) + γ \nabla_{{(H_{P})}_{r, n}}^{+} C_{L} (H_{P})\}},

본문에서 언급한 바와 같이, 본 논문의 목적은 Eq. (A2)와 같은 갱신식을 휴리스릭 곱셈 갱신 기법 대신 주요화최소화 기법으로 얻고자 하는 것이다. 불행히도 현재 $C_{L} (H_{P})$ 에 대한 주요화최소화 갱신식을 얻기가 어렵기 때문에, 본 논문에서는 이를 제외한 나머지 비용 함수에 대한 주요화최소화 갱신식을 얻고자 하였다. 이를 위하여 Eq. (A2)를 다음과 같이 간소화하자.

(A3)

{(H_{P})}_{r, n} \leftarrow \frac{\{U^{-} (H_{P}) + γ \nabla_{{(H_{P})}_{r, n}}^{-} C_{L} (H_{P}) {(H_{P})}_{r, n}\}}{U^{+} (H_{P}) + γ \nabla_{{(H_{P})}_{r, n}}^{+} C_{L} (H_{P})},

여기서 $U^{-} (H_{P})$ 와 $U^{+} (H_{P})$ 는 각각 𝛾=0일 때의 갱신식에 대한 분자 및 분모 항이다. 𝛾=0인 경우의 비용 함수는 다음과 같다.

(A4)

C_{E T S} (H_{P}) = C_{E} (V ∣ W H) + α C_{T}^{'} (H_{P}) + δ C_{S} (H_{P}) .

이 비용함수에 대해 주요화최소화 조건 Eq. (25)를 만족하는 보조함수 $G (H_{P} ∣ {\tilde{H}}_{P})$ 는 Eq, (26)과 같다. 이를 ${(H_{P})}_{r, n}$ 에 대해 미분하면,

(A5)

\nabla_{{(H_{P})}_{r, n}} G (H_{P} ∣ {\tilde{H}}_{P}) = \nabla_{{(H_{P})}_{r, n}} G_{K L} (H_{P} ∣ {\tilde{H}}_{P}) + α \nabla_{{(H_{P})}_{r, n}} C_{T}^{'} (H_{P}) + δ \nabla_{{(H_{P})}_{, n}, n} C_{S} (H_{P})

이고, 각 항은 다음과 같다.

(A6)

\nabla_{{(H_{P})}_{r, n}} G_{K L} (H_{P} ∣ {\tilde{H}}_{P}) = {(W_{P}^{T} 1_{K \times N})}_{r, n} - {(W_{P} {}^{T}{\frac{V}{WH}})}_{r, n} \frac{{({\tilde{H}}_{P})}_{r, n}}{{(H_{P})}_{r, n}} .

(A7)

\nabla_{{(H_{P})}_{r, n}} C_{T}^{'} (H_{P}) = 2 - \frac{{(H_{P})}_{r, n - 1} {(H_{P})}_{r, n + 1}}{{(H_{P})}_{r, n}} .

(A8)

\nabla_{{(H_{P})}_{r, n}} C_{S} (H_{P}) = {||w_{P, r}||}_{1} .

본문과 같이 사전 설계 단계에서 ${||w_{P, r}||}_{1} = 1$ 과 같이 설계하였다고 가정하고 위 식을 정리하면,

(A9)

\nabla_{{(H_{P})}_{r, n}} G (H_{P} ∣ {\tilde{H}}_{P}) = {(W_{P} {}^{T}1_{K \times N})}_{r, n} + 2 α + δ - \frac{1}{{(H_{P})}_{r, n}} [{(W_{P} {}^{T}{\frac{V}{WH}})}_{r, n} {({\tilde{H}}_{P})}_{r, n} + α \{{(H_{P})}_{r, n - 1} + {(H_{P})}_{r, n + 1}\}]

과 같이 정리할 수 있다. 보조 함수 $G (H_{P} ∣ {\tilde{H}}_{P})$ 를 최소화하는 최적해 ${({\hat{H}}_{P})}_{r, n}$ 는 보조 함수의 미분 Eq. (A9)를 0으로 만드므로,

(A10)

{(W_{P} {}^{T}1_{I \times N})}_{r, n} + 2 α + δ = \frac{1}{{({\hat{H}}_{P})}_{r, n}} [{(W_{P} {}^{T}{\frac{V}{WH}})}_{r, n} {({\tilde{H}}_{P})}_{r, n} + α \{{(H_{P})}_{r, n - 1} + {(H_{P})}_{r, n + 1}\}]

이고, 이를 ${({\hat{H}}_{P})}_{r, n}$ 에 대해 정리하면 다음과 같은 식을 얻는다.

(A11)

{({\hat{H}}_{P})}_{r, n} = \frac{[{(W_{P} {}^{T}{\frac{V}{WH}})}_{r, n} {({\tilde{H}}_{P})}_{r, n} + α \{{(H_{P})}_{r, n - 1} + {(H_{P})}_{r, n + 1}\}]}{{(W_{P} {}^{T}1_{K \times N})}_{r, n} + 2 α + δ} .

Eqs. (31), (32), (33), (34), (35)의 정의를 이용하여 Eq. (A11)을 간소화 하고, 이를 갱신식의 형태로 다시 표현하면

(A12)

{(H_{P})}_{r, n} \leftarrow \frac{P_{1} + α P_{2}}{M_{1} + 2 α + δ}

와 같고, 이는 𝛾=0일 때의 갱신식에 해당하므로 위 식의 좌변의 분자항 및 분모항은 각각 Eq. (A3)의 $U^{-} (H_{P})$ 와 $U^{+} (H_{P})$ 에 해당한다. 따라서 Eq. (A12)의 분자항 및 분모항을 각각 Eq. (A3)에 대입하면 Eq. (36)의 식을 얻을 수 있다.

References

L. E. Kinsler, A. R. Frey, A. B. Coppens, and J. V. Sanders, Fundamentals of Acoustics, 4th ed (John Wiley & Sons, Hoboken, 2000), pp. 456-464.

T. Collins and P. Atkins, "Doppler-sensitive active sonar pulse designs for reverberation processing," IEE Proc. Radar Sonar Navig. 145, 347-353 (1998).

10.1049/ip-rsn:19982434

S. Ward, "The use of sinusoidal frequency modulated pulses for low Doppler detection," Proc. IEEE Oceans, 2147-2151 (2001).

10.1109/OCEANS.2001.968331

G. Kim, K. Lee, K. Yoon, and S. Lee, "A study on pulse train waveforms for high duty cycle sonar systems: optimization scheme and relationship between orthogonality and bandwidth," IEEE Access, 9, 119800-119817 (2021).

10.1109/ACCESS.2021.3107907

S. Kay and J. Salisbury, "Improved active sonar detection using autoregressive preshiteners," J. Acoust. Soc. Am. 87, 1603-1611 (1990).

10.1121/1.399408

G. Ginolhac and G. Jourdain, "Principal component inverse algorihtm for detection in the presence of reverberation," IEEE J. Oean. Eng. 27, 310-321 (2002).

10.1109/JOE.2002.1002486

W. Li, Q. Zhang, X. Ma, and C. Hou, "Active sonar detection in reverberation via signal subspace extraction algorithm," EURASIP J. Wirel. Commun. Netw. 2010, 1-10 (2010).

10.1155/2010/981045

D. D. Lee and H. S. Seung, "Learning the parts of objects by non-negative matrix factorization," Nature, 401, 788-791 (1999).

10.1038/4456510548103

D. D. Lee and H. S. Seung, "Algorithms for non-negative matrix factorization," Proc. NeurIPS, 556-562 (2001).

P. Smaragdis and J. C. Brown, "Non-negative matrix factorization for polyphonic music transcription," Proc. IEEE WASPAA, 177-180 (2003).

10.1109/ASPAA.2003.1285860

M. Sun, Y. Li, J. F. Gemmeke, and X. Zhang, "Speech enhancement under low SNR conditions via noise estimation using sparse and low-rank NMF with Kullback-Leibler divergence," IEEE/ACM Trans. Audio. Speech. Lang. Process. 23, 1233-1242 (2015).

10.1109/TASLP.2015.2427520

S. Lee and J.-S. Lim, "Reverberation suppression using non-negative matrix factorization to detect low-doppler target with continuous wave active sonar," EURASIP J. Adv. Signal Process. 2019, 1-18 (2019).

10.1186/s13634-019-0608-6

G. Kim, K. Lee, and S. Lee, "Linear frequency modulated reverberation suppression using non-negative matrix factorization methods, dechirping trasnformation, and modulo operation," IEEE Access, 8, 110720-110737 (2020).

10.1109/ACCESS.2020.3001865

A. Marmin, J. H. de Morais Goulart, and C. Fevotte, "Majorization-minimization for sparse nonnegative matrix factorization with the beta-divergence," IEEE Trans. Signal. Process. 71, 1435-1447 (2023).

10.1109/TSP.2023.3266939

S. Lee and G. Kim, "A Study on the active sonar reverberation suppression method based on non-negative matrix factorization with beta-divergence function" (in Korean), J. Acoust. Soc. Kr. 43, 369-382 (2024).

D. A. Abaham and A. P. Lyons, "Simulation of non-rayleigh reverberation and clutter," IEEE J. Oceanic Eng. 29, 347-362 (2004).

10.1109/JOE.2004.828202

The Journal of the Acoustical Society of KoreaISSN:1225-4428(Print) 2287-3775(Online)한국음향학회

Preview

Improvement of non-negative matrix factorization-based active sonar reverberation suppression method using L1-norm and majorization-minimization

ABSTRACT

MAIN

(1)

(2)

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

Fig. 1.

An illustrative diagram for NMF-based active sonar reverberation method.[15]

(11)

(12)

(13)

(14)

(15)

(16)

(17)

(18)

(19)

(20)

(21)

(22)

(23)

(24)

(25)

(26)

(27)

(28)

(29)

(30)

(31)

(32)

(33)

(34)

(35)

(36)

(37)

(38)

(39)

Table 1.

Summary of the proposed algorithm.

(40)

Fig. 2.

(Color available online) Performance comparisons in various input SRR conditions.

Fig. 3.

(Color available online) SNR performance comparisons in various parameter values for temporal continuity constraints with input SRR condition of –17 dB,

Fig. 4.

(Color available online) SNR performance comparisons in various parameter values for temporal length limitation constraint with input SRR condition of –17 dB.

Fig. 5.

(Color available online) SNR performance comparisons for various parameter values for sparsity constraint with input SRR condition of –17 dB.

Fig. 6.

(Color available online) Magnitude spectrograms of output signals for (a) 𝛿=1 and (b) 𝛿=10 cases with input SRR condition of –17 dB.

Fig. 7.

(Color available online) SNR performance comparisons according to the number of iterations for –17 dB input SRR condition.

Fig. 8.

(Color available online) Calculation time comparisons for NMF iterations between the conventional and proposed methods.

Acknowledgements

부 록

(A1)

(A2)

(A3)

(A4)

(A5)

(A6)

(A7)

(A8)

(A9)

(A10)

(A11)

(A12)

References

An illustrative diagram for NMF-based active sonar reverberation method.^[15]