具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明的目的是提供一种生物腔体光声内窥成像方法及系统，考虑到生物组织的复杂性，在构建光吸收能量分布图像过程中结合声速不同的条件，使得构建分图像更加准确。

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。

一种生物腔体光声内窥成像方法，所述方法包括：

获取生物腔体的位置r处的光吸收能量初始值A₀(r)和声速初始值c_s,0(r)。

本实施例中，光吸收能量初始值取A₀(r)＝0，声速初始值取c_s,0(r)＝1600m/s。

分别计算生物腔体的位置r处的第k次的光吸收能量A_k(r)和第k+1次的光吸收能量A_k+1(r)；k为迭代次数，k≥1；

具体的，本实施例根据公式A_k(r)＝A_k-1(r)-W^-1(A_k-1(r))G₁'(r,c_s,k-1(r),A_k-1(r))计算所述第k次的光吸收能量A_k(r)和所述第k+1次的光吸收能量A_k+1(r)；

其中，

argmin[·]是使g(A_k(r))+λΦ_TV(A_k(r))取最小值时的A_k(r)；

W(A_k-1(r))是所述第一光声关系函数的近似黑塞矩阵；

W^-1(A_k-1(r))是W(A_k-1(r))的逆矩阵；

A是光吸收能量分布矩阵；A≥0指的是光吸收能量分布矩阵中大于或等于0的光吸收能量；

g(A_k(r))＝||p_m(r)-H(c_s,k-1(r))·A_k(r)||²是第二光声关系函数；

p_m(r)是位置r处光声信号的测量值；

H(c_s,k-1(r))是与声速有关的第一算子；

p(r)是位置r处光声信号的理论值；

||·||是2-范数；

Φ_TV(A_k(r))是TV正则化项，

λ是TV正则化参数；

L_k是关于A_k(r)方差特性的稀疏矩阵；

η＞0是常数；

G₁'(r,A_k-1(r),c_s,k-1(r))是所述一光声关系函数的梯度。

在具体处理过程中，采用蒙特卡罗模拟的方法，仿真得到在短脉冲激光照射下腔体壁组织中的光能量沉积。然后，采用时域有限差分算法求解极坐标系下的离散光声波动方程，得到组织产生的光声信号的理论值：

其中，(i,j)是腔体横截面上的点r在θ-l极坐标系中的坐标；Δθ和Δl分别是θ轴和l轴上的单位长度；Δt是离散时间间距；n是离散时刻；p⁽ⁿ⁺¹⁾(i,j)是位置为(i,j)的质点在n+1时刻产生的光声信号的理论值；和分别是位置为(i,j)的质点在时刻n沿θ方向和l方向的振动速度；c_s(i,j)是超声波在位置(i,j)；处的传播速度；β是组织的等压膨胀系数；C_p是组织的比热容；ρ₀是组织的密度；I⁽ⁿ⁾是激光脉冲函数在时刻n的值；A(i,j)是位置(i,j)处的光吸收能量。

分别计算生物腔体的位置r处的第k次的声速c_s,k(r)和第k+1次的声速c_s,k+1(r)；

本实施例中，

对第三光声关系函数进行k迭代后，得到位置r处声速为

经过求解，

根据公式分别计算所述第k次的声速c_s,k(r)和所述第k+1次的声速c_s,k+1(r)。

其中，

▽f(q_k)是f(q_k)的梯度；

f(q_k)＝||p_m(r)-H(q_k)·A_k-1(r)||²；

γ_k-1是第k-1次迭代后的权重参数，γ_k-2是第k-2次迭代后的权重参数；

第k次迭代后的权重参数γ_k为

H(q_k)是与声速有关的第二算子；

D是f(c_s,k-1(r))导数的普希茨常数；

f(c_s,k-1(r))＝||p_m(r)-H(c_s,k-1(r))·A_k-1(r)||²；

c_s声速分布矩阵，c_s≥0是声速分布矩阵中大于或等于0的声速。

计算所述第k次的光吸收能量A_k(r)与所述第k+1次的光吸收能量A_k+1(r)的光吸收能量绝对差ε_A,k，以及所述第k次的声速c_s,k(r)与所述第k+1次的声速c_s,k+1(r)的声速绝对差所述收敛条件为ε_A,k＜ε_A且其中ε_A为光吸收能量的收敛公差，为声速的收敛公差。

本实施例中，光吸收能量的收敛公差ε_A取0.01，声速的收敛公差取0.01。

在迭代计算的过程中，如果所有迭代值均不满足收敛条件，则可以通过改变光吸收能量的收敛公差ε_A与声速的收敛公差的取值，调整收敛条件，以使迭代过程完整。

判断所述光吸收能量绝对差ε_A,k和所述声速绝对差是否满足收敛条件，得到收敛判断结果；

若所述收敛判断结果为是，则输出位置r处的第k+1次光吸收能量A_k+1(r)和声速c_s,k+1(r)，分别构建光吸收能量分布图和声速分布图；

若所述收敛判断结果为否，则更新迭代次数，返回至所述分别计算生物腔体的位置r处的第k次的光吸收能量A_k(r)和第k+1次的光吸收能量A_k+1(r)。

所述构建光吸收能量分布图像和声速分布图像具体包括：

采用如下公式对极坐标下的光吸收能量分布矩阵A进行归一和灰度处理：

B(i,j)是对A(i,j)进行归一化和灰度化处理后的值；

i是生物腔体位置r处在极坐标系下的横坐标，j是生物腔体位置r处在极标系下的纵坐标；

A(i,j)是所述光吸收能量分布矩阵A中第i行第j列的元素值；

min(A)是所述光吸收能量分布矩阵A中元素的最小值，max(A)是所述光吸收能量分布矩阵A中元素的最大值；

采用如下公式对极坐标下的声速分布矩阵c_s进行归一和灰度处理；

E(i,j)是对c_s(i,j)进行归一化和灰度化处理后的值；

c_s(i,j)是所述声速分布矩阵c_s中第i行第j列的元素值；

min(c_s)是所述声速分布矩阵c_s中的最小值；max(c_s)是所述声速分布矩阵c_s中的最大值；

根据公式B′(x,y)＝B(jcosi,jsini)将所述B(i,j)转化为平面直角坐标中的B′(x,y)；

根据公式E′(x,y)＝E(jcosi,jsini)将所述E(i,j)转化为平面直角坐标中的E′(x,y)；

其中，i∈[0,360]，j∈[0,d]；

d是θ-l极坐标系中极径的最大值；

B′(x,y)是所述光吸收能量分布图中直角坐标系下点(x,y)的灰度值；E′(x,y)是所述声速分布图中直角坐标系下点(x,y)的灰度值。

本发明还公开了一种生物腔体光声内窥成像系统，所述系统包括：

初始值获取模块，用于获取生物腔体的位置r处的光吸收能量初始值A₀(r)和声速初始值c_s,0(r)。

本实施例中，光吸收能量初始值取A₀(r)＝0，声速初始值取c_s,0(r)＝1600m/s。

光吸收能量计算模块，用于分别计算生物腔体的位置r处的第k次的光吸收能量A_k(r)和第k+1次的光吸收能量A_k+1(r)；k为迭代次数；k≥1。

本实施例中，根据公式A_k(r)＝A_k-1(r)-W^-1(A_k-1(r))G₁'(r,c_s,k-1(r),A_k-1(r))计算所述第k次的光吸收能量A_k(r)和所述第k+1次的光吸收能量A_k+1(r)。

其中，

argmin[·]是使g(A_k(r))+λΦ_TV(A_k(r))取最小值时的A_k(r)；

W(A_k-1(r))是所述第一光声关系函数的近似黑塞矩阵；

W^-1(A_k-1(r))是W(A_k-1(r))的逆矩阵；

A是光吸收能量分布矩阵；A≥0指的是光吸收能量分布矩阵中大于或等于0的光吸收能量；

g(A_k(r))＝||p_m(r)-H(c_s,k-1(r))·A_k(r)||²是第二光声关系函数；

p_m(r)是位置r处光声信号的测量值；

H(c_s,k-1(r))是与声速有关的第一算子；

||·||是2-范数；

Φ_TV(A_k(r))是TV正则化项，

λ是TV正则化参数；

L_k是关于A_k(r)方差特性的稀疏矩阵；

η＞0是常数；

G₁'(r,A_k-1(r),c_s,k-1(r))是所述一光声关系函数的梯度。

在具体处理过程中，采用蒙特卡罗模拟的方法，仿真得到在短脉冲激光照射下腔体壁组织中的光能量沉积。然后，采用时域有限差分算法求解极坐标系下的离散光声波动方程，得到组织产生的光声信号的理论值：

即

其中，(i,j)是腔体横截面上的点r在θ-l极坐标系中的坐标；Δθ和Δl分别是θ轴和l轴上的单位长度；Δt是离散时间间距；n是离散时刻；p⁽ⁿ⁺¹⁾(i,j)是位置为(i,j)的质点在n+1时刻产生的光声信号的理论值；和分别是位置为(i,j)的质点在时刻n沿θ方向和l方向的振动速度；c_s(i,j)是超声波在位置(i,j)；处的传播速度；β是组织的等压膨胀系数；C_p是组织的比热容；ρ₀是组织的密度；I⁽ⁿ⁾是激光脉冲函数在时刻n的值；A(i,j)是位置(i,j)处的光吸收能量。

声速计算模块，用于分别计算生物腔体的位置r处的第k次的声速c_s,k(r)和第k+1次的声速c_s,k+1(r)。

本实施例中，

对第三光声关系函数进行k迭代后，得到位置r处声速为

经过求解，

根据公式分别计算所述第k次的声速c_s,k(r)和所述第k+1次的声速c_s,k+1(r)。

其中，

▽f(q_k)是f(q_k)的梯度；

f(c_s,k-1(r))＝||p_m(r)-H(c_s,k-1(r))·A_k-1(r)||²；

γ_k-1是第k-1次迭代后的权重参数，γ_k-2是第k-2次迭代后的权重参数；

第k次迭代后的权重参数γ_k为

H(q_k)是与声速有关的第二算子；

D是f(c_s,k-1(r))导数的普希茨常数；

f(c_s,k-1(r))＝||p_m(r)-H(c_s,k-1(r))·A_k-1(r)||²。

c_s声速分布矩阵，c_s≥0是声速分布矩阵中大于或等于0的声速。

绝对差计算模块，用于计算所述第k次的光吸收能量A_k(r)与所述第k+1次的光吸收能量A_k+1(r)的光吸收能量绝对差ε_A,k，以及所述第k次的声速c_s,k(r)与所述第k+1次的声速c_s,k+1(r)的声速绝对差所述收敛条件为ε_A,k＜ε_A且

其中ε_A为光吸收能量的收敛公差，为声速的收敛公差。

具体的，光吸收能量的收敛公差ε_A取0.01，声速的收敛公差取0.01。

在迭代计算的过程中，如果所有迭代值均不满足收敛条件，则可以通过改变光吸收能量的收敛公差ε_A与声速的收敛公差的取值，调整收敛条件，以使迭代过程完整。

收敛判断结果模块，用于判断所述光吸收能量绝对差ε_A,k和所述声速绝对差是否满足收敛条件，得到收敛判断结果；

图像构建模块，与所述收敛判断结果模块连接，用于在所述收敛判断结果模块为是时，则根据位置r处的第k+1次光吸收能量A_k+1(r)和声速c_s,k+1(r)分别构建光吸收能量分布图和声速分布图；

迭代更新模块，分别与所述收敛判断结果模块和所述光吸收能量计算模块连接，用于在所述收敛判断结果模块为否时，则更新迭代次数，返回至光吸收能量计算模块。

所述图像构建模块具体包括：

光吸收能量分布矩阵处理单元，用于采用如下公式对极坐标下的光吸收能量分布矩阵A进行归一和灰度处理：

B(i,j)是对A(i,j)进行归一化和灰度化处理后的值；

i是生物腔体位置r处在极坐标系下的横坐标，j是生物腔体位置r处在极标系下的纵坐标；

A(i,j)是所述光吸收能量分布矩阵A中第i行第j列的元素值；

min(A)是所述光吸收能量分布矩阵A中元素的最小值，max(A)是所述光吸收能量分布矩阵A中元素的最大值；

声速分布矩阵处理单元，用于采用如下公式对极坐标下的声速分布矩阵c_s进行归一和灰度处理；

E(i,j)是对c_s(i,j)进行归一化和灰度化处理后的值；

c_s(i,j)是所述声速分布矩阵c_s中第i行第j列的元素值；

min(c_s)是所述声速分布矩阵c_s中的最小值；max(c_s)是所述声速分布矩阵c_s中的最大值；

光吸收能量分布图灰度值转换单元，用于根据公式B′(x,y)＝B(jcosi,jsini)将所述B(i,j)转化为平面直角坐标中的B′(x,y)；

声速分布图像灰度值转换单元，用于根据公式E′(x,y)＝E(jcosi,jsini)将所述E(i,j)转化为平面直角坐标中的E′(x,y)；

其中，i∈[0,360]，j∈[0,d]；

d是θ-l极坐标系中极径的最大值；

B′(x,y)是所述光吸收能量分布图中直角坐标系下点(x,y)的灰度值；E′(x,y)是所述声速分布图中直角坐标系下点(x,y)的灰度值。

对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。

本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。