CN112823283A - 通过使用超声非侵入性地表征非均匀介质的方法和系统 - Google Patents
通过使用超声非侵入性地表征非均匀介质的方法和系统 Download PDFInfo
- Publication number
- CN112823283A CN112823283A CN201980061452.1A CN201980061452A CN112823283A CN 112823283 A CN112823283 A CN 112823283A CN 201980061452 A CN201980061452 A CN 201980061452A CN 112823283 A CN112823283 A CN 112823283A
- Authority
- CN
- China
- Prior art keywords
- matrix
- basis
- distortion
- reflection
- determining
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000002604 ultrasonography Methods 0.000 title claims abstract description 99
- 238000000034 method Methods 0.000 title claims description 76
- 239000011159 matrix material Substances 0.000 claims abstract description 398
- 230000004075 alteration Effects 0.000 claims abstract description 122
- 238000013507 mapping Methods 0.000 claims abstract description 30
- 229940050561 matrix product Drugs 0.000 claims abstract description 16
- 230000009466 transformation Effects 0.000 claims abstract description 12
- 230000036962 time dependent Effects 0.000 claims abstract description 4
- 238000012937 correction Methods 0.000 claims description 69
- 230000003595 spectral effect Effects 0.000 claims description 24
- 238000012512 characterization method Methods 0.000 claims description 21
- 238000000354 decomposition reaction Methods 0.000 claims description 12
- 230000001427 coherent effect Effects 0.000 claims description 4
- 239000013598 vector Substances 0.000 description 31
- 238000009792 diffusion process Methods 0.000 description 22
- 239000000523 sample Substances 0.000 description 19
- 239000010410 layer Substances 0.000 description 16
- 230000005540 biological transmission Effects 0.000 description 13
- 229920005372 Plexiglas® Polymers 0.000 description 12
- 238000001914 filtration Methods 0.000 description 12
- VVQNEPGJFQJSBK-UHFFFAOYSA-N Methyl methacrylate Chemical compound COC(=O)C(C)=C VVQNEPGJFQJSBK-UHFFFAOYSA-N 0.000 description 9
- 238000003384 imaging method Methods 0.000 description 9
- 238000005259 measurement Methods 0.000 description 9
- 238000012285 ultrasound imaging Methods 0.000 description 8
- 230000008901 benefit Effects 0.000 description 7
- 230000004044 response Effects 0.000 description 7
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 7
- 238000013459 approach Methods 0.000 description 6
- 238000011835 investigation Methods 0.000 description 6
- 238000002310 reflectometry Methods 0.000 description 6
- 210000001519 tissue Anatomy 0.000 description 6
- 238000003491 array Methods 0.000 description 5
- 238000004364 calculation method Methods 0.000 description 5
- 238000010586 diagram Methods 0.000 description 5
- 238000002592 echocardiography Methods 0.000 description 5
- 238000002474 experimental method Methods 0.000 description 5
- 125000001475 halogen functional group Chemical group 0.000 description 5
- 230000008569 process Effects 0.000 description 5
- 238000004458 analytical method Methods 0.000 description 4
- 230000001934 delay Effects 0.000 description 4
- 238000012545 processing Methods 0.000 description 4
- 238000013519 translation Methods 0.000 description 4
- 238000011144 upstream manufacturing Methods 0.000 description 4
- 208000004434 Calcinosis Diseases 0.000 description 3
- 230000003044 adaptive effect Effects 0.000 description 3
- 238000002059 diagnostic imaging Methods 0.000 description 3
- 239000006185 dispersion Substances 0.000 description 3
- 210000004185 liver Anatomy 0.000 description 3
- 238000010606 normalization Methods 0.000 description 3
- 239000004926 polymethyl methacrylate Substances 0.000 description 3
- 210000001747 pupil Anatomy 0.000 description 3
- 238000011002 quantification Methods 0.000 description 3
- 230000009467 reduction Effects 0.000 description 3
- 210000000481 breast Anatomy 0.000 description 2
- 230000008859 change Effects 0.000 description 2
- 238000006243 chemical reaction Methods 0.000 description 2
- 230000004069 differentiation Effects 0.000 description 2
- 230000006872 improvement Effects 0.000 description 2
- 238000001727 in vivo Methods 0.000 description 2
- 238000012804 iterative process Methods 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 210000003205 muscle Anatomy 0.000 description 2
- 230000003287 optical effect Effects 0.000 description 2
- 210000000056 organ Anatomy 0.000 description 2
- 210000004872 soft tissue Anatomy 0.000 description 2
- 239000002344 surface layer Substances 0.000 description 2
- 238000003325 tomography Methods 0.000 description 2
- 230000001131 transforming effect Effects 0.000 description 2
- 238000002834 transmittance Methods 0.000 description 2
- 238000009825 accumulation Methods 0.000 description 1
- 210000000577 adipose tissue Anatomy 0.000 description 1
- 238000004422 calculation algorithm Methods 0.000 description 1
- 230000015556 catabolic process Effects 0.000 description 1
- 238000004590 computer program Methods 0.000 description 1
- 230000021615 conjugation Effects 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 230000007423 decrease Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000006731 degradation reaction Methods 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 229940079593 drug Drugs 0.000 description 1
- 239000003814 drug Substances 0.000 description 1
- 230000009977 dual effect Effects 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 230000005284 excitation Effects 0.000 description 1
- 230000005621 ferroelectricity Effects 0.000 description 1
- 230000010006 flight Effects 0.000 description 1
- 238000003702 image correction Methods 0.000 description 1
- 230000033001 locomotion Effects 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 238000009659 non-destructive testing Methods 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 230000010363 phase shift Effects 0.000 description 1
- 238000012805 post-processing Methods 0.000 description 1
- 238000000513 principal component analysis Methods 0.000 description 1
- 230000001902 propagating effect Effects 0.000 description 1
- 238000004451 qualitative analysis Methods 0.000 description 1
- 238000012552 review Methods 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 238000001228 spectrum Methods 0.000 description 1
- 238000010183 spectrum analysis Methods 0.000 description 1
- 230000007704 transition Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N29/00—Investigating or analysing materials by the use of ultrasonic, sonic or infrasonic waves; Visualisation of the interior of objects by transmitting ultrasonic or sonic waves through the object
- G01N29/44—Processing the detected response signal, e.g. electronic circuits specially adapted therefor
- G01N29/46—Processing the detected response signal, e.g. electronic circuits specially adapted therefor by spectral analysis, e.g. Fourier analysis or wavelet analysis
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N29/00—Investigating or analysing materials by the use of ultrasonic, sonic or infrasonic waves; Visualisation of the interior of objects by transmitting ultrasonic or sonic waves through the object
- G01N29/04—Analysing solids
- G01N29/06—Visualisation of the interior, e.g. acoustic microscopy
- G01N29/0654—Imaging
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N29/00—Investigating or analysing materials by the use of ultrasonic, sonic or infrasonic waves; Visualisation of the interior of objects by transmitting ultrasonic or sonic waves through the object
- G01N29/04—Analysing solids
- G01N29/06—Visualisation of the interior, e.g. acoustic microscopy
- G01N29/0654—Imaging
- G01N29/0672—Imaging by acoustic tomography
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N29/00—Investigating or analysing materials by the use of ultrasonic, sonic or infrasonic waves; Visualisation of the interior of objects by transmitting ultrasonic or sonic waves through the object
- G01N29/44—Processing the detected response signal, e.g. electronic circuits specially adapted therefor
- G01N29/4463—Signal correction, e.g. distance amplitude correction [DAC], distance gain size [DGS], noise filtering
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S15/00—Systems using the reflection or reradiation of acoustic waves, e.g. sonar systems
- G01S15/88—Sonar systems specially adapted for specific applications
- G01S15/89—Sonar systems specially adapted for specific applications for mapping or imaging
- G01S15/8906—Short-range imaging systems; Acoustic microscope systems using pulse-echo techniques
- G01S15/8909—Short-range imaging systems; Acoustic microscope systems using pulse-echo techniques using a static transducer configuration
- G01S15/8915—Short-range imaging systems; Acoustic microscope systems using pulse-echo techniques using a static transducer configuration using a transducer array
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S15/00—Systems using the reflection or reradiation of acoustic waves, e.g. sonar systems
- G01S15/88—Sonar systems specially adapted for specific applications
- G01S15/89—Sonar systems specially adapted for specific applications for mapping or imaging
- G01S15/8906—Short-range imaging systems; Acoustic microscope systems using pulse-echo techniques
- G01S15/8977—Short-range imaging systems; Acoustic microscope systems using pulse-echo techniques using special techniques for image reconstruction, e.g. FFT, geometrical transformations, spatial deconvolution, time deconvolution
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
- G01S7/52—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S15/00
- G01S7/52017—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S15/00 particularly adapted to short-range imaging
- G01S7/52046—Techniques for image enhancement involving transmitter or receiver
- G01S7/52049—Techniques for image enhancement involving transmitter or receiver using correction of medium-induced phase aberration
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N2291/00—Indexing codes associated with group G01N29/00
- G01N2291/02—Indexing codes associated with the analysed material
- G01N2291/024—Mixtures
- G01N2291/02475—Tissue characterisation
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N2291/00—Indexing codes associated with group G01N29/00
- G01N2291/10—Number of transducers
- G01N2291/106—Number of transducers one or more transducer arrays
Landscapes
- Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Immunology (AREA)
- Pathology (AREA)
- Chemical & Material Sciences (AREA)
- Analytical Chemistry (AREA)
- Biochemistry (AREA)
- General Health & Medical Sciences (AREA)
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Acoustics & Sound (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Signal Processing (AREA)
- Computer Networks & Wireless Communication (AREA)
- Mathematical Physics (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Ultra Sonic Daignosis Equipment (AREA)
- Investigating Or Analyzing Materials By The Use Of Ultrasonic Waves (AREA)
Abstract
在第一方面，本说明书涉及一种通过使用超声非侵入性地表征非均匀介质的系统，该系统包括：至少一个第一换能器阵列(10)，其被配置为在所述非均匀介质的区域内产生一系列入射超声波并记录由所述区域反向散射的超声波随时间的变化；以及计算单元(42)，该计算单元与所述至少一个第一换能器阵列耦合并被配置为：记录在输入处的发射基和输出处的换能器基之间定义的实验性反射矩阵；从所述实验性反射矩阵确定在输入和输出处以聚焦基定义的至少一个第一宽带反射矩阵；确定在所述聚焦基和观察基之间定义的第一畸变矩阵，所述第一畸变矩阵是从在所述聚焦基和用于校正像差的基之间定义的所述第一宽带反射矩阵与针对模型介质在与上述相同的基之间定义的参考反射矩阵的相位共轭矩阵的逐项矩阵乘积直接获得或在基变换之后获得；从所述第一畸变矩阵确定所述非均匀介质的物理参数的至少一个映射。
Description
技术领域
本说明书涉及对非均匀介质进行非侵入性超声表征的方法和系统，尤其适用于医学成像或非破坏性测试，并且更普遍而言，适用于可以使用超声成像的所有领域。
背景技术
在声成像领域，试图用超声波主动探测未知介质来对其进行表征。这尤其是医学成像中回波描记术的原理。
声成像系统的分辨率可以被定义为辨别对象的微小细节的能力。原则上，声成像系统受衍射限制，并且其理论分辨率由λ/2(其中λ是声在介质中的波长)或检测器的有限角孔径给出。然而，实际上，分辨率通常由于声速在传播介质中的变化而降低。
具体而言，通常在声成像中，介质被认为是均匀的，其具有恒定的声速c0。然而，并不总是能够满足均匀介质的假设。例如，在对肝脏进行回波描记术的情况下，探头被放置在患者的肋骨之间。声波在到达目标器官之前先经过一系列脂肪层和肌肉层。软组织具有不同的机械性能。因此，声速远非均匀，而是位于针对脂肪组织的1450m/s和针对肝脏的1600m/s之间。声速的变化会导致波根据其传播经过的位置而有不同的相移。这导致声波波前的像差，该像差导致得到的超声图像产生畸变，并因此导致该超声图像的分辨率和对比度下降。像差可能会损及医疗检查。
在过去的几十年中，已经开发出各种方法来校正像差(这主要发生在医学成像领域)，以便在超声图像的分辨率由于介质的非均匀性而降低的情况下获得具有最佳分辨率的超声图像。然而，在复数像差的情况下，波前畸变会根据图像的区域而变化。符合等晕域条件的这些区域通常无法在理论上确定，因此构成了应用现有技术的基本限制。
如图1A至图1C所示，传统的回波描记方法使用压电换能器11的阵列10，这些压电换能器11可以独立地发射和接收超声脉冲。每个换能器的位置由向量u标识。当将这样的阵列面对要研究的介质20放置时，可以以各种方式对介质20进行声穿透。
对要研究的介质进行声穿透的第一种方法是从阵列中的一个换能器发射超声脉冲，该换能器的位置通过uin标识(图1A，左图)。对于一维(或二维)换能器阵列，这会导致发散的圆柱状(或球状)入射波。该入射波由介质20中的散射体21反射，并且每个换能器11将记录反向散射场随时间的变化(图1A，右图)。通过连续使用每个换能器作为波源来重复此操作，可以测量每个换能器之间的所有脉冲响应R(uout,uin,t)，其中uout表示检测器的位置。这些响应形成以换能器基表示的反射矩阵Ruu(t)。这种测量的优点在于，该矩阵包含有关所研究介质的所有信息，然后例如能够对其应用一组矩阵运算，以对介质进行成像。然而，这种获取方式要假定介质在整个测量过程中保持固定，这在体内实验的情况下可能是非常有限的。另外，单个压电元件发射的能量较低，这可能导致较差的信噪比。
已知还有其他方法可以对要研究的介质进行声穿透，这些技术称为“波束成形”，指的是通过应用延迟线对波进行相干控制。在这些技术中，已知有聚焦发射，这是常规的回波描记术广泛使用的一项技术。如图1B(左图)所示，仅需对换能器11施加适当的一组延迟以校正行进时间，从而使所有脉冲一起到达目标焦点rin。由于衍射的物理限制，超声集中在由超声探头的孔所界定的区域中。注意，在回波描记术中，也要在接收时进行聚焦。然后，如图1B的右图所示，对由阵列10中的元件11捕获的所有回波进行处理，以对接收透镜的效果进行模拟。换能器接收到的信号通过在时间上进行移位而恢复同相。这些延迟与向发射施加的延迟相同。在发射阶段，所有信号都在点rin发生干涉。在接收时，来自同一点rout＝rin的信号通过在弹道时间t＝2||uout-rin||/c的信号加和而发生电子干涉。该加和给出了在接收时聚焦的最终结果。图1B中所示的方法(被称为共聚焦方法，具有发射和接收的双重聚焦)可以对介质的反射直接成像，其具有受衍射限制的横向分辨率、仅受初始脉冲持续时间限制的出色的轴向分辨率以及出色的对比度。然而，该方法很耗时，因为从任何理论上讲，该方法都必须聚焦发射在介质中的每个点。
最近开发的另一种波束成形技术包括用一系列平面波对介质进行声穿透。图1C示出了这种“平面波”回波描记法的原理，该原理例如在G.Montaldo等人的文章《Coherentplane-wave compounding for very high frame rate ultrasonography and transientelastography》(IEEE Trans.Ultrason.,Ferroelect.Freq.Control 56 489-506,2009)中进行了描述。在发射时向每个信号施加延迟(图1C，左图)，以形成相对于换能器阵列10以角度θin倾斜的波前。接收时(图1C，右图)，所有传感器uout针对一系列入射平面波测量由介质反向散射的场R(uout,θin,t)，该入射平面波的入射角为变化的θin。所有这些响应形成在输入处的傅立叶基(或平面波基)和输出处的换能器基之间定义的反射矩阵Ruθ(t)。一旦记录了该矩阵，就在信号以相干方式加和前对其进行时间移位，以针对超声图像的每个点rin使发射和接收的数据在数值上集中。与标准(聚焦发射)回波描记术相比，该方法形成超声图像所需的获取次数更少，对于为了获得相同水平的超声图像的对比度和分辨率来说，获取次数也更少。
受1950年代天文学中使用的光学像差校正技术的启发，为了校正超声成像中的相位畸变，引入了自适应聚焦技术(例如，参见D.J.Philips,S.W.Smith,O.T.Von Ramm andF.L.Thurstone,"A phase compensation technique for B-Modeechoencephalography".Ultrasound in Medicine,1,345-404,1975)。
如图2所示，该校正技术适用于所有常规的超声成像方法(图1A至图1C)。该技术假设存在薄的畸变层22，其在换能器11上引起简单的时间偏移。因此，可以通过在波束成形中通常使用的延迟中添加其他延迟定律来补偿这些时间偏移。如图2所示，该校正技术包括在不考虑像差的情况下介质20的激发(25)、像差定律的估计(26–27)以及在考虑畸变器的情况下新的发射/接收(28–29)。通过使用发射器和检测器之间的反馈回路，可以实时迭代系统，以找到用于校正像差的最佳时移。已经提出了几种方法，用于基于超声散斑(也就是说，由欠分辨散射体的随机分布进行的超声反射所产生的信号)来计算要施加到阵列中的每个换能器的延迟。一种方法是基于相邻换能器接收的信号之间的相关性(例如，参见M.O'Donnell and S.W.Flax，"Phase aberration measurements in medical ultrasound:Human studies",Ultrasonic Imaging 10,1-11,1988)。另一种方法是基于超声散斑的对比度优化((Nock,G.E.Trahey and S.W.Smith"Phase aberration correction inmedical ultrasound using speckle brightness as a quality factor",J.Acoust.Soc.Am,Vol.85,5,pp 1819-1833,1989)。然而，这些校正技术是基于薄的畸变层22的假设，该薄的畸变层22位于尽可能靠近换能器阵列10的位置。因此，对于体积性畸变器，这些方法仅在单个等晕域上有效。
等晕域的概念如图3A和图3B所示。在换能器阵列10与焦平面之间存在畸变器22。在图3A的示例中，向换能器阵列10的换能器11施加延迟，使得可以校正在畸变器22的输出处的声波前部(圆柱形波)。这导致在焦平面的给定点rin处产生衍射受限的焦点。如果使该经过相同校正的波前稍微变换角度(图3A)(或平移)以允许聚焦在│rin–r'in│<lc的点r'in上(其中lc是畸变器的参数特征，称为畸变器的相干长度)，那么来自畸变器的波前保持被校正，并且在r'in处的聚焦始终衍射受限。rin和r'in点属于声穿透区中的相同等晕域。然而，如果使校正后的入射波前更进一步变换角度(图3B)以允许聚焦在│rin–r'in│>lc的点r'in上，则来自畸变器的波前会畸变，并且所施加的校正不能使r'in处进行衍射受限的聚焦。在这种情况下，点rin和r'in属于声穿透区域的不同等晕域。
除了校正场限于单个等晕域外，由于波场的每个频谱分量可能会经历不同的相位畸变，上述像差校正技术通常不适用于超宽带信号(例如用于超声成像的信号)。为了对超声成像中的像差进行超宽带校正，已经开发出更复杂的方法，例如从时间反转概念衍生的逆滤波或匹配滤波方法。这些方法包括通过合适的时空滤波器对发射和接收的信号进行卷积，该时空滤波器使得可以优化波在介质中的聚焦并校正其像差。然而，要了解此过滤器，要求存在强大的散射体(C.Prada and M.Fink,"Eigenmodes of the time reversaloperator",Wave Motion,20,151-163,1994)或位于要成像的介质内的尤其具有侵入性的检测器阵列(例如参见M.Tanter,J.-L.Thomas and M.Fink."Time reversal and theinverse filter",J.Acoust.Soc.Am.108,223-234,2000)。
为了迭代地校正超声散斑中的像差，一种最新的方法直接受益于反射矩阵及其相关矩阵。因此，在J-L.Robert等人的文章("Green's function estimation in speckleusing the decomposition of the time reversal operator:Application toaberration correction in medical imaging",J.Acoust.Soc.Am.,123,2(2008))中，测量在输入处的聚焦发射基与输出处的换能器基之间的反射矩阵Rur。可以将聚焦在介质的点rin处的发射集合看作放置在介质的该点内的虚拟换能器阵列，它们各自的大小由每个焦点给出。因此，该反射矩阵包含在所研究介质内的虚拟阵列(rin)中的每个换能器与实际阵列中的每个换能器(uout)之间的脉冲响应R(uout,rin,t)的集合。通过将虚拟阵列的大小限制为一个等晕域，已显示出相关矩阵
(其中
是Rur的转共轭矩阵(或共轭矩阵的转置矩阵))的第一特征向量给出了使得可以聚焦在虚拟散射体(假星)上的波前，该虚拟散射体的大小由初始焦点给出。通过使该波前变换角度，可以聚焦在该等晕域的每个点上。然后可以使用该新的聚焦发射基对该过程进行重新迭代，重建新的反射矩阵Rur，并获得更精确的像差校正定律等。这种方法的缺点之一是需要重新迭代该过程，这会大大降低帧率。另外，像现有技术的其他技术一样，图像校正场再次限于单个等晕域。
在本说明书的上下文中，提出了一种原创方法，该方法既不基于假星的产生，也不基于迭代方法。所提出的方法还使得能够同时校正散斑图像的各个等晕域中的像差，还可以包括更多会引起入射波的镜面反射的确定性对象(界面、大散射体等)。
发明内容
根据第一方面，本说明书涉及一种用于对非均匀介质进行非侵入性超声表征的方法，该方法包括：
-通过至少一个第一换能器阵列在所述非均匀介质的区域中产生一系列入射超声波的步骤；
-记录在输入处的发射基和输出处的换能器基之间定义的实验性反射矩阵的步骤；
-基于所述实验性反射矩阵确定至少一个第一宽带反射矩阵的步骤，所述第一宽带反射矩阵在输入处和输出处以聚焦基在第一谱带中定义；
-确定在所述聚焦基和观察基之间在所述第一谱带中定义的第一畸变矩阵的步骤，所述第一畸变矩阵是由在所述聚焦基和像差校正基之间定义的所述第一宽带反射矩阵与针对模型介质在与上述相同的基之间定义的参考反射矩阵的相位共轭矩阵的逐项矩阵乘积直接获得或在基变换之后获得；
-基于所述第一畸变矩阵确定所述非均匀介质中的所述声穿透区的物理参数的至少一个映射的步骤。
如将在说明书的其余部分更详细地描述的，所述非均匀介质的物理参数可以包括例如与所述声穿透区中包含的至少一个等晕域相关联的像差定律、表征所述介质的局部反射的参数、表征在所述介质中的声速的参数、多重散射率。
在本说明书的含义内的非均匀介质包括任何包含散射体的介质，该散射体相对于周围介质表现出声阻抗的中断，所述周围介质也可能表现出在空间上不均匀的声速(例如多层介质)。在本说明书中描述的方法和系统尤其使得可以对这种非均匀介质进行非侵入性深度超声表征。
超声波波前通过非均匀介质传播的像差对应于来自该介质的波前与在理想情况下预期的波前之间的差异。这些像差可以例如与声速在所述介质中的空间波动相关联。例如，在回波描记术的情况下，所述声波在到达目标器官之前穿过各种类型的组织。由于软组织具有不同的机械性能，所述声速远非均匀。这导致声波波前的像差，该像差导致所得到的超声图像产生畸变，并因此导致其分辨率和对比度下降。这些像差可能会损及医疗检查。
最佳像差校正基是使所述声穿透区中包含的等晕域的尺寸最大化的基。例如，在薄的畸变器位于所述换能器阵列附近的情况下，像差校正基可以是换能器基，或者在平移不变的多层介质的情况下，像差校正基可以是平面波基。在空间不均匀的体积性畸变器的一般情况下，无论所考虑的像差校正基是什么，都存在多个等晕域。因此研究所述畸变矩阵使得可以确定与这些区域中的每一个相关联的像差定律。
所述聚焦基可以包含在相对于所述换能器阵列位于所述介质的给定深度处的单个焦平面中，或者可以通过同时考虑位于不同深度的焦点，以体积来定义。
记录了所述实验性反射矩阵的发射基可以等同于所述换能器基、所述平面波基(“傅立叶基”)或任何聚焦基。
可以在各个基之间定义反射矩阵或所述畸变矩阵，通过在时域中应用延迟线或通过在频域中与传播矩阵进行矩阵乘积来执行基变换。在说明书的其余部分，术语“聚焦”反射矩阵或“聚焦”畸变矩阵分别是指在输入处和输出处以聚焦基定义的反射矩阵或畸变矩阵。
记录的所述实验性反射矩阵是“实数”矩阵(也就是说，是由时域中的实系数组成的矩阵)，每个换能器记录的电信号是实数。然而，所述宽带聚焦反射矩阵是“复数”矩阵(也就是说，由复数分析信号形成的系数组成的矩阵)。可以通过仅保留所记录的电信号频谱的正频率来获得形成所述宽带聚焦反射矩阵的复数分析信号。这些复数分析信号的幅角与所述介质反射的波的相位相吻合。
模型介质或参考介质是例如声速在其中恒定的均匀介质。在理论上可以针对该参考介质建立作为传播矩阵的参考反射矩阵(或更简单而言，“参考矩阵”)(也就是说，针对该参考介质用每个焦平面上的虚拟平面镜获得的超声波的反射矩阵)。根据传播介质的理论知识程度，所述参考介质可以采用更复杂的形式(例如多层介质等)。
在聚焦基和像差校正基之间定义的第一宽带反射矩阵与针对模型介质定义的参考反射矩阵的相位共轭矩阵的逐项矩阵乘积或“哈达玛积”，等于从所述第一宽带反射矩阵中的每个元素的相位减去所述参考反射矩阵中的相应元素的相位。具体而言，一个矩阵或第一矩阵的相位共轭(在本说明书中更简称为“共轭”)的元素是复数，该复数的模与所述第一矩阵的相同，但相位相反。因此，从所述第一宽带反射矩阵的每个元素的相位减去预期的弹道分量(由所述模型介质定义)，这使得可以针对每个输入聚焦发射隔离由所述介质反射的场的畸变分量。申请人已经表明，分析所述畸变矩阵使得尤其可以辨别在所述声穿透区中包含的等晕域，并且可以在观察平面中确定与每个等晕域相关联的像差定律。
根据一个或多个示例性实施例，所述观察基是所述像差校正基。可以试图确定所述聚焦基和所述像差校正基之间的第一畸变矩阵，例如以量化和/或校正像差。
在这种情况下，所述第一畸变矩阵的确定可以包括：
-在输入处或输出处对所述第一宽带聚焦反射矩阵进行基变换，以获得在所述聚焦基和所述像差校正基之间确定的所述第一宽带反射矩阵；并且
-对所述第一宽带反射矩阵与在相同的基之间定义的参考反射矩阵的相位共轭矩阵进行逐项矩阵乘积。
根据一个或多个示例性实施例，所述观察基是聚焦基，并且参考“聚焦的”第一畸变矩阵。可以试图确定“聚焦”畸变矩阵(也就是说，在输入处和输出处以聚焦基定义的矩阵)，例如以确定扩散函数(或“点扩散函数”，PSF)的映射，扩散函数被定义为在输入处和输出处的聚焦基之间的空间脉冲响应。基于所述PSF，例如可以对在所研究介质的声穿透区中的聚焦质量(即像差水平)、声速或多重散射率进行局部量化。
可以通过从在所述聚焦基和所述像差校正基之间定义的所述第一畸变矩阵的基变换来获得所述聚焦的第一畸变矩阵。
或者，可以基于所述第一宽带聚焦反射矩阵，借助于所述第一宽带聚焦反射矩阵的每个行和/或列与以相同的基定义的参考反射矩阵的相同行和/或列之间的空间相关性来直接构建所述聚焦的第一畸变矩阵。基于所述聚焦的第一畸变矩阵，仅通过基变换就可以确定在所述聚焦基和所述像差校正基之间定义的所述第一畸变矩阵。
根据一个或多个示例性实施例，所述第一宽带聚焦反射矩阵的确定包括：
-在输入处和输出处对以频域表示的所述实验性反射矩阵进行基变换，以形成聚焦反射矩阵；
-在所述聚焦反射矩阵的所述第一谱带上进行相干加和，以形成所述第一宽带聚焦反射矩阵。
或者，可以通过将时移应用于所述实验性反射矩阵来获得所述第一宽带聚焦反射矩阵。然而，这种方法不能分析不同谱带上的实验数据，这在像差的频率分散的情况下可能是不利的。
根据一个或多个示例性实施例，所述非均匀介质的物理参数的至少一个映射的确定包括：
-确定所述第一畸变矩阵在所述聚焦基中的不变量，以在所述聚焦基中识别至少一个第一等晕域；
-对于识别出的每个第一等晕域，确定至少一个第一像差定律在所述校正基中的映射。
为了基于所述第一畸变矩阵来确定所述等晕域和相关的像差定律，可以采用几种方法。根据一个或多个示例性实施例，对这些不变量的确定包括所述第一畸变矩阵的奇异值分解、所述第一畸变矩阵归一化后(也就是说，其每个元素的模已经被归一化，但是其相位被保留)的奇异值分解或者所述第一畸变矩阵的相关矩阵归一化后(也就是说，所述第一畸变矩阵的相关矩阵的每个元素的模已经被归一化)的奇异值分解。
特别假定存在非均匀介质进行的镜面反射，申请人已经表明，可以通过所述第一畸变矩阵的奇异向量的线性组合来获得非均匀介质关于整个声穿透区的局部反射率的映射。
申请人还已经表明，所述畸变矩阵的奇异值分解使得可以过滤信号子空间(以其行和/或其列之间具有实质相关性为特征的矩阵)中的噪声子空间(其行和列之间不具有相关性的随机矩阵)，所述噪声子空间既包含实验噪声，又包含由聚焦基上游发生的多重散射事件引起的反射场的非相干贡献。
根据一个或多个示例性实施例，超声表征方法还包括在所述观察基中确定由所述一个或多个第一像差定律校正的第一反射矩阵。尤其是在假定样品发生随机散射式反射时，这使得可以确定所述介质在像差校正后的局部反射率的映射或“图像”。
根据一个或多个示例性实施例，基于所述声穿透区的所述校正后的第一反射矩阵，可以确定像差校正后的畸变矩阵，该校正后的畸变矩阵在所述像差校正基中与在所述像差校正基中确定的所述校正后的反射矩阵和所述参考反射矩阵的相位共轭矩阵的逐项矩阵乘积相对应。所述像差校正后的畸变矩阵使得可以通过确定其在聚焦基中的不变量，在识别出的所述第一等晕域中完善校正。还使得可以在聚焦基中识别至少一个第二等晕域，并且对于识别出的每个第二等晕域，确定第二像差定律在所述像差校正基中的映射。因此，取决于声穿透区中包含的等晕域数量，可以根据需要多次迭代该方法，以获得介质在像差校正后的局部反射率的映射或“图像”。当假定随机散射或中间反射(也就是说，镜面反射和随机散射性反射的混合反射)时，该迭代过程尤其适用。
根据一个或多个示例性实施例，对所述非均匀介质的物理参数的至少一个映射的确定包括针对识别出的至少一个第一等晕域，在所述聚焦基中确定第一扩散函数(PSF)。扩散函数对应于在像差校正平面中测量的像差定律的空间傅里叶变换。其在每个等晕域内在空间上都是不变的。根据一个示例，所述扩散函数可以从所述聚焦畸变矩阵(也就是说，在输入处和输出处以聚焦基所定义的)的奇异值分解获得。
根据一个或多个示例性实施例，所述非均匀介质的物理参数的至少一个映射的确定包括确定在所述声穿透区中的声速。这涉及例如研究所述畸变矩阵的第一奇异值作为速度模型的函数的变化，所述速度模型用于对参考矩阵进行建模。最佳速度模型是使畸变矩阵的第一奇异值最大化的模型。一种可选方案是研究聚焦畸变矩阵。具体而言，这使得可以研究扩散函数作为用于计算参考矩阵的速度模型的函数的变化。最佳速度模型是使波长归一化后的扩散函数的宽度最小化的模型。
根据一个或多个示例性实施例，所述非均匀介质的物理参数的至少一个映射的确定包括确定在所述声穿透区中的多重散射率。例如，这涉及研究像差已经被预先校正的聚焦畸变矩阵。超出焦点后获得的非相干强度(多重散射)的平均值与在焦点处的强度(所述畸变矩阵的中心行)的比值可以给出局部多重散射率。
根据一个或多个示例性实施例，所述表征方法还包括识别和/或消除反射场的镜面分量和/或在换能器阵列与非均匀介质的各个界面或层之间产生的多重反射。为此，可以在输入处和输出处将畸变矩阵投影到傅立叶基上。在此基础上，所述场的镜面反射分量和多重反射分量看起来是精确的反射角和入射角对。因此，它们可以很容易被过滤，仅保留所述反射场的随机散射分量(散斑)。然后，对随机散射分量的这种区分允许直接获得像差定律，该像差定律将被应用在输入处和输出处以校正所述反射矩阵并获得所述介质的最佳图像，这在镜面分量占主导时是不可能的。另外，在傅立叶平面中过滤所述畸变矩阵使得可以消除多重反射(探头与人体的表面层之间的多次回波)，该多重反射污染例如在浅深度处的超声图像。
根据一个或多个示例性实施例，所述表征方法还包括：
-基于所述实验性反射矩阵确定至少一个第二宽带聚焦反射矩阵的步骤，所述第二宽带聚焦反射矩阵在输入处和输出处以聚焦基在与所述第一谱带不同的第二谱带中定义；
-确定在所述聚焦基和所述观察基之间在所述第二谱带中定义的至少一个第二畸变矩阵的步骤，所述第二畸变矩阵是由在所述聚焦基和所述像差校正基之间定义的所述第二宽带反射矩阵与针对模型介质在与上述相同的基之间定义的参考反射矩阵的相位共轭矩阵的逐项矩阵乘积直接获得或在基变换之后获得；
在不同谱带上定义的畸变矩阵可用于确定相同数量的校正后的反射矩阵，然后可以将所述校正后的反射矩阵相加以形成单个宽带反射矩阵。该频谱分析使得可以校正由声速在所述介质中的频率分散而引起的色差。在不同谱带上获得的所述多个畸变矩阵也可用于确定所述声速或所述多重散射率的频谱映射。
根据第二方面，本说明书涉及一种用于对非均匀介质进行非侵入性超声表征的系统，该系统被配置为执行如上所述的所有示例性超声表征方法。根据第二方面的超声表征系统包括：
-至少一个第一换能器阵列，其被配置为在所述非均匀介质的区域中产生一系列入射超声波，并记录由所述区域反向散射的超声波随时间的变化；
-计算单元，其与所述至少一个第一换能器阵列耦合，并被配置为：
ο记录在输入处的发射基和输出处的换能器基之间定义的实验性反射矩阵；
ο基于所述实验性反射矩阵确定至少一个第一宽带聚焦反射矩阵，所述第一宽带聚焦反射矩阵在输入处和输出处以聚焦基在第一谱带中定义；
ο确定在所述聚焦基和观察基之间在所述谱带中定义的第一畸变矩阵，所述第一畸变矩阵是由在所述聚焦基和所述校正基之间确定的所述第一宽带反射矩阵与针对模型介质在与上述相同的基之间定义的第一参考反射矩阵的相位共轭矩阵的逐项矩阵乘积直接获得或在基变换之后获得；
ο基于所述第一畸变矩阵确定所述非均匀介质的物理参数的至少一个映射。
根据本说明书的表征系统可以包括既能发射和又能接收的至少一个换能器阵列，或者包括多个换能器阵列，其中一些专用于超声波的发射而其他的专用于超声波的接收。
附图说明
通过阅读下面的详细描述并参照附图，上述呈现的其他优点和技术特征将变得显而易见。
-图1A至图1C(已描述)示出了用于超声成像和定量化的已知发射/接收机制；
-图2(已描述)示出了根据现有技术基于自适应聚焦来校正超声成像中的像差的技术；
-图3A和图3B(已描述)示意性地示出了等晕域的概念；
-图4示出了用于实施根据本说明书的超声表征方法的示例性超声表征系统；
-图5示出了根据本说明书的方法的一个步骤，旨在根据一个示例性实施例，基于实验性反射矩阵来确定第一宽带聚焦反射矩阵；
-图6示出了为了以聚焦基记录的反射矩阵的投影而进行的基变换矩阵运算；
-图7A、图7B示出了根据本说明书的方法的一个步骤，旨在根据一个示例性实施例，基于宽带聚焦反射矩阵来确定第一畸变矩阵；
-图8示出了在没有像差以及具有像差的情况下体模的超声图像以及相应的聚焦畸变矩阵，该聚焦畸变矩阵在给定深度处确定并且对应于图像中的线。
-图9A和图9B分别示出了在如图8所示的具有像差的情况下与体模相对应的畸变矩阵的奇异值分布(图9A)，以及相关联的第一特征向量在傅立叶平面中的相位(图9B)；
-图10示出了畸变矩阵的示例性用途(体模上的像差水层产生单个等晕域的情况)，更精确而言，示出了未经/经过校正的超声图像、在给定深度处确定并对应于图像中的线的相应聚焦畸变矩阵以及未经/经过校正的扩散函数的表示，该扩散函数是在声穿透区上取平均值并且由聚焦畸变矩阵获得；
图11示出了畸变矩阵的示例性用途(在产生多个等晕域的人体乳房体模上)，并且更精确而言，示出了常规超声图像和使用畸变矩阵的第一向量和畸变矩阵的第二向量校正像差之后的超声图像；
图12A至图12C示出了畸变矩阵用于过滤多重反射的示例性用途，更精确而言，图12A是描述实验装置的示意图(超声探头和体模之间放置有机玻璃板)，图12B是针对均匀速度模型(c＝1800m/s)在从图像中去除镜面反射分量和多重反射前后的超声图像、聚焦畸变矩阵在过滤前后针对超声图像中的一条线的示例、投影到平面波基中的相同畸变矩阵，图12C是说明镜面反射分量和多重反射分量在平面波基中的特征的图。
图13示出了畸变矩阵用于声速的定量测量的示例性用途，并且更精确而言，示出了所研究体模的超声图像，聚焦畸变矩阵针对图像中的线以及均匀速度模型(c＝1540m/s)的示例、在输入处的焦平面和傅立叶平面之间的对应畸变矩阵、通过波长归一化的扩散函数的宽度随波速c变化的函数、畸变矩阵的归一化后的第一奇异值随波速c变化的函数；
图14示出了畸变矩阵用于建立在介质中的声速的映射的示例性用途，并且更精确而言，示出了体模针对均匀速度模型(c＝1540m/s)且穿过有机玻璃板的情况下得到的常规超声图像，以及过滤多重反射并针对图像的每一行优化速度模型后的相同图像、声速在深度范围内的分布；
图15示出了畸变矩阵用于多重散射的局部量化的示例性用途，并且更具体而言，示出了标准的体模超声图像、解释了在畸变矩阵中出现多重散射的情况的图、针对超声图像中的线的聚焦畸变矩阵、从前述矩阵推导的显示与多重散射分量相关联的非相干背景的扩散函数示例、超声图像的每个像素中的多重散射率；
将参照附图描述各个实施例，在这些实施例中，相似或相同的元件采用相同的附图标记。
具体实施方式
在下面的详细描述中，仅详细描述了一些实施例以确保描述的清楚性，但是这些示例并不旨在限制从本说明书中所示出的原理的一般范围。
本说明书中描述的各个实施例和方面可以以多种方式进行组合或简化。特别而言，除非另有说明，否则各种方法中的步骤可以重复、反向或并行执行。
图4示出了根据本说明书的用于实现对非均匀介质20的超声表征方法的示例性超声表征系统40。系统40包括换能器11的至少一个第一阵列10，例如线性或二维阵列。换能器是例如超声压电换能器，其可以采用与介质20接触的刚性棒的常规形式。换能器阵列形成例如探头设备41的一部分。换能器阵列与计算单元42连接，该计算单元42本身可以与显示设备43连接；计算单元向/从每个换能器11传输并记录电信号。然后，超声换能器将这些电信号转换成超声波，反之亦然。计算单元42被配置用于实施计算或处理步骤，尤其是用于实施根据本说明书的方法的步骤。
在图4中，如同在本说明书的其余部分，参考了用于发射和接收的换能器阵列，应当理解，在更一般的情况下，可以同时使用多个换能器阵列。这些换能器阵列可以是既发射也接收，或者仅其中某些用于发射而其他用于接收。
在本说明书中，当参考尤其用于实施方法步骤的计算或处理步骤时，应理解，每个计算或处理步骤可以由软件、硬件、固件、微代码或这些技术的任何适当组合来实现。当使用软件时，每个计算或处理步骤可以通过计算机程序指令或软件代码来实施。这些指令可以存储在或传输至计算机(或计算单元)可读的存储介质和/或可以由计算机(或计算单元)执行，以实施这些计算或处理步骤。
本说明书中所示的各种实验结果(图5至图15)取自两个超声体模，即再现生物组织的机械特性(c＝1540m/s)并产生生物组织的超声散斑特征的合成样品。第一体模使得可以量化成像系统在分辨率和对比度方面的性能。具体而言，该第一体模包括位于各种深度和横向位置的近点状目标以及直径为5mm的非均匀回声内含物(图5至图10以及图12至图15)。第二体模再现了乳腺组织的形状和特性，并显示出明亮的回声，模仿了同一组织内部的微钙化现象(图11)。
本说明书描述了用于对非均匀样品进行非侵入性超声表征的方法和系统。这些方法和系统基于至少一个第一矩阵(在本说明书的其余部分被称为“畸变矩阵”)的确定。根据一些示例，图5至图7示出了用于确定畸变矩阵的方法的步骤。
如图5所示，步骤51，非侵入性超声表征方法包括通过换能器11(图4)的阵列10在非均匀介质的区域中产生一系列入射超声波的步骤以及记录在输入处的发射基i和输出处的换能器基u之间定义的实验性反射矩阵Rui(t)的步骤。如图1A至图1C所述，输入处的发射基i例如是换能器基u、平面波基θ或聚焦基r。在图5的示例中，例如是在输入处的平面波基和输出处的换能器基之间记录的实验性反射矩阵Ruθ(t)。
所述方法还包括基于所述实验性反射矩阵Rui(t)确定至少一个第一宽带反射矩阵Rrr(Δω1)的步骤，所述第一宽带反射矩阵Rrr(Δω1)在输入处和输出处以聚焦基r在第一谱带Δω1中定义；根据一个示例性实施例，通过图5和图6示出了在构建“畸变矩阵”之前的该步骤。
在图5的示例中，离散傅立叶变换操作52使得可以在频域Ruθ(ω)中表达实验性反射矩阵。对实验性反射矩阵进行“宽带聚焦”的步骤53包括将反射矩阵Ruθ(ω)在输入处和输出处投影到焦点基中，以获得聚焦反射矩阵Rrr(ω)，然后在谱带Δω1中上进行相干加和。然后获得宽带聚焦反射矩阵，其中示出了在给定深度z处的一个示例(图像54)。
图6示出了将反射矩阵Ruθ(ω)在输入处和输出处投影到焦点基r中。举例来说，图像63表示在第一体模中获取的矩阵Ruθ(ω)在频率为f＝3.48MHz处的相位，该第一体模包括近点状目标以及直径为5mm的非均匀回声内含物，该第一体模的超声图像如图5所示(图像56)。这相当于在发射(61)和接收(62)时执行波束成形，以对来自同一散射体的回声进行相干加和。在数学上，可以使用传播矩阵来描述这些操作，所述传播矩阵描述波在介质中的向外和返回路径。因此，通过矩阵G′0(ω)＝[G′0(rin,θin,ω)],，对从入射平面波基到聚焦基的转变进行建模，所述矩阵的元素由下式给出：
G′0(rin,θin,ω)＝exp[jk(zincosθin+xinsinθin)] (1)
其中rin＝(xin,zin)是介质中的焦点的位置向量，θin是发射波的入射角，并且k＝ω/c是波数。从介质反射到换能器阵列的回声路径可以通过传播矩阵G0(ω)＝[G0(rout,uout,ω)]来描述，该矩阵的元素由下式给出：
其中rout＝(xout,zout)是介质中的焦点的位置向量，并且uout是换能器的位置向量。上式是Rayleigh-Sommerfeld积分的矩阵平移。其第一项是波动方程的二维格林函数，或者在这里是其渐近形式。右边项是换能器阵列的法线与径向向量(uout-rout)之间的方向余弦。
发射和接收时波束成形的运算可以解释为回声向其相关联的散射体的反向传播。在频域中，这理解为传播矩阵的相位共轭。因此，可以通过以下矩阵乘积从实验性反射矩阵推导出聚焦反射矩阵Rrr：
其中符号×表示矩阵乘积，指数
表示共轭转置矩阵。当然，也可以基于在输入处以另一个基(例如换能器基座)表达的实验性反射矩阵来将反射矩阵在输入和输出处投影到焦点基中。对此，在式(3)中使用传播矩阵G0来代替传播矩阵G′0就足够了。
下一步是对在谱带Δω1上获得的反射矩阵进行加和，以获得宽带聚焦反射矩阵
τ是与以聚焦基表示的虚拟换能器阵列相关联的回声时间。宽带聚焦反射矩阵中的每个系数R(rout,rin,τ)对应于在rin处的虚拟换能器发射脉冲后立即由在rout处的虚拟换能器记录的分析信号(图55，图5)。弹道时间τ＝0时的矩阵Rrr(τ＝0)的对角线对应于发射和接收时聚焦在同一点上的共聚焦信号。因此，如果聚焦基被限制在距离实际换能器阵列的距离为z的焦平面上，将直接给出超声图像中的一条线(图像56)。然而，如果聚焦基包含要成像的体积中的所有焦点，则Rrr(τ＝0)的对角线直接给出整个超声图像。在后一种情况下，反射矩阵是分块矩阵，仅其对角线块为非零。这些对应于介质每个深度z处的聚焦反射矩阵。在本说明书的其余部分，为了简化表示，将宽带反射矩阵简称为Rrr，不再提及时间τ，因为始终将其设置为弹道时间τ＝0。然而，对于某些应用，尤其是对于高阶像差的校正、多重散射的研究或声速的测量，也可以在时间τ≠0时检查反射矩阵和畸变矩阵。
图7A、图7B示出了根据本说明书的方法的一个步骤，旨在根据一个示例性实施例，基于宽带聚焦反射矩阵Rrr来确定第一畸变矩阵D。
为了估计波前的像差分量，首先，选择用于校正像差的最佳基。例如，在薄的畸变器位于换能器阵列附近的情况下，换能器基看来是最合适的。在所述表征方法的该说明性示例中所选择的情况下，非均匀介质是平移不变的多层介质，平面波基似乎更合适。这例如是对肝脏采用回波描记术的情况，其中声波必须穿过一系列的肌肉层和脂肪组织层。在下文中考虑该示例，但是本说明书绝不限于从傅立叶基对像差进行校正。
图像71(图7A)示出了在存在像差的情况下记录的宽带聚焦反射矩阵Rrr的示例，该像差是通过在换能器阵列和体模(c＝1540m/s)之间插入10mm厚的水层(c＝1480m/s)生成的。首先，将反射矩阵Rrr在输出(或输入)处的基变更为平面波基，使得可以形成宽带反射矩阵Rθr，使得
Rθr＝tG′0(ωc)×Rrr (5)
其中考虑了谱带Δν1的中心频率ωc处的矩阵G′0。在物理上，由此获得的矩阵Rθr的每一行对应于针对每个声穿透点rin在远场中反向散射的波前。每个波前都有两个分量：几何分量，其与换能器平面和焦平面之间的波的衍射关联；以及畸变分量，其与波前向外和返回时所经历的像差关联。为了分离这两个分量，将矩阵Rθr与会在模型介质(例如没有像差(仅具有几何分量)的均匀介质)中获得的参考矩阵进行比较。在没有像差的均匀介质的情况下，该参考矩阵不过是在均匀介质中的传播矩阵G′0。然后从在聚焦基和像差校正基之间确定的所述第一宽带反射矩阵Rθr与在上述相同的基之间在中心频率ωc处定义的参考反射矩阵
之间的哈达玛矩阵乘积或逐项矩阵乘积来创建新矩阵(在本说明书中称为畸变矩阵)：
其中符号o表示逐项矩阵乘积(哈达玛矩阵乘积)，指数
表示共轭转置矩阵。回想一下，矩阵tG′0的相位共轭矩阵
的元素是复数，该复数的模与tG′0中的模相同，但相位相反。就矩阵系数而言，该式子写为：
上面执行的运算由图7B示出。包括从测得的波前减去确定的几何分量，从而隔离其畸变分量。矩阵Dθr的每一列都包含针对每个焦点rin从远场看到的波前的畸变分量。因此，在实践中，为了隔离反射的波前的畸变分量，从反射场的相位中减去在没有像差的理想情况下(对于模型介质而言)获得的场的相位。
在模型介质不均匀的情况下，参考反射矩阵G′0的表达式变得比式(1)给出的更复杂。例如，在多层介质的情况下，可以通过与波在每一层中的传播相关联的不同传输矩阵之间的矩阵乘积来分析地计算矩阵G′0。对于更复杂的模型介质(通常不是平移不变的介质，例如弯曲界面)，可以通过数值模拟或半分析计算来确定矩阵G′0。注意，如果参考反射矩阵的系数具有不均一的模，则可以通过对这些系数中的每一个施加恒定的模但保持其相位来对所述参考矩阵进行归一化。
图7A中的图像73显示了从反射矩阵推导出的畸变矩阵Dθr的相位，其模在图像71示出。后者对应于在存在像差水层的情况下在体模上产生的超声图像83(图8)中的线。图像74示出了针对焦平面的每个点隔离出的存在像差的波前，及其存储在畸变矩阵中的方式。
也可以在焦平面中研究畸变矩阵。可以通过以下基变换，基于在观察平面中表示的畸变矩阵Dθr来获得“焦平面”畸变矩阵Drr：
其中
是归一化后的传播矩阵，其元素由下式给出
因此，式(8)就矩阵系数而言重写如下：
其中，kc是中心频率ωc处的波数。在像差校正基是傅立叶基的当前情况下，仅通过Dθr在输出处的空间傅立叶变换来简单地推导出聚焦畸变矩阵Drr。然而，应该记住，如果校正平面不再是傅立叶平面，则两个矩阵之间的关系将不同[式(8)]。
Drr的列(其模在图像76中示出)对应于在重新集中在焦点rin上的图像平面中反射的场。因此，矩阵Drr给出了在超声图像的每个点处的反射扩散函数的变化，该扩散函数是成像系统(也就是说，重新集中在焦点上的聚焦点)的空间脉冲响应。如下所述，扩散函数使得可以量化和表征像差以及由样品在焦平面上游引起的多重散射。
作为说明，图8通过回溯图像80的实验图示出了在没有像差的情况下获得的体模的超声图像81。图像82示出了在给定深度处确定的并对应于图像81中的线的聚焦畸变矩阵。图像83是在存在像差(像差水层22，图85)的情况下获得的体模的超声图像，并且图像84示出了在给定深度处确定的并对应于图像83中的线的聚焦畸变矩阵。在没有像差的实验中，由介质反射的所有能量都出现在畸变矩阵的中心行(rin＝rout)：扩散函数的宽度δ接近于衍射施加的极限(δ～λz/D，其中D是换能器阵列的大小)，这解释了超声图像在分辨率和对比度方面的出色质量。然而，在存在像差水层的情况下，反射的能量会扩展到中央行的外部。空间脉冲响应的降级意味着超声图像的分辨率显著降低。该示例显示了“焦平面”畸变矩阵如何能够独立于超声图像的定性分析来量化在介质中聚焦的质量。
为了进一步分析畸变矩阵Dθr，有利的是，可以一方面考虑样品的反射主要是镜面反射的情况，另一方面考虑反射是随机散射的情况。实际上，通常是结合了随机散射性反射和镜面反射的中间情况。在这两种情况下，畸变矩阵Dθr在其列之间都显示出相关性。这些相关性对应于针对来自相同等晕域的波前的相同畸变模式的重复。对于引起例如散斑类型(未分辨散射体的随机分布)的随机散射性反射的样品，可以研究畸变矩阵Dθr的归一化后的相关矩阵，以便更高效地提取这些相同的相关性。无论哪种情况，如以下将更详细描述的那样，例如通过奇异值分解(SVD)搜索矩阵Dθr中的不变量，可以针对焦平面中的每个点提取畸变器的复数透射率，并因此可以最佳地校正声穿透区中包含的每个等晕域中的像差。
假定样品引起随机散射性反射，下面描述使用畸变矩阵来确定一个或多个像差定律在观察基中的映射(尤其是建立声穿透区的局部反射的映射)的示例。
为此，寻求畸变矩阵中的不变量(换言之，在声穿透区的等晕域内在空间上不变的像差定律)。本领域技术人员已知各种方法用于搜索这种矩阵中的不变量，例如奇异值分解(或“SVD”)或主成分分析(“PCA”)。
对于提取矩阵的行或列之间的相关性来说，奇异值分解是强大工具。数学上，大小为N×M的矩阵Dθr的SVD书写如下：
U和V是大小为N×M的单位矩阵，其列Ui和Vi对应于输出特征向量和输入特征向量。每个输出特征向量Ui以由角度θout标识的傅里叶平面来定义。因此，每个输入特征向量Vi以由向量rin标识的焦平面来定义。Σ是大小为N×M的矩阵，只有其对角线元素不为零：
矩阵Σ的对角线元素是矩阵Dθr的实数奇异值σi，该奇异值σi为正数并按降序排列：
因此，矩阵Dθr的系数(θout,rin)被写作以下总和：
其中L＝min(M,N)。
SVD主要将矩阵分解为两个子空间：信号子空间(以其行和/或列之间具有实质相关性为特征的矩阵)和噪声子空间(在其行和列之间不具有相关性的随机矩阵)。信号子空间与最大的奇异值相关联，而噪声子空间与最小的奇异值相关联。一方面，D的SVD因此可以过滤噪声子空间，该噪声子空间既包含实验噪声也包含由多重散射事件引起的反射场的非相干贡献。另一方面，信号子空间的每个奇异状态将使得可以根据输出特征向量Ui来提取波针对图像的每个等晕域在像差校正基上经历的畸变。特别而言，在单个等晕域的情况下，可以示出第一输出特征向量U1直接给出了校正平面和焦平面之间的畸变器的复数透射率A：
U1(θout)＝A(θout) (15)
下面的图示出了畸变矩阵的各种应用，无论该畸变矩阵是以焦点基定义还是在焦点基与像差校正基之间定义。
图9A和图9B示出了将SVD应用于图8的图像84中所示的畸变矩阵Dθr的结果。Dθr的奇异值σi的谱图(图9A)示出了第一奇异值σ1的优势。这表明存在单个等晕域。相关联的第一特征向量U1(图9B)显示出抛物线相位，这是散焦缺陷的特征，可以通过其中声速(1490m/s)与在体模中的声速(1540m/s)有很大不同的水层的存来解释。
为了校正宽带反射矩阵Rrr(图7A中的71)，在此示例中，将矩阵Rrr投影到像差校正平面(此处为傅里叶平面)中：
Rθθ＝tG′0(ωc)×Rrr×G′0 (16)
然后通过在矩阵Rθθ的输入和输出处应用U1的相位共轭来校正像差：
其中R′θθ是校正后的反射矩阵，函数arg{…}表示花括号内的向量的相位。关系式(18)转换为矩阵系数，如下所示：
在物理上，相位共轭运算包括重新发射由与所测得的畸变的相位相反的相位调制的波前，以及在输出处施加类似的校正。因此，该操作使得可以完美补偿波在其向外和返回路径上积累的相位畸变。然后可以经由在输入和输出处从瞳孔平面到焦平面的基变换从R′θθ推导出校正矩阵R′rr。可以从矩阵R′rr的对角线(rin＝rout)推导出校正后的共焦图像I′：
I′(rin)＝R′(rin,rin) (19)
因此，I′(rin)是对样品反射率ρ(rin)的可靠估计。
图10通过比较像差的超声图像101和根据式(17)描述的校正获得的超声图像103来示出该校正的益处。特别是借助于在校正像差后看起来对比度更强并且分辨率更好的回声目标，可以观察到这两个图像之间的定性增益。
然而，在回波描记术中，大多数时候没有回声目标来判断图像的分辨率。为了量化后者，可以使用聚焦畸变矩阵。从初始的反射矩阵推导出的Drr以特定深度在图像102中呈现。从校正后的反射矩阵推导出的D′rr以相同的深度在图像104中呈现。与第一种情况相比，在第二种情况下，能量更加集中在聚焦畸变矩阵的中心行周围。
基于聚焦畸变矩阵，还可以确定扩散函数。这对应于在像差校正平面中测量的像差定律的空间傅里叶变换。其在每个等晕域内在空间上都是不变的。该扩散函数也可以从聚焦畸变矩阵的奇异值分解获得。如果矩阵G′0是酉矩阵(傅立叶基和聚焦基彼此共轭的情况)，则矩阵Drr具有与Dθr类似的奇异值分解。输出特征向量(被称为Wi)通过简单的傅立叶变换与特征向量Ui直接关联：
图10的图表105比较了由此获得的扩散函数在校正前后的模。在分辨率方面的增益是明显的，并且在校正像差之后达到了衍射极限。
在进行的实验中，测得的像差不依赖于频率，因为在所研究的频带内，声速在体模和水中的变化很小。实际上，情况并非总是如此，并且畸变矩阵的频率分析可能被证明是明智的，以便确定能最佳地校正像差的空频自适应滤波器(或空时等效)。具体而言，可以在以不同中心频率ωc,i为中心的不同谱带Δωi上研究畸变矩阵D。可以在这些谱带中的每一个上提取像差定律U1(ωc,i)，并且可以从其推导校正后的反射矩阵R′rr(ωc,i)[式(17)]。在每个谱带Δωi独立校正的反射矩阵的相干加和使得可以获得超宽带校正反射矩阵R′rr：
该像差频率校正在波束成形中等效于由空时滤波器发射和测量的聚焦信号的卷积U1(τ)，其中
到目前为止，已经利用了平移不变畸变层的情况，因此仅生成一个等晕域。现在描述畸变矩阵的第二示例性用途，其出于非平移不变的体积性畸变器的情况下，从而产生多个等晕域。
在这种情况下，倾向于采用迭代方法(后处理)来校正像差。该过程的步骤#0包括通过畸变矩阵Dur的第一特征向量U1的相位共轭来在输出处校正宽带反射矩阵Rθr：
其中
是所得的校正后的矩阵。该初始步骤使得可以在整个图像上对像差进行整体校正。
步骤#1在于重新计算从
推导出的新畸变矩阵
该
是这次在输入处的傅立叶平面和输出处的焦平面之间定义的。对象的随机性质意味着在傅立叶平面上研究
的相关矩阵(即
)更有意义。该矩阵B(1)的相位的奇异值分解(exp[jargB(1)}])给出了一组新的特征向量
这些特征向量与图像的每个等晕域相关联。因此，这次可以在输入处校正相应的反射矩阵：
接下来的步骤包括通过交替校正在输出处(偶数迭代)和输入处(奇数迭代)的残留像差来再现相同的过程。然而，在每个步骤中，仍然使用在瞳孔平面中的畸变矩阵的相关矩阵的相位(exp[jargB(n)}])的第一特征向量
进行校正。具体而言，在步骤#1中进行等晕域的选择。根据校正是在输出处还是在输入处进行，在瞳孔平面中的相关矩阵B(n)为：对于偶数n(输出)，
对于奇数n(输入)，
在所述过程的每个步骤中，可以从在焦平面中表示的反射矩阵
的对角线推断出视场图像。在实践中，只需进行几次迭代就可以针对所选等晕域获得最佳校正。然后可以通过组合针对每个等晕域确定的校正来获得整个视场的图像。
图11通过在乳房体模上进行实验来示出根据本说明书的方法的益处，所述乳房体模产生多个等晕域并且在20至25mm之间的深度处包含微钙化。图11示出了常规超声图像111、使用归一化后的相关矩阵exp[jarg{B(1)}]的第一特征向量
校正像差之后的超声图像112以及关于前述超声图像的一部分的特写113。还示出了在使用第二向量
校正像差之后的超声图像116，并且图像117中是关于前述超声图像的一部分的特写。从这些图像中可以看出，第一特征向量
与包含微钙化(在20到25mm之间的深度处)的等晕域相关联。具体而言，在该深度，超声图像(112)对比度更高且分辨更好。相反，第二特征向量
可以在较浅的深度(10到15mm之间)最佳地校正图像(117)。
对样品中每个点的像差定律的了解，不仅可用于成像，而且可用于将超声波物理性地聚焦在样品的任意点上。例如，在图2的实验设备上，可以经由空时滤波器向换能器阵列施加从畸变矩阵的特征向量推导出的初始聚焦信号的卷积。在存在单个等晕域的情况下，将考虑第一特征向量U1[式(22)]的傅立叶变换。注意，如果像差定律是频率稳定的，则空时滤波器将仅包含应用与换能器基中的特征向量U1的相位成正比的延迟定律Δt
Δt(u)＝-arg{U1(u)}/ωc (25)
在存在多个等晕域的情况下，通过迭代过程确定像差定律(请参见先前段落)。根据向量W，将样品在输入处获得的所有像差校正汇总在一起：
在镜面反射模式下，并且在超声图像包含在单个等晕域中的情况下，可以示出畸变矩阵的第一输出特征向量U1给出了在向外以及返回时由畸变器累积引起的畸变：
U1(θout)＝A(θout)A(θin)δ(θout+θin) (27)
其中A(θout)和A(θin)是波前在向外以及返回时所经历的并且投影到傅立叶平面中的畸变。前式中的狄拉克分布δ反映了这样一个事实，在镜面反射模式下，入射角为θin的平面波在输出处被反射为角度为θout＝-θin的平面波。此外，可以示出输入特征向量V1关于其部分可以直接提供对象的反射率ρ：
V1(rin)＝ρ(rin) (28)
通过在输入或输出处仅施加一次特征向量U1的相位共轭应，可以对反射矩阵进行像差校正：
R′ur＝exp(-j×arg{U1})oRur (29)
在超声图像包含P>1个等晕域的情况下，畸变矩阵的秩为P。输入特征向量Vi分解焦平面中越过不同等晕域的对象，并且与傅立叶平面中的不同像差定律Ui相关联。通过相关特征值
加权的特征向量Vi的模的线性组合，最终给出校正由上游引起的像差后的视场图像
然而，反射器永远不会是完美的镜面反射，因此其随机散射分量将无法通过式(29)进行最佳校正。此外，镜面反射器可能在其界面之间引起多重反射，这可能会使超声图像模糊。因此，必须能够分离由介质反射的场的镜面反射分量和随机散射分量。根据一个或多个示例性实施例，表征方法因此还包括识别和/或消除反射场的镜面分量和/或在换能器阵列与非均匀介质的各个界面或层之间产生的多重反射。为此，可以在输入处和输出处将畸变矩阵投影到像差校正基上。可以通过以下基变换，基于在观察平面中表示的畸变矩阵Dθr来获得“焦平面”畸变矩阵Dθθ：
Dθθ＝Dθr tG′0(ωc) (31)
因此，前式就矩阵系数而言重写如下：
如图13C所示，对于精确的反射角和入射角对θin和θout，出现镜面反射分量和多重反射分量，使得
θin+θout＝Δθ (33)
其中Δθ对应于镜面反射器和换能器阵列之间的角度。如果镜面反射器平行于换能器阵列，则Δθ为零。因此，可以容易地滤除镜面反射分量和多重反射分量，并且仅保留反射场的随机散射分量(散斑)，以便精确确定由有机玻璃板引起的像差并进行校正。然后，对随机散射分量的这种区分使得可以直接获得像差定律，该像差定律要被用在输入处和输出处，以校正所述反射矩阵并获得所述介质的最佳图像[式(17)]，这在镜面分量占主导时是不可能的[式(29)]。另外，在傅立叶平面中过滤所述畸变矩阵使得可以消除多重反射(探头与人体的表面层之间的多次回波)，该多重反射污染例如在浅深度处的超声图像。
图12A至图12C示出了使用畸变矩阵对多重反射进行过滤以及该畸变矩阵后续在像差校正和声速映射中的应用。图12A示出了实验装置的图，即有机玻璃板121放置在超声探头41和体模20之间。相应的超声图像122是使用均匀速度模型(c＝1800m/s)构造的。该图像由于有机玻璃板引起的多重反射和像差而大大降级。图12B示出了针对超声图像的深度(z＝25mm)的畸变矩阵的各种形式。图像124示出了其中扩散函数因有机玻璃板的存在而特别降级的聚焦畸变矩阵的模。图像125示出了傅立叶平面中的畸变矩阵。对于θin+θout＝0，强镜面分量很明显。图像126示出了应用镜面滤波器后的相同畸变矩阵，所述镜面滤波器在于消除|θin+θout|＜1/60rad的场的所有分量。图像127示出了过滤镜面反射分量后的聚焦畸变矩阵。与图像126(在滤波之前)的比较示出了畸变矩阵虽然仍然具有强烈像差，却具有典型的随机散射模式(散斑)的随机外观。如此过滤后的畸变矩阵的中心行给出了图12A的超声图像123。与原始超声图像122的比较表明，所施加的滤波器已成功地消除了由有机玻璃板引起的多重反射。即使图像123仍然受到像差的影响，但体模中的每个目标现在都可见，而初始图像122中却不是这种情况。
根据一个或多个示例性实施例，所述非均匀介质的物理参数的至少一个映射的确定包括确定在所述声穿透区中的声速。这涉及例如研究所述畸变矩阵的第一奇异值作为速度模型的函数的变化，该速度模型用于对参考矩阵进行建模。最佳速度模型是使畸变矩阵的第一奇异值最大化的模型。一种可选方案是研究聚焦畸变矩阵。具体而言，这使得可以研究扩散函数作为用于计算参考矩阵的速度模型的函数的变化。最佳速度模型是使波长归一化后的扩散函数宽度最小化的模型。
因此，图13示出了畸变矩阵用于声速的定量测量的示例性用途。图13示出了已经在图像81(图8)中示出的体模在没有像差的情况下的超声图像131。图像132示出了聚焦畸变矩阵针对图像中的线以及均匀速度模型(c＝1540m/s)的示例，图像133示出了在输入处的焦平面和傅立叶平面之间的对应的畸变矩阵。图像134示出了由波长λ归一化的扩散函数的宽度δ随波速c变化的函数，图像135示出了畸变矩阵的归一化后的第一奇异值(
)随波速c变化的函数。
图像134中的比率δ/λ清楚地显示为至少大约c＝1539m/s，而体模制作商给出的速度为1540m/s。适当的速度模型意味着超声波在发射和接收时的最佳聚焦，以及相应地有更好的分辨率的扩散函数。因此，聚焦畸变矩阵为定量测量体模中的声速提供了相关标准。类似地，图表135中的畸变矩阵的第一奇异值
在速度c＝1542m/s附近显示为最大值。同样，这里测得的速度接近制造商给出的速度。具体而言，用于计算畸变矩阵的速度模型越接近现实，该畸变矩阵就越接近秩为1的奇异矩阵。该实验示出了畸变矩阵如何为声速的定量测量提供各种标准。为了进行概念验证，这里仅限于均匀介质的情况，但是，当然可以使用畸变矩阵及其参数δ/λ或
作为信息反馈，以构建介质中声速的局部映射。第一步将在图14中示出。
除了介质的反射率或样品每个点的像差定律知识外，畸变矩阵还可以对所研究样品的光学指标进行3D层析成像。从“焦平面”畸变矩阵Drr，可以获得样品中的扩散函数以及等晕域。不是直接校正图像，该想法是改变参考介质，畸变矩阵是基于该参考介质从反射矩阵计算的。这导致了对畸变矩阵的迭代方法，其中使参考介质沿更复杂的模型(例如多层介质)的方向变化，以减小扩散函数的空间范围并增大等晕域的大小。通过越来越深入样品，可以逐步重建声速的三维映射，同时获得介质反射率的图像，参考介质越接近现实，所述图像就越可靠。
图14示出了畸变矩阵用于对媒介中的声速建立映射的示例性用途。图14A示出了在均匀速度模型(c＝1540m/s)且穿过有机玻璃板的情况下，体模的常规超声图像141。由于有机玻璃板引起的像差和多重反射，该超声图像质量较差。相比之下，在超声图像142中成功地示出体模的反射率。后者的构建首先是通过过滤图像中的多重反射(图12)，然后通过针对图像中的每一行优化速度模型(图13)。这些操作使得可以获得仅受衍射限制的最佳分辨率。由于镜面反射分量的过滤，仅目标的低空间频率很难重新组合。图143示出了针对在深度上积分的声速而获得的曲线。该速度在进入有机玻璃层时(z＝14mm)急剧增加，然后随着离开有机玻璃层(z＝17mm)而逐渐降低。然后，一种算法使得可以解决反向问题，该反向问题包括从该积分的速度分布图推断出声速随深度变化的精确分布图。在有机玻璃中测得的声速为2700m/s，如所预计的，在体模中测得的声速为1540m/s。
根据一个或多个示例性实施例，所述非均匀介质的物理参数的至少一个映射的确定包括确定在所述声穿透区中的多重散射率。例如，这涉及研究像差已经被预先校正的聚焦畸变矩阵。超出焦点后获得的非相干强度(多重散射)的平均值与焦点处的强度(所述畸变矩阵的中心行)的比值可以给出局部多重散射率。
图15示出了畸变矩阵对于多重散射的局部量化的优势。图15显示了标准体模超声图像151。图152解释了在畸变矩阵中表现为非相干场的多个散射路径的本质。这些发生在焦平面的上游，因为它们与在焦平面中单一散射的波具有相同的飞行次数。聚焦畸变矩阵在图像153中示出，并且对应于超声图像中的线(z＝75mm)。其中心行(rin＝rout)清楚地显示了单一散射贡献的超强度特征。然而，无论距离|rin-rout|如何，都具有均匀强度的非相干背景。该非相干背景可能是由超声测量中的噪声引起的，或者可能与多重散射有关。为了知道其起源，可以检查波的空间互易性。如果R(rin,rout)＝R(rout,rin)，则该不相干背景与物理起源的信号有关，因此为多重散射。否则，只能是噪音。
聚焦畸变矩阵Drr提供成像系统在视场的任何点处的扩散函数154。可以基于非相干背景的水平(|rout-rin|＞＞δ)与Drr中心行上(rout＝rin)的信号水平之比来局部测量多重散射率γ(rin)：
其中符号<…>表示除中心元素以外(|rout-rin|＞＞δ)沿Drr的每一行的平均值。IM(rin)＝＜|D(rout≠rin,rin)|2＞是点rin处的多重散射强度(非相干背景)。IS(rin)是点rin处的单一散射强度。式(34)指示在中心行(rout＝rin)上，存在以下关系：|D(rout＝rin,rin)|2＝IS(rin)+2IM(rin)。因子2来自相干反向散射现象，该现象由rin中虚源上的反向散射波由于多重散射贡献引起的2倍超强度组成。该现象是由沿着相同轨迹但方向相反(反向路径)的多重散射波之间的相长干涉引起的。[A.Aubry et al.,Appl.Phys.Lett.92,124101,2008]。
图15的图像155示出了在超声图像151的每个点处测得的多重散射率。在较浅的深度，该程度不可忽略，但与电子噪声(信号的非互易性)相关。另一方面，该程度在很深处的增加显然与多重散射现象有关。特别而言，在超过60mm的深度，多重散射率超过50％。
该测量结果已通过先前的体内研究得到证实[A.Aubry et al.,J.Acous.Soc.Am.129,225-233,2011]。本发明的贡献在于聚焦畸变矩阵使得可以局部测量多重散射率，而在此之前，只能获得在每个回波时间或深度处的平均值。[A.Aubry et al.,J.Acous.Soc.Am.129,225-233,2011]。
通过花费比弹道时间[式(4)]更长的时间，扩散函数还可以跟着介质中扩散晕的增长，并从中得到多重散射波的传输参数的局部测量(散射系数或传输平均自由程)。研究短时间内的扩散晕可以获得比现有技术在随机散射层析成像中所获得的更好的空间分辨率[A.Aubry et al.,Appl.Phys.Lett.92,124101,2008]。
尽管通过多个详细的示例性实施例进行了描述，但是用于非侵入性超声表征的方法和系统包括不同的变型、修改和改进，这对于本领域技术人员而言将是显而易见的，应理解，这些不同的变型、修改和改进都落入由所附权利要求书限定的本发明范围内。
Claims (16)
1.一种用于对非均匀介质进行非侵入性超声表征的方法，该方法包括：
-通过至少一个第一换能器阵列在所述非均匀介质的区域中产生一系列入射超声波的步骤；
-记录在输入处的发射基(i)和输出处的换能器基(u)之间定义的实验性反射矩阵(Rui(t))的步骤；
-基于所述实验性反射矩阵确定至少一个第一宽带反射矩阵(Rrr(Δω1))的步骤，所述第一宽带反射矩阵在输入和输出处以聚焦基在第一谱带(Δω1)中定义；
-确定在所述聚焦基和观察基之间在所述第一谱带中定义的第一畸变矩阵(Dθr,Drr)的步骤，所述第一畸变矩阵是由在所述聚焦基和像差校正基之间定义的所述第一宽带反射矩阵(Rθr(Δω1))与针对模型介质在与上述相同的基之间定义的参考反射矩阵的相位共轭矩阵的逐项矩阵乘积直接获得或在基变换之后获得；
-基于所述第一畸变矩阵确定所述非均匀介质的物理参数的至少一个映射的步骤。
2.根据权利要求1所述的方法，其中，所述第一宽带反射矩阵(Rrr(Δω1))的确定包括：
-在输入和输出处对在频域中表示的所述实验性反射矩阵(Rθu(ω))进行基变换，以形成在输入和输出处以聚焦基在所述频域中定义的反射矩阵(Rrr(ω))；
-对在输入和输出处以聚焦基定义的所述反射矩阵的所述第一给定谱带(Δω1)进行相干加和，以形成在输入和输出处以聚焦基定义的所述第一宽带反射矩阵(Rrr(Δω1))。
3.根据前述权利要求中的任一项所述的方法，其中，所述第一畸变矩阵的确定包括：
-在输入处或输出处对所述第一宽带聚焦反射矩阵(Rrr)进行基变换，以获得在所述聚焦基和所述像差校正基之间确定的所述第一宽带反射矩阵(Rθr)；
-对在所述聚焦基和所述像差校正基之间定义的所述第一宽带反射矩阵(Rθr)与在上述相同的基之间定义的所述参考反射矩阵的相位共轭矩阵进行逐项矩阵乘积。
4.根据前述权利要求中的任一项所述的方法，其中：
-所述观察基是聚焦基；
-所述第一畸变矩阵在输入处和输出处以聚焦基定义，并且从在输入处和输出处以聚焦基定义的所述第一宽带反射矩阵(Rrr)直接获得，或者通过在所述聚焦基和所述像差校正基之间定义的所述第一畸变矩阵进行基变换获得。
5.根据权利要求4所述的方法，其中，所述非均匀介质的物理参数的至少一个映射的确定包括，从在输入处和输出处以聚焦基定义的所述第一畸变矩阵并且针对识别出的至少一个第一等晕域确定至少一个第一点扩散函数(PSF)。
6.根据前述权利要求中的任一项所述的方法，其中，所述非均匀介质的物理参数的至少一个映射的确定包括：
-确定所述第一畸变矩阵在所述聚焦基中的不变量，以在所述聚焦基中识别至少一个第一等晕域；
-对于识别出的每个第一等晕域，确定至少一个第一像差定律在所述校正基中的映射。
7.根据权利要求6所述的方法，其中，所述第一畸变矩阵的不变量的确定包括矩阵组中的至少一个矩阵的奇异值分解，所述矩阵组包括所述第一畸变矩阵、归一化后的所述第一畸变矩阵、所述第一畸变矩阵的归一化后的相关矩阵。
8.根据权利要求6和7中的任一项所述的方法，其中，所述非均匀介质的物理参数的至少一个映射的确定还包括在所述观察基中确定由所述一个或多个第一像差定律校正的第一宽带反射矩阵。
9.根据权利要求8所述的方法，还包括：
-确定校正后的畸变矩阵，该校正后的畸变矩阵通过所述校正后的第一宽带反射矩阵与所述参考反射矩阵的相位共轭矩阵的逐项矩阵乘积在所述聚焦基和所述像差校正基之间定义。
10.根据权利要求9所述的方法，还包括：
确定所述校正后的畸变矩阵在所述聚焦基中的不变量，以在所述焦平面中识别至少一个第二等晕域；
-对于识别出的每个第二等晕域，确定至少一个第二像差定律在所述校正基中的映射。
11.根据前述权利要求中的任一项所述的方法，还包括：
-基于所述实验性反射矩阵确定至少一个第二宽带反射矩阵(Rrr(Δω2))的步骤，所述第二宽带反射矩阵在输入处和输出处以聚焦基在与所述第一谱带不同的第二谱带中定义；
-确定在所述聚焦基和所述观察基之间在所述第二谱带中定义的至少一个第二畸变矩阵(Dθr(Δω2),Drr(Δω2))的步骤，所述第二畸变矩阵是由在所述聚焦基和所述像差校正基之间定义的所述第二宽带反射矩阵(Dθr(Δω2))与针对模型介质在与上述相同的基之间定义的参考反射矩阵的相位共轭矩阵的逐项矩阵乘积直接获得或在基变换之后获得。
12.根据权利要求11所述的方法，还包括：
-确定所述至少一个第二畸变矩阵(Dθr(Δω2))在所述聚焦基中的不变量，以在所述焦平面中识别至少一个第一等晕域；
-对于识别出的每个等晕域，确定至少一个第一像差定律在所述校正基中的映射；
-在所述观察基中，确定由所述一个或多个第一像差定律校正的至少一个第二宽带反射矩阵；
-将所述校正后的第一宽带反射矩阵与所述校正后的至少一个第二宽带反射矩阵相加，以形成校正后的超宽带反射矩阵。
13.根据前述权利要求中的任一项所述的方法，其中，所述非均匀介质的物理参数的至少一个映射的确定包括确定声穿透区中的声速。
14.根据前述权利要求中的任一项所述的方法，其中，所述非均匀介质的物理参数的至少一个映射的确定包括确定声穿透区中的多重散射率。
15.根据前述权利要求中的任一项所述的方法，还包括基于所述第一畸变矩阵来识别和/或消除镜面反射分量和/或多重反射分量，所述第一畸变矩阵在输入和输出处以平面波基表示。
16.一种用于对非均匀介质(20)进行非侵入性超声表征的系统(40)，包括：
-至少一个第一换能器阵列(10)，其被配置为在所述非均匀介质的区域中产生一系列入射超声波，并记录由所述区域反向散射的超声波随时间的变化；
-计算单元(42)，其与所述至少一个第一换能器阵列耦合，并被配置为：
○记录在输入处的发射基(i)和输出处的换能器基(u)之间定义的实验性反射矩阵(Rui(t))；
○基于所述实验性反射矩阵确定至少一个第一宽带反射矩阵(Rrr(Δω1))，所述第一宽带反射矩阵在输入和输出处以聚焦基在第一给定谱带(Δω1)中定义；
○确定在所述聚焦基和观察基之间在所述第一谱带中定义的第一畸变矩阵(Dθr，Drr)，所述第一畸变矩阵是由在所述聚焦基和像差校正基之间定义的所述第一宽带反射矩阵(Rθr(Δω1))与针对模型介质在与上述相同的基之间定义的参考反射矩阵的相位共轭矩阵的逐项矩阵乘积直接获得或在基变换之后获得；
○基于所述第一畸变矩阵确定所述非均匀介质的物理参数的至少一个映射。
Applications Claiming Priority (3)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| FR1856720A FR3084166B1 (fr) | 2018-07-19 | 2018-07-19 | Procedes et systemes de caracterisation ultrasonore non invasive d'un milieu heterogene |
| FR1856720 | 2018-07-19 | ||
| PCT/EP2019/069160 WO2020016250A1 (fr) | 2018-07-19 | 2019-07-16 | Procédés et systèmes de caractérisation ultrasonore non invasive d'un milieu hétérogène |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN112823283A true CN112823283A (zh) | 2021-05-18 |
| CN112823283B CN112823283B (zh) | 2024-08-02 |
Family
ID=63684140
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201980061452.1A Active CN112823283B (zh) | 2018-07-19 | 2019-07-16 | 通过使用超声非侵入性地表征非均匀介质的方法和系统 |
Country Status (6)
| Country | Link |
|---|---|
| US (1) | US11346819B2 (zh) |
| EP (1) | EP3824280B1 (zh) |
| KR (1) | KR20210042907A (zh) |
| CN (1) | CN112823283B (zh) |
| FR (1) | FR3084166B1 (zh) |
| WO (1) | WO2020016250A1 (zh) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114778691A (zh) * | 2022-06-16 | 2022-07-22 | 南京航空航天大学 | 一种变阵列形式的超声导波定量化成像方法 |
| CN115619684A (zh) * | 2022-11-08 | 2023-01-17 | 清华大学 | 基于编码摄像的非侵入式激光扫描的显微散射成像方法 |
Families Citing this family (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111494816A (zh) * | 2020-02-28 | 2020-08-07 | 南北兄弟药业投资有限公司 | 一种超声精准自适应聚焦系统和方法 |
| FR3108984B1 (fr) * | 2020-04-02 | 2022-12-23 | Safran | Contrôle non destructif d’une pièce mécanique en un matériau polycristallin |
| FR3114157B1 (fr) | 2020-09-15 | 2022-07-29 | Supersonic Imagine | Procédé et système de caractérisation ultrasonore d’un milieu |
| FR3114156A1 (fr) * | 2020-09-15 | 2022-03-18 | Supersonic Imagine | Procédé et système de caractérisation ultrasonore d’un milieu |
| FR3114155B1 (fr) * | 2020-09-15 | 2022-07-29 | Supersonic Imagine | Procédé et système de caractérisation ultrasonore d’un milieu |
| FR3114159A1 (fr) * | 2020-09-15 | 2022-03-18 | Supersonic Imagine | Procédé et système de caractérisation ultrasonore d’un milieu |
| FR3114158B1 (fr) * | 2020-09-15 | 2022-07-29 | Supersonic Imagine | Procédé et système de caractérisation ultrasonore d’un milieu |
| US11637043B2 (en) * | 2020-11-03 | 2023-04-25 | Applied Materials, Inc. | Analyzing in-plane distortion |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2002032316A1 (fr) * | 2000-10-20 | 2002-04-25 | Centre National De La Recherche Scientifique - Cnrs | Procede et dispositif non invasif de focalisation d'ondes acoustiques |
| WO2005008280A1 (en) * | 2003-07-17 | 2005-01-27 | Angelsen Bjoern A J | Corrections for wavefront aberrations in ultrasound imaging |
| CN101631591A (zh) * | 2007-02-21 | 2010-01-20 | 超声成像公司 | 一种用于优化经过引入像差的元件的波的聚焦的方法 |
| CN102576527A (zh) * | 2009-09-03 | 2012-07-11 | 皇家飞利浦电子股份有限公司 | 对经颅超声畸变的基于对侧阵列的校正 |
| CN107003403A (zh) * | 2014-09-26 | 2017-08-01 | 国家科学研究中心 | 声成像的方法和设备 |
Family Cites Families (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP4564286B2 (ja) * | 2004-06-14 | 2010-10-20 | 株式会社東芝 | ３次元超音波画像化装置 |
| JP5373308B2 (ja) * | 2008-03-31 | 2013-12-18 | 富士フイルム株式会社 | 超音波撮像装置及び超音波撮像方法 |
| JP6168802B2 (ja) * | 2012-04-12 | 2017-07-26 | キヤノン株式会社 | 処理装置、処理方法、及びプログラム |
-
2018
- 2018-07-19 FR FR1856720A patent/FR3084166B1/fr active Active
-
2019
- 2019-07-16 EP EP19745574.4A patent/EP3824280B1/fr active Active
- 2019-07-16 KR KR1020217003587A patent/KR20210042907A/ko not_active Application Discontinuation
- 2019-07-16 CN CN201980061452.1A patent/CN112823283B/zh active Active
- 2019-07-16 WO PCT/EP2019/069160 patent/WO2020016250A1/fr active Application Filing
- 2019-07-16 US US17/261,095 patent/US11346819B2/en active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2002032316A1 (fr) * | 2000-10-20 | 2002-04-25 | Centre National De La Recherche Scientifique - Cnrs | Procede et dispositif non invasif de focalisation d'ondes acoustiques |
| WO2005008280A1 (en) * | 2003-07-17 | 2005-01-27 | Angelsen Bjoern A J | Corrections for wavefront aberrations in ultrasound imaging |
| CN101631591A (zh) * | 2007-02-21 | 2010-01-20 | 超声成像公司 | 一种用于优化经过引入像差的元件的波的聚焦的方法 |
| CN102576527A (zh) * | 2009-09-03 | 2012-07-11 | 皇家飞利浦电子股份有限公司 | 对经颅超声畸变的基于对侧阵列的校正 |
| CN107003403A (zh) * | 2014-09-26 | 2017-08-01 | 国家科学研究中心 | 声成像的方法和设备 |
Non-Patent Citations (3)
| Title |
|---|
| AUBRY, ALEXANDRE等: "Detection and imaging in a random medium: A matrix method to overcome multiple scattering and aberration", 《JOURNAL OF APPLIED PHYSICS》, vol. 106, no. 4, XP012123911, DOI: 10.1063/1.3200962 * |
| JEAN-LUC ROBERT等: "Green’s function estimation in speckle using the decomposition of the time reversal operator: Application to aberration correction in medical imaging", 《THE JOURNAL OF THE ACOUSTICAL SOCIETY OF AMERICA》, vol. 123, no. 2, pages 866 - 877, XP055569717, DOI: 10.1121/1.2816562 * |
| 朱新杰;邓明晰;都东;韩赞东;: "接收孔径对超声导波合成孔径阵列成像检测的影响分析", 机械工程学报, no. 12, pages 149 - 156 * |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114778691A (zh) * | 2022-06-16 | 2022-07-22 | 南京航空航天大学 | 一种变阵列形式的超声导波定量化成像方法 |
| CN114778691B (zh) * | 2022-06-16 | 2022-09-06 | 南京航空航天大学 | 一种变阵列形式的超声导波定量化成像方法 |
| CN115619684A (zh) * | 2022-11-08 | 2023-01-17 | 清华大学 | 基于编码摄像的非侵入式激光扫描的显微散射成像方法 |
| CN115619684B (zh) * | 2022-11-08 | 2023-09-05 | 清华大学 | 基于编码摄像的非侵入式激光扫描的显微散射成像方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| FR3084166A1 (fr) | 2020-01-24 |
| EP3824280A1 (fr) | 2021-05-26 |
| WO2020016250A1 (fr) | 2020-01-23 |
| US11346819B2 (en) | 2022-05-31 |
| KR20210042907A (ko) | 2021-04-20 |
| FR3084166B1 (fr) | 2020-10-16 |
| US20220003721A1 (en) | 2022-01-06 |
| CN112823283B (zh) | 2024-08-02 |
| EP3824280B1 (fr) | 2022-08-03 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN112823283B (zh) | 通过使用超声非侵入性地表征非均匀介质的方法和系统 | |
| JP6408297B2 (ja) | ビームフォーミング方法、計測イメージング装置、及び、通信装置 | |
| KR101942595B1 (ko) | 영상 획득속도 최적화를 구비한 영상장치 | |
| US11776526B2 (en) | Method and system for ultrasonic characterization of a medium | |
| US11761928B2 (en) | Method and system for ultrasonic characterization of a medium | |
| US11768181B2 (en) | Method and system for ultrasonic characterization of a medium | |
| CN113424073B (zh) | 材料非线性体积弹性的超声估算 | |
| Szasz | Advanced beamforming techniques in ultrasound imaging and the associated inverse problems | |
| US20220082693A1 (en) | Method and system for ultrasonic characterization of a medium | |
| US20240036004A1 (en) | Method and system for ultrasonic characterization of a medium | |
| Montero | New developments on quantitative imaging using ultrasonic waves | |
| JP7515563B2 (ja) | 超音波を使用して不均一な媒体を非侵襲的に特性化するための方法及びシステム |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |