CN116011172B - 一种基于声散射理论的启发式校正数值声场误差的方法及其存储介质 - Google Patents
一种基于声散射理论的启发式校正数值声场误差的方法及其存储介质 Download PDFInfo
- Publication number
- CN116011172B CN116011172B CN202211484315.9A CN202211484315A CN116011172B CN 116011172 B CN116011172 B CN 116011172B CN 202211484315 A CN202211484315 A CN 202211484315A CN 116011172 B CN116011172 B CN 116011172B
- Authority
- CN
- China
- Prior art keywords
- solution
- numerical
- sound
- parameters
- acoustic
- 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.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 56
- 238000003860 storage Methods 0.000 title claims description 4
- 238000004364 calculation method Methods 0.000 claims abstract description 47
- 230000002068 genetic effect Effects 0.000 claims abstract description 31
- 210000000349 chromosome Anatomy 0.000 claims abstract description 30
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 25
- 238000004458 analytical method Methods 0.000 claims abstract description 21
- 230000035772 mutation Effects 0.000 claims abstract description 16
- 238000004088 simulation Methods 0.000 claims description 40
- 230000005428 wave function Effects 0.000 claims description 17
- 230000008569 process Effects 0.000 claims description 11
- 238000001514 detection method Methods 0.000 claims description 10
- 239000000463 material Substances 0.000 claims description 7
- 238000006243 chemical reaction Methods 0.000 claims description 6
- 238000010606 normalization Methods 0.000 claims description 5
- 238000012937 correction Methods 0.000 claims description 4
- 238000005070 sampling Methods 0.000 claims description 4
- 238000010187 selection method Methods 0.000 claims description 4
- 230000005540 biological transmission Effects 0.000 claims description 3
- 230000005855 radiation Effects 0.000 claims description 3
- 238000002474 experimental method Methods 0.000 description 21
- 238000012795 verification Methods 0.000 description 6
- 238000013461 design Methods 0.000 description 4
- 238000005259 measurement Methods 0.000 description 4
- 238000002604 ultrasonography Methods 0.000 description 4
- 230000007547 defect Effects 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 3
- 230000002159 abnormal effect Effects 0.000 description 2
- 230000005856 abnormality Effects 0.000 description 2
- 238000003384 imaging method Methods 0.000 description 2
- 230000006872 improvement Effects 0.000 description 2
- 230000014759 maintenance of location Effects 0.000 description 2
- 230000010355 oscillation Effects 0.000 description 2
- 206010006187 Breast cancer Diseases 0.000 description 1
- 208000026310 Breast neoplasm Diseases 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 238000010835 comparative analysis Methods 0.000 description 1
- 239000013078 crystal Substances 0.000 description 1
- 238000005520 cutting process Methods 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000002059 diagnostic imaging Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 239000000284 extract Substances 0.000 description 1
- 238000012905 input function Methods 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 230000005865 ionizing radiation Effects 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 238000012216 screening Methods 0.000 description 1
- 238000004611 spectroscopical analysis Methods 0.000 description 1
- 230000009897 systematic effect Effects 0.000 description 1
- 210000001519 tissue Anatomy 0.000 description 1
Classifications
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation
Landscapes
- Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
Abstract
一种基于声散射理论的启发式校正数值声场误差的方法，包括：确定声波传播数值模型的超参数及取值范围；设置遗传算法参数，并初始化染色体信息得到初始种群；建立每个染色体所对应的声波传播的数值模型，通过计算得到所述模型的数值解；求解二维亥姆霍兹方程实现对应声散射问题的解析解计算；将得到的声波传播的数值解与解析解进行比较，进行振幅标准化以及相对平均误差分析，获得种群中所有染色体的适应度、最大适应度以及最大适应度值对应的染色体；判断遗传算法是否达到终止条件，若是，则停止迭代，输出参数配置结果；若否，进行选择、交叉、变异等遗传操作产生新一代种群，返回步骤3继续迭代，直至结果满足迭代停止条件。
Description
技术领域
本发明属于超声仿真技术，具体涉及一种数值声场误差的校正方法。
背景技术
超声CT作为一种新型医学影像技术，通过控制超声换能器阵列发射和接受超声波，实现断层扫描功能，具有无创无痛、无电离辐射、低成本等优点，在乳腺癌的早期筛查方面具有重要意义。在超声CT设备开发及应用中，超声波仿真提供了一个重要的途径来帮助我们理解超声与组织的相互作用，创新超声CT技术，并以更加高效和经济的方式进行超声换能器的设计。超声波仿真对实际声波传播现象进行数学建模，运用数值模型模拟实际的物理系统。声波传播数值模型的合理建立及精确求解，可以方便地比较各超声成像方法的特征并作出改进，此外还可以指导超声换能器的设计，找到最佳发射频率、探测距离等超参数，从而获得最佳的成像性能。迄今为止，各种声波传播数值模型的建立及求解方法已经被提出，一些比较常见的方法有伪谱法、有限差分法、有限元法等，同时基于这些方法开发了各种各样超声仿真软件，例如k-Wave、Field II、OnScale、FOCUS、Simsonic、Comsol等。
在利用超声仿真软件进行声波传播模拟的过程中，需要将物理模型离散为若干个计算网格，把空间和时间上连续分布的物理量用有限个离散节点上的值来表示，方便计算机进行仿真计算。因此需要对真实的超声系统进行抽象、假设、简化以及超参数配置。超参数与普通参数不同，普通参数是在对数值模型进行求解计算过程中生成的，其值可从仿真数据中获得，而超参数的值无法从数据中获取，一般根据已有的经验确定，是在对数值模型进行求解前配置的参数。在超声仿真中，超参数包括硬件系统可调控参数如发射信号中心频率、换能器探测距离等，以及仿真设置参数如时间步长、网格数、网格尺寸等其他手动设置的参数。在使用超声仿真软件时，如何高效准确的选取仿真超参数，并保证仿真结果即数值解的准确性和鲁棒性一直是一个难题。
现有技术往往通过网格无关性验证来确定网格参数，即逐步细化网格，比较不同网格数量下的计算结果，选取与网格尺寸无关解的最大网格尺寸作为网格划分方案。对于大规模声波传播数值计算而言，这种方法需要手动加密网格并且计算耗时较长，且网格无关解并不一定代表结果准确符合实际，存在一定缺陷。现有相关技术方案(公开号:CN113671044A)利用Comsol软件，通过设计超声换能器仿真正交实验，对换能器中心频率、阵元宽度、阵元数及阵元间隙四个因素的不同取值进行组合，然后采用这些参数组合进行仿真，根据不同组合对声场的影响来选取换能器最佳检测参数。但是，该技术方案一方面需要结合软件使用者的经验来配置仿真网格、边界条件等参数，对使用者的要求较高，另一方面仿真正交实验所选取的最优检测参数只能是仿真时所用的其中一种参数组合，最终优选结果不会超越实验所用组合的范围，具有一定的局限性，且该方案在仿真结束后设计对应的物理实验来验证仿真结果的准确性，存在工作量大、效率低下的问题。
现有技术的缺点可以归纳为：
1)通过网格无关性验证来确定网格参数不能保证仿真结果符合物理实际。
理论上来说，仿真网格划分的越细密，数值解的精度越高。但在实际应用中，网格尺寸越大与实际物理结果偏差越大，尺寸过小，计算时间成本急剧增加但是精度提升有限。因此，实际应用中应选择满足计算精度的网格，用尽可能低的计算成本获取更精确的结果，最大化的利用计算资源。现有技术往往通过网格无关性验证来确定网格参数，即逐步细化网格，比较不同网格数量下的计算结果，当网格疏密对结果影响不大时，该结果即为网格无关解，选取与网格参数无关解的最大网格尺寸作为网格划分方案。对于大规模声波传播数值计算而言，这种方法需要手动加密网格并且计算耗时较长，且网格无关解并不一定代表仿真结果符合物理实际，存在一定缺陷。
2)基于物理实验的验证方法存在工作量大、效率低、成本高问题。
现有技术方案在验证仿真结果是否准确，一般通过物理实验与超声仿真相结合的方式来量化仿真结果的误差。验证过程通常需要将超声仿真模型的计算结果与对应的实际物理实验的结果相对照，并基于仿真与实验结果是否一致来判断结果的准确性。进行超声实验通常需要准备硬件设备及实验耗材，搭建物理实验平台，而新型超声换能器的制造成本昂贵，同时需要进行大量的实验测量，因此该方法存在工作量大、效率低、成本高的问题。
3)超参数的选择强烈依赖于建模方法的选择、物理场景以及先验信息。
实现高效准确的选取超参数并保证仿真结果的准确可靠，严重依赖于建模方法的选择、物理场景以及先验信息。在使用仿真软件建立声波传播数值模型时，需要根据超声换能器实际运行的物理情景，选择合适的建模方法，根据先验信息为模型施加对应的材料属性、边界条件、发射的超声波信号等，并综合考虑运算精度和时间，对模型进行网格划分。而且超声系统包含多种超参数，包括硬件系统可调控参数如发射信号中心频率、换能器探测距离等，以及仿真设置参数如时间步长、网格数、网格尺寸等，参数的系统性误差会导致仿真结果的不确定性，需要对换能器本身、超声知识、仿真软件相当熟悉才能够保证结果的可靠性与准确性，否则将会产生仿真数据不准确、换能器设计优化工作延缓等问题。
4)采取仿真正交实验法配置超参数具有一定的局限性。
正交实验法是利用正交表来研究多因素和多水平实验的一种设计方法，实验中对结果有影响的变量称为因素，各因素可能的取值称为水平。现有相关技术方案(公开号:CN113671044A)利用Comsol软件，通过设计超声换能器仿真正交实验，对换能器中心频率、阵元宽度、阵元数及阵元间隙四个因素的不同取值进行组合，然后采用这些不同因素的水平组合进行仿真，根据不同组合对声场的影响来配置换能器最佳检测参数。但是进行仿真正交实验所得到的最佳检测参数只能是实验中所用的某一种水平组合，最终优选结果不会超越实验所用组合的范围，具有一定的局限性。
发明内容
本发明克服现有技术的上述缺点，提供一种基于声散射理论的启发式校正数值声场误差的方法。本发明采用声波仿真软件配置并求解声散射问题的声波传播数值模型得到数值解，并基于声散射理论实现声散射问题的解析解求解，对数值解和解析解振幅标准化后再进行相对误差分析，并结合遗传算法自动确定高精度数值模型超参数的配置方案，从而实现数值声场的误差校正。
本发明解决的技术问题1：针对所述问题1)、2)，本发明旨在进行超参数的精确配置，保证数值解的准确且符合实际，同时降低工作量与成本，提升效率。本发明采用声波仿真软件配置并求解声散射问题的声波传播数值模型得到数值解，并基于声散射理论得到声散射问题的解析解，采用相对平均误差验证数值解的精确性。由于解析解与数值解计算均只需在计算机上完成，不需要额外的实验设备和实验测量，从而有效降低实际物理实验带来的高成本与高工作量。
本发明解决的技术问题2：针对所述问题3)、4)，本发明旨在降低声波传播数值模型的超参数配置对建模方法、物理场景以及先验信息的依赖，提高最优参数配置方案的搜索能力。本发明结合遗传算法实现超参数的自动化配置，通过遗传操作(选择、交叉、变异)来增加获得最优超参数配置方案的概率。
本发明提出一种基于声散射理论的启发式校正数值声场误差的方法，包括如下步骤：
步骤1：确定需要进行配置的声波传播数值模型的超参数及其取值范围；
步骤2：设置遗传算法参数，并在超参数取值范围内初始化染色体信息得到初始种群；
步骤3：采用声波仿真软件建立声散射问题中每个染色体所对应的声波传播的数值模型，包括计算网格、声波传播的材料介质属性、换能器的检测参数和位置参数等，通过计算得到所述模型的数值解，并将数值解从时域转化为频域；
步骤4：基于声散射理论求解二维亥姆霍兹(Helmholtz)方程实现对应声散射问题的解析解计算；确定解析解的截断误差范围来将解析解的精度控制在所期望范围内，将其截断为有限和的方式近似逼近解析解；
步骤5：将得到的声波传播的数值解与解析解进行比较，进行振幅标准化以及相对平均误差分析，获得种群中所有染色体的适应度、最大适应度以及最大适应度值对应的染色体；
步骤6：判断遗传算法是否达到终止条件，若是，则停止迭代，输出参数配置结果；若否，进行选择、交叉、变异等遗传操作产生新一代种群，返回步骤3继续迭代，直至结果满足迭代停止条件。
作为优选，步骤1中，声波传播数值模型的超参数包括但不限于：发射信号中心频率fc，换能器与散射体的距离d，仿真时间步长dt，仿真网格间距dx、dy，网格数Nx、Ny等。
作为优选，步骤2中，设置遗传算法参数，包括：种群大小np、最大迭代次数niter、选择率ps、交叉率pc、变异率pv。遗传算法的初始种群通过在超参数取值范围内进行随机抽样对染色体进行初始化来获得。
作为优选，步骤3中，采用声波仿真软件配置及求解声散射问题的声波传播的数值模型得到数值解，这些仿真软件包括但不限于k-Wave、Field II、OnScale、FOCUS、Simsonic、Comsol等，通过快速傅里叶变换将数值解从时域解转换为频域解。
作为优选，步骤3中，声波传播的数值模型，包括计算网格、声波传播的材料介质属性、换能器的检测参数和位置参数。
作为优选，所述声散射问题对二维空间下的圆形散射体进行求解，但本方法不仅限于二维，可以推广到三维圆柱/球形仿体。所述二维声波散射问题中，超声换能器为N个阵元组成的阵列，阵元通过电声能量转换实现超声波的发射和接收，发射阵元产生的声波通过圆形散射体(记为B)并产生散射波，随后这些散射波被接收阵元所接收。为便于计算，阵元被简化为点阵元，其在二维空间中的位置用坐标向量zi,(i＝1,...,N)来表示。所述超声换能器除直线形外、还可以是曲线、环形等各种形状换能器，同时可以位于二维空间中的任意位置。圆形散射体B是声速为c、半径为rB的均质介质。二维空间中的介质声速为c0，该空间的坐标原点o位于圆形散射体的圆心。
作为优选，步骤4中，所述解析解的计算包括解析解计算模型的建立以及解析解截断误差范围的确定。根据二维声散射问题的解析解的数学性质，该问题可以用亥姆霍兹(Helmholtz)方程来描述，数学表达式如下：
其中z表示二维空间中某一点的位置坐标向量，ω表示角频率，Φ(z,ω)表示z处的全波函数，K表示波数，Q(z,ω)表示源项。
方程(1)满足远场条件与辐射条件：
求解该方程可得所述声散射问题的解析解。在解析解的计算过程中，全波函数Φ(z)可以分为入射波函数Φinc和散射波函数Φsct，散射体B内的一点z∈B的全波函数Φ(z)与散射体B外的一点的散射波函数Φsct(z)为级数，通过确定解析解的截断误差范围来将解析解的精度控制在所期望范围内，将级数截断为有限和的方式近似逼近解析解。
作为优选，步骤5中，所述振幅标准化操作首先利用傅里叶级数分别对数值解与解析解结果进行拟合，然后基于其均值和标准差对拟合后的结果进行标准化(standardization)。将标准化后的数值解与解析解进行比较，计算二者的相对平均误差(Mean Relative Error,MRE)，具体为：
其中，Φ(i) e表示第i个阵元位置zi处的解析解，Φ(i) n表示第i个阵元位置zi处的数值解，N为总的阵元数，下标e和n分别表示解析解(exact solution)与数值解(numericalsolution)。所述适应度函数为相对平均误差的倒数，具体为：
作为优选，步骤6中，所述遗传算法采用二值码串来表示每个超参数的值，对二值码串进行选择、交叉、变异操作，具体为：采用轮盘赌选择法根据选择率选择部分个体，利用精英保留策略，将上一迭代过程中种群的最优个体保留至下一代；采用单点交叉根据交叉概率选择个体进行交叉操作，生成新个体；采用基本位变异根据变异概率选择个体进行变异操作，生成新个体；最终所有新生成个体组成新一代种群。
一种计算机可读存储介质，其特征在于，其上存储有程序，该程序被处理器执行时，实现一种基于声散射理论的启发式校正数值声场误差的方法。
本发明的关键点是：
1.本发明提出一种基于声散射理论的启发式校正数值声场误差的方法，该方法不局限于高精度数值声场网格参数的配置，同时也能够满足换能器中心频率、换能器与散射体的位置、网格参数、仿真时间步、边界条件其他超参数的配置需求；
2.本发明不仅可以满足声波传播数值模型的单个超参数的配置需求，也可以同时对多个超参数进行配置；
3.采用超声仿真软件建立二维圆形超声散射问题的计算模型，超声仿真软件不局限于k-Wave这一种，可以采用其他仿真软件如：Field II、OnScale、FOCUS、Simsonic、Comsol等；
4.声波传播的数值模型建模过程中，对超声换能器的复杂外形和空间位置没有限制，所述超声换能器除直线形外、还可以是曲线、环形等各种形状换能器，同时可以位于二维空间中的任意位置；
5.本发明采用二维圆形声散射体进行计算，但是本方法不仅限于二维，可以推广到三维圆柱/球形散射体。
6.本发明基于声散射理论计算解析解，所采取的解析解级数截断方案能够令解析解准确的收敛于理论结果；
7.本发明基于声散射理论，通过计算解析解与数值解之间的相对平均误差指标，能够评价仿真结果的准确性；
8.本发明在对数值解和解析解计算结果进行处理时，采用了傅里叶级数拟合数据，然后对结果进行标准化，剔除其中的异常点，避免了用单个数据进行缩放时的数据异常造成的误差；
9.本发明结合遗传算法，实现超参数的自动化配置，降低声波传播数值模型的超参数配置对仿真软件使用者经验的依赖性；
10.由于解析解与数值解计算均只需在计算机上完成，不需要额外的实验设备和实验测量，从而有效降低实际物理实验带来的高成本与高工作量。
本发明的优点是：
1.本发明提出一种基于声散射理论的启发式校正数值声场误差的方法，基于声散射理论实现数值解与解析解的对比分析，能够确定保证计算准确性的网格划分方案。且该方法不局限于网格的配置，同时也能够满足其他超参数的配置需求；
2.本发明结合遗传算法实现超参数的自动化配置，降低声波传播数值模型的超参数配置对仿真软件使用者经验的依赖性；
3.本发明结合遗传算法，通过遗传操作(选择、交叉、变异)来增加获得最优超参数配置方案的概率；
4.本发明基于声散射理论，通过计算解析解与数值解之间的相对平均误差指标，避免了使用仿真软件进行交叉验证存在的无法评估数值解与真实值之间的差距的问题；
5.由于解析解与数值解计算均只需在计算机上完成，不需要额外的实验设备和实验测量，从而有效降低实际物理实验带来的高成本与高工作量。
附图说明
图1是本发明的二维圆形超声散射问题示意图。
图2是本发明的|Φn(z)|与|Φn sct(z)|随n的变化曲线图。
图3是本发明的|Φn(z)|＝δ1时，n与|z|之间的关系示意图。
图4是本发明方法流程示意图。
图5是本发明实施例中声传播数值模型建立流程图。
具体实施方式
本发明详细实现流程如附图4所示，下面结合具体实施例对本发明进行详细说明。
步骤1：确定需要进行配置的声波传播数值模型的超参数及取值范围；在本实施例中，超参数为发射信号中心频率fc，取值范围[a,b]
步骤2：设置遗传算法参数种群大小np、最大迭代次数niter、选择率ps、交叉率pc、变异率pv，并在超参数fc的取值范围[a,b]内采用随机抽样法初始化染色体信息得到初始种群
步骤3：根据上一步得到的染色体参数，采用k-Wave调用4个输入函数包括计算网格Kgrid、材料介质Medium、声源项的位置和特征Source、传感器位置和特征Sensor来建立声散射问题的声波传播的数值模型；通过调用k-Wave的仿真函数之一kspaceFirstOrder2D进行计算得到所述模型的数值解Φ(z,t)n，该结果为时域结果，通过快速傅里叶变换将数值解从时域解Φ(zi,t)n转换为频域解Φ(zi,ω)n，并提取中心角频率处的幅值|Φ(zi,ωc)n|；
步骤4：根据步骤3建立的数值模型，建立解析解的计算模型，所需参数包括：超声换能器的阵元数N、发射信号中心频率fc、散射体半径rB、背景声速c0、散射体声速c、换能器各阵元的位置坐标向量zi,(i＝1,...,N)等。确定解析解的级数截断项n之后，完成解析解Φ(zi,ωc)e的计算；
步骤5：将得到的声波传播的数值解与解析解进行比较，首先对数值解与解析解进行振幅标准化然后再进行相对平均误差分析，并根据获得种群中所有染色体的适应度/>并找到F中的最大适应度Fitnessbest＝max(F)以及最大适应度值对应的染色体fcbest，max(·)表示为取最大值操作函数；
步骤6：判断迭代次数是否达到niter，若是，则停止迭代，输出参数最优结果fcbest；若否，采用二进制编码将每个染色体的值用二值码串来表示，并根据选择率ps，采用轮盘赌选择法选择部分个体，同时采用精英保留策略，将种群中的最优个体保留至下一代，然后采用单点交叉法，根据交叉率pc，选择个体进行交叉操作，生成新个体，采用基本位变异法，根据变异率pv，选择个体进行变异操作，生成新个体，最后将经过选择、交叉、变异等遗传操作产生的个体进行合并生成新一代种群，迭代次数+1，返回步骤3继续迭代，直至迭代次数达到niter，输出最优结果fcbest，即为最优超参数配置方案。
技术方案实现原理是：
本发明采用二维声波散射问题的解析解与数值解进行相对误差分析，结合遗传算法，达到声波传播数值模型的超参数自动配置的目的。具体技术方案实现原理及技术手段如下：
(1)二维声波散射问题的定义。
如附图1所示，在二维声波散射问题中，超声换能器为N个阵元组成的压电晶体阵列，阵列上的阵元通过压电效应实现电声能量转换进而实现超声波的发射和接收，发射阵元产生的声波通过圆形散射体(记为B)并产生散射波，随后这些散射波被接收阵元所接收。在该问题中，阵元被简化为点阵元，其在二维空间中的位置用坐标向量zi,(i＝1,...,N)来表示。所述超声换能器除图1中所示直线形外、还可以是曲线、环形等各种形状换能器，同时可以位于二维空间中的任意位置。圆形散射体B是声速为c、半径为rB的均质介质。二维空间中的介质声速为c0，该空间的坐标原点o位于圆形散射体的圆心。
该问题可以用亥姆霍兹(Helmholtz)方程来描述，数学表达式如下：
其中z表示二维空间中某一点的位置坐标向量，ω表示角频率，Φ(z,ω)表示z处的全波函数，K表示波数，Q(z,ω)表示源项。
方程(1)满足远场条件与辐射条件：
圆形散射体B的超声波场对于上述方程有一个解析解Φ(z,ω)。因此可以通过比较换能器中心频率下所有阵元中心处(坐标向量记为zi,i＝1,...,N)的解析解与数值解，当两者的数值误差满足要求时，则表明仿真参数的满足对精确性的要求。
(2)二维声波散射问题的解析解。
设置解析解的计算所需参数，这些参数包括超声换能器的阵元数N、发射信号中心频率fc、散射体半径rB、背景声速c0、散射体声速c、换能器各阵元的位置坐标向量zi,(i＝1,...,N)等。
下面介绍解析解Φ(z,ω)的计算，其中ω＝2πfc，Φ(z,ω)简写为Φ(z)：
在角频率为ω的二维点波源场中，圆形散射体B的波数为K，二维空间介质的波数为K0，波数K的计算公式如下
其中，fc表示中心频率，c表示声速；
对于散射体B外的计算域中的一点全波函数Φ(z)可以分为入射波函数Φinc和散射波函数Φsct，具体计算公式如下：
Φ(z)＝Φinc(z)+Φsct(z) (7)
对于散射体B内的计算域中的一点z∈B，全波函数用如下公式表示
其中αn、βn为与散射体边界条件相关的系数，根据散射体B的边界内外具有相同的压力和压力梯度确定，ψ为入射波的传播角，θ为散射波的传播角。
其中，对于球面入射波，系数γn由下式确定：
z0为入射点声源位置坐标向量，i为虚数单位，e为自然常数。
Jn、分别记为第一类贝塞尔(Bessel)函数与汉克尔(Hankel)函数。
(3)解析解的截断误差范围的确定。
在二维声波散射问题的解析解的计算过程中，散射体B内的一点z∈B全波函数Φ(z)与散射体B外的一点散射波函数Φsct(z)为级数，需要通过无穷项求和的过程才能得到的结果，而计算机只能完成有限次运算，因此需通过将级数截断为有限和的方式近似逼近解析解。这样就产生了有限过程代替无限过程的误差，即截断误差。我们需要确定解析解的截断误差范围来控制解析解的精度。对于二维空间内一点，记第n项为Φn(z)、Φn sct(z)，具体表达式如下：
Φn(z)＝βnJn(K|z|)ein(θ-ψ),z∈B (14)
由于Φn(z)、Φn sct(z)均为复函数，我们对其取模以便于观察，如附图2所示|Φn(z)|、|Φn sct(z)|的值在一定范围内存在震荡现象，超出该震荡范围后逐渐衰减至0，为了控制截断误差范围，需要控制截断项在衰减区范围内，例如，取截断项n使|Φn(z)|＜δ1、|Φn sct(z)|＜δ2，其中δ1、δ2为误差范围，考虑实际问题对计算精度的需求进行设置，取值越小精度越高，计算也成本更高，一般δ1、δ2取10-8已足够精确；
具体的，在散射体径向方向取一条观测线，这条线的长度为散射体中心点到二维空间边界的距离，然后将这条线分成若干等间距的点，对于散射体B内部和外部的点，计算每个点处的|Φn(z)|、|Φn sct(z)|的值随n的变化，这里以为例|Φn(z)|令|Φn(z)|＝δ1，得到n和z的值，通过直线拟合的方式获得n与z之间的关系，如附图3所示给出了5MHz频率下的拟合结果，取截断项n在该直线上方范围内控制解析解的精度在期望范围内，即n＞a|z|+b，满足|Φn(z)|＜δ1，其中a为直线斜率，b为截距。|Φn sct(z)|的处理方式同|Φn(z)|，确定级数截断项n之后，根据(2)中给出的解析解计算公式即可完成解析解Φ(zi,ωc)e的计算，i表示第i个阵元(i＝1,...,N)，ωc＝2πfc为中心角频率；
(4)声波传播数值模型建立及求解。
采用超声仿真软件(如k-Wave、Field II、OnScale、FOCUS、Simsonic、Comsol等)建立二维声波散射问题的数值模型，超声换能器的阵元数N、发射信号中心频率fc、散射体半径rB、背景声速c0、散射体声速c、换能器各阵元的位置坐标向量zi,(i＝1...N)等与解析解计算模型保持一致，同时，需要输入时间步长dt、二维空间的网格间距dx、dy、网格数Nx、Ny等仿真参数(与所采用的仿真软件相关，不同仿真软件所需的仿真参数有所区别)，根据这些参数定义仿真计算网格、声传播介质材料参数、超声换能器的位置和检测参数等，通过计算得到二维声波散射问题的数值解。
(5)将数值解从时域数据转换为频域数据。
由于数值解多为时域信号Φ(z,t)n，而解析解得到的是频域信号Φ(z,ωc)e，因此，需要通过快速傅里叶变换时域数值解Φ(zi,t)n转换为频域结果Φ(zi,ω)n，具体为
Φ(zi,ω)n＝FFT(Φ(zi,t)n),(i＝1,...,N) (16)
其中，FFT()表示快速傅里叶变换函数，经过FFT变换后，提取中心角频率处的幅值|Φ(zi,ωc)n|。
(6)振幅标准化。
由于数值解与解析解的振幅不同，因此需要对振幅进行标准化处理，以便下一步进行相对平均误差的计算。为了剔除了其中的异常点，避免了用单个数据进行缩放时的数据异常造成的误差，提高计算的准确性，首先利用傅里叶级数分别对数值解与解析解结果进行拟合，然后对二者拟合后的结果进行标准化(standardization)后得到Φ(i) n，Φ(i) e，i＝1,...,N，即对拟合过后的换能器所有阵元的声场数据的进行幅值标准化，以便下一步进行进行相对误差分析。
为了实现对振幅的标准化处理，首先采用傅里叶拟合数值解与解析解结果，将数值解与解析解看作阵元数与幅值的关系。对于拟合函数选取有两个重要的考量：拟合精确度尽可能高、拟合函数的表达要尽可能简洁。经过综合考量，八阶傅里叶级数可以很好的满足这两个要求，拟合公式为
式中，a0、an、bn、w为待定系数，i取值为i＝1,...,N对应第1,...,N个阵元，f(i)为拟合后的结果，对应仿真结果|Φ(zi,ωc)n|,(i＝1,...,N)或解析解结果|Φ(zi,ωc)e|,(i＝1,...,N)。
对拟合后的结果进行标准化，公式如下：
下标n,e分别表示数值解与解析解结果，μ为f(i),(i＝1,...,N)的平均值，σ为f(i)的标准差。
(7)相对平均误差计算。
将标准化后的数值解与解析解进行比较，计算二者的相对平均误差，具体公式如下：
其中，Φ(i) e表示第i个阵元位置zi处的解析解，Φ(i) n表示第i个阵元位置zi处的数值解，N为总的阵元数。
(8)采用遗传算法实现声波传播数值模型的超参数自动配置。
采用遗传算法实现声波传播数值模型的超参数的自动配置，遗传算法模拟自然界生物进化机制，首先随机产生一组初始解(称为种群)，种群中的个体代表一个解，称为染色体，每个染色体的优劣用适应度来评价，这些染色体在后续迭代过程中不断进化。新一代种群形成过程中，通过选择操作将父代中高适应度值的个体以更高的概率遗传到子代，通过交叉和变异操作，保证子代的多样性。经过若干次迭代后，算法收敛于最优染色体，即最优解。
采用遗传算法实现声波传播数值模型的超参数配置的具体流程如下：
1)确定需要进行配置的声波传播数值模型的超参数及取值范围；
2)设置遗传算法参数，包括：种群大小np、最大迭代次数niter、选择率ps、交叉率pc、变异率pv。通过在超参数取值范围内进行随机抽样对染色体(超参数)进行初始化来获得大小为np的初始种群；
3)将上一步得到的染色体参数通过声波仿真软件建立每个染色体对应的声波传播数值模型，求解得到数值解；
4)基于声散射理论计算对应染色体的解析解。
5)将3)中所得数值解与解析解进行比较，进行振幅标准化以及相对平均误差分析，获得种群中所有染色体的适应度、最大适应度以及最大适应度值对应的染色体；适应度函数为相对平均误差的倒数，具体为：
6)判断迭代次数是否达到niter，若是，则停止迭代，输出最优结果；若否，进行下一步；
7)采用二进制编码将每个染色体的值用二值码串来表示，根据选择率ps，采用轮盘赌选择法选择部分个体，同时采用精英保留策略，将种群中的最优个体保留至下一代；
8)采用单点交叉法，根据交叉率pc，选择个体进行交叉操作，生成新个体；
9)采用基本位变异法，根据变异率pv，选择个体进行变异操作，生成新个体；
10)将经过遗传操作(选择、交叉、变异)的个体进行合并，生成新的种群，迭代次数+1，返回3)继续迭代，直至迭代次数达到niter，输出最优结果。
Claims (6)
1.一种基于声散射理论的启发式校正数值声场误差的方法，包括如下步骤：
步骤1：确定需要进行配置的声波传播数值模型的超参数及取值范围；
步骤2：设置遗传算法参数，并在超参数取值范围内初始化染色体信息得到初始种群；
步骤3：采用声波仿真软件建立声散射问题中每个染色体所对应的声波传播的数值模型，包括计算网格、声波传播的材料介质属性、换能器的检测参数和位置参数，声波仿真软件包括k-Wave、Field II、OnScale、FOCUS、Simsonic、Comsol，通过计算得到所述模型的数值解，并通过快速傅里叶变换将数值解从时域转化为频域；具体步骤包括：对二维空间下的圆形散射体进行求解；二维声波散射问题中，超声换能器为N个阵元组成的阵列，阵元通过电声能量转换实现超声波的发射和接收，发射阵元产生的声波通过圆形散射体B并产生散射波，随后这些散射波被接收阵元所接收；为便于计算，阵元被简化为点阵元，其在二维空间中的位置用坐标向量zi来表示，其中i＝1,...,N；所述超声换能器包括直线形、曲线、环形形状换能器，同时位于二维空间中的任意位置；圆形散射体B是声速为c、半径为rB的均质介质；二维空间中的介质声速为c0，该空间的坐标原点o位于圆形散射体的圆心；
步骤4：基于声散射理论求解二维亥姆霍兹(Helmholtz)方程实现对应声散射问题的解析解计算；确定解析解的截断误差范围来将解析解的精度控制在所期望范围内，将其截断为有限和的方式近似逼近解析解；所述解析解的计算包括解析解计算模型的建立以及解析解截断误差范围的确定；根据二维声散射问题的解析解的数学性质，该问题用亥姆霍兹(Helmholtz)方程来描述，数学表达式如下：
其中z表示二维空间中某一点的位置坐标向量，ω表示角频率，Φ(z,ω)表示z处的全波函数，K表示波数，Q(z,ω)表示源项；
方程(1)满足远场条件与辐射条件：
求解该方程可得所述声散射问题的解析解；在解析解的计算过程中，全波函数Φ(z)分为入射波函数Φinc和散射波函数Φsct，散射体B内的一点z∈B全波函数Φ(z)与散射体B外的一点散射波函数Φsct为级数，通过确定解析解的截断误差范围来将解析解的精度控制在所期望范围内，将级数截断为有限和的方式近似逼近解析解；
步骤5：将得到的声波传播的数值解与解析解进行比较，进行振幅标准化以及相对平均误差分析，获得种群中所有染色体的适应度、最大适应度以及最大适应度值对应的染色体；
步骤6：判断遗传算法是否达到终止条件，若是，则停止迭代，输出参数配置结果；若否，进行选择、交叉、变异的遗传操作产生新一代种群，返回步骤3继续迭代，直至结果满足迭代停止条件。
2.如权利要求1所述的一种基于声散射理论的启发式校正数值声场误差的方法，其特征在于，步骤1中，声波传播数值模型的超参数包括：发射信号中心频率fc，换能器与散射体的距离d，仿真时间步长dt，仿真网格间距dx、dy，网格数Nx、Ny。
3.如权利要求1所述的一种基于声散射理论的启发式校正数值声场误差的方法，其特征在于，步骤2中，设置遗传算法参数，包括：种群大小np、最大迭代次数niter、选择率ps、交叉率pc、变异率pv；遗传算法的初始种群通过在超参数取值范围内进行随机抽样对染色体进行初始化来获得。
4.如权利要求1所述的一种基于声散射理论的启发式校正数值声场误差的方法，其特征在于，步骤5中，所述振幅标准化操作首先利用傅里叶级数分别对数值解与解析解结果进行拟合，然后基于其均值和标准差对拟合后的结果进行标准化(standardization)；将标准化后的数值解与解析解进行比较，计算二者的相对平均误差(Mean Relative Error,MRE)，具体为：
其中，Φ(i) e表示第i个阵元位置zi处的解析解，Φ(i) n表示第i个阵元位置zi处的数值解，N为总的阵元数，下标e和n分别表示解析解(exact solution)与数值解(numericalsolution)；所述适应度函数为相对平均误差的倒数，具体为：
5.如权利要求1所述的一种基于声散射理论的启发式校正数值声场误差的方法，其特征在于，步骤6中，所述遗传算法采用二值码串来表示每个超参数的值，对二值码串进行选择、交叉、变异操作，具体为：采用轮盘赌选择法根据选择率选择部分个体，利用精英保留策略，将上一迭代过程中种群的最优个体保留至下一代；采用单点交叉根据交叉概率选择个体进行交叉操作，生成新个体；采用基本位变异根据变异概率选择个体进行变异操作，生成新个体；最终所有新生成个体组成新一代种群。
6.一种计算机可读存储介质，其特征在于，其上存储有程序，该程序被处理器执行时，实现权利要求1-5中任一项所述的一种基于声散射理论的启发式校正数值声场误差的方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211484315.9A CN116011172B (zh) | 2022-11-24 | 2022-11-24 | 一种基于声散射理论的启发式校正数值声场误差的方法及其存储介质 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211484315.9A CN116011172B (zh) | 2022-11-24 | 2022-11-24 | 一种基于声散射理论的启发式校正数值声场误差的方法及其存储介质 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116011172A CN116011172A (zh) | 2023-04-25 |
| CN116011172B true CN116011172B (zh) | 2023-12-01 |
Family
ID=86023762
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202211484315.9A Active CN116011172B (zh) | 2022-11-24 | 2022-11-24 | 一种基于声散射理论的启发式校正数值声场误差的方法及其存储介质 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116011172B (zh) |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107409256A (zh) * | 2015-02-16 | 2017-11-28 | 歌拉利旺株式会社 | 声场校正装置、声场校正方法和声场校正程序 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US9513372B2 (en) * | 2013-01-22 | 2016-12-06 | Schlumberger Technology Corporation | Automatic processing of ultrasonic data |
-
2022
- 2022-11-24 CN CN202211484315.9A patent/CN116011172B/zh active Active
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107409256A (zh) * | 2015-02-16 | 2017-11-28 | 歌拉利旺株式会社 | 声场校正装置、声场校正方法和声场校正程序 |
Non-Patent Citations (2)
| Title |
|---|
| Binaural Sound Localization Based on Detection of Multi-Band Zero-Crossing Points;LI Cong-qing等;2009 Second International Conference on Intelligent Networks and Intelligent Systems;第393-396页 * |
| 平面波经颅超声成像相位校正及散斑跟踪;李倩岩等;应用声学;第41卷(第1期);第132-142页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116011172A (zh) | 2023-04-25 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Gu et al. | Modeling of wave propagation for medical ultrasound: a review | |
| Jensen | A multi-threaded version of Field II | |
| Douglas Mast et al. | Simplified expansions for radiation from a baffled circular piston | |
| EP2519321A2 (en) | High intensity focused ultrasound transducer optimization | |
| CN115324564B (zh) | 固井质量检测方法、装置、计算设备及存储介质 | |
| Macaulay et al. | Accuracy of the Kirchhoff-approximation and Kirchhoff-ray-mode fish swimbladder acoustic scattering models | |
| CN116011172B (zh) | 一种基于声散射理论的启发式校正数值声场误差的方法及其存储介质 | |
| Flores et al. | Experimental study to evaluate the generation of reverberant shear wave fields (R-SWF) in homogenous media | |
| Kabanikhin et al. | On the problem of modeling the acoustic radiation pattern of source for the 2D first-order system of hyperbolic equations | |
| Soneson et al. | Gaussian representation of high-intensity focused ultrasound beams | |
| Jansen et al. | Enhanced radon domain beamforming using deep-learning-based plane wave compounding | |
| Hu et al. | Assessment of homodyned K distribution modeling ultrasonic speckles from scatterers with varying spatial organizations | |
| Suvorov et al. | Numerical simulation of acoustic emission using acoustic contact elements | |
| US11464463B2 (en) | Elastography based on x-ray computed tomography and sound wave integration | |
| Duran et al. | Gpu accelerated acoustic field determination for a continuously excited circular ultrasonic transducer | |
| Wang et al. | An easily-achieved time-domain beamformer for ultrafast ultrasound imaging based on compressive sensing | |
| CN114065630A (zh) | 基于遗传算法的不确定参数聚焦匹配场声源功率估计方法 | |
| CN101739503A (zh) | 一种预测光在生物组织中传播模型的实现方法 | |
| Huijssen et al. | Simulations of the nonlinear acoustic pressure field without using the parabolic approximation | |
| Le Moign et al. | Optimization of virtual sources distribution in 3D echography | |
| Huang et al. | Auto-tuning numerical method for acoustic wave simulation using analytical solution | |
| Pedersen et al. | Finite-amplitude sound propagation effects in volume backscattering measurements for fish abundance estimation | |
| Wang et al. | Finite difference-embedded UNet for solving transcranial ultrasound frequency-domain wavefield | |
| Ochmann et al. | Construction of analytical solutions for the error estimation of acoustical boundary element solvers | |
| Martellotta | Caveats and pitfalls in acoustic simulation of non-existing buildings |
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 |