음원의 위치 및 강도 추정을 위한 압축센싱 기법은 등가음원법(equivalent source method)을 이용한 정식화에서 출발한다. 즉, Fig. 1에 나타낸 바와 같이, 어떤 관심 공간 내에 존재하는 음원을 N개의 격자로 이산화, 혹은 N개의 잠재적인 단극음원으로 분포시켜 이들의 조합으로 표현할 수 있다고 가정하자. 이때 M개의 배열 센서(마이크로폰 혹은 하이드로폰)를 이용하여 음장을 계측할 경우, m번째 센서에서 계측된 압력값 y_m은 그린함수를 통해 주파수 영역에서 식 (1)과 같이 나타낼 수 있다.

Fig. 1 
Distribution of N equivalent monopole sources in the interested domain (n-th equivalent source with strength x_n is located at each grid point, m-th sensor in the array receives the signal y_m)

여기서 x_n [Pa⋅m]은 n번째 잠재적 단극음원의 강도이고, r_mn [m]은 n번째 잠재적 단극음원과 m번째 센서의 거리이며, k (=2πf/c) [1/m]는 주파수 f [Hz]에 해당하는 파수를 나타낸다(단, c(=1500) [m/s]는 음속을 의미함). 식 (1)을 M개 센서에 대하여 확장하면 관측벡터(measurement vector) y∈ℂ^M은 식 (2)와 같이 나타낼 수 있다.

여기서 A∈ℂ^M×N은 관측행렬(measurement matrix)이라고 한다.

대부분의 경우, 관심 공간 내 잠재적 단극음원의 개수 N은 센서 개수 M보다 크다(M < N). 따라서 식 (2)의 관측행렬 A는 미지수의 개수(= N)가 식의 개수(= M)보다 많은 부족결정 시스템(underdetermined system)을 구성하며 무한개의 해를 가지게 된다. 그러나 실 음원의 개수 s가 잠재적 단극음원의 개수 N보다 매우 작다면(s ≪ N), 구하고자 하는 해 x∈ℂ^N는 대부분의 원소가 0 또는 0에 가까운 값을 가지는 희소벡터(sparse vector)가 되고, 식 (3)과 같은 l₀-norm 최소화 문제를 풀어 해를 구할(또는 복구할) 수 있다.

여기서 l₀-norm ‖x‖₀은 0이 아닌 원소의 개수를 나타낸다.

그러나 식 (3)은 다항 시간(polynomial-time) 내에 해를 구할 수 없는 NP-hard 문제이기에 사실상 계산이 불가능함이 알려져 있다⁽⁷⁾. 따라서 식 (4)와 같이 볼록 완화(convex relaxation)을 통한 l₁-norm 최소화 문제로 근사하여 해를 추정하게 되며, 이를 기저 추적(basis pursuit)이라고도 한다⁽⁸⁾.

여기서 l₁-norm은 ‖x‖₁=∑|x_n|으로, 원소의 절대값 합을 나타낸다.

식 (4)로 주어지는 l₁-norm 최소화 문제의 의미를 살피기 위해 다음과 같이 통상적으로 사용되는 l₂-norm 최소화 문제와 비교해 보자.

여기서 l₂-norm은 ‖x‖₂=(∑x_n²)^1/2이다.

N = 2, M = 1인 y = Ax 문제에 대한 l₁-norm 최소화 문제(식 (4))와 l₂-norm 최소화 문제(식 (5))의 예를 Fig. 2에 나타내었다. 식 (4), (5)는 목적함수 정의에 사용된 norm에 따라 각각 마름모꼴의 등고선(Fig. 2(a))과 원형의 등고선(Fig. 2(b))으로 나타나며, 제약조건은 빨간색 직선으로 표시된다. 만약 제약조건이 없다면 각각의 최소화 문제는 원점을 해로 가질 것이다. 그러나 식 (4), (5)는 제약조건이 존재하므로, 이를 만족하면서 목적함수가 최소값을 가지는 지점, 즉 가장 작은 값을 가지는 등고선과 제약조건의 접점이 해 x가 된다.

Fig. 2 
Comparison of N = 2 and M = 1

각 최소화 문제의 해를 살펴보자. 식 (5)의 해는 x₁, x₂ 모두 0이 아닌 값을 가졌으며, 이는 l₂-norm 최소화 문제의 해가 희소치가 ‘2’임을 나타낸다. 반면, 식 (4)의 해는 x₁축 위에 존재하므로, l₁-norm 최소화 문제의 해는 희소치 s가 ‘1’이고, l₂-norm 최소화 문제를 통해 얻은 희소치보다 작음을 알 수 있다. 따라서 식 (4)의 l₁-norm 최소화 문제를 이용하면 희소치 s의 최소화가 가능함을 설명할 수 있다.

단, l₁-norm 최소화 문제(식 (4))의 해와 l₀-norm 최소화 문제(식 (3))의 해가 일치하는가 하는 유일성 문제가 발생한다. 이 일치성을 보증하는 방법은 두 가지가 있다. 하나는 잠재적 단극음원의 개수 N과 센서의 개수 M이 식 (6)의 관계식을 만족하는 것이다.

식 (6)의 조건을 충족하면 l₁-norm 최소화 문제의 해와 l₀-norm 최소화 문제의 해가 높은 확률로 일치함이 알려져 있다⁽⁹⁾. 다른 하나는 제한등방성이라는 조건을 만족하는 것으로, 관련 내용은 다음 절에서 설명하겠다.

실제 환경에서의 계측은 잡음을 포함하므로, 이를 고려한 식은 식 (7)과 같다.

여기서 e=[e₁, ..., e_n]^T∈ℂ^N은 서로 상관관계가 없는 원소로 구성된 가우시안(gaussian) 잡음벡터이며, 위 첨자 T는 전치(transpose)를 의미한다. 이에 따라 식 (3), (4) 또한 각각 식 (8), (9)로 수정된다.

여기서 ε은 잡음벡터 e의 l₂-norm 크기를 나타낸다.

따라서 식 (4), (9)의 해를 선형계획법(linear programming)이나 이차계획법(quadratic programming) 등을 통해 구하거나, Matlab CVX toolbox⁽¹⁰⁾와 같은 공개된 코드를 이용하여 구할 수 있다. 그러나 식 (9)의 해는 ε값에 민감하게 반응하기에 과적합(overfitting)과 같은 문제가 발생할 수 있으며, 이 연구에서는 Choo, Y. et al.⁽⁵⁾이 제안한 식 (10)의 제약조건을 가지는 l₁-norm 최소화 문제를 사용함으로써 이 문제를 해결하고자 하였다.

여기서 ε`(=y^Te)은 잡음벡터 e와 관측벡터 y의 내적(inner product) 절대값이다.

식 (10)은 관측벡터 y와 상관관계가 가장 큰 원소 하나를 가지는 해, 다시 말해서 희소치 s가 ‘1’인 x를 구하는 특징이 있으며, ε`값 변화에 추정치가 크게 달라지지 않는 강인성을 가지고 있다.

단, 식 (10)을 비롯하여 앞에서 나타낸 최소화 문제들은 주파수 영역에서 정의되었음을 강조하며, 이에 따라 해석 주파수 대역 선정에 대한 문제가 남는다. 이 연구에서는 파워 스펙트럼 밀도의 변화 관찰 및 스펙트럴 커토시스 해석을 통해 최적 주파수 대역을 선정하였다. 그 후, 선정된 주파수 대역에 대한 식 (10)의 반복 적용 및 비상관 평균화(incoherent averaging)를 통해 최종해, 즉 위치와 강도를 얻는다.

관측행렬 A는 가정된 음원 모형에서 도출되며, 해 x는 A와 밀접한 관련이 있다. RIP 조건은 관측행렬 A가 얼마나 직교행렬에 가까운지 판단하는 척도로, 제한등방상수(restricted isometry constant, RIC) δ_s(0 ≤ δ_s ≤ 1)를 통해 추정할 수 있다^(11,12). δ_s는 s개의 0이 아닌 원소를 가지는 모든 s-희소벡터 x에 대하여 식 (11)을 만족하는 값 가운데 가장 작은 것을 가리킨다.

여기서 A는 각각의 열(column)이 l₂-norm으로 정규화(normalized)되어 있다고 가정한다. 만약 δ_s가 특정한 기준보다 낮으면 A는 s-차수(s-th order)의 RIP 조건을 만족한다고 할 수 있다.

RIP의 핵심 개념은 s개의 열로 구성된 A의 부분집합이 에너지의 손실 또는 이득 없이 정규직교변환(orthogonal transformation)처럼 기능함⁽¹³⁾에 있으며, 일반적으로 Candes⁽¹²⁾가 제시한 δ₂s ≤ 0.414를 RIP 충족을 위한 필요조건으로 본다. 여기서 RIC의 아래 첨자가 2s로 바뀐 것에 주의해야 하며, s-희소벡터 x를 구하기 위해서는 δ_s가 아닌 2배 차수에 대한 δ₂s 가 Candes criteria를 만족하는지 확인해야 함을 강조한다. 이때 δ₂s 는 식 (12)의 δ를 이용하여 통계적으로 추정이 가능하다.

식 (12)는 식 (11)에 δ_s에 대하여 정리한 후, 부등식을 등식으로 수정함으로써 얻을 수 있으며, δ는 몬테카를로법(monte carlo method)⁽¹⁴⁾을 통해 분포를 구할 수 있다. 이때 식 (12)의 x는 s-희소벡터가 아닌 2s-희소벡터이다.

Fig. 3에 δ의 분포 예시를 나타내었으며, 그 분포는 특정한 확률밀도함수 P(δ)로 근사할 수 있다. 만약 근사된 P(δ)가 식 (13)의 정규분포를 따른다면, 식 (14)를 이용하여 제한등방상수 δ_2s를 추정할 수 있다⁽¹⁵⁾.

Fig. 3 
Histogram using 100 000 random draws for the estimation of δ_2s(s=3) (M = 24, N = 441 and f = 8 kHz)

여기서 E(δ)는 δ의 평균이며, σ는 표준편차를 나타낸다.

정리하자면 식 (14)를 통한 δ₂s 의 추정치는 근사된 확률밀도함수의 누적분포가 99.9 %에 해당하는 값임을 알 수 있다. 따라서 식 (13)의 P(δ)가 정규분포로 근사되지 않더라도 그것의 누적분포를 통해 δ₂s를 추정할 수 있다.

제안된 추정법을 캐비테이션 소음 실측 결과에 적용하고자, 선박해양플랜트연구소(KRISO)의 캐비테이션 터널에서 KVLCC2 모형선(Fig. 4(a))을 이용한 초생 캐비테이션 발생 실험을 수행하였다. 날개 수가 4개인 프로펠러가 장착된 모형선을 사용하였으며, 실선 KVLCC2와의 비율 λ는 1:39.44이다.

Fig. 4 
Respect to the frequency 0 kHz ~ 120 kHz

캐비테이션 소음신호 계측에는 9개의 청음센서를 사용하였으며, 설치공간의 제약으로 Fig. 4(a)와 같이 프로펠러 상부 선체 표면에 센서를 매립하였다. Fig. 4(b)에는 RIP 조건 분석(Fig. 4(c))을 통하여 도출한 센서 배열의 좌표를 도시하였다. 프로펠러 중심을 원점으로 한 주변 0.15×0.3×0.3 m³의 영역을 가로, 세로, 높이 10 mm 간격의 15 376개 격자로 분할하여 관측 공간(Fig. 4(a))으로 선정하였고, 이를 바탕으로 M = 9, N = 15 376의 크기(size)를 가지는 관측행렬 A를 RIP 조건 분석에 사용하였다. 캐비테이션의 광대역 특성⁽¹⁶⁾ 및 샘플링 주파수를 고려하여 0 kHz ~ 120 kHz의 주파수 대역에서 수행하였으며, 분석 결과로부터 20 kHz ~ 120 kHz 주파수 범위에서 센서 배열이 추정 가능한 캐비테이션 소음원의 개수는 최대 2개임을 알 수 있었다.

터널 내 압력 조건에 따라 4개의 실험조건(condition 1 ~ 4)을 선정하였고, 각각의 압력 조건은 Table 1과 같다. 각 조건에 해당하는 캐비테이션 수(cavitation number, σ_n)는 식 (15)로 정의된다.

여기서 n [Hz]은 프로펠러의 초당 회전수, D [m]는 프로펠러의 직경, P [Pa]는 터널 내 압력, P_n[Pa]은 포화수증기압을 나타낸다.

Fig. 5는 고속카메라를 이용한 실험조건별 캐비테이션 관찰 결과로, 선수에서 선미를 본 방향으로 촬영하였다. Condition 1은 육안으로 거의 관측되지 않을 만큼 캐비테이션이 없는 비공동 상황을 나타내며(Fig. 5(a)), condition 2에서는 Fig. 5(b)와 같이 프로펠러 날개끝 부근에서 보텍스 캐비테이션이 나타났으나 간헐적 관찰만 가능하였다. 이때 보텍스 캐비테이션은 날개 끝단에서 생성되지 않고 프로펠러 날개가 12시 방향을 통과한 후 반시계방향으로 20° ~ 30° 회전했을 때 불규칙하게 발생하였으며, 길게 늘어지지 않고 짧게 나타났다가 바로 사라지는 경향을 띠었다. Condition 3, 4에서는 공통적으로 날개끝에서 길게 늘어진 형태의 보텍스 캐비테이션과 날개 끝단을 덮는 형태의 시트 캐비테이션(sheet cavitation)이 관측되었다(Fig. 5(c), (d)).

Fig. 5 
Cavitation pattern for different test conditions (captured using high-speed camera)

Fig. 6에는 계측된 소음신호 중 음압 레벨이 비교적 양호한 F2 채널의 실험 조건별 파워 스펙트럼 밀도를 도시하였으며, 주파수 분해능 64 Hz 및 해닝 창문함수(hanning window function)을 적용하여 75 %의 오버랩(overlapping)으로 641회의 앙상블 평균(ensemble average)을 취해 계산하였다.

Fig. 6 
Power spectral density of cavitation noise for different test conditions (channel: F2)

Fig. 5에서 터널 내 압력이 감소할수록 더 많은 캐비테이션 관측이 가능했던 점을 미루어, Fig. 6의 파워 스펙트럼 밀도 역시 눈에 띄는 차이가 있을 것으로 예상하였다. 그러나 condition 1에서 4로 터널 내 압력이 감소해도 음압 레벨(sound pressure level)의 미세한 증가만 있을 뿐이었고, 20 kHz ~ 60 kHz 대역과 100 kHz ~ 120 kHz 대역에서 변조 및 고조파 성분이 발생하여 낮은 신호 대 잡음비(signal to noise ratio, SNR)의 신호가 나타났다. 즉, 2장에서 언급한 것처럼 주파수 대역의 선정이 필요하나, 파워 스펙트럼 밀도의 변화 관찰만으로는 판단하기 어려움을 알 수 있다. 따라서 이 연구에서는 서론에서 언급한 바와 같이, 날개끝 보텍스 캐비테이션 소음이 순환정상성 신호로 분류될 수 있음에 주목하여 스펙트럴 커토시스 해석을 수행하였고, Fig. 7에 결과를 도시하였다. 초생 캐비테이션을 간헐적으로 관측할 수 있었던 condition 2의 시계열 자료를 사용하였으며, 해석 결과로부터 64 kHz ~ 81 kHz를 최적 주파수 대역으로 선정하였다.

Fig. 7 
Spectral kurtosis analysis of condition 2

캐비테이션은 일반적으로 생성과 소멸을 거치면서 짧은 시간 동안 큰 압력을 발생시키므로, 이때 수반되는 충격 신호를 캐비테이션 사건으로 볼 수 있다. 이에 이 연구에서는 Fig. 8과 같이 파워 스펙트럼 밀도 해석에 사용된 시계열 자료에서 캐비테이션 사건이 극히 일부분을 차지하는 점을 미루어, 개별 캐비테이션 사건의 시계열 자료를 사용하면 분석에 불필요한 잡음 성분이 억제 또는 제거되어 SNR 향상이 가능할 것으로 판단하였다. 캐비테이션 대부분이 프로펠러 날개가 전면 기준 반시계방향으로 약 20°~30° (단, 12시 방향을 0°로 시작) 지점을 통과할 때 관측되었음(Fig. 5)을 바탕으로 캐비테이션 사건의 시작-종료 사이의 지속 시간을 약 3 ms로 간주하였다. 이 지속 시간에 대하여 충격 신호를 포함하는 시계열 자료를 Fig. 9(a)와 같이 임의 선택하여 Fig. 9(b)에 도시하였고, Table 2에 시험 조건별 선택한 캐비테이션 사건 개수를 나타내었다. Fig. 5의 관찰 결과에서 알 수 있듯, condition 1에서 4로 터널 내 압력이 감소할수록 개별 캐비테이션 사건 개수는 증가하는 추세를 보였다.

Fig. 8 
Time series data of cavitation noise for different test conditions (channel: F2)

Fig. 9 
Example of selected individual cavitation event

Fig. 10에는 Fig. 9(b)에 도시한 캐비테이션 사건에 대한 파워 스펙트럼 밀도를 나타내었으며, 시계열 자료의 과도 신호(transient signal) 특성과 길이를 고려하여 주파수 분해능 333.3 Hz 및 사각 창문함수(rectangular window function)를 적용하였다. SNR은 전 주파수 대역에서 Fig. 6 대비 약 10 dB 정도 상승하였음을 알 수 있는데, 이는 광범위하게 나타나던 변조 및 고조파 성분이 캐비테이션 소음신호에 대부분 가려졌기(masking) 때문으로 판단된다.

Fig. 10 
Power spectral density of individual cavitation of Fig. 9(b)

개별 캐비테이션 소음의 발생위치 및 강도 추정을 위해 스펙트럴 커토시스 해석을 통해 선정한 65 kHz ~ 81 kHz의 관심 주파수 대역(Fig. 10의 흑색 점선)을 적용하였고, 333.3 Hz의 해상도로 분할하여 관측벡터 y를 추출하였다. 앞서 결정된 관측행렬 A와 함께 식 (10)을 관심 주파수 대역에 대해 반복 적용하였으며, 추정에는 CVX toolbox⁽¹⁰⁾를 사용하였다. 이때 식 (10)의 ε`는 주파수와 무관하게 0.01로 일정한 값을 설정하였다. 잡음의 성질을 정의할 수 없기 때문에 오차가 존재할 가능성이 다분하나 여러 ε`값에 대한 계산을 반복하였고, 0.01 이하의 ε`에서 추정치에 대한 오차는 최대 0.5 % 미만으로 큰 변화가 없음을 사전에 확인하였다.

Fig. 11에는 개별 캐비테이션 사건에 대한 시험 조건별 추정 결과를 나타내었다. 각각 위에서 아래를 바라본 방향(x-y 평면)과 선수에서 선미를 바라본 방향(y-z 평면)의 그림을 도시하였고, 각 그림 내부의 흑색 선과 원은 프로펠러 디스크를 나타낸다. 또한, 회색조의 점으로 캐비테이션 소음의 강도를 나타내었으며, 시험 조건별 강도의 최대값을 범위로 각각 지정하였다. 추정 위치는 프로펠러의 상부, 다시 말해서 12시 방향에 집중적으로 나타났고, 시험 조건과 상관없이 일정한 영역에서 군집을 이루었으며, 그 범위는 고속카메라 관찰 결과(Fig. 5)와 큰 상관관계를 띠고 있었다. 추정 강도는 Fig. 8의 음압 수준을 거리 역산한 추정 강도와 유사하게 나타났다. 여기서 Fig. 11(a), (b)의 추정 강도 범위가 약 20 kPa⋅m로 유사하였는데, 이를 통해 육안 관찰에서는 비공동 상황이었던 condition 1의 터널 내 압력에서도 condition 2와 같은 날개끝 보텍스 캐비테이션이 발생하였음을 알 수 있었다.

Fig. 11 
The results of localization- and source strength estimation of individual cavitations for different test conditions

Fig. 12는 개별 캐비테이션 사건에 대한 추정 강도를 dB 단위로 나타낸 결과이다(단, dB의 기준값은 추정 강도의 최대치를 이용하였음). 추정 강도는 각 시험 조건별로 대략 40 dB ~ 50 dB 정도의 산포를 가지나, 평균값의 추이는 터널 내 압력 감소에 따라, 다시 말해서 캐비테이션 양이 증가함에 따라 비례하는 특성을 보였다. 따라서 이 연구에서 제안하는 기법을 통해 캐비테이션 위치 및 강도의 정량적 추정이 가능할 것으로 기대한다.

Fig. 12 
The dB distribution and mean value of the source-strength of individual cavitations for different test conditions