一种生物腔体光声内窥成像方法及系统发明专利

一种生物腔体光声内窥成像方法及系统

摘要：本发明公开一种生物腔体光声内窥成像方法及系统，获取生物腔体的位置r处的光吸收能量初始值A0(r)和声速初始值cs,0(r)，分别计算位置r处的第k次的光吸收能量Ak(r)、第k+1次的光吸收能量Ak+1(r)、第k次的声速cs,k(r)和第k+1次的声速cs,k+1(r)；通过判断光吸收能量绝对差εA,k和声速绝对差是否满足收敛条件，决定是继续迭代还是直接构建光吸收能量分布图和声速分布图。本发明考虑生物组织的复杂性，综合考虑速度的差异性，结合光吸收能量和不同声速构建光声图像，减少由待测组织内的声速分布不均匀导致的误差，得到高质量的腔体横截面光吸收能量分布图。

专利类型：发明专利

专利申请号：CN201910806275.7

专利申请（专利权）人：华北电力大学(保定)
权利人地址：河北省保定市莲池区永华北大街619号

专利发明（设计）人： 孙正 ,段爽

主权利要求:一种生物腔体光声内窥成像方法及系统
技术领域
[0001] 本发明涉及医学成像技术领域，特别是涉及一种生物腔体光声内窥成像方法及系统。
背景技术
[0002] 生物光声内窥(Photoacousticendoscopy,PAE)成像是一种新型的医学功能成像方法，它将特制的成像导管直接插入目标腔体中，导管探头发出短脉冲激光照射周围组织，组织吸收光能量后受热膨胀，激发出超声波并向组织表面传播，即光声信号。位于成像导管顶端的超声探测器在腔内扫描周围组织的过程中接收来自各个方向的光声信号。将在不同测量位置收集的光声压时间序列送入计算机后，采用合适的图像重建算法即可反演重建出腔体横截面的光吸收能量分布图像，反映组织的形态结构和功能成分。
[0003] 考虑到生物组织的复杂性，超声波在具有不同成分的组织中传播时，其速度通常存在较大差异。而多数情况下，在进行PAE成像前通常无法预知组织中的声速分布情况，因此在重建图像的过程中，恒定声速的假设可能导致重建图像中存在严重的声学畸变、伪影、模糊以及目标错位等问题。
发明内容
[0004] 本发明的目的是提供一种生物腔体光声内窥成像方法及系统，减少由待测组织内的声速分布不均匀导致的误差，得到高质量的腔体横截面光吸收能量分布图像。
[0005] 为实现上述目的，本发明提供了如下方案：
[0006] 一种生物腔体光声内窥成像方法，所述方法包括：
[0007] 获取生物腔体的位置r处的光吸收能量初始值A0(r)和声速初始值cs,0(r)；
[0008] 分别计算生物腔体的位置r处的第k次的光吸收能量Ak(r)和第k+1次的光吸收能量Ak+1(r)；k为迭代次数，k≥1；
[0009] 分别计算生物腔体的位置r处的第k次的声速cs,k(r)和第k+1次的声速cs,k+1(r)；
[0010] 计算所述第k次的光吸收能量Ak(r)与所述第k+1次的光吸收能量Ak+1(r)的光吸收能量绝对差εA,k，以及所述第k次的声速cs,k(r)与所述第k+1次的声速cs,k+1(r)的声速绝对差[0011] 判断所述光吸收能量绝对差εA,k和所述声速绝对差 是否满足收敛条件，得到收敛判断结果；
[0012] 若所述收敛判断结果为是，则输出位置r处的第k+1次光吸收能量Ak+1(r)和声速cs,k+1(r)，分别构建光吸收能量分布图和声速分布图；
[0013] 若所述收敛判断结果为否，则更新迭代次数，返回至所述分别计算生物腔体的位置r处的第k次的光吸收能量Ak(r)和第k+1次的光吸收能量Ak+1(r)。
[0014] 可选的，根据公式Ak(r)＝Ak‑1(r)‑W‑1(Ak‑1(r))G1'(r,cs,k‑1(r),Ak‑1(r))计算所述第k次的光吸收能量Ak(r)和所述第k+1次的光吸收能量Ak+1(r)；
[0015] 其中，
[0016]
[0017] argmin[·]是使g(Ak(r))+λΦTV(Ak(r))取最小值时的Ak(r)；
[0018] W(Ak‑1(r))是所述第一光声关系函数的近似黑塞矩阵；
[0019] W‑1(Ak‑1(r))是W(Ak‑1(r))的逆矩阵；
[0020] A是光吸收能量分布矩阵；A≥0指的是光吸收能量分布矩阵中大于或等于0的光吸收能量；
[0021] g(Ak(r))＝||pm(r)‑H(cs,k‑1(r))·Ak(r)||2是第二光声关系函数；
[0022] pm(r)是位置r处光声信号的测量值；
[0023] H(cs,k‑1(r))是与声速有关的第一算子；
[0024] ||·||是2‑范数；
[0025] ΦTV(Ak(r))是TV正则化项，
[0026] λ是TV正则化参数；
[0027] Lk是关于Ak(r)方差特性的稀疏矩阵；
[0028] η＞0是常数；
[0029] G1'(r,Ak‑1(r),cs,k‑1(r))是所述一光声关系函数的梯度。
[0030] 可选的，根据公式
[0031] 分别计算所述第k次的声速cs,k(r)和所述第k+1次的声速cs,k+1(r)；
[0032] 其中，
[0033]
[0034] ▽f(qk)是f(qk)的梯度；
[0035] f(qk)＝||pm(r)‑H(qk)·Ak‑1(r)||2；
[0036]
[0037] γk‑1是第k‑1次迭代后的权重参数，γk‑2是第k‑2次迭代后的权重参数；
[0038] 第k次迭代后的权重参数γk为
[0039] H(qk)是与声速有关的第二算子；
[0040] D是f(cs,k‑1(r))导数的普希茨常数；
[0041] f(cs,k‑1(r))＝||pm(r)‑H(cs,k‑1(r))·Ak‑1(r)||2。
[0042] 可选的，所述收敛条件为εA,k＜εA且
[0043] 其中εA为光吸收能量的收敛公差， 为声速的收敛公差。
[0044] 可选的，所述构建光吸收能量分布图像和声速分布图像具体包括：
[0045] 采用如下公式对极坐标下的光吸收能量分布矩阵A进行归一和灰度处理：
[0046]
[0047] B(i,j)是对A(i,j)进行归一化和灰度化处理后的值；
[0048] i是生物腔体位置r处对应极坐标系下的横坐标，j是生物腔体位置r处对应极坐标系下的纵坐标；
[0049] A(i,j)是所述光吸收能量分布矩阵A中第i行第j列的元素值；
[0050] min(A)是所述光吸收能量分布矩阵A中元素的最小值，max(A)是所述光吸收能量分布矩阵A中元素的最大值；
[0051] 采用如下公式对极坐标下的声速分布矩阵cs进行归一和灰度处理；
[0052]
[0053] E(i,j)是对cs(i,j)进行归一化和灰度化处理后的值；
[0054] cs(i,j)是所述声速分布矩阵cs中第i行第j列的元素值；
[0055] min(cs)是所述声速分布矩阵cs中的最小值；max(cs)是所述声速分布矩阵cs中的最大值；
[0056] 根据公式B′(x,y)＝B(jcosi,jsini)将所述B(i,j)转化为平面直角坐标中的B′(x,y)；
[0057] 根据公式E′(x,y)＝E(jcosi,jsini)将所述E(i,j)转化为平面直角坐标中的E′(x,y)；
[0058] 其中，i∈[0,360]，j∈[0,d]；
[0059] d是θ‑l极坐标系中极径的最大值；
[0060] B′(x,y)是所述光吸收能量分布图中直角坐标系下点(x,y)的灰度值；E′(x,y)是所述声速分布图中直角坐标系下点(x,y)的灰度值。
[0061] 一种生物腔体光声内窥成像系统，所述系统包括：
[0062] 初始值获取模块，用于获取生物腔体的位置r处的光吸收能量初始值A0(r)和声速初始值cs,0(r)；
[0063] 光吸收能量计算模块，用于分别计算生物腔体的位置r处的第k次的光吸收能量Ak(r)和第k+1次的光吸收能量Ak+1(r)；k为迭代次数，k≥1；
[0064] 声速计算模块，用于分别计算生物腔体的位置r处的第k次的声速cs,k(r)和第k+1次的声速cs,k+1(r)；
[0065] 绝对差计算模块，用于计算所述第k次的光吸收能量Ak(r)与所述第k+1次的光吸收能量Ak+1(r)的光吸收能量绝对差εA,k，以及所述第k次的声速cs,k(r)与所述第k+1次的声速cs,k+1(r)的声速绝对差
[0066] 收敛判断结果模块，用于判断所述光吸收能量绝对差εA,k和所述声速绝对差 是否满足收敛条件，得到收敛判断结果；
[0067] 图像构建模块，与所述收敛判断结果模块连接，用于在所述收敛判断结果模块为是时，则根据位置r处的第k+1次光吸收能量Ak+1(r)和声速cs,k+1(r)分别构建光吸收能量分布图和声速分布图；
[0068] 迭代更新模块，分别与所述收敛判断结果模块和所述光吸收能量计算模块连接，用于在所述收敛判断结果模块为否时，则更新迭代次数，返回至光吸收能量计算模块。
[0069] 可选的，根据公式Ak(r)＝Ak‑1(r)‑W‑1(Ak‑1(r))G1'(r,cs,k‑1(r),Ak‑1(r))计算所述第k次的光吸收能量Ak(r)和所述第k+1次的光吸收能量Ak+1(r)；
[0070] 其中，
[0071]
[0072] argmin[·]是使g(Ak(r))+λΦTV(Ak(r))取最小值时的Ak(r)；
[0073] W(Ak‑1(r))是所述第一光声关系函数的近似黑塞矩阵；
[0074] W‑1(Ak‑1(r))是W(Ak‑1(r))的逆矩阵；
[0075] A是光吸收能量分布矩阵；A≥0指的是光吸收能量分布矩阵中大于或等于0的光吸收能量；
[0076] g(Ak(r))＝||pm(r)‑H(cs,k‑1(r))·Ak(r)||2是第二光声关系函数；
[0077] pm(r)是位置r处光声信号的测量值；
[0078] H(cs,k‑1(r))是与声速有关的第一算子；
[0079] ||·||是2‑范数；
[0080] ΦTV(Ak(r))是TV正则化项，
[0081] λ是TV正则化参数；
[0082] Lk是关于Ak(r)方差特性的稀疏矩阵；
[0083] η＞0是常数；
[0084] G1'(r,Ak‑1(r),cs,k‑1(r))是所述一光声关系函数的梯度。
[0085] 可选的，根据公式
[0086] 分别计算所述第k次的声速cs,k(r)和所述第k+1次迭代后的声速cs,k+1(r)；
[0087] 其中，
[0088]
[0089] ▽f(qk)是f(qk)的梯度；
[0090] f(qk)＝||pm(r)‑H(qk)·Ak‑1(r)||2；
[0091]
[0092] γk‑1是第k‑1次迭代后的权重参数，γk‑2是第k‑2次迭代后的权重参数；
[0093] 第k次迭代后的权重参数γk为
[0094] H(qk)是与声速有关的第二算子；
[0095] D是f(cs,k‑1(r))导数的普希茨常数；
[0096] f(cs,k‑1(r))＝||pm(r)‑H(cs,k‑1(r))·Ak‑1(r)||2。
[0097] 可选的，所述收敛条件为εA,k＜εA且
[0098] 其中εA为光吸收能量的收敛公差， 为声速的收敛公差。
[0099] 可选的，所述图像构建模块具体包括：
[0100] 光吸收能量分布矩阵处理单元，用于采用如下公式对极坐标下的光吸收能量分布矩阵A进行归一和灰度处理：
[0101]
[0102] B(i,j)是对A(i,j)进行归一化和灰度化处理后的值；
[0103] i是生物腔体位置r处对应极坐标系下的横坐标，j是生物腔体位置r处对应极坐标系下的纵坐标；
[0104] A(i,j)是所述光吸收能量分布矩阵A中第i行第j列的元素值；
[0105] min(A)是所述光吸收能量分布矩阵A中元素的最小值，max(A)是所述光吸收能量分布矩阵A中元素的最大值；
[0106] 声速分布矩阵处理单元，用于采用如下公式对极坐标下的声速分布矩阵cs进行归一和灰度处理；
[0107]
[0108] E(i,j)是对cs(i,j)进行归一化和灰度化处理后的值；
[0109] cs(i,j)是所述声速分布矩阵cs中第i行第j列的元素值；
[0110] min(cs)是所述声速分布矩阵cs中的最小值；max(cs)是所述声速分布矩阵cs中的最大值；
[0111] 光吸收能量分布图灰度值转换单元，用于根据公式B′(x,y)＝B(jcosi,jsini)将所述B(i,j)转化为平面直角坐标中的B′(x,y)；
[0112] 声速分布图灰度值转换单元，用于根据公式E′(x,y)＝E(jcosi,jsini)将所述E(i,j)转化为平面直角坐标中的E′(x,y)；
[0113] 其中，i∈[0,360]，j∈[0,d]；
[0114] d是θ‑l极坐标系中极径的最大值；
[0115] B′(x,y)是所述光吸收能量分布图中直角坐标系下点(x,y)的灰度值；E′(x,y)是所述声速分布图中直角坐标系下点(x,y)的灰度值。
[0116] 根据本发明提供的具体实施例，本发明公开了以下技术效果：本发明获取获取生物腔体的位置r处的光吸收能量初始值A0(r)和声速初始值cs,0(r)，然后分别计算位置r处的第k次的光吸收能量Ak(r)、第k+1次的光吸收能量Ak+1(r)、第k次的声速cs,k(r)和第k+1次的声速cs,k+1(r)；通过判断光吸收能量绝对差εA,k和声速绝对差 是否满足收敛条件，决定是继续迭代还是直接构建光吸收能量分布图像和声速分布图像。本发明考虑生物组织的复杂性，综合考虑速度的差异性，结合光吸收能量和不同声速重建光声图像，减少由待测组织内的声速分布不均匀导致的误差，得到高质量的腔体横截面光吸收能量分布图像。
附图说明
[0117] 为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
[0118] 图1为本发明的生物腔体光声内窥成像方法的示意图；
[0119] 图2为本发明的生物腔体光声内窥成像系统的模块结构图；
[0120] 图3是XOY直角坐标系中的生物腔体横截面示意图；
[0121] 图4是θ‑l极坐标系中的多层腔体壁组织示意图。
具体实施方式
[0122] 下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0123] 本发明的目的是提供一种生物腔体光声内窥成像方法及系统，考虑到生物组织的复杂性，在构建光吸收能量分布图像过程中结合声速不同的条件，使得构建分图像更加准确。
[0124] 为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。
[0125] 一种生物腔体光声内窥成像方法，所述方法包括：
[0126] 获取生物腔体的位置r处的光吸收能量初始值A0(r)和声速初始值cs,0(r)。
[0127] 本实施例中，光吸收能量初始值取A0(r)＝0，声速初始值取cs,0(r)＝1600m/s。
[0128] 分别计算生物腔体的位置r处的第k次的光吸收能量Ak(r)和第k+1次的光吸收能量Ak+1(r)；k为迭代次数，k≥1；
[0129] 具体的，本实施例根据公式Ak(r)＝Ak‑1(r)‑W‑1(Ak‑1(r))G1'(r,cs,k‑1(r),Ak‑1(r))计算所述第k次的光吸收能量Ak(r)和所述第k+1次的光吸收能量Ak+1(r)；
[0130] 其中，
[0131]
[0132] argmin[·]是使g(Ak(r))+λΦTV(Ak(r))取最小值时的Ak(r)；
[0133] W(Ak‑1(r))是所述第一光声关系函数的近似黑塞矩阵；
[0134] W‑1(Ak‑1(r))是W(Ak‑1(r))的逆矩阵；
[0135] A是光吸收能量分布矩阵；A≥0指的是光吸收能量分布矩阵中大于或等于0的光吸收能量；
[0136] g(Ak(r))＝||pm(r)‑H(cs,k‑1(r))·Ak(r)||2是第二光声关系函数；
[0137] pm(r)是位置r处光声信号的测量值；
[0138] H(cs,k‑1(r))是与声速有关的第一算子；
[0139]
[0140] p(r)是位置r处光声信号的理论值；
[0141] ||·||是2‑范数；
[0142] ΦTV(Ak(r))是TV正则化项，
[0143] λ是TV正则化参数；
[0144] Lk是关于Ak(r)方差特性的稀疏矩阵；
[0145] η＞0是常数；
[0146] G1'(r,Ak‑1(r),cs,k‑1(r))是所述一光声关系函数的梯度。
[0147] 在具体处理过程中，采用蒙特卡罗模拟的方法，仿真得到在短脉冲激光照射下腔体壁组织中的光能量沉积。然后，采用时域有限差分算法求解极坐标系下的离散光声波动方程，得到组织产生的光声信号的理论值：
[0148]
[0149] 其中，(i,j)是腔体横截面上的点r在θ‑l极坐标系中的坐标；Δθ和Δl分别是θ轴(n+1)和l轴上的单位长度；Δt是离散时间间距；n是离散时刻；p (i,j)是位置为(i,j)的质点在n+1时刻产生的光声信号的理论值； 和 分别是位置为(i,j)的质点在时刻
n沿θ方向和l方向的振动速度；cs(i,j)是超声波在位置(i,j)；处的传播速度；β是组织的等(n)
压膨胀系数；Cp是组织的比热容；ρ0是组织的密度；I 是激光脉冲函数在时刻n的值；A(i,j)是位置(i,j)处的光吸收能量。
[0150] 分别计算生物腔体的位置r处的第k次的声速cs,k(r)和第k+1次的声速cs,k+1(r)；
[0151] 本实施例中，
[0152]
[0153] 对第三光声关系函数进行k迭代后，得到位置r处声速为
[0154] 经过求解，
[0155]
[0156] 根据公式分别计算所述第k次的声速cs,k(r)和所述第k+1次的声速cs,k+1(r)。
[0157] 其中，
[0158]
[0159] ▽f(qk)是f(qk)的梯度；
[0160] f(qk)＝||pm(r)‑H(qk)·Ak‑1(r)||2；
[0161]
[0162] γk‑1是第k‑1次迭代后的权重参数，γk‑2是第k‑2次迭代后的权重参数；
[0163] 第k次迭代后的权重参数γk为
[0164] H(qk)是与声速有关的第二算子；
[0165] D是f(cs,k‑1(r))导数的普希茨常数；
[0166] f(cs,k‑1(r))＝||pm(r)‑H(cs,k‑1(r))·Ak‑1(r)||2；
[0167] cs声速分布矩阵，cs≥0是声速分布矩阵中大于或等于0的声速。
[0168] 计算所述第k次的光吸收能量Ak(r)与所述第k+1次的光吸收能量Ak+1(r)的光吸收能量绝对差εA,k，以及所述第k次的声速cs,k(r)与所述第k+1次的声速cs,k+1(r)的声速绝对差所述收敛条件为εA,k＜εA且 其中εA为光吸收能量的收敛公差， 为声速的收敛公差。
[0169] 本实施例中，光吸收能量的收敛公差εA取0.01，声速的收敛公差 取0.01。
[0170] 在迭代计算的过程中，如果所有迭代值均不满足收敛条件，则可以通过改变光吸收能量的收敛公差εA与声速的收敛公差 的取值，调整收敛条件，以使迭代过程完整。
[0171] 判断所述光吸收能量绝对差εA,k和所述声速绝对差 是否满足收敛条件，得到收敛判断结果；
[0172] 若所述收敛判断结果为是，则输出位置r处的第k+1次光吸收能量Ak+1(r)和声速cs,k+1(r)，分别构建光吸收能量分布图和声速分布图；
[0173] 若所述收敛判断结果为否，则更新迭代次数，返回至所述分别计算生物腔体的位置r处的第k次的光吸收能量Ak(r)和第k+1次的光吸收能量Ak+1(r)。
[0174] 所述构建光吸收能量分布图像和声速分布图像具体包括：
[0175] 采用如下公式对极坐标下的光吸收能量分布矩阵A进行归一和灰度处理：
[0176]
[0177] B(i,j)是对A(i,j)进行归一化和灰度化处理后的值；
[0178] i是生物腔体位置r处在极坐标系下的横坐标，j是生物腔体位置r处在极标系下的纵坐标；
[0179] A(i,j)是所述光吸收能量分布矩阵A中第i行第j列的元素值；
[0180] min(A)是所述光吸收能量分布矩阵A中元素的最小值，max(A)是所述光吸收能量分布矩阵A中元素的最大值；
[0181] 采用如下公式对极坐标下的声速分布矩阵cs进行归一和灰度处理；
[0182]
[0183] E(i,j)是对cs(i,j)进行归一化和灰度化处理后的值；
[0184] cs(i,j)是所述声速分布矩阵cs中第i行第j列的元素值；
[0185] min(cs)是所述声速分布矩阵cs中的最小值；max(cs)是所述声速分布矩阵cs中的最大值；
[0186] 根据公式B′(x,y)＝B(jcosi,jsini)将所述B(i,j)转化为平面直角坐标中的B′(x,y)；
[0187] 根据公式E′(x,y)＝E(jcosi,jsini)将所述E(i,j)转化为平面直角坐标中的E′(x,y)；
[0188] 其中，i∈[0,360]，j∈[0,d]；
[0189] d是θ‑l极坐标系中极径的最大值；
[0190] B′(x,y)是所述光吸收能量分布图中直角坐标系下点(x,y)的灰度值；E′(x,y)是所述声速分布图中直角坐标系下点(x,y)的灰度值。
[0191] 本发明还公开了一种生物腔体光声内窥成像系统，所述系统包括：
[0192] 初始值获取模块，用于获取生物腔体的位置r处的光吸收能量初始值A0(r)和声速初始值cs,0(r)。
[0193] 本实施例中，光吸收能量初始值取A0(r)＝0，声速初始值取cs,0(r)＝1600m/s。
[0194] 光吸收能量计算模块，用于分别计算生物腔体的位置r处的第k次的光吸收能量Ak(r)和第k+1次的光吸收能量Ak+1(r)；k为迭代次数；k≥1。
[0195] 本实施例中，根据公式Ak(r)＝Ak‑1(r)‑W‑1(Ak‑1(r))G1'(r,cs,k‑1(r),Ak‑1(r))计算所述第k次的光吸收能量Ak(r)和所述第k+1次的光吸收能量Ak+1(r)。
[0196] 其中，
[0197]
[0198] argmin[·]是使g(Ak(r))+λΦTV(Ak(r))取最小值时的Ak(r)；
[0199] W(Ak‑1(r))是所述第一光声关系函数的近似黑塞矩阵；
[0200] W‑1(Ak‑1(r))是W(Ak‑1(r))的逆矩阵；
[0201] A是光吸收能量分布矩阵；A≥0指的是光吸收能量分布矩阵中大于或等于0的光吸收能量；
[0202] g(Ak(r))＝||pm(r)‑H(cs,k‑1(r))·Ak(r)||2是第二光声关系函数；
[0203] pm(r)是位置r处光声信号的测量值；
[0204] H(cs,k‑1(r))是与声速有关的第一算子；
[0205] ||·||是2‑范数；
[0206] ΦTV(Ak(r))是TV正则化项，
[0207] λ是TV正则化参数；
[0208] Lk是关于Ak(r)方差特性的稀疏矩阵；
[0209] η＞0是常数；
[0210] G1'(r,Ak‑1(r),cs,k‑1(r))是所述一光声关系函数的梯度。
[0211] 在具体处理过程中，采用蒙特卡罗模拟的方法，仿真得到在短脉冲激光照射下腔体壁组织中的光能量沉积。然后，采用时域有限差分算法求解极坐标系下的离散光声波动方程，得到组织产生的光声信号的理论值：
[0212] 即
[0213] 其中，(i,j)是腔体横截面上的点r在θ‑l极坐标系中的坐标；Δθ和Δl分别是θ轴(n+1)和l轴上的单位长度；Δt是离散时间间距；n是离散时刻；p (i,j)是位置为(i,j)的质点在n+1时刻产生的光声信号的理论值； 和 分别是位置为(i,j)的质点在时刻
n沿θ方向和l方向的振动速度；cs(i,j)是超声波在位置(i,j)；处的传播速度；β是组织的等(n)
压膨胀系数；Cp是组织的比热容；ρ0是组织的密度；I 是激光脉冲函数在时刻n的值；A(i,j)是位置(i,j)处的光吸收能量。
[0214] 声速计算模块，用于分别计算生物腔体的位置r处的第k次的声速cs,k(r)和第k+1次的声速cs,k+1(r)。
[0215] 本实施例中，
[0216]
[0217] 对第三光声关系函数进行k迭代后，得到位置r处声速为
[0218] 经过求解，
[0219]
[0220] 根据公式 分别计算所述第k次的声速cs,k(r)和所述第k+1次的声速cs,k+1(r)。
[0221] 其中，
[0222]
[0223] ▽f(qk)是f(qk)的梯度；
[0224] f(cs,k‑1(r))＝||pm(r)‑H(cs,k‑1(r))·Ak‑1(r)||2；
[0225]
[0226] γk‑1是第k‑1次迭代后的权重参数，γk‑2是第k‑2次迭代后的权重参数；
[0227] 第k次迭代后的权重参数γk为
[0228] H(qk)是与声速有关的第二算子；
[0229] D是f(cs,k‑1(r))导数的普希茨常数；
[0230] f(cs,k‑1(r))＝||pm(r)‑H(cs,k‑1(r))·Ak‑1(r)||2。
[0231] cs声速分布矩阵，cs≥0是声速分布矩阵中大于或等于0的声速。
[0232] 绝对差计算模块，用于计算所述第k次的光吸收能量Ak(r)与所述第k+1次的光吸收能量Ak+1(r)的光吸收能量绝对差εA,k，以及所述第k次的声速cs,k(r)与所述第k+1次的声速cs,k+1(r)的声速绝对差 所述收敛条件为εA,k＜εA且
[0233] 其中εA为光吸收能量的收敛公差， 为声速的收敛公差。
[0234] 具体的，光吸收能量的收敛公差εA取0.01，声速的收敛公差 取0.01。
[0235] 在迭代计算的过程中，如果所有迭代值均不满足收敛条件，则可以通过改变光吸收能量的收敛公差εA与声速的收敛公差 的取值，调整收敛条件，以使迭代过程完整。
[0236] 收敛判断结果模块，用于判断所述光吸收能量绝对差εA,k和所述声速绝对差 是否满足收敛条件，得到收敛判断结果；
[0237] 图像构建模块，与所述收敛判断结果模块连接，用于在所述收敛判断结果模块为是时，则根据位置r处的第k+1次光吸收能量Ak+1(r)和声速cs,k+1(r)分别构建光吸收能量分布图和声速分布图；
[0238] 迭代更新模块，分别与所述收敛判断结果模块和所述光吸收能量计算模块连接，用于在所述收敛判断结果模块为否时，则更新迭代次数，返回至光吸收能量计算模块。
[0239] 所述图像构建模块具体包括：
[0240] 光吸收能量分布矩阵处理单元，用于采用如下公式对极坐标下的光吸收能量分布矩阵A进行归一和灰度处理：
[0241]
[0242] B(i,j)是对A(i,j)进行归一化和灰度化处理后的值；
[0243] i是生物腔体位置r处在极坐标系下的横坐标，j是生物腔体位置r处在极标系下的纵坐标；
[0244] A(i,j)是所述光吸收能量分布矩阵A中第i行第j列的元素值；
[0245] min(A)是所述光吸收能量分布矩阵A中元素的最小值，max(A)是所述光吸收能量分布矩阵A中元素的最大值；
[0246] 声速分布矩阵处理单元，用于采用如下公式对极坐标下的声速分布矩阵cs进行归一和灰度处理；
[0247]
[0248] E(i,j)是对cs(i,j)进行归一化和灰度化处理后的值；
[0249] cs(i,j)是所述声速分布矩阵cs中第i行第j列的元素值；
[0250] min(cs)是所述声速分布矩阵cs中的最小值；max(cs)是所述声速分布矩阵cs中的最大值；
[0251] 光吸收能量分布图灰度值转换单元，用于根据公式B′(x,y)＝B(jcosi,jsini)将所述B(i,j)转化为平面直角坐标中的B′(x,y)；
[0252] 声速分布图像灰度值转换单元，用于根据公式E′(x,y)＝E(jcosi,jsini)将所述E(i,j)转化为平面直角坐标中的E′(x,y)；
[0253] 其中，i∈[0,360]，j∈[0,d]；
[0254] d是θ‑l极坐标系中极径的最大值；
[0255] B′(x,y)是所述光吸收能量分布图中直角坐标系下点(x,y)的灰度值；E′(x,y)是所述声速分布图中直角坐标系下点(x,y)的灰度值。
[0256] 对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。
[0257] 本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。

专利地区：河北

专利申请日期：2019-08-29

专利公开日期：2022-12-27

专利公告号：CN112438702B