Passive Acoustic Mapping Using Compressive Sensing

Abstract

a processing means configured to: define a sample period and a sample space in which a signal can be sampled at each of a plurality of sample points within the sample space; define a grid of candidate bubble locations; define a random projection; sample each of the output signals at groups of sample points defined by the random projection; define a dictionary defining sample values ​​to be obtained for each signal of each group of sensor elements; define a vector having elements that identify candidate bubble locations; and perform a minimization operation to derive vector element values.

Description

[Technical Field]

[0001] The present invention relates to passive acoustic mapping, and in particular to passive ultrasound imaging of cavitation, for example in medical ultrasound systems. [Background technology]

[0002] Real-time treatment monitoring of therapeutic ultrasound—determining both the spatial extent and dose of therapy—is a key issue for widespread clinical adoption, whether the application is high-intensity focused ultrasound (HIFU) ablation of cancerous tissue or ultrasound-mediated targeted drug delivery. One method for achieving this is passive acoustic mapping (PAM), also known as passive cavitation mapping (PCM) or passive cavitation imaging (PCI). PAM is a technique that allows imaging of bubble activity without knowing the time-of-flight of the interrogating pulse that causes the bubble activity. This technique has been used in various applications of therapeutic ultrasound, including cavitation-mediated drug delivery, soft tissue thermal ablation, and transcranial therapy monitoring. Examples of such systems are described in WO 2010 / 052494 A1, which describes a basic method for PAM, and in WO 2014 / 041370 A1, which describes a method for PAM using a sparse transducer array.

[0003] Since its initial development, PAM has been extensively studied, and techniques such as optimal beamforming, planar projection, and non-Gaussian signal processing have been demonstrated to improve resolution or processing speed. However, while hardware developments such as high-channel-count ultrasound research platforms offer unprecedented data acquisition capabilities, real-time processing of growing data streams poses developmental challenges.

[0004] Compressive sensing relies on the fact that many real-world signals are sparse in some domains and therefore can be sampled at rates significantly lower than the Nyquist rate and then reconstructed by an optimization process (R. G. Baraniuk, "Compressive Sensing [Lecture Notes]," IEEE Signal Processing Magazine, Vol. 24, No. 4, pp. 118-121, July 2007, doi:10.1109 / MSP.2007.4286571). This technique has been widely studied for applications such as radar (R. Baraniuk and P. Steeghs, "Compressive Radar Imaging," 2007 IEEE Radar Conference, Apr. 2007, pp. 128-133, doi:10.1109 / RADAR.2007.374203.), MRI (M. Lustig, DL Donoho, JM Santos, and JM Pauly, "Compressed Sensing MRI," IEEE Signal Processing Magazine, vol. 25, No. 2, pp. 72-82, Mar. 2008, doi:10.1109 / MSP.2007.914728), and conventional diagnostic ultrasound imaging.Previous work has focused on matching and basis tracking (M. Gyoengy and C.M. Coviello, "Passive cavitation mapping with temporal sparsity constraint," The Journal of the Acoustical Society of America, vol. 130, No. 5, pp. 3489-3497, Nov. 2011, doi:10.1121 / 1.3626138.), common array processing (C.M. Coviello, R.J. Kozick, A. Hurrell, P.P. Smith, and C. Coussios, "Thin-film sparse boundary array design for passive acoustic mapping during ultrasound therapy," IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 59, No. 10, pp. 2322-2330, Oct. 2012, doi:10.1121 / 1.3626138.) doi:10.1109 / TUFFC.2012.2457.), and PAM based on sparse reconstruction of radio frequency data and sparse selection of array elements (C. Crake et al., "Passive acoustic mapping and B-mode ultrasound imaging utilizing compressed sensing for real-time monitoring of cavitation-enhanced drug delivery," The Journal of the Acoustical Society of America, vol. 143, No. 3, Art. No. 3, Mar. 2018, doi:10.1121 / 1.5036143). Summary of the Invention [Problem to be solved by the invention]

[0005] The present invention provides a system in which PAM operates directly on sparsely sampled data, sometimes referred to as compressed sensing PAM (CS-PAM). [Means for solving the problem]

[0006] The present invention provides a passive compressional wave imaging system for locating cavitation bubbles, comprising a plurality of sensor elements arranged in an array, each sensor element configured to generate an output signal; and processing means, wherein the processing means is configured to: define a sample period and a sample space in which a signal can be sampled at each of a plurality of sample points in the sample space; define a grid of candidate bubble locations; define a random projection for each of the output signals that identifies a different group of sample points; sample each of the output signals at the groups of sample points defined by the random projection to generate sample data; define a dictionary that defines sample values ​​to be obtained for the signal of each of the sensor elements at each of the respective groups of sample points over the sample period for a bubble at each of the candidate bubble locations; define a vector having elements that identify the candidate bubble locations; and perform a minimization operation to derive the vector element values, and therefore the bubble locations, from the sample data, basis, and random projection.

[0007] The vector defining the candidate bubble locations may be the same for all sensor elements or may be defined individually for each of the sensor elements. The random projection can be defined by at least one random projection matrix, which can be one random projection matrix defined as multiple sub-matrices, or two or more random projection matrices each defining at least one compression stage.

[0008] The sample period may be the time it takes for one still image to be generated. If the system is configured to generate video images, the sample period is determined by the frame rate of the video.

[0009] The sample space may be time.The processing means may be further configured to define a sample rate, which defines the sample points as a number of sample instants within a sample period at which samples may be taken.

[0010] The sample space may be frequency. The processing means may be further configured to define the sample points as a plurality of sample frequencies at which samples may be taken. The processing means may be configured to separate an analog component of the sensor output signal at each of the frequency sample points and sample each of the components over a sample period.

[0011] If the sample space is time, the processing means may be arranged to sample each of the sensor output signals and determine frequency component samples for each of the output signals from the time domain samples.

[0012] The processing means may be configured to form a linear combination of the sample data from the sensors to thereby form virtual sensor data corresponding to each of a plurality of virtual sensors. The random projection matrix may be further configured to define a linear combination. Alternatively, the linear combination may be defined by a further matrix.

[0013] The processing means may be configured to identify one of the sensor elements as a reference sensor element. The processing means may be configured to sample the output signal from the reference sensor element at each of the sample points to generate reference data. The processing means may be configured to define a delay time for each combination of one of the candidate bubble locations and one of the sensor elements. The processing means may be configured to determine a dictionary from the reference data and the delay times.

[0014] The present invention further provides a method for passive compressional wave imaging of cavitation bubbles, comprising the steps of defining a sample period and a plurality of sample points across a sample space; defining a grid of candidate bubble locations; defining a random projection for each of the output signals that identifies a different group of sample points; sampling each of the output signals at the groups of sample points defined by the random projection to generate sample data; defining a dictionary that defines sample values ​​to be obtained for each signal of each group of sensor elements at each of the sample points across the sample period for a bubble at each of the candidate bubble locations; defining a vector having elements that identify the candidate bubble locations; and performing a minimization operation to derive the vector element values, and therefore the bubble locations, from the sample data, basis, and random projection.

[0015] Thus, the present invention may provide a compressed sensing (CS) technique that can be used to perform PAM using sparsely sampled data, which can improve image quality and processing speed for a given dataset.

[0016] The system may further comprise any one or more features, in any combination, of the preferred embodiments of the invention described herein, by way of example only, with reference to the accompanying drawings. [Brief explanation of the drawings]

[0017] [Figure 1] 1 is a schematic diagram of a sensor array forming part of an embodiment of the present invention; [Figure 2] 1 is a schematic diagram of an imaging system according to an embodiment of the present invention; [Figure 3] FIG. 3 illustrates the generation of a plane acoustic wave from the system of FIG. 2. [Figure 4] FIG. 3 illustrates the generation of focused acoustic waves from the system of FIG. 2. [Figure 5] FIG. 10 is a functional block diagram of a system in accordance with another embodiment of the present invention in which sensor signals are sampled in the time domain, transformed to the frequency domain, and then compressed. [Figure 6] FIG. 10 is a functional block diagram of a system of a further embodiment of the present invention in which the sensor signal is sampled in the frequency domain. [Figure 7] FIG. 1 is a schematic diagram of an embodiment of the present invention having a mechanically controlled transducer array. [Figure 8] FIG. 1 is a schematic diagram of one embodiment of the present invention having a handheld transducer array. [Figure 9] FIG. 1 is a schematic diagram of an embodiment of the present invention having an internal probe with an attached transducer array. [Figure 10] FIG. 1 is a schematic diagram of an embodiment of the present invention in which different transducer arrays are mounted on separate modules. [Figure 11] FIG. 1 shows simulated results of the CS PAM method used to image two acoustic sources. [Figure 12] FIG. 1 shows the results of the CS PAM method used in vivo in a pig liver to image a single acoustic source. DETAILED DESCRIPTION OF THE INVENTION

[0018] We first describe the basic methodology of compressed sensing. Without loss of generality, consider the case of a single far-field acoustic source propagating from a location denoted as "x" towards a linear array of i={1,2,...,L} receiving elements or sensors, as shown in Figure 1. Consider the signal recorded at sensor i as:

[0019] ξ i (t)=w·s(t+Δ i (θ s )-R / c) (1) where w is the sensor coefficient, i.e., the conversion factor from the received acoustic signal detected as a variable pressure to the output electrical signal as a variable voltage, s is the signal generated by the acoustic source, R is the range, i.e., the distance from the array origin to the acoustic source, c is the speed of sound between the acoustic source and the sensor, and θ s is the angle of arrival for the acoustic wave from the acoustic source to the array origin (defined as the location of sensor 1 in Figure 1) measured relative to a line (or plane) through the sensor array, and Δ i (θ s ) is the relative time delay at sensor i. The relative time delay is parameterized by the angle of arrival as follows:

[0020]

number

[0021] This means that in the far field, θ s is assumed to be constant across all sensors. Each sensor signal is measured over N sampling periods. t for samples with sampling frequency F s =1 / T s When digitized by , the signal at each array element can be expressed by the following equation:

[0022]

number

[0023] Next, we introduce the angular space grid, B={θ1θ2… θ N} (where N defines the size and resolution of the grid) and a sparsity pattern vector b that selects the actual acoustic source angle from the angle space grid. Thus, the sparsity pattern vector b only has non-zero entries at the indices where the grid angle to be selected is located. This can be related to the received data as follows:

[0024] ξ i =Ψ i b (4) where Ψ i is a sparse basis that defines the sample values ​​for all sensor elements at each sample time for all possible acoustic source angles. i The jth column of

[0025]

number

[0026] where the subscript N t represents an arbitrary number of samples, and j is the index of the possible acoustic source angles. The problem in this context is to recover the sparse vector b (and therefore the true source angle) from the compressed measurements. Rather than measuring every signal at every sample time, we measure M samples or "compressed measurements" per sensor to reduce the data by taking linear projections on M random vectors, producing the following model:

[0027] β i =Φ i ξ i =Φ i Ψ i b, i=1,…,L (6) Here, (M×N t ) Random projection matrix Φ i teeth,

[0028]

number

[0029] and each line is defined by

[0030]

number

[0031] (m=1,...,M) is a random (e.g., independent and identically distributed Bernoulli or Gaussian) row vector that defines, for each array element, a random group of sample times selected from all sample times within the sample period. Φ i The particular choice of Ψ i The aim is to ensure small cross-correlation with

[0032] The goal is to solve from (L x M) measurements to find the angle of arrival, where L is the number of sensors and M is the number of samples taken of the signal from each sensor with respect to b. The l1 minimization problem to solve is written as follows:

[0033]

number

[0034] Here, Φ=diag{Φ1,…,Φ L}, Ψ=[Ψ1,…,Ψ L ], and β = [β1,…,β L ]. However, if the acoustic source signal s is unknown, the l1 minimization cannot be similarly constructed. In this case, assuming that the true acoustic source signal s is recorded in the sampled sensor data at each sensor, one of the sensors can be selected as a reference sensor, which should typically be sampled at the Nyquist rate. Compressed Sensing PAM To adapt the concept of compressed sensing far-field beamforming to near-field imaging, PAM, a voxel space grid is defined instead of an angular space grid. First, the coordinate r s =[x s ,z s Assume there is a single known bubble at a particular defined voxel (defined in terms of 2D imaging with a linear sensor array along the x-axis, on a horizontal or xz grid) with ]. This means that

[0035]

number

[0036] yields a signal model for sensor i in where p is the acoustic pressure at the source s and w i is the piezoelectric efficiency of the sensor element, and r i =[x i ,z i ] is the sensor position, and 1 / (4π·|r i -r s |) takes into account spherical diffusion. α i =w i / (4π·|r i -r s |), and then sampling this received data as previously shown in equation (3), we obtain the following for each sensor:

[0037]

number

[0038] Here, the defined voxel grid [x min ,z min ]×[x max ,z max Assume a sparse grid B of voxel locations (voxels and pixels are synonymous depending on whether the system is imaging in 2D or 3D, respectively) determined by:

[0039] B={ρ1ρ2… ρ N} (11) where N ultimately defines the resolution of the grid. If the true acoustic source is contained in the grid, i.e., the acoustic source voxel position r s is ρ j If it appears in the grid as one of the candidate locations, then its location should be recoverable using the same method. Similar to the far-field method described above, let b be the sparsity pattern vector, which is related to the received data by the following equation:

[0040] ξ i =Ψ i b (12) Again, Ψ i is a sparse basis for each sensor i=1,...,L, whose j-th column is a vector defined by the following equation:

[0041]

number

[0042] where j identifies a particular voxel location in the grid corresponding to the index j∈{1,…,N}, and for each sensor i=1,…,L, α i,j =w i / (4π·|r i -ρ j |).

[0043] The problem with CS-PAM is s The goal is to recover the σ, and therefore proceed by measuring random samples as linear projections from each sensor signal. This provides an equivalent form to equations (6)-(8), but now

[0044]

number

[0045] provides a sparse vector of likely voxel locations. Projections are described in more detail below. One thing to note is that the formulation assumes a known bubble and does not assume multiple unknown bubbles. In the presence of multiple bubbles, the formulation holds, provided there are limitations on the mutual coherence between the individual bubbles. As with the beamforming example, the bubble signal is unknown, so a reference sensor sampled and recorded at the Nyquist rate must be used, as described next.

[0046] In practice, the bubble signal as a function of time, p(t), is unknown, and therefore the basis cannot be defined as in (13). In this case, we follow the approach in A.C. Gurbuz, V. Cevher, and J.H. McClellan, "Bearing Estimation via Spatial Sparsity using Compressive Sensing," IEEE Transactions on Aerospace and Electronic Systems, vol. 48, No. 2, pp. 1358-1369, 2012, doi:10.1109 / TAES.2012.6178067, and use the subscript i RS Select a reference sensor (RS) with ∈{1,...,N};

[0047]

number

[0048] Let denote the signal measured at RS. The basis is formed by scaling and time-shifting the full-rate samples of RS as follows:

[0049]

number

[0050] The rationale for this approach is that when the sensor data does not contain noise from a single bubble, as in (9), and when the grid in (11) contains acoustic source locations (i.e., ρ j =r s ), (12) using (14) yields a sparse solution vector b with a single non-zero element at index j for all sensors i=1,...,L, except for RS. RS For RS with j=1,...,N,

[0051]

number

[0052] All columns of are identical and equal to samples of the RS signal, and therefore RS is excluded when solving for the sparse solution vector. CS-PAM in the frequency domain The time delays of the RS signals required to form the dictionary basis vectors in (14) can be efficiently obtained in the frequency domain as described below. The system can also use a two-level compressive sampling method across frequency bins and sensor data, which corresponds to a linear projection (i.e., compressive sampling) across space and time. The CS-PAM image is then obtained directly from the compressed data.

[0053] Sampled sensor data ξ i (t0),ξ i (t0+T s ),…,ξ i (t0+(N t -1)T s The discrete Fourier transform (DFT) of

[0054]

number

[0055] It is defined as follows. This is the frequency bin index k=0,…,N f -1 and the sensor index i=1,...,L, where the number of frequency bins N f ≧N t is usually chosen as an integer power of 2 so that the DFT can be efficiently computed with a Fast Fourier Transform (FFT) algorithm. The DFT frequency bin values ​​are s / N f Since the sensor data is a real value,

[0056]

number

[0057] Only k is unique, k=1,…,N f Regarding 2-1

[0058]

number

[0059] where: * denotes the complex conjugate. The subscripts k=0 and k=N f DFT bins with sampling rate F s (The components at DC and Nyquist frequencies in terms of samples / second.) In general, additional frequency bins may be excluded to achieve a filtering effect, for example, to take into account the frequency response of the sensor, band-stop filtering at frequency intervals around harmonics of the HIFU irradiation pulse, or other purposes. After this filtering, N b ≦N f Let / 2-1 denote the number of frequency bins retained for each sensor. Also, N b corresponds to the frequency bins to be retained {1,…,N f / 2-1}, and the elements of this set are

[0060]

number

[0061] where the frequencies to be retained in Hz are k=1,…,N b Regarding

[0062]

number

[0063] Finally, we denote the complex-valued frequency-domain samples of the filtered sensor data by the vector Ξ i (i=1,…,L), and each vector has dimension N b , and the elements from (15)

[0064]

number

[0065] (k=1,…,N b ) These vectors are used for subsequent processing to generate the CS-PAM image. Next, we form the frequency domain counterparts of the dictionary basis vectors in (14). We use the same notation Ψ for the frequency domain dictionary basis vectors. i where the matrix for each sensor i=1,...,L is of dimension N b The frequency domain RS data is a complex value with elements

[0066]

number

[0067] (k=1,…,N b ) vector

[0068]

number

[0069] where the elements in the frequency domain dictionary matrix are given by:

[0070]

number

[0071] where f k are the frequency components retained after filtering in Hz, for i∈{1,…,L} and i≠i RS , k=1,…,N b , j=1,...,N. Then, similar to (4), Ψ i The columns of provide a basis in which the frequency domain sensor data is described by a sparse vector b with non-zero elements, with indices corresponding to the bubble locations.

[0072] Ξ i =Ψ i b, i∈{1,…,L} and i≠i RS (17) The goal of CS-PAM is to recover b from the compressed data. As explained in the previous section, RS (16) in

[0073]

number

[0074] Since RS is a matrix with identical columns equal to ,RS is excluded in (17). Compressed sampling (CS) is done in two steps. First, the data for each sensor is b complex frequency samples to M F ≦N b The data is then reduced to L to M real frequency samples across the sensors. L ≦L values, total M=M F M LCompression across sensors is achieved by multiple linear combinations of the compressed data from each sensor. Each linear combination has random coefficients, which generate virtual sensor data corresponding to the (compressed) data from a small number of virtual sensors. The linear combinations are the same for each possible bubble position. This results in a compressed dictionary basis A' for each sensor. i The N t The overall compression ratio for time domain samples is:

[0075]

number

[0076] It has been found that CS-PAM images can be formed using only a few percent of the original data. The compression across frequency bins has complex-valued independent and identically distributed (iid) elements (M F ×N b ) Frequency compression matrix Φ' i and the compressed sensor data β' i and the compressed dictionary matrix A' i are generated, which form the compressed dictionary.

[0077] β' i =Re{Φ' i Ξ i}(M F ×1) A' i =Re{Φ' i Ψ i}(M F ×N), i∈{1,…,L}, i≠i RS (19) where Re{·} is the real part. The real part arises because the range (F s / 2,F s) must be compressed with the corresponding complex conjugate weights, and including these terms in the linear combination produces the real part in (19).

[0078] The compression across sensors is determined by the real-valued i.i.d. element φ m,i (M L ×(L−1)) sensor compression matrix Φ, yielding:

[0079]

number

[0080] Then, similar to (8), the l1 minimization problem to be solved is:

[0081]

number

[0082] where:

[0083]

number

[0084] is. A sparse solution vector that defines which of the grid of source locations contains the source

[0085]

number

[0086] is the voxel grid B={ρ1ρ2… ρ N}

[0087]

number

[0088] is used to generate a CS-PAM image by plotting the values ​​of . To allow for noise in the sensor data and other model mismatches, the following problem is solved instead:

[0089]

number

[0090] where ε is the error value. We now describe further examples of random projection matrices and therefore the methods used by the processing means to perform two-stage compression with a single matrix.

[0091] All frequency bins k=1,…,N f / 2-1 Elements

[0092]

number

[0093] vector Ξ' containing frequency domain data from all sensors with i (i=1,…,L). Then, Φ'' i Let us denote the complex-valued random matrix Φ' in (19). i and the same dimension (M F ×(N f / 2-1)), and the columns of Φ' are i column

[0094]

number

[0095] and Φ'' i The other columns of are zero, which filters the selected frequency bins. Then, the augmented dictionary basis vector matrix Ψ' is obtained as in (16). i is defined similarly, but has more rows and therefore dimensions ((N f / 2-1)×N), where f in (16)k teeth

[0096]

number

[0097] where k=1,…,N f / 2-1. In this case, (19) becomes β' i =Re{Φ'' i Ξ' i}(M F ×1) A' i =Re{Φ'' i Ψ' i}(M F ×N) (23) where the compression and filtering across frequency bins is done by the projection matrix Φ'' i Then, Φ(its element φ m,i is found in (20)) L × L), and an additional column of zeros is added to the column i corresponding to the index of the reference sensor. RS The projection matrix Φ' compresses across sensors, and the elements of Φ' are m,i It is expressed as:

[0098] The matrix A and vector β in (21), which correspond to the compressed dictionary and compressed sensor data, respectively, can be converted into a single projection matrix Φ as follows: T First,

[0099]

number

[0100] Define Next, β=Re{Φ T Ξ'} and A=Re{Φ T Ψ'} (25) Define Φ TThe block in

[0101]

number

[0102] is defined by Φ T The structure of Φ'' i shows that for each sensor i=1,...,L, we compress across frequency bins, where each row of the Φ′ matrix contains the coefficients of a linear combination that compresses across sensors. The formula in (26) performs two stages of compression using a single projection matrix.

[0103] To form the dictionary matrix in (16) and the matrix A from (19)-(21), RS must be sampled at the full Nyquist rate. However, the β in (20) m The compressed samples from the other sensors at are the continuous-time sensor signals ξ observed over the time interval t ∈ [t0,t0+T]. i (t), where the observation time or sample period is T=N t T s The matrix Φ' in (19) i is the dimension (M F ×N b ), and so these can be augmented with zero columns to give dimensions (M F ×(N f / 2-1)) with Φ'' i Φ' is formed. i The columns of Φ'' i column

[0104]

number

[0105] and Φ'' i The other columns of are zero. Then, i to form the inverse FFT of each row of Φ′′, a matrix containing real-valued time samples. iis formed, and then each row is interpolated to form a continuous time function Φ’’ m,i (t)(m = 1, …, M F ). Then, the compressed samples in (19) are obtained as follows.

[0106]

Number

[0107] Full Nyquist rate sampling of sensor signals other than RS is not necessary. Formulation of Multiple Measurement Vectors (MMV) The CS-PAM method described above can be called a single measurement vector (SMV) formulation. This is because the compressed measurements from all sensors and frequency bins are stacked into the vector β defined after (21), and correspondingly, the compressed dictionary is also stacked into the matrix A defined after (21). The linear model β = Ab in (21) has a total of N unknowns in the vector b (however, the sparse solution to (22)

[0108]

Number

[0109] has far fewer non-zero entries than N), and β has M F ·M L measurements, so β = Ab is a set of M F ·M L equations. N is the number of voxels in the grid defined in (11), and M F , M L should be recalled to characterize the compression over frequency bins and sensors respectively. The sparse recovery algorithm for solving (22) typically assumes that the model β = Ab is underdetermined (fewer equations than unknowns), which requires M F ·M L < N. Excessive compression may be required to meet this condition, and furthermore, the model Ξ in (17)i =Ψ i b is the same b vector, but the data Ξ i Let Ψ be the dictionary basis for that sensor. i This suggests describing the CS-PAM in terms of σ, which can be overly limiting. The Multiple Measurement Vector (MMV) formulation of CS-PAM addresses both of these issues and is described below. The MMV formulation was not considered in previous work by A.C. Gurbuz, V. Cevher, and J.H. McClellan, "Bearing Estimation via Spatial Sparsity using Compressive Sensing," IEEE Transactions on Aerospace and Electronic Systems, vol. 48, No. 2, pp. 1358-1369, 2012, doi:10.1109 / TAES.2012.6178067.

[0110] The MMV formulation assigns a different vector b i Modify (17) to include Ξ i =Ψ i b i , i∈{1,…,L} and i≠i RS (28) each b i The vector is sparse and has non-zero elements at the same positions, indicating the location of the voxels where energy is emitted. The standard MMV model in compressed sensing uses the same dictionary Ψ for each sensor, and therefore the model in (28) is called generalized MMV (GMMV) because each sensor has a different dictionary matrix defined by (16). Compression of the measurements and dictionary matrix across frequency bins is done similarly to (19), resulting in the following GMMV model for the compressed data:

[0111] β' i =A' i b i , i∈{1,…,L} and i≠i RS (29) MF When N, these equations for each sensor are underdetermined, which can be achieved more easily than in the case of the SMV model. In the GMMV formulation, no compression is performed across sensors, and thus the processing in (20) is omitted. This is because with a fixed or limited number of sensors, only a slight reduction in compression across sensors is possible. In contrast, for bubble imaging, most of the compression can be achieved across frequency bins. Here, the generalization of the l1 minimization problem in (22) is as follows.

[0112]

Number

[0113] Computationally efficient algorithms for solving (30) are described in (L. Li, X. Huang, and J. A. K. Suykens, "Signal recovery for jointly sparse vectors with different sensing matrices," Signal Processing, Volume 108, 2015, Pages 451-458, ISSN 0165-1684, https: / / doi.org / 10.1016 / j.sigpro.2014.10.010) and (R. Heckel and H. Bolcskei, "Joint sparsity with different measurement matrices," 2012 50th Annual Allerton Conference on Communication, Control, and Computing (Allerton), 2012, pp. 698-702, doi:10.1109 / Allerton.2012.6483286). The GMMV formulation can also be viewed as the block-sparse model discussed in (Y. C. Eldar, P. Kuppinger, and H. Bolcskei, "Block-sparse signals: Uncertainty relations and efficient recovery," IEEE Trans. Signal Process., vol. 58, No. 6, pp. 3042-3054, 2010). i When = b, the GMMV problem in (30) becomes the SMV problem in (22). Methods for determining voxel grid resolution and CS parameters The spacing (or resolution) of the voxel grid in (11) must be chosen appropriately. If the spacing is too large, the dictionary basis cannot adequately model the sensor measurements, and acoustic sources are lost in the CS-PAM image. If the spacing is too small, the dictionary basis has strong consistency, which limits the number of distinct acoustic sources that can be recovered in the CS-PAM image. In this section, we describe a method for selecting an appropriate resolution of the voxel grid in the lateral and axial directions.

[0114] Region of interest [x min ,z min ]×[x max ,z max ] at voxel position r s Consider an acoustic source located at . Assume that p(t) is the time waveform emitted by the acoustic source. If the sensor element position is r i (i=1,…,L), the time delay waveform at the sensor is

[0115]

number

[0116] In this discussion, we omit the spherical diffusion coefficient and the scaling factor. If p(t) is sampled over the measurement interval, as mentioned in the paragraph after (15), we can define Ξ as a function of dimension N after removing certain frequency bins. b Then, the DFT vector Ξ at each sensor is i has phase-shifted elements as follows:

[0117]

number

[0118] Next, the acoustic source position r s For a set of voxel locations ρ densely sampled (with small grid spacing) over a region near , the dictionary basis vector Ψ for each sensor i i(ρ) is formed. The index of the reference sensor is i RS If , according to (16) with the scaling factor omitted, these dictionary vectors are

[0119]

number

[0120] It has. The following dot product then evaluates the match between the sensor data in (31) and the dictionary basis vectors for voxel locations ρ in the vicinity of the acoustic source location:

[0121]

number

[0122] The inner product (33) is the acoustic source position ρ=r s achieves a maximum value at

[0123]

number

[0124] and f(ρ) ≦ f(r s ) The function f(ρ) is a function of r with a null. s The main lobe near r s The voxel grid spacing in CS-PAM in the axial and lateral directions is r s If the adjacent voxels are f(ρ)≒f(r s ) / 2. This process involves scanning various source locations r within the region of interest to find a suitable grid spacing for the entire region. s, or the grid spacing can be non-uniform across the region if desired. The axial resolution is determined by the number of frequency bins (larger bandwidth suggests better axial resolution), and the lateral resolution is determined by the array aperture (larger aperture suggests better lateral resolution). (32) is the power spectrum of the acoustic source signal at the selected frequency bin |[Ξ] k | 2 Note that these values ​​can be set to 1 for simplicity, or more realistically for CS-PAM, can be calculated from a model for the cavitation bubble waveform induced by HIFU irradiation.

[0125] This method reduces the dimension of the compressed sampling matrix (M F and M L ) can be extended to include compression across frequency bins and sensors. In this case, (19) and (20) can be used to select Ξ in (31). i and Ψ in (32) i (ρ) is calculated using m and A m (ρ), and then the dot product in (33) is transformed as follows:

[0126] f(ρ)=β T A(ρ) (34) where β and A(ρ) are formed by stacking operations in the lines after (21). The goal is to find the smallest possible M to maximize compression. F and M L and f(ρ) in (34) is selected as r s The goal is also to ensure that the main lobe structure for voxel locations ρ close to ρ is preserved. If the main lobe structure is lost, the compressed data will not retain enough information to locate the acoustic source.

[0127] Referring to FIG. 2, an ultrasound system 200 according to one embodiment of the present invention may include a geometrically focused ultrasound transducer 201 with an array of ultrasound transducer elements 202 positioned in an aperture 203 at the center of the transducer 201. Each transducer element 202 may be operable to generate and detect ultrasound waves. Thus, it may be usable in an active mode, which generates and detects ultrasound waves to produce a reflected (e.g., B-mode) ultrasound image, or in a passive mode, which detects ultrasound waves only. The array may be a linear or convex linear array extending primarily in a direction referred to as the x-direction, as shown in FIG. 2. The direction perpendicular to the x-direction, along the axis of the transducer, is referred to as the z-direction. Thus, the imaging surface of the array is the xz plane. The direction perpendicular to both the x- and z-directions is referred to as the y-direction.

[0128] The control unit 204 is configured to control the generation of ultrasound signals by each transducer element 202 and to receive detection signals from each transducer element 202. The control unit 204 may include a preamplifier and filter block 206 and an oscillator 208. The control unit 204 is also configured to control the transducer 201, for example, to control the power and frequency of ultrasound waves generated by the transducer 201, using a signal from the oscillator 208 to control the frequency of the ultrasound waves. Although the control unit 204 is described functionally, it should be understood that it may be implemented as a single processor or as two or more separate processors performing different functions within the system, such as control and analysis functions. The control unit is connected to a display screen 210 that can display data derived from the detection signals in an appropriate format. In this case, the therapy transducer 201 has a focus at a focal region 214, where it generates ultrasound waves of the highest intensity.

[0129] The configuration of Figure 2 can be implemented using a variety of components and systems, but in one embodiment, an ultrasound data acquisition system may be used that enables simultaneous raw radio frequency (RF) data or in-phase / quadrature (I / Q) data (which can be demodulated back to RF) from multiple individual elements 202 across a wide ultrasound bandwidth (e.g., 1-15 MHz). When the array is used in passive mode, pulse transmission can be switched off so that the array operates in receive only. In some modes, one group of transducer elements 202 is used in active mode and another group is used in passive mode, thus allowing active and passive detection to be used simultaneously. To make the system clinically applicable, a therapeutic ultrasound transducer 201 with a central aperture 203 for a linear (rectilinear or convex) detector array 202 may be used.

[0130] 3 and 4, the control unit 204 may be configured to control the transducer elements 202 as a phased array. For example, the control unit may be configured to generate a transmit signal for each transducer element 202 to control the frequency and timing, i.e., relative phase, of the vibrations produced by each transducer element 202, and therefore the ultrasound waves. Typically, the vibration frequency of each transducer element is controlled to be the same, and the phase or timing of each element is varied to steer the ultrasound waves generated by the array as a whole. The transducer elements may be configured to vibrate in phase with each other, which generates ultrasound waves with straight, parallel wavefronts 220 that all travel in the same direction, as shown in FIG. 2. This is suitable for anatomical ultrasound (B-mode or other reflection) imaging. As shown in FIG. 4, if the phase of the vibrations of the elements 202 are shifted so that the vibrations at the outer edges of the array are in phase with each other and the delay increases toward the center of the array, this generates ultrasound waves with curved wavefronts 222 that converge at the focal region 214.

[0131] It should be appreciated that rather than having a separate transducer 201 for focused ultrasound transmission, an array of ultrasound elements 204 can be used to generate ultrasound waves for both anatomical imaging and therapeutic FUS by rapidly switching between two different phase configurations. Furthermore, because both the anatomical ultrasound imaging and passive acoustic mapping described above require detection of ultrasound waves reflected or generated by the tissue being imaged, this detection can, in some embodiments of the present invention, be performed by the same transducer elements 202 used to transmit ultrasound. However, this requires further time division, and it is often preferable to have separate transducer arrays for detecting or receiving ultrasound waves. Thus, in the following description referring to transmit and receive arrays, these are typically separate arrays, but may alternatively be the same array.

[0132] Referring to FIG. 5, in an actual implementation of the system, the transducer elements 202 are divided into two groups, with the transducer elements of each group arranged in an array: a transmit array 202a and a receive array 202b, both of which are attached to the probe 230. Each of these arrays 202a, 202b may be a linear array, as described above. One of the sensor elements is designated as a reference sensor 202c, which may be located, for example, at the center of the receive array 202b. Thus, the separate geometric FUS transducer 201 of FIG. 1 is omitted. In a variation of this configuration, the transducer elements are provided in a therapy group and an imaging group. The therapy group of transducer elements may then be configured to generate focused ultrasound, and the imaging group of transducer elements may be configured to both transmit and receive reflected imaging ultrasound and to receive passive cavitation mapping signals. An imaging group may be configured to transmit flat, planar or focused ultrasound waves, for example using subgroups of transducer elements that are controlled to provide different focal points at a common focal depth, to provide good quality images at that depth.

[0133] The probe may further comprise a position sensor 232 suitable for any of a number of known position and orientation sensing systems. For example, the position sensor 232 may be a sensor of an electromagnetic motion tracking system and may comprise one or more markers that allow the position of the probe to be tracked by a stereo optical camera and an infrared camera or a laser tracker.

[0134] The control unit 204 includes a main processor 204a and an ultrasound front-end processor 204b. The ultrasound front-end 204b includes a transmit side and a receive side. The transmit side includes generators configured to generate ultrasound signals to be transmitted by the transmit array 202a. For example, these may include a therapy signal generator 240 and an imaging signal generator 242. Each of these signal generators 240, 242 is configured to output a transmit signal that is input to the individual transducer elements of the transmit array 202a, thereby transmitting ultrasound signals in the required formats for therapeutic FUS and B-mode ultrasound imaging, respectively, as shown in FIGS. 3 and 4. The ultrasound front-end 204b further includes two amplifiers 244, 246 configured to amplify the transmit signals from the signal generators 240, 242, respectively, and a multiplexer 248 configured to receive both amplified transmit signals and transmit each amplified transmit signal to the appropriate transmit array 202a at the appropriate time in a time-division multiplexed manner. A timing controller 264, which forms part of the main processor 204a, provides trigger or synchronization signals to coordinate the timing of the signal generators and multiplexers, as described in more detail below.

[0135] The receive side of the ultrasound front end 204b includes a multiplexer 250 configured to receive all detection signals from the receive array 202b and separate the anatomical ultrasound imaging signals (e.g., B-mode) from the passive acoustic mapping signals into the anatomical ultrasound imaging channel and the PAM channel. This separation is based on timing, since the received anatomical imaging signals are reflections of the transmitted anatomical imaging signals, and the received PAM signals are generated by cavitation caused by the transmitted FUS signals. The anatomical imaging signals are time-gain compensated by a TGC module 252, digitized by an ADC 254, and filtered by a digital filter 256. The PAM signals are filtered by an analog filter 258 to separate the high-frequency broadband signal used for PAM and amplified by a low-noise amplifier 260.

[0136] The receive side of the ultrasonic front end 204a further comprises another filter 270 configured to filter the output signal from the reference sensor 202c to isolate the high frequency wideband signal used for PAM, and a low noise amplifier 272 configured to amplify the output from the filter.

[0137] The filtered and amplified outputs from each of the main array of sensor elements and the reference sensor are all input to an ADC 290, which is configured to sample each signal at a series of sample points in the time domain at the Nyquist rate and output a digitized signal. An FFT module 292 receives the digitized signals and transforms them into the frequency domain for each sample period. A random projection module 261 receives the frequency components from each sensor signal and performs a selection of those frequency components for sensor data compression.

[0138] The FFT data from the reference sensor signals is input to a sparse dictionary and parameters unit 278, and the FFT data from the main sensor array is input to a random projection unit 294, which performs frequency component selection from each sensor signal similar to the random projection unit 261 in FIG. 5 except that it operates on digital signals.

[0139] This can be done in hardware using a filter component and an integrator that integrates each sensor data over each sample period. The integrated value is then sampled and digitized by the ADC 262 once per sample period.

[0140] A high-speed ADC 274 is configured to sample the filtered reference sensor signal at each sample frequency, allowing construction of a frequency domain basis. The main processor 204a includes a timing controller 264 configured to provide timing inputs to the ultrasound signal generators 240, 242, the multiplexers 248, 250, and the ADCs 254, 262, as described above. The timing inputs to the ultrasound signal generators 240, 242 are configured to trigger alternating generation of therapeutic FUS signals and anatomical ultrasound imaging signals over respective series of short periods. The multiplexer 248 uses the timing signals received from the timing controller 264 to multiplex the two types of ultrasound generation signals onto a control input to the transmitter array 202a, causing the transmitter array 202a to generate alternating pulses of therapeutic FUS and focused or unfocused pulsed imaging ultrasound. The receive multiplexer 250 is controlled by the timing controller 264 to switch signals from the receive array 202b between the B-mode channel and the PAM channel based on the transmission times of the two different types of ultrasound and known or expected delays between transmission and reception of the associated ultrasound signals. The timing controller 264 also provides timing inputs to the ADCs 254, 262 to control the sampling rate of the analog signals by the ADCs.

[0141] The main processor 204a further comprises an apodization unit 266 and an image reconstruction unit 268, which are configured to receive filtered digital signals from the B-mode channel of the ultrasound front end 204b and generate from those signals anatomical ultrasound images in the form of a 2D B-mode ultrasound imaging stream, which may include, for example, a series of time-stamped 2D image frames.

[0142] The main processor 204a further includes a sparse dictionary and parameter unit 278 and a PAM imaging unit 280. The sparse dictionary and parameter unit 278 is configured to receive the digitized frequency-domain reference signal from the FFT unit 292 and generate therefrom a sparse dictionary basis (compressed dictionary) defined in equation (19). The PAM imaging unit 280 is configured to receive the compressed data signal from the random projection unit 294 and the sparse dictionary from the sparse dictionary and parameter unit 278 and generate therefrom a PAM image, for example in the form of a 2D PAM imaging stream having one image frame per sample period.

[0143] The main processor 204a further comprises a local memory 274 and a data processor 276. The local memory 274 is configured to receive and store image streams from the anatomical and PAM imaging channels. It is also configured to receive and store 3D anatomical image data from a previous scan of the patient, for example, in the form of a set of 2D image slices obtained by a CT or MRI scan (which may be in DICOM format, for example). The data processor 276 is configured to process the data stored in the local memory to generate a composite image or image stream in which the PAM image of the FUS-induced cavitation is superimposed on a 3D anatomical image of the patient, which may be a pre-scanned 3D image or a 3D ultrasound image generated from compounded 2D B-mode ultrasound image slices. The processing method by which the composite image can be obtained is not relevant to the present invention and will not be described.

[0144] Referring to Figure 6, in a variation of the embodiment of Figure 5, the sensor signal is sampled in the frequency domain. Components in this embodiment corresponding to those in Figure 5 are given the same reference numerals and will not be described in detail. A random projection unit 261 is configured to receive the amplified sensor output signals from each sensor element 202b and filter and compress them as discussed above with reference to equation (27), and the integral in (27) is sampled by an ADC 262 at the image frame rate. These compressed samples are then input to a CS-PAM unit 280. A separate ADC 274 is configured to sample the RS signal data at the full Nyquist rate, an FFT unit 275 transforms the data to the frequency domain for all frequency bins, and the frequency-domain reference data is input to a sparse dictionary and parameter unit 278.

[0145] In a further embodiment, the sensor signals, including the reference signal, are sampled in the time domain, there is no FFT transformation to the frequency domain, the basis is constructed in the time domain as described above in equation (14), and the voxel locations are determined by minimization accordingly. A system for this method has components corresponding to those in Figure 6, except that the FFT unit 275 is omitted, the reduced sampling of the sensor signals is performed in the time domain, and the compressed dictionary is defined in the time domain.

[0146] 7-9, the systems of FIGS. 2-5 or 6 can be incorporated into a variety of different probe configurations. For example, referring to FIG. 7, the probe 230 may be connected to an electronically or manually controlled motion control system 600, such as a robotic arm or multi-axis motion system, configured to support the probe and control its movement in six degrees of freedom (three along orthogonal translational motion directions and three about orthogonal rotational axes) in response to inputs from an electronic control system or manual inputs. In this case, the position and orientation of the probe 230 may be determined from a position sensor 632 attached to the probe, as described above, or from the operation of the motion control system.

[0147] 8, the probe may be a handheld external probe 730 having a contact surface 700 configured to be placed in contact with the patient's skin, with the transmit array 702a and receive array 702b mounted within or adjacent to the contact surface 700. In this case, a position sensor 732 is mounted within or fixed to the body of the probe 730 and can record any manual movement of the probe.

[0148] 9, the probe may be a handheld intracorporeal or intraoperative probe 830 having a contact surface 800 configured to be placed in contact with a patient within the body, with a transmit array 802a and a receive array 802b mounted within or adjacent to the contact surface 800. Again, a position sensor 832 is mounted on the body of the probe 830 and can record any manual movement of the probe in a similar manner as described above.

[0149] Referring to FIG. 10 , rather than all being mounted on a common probe, the system may include separate probes, each supporting a different group of transducer elements. For example, the system may include an anatomical imaging probe 900 supporting an array 902 of transducer elements. The array 902 may be a linear or convex array and may be configured to transmit reflected imaging ultrasound and receive reflected ultrasound, for example, as shown in FIG. 2 , using the same group of transducer elements in the array 902, different groups of transducer elements in the array 902, or separate arrays of transducer elements. The system may include an FUS probe 904 having an array 906 of transducer elements configured to generate therapeutic FUS, as shown in FIG. 3 . The system may further include a PAM probe 908 having an array 910 of transducer elements configured for reception. Each probe 900, 904, 908 may further include position sensors 912, 914, 916 so that the position of the probe, and therefore the transducer array mounted on the probe, can be monitored. This system can operate similarly to the system described above, except that the calibration step 300 used to spatially align the PAM image with the anatomical ultrasound image is replaced by a transformation that changes over time and is determined for each image frame of the PAM image stream based on the relative positions of the two probes 900, 908.

[0150] In a further embodiment, the system includes a FUS ultrasound transmitter mounted on a probe or robotic arm and a movable probe having PAM imaging receive transducer elements and both the transmit and receive elements of the anatomical ultrasound imaging system, either as separate arrays or as a single common array, or as a transmit array for anatomical ultrasound imaging and a receive array for both PAM imaging and anatomical ultrasound imaging.

[0151] Simulations of the system in Figure 5 were performed using a variety of different PAM methods, specifically a curvilinear array with L = 128 elements and a voxel grid with N = 2626 locations (101 lateral × 26 axial, 40 mm × 40 mm area). For CS-PAM, a compression ratio of approximately 2% was used. CS-PAM uses only 2% of the processed data compared to conventional PAM. Other methods used in the simulations were standard beamforming (SB), standard Capon beamforming (SCB), and robust Capon beamforming (RCB).

[0152] CS-PAM:M F =300, M L For =32 simulations, the optimization is solved using perturbed orthogonal matching pursuit (POMP). For all simulations, the bubble signal is generated with a random pulse train bubble model, SNR 40 dB, and the resolution is evaluated using two bubble sources, and the performance is evaluated using bubble clouds, on and off-grid (on-grid shown in Figure 11), from which it can be seen that the simulation using CS PAM has significantly less artifacts than the other methods.

[0153] In vivo treatment of a pig liver using SonoTran particles was also performed, and the results are shown for a single region of cavitation in Figure 12. Again, CS PAM has significantly smaller artifacts than the other methods.

Claims

1. 1. A passive compressional wave imaging system for locating cavitation bubbles, comprising: a plurality of pressure sensor elements arranged in an array, each configured to generate an output signal; and processing means, the processing means comprising: defining a sample period and a sample space within which the output signal can be sampled at each of a plurality of sample points; Define a grid of candidate bubble locations; defining, for each of said output signals, a random projection that identifies a different group of said sample points; sampling each of the output signals at the group of sample points defined by the random projection to generate sample data; defining a dictionary defining sample values ​​to be taken for the output signal of each of the respective groups of the plurality of pressure sensor elements at each of the sample points over the sample period for a bubble at each of the candidate bubble locations; defining a vector having a plurality of elements, each of the plurality of elements having a value, the values ​​of the plurality of elements specifying at which of the candidate bubble locations a bubble is located; performing a minimization operation to derive the values ​​of the elements of the vector and the corresponding identified bubble locations from the sample data, the dictionary, and the random projections; The system is configured as follows.

2. The system of claim 1 , wherein the sample space is time.

3. 3. The system of claim 2, wherein the processing means is further configured to define a sample rate that defines the sample points as a number of sample times within the sample period at which samples can be taken.

4. The system of claim 1 , wherein the sample space is frequency.

5. 5. The system of claim 4, wherein the processing means is further configured to define the sample points as a plurality of sample frequencies at which samples can be taken.

6. 6. A system according to claim 4 or 5, wherein the processing means is arranged to separate an analogue component of the output signal at each of the sample points and to sample each of the analogue components over the sample period.

7. 3. The system of claim 2, wherein the processing means is configured to sample each of the output signals and determine frequency content samples for each of the output signals from time domain samples.

8. 8. The system of claim 1, wherein the processing means is configured to form a linear combination of the sample data from the plurality of pressure sensor elements, thereby forming virtual sensor data corresponding to a plurality of virtual sensors.

9. The system of claim 1 , wherein the random projection is defined by at least one matrix.

10. 10. The system of claim 1, wherein the processing means is configured to identify one of the plurality of pressure sensor elements as a reference sensor element, sample the output signal from the reference sensor element at each of the sample points to generate reference data, define a delay time for each combination of one of the candidate bubble locations and one of the plurality of pressure sensor elements, and determine the dictionary from the reference data and the delay time.

11. 1. A method for passive compression wave imaging of cavitation bubbles, comprising: defining a sample period and a number of sample points across a sample space; defining a grid of candidate bubble locations; defining a random projection for each of a plurality of pressure sensor element output signals, the random projection identifying a respective different group of said sample points; sampling each of the output signals at the group of sample points defined by the random projection to generate sample data; defining a dictionary defining sample values ​​to be obtained for the output signal of each of the plurality of pressure sensor elements at each of the respective groups of sample points over the sample period for a bubble at each of the candidate bubble locations; defining a vector having a plurality of elements, each of the plurality of elements having a value, the values ​​of the plurality of elements specifying at which of the candidate bubble locations a bubble is located; performing a minimization operation to derive the values ​​of the elements of the vector and the corresponding identified bubble locations from the sample data, the dictionary, and the random projection; A method comprising:

12. The method of claim 11 , wherein the sample space is time.

13. The method of claim 12 further comprising defining a sample rate that defines the sample points as a number of sample times within the sample period at which samples can be taken.

14. The method of claim 11 , wherein the sample space is frequency.

15. The method of claim 14 further comprising defining the sample points as a plurality of sample frequencies at which samples can be taken.

16. 16. The method of claim 14 or 15, further comprising isolating an analog component of the output signal at each of the sample points and sampling each of the analog components over the sample period.

17. The method of claim 12 , further comprising sampling each of the output signals and determining frequency component samples for each of the output signals from the time domain samples.

18. 18. The method of claim 11, further comprising forming a linear combination of the sample data from the plurality of pressure sensor elements, thereby forming virtual sensor data corresponding to a plurality of virtual sensors.

19. 19. The method of any one of claims 11 to 18, wherein the random projection is defined by at least one matrix.

20. 20. The method of claim 11, further comprising: identifying one of the plurality of pressure sensor elements as a reference sensor element; sampling the output signal from the reference sensor element at each of the sample points to generate reference data; defining a delay time for each combination of one of the candidate bubble locations and one of the plurality of pressure sensor elements; and determining the dictionary from the reference data and the delay time.

Images