CN104887266A

Movatterモバイル変換

Info

Publication number: CN104887266A
Application number: CN201510290576.0A
Authority: CN
Inventors: 万明习; 路舒宽
Original assignee: Xian Jiaotong University
Current assignee: Xian Jiaotong University
Priority date: 2015-05-29
Filing date: 2015-05-29
Publication date: 2015-09-09
Anticipated expiration: 2035-05-29
Also published as: CN104887266B

Abstract

The invention provides a method for small-area three-dimensional passive cavitation imaging and three-dimensional composite imaging based on an area array, wherein a two-dimensional passive cavitation imaging algorithm is expanded to a three-dimensional algorithm and a multi-stage continuous self-adaption trapped wave system is designed to research a mechanism of inertial oscillation of a cavitation micro bubble in a sound field, and a dynamic aperture receiving technique is used for further processing of three-dimensional data. Delayed array element receiving of the area array is then disposed to adjust an angle of a passive receiving plane, three-dimensional conversion is carried out on the three-dimensional data and then three-dimensional recombination is carried out on the data in order to realize the three-dimensional passive cavitation imaging. By the method for the small-area three-dimensional passive cavitation imaging and the three-dimensional composite imaging based on the area array provided by the invention, three-dimensional passive cavitation multi-angle compound imaging with high resolution, high signal-to-noise ratio and high efficiency can be obtained, and a selective monitoring manner can be provided for clinical operations by combing the passive cavitation imaging with other techniques.

Description

Translated fromChinese

基于面阵的小区域三维被动空化成像及三维复合成像方法Small area 3D passive cavitation imaging and 3D composite imaging method based on area array

技术领域technical field

本发明涉及超声成像技术领域，特别涉及一种适用于高分辨率、高信噪比和高效率的三维被动空化成像方法。The invention relates to the technical field of ultrasonic imaging, in particular to a three-dimensional passive cavitation imaging method suitable for high resolution, high signal-to-noise ratio and high efficiency.

背景技术Background technique

高强度聚焦超声(High Intensity Focused Ultrasound,HIFU)换能器在聚焦声场作用下会在焦域附近产生空化效应，空化效应是指液体中的不稳定结构被破坏产生空腔(即微米级尺寸的气泡)并在外界驱动作用下产生的膨胀、收缩、塌陷的动力学过程，按其物理机制又可分为惯性空化和非惯性空化。惯性空化是指在负压相的作用下微泡急剧膨胀，在正压相的作用下则急剧收缩直至塌陷，这一过程中微泡在很小的区域内集聚了很大的能量并突然释放，因此会产生局部的高温高压、化学自由基、高速微射流和冲击波等物理和化学现象，可以对正常的细胞结构和生物体内酶的活性造成很大的破坏，同时可以对肿瘤细胞和癌变组织进行有效的杀伤从而达到用HIFU治疗的作用。随着超声造影剂研制的进展，空化效应与超声造影剂的协同作用，在基因治疗、药物传输、溶栓治栓、炎症和肿瘤的靶向显影与治疗等方面展现了巨大的应用。The High Intensity Focused Ultrasound (HIFU) transducer will produce cavitation effect near the focal region under the action of the focused sound field. The dynamic process of expansion, contraction, and collapse under the external drive can be divided into inertial cavitation and non-inertial cavitation according to its physical mechanism. Inertial cavitation refers to the rapid expansion of microbubbles under the action of negative pressure phase, and sharp contraction and collapse under the action of positive pressure phase. In this process, microbubbles accumulate a lot of energy in a small area and suddenly collapse. Therefore, local physical and chemical phenomena such as high temperature and high pressure, chemical free radicals, high-speed micro jets and shock waves will be generated, which can cause great damage to normal cell structure and enzyme activity in organisms, and can also damage tumor cells and cancer The tissue is effectively killed to achieve the effect of HIFU treatment. With the progress in the development of ultrasound contrast agents, the synergistic effect of cavitation and ultrasound contrast agents has shown great applications in gene therapy, drug delivery, thrombolysis and thrombolysis, and targeted development and treatment of inflammation and tumors.

基于空化微泡的声学成像技术大致有被动成像技术和主动成像技术。其中被动成像技术是利用一个换能器被动接收由于空化微泡所导致的声发射和声散射。HIFU治疗中使用核磁共振技术可实时准确监控组织温度，可防止目标区域温度过高造成不可逆的损伤，但是声学成像技术对于HIFU作用中空化瞬态过程的监控更加具备优势；作为一种被动的方法，它比磁共振成像技术可以提供更高的帧频，而且在以后的临床应用中由于廉价更加容易被医疗机构所接受。现有的被动空化成像(Passive Cavitation Imaging,PCI)技术主要是在时间域应用时间曝光声学(Time Exposure Acoustic,TEA)和通过快速的傅氏变换将接收数据转换到频率域，然后应用TEA原理进行成像。此外，稀疏解卷积、最小方差无失真响应(Minimum Variance Distortionless Response,MVDR)以及鲁棒的Capon自适应波束合成等技术也应用到了其中。Acoustic imaging techniques based on cavitation microbubbles generally include passive imaging techniques and active imaging techniques. Among them, the passive imaging technology uses a transducer to passively receive acoustic emission and acoustic scattering caused by cavitation microbubbles. The use of nuclear magnetic resonance technology in HIFU treatment can accurately monitor tissue temperature in real time and prevent irreversible damage caused by excessive temperature in the target area. However, acoustic imaging technology has more advantages in monitoring the transient process of cavitation during HIFU treatment; as a passive method , it can provide a higher frame rate than magnetic resonance imaging technology, and it is easier to be accepted by medical institutions because of its low cost in future clinical applications. The existing Passive Cavitation Imaging (PCI) technology mainly applies Time Exposure Acoustic (TEA) in the time domain and converts the received data to the frequency domain through fast Fourier transform, and then applies the principle of TEA for imaging. In addition, techniques such as sparse deconvolution, Minimum Variance Distortionless Response (MVDR) and robust Capon adaptive beamforming are also applied.

然而这些二维的PCI算法只能得到一个平面方向的空化成像，在应用中只能从一个角度反应组织或者器官的变化，临床医生无法选择最优的成像平面来诊断。例如在监控HIFU治疗肿瘤的过程中，操作者必须有丰富的经验，而且二维的成像受操作者的主观因素影响很大。目前超声三维重建技术是对一组二维图像进行目标定位和空体素补差后进行三维重建，提供空间三维尺度上更丰富的信息量，是在二维成像基础上所做的后续处理；基于像素的三维重建算法(Pixel-BasedMethods,PBM)是在已经获取二维超声图像的基础上，将图像上像素的坐标转换到三维区域中的某个体素并对其中未赋值的体素插值，但是这种方法计算速度较慢；基于体素的三维重建算法(Voxel-Based Methods,VBM)首先定义一个三维标准晶格并对其中的每一个体素进行插值，插值方法有中值滤波算法、距离加权算法、线性插值算法、基于探头轨迹的算法等等。但是这两种方法都无法实现三维实时成像，而且当二维超声图像分辨率较低时，还有可能会产生比较大的图像伪影和误差。因此，一幅三维超声成像质量的差异主要还是由二维成像来决定，它只是从二维成像到三维成像的一种量变，无法达成真正意义上的“三维”，在临床上也无法大范围的使用。However, these two-dimensional PCI algorithms can only obtain cavitation imaging in one plane direction, and can only reflect changes in tissues or organs from one angle in application, and clinicians cannot choose the optimal imaging plane for diagnosis. For example, in the process of monitoring HIFU treatment of tumors, the operator must have rich experience, and the two-dimensional imaging is greatly affected by the subjective factors of the operator. At present, ultrasonic 3D reconstruction technology is to perform 3D reconstruction on a group of 2D images after target positioning and empty voxel compensation, which provides more abundant information on the spatial 3D scale, and is a follow-up processing based on 2D imaging; The pixel-based three-dimensional reconstruction algorithm (Pixel-Based Methods, PBM) is to convert the coordinates of the pixels on the image to a certain voxel in the three-dimensional area and interpolate the unassigned voxels on the basis of the acquired two-dimensional ultrasound image. The calculation speed of this method is slow; the voxel-based 3D reconstruction algorithm (Voxel-Based Methods, VBM) first defines a 3D standard lattice and interpolates each voxel in it. The interpolation methods include median filter algorithm, distance Weighting algorithms, linear interpolation algorithms, probe trajectory based algorithms, etc. However, these two methods cannot realize three-dimensional real-time imaging, and when the resolution of the two-dimensional ultrasound image is low, relatively large image artifacts and errors may be generated. Therefore, the difference in the quality of a three-dimensional ultrasound image is mainly determined by two-dimensional imaging, which is only a quantitative change from two-dimensional imaging to three-dimensional imaging. usage of.

在散乱数据插值方面应用较多的径向基函数(RBF)插值方法是样条插值的一种近似，可以得到平滑的三维重建数据。径向基函数可以保证无序三维数据相应的线性方程组矩阵可逆，因此不会出现一般FBM算法解方程不稳定的问题。不过用径向基函数实现三维重建要考虑到过大的计算量问题。传统单一的RBF插值方法对原始三维体数据不做任何处理，虽然这样准确度非常高，但是三维重建速度大大减慢，在临床上根本无法实现。所以考虑在插值之前对数据进行分割预处理，从而提高计算效率。另外，B样条插值对于均匀的矩形网格数据插值效果较好，由于二维面阵采集到的空化数据本身具备这个特性，因此可以将径向基函数插值和B样条插值方法结合起来用于被动空化成像的三维重建上。如此可以在保证计算稳定性的前提下缩短计算时间并使重建图像尽可能平滑。Radial basis function (RBF) interpolation method, which is widely used in interpolation of scattered data, is an approximation of spline interpolation, which can obtain smooth 3D reconstruction data. Radial basis functions can ensure that the matrix of linear equations corresponding to unordered three-dimensional data is invertible, so there will be no problem of unstable solution equations of general FBM algorithms. However, the implementation of 3D reconstruction with radial basis functions has to take into account the problem of excessive calculation. The traditional single RBF interpolation method does not do any processing on the original 3D volume data. Although the accuracy is very high, the 3D reconstruction speed is greatly slowed down, which is impossible to achieve clinically. So consider preprocessing the data before interpolation to improve computational efficiency. In addition, B-spline interpolation is better for uniform rectangular grid data interpolation. Since the cavitation data collected by the two-dimensional array has this characteristic, it is possible to combine radial basis function interpolation and B-spline interpolation methods. For 3D reconstruction of passive cavitation imaging. In this way, the calculation time can be shortened and the reconstructed image can be as smooth as possible under the premise of ensuring calculation stability.

鉴于以上内容，很有必要提出一种针对空化微泡的三维被动成像和三维复合成像方法。In view of the above, it is necessary to propose a 3D passive imaging and 3D composite imaging method for cavitation microbubbles.

发明内容Contents of the invention

针对二维PCI在成像角度、成像质量等方面的限制，本发明的目的在于提供一种基于面阵的小区域三维被动空化成像及三维复合成像方法。Aiming at the limitations of 2D PCI on imaging angle, imaging quality, etc., the purpose of the present invention is to provide a small area 3D passive cavitation imaging and 3D composite imaging method based on area array.

为实现上述目的，本发明采用了以下技术方案：To achieve the above object, the present invention adopts the following technical solutions:

步骤一：基于谱估计理论设计自适应滤波系统，对面阵采集到的三维数据进行自适应滤波，得到表征惯性空化的宽带噪声或得到表征稳态空化的次谐波，通过计算散射点与阵元之间的距离和对信号丢失和接收灵敏度的补偿系数，采用被动成像算法计算时间累计区间内的源强度三维数据，并应用动态孔径接收技术和幅度变迹技术对源强度三维数据进行处理，得到惯性空化或稳态空化的三维体数据；Step 1: Design an adaptive filtering system based on spectral estimation theory, and perform adaptive filtering on the 3D data collected by the area array to obtain broadband noise representing inertial cavitation or subharmonics representing steady-state cavitation. The distance between array elements and the compensation coefficient for signal loss and receiving sensitivity, the passive imaging algorithm is used to calculate the three-dimensional source intensity data in the time accumulation interval, and the dynamic aperture receiving technology and amplitude apodization technology are used to process the three-dimensional source intensity data , to obtain the three-dimensional volume data of inertial cavitation or steady-state cavitation;

步骤二：重复设置面阵阵元时间延时以调整成像角度，在每次调整成像角度后重新采集三维数据并按照步骤一处理，并将经过步骤一处理后的三维数据转换到参考坐标系下，对多次转换到参考坐标系下的三维数据处理结果取平均得到用于三维重建的预处理数据；Step 2: Repeatedly set the time delay of the area array element to adjust the imaging angle, re-acquire the 3D data after each adjustment of the imaging angle and process it according to step 1, and convert the 3D data processed in step 1 to the reference coordinate system , taking the average of the three-dimensional data processing results converted to the reference coordinate system multiple times to obtain preprocessed data for three-dimensional reconstruction;

步骤三：对所述预处理数据分割子区域并设置径向基函数和窗函数，然后结合径向基函数插值和B样条插值方法对所述预处理数据进行三维重建，得到三维空化成像。Step 3: Divide the preprocessed data into sub-regions and set radial basis function and window function, and then perform three-dimensional reconstruction on the preprocessed data by combining radial basis function interpolation and B-spline interpolation methods to obtain three-dimensional cavitation imaging .

步骤一中所述自适应滤波具体包括以下步骤：The adaptive filtering described in step 1 specifically includes the following steps:

(2.1)对面阵第一个阵元通道的信号进行快速傅里叶变换，得到空化信号的功率谱估计，设置两路参考信号，第一路参考信号为正弦信号x₁，第二路参考信号为x₁经过一个正交相移网络后的信号x₂，其中正弦信号的相位均为0，幅度C均为1，频率f为第一个谐波频率：(2.1) Perform fast Fourier transform on the signal of the first element channel of the area array to obtain the power spectrum estimation of the cavitation signal, and set two reference signals. The first reference signal is a sinusoidal signal x₁ , and the second reference signal The signal is the signal x₂ after x₁ passes through a quadrature phase shift network, in which the phase of the sinusoidal signal is 0, the amplitude C is 1, and the frequency f is the first harmonic frequency:

x₁＝C cos(2πft),x₂＝C sin(2πft) (5)x₁ ＝C cos(2πft), x₂ ＝C sin(2πft) (5)

(2.2)设置自适应陷波器中的权阶数为2，权系数的初始值分别为w₁＝0.1、w₂＝0.1，对步骤(2.1)中所述两路参考信号进行加权得到加权信号y＝w₁x₁(i)+w₂x₂(i)，计算误差信号第i个采样点对应的数值e(i)＝x(i)-y，其中x(i)为第一个阵元通道的信号，x₁(i)为第一路参考信号第i个采样点对应数值，x₂(i)为第二路参考信号第i个采样点对应数值；(2.2) Set the weight order number in the adaptive notch filter to be 2, the initial values of the weight coefficients are respectively w₁ =0.1, w₂ =0.1, and weight the two reference signals described in step (2.1) to obtain weighted Signal y=w₁ x₁ (i)+w₂ x₂ (i), calculate the value e(i)=x(i)-y corresponding to the i-th sampling point of the error signal, where x(i) is the first The signal of the array element channel, x₁ (i) is the value corresponding to the i-th sampling point of the first reference signal, and x₂ (i) is the corresponding value of the i-th sampling point of the second reference signal;

(2.3)设置自适应陷波的收敛系数μ，根据最陡下降法产生两个权系数的迭代公式：(2.3) Set the convergence coefficient μ of the adaptive notch, and generate the iterative formula of two weight coefficients according to the steepest descent method:

w₁＝w₁+2μe(i)x₁(i),w₂＝w₂+2μe(i)x₂(i) (6)w₁ =w₁ +2μe(i)x₁ (i),w₂ =w₂ +2μe(i)x₂ (i) (6)

(2.4)对于第一个阵元通道的所有采样点x(1),x(2),...,x(N)，重复步骤(2.2)和(2.3)，N为采样点数，得到误差信号e，误差信号即为滤除了第一个频率分量的信号；(2.4) For all sampling points x(1), x(2),...,x(N) of the first array element channel, repeat steps (2.2) and (2.3), N is the number of sampling points, and get the error Signal e, the error signal is the signal that has filtered out the first frequency component;

(2.5)级联下一个自适应陷波器，改变陷波频率，重复步骤(2.1)～(2.4)的过程，得到第一个阵元通道将高谐波和超高次谐波滤除掉之后的宽带噪声信号，单个自适应陷波器的传递函数为：(2.5) Cascade the next adaptive notch filter, change the notch frequency, repeat the process of steps (2.1) to (2.4), and obtain the first array element channel to filter out high harmonics and ultra-high harmonics After broadband noise signal, the transfer function of a single adaptive notch filter is:

$H h ((z z)) = = \frac{{z z}^{22} - - 22 z z cos cos wT wxya + + 11}{{z z}^{22} - - 22 z z cos cos wT wxya ((11 - - {μC μC}^{22})) + + ((11 - - 22 μ μ {C C}^{22}))} - - - - - - ((77))$

其中，T表示采样周期，z表示z变换的复参量，w表示角频率；Among them, T represents the sampling period, z represents the complex parameter of z transformation, and w represents the angular frequency;

(2.6)重复以上步骤，使面阵每个阵元采集到的信号经过多级连自适应陷波系统进行滤波处理，得到一组反映微泡惯性振动机制的三维数据。(2.6) Repeat the above steps, so that the signals collected by each element of the area array are filtered by a multi-stage adaptive notch system to obtain a set of three-dimensional data reflecting the inertial vibration mechanism of the microbubbles.

步骤一中所述采用被动成像算法计算时间累计区间内的源强度三维数据，具体包括以下步骤：In Step 1, the passive imaging algorithm is used to calculate the three-dimensional source intensity data within the time accumulation interval, which specifically includes the following steps:

(3.1)通过面阵被动接收到的射频压力信号p_i,j(x,y,z,t)经过滤波处理后得到的信号为：(3.1) The RF pressure signal p_i,j (x,y,z,t) received passively by the area array is filtered and processed for:

${\overset{^^}{p p}}_{i i,, j j} ((x x,, y the y,, z z,, t t)) = = {p p}_{i i,, j j} ((x x,, y the y,, z z,, t t)) * * F f ((t t)) - - - - - - ((88))$

其中F(t)为滤波函数，与p_i,j(x,y,z,t)在时域上做卷积运算；Among them, F(t) is a filter function, which is convolved with p_i,j (x,y,z,t) in the time domain;

(3.2)对做时间延时和补偿：(3.2) Yes Do time delay and compensation:

${F f}_{i i,, j j} ((x x,, y the y,, z z,, τ τ)) = = {λ λ}_{i i,, j j} \times \times {\overset{^^}{p p}}_{i i,, j j} ((x x,, y the y,, z z,, t t - - {d d}_{i i,, j j} / / c c + + τ τ)) - - - - - - ((99))$

其中F_i,j(x,y,z,τ)是后向传播信号，λ_i,j＝4πd_i,j为补偿目标源信号几何传播产生的信号丢失及补偿被动接收阵元灵敏度的系数，c表示声传播速度，τ表示时间延时；where F_i,j (x,y,z,τ) is the backward propagating signal, λ_i,j =4πd_i,j is the coefficient for compensating the signal loss caused by the geometric propagation of the target source signal and compensating the sensitivity of the passive receiving array element, c represents the speed of sound propagation, τ represents the time delay;

(3.3)计算从散射点(x,y,z)到面阵阵元(i,j)的传播距离d_i,j：(3.3) Calculate the propagation distance d_i,j from the scattering point (x,y,z) to the surface array element (i,j):

${d d}_{i i,, j j} = = \sqrt{{(({x x}_{i i} - - x x))}^{22} + + {(({y the y}_{i i} - - y the y))}^{22} + + {z z}^{22}} - - - - - - ((1010))$

x_i、y_j分别是面阵阵元的横、纵坐标；x_i and y_j are the abscissa and ordinate of the planar array element respectively;

(3.4)估计空化区域散射点(x,y,z)处的源强度E：(3.4) Estimate the source intensity E at the scattering point (x, y, z) in the cavitation region:

${\overset{&OverBar; &OverBar;}{E E.}}_{Source source} ((x x,, y the y,, z z)) = = \frac{22}{{t t}_{22} - - {t t}_{11}} {&Integral; &Integral;}_{{t t}_{11}}^{{t t}_{22}} {{{Σ Σ}_{\underset{{i i}_{11} < < {i i}_{22}}{{i i}_{11},, {i i}_{22} = = 11}}^{n no} {Σ Σ}_{\underset{{j j}_{11} < < {j j}_{22}}{{j j}_{11},, {j j}_{22} = = 11}}^{m m} {F f}_{{i i}_{11},, {j j}_{11}} ((x x,, y the y,, z z,, τ τ)) {F f}_{{i i}_{22},, {j j}_{22}} ((x x,, y the y,, z z,, τ τ))}} dτ dτ - - - - - - ((1111))$

其中，表示源强度的估计，n和m分别为面阵的横向阵元个数和纵向阵元个数，t₁表示时间累计区间下限，t₂表示时间累计区间上限。in, Indicates the estimation of source intensity, n and m are the number of horizontal array elements and vertical array elements respectively, t₁ indicates the lower limit of the time accumulation interval, and t₂ indicates the upper limit of the time accumulation interval.

步骤一中所述应用动态孔径接收技术和幅度变迹技术对源强度三维数据进行处理，具体包括以下步骤：The application of dynamic aperture receiving technology and amplitude apodization technology in step 1 to process the three-dimensional source intensity data includes the following steps:

(4.1)对于面阵每一个阵元(i,j)上的采样点k，确定出采样点的坐标位置(X,Y,Z)，Z＝(k+nSampleOffset)×sampleSpacing，X＝i×pitch，Y＝j×pitch，nSampleOffset表示采样偏移量，sampleSpacing表示采样间隔，pitch表示阵元间隔；(4.1) For the sampling point k on each array element (i, j) of the surface array, determine the coordinate position (X, Y, Z) of the sampling point, Z=(k+nSampleOffset)×sampleSpacing, X=i× pitch, Y=j×pitch, nSampleOffset represents the sampling offset, sampleSpacing represents the sampling interval, and pitch represents the element interval;

(4.2)在焦距与面阵换能器有效直径之比Fnum的基础上计算横向孔径A_x＝Z/(2×Fnum)以及横向半孔径halfA_x＝[A_x/2]_-，并且判断halfA_x是否超过了面阵横向最大孔径maxA_x的一半，如果超过，则横向半孔径halfA_x＝maxA_x/2，从而得到横向孔径下标x＝-halfA_x～halfA_x；在Fnum的基础上计算纵向孔径A_y＝Z/(2×Fnum)以及纵向半孔径halfA_y＝[A_y/2]_-，并且判断halfA_y是否超过了面阵纵向最大孔径maxA_y的一半，如果超过，则纵向半孔径halfA_y＝maxA_y/2，从而得到纵向孔径下标y＝-halfA_y～halfA_y；(4.2) Calculate the lateral aperture A_x =Z/(2×Fnum) and the lateral half aperture halfA_x =[A_x /2]_- on the basis of the ratio Fnum of the focal length and the effective diameter of the area array transducer, and judge halfA Whether_x exceeds half of the maximum lateral aperture maxA_x of the array, if it exceeds, then the lateral half aperture halfA_x = maxA_x /2, thus obtaining the lateral aperture subscript x = -halfA_x ~halfA_x ; calculated on the basis of Fnum Longitudinal aperture A_y =Z/(2×Fnum) and longitudinal semi-aperture halfA_y =[A_y /2]_- , and judge whether halfA_y exceeds half of the maximum longitudinal aperture maxA_y of the array, if so, the longitudinal half Aperture halfA_y = maxA_y /2, thus obtaining the longitudinal aperture subscript y = -halfA_y ~halfA_y ;

(4.3)在采样深度和横纵向孔径下标x,y的基础上计算时间延时：(4.3) Calculate the time delay based on the sampling depth and the horizontal and vertical aperture subscripts x, y:

$delays delays = = {{Z Z + + \sqrt{[[{Z Z}^{22} + + {((X x - - ((i i + + x x)) \times \times pitch pitch))}^{22} + + {((Y Y - - ((j j + + y the y)) \times \times pitch pitch))}^{22}]]}}} / / c c \times \times Fs Fs - - - - - - ((1212))$

其中，delays表示时间延时，c表示声传播速度，Fs表示采样频率；Among them, delays represents the time delay, c represents the speed of sound propagation, and Fs represents the sampling frequency;

(4.4)设置hanning窗进行幅度遍迹，hanning窗的表达式为：(4.4) Set the hanning window for amplitude traversal, the expression of the hanning window is:

$win win ((x x,, y the y)) = = 0.25 0.25 [[11 - - cos cos ((22 π π \frac{x x}{N N - - 11}))]] [[11 - - cos cos ((22 π π \frac{y the y}{N N - - 11}))]],,$

x＝-halfA_x,...,halfA_x,y＝-halfA_y,...,halfA_y (13)x＝-halfA_x ,...,halfA_x ,y＝-halfA_y ,...,halfA_y (13)

其中，N表示单个阵元通道的采样点数；Among them, N represents the number of sampling points of a single array element channel;

(4.5)在已得到时间延时delays和横纵向孔径下标的基础上，从三维数据每条线(x,y)中寻找相应的值填充到预先初始化的数组chnls(x,y)中，从而得到幅度变迹处理后位置(X,Y,Z)的值：(4.5) On the basis of the obtained time delay delays and horizontal and vertical aperture subscripts, from Find the corresponding value in each line (x, y) of the three-dimensional data and fill it into the pre-initialized array chnls(x, y), so as to obtain the value of the position (X, Y, Z) after amplitude apodization processing:

$((X x,, Y Y,, Z Z)) = = {Σ Σ}_{x x = = - - half half {A A}_{x x}}^{{halfA halfA}_{x x}} {Σ Σ}_{y the y = = - - {halfA halfA}_{y the y}}^{{halfA halfA}_{y the y}} chnls chnls ((x x,, y the y)) \times \times win win ((x x,, y the y)) - - - - - - ((1414)) . .$

步骤二中所述将经过步骤一处理后的三维数据转换到参考坐标系下，对多次转换到参考坐标系下的三维数据处理结果取平均得到用于三维重建的预处理数据，具体包括以下步骤：In Step 2, the 3D data processed in Step 1 is converted to the reference coordinate system, and the processing results of the 3D data converted to the reference coordinate system are averaged to obtain the preprocessed data for 3D reconstruction, specifically including the following step:

以下假设参考坐标系为HIFU照射轴向方向的垂直方向，参考坐标系坐标为(x,y,z)，旋转坐标系坐标为(x',y',z')，旋转坐标系的方向向量为(v'_x,v'_y,v'_z)，平移矩阵用T_arb表示，旋转矩阵用R_arb表示，转换矩阵用M_arb表示；The following assumes that the reference coordinate system is the vertical direction of the HIFU irradiation axial direction, the coordinates of the reference coordinate system are (x, y, z), the coordinates of the rotating coordinate system are (x', y', z'), and the direction vector of the rotating coordinate system is (v'_x , v'_y , v'_z ), the translation matrix is represented by T_arb , the rotation matrix is represented by R_arb , and the transformation matrix is represented by M_arb ;

(5.1)计算平移矩阵：(5.1) Calculate the translation matrix:

${T T}_{arb arb} = = [\begin{matrix} 11 & 00 & 00 & 00 \\ 00 & 11 & 00 & 00 \\ 00 & 00 & 11 & 00 \\ {Δ Δ}_{x x} & {Δ Δ}_{y the y} & {Δ Δ}_{z z} & 11 \end{matrix}] - - - - - - ((1515))$

其中，Δ_x、Δ_y、Δ_z分别表示参考坐标系原点和旋转坐标系原点在三个坐标轴之间的平移量；Among them, Δ_x , Δ_y , and Δ_z represent the translation between the origin of the reference coordinate system and the origin of the rotating coordinate system between the three coordinate axes;

根据坐标系的方向向量构造的旋转矩阵为：The rotation matrix constructed according to the direction vector of the coordinate system is:

${R R}_{arb arb} = = [\begin{matrix} {v v}_{{x x}_{11}}^{' '} & {v v}_{{y the y}_{11}}^{' '} & {v v}_{{z z}_{11}}^{' '} & 00 \\ {v v}_{{x x}_{22}}^{' '} & {v v}_{{y the y}_{22}}^{' '} & {v v}_{{z z}_{22}}^{' '} & 00 \\ {v v}_{{x x}_{33}}^{' '} & {v v}_{{y the y}_{33}}^{' '} & {v v}_{{z z}_{33}}^{' '} & 00 \\ 00 & 00 & 00 & 11 \end{matrix}] - - - - - - ((1616))$

(v'_xi,v'_yi,v'_zi)为旋转坐标系的方向向量(v'_x,v'_y,v'_z)分解到参考坐标系上所得到的分解向量，i＝1～3，即(v'_xi , v'_yi , v'_zi ) is the decomposition vector obtained by decomposing the direction vector (v'_x , v'_y , v'_z ) of the rotating coordinate system into the reference coordinate system, i=1～3 ,Right now

${v v}_{x x}^{' '} = = [[{v v}_{{x x}_{11}}^{' '} (({v v}_{x x},, {v v}_{y the y},, {v v}_{z z})),, {v v}_{{x x}_{22}}^{' '} (({v v}_{x x},, {v v}_{y the y},, {v v}_{z z})),, {v v}_{{x x}_{33}}^{' '} (({v v}_{x x},, {v v}_{y the y},, {v v}_{z z}))]]$

${v v}_{y the y}^{' '} = = [[{v v}_{{y the y}_{11}}^{' '} (({v v}_{x x},, {v v}_{y the y},, {v v}_{z z})),, {v v}_{{y the y}_{22}}^{' '} (({v v}_{x x},, {v v}_{y the y},, {v v}_{z z})),, {v v}_{{y the y}_{33}}^{' '} (({v v}_{x x},, {v v}_{y the y},, {v v}_{z z}))]]$

${v v}_{z z}^{' '} = = [[{v v}_{{z z}_{11}}^{' '} (({v v}_{x x},, {v v}_{y the y},, {v v}_{z z})),, {v v}_{{z z}_{22}}^{' '} (({v v}_{x x},, {v v}_{y the y},, {v v}_{z z})),, {v v}_{{z z}_{33}}^{' '} (({v v}_{x x},, {v v}_{y the y},, {v v}_{z z}))]]$

参考坐标系经过平移和旋转后得坐标转换矩阵：The coordinate transformation matrix obtained after translation and rotation of the reference coordinate system:

M_arb＝T_arb(Δ_x,Δ_y,Δ_z)R_arb(v'_x,v'_y,v'_z) (17)M_arb ＝T_arb (Δ_x ,Δ_y ,Δ_z )R_arb (v'_x ,v'_y ,v'_z ) (17)

将旋转坐标系下的数据对应到参考坐标系下的数据，求解方程(18)，得用于三维重建的预处理三维体数据E_pre3D(x,y,z)：Corresponding the data in the rotating coordinate system to the data in the reference coordinate system, and solving equation (18), the preprocessed 3D volume data E_pre3D (x, y, z) for 3D reconstruction can be obtained:

${E E.}_{pre pre 33 D D.} ((x x,, y the y,, z z)) = = [[{x x}^{' '},, {y the y}^{' '},, {z z}^{' '} 11]] {R R}_{arb arb}^{- - 11} (({v v}_{x x}^{' '},, {v v}_{y the y}^{' '} {v v}_{z z}^{' '})) {T T}_{arb arb}^{- - 11} (({Δ Δ}_{x x},, {Δ Δ}_{y the y},, {Δ Δ}_{z z})) - - - - - - ((1818))$

(5.2)用来表示第q个角度下得到的预处理三维体数据，重复Q次设置接收阵元时间延时，并得到Q个角度下的预处理三维体数据，然后对Q个角度下的预处理三维体数据进行平均，则可得到Q次平均后的预处理三维体数据：(5.2) with to represent the preprocessed 3D volume data obtained at the qth angle, repeat Q times to set the time delay of receiving array elements, and obtain the preprocessed 3D volume data at Q angles, and then perform the preprocessing 3D volume data at Q angles If the data is averaged, the preprocessed three-dimensional volume data after Q averages can be obtained:

${\overset{&OverBar; &OverBar;}{E E.}}_{pre pre 33 D D.}^{Q Q} ((x x,, y the y,, z z)) = = \frac{11}{Q Q} {Σ Σ}_{q q = = 11}^{Q Q} {E E.}_{pre pre 33 D D.}^{((q q))} ((x x,, y the y,, z z)) - - - - - - ((1919))$

(5.3)第Q+1次平均后的预处理三维体数据为：(5.3) The preprocessed three-dimensional volume data after the Q+1 average is:

${\overset{&OverBar; &OverBar;}{E E.}}_{pre pre 33 D D.}^{Q Q + + 11} ((x x,, y the y,, z z)) = = \frac{Q Q}{Q Q + + 11} {\overset{&OverBar; &OverBar;}{E E.}}_{pre pre 33 D D.}^{Q Q} ((x x,, y the y,, z z)) + + \frac{11}{Q Q + + 11} {E E.}_{pre pre 33 D D.}^{((Q Q + + 11))} ((x x,, y the y,, z z)) - - - - - - ((2020))$

其中为Q次平均后的预处理三维体数据，为第Q+1个角度下的预处理三维体数据。in is the preprocessed 3D volume data after Q averages, is the preprocessed three-dimensional volume data at the Q+1th angle.

步骤三中所述对所述预处理数据分割子区域并设置径向基函数和窗函数，具体包括以下步骤：As described in step 3, the sub-regions are divided into the preprocessed data and the radial basis function and window function are set, which specifically includes the following steps:

(6.1)假设体数据点子区域C_k被数据盒D_k所包围，|C_k|为体数据点子区域点集内的点数，假设S_max和s为子区域可允许的最大数据点数和最小尺寸；(6.1) Assume that the volume data point sub-area C_k is surrounded by the data box D_k , |C_k | is the number of points in the volume data point sub-area point set, assuming S_max and s are the maximum number of data points and the minimum size allowed by the sub-area ;

(6.2)如果满足条件：|C_k|>S_max且沿着x轴方向的尺寸大于最小边长，则计算C_k内所有点的重心，记为(x_w,y_w,z_w)，然后用平面x＝x_w分割点集，如此往复循环直到子区域内的数据点数|C_k|和尺寸满足该条件；(6.2) If the condition is met: |C_k |>S_max and the size along the x-axis direction is greater than the minimum side length, then calculate the center of gravity of all points in C_k , recorded as (x_w , y_w , z_w ), Then use the plane x=_xw to divide the point set, so reciprocate until the number of data points | C_k | and the size in the sub-region meet the condition;

(6.3)然后把相邻的两个数据盒做并集运算，使两个数据盒相交，并减去公共的部分，假如相并之后的数据盒不满足条件：|C_k|>S_max且沿着y轴方向的尺寸大于最小边长，则用平面y＝y_w分割点集直到满足该条件；(6.3) Then perform a union operation on two adjacent data boxes, make the two data boxes intersect, and subtract the common part, if the combined data box does not meet the condition: |C_k |>S_max and The size along the y-axis direction is greater than the minimum side length, then use the plane y=y_w to divide the point set until the condition is met;

(6.4)接着继续把相邻的两个数据盒做并集运算，使两个数据盒相交，并减去公共的部分，假如相并之后的数据盒不满足条件：|C_k|>S_max且沿着z轴方向的尺寸大于最小边长，则用平面z＝z_w分割点集直到满足该条件；(6.4) Then continue to do the union operation of two adjacent data boxes, make the two data boxes intersect, and subtract the common part, if the combined data box does not meet the condition: |C_k |>S_max And the size along the z-axis direction is greater than the minimum side length, then use the plane z=z_w to divide the point set until the condition is met;

(6.5)设置第k个子区域径向基插值时所需的窗函数：(6.5) Set the window function required for radial basis interpolation in the kth sub-region:

${ω ω}_{k k} ((r r)) = = ((11 - - \frac{| | 22 x x - - (({x x}_{k k max max} + + {x x}_{k k min min})) | |}{{x x}_{k k max max} - - {x x}_{k k min min}})) ((11 - - \frac{| | 22 y the y - - (({y the y}_{k k max max} + + {y the y}_{k k min min})) | |}{{y the y}_{k k max max} - - {y the y}_{k k min min}})) ((11 - - \frac{| | 22 z z (({z z}_{k k max max} + + {z z}_{k k min min})) | |}{{z z}_{k k max max} - - {z z}_{k k min min}})) - - - - - - ((21 twenty one))$

其中x_kmin,x_kmax,y_kmin,y_kmax,z_kmin,z_kmax,分别为数据盒D_k在三维尺度x,y,z上边界的下限和上限；Among them, x_kmin , x_kmax , y_kmin , y_kmax , z_kmin , z_kmax are respectively the lower limit and upper limit of the upper boundary of the data box D_k on the three-dimensional scale x, y, z;

(6.6)假设第k个子区域数据点集C_k＝{r_i(i＝1,2,...,|C_k|)}，所对应的函数值为f(X_i)(i＝1,2,...,|C_k|)，则得方程组：(6.6) Assuming the kth sub-region data point set C_k ={r_i (i=1,2,...,|C_k |)}, the corresponding function value is f(X_i )(i=1 ,2,...,|C_k |), then the system of equations is obtained:

$[\begin{matrix} {a a}_{11} & {a a}_{22} & . . . . . . & {a a}_{| | {C C}_{k k} | |} \end{matrix}] [\begin{matrix} φ φ ((| | | | {X x}_{i i} - - {X x}_{11} | | | |)) \\ φ φ ((| | | | {X x}_{i i} - - {X x}_{11} | | | |)) \\ . . . . . . \\ φ φ ((| | | | {X x}_{i i} - - {X x}_{j j} | | | |)) \end{matrix}] = = f f (({X x}_{i i})) ((i i = = 1,2 1,2,, . . . . . .,, | | {C C}_{k k} | |)) - - - - - - ((22 twenty two))$

(6.7)解方程组得到的值，从而计算出所有子区域的径向插值曲面方程，利用此曲面方程求得网格点对应的函数值，从而取得B样条插值的点集B_ij(i＝1,2,...,M,j＝1,2,...,N)。(6.7) Solve the system of equations to get value, so as to calculate the radial interpolation surface equation of all sub-regions, and use this surface equation to obtain the function value corresponding to the grid point, so as to obtain the B-spline interpolation point set B_ij (i=1,2,.. .,M,j=1,2,...,N).

步骤三中所述结合径向基函数插值和B样条插值方法对所述预处理数据进行三维重建，具体包括以下步骤：The three-dimensional reconstruction of the preprocessed data in combination with radial basis function interpolation and B-spline interpolation methods described in step three specifically includes the following steps:

(7.1)假设经过处理得到的预处理三维体数据点集合为C，设C₀为集合C在XOY平面上的投影点集，并划定出C₀的边界；(7.1) Assuming that the preprocessed three-dimensional volume data point set obtained through processing is C, let C be the projection point set of set C on the XOY plane, and delineate the boundary of C₀_;

(7.2)将投影点集C₀分割成M×N个规则的小矩形网格，相邻矩形的交点则对应着B样条插值的点在XOY平面上的投影；(7.2) Divide the projection point set_C0 into M×N regular small rectangular grids, and the intersection points of adjacent rectangles correspond to the projection of the B-spline interpolation points on the XOY plane;

(7.3)将集合C划分成L个子区域，每一块子区域C_k数据点的个数为|C_k|个，k＝1,2,...,L；(7.3) Divide the set C into L sub-regions, the number of data points in each sub-region C_k is |C_k | pieces, k=1,2,...,L;

(7.4)设置不同的径向基函数rbf_k对每一个子区域C_k内的数据点分别插值，得到插值后的三维子区域；(7.4) Different radial basis functions rbf_k are set to interpolate the data points in each sub-area C_k respectively to obtain a three-dimensional sub-area after interpolation;

(7.5)根据每个子区域插值函数rbf_k分别求出C₀划分后网格点对应的函数值，组成了B样条插值的点集B_ij(i＝1,2,...,M,j＝1,2,...,N)；(7.5) According to the interpolation function rbf_k of each sub-area, the function values corresponding to the grid points after_C0 division are respectively obtained, and the point set B_ij of B-spline interpolation (i=1,2,...,M, j=1,2,...,N);

(7.6)对点集B_ij做B样条插值运算，生成B样条插值后的三维成像。(7.6) Perform B-spline interpolation operation on point set B_ij to generate 3D imaging after B-spline interpolation.

本发明的有益效果体现在：本发明是一种在空间中对空化微泡进行被动成像和三维复合成像的方法。其将一般的被动成像算法由二维延伸到三维，并针对FBM重建算法在过大计算量和方程求解不稳定等方面的问题，引入了区域分割预处理和B样条插值方法。综合以上方法提出的三维被动空化成像方法对于提高三维成像的分辨率、信噪比和效率有着重要的意义。The beneficial effect of the present invention is embodied in that the present invention is a method for passive imaging and three-dimensional composite imaging of cavitation microbubbles in space. It extends the general passive imaging algorithm from two-dimensional to three-dimensional, and introduces the region segmentation preprocessing and B-spline interpolation method to solve the problems of the FBM reconstruction algorithm in the aspects of excessive calculation and unstable equation solution. Combining the above methods, the 3D passive cavitation imaging method proposed is of great significance for improving the resolution, signal-to-noise ratio and efficiency of 3D imaging.

本发明具有以下优点：The present invention has the following advantages:

1、三维被动空化成像方法拓展了二维的PCI成像方法，引入的三维复合成像可以使图像的信噪比进一步提高，数据的实时更新避免了过大数据量无法存储的问题，可以得到高分辨率、高信噪比和高效率的三维空化成像，为临床提供了实时三维监控的可能。1. The three-dimensional passive cavitation imaging method expands the two-dimensional PCI imaging method, and the introduction of three-dimensional composite imaging can further improve the signal-to-noise ratio of the image. Three-dimensional cavitation imaging with high resolution, high signal-to-noise ratio and high efficiency provides the possibility of real-time three-dimensional monitoring for clinical practice.

2、根据自适应滤波理论设计的多级连自适应陷波系统相比一般的数字带阻滤波器更容易控制带宽和品质因数，保留的宽带噪声更加可信，提供了研究惯性空化机制的方法。2. The multi-stage adaptive notch system designed according to the adaptive filtering theory is easier to control the bandwidth and quality factor than the general digital band-stop filter, and the retained broadband noise is more reliable, which provides a basis for studying the inertial cavitation mechanism method.

3、空化三维体数据的重建引入了子区域分割预处理，极大地降低了计算量，提高了计算效率，结合径向基函数插值和B样条插值使得三维成像更加稳定、图像更加平滑，且可得到实时更新的三维被动空化图像。3. The reconstruction of cavitation 3D volume data introduces sub-region segmentation preprocessing, which greatly reduces the amount of calculation and improves calculation efficiency. Combined with radial basis function interpolation and B-spline interpolation, the 3D imaging is more stable and the image is smoother. And real-time updated three-dimensional passive cavitation images can be obtained.

4、三维被动成像方法不仅可以应用于空化，也可以应用于造影微泡的监控、基因治疗、药物传递、组织消融等诸多方面。4. The three-dimensional passive imaging method can be applied not only to cavitation, but also to the monitoring of contrast microbubbles, gene therapy, drug delivery, tissue ablation and many other aspects.

附图说明Description of drawings

图1为三维被动空化成像方法中多级连自适应陷波系统框图；Figure 1 is a block diagram of a multi-stage adaptive notch system in a three-dimensional passive cavitation imaging method;

图2为三维被动空化成像方法中成像算法延伸框图；Fig. 2 is an extended block diagram of the imaging algorithm in the three-dimensional passive cavitation imaging method;

图3为三维被动空化成像方法中三维复合方法示意图；Fig. 3 is a schematic diagram of the three-dimensional composite method in the three-dimensional passive cavitation imaging method;

图4为三维被动空化成像方法中B样条插值三维重建仿真图；Fig. 4 is a B-spline interpolation three-dimensional reconstruction simulation diagram in the three-dimensional passive cavitation imaging method;

图5为三维被动空化成像方法中径向基函数插值三维重建仿真实例图；Fig. 5 is an example diagram of radial basis function interpolation 3D reconstruction simulation in the 3D passive cavitation imaging method;

图6为三维被动空化成像方法中三维重建程序流程图。Fig. 6 is a flow chart of the three-dimensional reconstruction program in the three-dimensional passive cavitation imaging method.

具体实施方式Detailed ways

下面结合附图和实施例对本发明做详细描述。The present invention will be described in detail below in conjunction with the accompanying drawings and embodiments.

本发明涉及到的三维被动空化成像方法包括以下具体步骤：The three-dimensional passive cavitation imaging method involved in the present invention includes the following specific steps:

步骤一：基于谱估计理论求出信号的功率谱估计，根据面阵采集到的三维数据的最大幅值设定阈值，以明确要滤除的频率分量从而减少计算时间。设计多级连自适应陷波系统，对面阵接收的n×m个通道逐次滤除掉其中的高次谐波和超高次谐波分量。Step 1: Based on the spectrum estimation theory, the power spectrum estimation of the signal is obtained, and the threshold is set according to the maximum amplitude of the three-dimensional data collected by the area array, so as to clarify the frequency components to be filtered and reduce the calculation time. A multi-stage adaptive notch system is designed to successively filter out the high-order harmonic and ultra-high-order harmonic components of the n×m channels received by the area array.

步骤二：计算散射点与阵元之间的距离和对信号丢失和接收灵敏度的补偿系数，采用被动成像算法计算时间累计区间内的源强度三维数据。Step 2: Calculate the distance between the scattering point and the array element and the compensation coefficient for signal loss and receiving sensitivity, and use the passive imaging algorithm to calculate the three-dimensional source intensity data in the time accumulation interval.

步骤三：提取步骤二中得到的源强度三维数据，应用动态孔径接收技术和幅度变迹技术对源强度三维数据进行处理，得到惯性空化的三维体数据。Step 3: Extract the 3D source intensity data obtained in Step 2, and apply dynamic aperture receiving technology and amplitude apodization technology to process the 3D source intensity data to obtain 3D volume data of inertial cavitation.

步骤四：重复设置面阵阵元时间延时(5～10次)以调整成像角度(5°～20°)，重新采集三维数据并按照步骤一至步骤三所述方法处理，并将经过处理后的三维数据转换到参考坐标系下，多次平均得到用于三维重建的预处理数据；Step 4: Repeatedly set the time delay of the area array element (5 to 10 times) to adjust the imaging angle (5° to 20°), re-acquire the 3D data and process it according to the method described in Step 1 to Step 3, and the processed Transform the 3D data into the reference coordinate system, and average multiple times to obtain the preprocessed data for 3D reconstruction;

步骤五：对步骤四得到的预处理数据进行子区域的分割以提高三维重建速度，选择径向基函数对子区域插值并计算B样条插值点集。Step 5: Segment the preprocessed data obtained in step 4 into sub-regions to increase the speed of 3D reconstruction, select the radial basis function to interpolate the sub-regions and calculate the B-spline interpolation point set.

步骤六：在步骤五子区域分割和B样条插值点集计算完成的基础上，做B样条插值运算，生成B样条插值后的三维重建成像。Step 6: On the basis of the sub-region segmentation and B-spline interpolation point set calculation in step 5, perform B-spline interpolation operation to generate a 3D reconstruction image after B-spline interpolation.

具体说明如下：The specific instructions are as follows:

步骤一、求出信号功率谱估计并且确定陷波频率的大小及个数，基于LMS算法设计多级连自适应陷波系统，遍历采集到的每一个通道的射频数据以及遍历每一个陷波频率进行自适应陷波处理，其具体步骤如下：Step 1. Find the signal power spectrum estimate and determine the size and number of notch frequencies, design a multi-stage adaptive notch system based on the LMS algorithm, traverse the collected RF data of each channel and traverse each notch frequency Perform adaptive notch processing, the specific steps are as follows:

参阅图1所示，描述如下：Referring to Figure 1, the description is as follows:

(1)对第一个通道的信号进行快速傅里叶变换，得到信号的功率谱估计，明确高次谐波信号和超高次谐波信号。(1) Perform fast Fourier transform on the signal of the first channel to obtain the power spectrum estimation of the signal, and clarify the high-order harmonic signal and ultra-high-order harmonic signal.

(2)设定阈值为thrs＝max{x(1),x(2),...,x(N)}×r，认为幅值大于thrs的频率分量都是需要滤除的，这一步的目的是为了加快计算速度。其中r为经验值，大约在1/30～1/20之间。(2) Set the threshold as thrs=max{x(1),x(2),...,x(N)}×r, and consider that the frequency components with amplitude greater than thrs need to be filtered out. This step The purpose is to speed up the calculation. Among them, r is an empirical value, which is about 1/30 to 1/20.

(3)根据设定的阈值确定出接下来自适应滤波要滤除的频率成分，要滤除的频率分量可能为nf₀/2(n＝1,2,...)其中的若干个。(3) Determine the frequency components to be filtered out in the next adaptive filtering according to the set threshold, and the frequency components to be filtered out may be several of nf₀ /2 (n=1, 2, . . . ).

(4)对第一个通道的信号进行快速傅里叶变换，得到空化信号的功率谱估计，设置参考信号为基本的正弦信号，正弦信号的相位均为0，幅度C均为1，频率f为第一个谐波频率。参考信号分为两路，一路是正弦信号本身x₁，另外一路是正弦信号经过一个正交相移网络后的信号x₂。表达式分别为：(4) Perform fast Fourier transform on the signal of the first channel to obtain the power spectrum estimation of the cavitation signal, set the reference signal as a basic sinusoidal signal, the phase of the sinusoidal signal is 0, the amplitude C is 1, and the frequency f is the first harmonic frequency. The reference signal is divided into two channels, one is the sinusoidal signal itself x₁ , and the other is the signal x₂ of the sinusoidal signal after passing through a quadrature phase shift network. The expressions are:

x₁＝C cos(2πft),x₂＝C sin(2πft) (9)x₁ ＝C cos(2πft), x₂ ＝C sin(2πft) (9)

(5)设置自适应陷波器中的权阶数为2，权系数的初始值分别为w₁＝0.1、w₂＝0.1，对(4)中两路信号进行加权得到加权信号y＝w₁x₁(i)+w₂x₂(i)，计算误差信号第i个采样点对应的数值e(i)＝x(i)-y，其中x(i)为第一个阵元通道的原始信号，x₁(i)为第一路参考信号第i个采样点对应数值，x₂(i)为第二路参考信号第i个采样点对应数值；(5) Set the weight order in the adaptive notch filter to 2, and the initial values of the weight coefficients are respectively w₁ =0.1, w₂ =0.1, and weight the two signals in (4) to obtain a weighted signal y=w₁ x₁ (i)+w₂ x₂ (i), calculate the value e(i)=x(i)-y corresponding to the i-th sampling point of the error signal, where x(i) is the first array element channel The original signal of , x₁ (i) is the value corresponding to the i-th sampling point of the first reference signal, and x₂ (i) is the corresponding value of the i-th sampling point of the second reference signal;

(6)设置自适应陷波的收敛系数μ，根据Wirdow-Hoff提出的最陡下降法产生两个权系数的迭代公式(6) Set the convergence coefficient μ of the adaptive notch, and generate an iterative formula of two weight coefficients according to the steepest descent method proposed by Wirdow-Hoff

w₁＝w₁+2μe(i)x₁(i),w₂＝w₂+2μe(i)x₂(i) (10)w₁ =w₁ +2μe(i)x₁ (i),w₂ =w₂ +2μe(i)x₂ (i) (10)

(7)对于第一个通道的所有采样点x(1),x(2),...,x(N)，重复步骤(5)和(6)的过程，其中N为采样点数，得到误差信号e，误差信号即为滤除了第一个频率分量的信号。(7) For all sampling points x(1), x(2),...,x(N) of the first channel, repeat steps (5) and (6), where N is the number of sampling points, and get Error signal e, the error signal is the signal from which the first frequency component has been filtered.

(8)在(7)的基础上级联下一个自适应陷波器，改变陷波频率，重复步骤(4)～(7)的过程。得到第一个通道将高谐波和超高次谐波滤除掉之后的宽带噪声信号，单个自适应陷波器的传递函数为：(8) On the basis of (7), cascade the next adaptive notch filter, change the notch frequency, and repeat the process of steps (4) to (7). Obtain the broadband noise signal after the first channel filters out high harmonics and ultra-high harmonics, and the transfer function of a single adaptive notch filter is:

$H h ((z z)) = = \frac{{z z}^{22} - - 22 z z cos cos wT wxya + + 11}{{z z}^{22} - - 22 z z cos cos wT wxya ((11 - - {μC μC}^{22})) + + ((11 - - 22 μ μ {C C}^{22}))} - - - - - - ((1111))$

(9)重复步骤(4)～(8)，使面阵每个阵元采集到的射频数据经过多级连自适应陷波系统，得到一组反映微泡惯性振动机制的射频数据。(9) Repeat steps (4) to (8), so that the radio frequency data collected by each array element of the area array passes through the multi-stage adaptive notch system to obtain a set of radio frequency data reflecting the inertial vibration mechanism of the microbubbles.

步骤二、计算散射点与阵元之间的距离和对信号丢失和接收灵敏度的补偿系数，采用被动成像算法计算时间累计区间内的源强度三维数据：Step 2. Calculate the distance between the scattering point and the array element and the compensation coefficient for signal loss and receiving sensitivity, and use the passive imaging algorithm to calculate the three-dimensional source intensity data in the time accumulation interval:

(1)假设有一个零时间平均值的空化源，根据线性波动理论，可得产生的压力场p(x,y,z,t)，即：(1) Assuming that there is a cavitation source with zero time average value, according to the linear wave theory, the resulting pressure field p(x,y,z,t) can be obtained, namely:

${&dtri; &dtri;}^{22} p p - - \frac{11}{{c c}^{22}} \frac{{&PartialD; &PartialD;}^{22} p p}{{&PartialD; &PartialD; t t}^{22}} = = - - s the s ((x x,, y the y,, z z,, t t)) - - - - - - ((1212))$

(2)通过面阵被动接收到的射频压力信号记为p_i,j，这里可以对射频信号进行相关的滤波处理以满足不同条件下的成像需求，即(2) The radio frequency pressure signal passively received by the area array is denoted as p_i,j , where relevant filtering can be performed on the radio frequency signal to meet the imaging requirements under different conditions, namely

${\overset{^^}{p p}}_{i i,, j j} ((x x,, y the y,, z z,, t t)) = = {p p}_{i i,, j j} ((x x,, y the y,, z z,, t t)) * * F f ((t t)) - - - - - - ((1313))$

其中F(t)为滤波函数，与p_i,j(x,y,z,t)在时域上做卷积运算。这个滤波器的设计可以是陷波器以滤除掉谐波成分来呈现惯性的空化，也可以是带通滤波器以提取次谐波成分来得到非惯性成像。Among them, F(t) is a filter function, which performs convolution operation with p_i,j (x,y,z,t) in the time domain. The design of this filter can be a notch filter to filter out harmonic components to present inertial cavitation, or a bandpass filter to extract sub-harmonic components to obtain non-inertial imaging.

(3)对做时间延时和补偿(3) yes Do time delay and compensation

${F f}_{i i,, j j} ((x x,, y the y,, z z,, τ τ)) = = {λ λ}_{i i,, j j} \times \times {\overset{^^}{p p}}_{i i,, j j} ((x x,, y the y,, z z,, t t - - {d d}_{i i,, j j} / / c c + + τ τ)) - - - - - - ((1414))$

其中F_i,j(x,y,z,τ)是后向传播信号，λ_i,j＝4πd_i,j为补偿目标源信号几何传播产生的信号丢失及补偿被动接收阵元灵敏度的系数，c表示声传播速度，τ时间延时；where F_i,j (x,y,z,τ) is the backward propagating signal, λ_i,j =4πd_i,j is the coefficient for compensating the signal loss caused by the geometric propagation of the target source signal and compensating the sensitivity of the passive receiving array element, c represents the speed of sound propagation, τ time delay;

(4)计算从散射点(x,y,z)到成像面阵阵元(i,j)的传播距离d_i,j：(4) Calculate the propagation distance d_i,j from the scattering point (x,y,z) to the imaging surface array element (i,j):

${d d}_{i i,, j j} = = \sqrt{{(({x x}_{i i} - - x x))}^{22} + + {(({y the y}_{i i} - - y the y))}^{22} + + {z z}^{22}} - - - - - - ((1515))$

(5)估计空化区域位置(x,y,z)处的源强度E。空化产生的信号是均值为零的随机变量，对其进行统计分析可得E{p_i,j}＝0。由于被动成像数据在压力场中是线性的，在“幅值”意义上得到成像显然无法实现，因此考虑在“强度”意义上考虑对空化区域的成像。然而，这会导致在一段时间的积累上产生很大的直流成分，因此作为对空化成像的估计，需要将直流成分去除掉。设置时间积累区间并计算后向传播阵列中每个阵元的射频压力信号，得到源强度为：(5) Estimate the source intensity E at the location (x, y, z) of the cavitation region. The signal generated by cavitation is a random variable with a mean value of zero. Statistical analysis on it can give E{p_i,j }=0. Since the passive imaging data is linear in the pressure field, imaging in the sense of "amplitude" is obviously impossible, so consider imaging the cavitation region in the sense of "intensity". However, this leads to a large DC component accumulated over a period of time, so as an estimate for cavitation imaging, the DC component needs to be removed. Set the time accumulation interval and calculate the RF pressure signal of each element in the backpropagation array, and the source strength is obtained as:

${\overset{&OverBar; &OverBar;}{E E.}}_{Source source} ((x x,, y the y,, z z)) = = \frac{11}{{t t}_{22} - - {t t}_{11}} {&Integral; &Integral;}_{{t t}_{11}}^{{t t}_{22}} {{{[[{Σ Σ}_{i i = = 11}^{n no} {Σ Σ}_{j j = = 11}^{m m} {F f}_{i i,, j j} ((x x,, y the y,, z z,, τ τ))]]}^{22} - - {Σ Σ}_{i i = = 11}^{n no} {Σ Σ}_{j j = = 11}^{m m} {F f}_{i i,, j j}^{22} ((x x,, y the y,, z z,, τ τ))}} dτ dτ - - - - - - ((1616))$

其中，表示源强度的估计，n和m分别为面阵的横向阵元个数和纵向阵元个数，t₁表示时间累计区间下限，t₂表示时间累计区间上限，第一项保留了空化成像中相干源的相似性从而可以得到焦点以外延伸区域源强度的计算，第二项则是滤除了成像中的直流成分。in, Indicates the estimation of the source intensity, n and m are the number of transverse array elements and longitudinal array elements respectively, t₁ indicates the lower limit of the time accumulation interval, t₂ indicates the upper limit of the time accumulation interval, the first item retains the cavitation imaging The similarity of the coherent source in the center can obtain the calculation of the source intensity in the extended area beyond the focal point, and the second term is to filter out the DC component in the imaging.

(6)降低计算量，使用恒等式：(6) To reduce the amount of calculation, use the identity:

${(({Σ Σ}_{i i = = 11}^{n no} {k k}_{i i}))}^{22} - - {Σ Σ}_{i i = = 11}^{n no} {k k}_{i i}^{22} = = 22 {\underset{i i < < j j}{Σ Σ}}_{i i,, j j = = 11}^{n no} {k k}_{i i} {k k}_{j j} - - - - - - ((1717))$

将上式修改为：Modify the above formula to:

${\overset{&OverBar; &OverBar;}{E E.}}_{Source source} ((x x,, y the y,, z z)) = = \frac{22}{{t t}_{22} - - {t t}_{11}} {&Integral; &Integral;}_{{t t}_{11}}^{{t t}_{22}} {{{Σ Σ}_{\underset{{i i}_{11} < < {i i}_{22}}{{i i}_{11},, {i i}_{22} = = 11}}^{n no} {Σ Σ}_{\underset{{j j}_{11} < < {j j}_{22}}{{j j}_{11},, {j j}_{22} = = 11}}^{m m} {F f}_{{i i}_{11},, {j j}_{11}} ((x x,, y the y,, z z,, τ τ)) {F f}_{{i i}_{22},, {j j}_{22}} ((x x,, y the y,, z z,, τ τ))}} dτ dτ - - - - - - ((1818))$

(7)被动空化算法的演变示意图参阅图2所示。(7) The schematic diagram of the evolution of the passive cavitation algorithm is shown in Figure 2.

步骤三、在步骤二处理的基础上应用动态孔径接收技术和幅度变迹技术对源强度三维数据进行处理，其具体步骤为如下：Step 3. Apply dynamic aperture receiving technology and amplitude apodization technology to process the three-dimensional source intensity data on the basis of the processing in step 2. The specific steps are as follows:

(1)设置数据尺寸、采样点数、阵元间隔pitch，声传播速度、采样间隔sampleSpacing和采样偏移量nSampleOffset以便下一步计算；(1) Set the data size, number of sampling points, array element pitch, sound propagation velocity, sampling interval sampleSpacing and sampling offset nSampleOffset for the next calculation;

(2)对于每一条线(i,j)上的采样点k，确定出采样点的坐标位置(X,Y,Z)。其中Z＝(k+nSampleOffset)×sampleSpacing，X＝i×pitch，Y＝j×pitch；(2) For the sampling point k on each line (i, j), determine the coordinate position (X, Y, Z) of the sampling point. Where Z=(k+nSampleOffset)×sampleSpacing, X=i×pitch, Y=j×pitch;

(3)在Fnum(焦距与换能器有效直径之比)的基础上计算横向孔径A_x＝Z/(2×Fnum)以及横向半孔径halfA_x＝[A_x/2]_-，[x]_-表示取靠近x的最小整数，并且判断halfA_x是否超过了面阵横向最大孔径maxA_x的一半，如果超过，则横向半孔径halfA_x＝maxA_x/2，从而得到横向孔径下标x＝-halfA_x～halfA_x；在Fnum的基础上计算纵向孔径A_y＝Z/(2×Fnum)以及纵向半孔径halfA_y＝[A_y/2]_-，并且判断halfA_y是否超过了面阵纵向最大孔径maxA_y的一半，如果超过，则纵向半孔径halfA_y＝maxA_y/2，从而得到纵向孔径下标y＝-halfA_y～halfA_y；(3) Calculate the transverse aperture A_x =Z/(2×Fnum) and the transverse semi-aperture halfA_x =[A_x /2]_- on the basis of Fnum (the ratio of the focal length to the effective diameter of the transducer), [x]_- means to take the smallest integer close to x, and judge whether halfA_x exceeds half of the maximum lateral aperture maxA_x of the array, if so, then the lateral half aperture halfA_x = maxA_x /2, thus obtaining the lateral aperture subscript x = - halfA_x ～halfA_x ; Calculate the longitudinal aperture A_y = Z/(2×Fnum) and the longitudinal semi-aperture halfA_y = [A_y /2]_- on the basis of Fnum, and judge whether halfA_y exceeds the longitudinal maximum of the area array Half of the aperture maxA_y , if it exceeds, then the longitudinal half aperture halfA_y = maxA_y /2, thus obtaining the longitudinal aperture subscript y = -halfA_y ~halfA_y ;

(4)在采样深度和横纵向孔径下标x,y的基础上计算延时：(4) Calculate the delay on the basis of the sampling depth and the horizontal and vertical aperture subscripts x, y:

$delays delays = = {{Z Z + + \sqrt{[[{Z Z}^{22} + + {((X x - - ((i i + + x x)) \times \times pitch pitch))}^{22} + + {((Y Y - - ((j j + + y the y)) \times \times pitch pitch))}^{22}]]}}} / / c c \times \times Fs Fs - - - - - - ((1919))$

其中，delays表示时间延时，c表示声传播速度，Fs表示采样频率，这里的延时用采样点来表示方便后续的计算；Among them, delays represents the time delay, c represents the speed of sound propagation, Fs represents the sampling frequency, and the delay here is represented by sampling points to facilitate subsequent calculations;

(5)幅度变迹技术是一种控制发射和接收声场分布的手段。当发射子阵中各个阵元施加相同幅度的激励信号，就形成了声场中的等幅度相干叠加。一般地，幅度变迹可以使子阵中中心阵元的激励信号幅度强，而两旁位置阵元的激励信号幅度逐渐减弱。常用的幅度变迹函数有hanning函数、hamming函数、Blackman函数等。选择hanning窗进行幅度遍迹，hanning窗的表达式为：(5) Amplitude apodization technology is a means to control the distribution of transmitting and receiving sound fields. When excitation signals of the same amplitude are applied to each array element in the transmitting sub-array, an equal-amplitude coherent superposition in the sound field is formed. Generally, amplitude apodization can make the excitation signal amplitude of the central array element in the sub-array strong, while the excitation signal amplitude of the array elements on both sides gradually weakens. Commonly used amplitude apodization functions include hanning function, hamming function, Blackman function and so on. Select the hanning window for amplitude traversal, the expression of the hanning window is:

x＝-halfA_x,...,halfA_x,y＝-halfA_y,...,halfA_y (20)x＝-halfA_x ,...,halfA_x ,y＝-halfA_y ,...,halfA_y (20)

(6)在已得到时间延时delays和横纵向孔径下标的基础上，从三维数据每条线(x,y)中寻找相应的值填充到预先初始化的数组chnls(x,y)中，从而得到幅度变迹处理后位置(X,Y,Z)的值：(6) On the basis of the obtained time delay delays and horizontal and vertical aperture subscripts, from Find the corresponding value in each line (x, y) of the three-dimensional data and fill it into the pre-initialized array chnls(x, y), so as to obtain the value of the position (X, Y, Z) after amplitude apodization processing:

(7)计算动态孔径接收处理过后此采样点的值(7) Calculate the value of this sampling point after the dynamic aperture receiving process

$((X x,, Y Y,, Z Z)) = = {Σ Σ}_{x x = = - - half half {A A}_{x x}}^{{halfA halfA}_{x x}} {Σ Σ}_{y the y = = - - {halfA halfA}_{y the y}}^{{halfA halfA}_{y the y}} chnls chnls ((x x,, y the y)) \times \times win win ((x x,, y the y)) - - - - - - ((21 twenty one))$

步骤四、为得到三维复合成像，需要设置面阵阵元时间延时，得到新的体数据，并按照不同的解决方案将新坐标系下的体数据转换到参考坐标下，多次平均得到用于三维重建的预处理数据，其具体步骤如下：Step 4. In order to obtain 3D composite imaging, it is necessary to set the time delay of the area array elements to obtain new volume data, and convert the volume data in the new coordinate system to the reference coordinates according to different solutions, and obtain the used volume data by averaging multiple times. The specific steps for preprocessing data for 3D reconstruction are as follows:

参阅图3所示，具体描述如下：Referring to Figure 3, the specific description is as follows:

以下假设参考坐标系为HIFU照射轴向方向的垂直方向，参考坐标系坐标为(x,y,z)，旋转坐标系坐标为(x',y',z')，参考坐标系的方向向量为(v_x,v_y,v_z)，旋转坐标系的方向向量为(v'_x,v'_y,v'_z)，平移矩阵用T表示，旋转矩阵用R表示，转换矩阵用M表示，面阵的横纵向孔径大小为A，被动接收平面为S，面阵阵元间隔为d。The following assumes that the reference coordinate system is the vertical direction of the HIFU irradiation axial direction, the coordinates of the reference coordinate system are (x, y, z), the coordinates of the rotating coordinate system are (x', y', z'), and the direction vector of the reference coordinate system is (v_x , v_y , v_z ), the direction vector of the rotating coordinate system is (v'_x , v'_y , v'_z ), the translation matrix is represented by T, the rotation matrix is represented by R, and the transformation matrix is represented by M , the horizontal and vertical aperture size of the area array is A, the passive receiving plane is S, and the element spacing of the area array is d.

(1)若被动接收平面S_x与XOY平面以y轴正方向(负方向)成α角，相应的时间延时为则旋转坐标系的原点和参考坐标系的原点之间存在平移量，平移矩阵为：(1) If the passive receiving plane S_x and the XOY plane form an angle α in the positive direction (negative direction) of the y-axis, the corresponding time delay is Then there is a translation amount between the origin of the rotating coordinate system and the origin of the reference coordinate system, and the translation matrix is:

${T T}_{x x} = = [\begin{matrix} 11 & 00 & 00 & 00 \\ 00 & 11 & 00 & 00 \\ 00 & 00 & 11 & 00 \\ 00 & 00 & {t t}_{x x} & 11 \end{matrix}] - - - - - - ((22 twenty two))$

其中以y轴正方向成α角，t_x＝d sinα；以y轴负方向成α角，t_x＝Ad sinα。Where the positive direction of the y-axis forms an angle α, t_x =d sin α; the negative direction of the y-axis forms an angle α, t_x =Ad sin α.

旋转点的x坐标不变，y,z的坐标相当于在YOZ平面内做正α旋转，旋转矩阵为：The x coordinate of the rotation point remains unchanged, and the coordinates of y and z are equivalent to positive α rotation in the YOZ plane. The rotation matrix is:

${R R}_{x x} = = [\begin{matrix} 11 & 00 & 00 & 00 \\ 00 & cos cos α α & sin sin α α & 00 \\ 00 & - - sin sin α α & cos cos α α & 00 \\ 00 & 00 & 00 & 11 \end{matrix}] - - - - - - ((23 twenty three))$

经过平移和旋转后可得坐标转换矩阵：After translation and rotation, the coordinate transformation matrix can be obtained:

M_x(α)＝T_x(α)R_x(α) (24)M_x (α) = T_x (α) R_x (α) (24)

将旋转坐标系下的数据对应到参考坐标系下的数据，求解下述方程可得：Correspond the data in the rotating coordinate system to the data in the reference coordinate system, and solve the following equation to get:

$[[x x,, y the y,, z z,, 11]] = = [[{x x}^{' '},, {y the y}^{' '},, {z z}^{' '},, 11]] {R R}_{x x}^{- - 11} ((α α)) {T T}_{x x}^{- - 11} ((α α)) - - - - - - ((2525))$

(2)若被动接收平面S_y与XOY平面以x轴正方向(负方向)成β角，相应的时间延时为则旋转坐标系的原点和参考坐标系的原点之间存在平移量，平移矩阵为：(2) If the passive receiving plane S_y and the XOY plane form an angle β with the positive direction (negative direction) of the x-axis, the corresponding time delay is Then there is a translation amount between the origin of the rotating coordinate system and the origin of the reference coordinate system, and the translation matrix is:

${T T}_{y the y} = = [\begin{matrix} 11 & 00 & 00 & 00 \\ 00 & 11 & 00 & 00 \\ 00 & 00 & 11 & 00 \\ 00 & 00 & {t t}_{y the y} & 11 \end{matrix}] - - - - - - ((2626))$

其中以x轴正方向成β角，t_y＝d sinβ；以x轴负方向成β角，t_y＝Ad sinβ。Where the positive direction of the x-axis forms an angle of β,_ty = d sin β; the negative direction of the x-axis forms an angle of β,_ty = Ad sin β.

旋转点的y坐标不变，x,z的坐标相当于在XOZ平面内做正β旋转，旋转矩阵为：The y coordinate of the rotation point remains unchanged, and the x and z coordinates are equivalent to positive β rotation in the XOZ plane. The rotation matrix is:

${R R}_{y the y} = = [\begin{matrix} cos cos β β & 00 & - - sin sin β β & 00 \\ 00 & 11 & 00 & 00 \\ sin sin β β & 00 & cos cos β β & 00 \\ 00 & 00 & 00 & 11 \end{matrix}] - - - - - - ((2727))$

M_y(β)＝T_y(β)R_y(β) (28)M_y (β)=T_y (β)R_y (β) (28)

将旋转坐标系下的数据对应到参考坐标系下的数据，求解方程可得：Correspond the data in the rotating coordinate system to the data in the reference coordinate system, and solve the equation to get:

$[[x x,, y the y,, z z,, 11]] = = [[{x x}^{' '},, {y the y}^{' '},, {z z}^{' '},, 11]] {R R}_{y the y}^{- - 11} ((β β)) {T T}_{y the y}^{- - 11} ((β β)) - - - - - - ((2929))$

(3)若被动接受平面绕任意轴旋转，在已知相对于参考坐标系的平移量之后，首先让旋转坐标系的原点经过平移后与参考坐标系的原点重合，平移矩阵为：(3) If the plane is passively accepted to rotate around any axis, after the translation relative to the reference coordinate system is known, first let the origin of the rotating coordinate system coincide with the origin of the reference coordinate system after being translated, and the translation matrix is:

${T T}_{arb arb} = = [\begin{matrix} 11 & 00 & 00 & 00 \\ 00 & 11 & 00 & 00 \\ 00 & 00 & 11 & 00 \\ {Δ Δ}_{x x} & {Δ Δ}_{y the y} & {Δ Δ}_{z z} & 11 \end{matrix}] - - - - - - ((3030))$

${R R}_{arb arb} = = [\begin{matrix} {v v}_{{x x}_{11}}^{' '} & {v v}_{{y the y}_{11}}^{' '} & {v v}_{{z z}_{11}}^{' '} & 00 \\ {v v}_{{x x}_{22}}^{' '} & {v v}_{{y the y}_{22}}^{' '} & {v v}_{{z z}_{22}}^{' '} & 00 \\ {v v}_{{x x}_{33}}^{' '} & {v v}_{{y the y}_{33}}^{' '} & {v v}_{{z z}_{33}}^{' '} & 00 \\ 00 & 00 & 00 & 11 \end{matrix}] - - - - - - ((3131))$

(v'_xi,v'_yi,v'_zi)(i＝1～3)为旋转坐标系方向向量(v'_x,v'_y,v'_z)分解到参考坐标系上所得到的分解向量；(v'_xi , v'_yi , v'_zi ) (i=1～3) is the decomposition vector obtained by decomposing the direction vector (v'_x , v'_y , v'_z ) of the rotating coordinate system into the reference coordinate system ;

旋转坐标系经过平移和旋转后得坐标转换矩阵：After the rotating coordinate system is translated and rotated, the coordinate transformation matrix is obtained:

M_arb＝T_arb(Δ_x,Δ_y,Δ_z)R_arb(v'_x,v'_y,v'_z) (32)M_arb ＝T_arb (Δ_x ,Δ_y ,Δ_z )R_arb (v'_x ,v'_y ,v'_z ) (32)

将旋转坐标系下的数据对应到参考坐标系下的数据，求解下面方程得用于三维重建的预处理三维体数据E_pre3D(x,y,z)：Correspond the data in the rotating coordinate system to the data in the reference coordinate system, and solve the following equation to get the preprocessed 3D volume data E_pre3D (x,y,z) for 3D reconstruction:

${E E.}_{pre pre 33 D D.} ((x x,, y the y,, z z)) = = [[{x x}^{' '},, {y the y}^{' '},, {z z}^{' '} 11]] {R R}_{arb arb}^{- - 11} (({v v}_{x x}^{' '},, {v v}_{y the y}^{' '} {v v}_{z z}^{' '})) {T T}_{arb arb}^{- - 11} (({Δ Δ}_{x x},, {Δ Δ}_{y the y},, {Δ Δ}_{z z})) - - - - - - ((3333))$

(4)用来表示第q个角度下得到的预处理三维体数据，重复多次设置接收阵元时间延时、采集数据并对采集数据做被动成像算法、动态孔径接收和幅度变迹技术的处理，得到多个角度下的预处理三维体数据，然后对多个角度下的预处理三维体数据进行平均，则可得到多次平均后的预处理三维体数据：(4) use to represent the preprocessed 3D volume data obtained at the qth angle, repeatedly setting the time delay of receiving array elements, collecting data, and processing the collected data with passive imaging algorithm, dynamic aperture receiving and amplitude apodization technology, and multiple The preprocessed 3D volume data under one angle, and then average the preprocessed 3D volume data under multiple angles, then the preprocessed 3D volume data after multiple averages can be obtained:

${\overset{&OverBar; &OverBar;}{E E.}}_{pre pre 33 D D.}^{Q Q} ((x x,, y the y,, z z)) = = \frac{11}{Q Q} {Σ Σ}_{q q = = 11}^{Q Q} {E E.}_{pre pre 33 D D.}^{((q q))} ((x x,, y the y,, z z)) - - - - - - ((3434))$

(5)由于大量的数据存储却难以避免，因此为了得到第Q+1次平均后的预处理三维体数据以对三维图像实时更新。采用下面公式(5) Since a large amount of data storage is unavoidable, the 3D image is updated in real time in order to obtain the preprocessed 3D volume data after the Q+1th average. Use the following formula

${\overset{&OverBar; &OverBar;}{E E.}}_{pre pre 33 D D.}^{Q Q + + 11} ((x x,, y the y,, z z)) = = \frac{Q Q}{Q Q + + 11} {\overset{&OverBar; &OverBar;}{E E.}}_{pre pre 33 D D.}^{Q Q} ((x x,, y the y,, z z)) + + \frac{11}{Q Q + + 11} {E E.}_{pre pre 33 D D.}^{((Q Q + + 11))} ((x x,, y the y,, z z)) - - - - - - ((3535))$

其中为Q次预处理三维体数据的平均，为第Q+1次的预处理三维体数据。in is the average of Q preprocessed 3D volume data, It is the preprocessed three-dimensional volume data of the Q+1 time.

步骤五、在步骤四得到三维重建前预处理的三维体数据的基础上分割子区域并计算B样条插值点集，其具体步骤如下：Step five, on the basis of the preprocessed three-dimensional volume data obtained in step four before three-dimensional reconstruction, segment the sub-regions and calculate the B-spline interpolation point set, the specific steps are as follows:

参阅图6所示，具体描述如下：Referring to Figure 6, the specific description is as follows:

(1)首先分割的子区域需要尽可能满足以下条件：(1) The sub-regions to be divided first need to meet the following conditions as much as possible:

1)各PCI子区域内的数据点个数差异不能太大；1) The difference in the number of data points in each PCI sub-area should not be too large;

2)各PCI子区域内的数据点投影对应的函数值差异不能太大；2) The difference in the function value corresponding to the data point projection in each PCI sub-area should not be too large;

3)各PCI子区域的边界上尽量保留原始数据点；3) Try to keep the original data points on the boundary of each PCI sub-area;

(2)假设体数据点子区域C_k被数据盒D_k所包围，|C_k|为体数据点子区域点集内的点数，假设S_max和s为子区域可允许的最大数据点数和最小尺寸；(2) Assume that the volume data point sub-area C_k is surrounded by the data box D_k , |C_k | is the number of points in the volume data point sub-area point set, assuming S_max and s are the maximum number of data points and the minimum size allowed by the sub-area ;

(3)如果|C_k|>S_max而且沿着x轴方向的尺寸大于最小边长，则计算C_k内所有点的重心，记为(x_w,y_w,z_w)，然后用平面x＝x_w分割点集，如此往复循环直到子区域内的数据点数|C_k|和尺寸满足“|C_k|>S_max而且沿着x轴方向的尺寸大于最小边长”；(3) If |C_k |>S_max and the size along the x-axis direction is greater than the minimum side length, then calculate the center of gravity of all points in C_k , denoted as (x_w , y_w , z_w ), and then use the plane x=x_w divides the point set, and so on and so forth until the number of data points in the sub-region |C_k | and the size meet "|C_k |>S_max and the size along the x-axis direction is greater than the minimum side length";

(4)然后把相邻的两个数据盒做并集运算，使两个数据盒相交，并减去公共的部分，假如相并之后的数据盒不满足“|C_k|>S_max而且沿着y轴方向的尺寸大于最小边长”这一条件，则用平面y＝y_w分割点集直到满足这一条件；(4) Then perform a union operation on two adjacent data boxes, make the two data boxes intersect, and subtract the common part, if the combined data box does not satisfy "|C_k |>S_max and along If the size in the y-axis direction is greater than the minimum side length", then use the plane y=y_w to divide the point set until this condition is met;

(5)接着继续把相邻的两个数据盒做并集运算，使两个数据盒相交，并减去公共的部分，假如相并之后的数据盒不满足“|C_k|>S_max而且沿着z轴方向的尺寸大于最小边长”这一条件，则用平面z＝z_w分割点集直到满足这一条件；(5) Then continue to perform the union operation on the two adjacent data boxes, make the two data boxes intersect, and subtract the common part, if the combined data box does not satisfy "|C_k |>S_max and The size along the z-axis direction is greater than the minimum side length", then use the plane z=z_w to divide the point set until this condition is satisfied;

(6)常用的径向基函数主要有：(6) The commonly used radial basis functions mainly include:

线性函数：φ(r)＝rLinear function: φ(r)=r

三次函数：φ(r)＝r³Cubic function: φ(r)=r³

Kriging方法的Gauss分布函数：Gauss distribution function of Kriging method:

Kriging方法的Markov分布函数：φ(r)＝e^-ar和其他概率分布函数Markov distribution function of Kriging method: φ(r)=e^-ar and other probability distribution functions

Hardy的multiquadric函数：φ(r)＝(c²+r²)^β以及逆multiquadric函数：φ(r)＝(c²+r²)^-βHardy's multiquadric function: φ(r)=(c² +r² )^β and inverse multiquadric function: φ(r)=(c² +r² )^-β

Duchon的薄板样条：其中d为空间维数Duchon's thin plate spline: where d is the dimension of the space

紧支柱正定径向基函数：这类函数Φ(x)＝φ(||x||)主要有如下形式的多项式：Positive definite radial basis functions of tight struts: This type of function Φ(x)=φ(||x||) mainly has polynomials of the following form:

$φ (r) = \{\begin{matrix} p (r), & 0 \leq r \leq 1 \\ 0, & r > 1 \end{matrix},$ 其中 $p (r) = Σ_{j = 0}^{N} c_{j} r^{j},$ N为次数。 $φ (r) = \{\begin{matrix} p (r), & 0 \leq r \leq 1 \\ 0, & r > 1 \end{matrix},$ in $p (r) = Σ_{j = 0}^{N} c_{j} r^{j},$ N is the number of times.

(7)如果将所有插值函数直接求和，会导致三维重建的错误。因此设置窗函数(7) If all interpolation functions are directly summed, it will lead to errors in 3D reconstruction. Therefore, setting the window function

${ω ω}_{k k} ((r r)) = = ((11 - - \frac{| | 22 x x - - (({x x}_{k k max max} + + {x x}_{k k min min})) | |}{{x x}_{k k max max} - - {x x}_{k k min min}})) ((11 - - \frac{| | 22 y the y - - (({y the y}_{k k max max} + + {y the y}_{k k min min})) | |}{{y the y}_{k k max max} - - {y the y}_{k k min min}})) ((11 - - \frac{| | 22 z z (({z z}_{k k max max} + + {z z}_{k k min min})) | |}{{z z}_{k k max max} - - {z z}_{k k min min}})) - - - - - - ((3838))$

(8)假设第k个子区域数据点集C_k＝{r_i(i＝1,2,...,|C_k|)}，所对应的函数值为f(X_i)(i＝1,2,...,|C_k|)，则可得方程组(8) Suppose the kth sub-region data point set C_k ={r_i (i=1,2,...,|C_k |)}, the corresponding function value is f(X_i )(i=1 ,2,...,|C_k |), then the system of equations can be obtained

$[\begin{matrix} {a a}_{11} & {a a}_{22} & . . . . . . & {a a}_{| | {C C}_{k k} | |} \end{matrix}] [\begin{matrix} φ φ ((| | | | {X x}_{i i} - - {X x}_{11} | | | |)) \\ φ φ ((| | | | {X x}_{i i} - - {X x}_{11} | | | |)) \\ . . . . . . \\ φ φ ((| | | | {X x}_{i i} - - {X x}_{j j} | | | |)) \end{matrix}] = = f f (({X x}_{i i})) ((i i = = 1,2 1,2,, . . . . . .,, | | {C C}_{k k} | |)) - - - - - - ((3939))$

(9)解方程组得到的值，求解每一个子区域的径向插值方程并计算网格点对应的函数值即B样条插值的点集B_ij(i＝1,2,...,M,j＝1,2,...,N)。(9) Solve the system of equations to get value, solve the radial interpolation equation of each sub-region and calculate the function value corresponding to the grid point, which is the point set B_ij of B-spline interpolation (i=1,2,...,M,j=1,2 ,...,N).

(10)B样条插值三维重建仿真实例参阅图4所示。(10) Refer to Figure 4 for an example of B-spline interpolation 3D reconstruction simulation.

步骤六、结合径向基函数插值和B样条插值对三维体数据进行三维重建，并对三维图像实时更新。其具体步骤如下：Step 6: Combining radial basis function interpolation and B-spline interpolation to perform 3D reconstruction on 3D volume data, and update the 3D image in real time. The specific steps are as follows:

(1)假设经过处理得到的预处理三维体数据点集合为C，设C₀为PCI点集合C在XOY平面上的投影点集，并划定出C₀的边界；(1) Assuming that the preprocessed three-dimensional volume data point set obtained through processing is C, let C be the projection point set of PCI point set C on the XOY plane, and delineate the boundary of C₀_;

(2)将投影点集C₀分割成M×N个规则的小矩形网格，相邻矩形的交点则对应着B样条插值的点在XOY平面上的投影；此处N与前文采样点数N无关。(2) Divide the projected point set C₀ into M×N regular small rectangular grids, and the intersection points of adjacent rectangles correspond to the projections of the B-spline interpolation points on the XOY plane; here N is the same as the number of sampling points above N is irrelevant.

(3)将集合C划分成L个子区域，每一块子区域C_k数据点的个数为|C_k|(k＝1,2,...,L)个；(3) Divide the set C into L sub-regions, and the number of data points in each sub-region C_k is |C_k |(k=1,2,...,L);

(4)设置不同的径向基函数rbf_k(k＝1,2,...,L)对每一个子区域C_k(k＝1,2,...,L)内的数据点分别插值，得到插值后的三维子区域；(4) Set different radial basis functions rbf_k (k=1,2,...,L) for the data points in each sub-region C_k (k=1,2,...,L) Interpolation to obtain the interpolated three-dimensional sub-region;

(5)根据每个子区域插值函数rbf_k(k＝1,2,...,L)分别求出C₀划分后网格点对应的函数值，组成了B样条插值的点集B_ij(i＝1,2,...,M,j＝1,2,...,N)；(5) According to the interpolation function rbf_k (k=1,2,...,L) of each sub-region, calculate the function value corresponding to the grid point after C₀ is divided, and form the point set B_ij of B-spline interpolation (i=1,2,...,M,j=1,2,...,N);

(6)对点集B_ij做B样条插值运算，生成B样条插值后的三维成像。(6) Perform B-spline interpolation operation on point set B_ij to generate 3D imaging after B-spline interpolation.

(7)径向基函数插值三维重建仿真实例参阅图5所示。(7) Refer to Figure 5 for an example of radial basis function interpolation 3D reconstruction simulation.

被动空化成像(PCI)方法第一次与基因治疗和药物传送等结合并取得了较好的结果这一事实说明PCI技术在监控治疗领域拥有极广的前景。此方法中所涉多角度复合成像和实时更新方法更可以满足临床上需求，将三维空化被动成像技术与HIFU测温结合起来可以用来实时监控HIFU治疗过程中组织消融的程度和预测HIFU焦前、焦点和焦后产生的损伤大小。还有，组织消融中存在稳态过程，因此可在步骤四中多次设置复合角度以实现更加精确的监控成像。在药物释放和传送、基因治疗、经颅超声溶栓和血脑屏障破坏等方面，此方法可以提供非侵入的、更有效的、可测量的定位和监控手段。The fact that the passive cavitation imaging (PCI) method has been combined with gene therapy and drug delivery for the first time and achieved good results shows that PCI technology has a very broad prospect in the field of monitoring and treatment. The multi-angle composite imaging and real-time update method involved in this method can better meet the clinical needs. The combination of three-dimensional cavitation passive imaging technology and HIFU temperature measurement can be used to monitor the degree of tissue ablation during HIFU treatment in real time and predict the focus of HIFU. The size of the damage produced before, after focus and after focus. In addition, there is a steady-state process in tissue ablation, so the composite angle can be set multiple times in step 4 to achieve more accurate monitoring imaging. This method can provide non-invasive, more effective, measurable localization and monitoring means in drug release and delivery, gene therapy, transcranial ultrasound thrombolysis and blood-brain barrier disruption.