基于广义逆算法的冲击波超压场重建方法*

郭亚丽，韩 焱，王黎明

(中北大学电子测试技术国家重点实验室，山西 太原 030051)

爆炸技术越来越多地应用于国防及民用的各个领域[1]，冲击波超压是爆炸产生杀伤和破坏作用的主要因素，一直是空气中爆炸规律研究的重要课题[2]。目前，虽然已总结了一些计算冲击波超压的经验公式[3]，但在实验研究时，由于测试环境、测试条件以及装药形状不完全相同，直接利用经验公式进行计算是不准确的。针对目前局部布点测试不能全面了解冲击波传播过程的不足以及超压经验公式的局限性，本文中，拟采用网络化测试技术获取冲击波信号，以计算机层析成像技术为基础，利用加权广义逆反演算法对爆炸冲击波速度场进行反演，并根据冲击波速度与峰值超压的关系，得到峰值超压场的分布。

1 冲击波超压场重建原理

炸药在空气中爆炸，形成空气冲击波。将测试区域划分为规则的网格单元，假设自由空气中的理想冲击波不沿网格单元的边界传播，在每个网格单元内近似沿某条直线传播，每个探测器 对应一条射线。冲击波到达探测器 阵元的时间为走时，冲击波在传输过程的走时是速度和几何路径的函数。根据走时层析成像原理，有[4]

DS=T

(1)

式中：T=(t1,t2,…,tm)T为各条射线走时的m维列向量,S=(s1,s2,…,sn)T为待求离散单元慢度值，为n维未知的列向量；D为m×n阶稀疏矩阵，其元素为dij，即第i条射线穿过第j个网格的射线长度。

2 观测系统优化布局

2.1 测试区域网格划分

将测试区域划分为若干规则的网格单元，每个网格单元内波速视为均匀。网格划分越密，层析成像的分辨率越高，但方程的未知数越多，解的不确定性也越高。网格划分应结合测试区域大小、先验信息、实验所用探测器 数目合理划分。

2.2 探测器优化布局及判断指标

探测器布设时，应使射线覆盖面广、分布均匀；减少射线路径矩阵中零元素的个数、降低其条件数。采用射线密度、射线正交性及矩阵D的条件数作为探测器优化布局的判断指标。

射线密度代表通过各网格像元的射线数目，而射线正交性是通过各网格像元的射线之间夹角的最大正弦值来度量，射线密度小和正交性差的区域，图像误差大；反之，则结果比较可靠。在射线总数一定时，根据探测目标形状，模型分布特点等合理分布射线[5]。

3 基于加权广义逆的冲击波场反演算法

从矩阵的观点求方程组(1), 也就是求系数矩阵D的逆矩阵的问题。但是对于爆炸层析成像问题，系数矩阵D一般是奇异矩阵, 不存在通常意义下的逆矩阵, 因此需要采用广义逆理论进行求解。对于冲击波层析成像，除采用稳定性好的广义逆反演算法，还应对反演矩阵进行加权处理形成加权广义逆反演，来增加反演过程中的信息量。

在广义逆反演时，既可对观测数据加权，也可以对模型参数加权。设A为m×n阶矩阵，P、Q分别为m×m、n×n的正定矩阵[6]，若有n×m阶矩阵X满足

AXA=A，XAX=X， (PAX)T=PAX， (QXA)T=QXA

(2)

则称X为A的加权广义逆A+，A+可以表示为A+=Q-1(PAQ-1)+P。

采用自然权矩阵，设P、Q均为对角阵，对称正定，矩阵P的对角线元素为当前模型的各射线的走时；矩阵Q的对角线元素为各射线在第j个网格上的总贡献，表示为射线穿越任何网格单元的长度和与相应单元速度乘积，J.G.Berryman[7]称其为射线分布矩阵。

这种权矩阵的最大特点是物理意义明确,表征了以下基于物理原理的判断：(1)冲击波作为信息载体,随传播距离的增加,衰减幅度增大，波形信号减弱，故射线路径较短者,其信噪比应较高,也较接近真实路径,故对小走时测值加重权。(2)一般射线越密集,重建图像就越精确。基于这样的考虑,对众多射线穿过的网格单元加重权。于是，加权广义逆走时层析成像算法写为S=Q-1(PDQ-1)+PT。

图1 爆炸场速度分布模型Fig.1 The model for velocity distribution in explosion field

图2 网格及探测器分布示意图 Fig.2 Schematic distribution of grids and detectors

4 数值模拟及实验结果

4.1 数值模拟

设置10 m×10 m自由场空气爆炸冲击波速度分布模型,见图1。根据实验条件，将测试区域划分为10×10的网格单元，炸点位于测试区域中心，见图2。根据对称性原理，对1/4区域进行速度场反演。选择的探测器数为13，选取2种不同的布设方式,见图3。图3(a)中探测器节

图3 探测器的2种不同的布设方式Fig.3 Two different distributions of detectors

图4 不同的探测器布设方式下的射线密度分布Fig.4 Half-line density distributions with two different distributions of detectors

图5 不同的探测器布设方式下的射线正交性分布Fig.5 Half-line orthogonality distribution with two different distributions of detectors

图6 速度分布反演的初始模型Fig.6 The initial model for velocity inversion

图7 各个网格单元速度相对误差Fig.7 The relative error of velocity in each grid

点基本呈均匀分布，其坐标分别为:S1(0,4.9)、S2(0,4.15)、S3(0,3.32)、S4(0,2.5)、S5(0,1.66)、S6(0,0.83)、S7(0,0)、S8(0.83,0)、S9(1.66,0)、S10(2.5,0)、S11(3.32,0)、S12(4.15,0)、S13(4.5,0)。图3(b)为按照探测器布设指标优化布局，各节点坐标为:S1(0,4.9)、S2(0,3.8)、S3(0,3)、S4(0,2)、S5(0,1.3)、S6(0,1)、S7(0,0.5)、S8(1.2,0)、S9(2,0)、S10(2.5,0)、S11(3,0)、S12(3.5,0)、S13(4.9,0)。

2种探测器布局方式下,测试区域射线密度及正交性分布见图4～5。由图4～5可以看出:射线总数一定时，探测器优化布局方式下射线分布更合理，正交性更好；同时，矩阵D的条件数为35.61，远小于均匀布局方式下的2.74×1016。可见，探测器优化布局使得方程更稳定。采用相同初始模型(图6)和加权广义逆算法，对2种布局下的爆炸场速度进行反演重建，相对于真实模型(图1)，探测器均匀布局和优化布局的网格平均相对误差分别为8.47%和3.32%，各个网格ng单元速度相对误差εv见图7，反演结果的二维图像见图8。由反演误差及图像分布均可得到：采用探测器优化布局得到的反演结果误差更小，更接近真实模型。

图8 不同的探测器布局方式下的反演结果Fig.8 Inversion results with two different distributions of detectors

图9 实验中探测器布局Fig.9 Detector distribution in experiment

图10 冲击波峰值超压场的重建结果 Fig.10 The reconstruction results for peak overpressure of shock wave

4.2 实验及重建结果

实验时,使用50 kg散装TNT炸药，炸药密度为0.80 g/cm3，形状接近柱形，吊离地面1.7m，进行空中爆炸。对32 m×32 m范围内自由冲击波场进行二维重建，爆炸点居于测试区域中心，使用了30个超压探测器 ，且探测器 与炸点在同一平面，可近似认为是在无限空气中的爆炸。根据有效的射线总数和探测区域的大小以及射线密度将测试区域划分为长、宽各64×64个大小相同的网格，探测器 阵元分布在网格边界和网格内个别节点上，如图9所示。

采用本文算法对爆炸冲击波峰值超压场进行重建，同时，采用如下经验公式计算峰值超压：

(6)

将某次实验测得的冲击波峰值超压pm,e与重建结果pm,r及经验公式计算结果进行比较，见表1。

通过多次实验数据与重建结果的比较，得到在近场区4 m以内，重建结果与实测值的偏差较大，接近10%，这是由于实验时，近场区干扰较大，实验数据本身也存在一定误差，而且测试点数较少，导致重建结果与实测值的偏差增大，4 m以外的测试区域内，采用本文算法重建得到的峰值超压与实测值较接近，相对偏差在5%以内，优于经验公式计算结果。

实验时，若增加超压探测器 数目，重建精度会大大提高。峰值超压场二维重建结果如图10所示，从图中可以看出，测试区域内冲击波峰值超压随距离的分布情况，距离爆炸中心5 m以内的区域超压衰减及其分布形状较明显，且其分布形式与装药形状基本一致；而5 m以外的区域由于冲击波超压衰减幅度较大，故分布不明显。

5 结 论

针对冲击波超压经验公式的局限性，采用网络化测试技术获取冲击波信号，利用不完全数据重建技术得到一定区域内爆炸冲击波速度场及峰值超压场分布。在一定范围内，重建结果与实测值较接近。然而在爆炸近场区测试时，对实验仪器性能要求较高，所测数据本身存在一定误差，且所测点数较少，故重建结果与实测值的偏差较大。本文中给出了一种冲击波速度和峰值超压的计算方法，但还需要更多的实验数据进行验证。

[1] 张守中.爆炸基本原理[M].北京:国防工业出版社,1988,397-530.

[2] W E贝克.空中爆炸[M〗.江科,译.北京:原子能出版社,1982:40-85.

[3] 李翼祺,马素贞.爆炸力学[M].北京:科学出版社,1992,259-262.

[4] Kepler W F, Bond L J, Frangopol D M.Improved assessment of mass concrete dams using acoustic travel time tomography: II: Aapplication[J].Construction and Building Materials, 2000,14(3):147-156.

[5] 裴正林,余钦范,狄帮让.井间地震层析成像分辨率研究[J].物探与化探,2006,26(3):218-224.Pei Zheng-lin, Yu Qin-fan, Di Bang-rang.Resolution in crosshole seismic tomography[J].Geophysical and Geochemical Exploration, 2006,26(3):218-224.

[6] 王振宇,刘国华,梁国钱.基于广义逆的层析成像反演方法研究[J].浙江大学学报:工学版,2005,39(1):1-5.Wang Zhen-yu, Liu Guo-hua, Liang Guo-qian.Study on inversion methods for travel time tomography based on generalized inverse theory[J].Journal of Zhejiang University: Engineering Science, 2005,39(1):1-5.

[7] Berryman J G.Fermat’s principle and nonlinear travel time tomography[J].Physical Review Letter, 1989,62(25):2953-2956.