I. 서 론

고강도 집속형 초음파나 충격파와 같은 고출력 초음파를 인체와 같은 매질에 조사할 때 수동 공동이 발생될 수 있다. 초음파가 진행하는 매질의 압축과 인장을 반복하는 초음파의 특성에 의해 매질에 조사된 초음파는 압력 변동을 유발하며, 주로 음압에 의해 기포군이 발생한다. 이렇게 발생된 기포는 어느 시점에서 관성적으로 급격히 수축하여 파열하면서 수동 초음파 공동 현상을 일으킨다. 기포의 파열에 의한 초음파 공동은 이차적인 초음파 또는 충격파를 발생시킨다.

이와 같은 고강도 초음파에 의한 기포의 발생, 기포의 동적 거동 및 관성적인 파열 등의 일련의 과정에서 방출되는 이차적인 초음파는 다양한 초음파의 효과(세포 파괴, 결석 분쇄 등)를 유발하기 때문에, 초음파 공동 현상으로 유발된 초음파 신호를 관찰하고 추적하는 것은 매우 중요하다.^[1],^[2],^[3],^[4],^[5],^[6] 그러나 수동 공동을 영상화할 때 초음파 영상을 위한 송신을 하지 않고 수신 집속 방법만 적용 가능하다. 배열 트랜스듀서에서 일반적으로 많이 사용되는 지연 후 더해주는(Delay-And-Sum, DAS) 빔포밍 방법은 수신 과정에서만 적용되므로 낮은 해상도와 영상 영역에 넓게 퍼지는 큰 부엽에 의하여 수동 공동을 관찰하기 어려운 문제점이 있다. 초음파 의료 영상 시스템의 집속 과정은 Fig. 1과 같다. 수신 집속 시스템에서, 수신된 각각의 채널의 신호에 수신 집속 시간 지연을 인가하면 영상점(Fig. 1의 a 점)에서 돌아온 신호는 모든 수신 채널에 동시에 도착하는 것처럼 시간이 정렬된다. 따라서 모든 수신 채널의 신호를 더하면 큰 신호를 얻는다. 그러나 영상점 밖(Fig. 1의 b 점)에서 돌아온 신호는 모든 수신 채널에 동시에 도착하지 않게 나타나므로 집속 시간을 인가한 후에도 시간축에서 정렬되지 않는다. 따라서 모든 채널 신호를 더하면 영상점보다 작은 신호가 나타난다. 영상점 밖의 신호는 집속 과정에서 완전히 제거되어야 하지만, 제거되지 않고 남으면 부엽이 되어 영상의 해상도와 대조도를 저하시킨다. 초음파 의료 영상 시스템에서 부엽을 억제하기 위하여 널리 사용되는 방법은 절족화(apodization) 방법이다. 이 방법은 해밍(Hamming) 윈도우와 같은 창함수를 송수신 채널 데이터에 곱하여 부엽을 억제시키지만 해상도를 저하시키는 단점이 있다.^[7],^[8],^[9]

Fig. 1.

Architecture of delay-and-sum beamforming system using array transducer.

이를 극복하기 위하여 모든 영상점마다 최적의 절족 윈도우를 계산하여 적용하는 최소 분산 빔포밍(Minimum Variance Beamforming, MVB)이 의료용 초음파 영상에 적용되고 있다.^[10],^[11] 이 방법은 해상도를 향상시키면서 부엽 억제 효과가 크므로 많은 연구가 이루어지고 있다. 수동 공동 영상(Passive Cavitation Imaging, PCI)의 경우 송신 집속을 적용할 수 없으므로 부엽의 영향이 더 크게 나타난다. 기포의 파열에서 발생하는 초음파 신호는 하나의 점발생원에서 발생하는 신호로 가정할 수 있으므로, 특히 MVB는 수동 공동 영상의 해상도 향상에 효과적으로 적용될 수 있다.

II. 수동 공동 영상에서 집속 신호처리

DAS 집속 방법에서 t 시간에 영상점 px에서 수신 집속된 신호 s_px(t)는 다음과 같다.

여기서 w_n는 절족화 윈도우이며 n번째 채널의 수신신호에 수신 집속 시간 지연을 인가한다. N_ch는 수신 채널의 개수이다. MVB에서는 최적의 절족화 윈도우를 각각 영상점마다 계산하여 적용한다.

초음파 의료 영상 시스템에서 축방향 해상도를 얻기 위하여 시간적으로 짧은 펄스 형태의 초음파를 송수신한다. 그러나 PCI에서는 기포의 파열에서 연속적으로 발생하는 초음파를 수신 집속만하여 영상을 구성하므로 축방향 공간 해상도가 떨어질 뿐 아니라 측방향으로 부엽의 영향이 크게 나타난다. 축방향 해상도는 수동 공동에서 발생하는 초음파의 주파수 스펙트럼에 의하여 결정되며, 부엽의 경우 집속 방법에 의존한다. 따라서 부엽의 영향을 줄이기 위한 방법이 필요하다. 본 논문에서는 수신 집속 과정에서 영상점에서 채널 신호 크기의 분포를 나타내는 파라미터인 무게중심과 평탄도를 계산하여 부엽의 영향을 추정하였다.

III. 수신 채널 데이터를 이용한 무게중심과 평탄도의 정의

초음파 수신 집속 시스템에서 수신 채널의 신호에 집속 지연을 인가하여 진행 거리의 시간 차이를 보상하여 정렬시킨다. Fig. 2(a)는 영상점에서 돌아온 초음파 신호에 집속 지연을 인가한 후의 수신 채널의 신호 모양이며, 모든 채널에 같은 시간에 수신되므로 전체 채널에 걸쳐 균일한 크기의 신호가 나타난다. Fig. 2(b)는 영상점 밖에서 돌아온 초음파 신호이다. 초음파 펄스 신호의 시간 길이가 짧은 경우 입사각에 따라 전체 수신 채널에서 각각 다른 시간에 초음파 신호가 도착한다. 따라서 영상점과 영상점 밖에서 온 신호는 신호가 돌아오는 입사방향에 따라서 수신 채널에 나타나는 신호의 분포가 다르므로, 모든 수신 채널에서 수신 신호가 존재하는 분포를 계산하면 영상점과 영상점 밖의 신호를 구별할 수 있게 된다. 수신 채널에 나타나는 신호의 분포를 나타내는 파라미터로 무게중심과 평탄도를 계산하였다. 무게중심 Centroid는 전체 수신 채널에서 각각의 수신 채널 신호 크기의 무게중심의 위치로 계산된다.

여기서 r(ch)는 ch번째 채널에 수신된 신호이며 그것의 절댓값인 $\left|r\left(ch\right)\right|$를 사용한다. N_ch는 전체 수신 채널의 개수이다. Fig. 2(a)의 경우 모든 채널에서 신호 크기가 같으므로 무게중심을 구하면 수신 채널에서 가운데 채널의 위치에 나타난다. 그러나 Fig. 2(b)의 부엽 신호의 경우 수신 채널의 중앙에서 벗어난 위치에 무게중심이 위치한다. 따라서 무게중심이 수신 채널의 가운데 채널에 가까울수록 영상점에서 돌아온 신호의 가능성이 높다.

Fig. 2.

Waveforms of echoes across transducer array: (a) from imaging point and (b) not from imaging point.

평탄도는 모든 수신 채널 신호 크기의 분산으로부터 구해지며 1에 가까울수록 모든 채널에 비슷한 크기의 신호가 수신되었음을 의미한다. 영상점에서의 평탄도 f(x,z)는 다음과 같이 정의한다.

여기서 var(r(ch))는 영상점 (x,z)에서 수신 채널 데이터의 분산값이다. 한 프레임의 모든 영상점에서 분산값을 모두 구한 뒤에 최댓값으로 정규화하여 0≤ f(x,z) ≤ 1의 값을 가지도록 하였다.

무게중심을 이용하여 영상의 부엽을 억제할 때는 평탄도의 위치에 따라 가중값으로 변환하여 사용한다. Fig. 3은 64채널 시스템에서 평탄도에 따른 가중값의 예이다. 실선은 삼각형 가중값이며, 파선은 Hamming 윈도우의 4제곱을 사용하여 만든 가중값이다. 채널 데이터에서 계산한 평탄도가 중앙의 채널에 위치하면 평탄도 가중값은 1, 채널의 양 끝에 있으면 0의 값을 가지도록 정규화하였다. 따라서 가중값이 1에 가까우면 수신 채널 신호가 영상점에서 돌아온 주엽 신호를 의미하며, 1보다 작아지면 영상점 밖에서 돌아온 부엽 신호를 의미한다. 초음파 영상에서 무게중심 가중값과 평탄도를 이용하여 부엽을 억제하는 방법은 다음과 같다.

여기서 Pixel(x,z)는 입력 영상값, α는 무게중심과 평탄도를 이용한 가중값, Pixel_weighting(x,z)은 가중값 처리 후의 영상값이다. 가중값 α를 구하는 한 가지 방법은 무게중심 가중값과 평탄도를 서로 곱하는 것이다. 인체와 같이 많은 반사체가 존재하는 스페클 영상에서는 모든 수신 채널에 비슷한 크기의 신호가 수신되므로, 무게중심을 이용하여 주엽과 부엽 신호를 구별하기 어렵다. 그러나 제안한 방법은 신호의 발생원이 점발생원인 경우 부엽을 억제하는 데 유용하다. 기포의 파열에서 발생하는 초음파 신호는 점발생원으로 가정할 수 있으므로 제안한 방법을 유용하게 적용할 수 있다.

Fig. 3.

Weighting values as a function of centroid in a 64 channel imaging system: triangular weight (solid) and weight corresponding to the fourth power of the Hamming window (dashed).

IV. 수동 공동 영상의 컴퓨터 시뮬레이션

다음의 시뮬레이션은 공동에 의하여 발생한 초음파를 이용하여 공동을 영상화하였다. 시뮬레이션 조건은 공동에서 발생하는 신호를 점발생원으로 모델링하였고, 중심주파수 5 MHz를 가지는 펄스 형태의 초음파 신호가 30 mm 깊이에서 발생한다고 가정하였다. 초음파 신호의 수신은 5 MHz 중심주파수를 가지는 64채널의 선형 트랜스듀서를 사용하고, 소자의 폭은 0.6 mm이며, 수신 신호는 80 MHz의 샘플링 주파수로 샘플링하여 무선주파수(Radio Frequency, RF) 데이터를 생성하였다. 프로그램은 MATLAB으로 작성하고 IBM PC에서 실행하였다. Fig. 4는 한 점의 공동에서 돌아온 RF의 수신 채널 신호의 모양이며 발생원의 위치에 따라 신호의 분포가 원호 모양을 나타낸다.

Fig. 4.

Temporal evolution of channel RF waveform due to cavitation.

PCI의 구현은 DAS 빔포밍과 MVB를 적용하여 공동군을 영상화하였다. Fig. 5(a)는 DAS로 얻은 영상이며, Fig. 5(b)는 MVB로 얻은 영상이다. 영상의 밝기는 MATLAB에서 그레이 컬러맵을 사용하여 사상하였고 DAS 영상은 선형 밝기로 표시하였고, MVB는 부엽이 DAS보다 더 적은 영역에서 작은 크기로 나타나므로 관찰을 위하여 영상의 밝기를 증가시켜서 40 dB 로그 압축을 적용하였다. 수신 집속만으로 공간 해상도를 얻어야 하므로 DAS 영상은 넓은 영역에 걸쳐 부엽이 나타난다. 의료용 초음파 영상 시스템에서와 같이 펄스 신호의 시간 게이팅에 의한 축방향 해상도를 얻을 수 없으므로 부엽은 축방향으로 길게 나타난다. RF 데이터가 점발생원이므로 MVB 영상의 경우 공간 해상도가 매우 좋게 나타난다. MVB의 계산은 Reference [11]에서 서술한 수식을 사용하여 64채널 시스템에서 공간 평활화를 적용하지 않고, 대각선 하중 계수는 0.01로 적용하여 계산하였다.

Fig. 5.

Images of a point like bubble obtained using (a) DAS without compression and (b) MVB with 40 dB logarithmic compression.

Fig. 6(a)는 Fig. 5의 영상 영역에서 계산한 무게중심 영상이며 Fig. 3의 Hamming 윈도우 값의 4제곱을 사용하였다. 수동공동이 있는 위치에서 축방향으로 큰 값이 나타나서 축방향으로의 해상도 향상은 기대할 수 없다. 그러나 부엽이 나타나는 ‘X’ 형태의 영역에는 작은 가중값을 나타낸다. Fig. 6(b)는 평탄도의 모습인데 수동 공동이 있는 영상점 주위에서 큰 값이 나타난다. Fig. 6(c)는 Fig. 6(a)와 Fig. 6(b)를 서로 곱하여 얻은 가중값이며, Fig. 6(b)에 비하여 부엽 영역에서 더 작은 값을 가진다. Fig. 6(c)의 가중값을 사용하여 영상에서 부엽을 억제하는 신호처리에 사용된다.

Fig. 6.

2-D maps of (a) centroid weighting, (b) flatness, and (c) multiplication of centroid weighting and flatness.

Fig. 7은 DAS와 MVB 영상에 가중값을 곱하여 얻었으며, (a)는 weighted DAS 영상이고 (b)는 weighted MVB 영상이다. 부엽을 관찰하기 위하여 두 영상 모두 40 dB 로그 압축을 적용하였으며 부엽이 나타나는 영상 영역이 줄어들었다. Figs. 8과 9에서는 측방향 음장 특성과 축방향 음장 특성을 정량적으로 비교하였다. 실선은 DAS, 점선은 weighted DAS, 파선은 MVB, 일점쇄선은 weighted MVB이다. 가중값 처리된 영상에서 측방향 음장에서 부엽이 10 dB 정도 감소하고 있으며 축방향 음장의 경우에도 영상점에서 멀어질수록 부엽이 많이 감소한다.

Fig. 7.

Images obtained using (a) weighted DAS and (b) weighted MVB, both with a 40 dB log compression.

Fig. 8.

Lateral field responses: DAS (solid), weighted DAS (dotted), MVB (dashed), and weighted MVB (dash-dotted).

Fig. 9.

Axial field responses: DAS (solid), weighted DAS (dotted), MVB (dashed), and weighted MVB (dash-dotted).

V. 실 험

Fig. 10은 집속형 초음파에 의해 발생된 기포의 파열에 의해 유도된 초음파 분포를 영상화하는 수동 공동 영상을 얻기 위해 구성한 실험 장치의 구성을 보여준다. 집속형 초음파 발생을 위한 초음파 트랜스듀서는 물을 저장하는 수조의 한쪽 벽에 고정되며, 초음파 신호를 수신할 64채널 선형 배열의 영상 프로브는 초음파 트랜스듀서의 축 방향과 수직을 이루도록 32 mm의 거리에 설치된다. Fig. 11(a)는 충격파 출력용 트랜스듀서와 영상용 트랜스듀서의 배치를 나타낸다. Fig. 11(b)는 초음파 전달 매질인 물속에서 집속형 초음파에 의해 발생한 기포군에 대한 초음파 B-모드 영상을 보여주며, B-모드 영상의 크기는 60 mm × 38.4 mm이다. 고출력의 초음파를 발생시키는 장비는 Shinewave-Sonic(HnT Medical, Korea)을 사용하여 충격파를 생성시켰다. 영상 데이터 획득은 EQB-12R(Alpinion, Korea)에서 5 MHz 중심주파수, 0.3 mm 피치, 128채널 선형 프로브에서 소자를 1개씩 건너뛰어서 전체 64개의 소자를 사용하였다. 따라서 배열 소자의 피치를 0.60 mm로 하여 64채널 데이터를 얻었다. RF 데이터는 40 MHz의 샘플링 주파수로 얻어서 메모리에 저장한 다음 PC로 옮겨서 2배 보간하여 80 MHz의 샘플링 주파수로 변환한 뒤에 신호처리를 하여 영상을 얻었다. 신호처리에 사용한 프로그램은 MATLAB으로 작성하고 IBM PC에서 계산하였다.

Fig. 10.

Experimental setup.

Fig. 11.

Experimental setup: (a) relative placement of shock wave generation and imaging probes for alignment and (b) B-mode image of cavitation due to shock wave.

집속형 초음파에 의해 초점 영역에서 발생된 기포는 시간이 지나면서 점차 커지다가 어느 순간 급속히 빠른 속도로 관성적으로 수축하여 파열되는 것이 초음파 영상에서 확인되었으며, 이 과정에서 수반된 기포의 동적 거동으로 이차적인 초음파 또는 충격파가 방출되는 것이 확인되었다. Fig. 12는 기포 붕괴에서 발생한 공동에 의한 초음파를 수신한 RF 데이터의 일부이다. 공동 신호는 약 100 μs의 시간 구간 동안 주로 발생하였다. 이 RF 데이터를 이용하여 DAS, MVB, weighted DAS, weighted MVB 영상을 계산하였다. 영상은 깊이 정보 없이 연속적으로 수신되므로, 영상을 얻고자 하는 시간에 수신된 채널 데이터에서 영상점에 대응되는 수신 집속 시간 지연을 인가하여 수신 집속을 하여 얻은 신호의 크기를 영상점에 사상하였다. 영상의 크기는 30 mm 초점 깊이를 중심으로 30 mm × 40 mm 크기의 영상을 계산하고 영상의 밝기를 MATLAB에서 그레이 맵을 이용하여 표시하였다. 영상의 노이즈를 줄이기 위하여 RF 데이터에 패스밴드가 3 MHz ~ 8 MHz인 대역통과필터를 적용하였다. Fig. 13에서 (a)는 한 개의 공동에 의한 RF 데이터, (b)는 DAS 영상, (c)는 MVB 영상이다. (d)는 두 개의 공동 신호가 0.25 μs 시간 차이를 가지고 수신된 RF 데이터를 이용하여 영상을 합성하였고 (e)는 DAS 영상, (f)는 MVB 영상이다. 두 경우 모두 MVB 영상의 부엽이 줄어들었으며, MVB의 경우, 실험 데이터의 노이즈 때문에 시뮬레이션에서보다 부엽의 억제 효과가 조금 떨어진다. 그러나 두 개의 공동의 경우 다른 시간에 수신된 신호이므로 RF 데이터에서 흰색 가로선으로 표시한 시간에 구한 영상에서 하나의 공동이 더 뚜렷하게 영상에 나타났고 다른 하나는 밝기가 줄어들었다.

Fig. 12.

Part of 64 channels of RF data.

Fig. 13.

(a) RF echo from single bubble collapse, (b) DAS image, (c) MVB image, (d) RF echo from two adjacent bubble collapse, (e) DAS image, and (f) MVB image.

Fig. 14에서 (a)와 (b)는 단일 기포에서의 무게중심 가중값과 평탄도이며, (c)는 이에 따른 가중값이다. (d)와 (e)는 각각 두 개의 기포에서의 무게중심 가중값과 평탄도이며, (f)는 이에 따른 가중값이다. 공동이 있는 주위의 부엽이 있는 영역에서 가중값은 작은 값을 가진다.

Fig. 14.

Plots of (a) centroid weighting from single bubble collapse data, (b) uniformity of single bubble collapse image, (c) weighting factor obtained by multiplying together centroid weighting and flatness of single bubble collapse image, (d) centroid weighting from two adjacent bubble collapse data, (e) flatness of two bubble collapse image, and (f) weighting factor obtained by multiplying together centroid weighting and flatness of two bubble collapse image.

Fig. 15는 Fig. 13의 영상에 Fig. 14의 가중값을 곱한 것이다. Fig. 15에서 (a)와 (b)는 각각 단일 공동 영상에서 가중값을 곱한 DAS와 MVB이며 (c)와 (d)는 각각 두 개의 공동이 있는 영상에서 가중값을 곱한 DAS와 MVB이다. 공동을 구별하는 해상도는 모두 향상되고 있지만 배경의 노이즈가 증가한다. 특히 두 개의 공동이 가중값 영상에서 잘 구별되고 있다. 가중값 영상의 경우 수신 채널에 밝은 주엽 신호가 없으면 오히려 배경 노이즈를 증가시키는 경향이 있다. 수동 공동 영상에서 공동의 유무가 중요한 정보이므로, 일정 크기 이하의 영상 신호를 제거하는 문턱값 처리(thresholding)를 함으로써 배경의 노이즈를 억제할 수 있다.

Fig. 15.

B-mode images of (a) weighted single bubble collapse using DAS, (b) weighted single bubble collapse using MVB, (c) weighted two bubble collapse using DAS, and (d) weighted two bubble collapse using MVB.

VI. 결 론

DAS 방법으로 수신 집속을 한 수동 공동 영상은 낮은 해상도와 큰 부엽 때문에 작은 신호 크기의 공동을 구별하기 어려운 문제가 있다. 주엽과 부엽 신호의 채널 데이터의 분포 특성을 이용하여 무게중심과 평탄도를 정의하고 이 값을 이용한 가중값을 영상에 적용하였다. 무게중심은 주변의 노이즈를 증폭하는 문제점이있지만 평탄도와 결합하여 주엽에 가까운 영역에서만 큰 가중값을 가지도록 하였다. 실험을 통해 한 개와 두 개의 공동이 나타나는 영상 데이터에 제안한 방법의 가중값을 곱하여 부엽을 줄이고 동시에 해상도가 증가함을 확인하였다.