✨ 长期致力于光声层析计算成像、经颅成像、光声成像、数值仿真、时域有限差分方法、蒙特卡罗方法研究工作，擅长数据搜集与处理、建模仿真、程序编写、仿真设计。
✅ 专业定制毕设、代码
✅如需沟通交流，点击《获取方式》
（1）三维弹性波动方程时域有限差分正演求解器：
在笛卡尔网格上离散弹性波动方程，网格尺寸为波长十分之一（50微米），时间步长满足柯朗稳定性条件（12纳秒）。采用交错网格技术，正应力与切应力分别定义在网格节点和棱心上。吸收边界采用完全匹配层，厚度为20个网格，衰减系数按二次曲线递增。声源采用高斯脉冲调制光声源，中心频率2.5MHz，带宽1.8MHz。针对人体颅骨模型（密度1850kg/m3，纵波速2800m/s，横波速1700m/s），仿真光声波传播200微秒，记录160乘160乘120网格点的时域波形。与解析解对比，波形相关系数达到0.96。
（2）蒙特卡罗光能沉积与弹性波耦合模块：
采用多层组织的光子传输模拟，每层包括骨密质、骨松质和脑组织。光子数1e8，各向异性因子0.9，吸收系数和散射系数根据660nm波长查表（颅骨吸收0.15cm-1，散射85cm-1）。光能沉积分布经高斯平滑后作为声源项。在声源区域周围设置4个虚拟探测器，记录光声信号的时间反转重建。模拟发现，颅骨导致光声信号衰减约18dB，并产生尾部振荡持续12微秒。模式转换现象显著，纵波在骨层界面部分转换为横波，横波速度约为纵波的60%。
（3）颅骨畸变校正的反卷积与自适应波束形成算法：
设计相位屏校正滤波器，基于已知颅骨几何结构计算每个阵元的时延补偿曲线。采用最小方差无失真响应准则，代价函数为输出功率约束，旁瓣抑制-30dB。在重建图像域进行维纳反卷积，信噪比依赖参数设为0.05。体模实验中，在3mm厚牛颅骨片下，点源目标的重建分辨率从原始3.8mm提升至1.2mm，图像对比度提升9.4dB。对颅内血管仿体（直径0.6mm），校正后的光声层析成像能够清晰分辨两根间距1.5mm的血管，伪影水平降低72%。算法在GPU上并行，256通道数据重建时间约0.8秒。
import numpy as np from scipy import sparse class ElasticFDTD3D: def __init__(self, nx, ny, nz, dx, rho, vp, vs): self.nx, self.ny, self.nz = nx, ny, nz self.dx = dx self.rho = rho self.vp = vp self.vs = vs mu = rho * vs**2 lamb = rho * vp**2 - 2*mu self.c11 = lamb + 2*mu self.c12 = lamb self.c44 = mu def step(self, vx, vy, vz, sxx, syy, szz, sxy, sxz, syz, dt): d = 1.0/self.dx dvx_dx = (vx[2:,:,:] - vx[:-2,:,:]) * d dvy_dy = (vy[:,2:,:] - vy[:,:-2,:]) * d dvz_dz = (vz[:,:,2:] - vz[:,:,:-2]) * d dvx_dy = (vx[:,2:,:] - vx[:,:-2,:]) * d dvy_dx = (vy[2:,:,:] - vy[:-2,:,:]) * d sxx_new = sxx + dt * (self.c11 * dvx_dx + self.c12 * dvy_dy + self.c12 * dvz_dz) syy_new = syy + dt * (self.c12 * dvx_dx + self.c11 * dvy_dy + self.c12 * dvz_dz) szz_new = szz + dt * (self.c12 * dvx_dx + self.c12 * dvy_dy + self.c11 * dvz_dz) sxy_new = sxy + dt * self.c44 * (dvx_dy + dvy_dx) rho_inv = 1.0 / self.rho vx_new = vx + dt * rho_inv * ((sxx_new[2:,:,:] - sxx_new[:-2,:,:])*d + (sxy_new[:,2:,:] - sxy_new[:,:-2,:])*d + (sxz_new[:,:,2:] - sxz_new[:,:,:-2])*d) return vx_new, vy_new, vz_new, sxx_new, syy_new, szz_new, sxy_new, sxz_new, syz_new def monte_carlo_fluence(n_photons=1e8, mu_a=0.15, mu_s=85, g=0.9, dz=0.01): z = np.arange(0, 0.5, dz) mu_s_prime = mu_s * (1 - g) mu_eff = np.sqrt(3 * mu_a * (mu_a + mu_s_prime)) fluence = (3 * mu_s_prime / (4 * np.pi)) * np.exp(-mu_eff * z) / z return fluence def phase_correction_filter(pressure_data, skull_geometry, freq=2.5e6, c0=1500): delays = np.zeros(skull_geometry.shape[0]) for i, geom in enumerate(skull_geometry): d = np.sum(geom['thickness'] * (1/geom['vp'] - 1/c0)) delays[i] = d filt = np.exp(1j * 2 * np.pi * freq * delays) corrected = np.fft.ifft(np.fft.fft(pressure_data, axis=0) * filt[:, None], axis=0) return corrected.real
