爆破振动信号波形预测模型的优化研究.pdf

第2 8 卷第1 期 爆破 V 0 1 ．2 8N o ．1 2 0 1 1 年3 月 B L A S T I N G ， M a r ．2 0 1 1 D O I 1 0 ．3 9 6 3 ／j ．i s s n ．1 0 0 1 4 8 7 X ．2 0 11 ．0 1 ．0 0 6 爆破振动信号波形预测模型的优化研究木 燕永峰1 ’2 ，陈士海1 , 2 ，张安康1 ’2 ，徐琪1 ’2 1 ．山东科技大学土木建筑学院，青岛2 6 6 5 1 0 ； 2 ．山东科技大学山东省土木工程防灾减灾重点实验室，青岛2 6 6 5 1 0 摘要依托临沂引水隧洞竖井爆破工程，通过M a t l a b 语言编程对实测数据进行小波包分析，并利用修正 萨道夫斯基公式 简称修正公式 对各个频带的信号进行了非线性预测分析，最终整合得到爆区内任意点的 爆破振动信号波形的预测模型。分析发现，修正萨道夫斯基公式的预测精度明显高于传统的经验公式，且基 于修正萨道夫斯基公式构建的多频带小泼包系数波形预测模型的预测精度较高 幅值预测误差不足2 ％，主 频预测误差6 ．5 ％ ，较单一衰减因子分析模型稳定。 关键词爆破振动；波形预测；小波包分析；修正萨道夫斯基公式 中图分类号0 3 8 1 文献标识码A文章编号 1 0 0 1 4 8 7 X 2 0 1 1 0 1 0 0 2 3 0 5 O p t i m i z a t i o nR e s e a r c ho nW a v e f o r mP r e d i c t i o n M o d e lo fB l a s t i n gV i b r a t i o nS i g n a l Y A NY o n g - f e n 9 1 ”，C H E NS h i ．h a i l ”，Z H A N GA n ，k a n 9 1 ’，X UQ i l ’2 1 ．C o l l e g eo fA r c h i t e c t u r ea n dC i v i lE n g i n e e r i n g ，S h a n d o n gU n i v e r s i t yo fS c i e n c ea n d T e c h n o l o g y ，Q i n g d a o2 6 6 5 1 0 ，C h i n a ；2 ．S h a n d o n gP r o v i n c i a lK e yL a b o r a t o r yo fC i d l E n g i n e e r i n gD i s a s t e rP r e v e n t i o na n dM i t i g a t i o n ，S D U S T ，Q i n g d a o2 6 6 510 ，C h i n a A b s t r a c t B a s e do nt h ev e r t i c a ls h a f tb l a s t i n gp r o j e c to ft h ed i v e r s i o nt u n n e li nL i n y iC i t y 。t I l ew a v e l e tp a c k e ta ． n a l y s i sp r o g r a m m m e dw i t hM a t l a bW a sc a r r i e do u tt oa n a l y z et h eo r i g i n a lm e a s u r e dd a t aa n dt h es i s r l a l si nv a r i o u s b a n d sw e r ep r e d i c t e dr e s p e c t i v e l yw i t hm e d i f i e dS a d a o v s kF o r m u l as o l v e db yn o n l i n e a rr e g r e s s i o na n a l y s i s ，a n dt h e n t h eb l a s t i n gv i b r a t i o nw a v e f o r mp r e d i c t i o nm o d e lo fa n yp o i n ti nt h eb l a s ta r e aw a so b t a i n e de v e n t u a l l y ．T h ea n a l y s e s s h o wt h a tt h ep r e d i c t i o np r e c i s i o no fm o d i f i e dS a d a o v s kF o r m u l ai se v i d e n t l yh i g h e rt h a nt h a to ft r a d i t i o n a le m p i r i c M f o r m u l a ．A n dt h ep r e d i c t i o np r e c i s i o no fm u l t i ．b a n dw a v e l e tc o e f f i c i e n tw a v e f o t i np r e d i c t i o nm o d e lb a s e do nt h em o d i ． f l e dS a d a o v s kF o r m u l ai sr e l a t i v e l yh i g h e r T h ep r e d i c t i o ne l r o ro fa m p l i t u d ei sl e s st h a n2 ％a n dt h a t o fd o m i n a n t f r e q u e n c yi so n l y6 ．5 ％ ，w h i c hi sm o r es t a b l ea n dr e l i a b l et h a nt h es i n g l ea t t e n u a t i o nf a c t o ra n a l y s i sm e t h o d ． K e yw o r d s b l a s t i n gv i b r a t i o n ；w a v e f o r mp r e d i c t i o n ；w a v e l e tp a c k e ta n a l y s i s ；m o d i f i e dS o t a o v s kf o r m u l a 一直以来，爆破振动预测是各国爆破振动研究 者致力解决的重大问题，很多学者从理论与经验的 不同角度进行了有益的尝试与探索。综合国内外在 收藕日期2 0 1 0 0 7 一0 2 作者简介燕永峰 1 9 8 6 ，男，硕士研究生，主要从事爆破震动预测 方面的研究， E m a i l s k d y y f 1 6 3 ．C O r n 。 基金项目国家自然科学基金 N o ．5 0 8 7 8 1 2 3 和N o ．5 0 7 7 8 1 0 7 ；山东 省“泰山学者”建设工程基金；山东省土木工程防灾减灾 重点实验室基金 这方面取得的成果，爆破预测的发展可分为爆破振 动强度的预测和爆破振动波形的预测两个阶段。早 期的预测一般是围绕着爆破振动强度进行的，也就 是对爆破振动峰值速度 P P V 的预测，典型的例子 是前苏联的萨道夫斯基经验公式；后期，随着对爆破 振动特征及其破坏机理研究的深入，人们开始着眼 于爆破振动信号波形的预测。爆破振动波形预测的 实现便于施工设计人员在事前准确了解爆破过程中 万方数据 爆破2 0 1 1 年3 月 产生的振动信号各频率成分的幅值、主频、相位分 布、能量大小和所占能量比例，进而为防灾减灾工作 提供了明确的目标和准确的依据，这方面的研究已 取得不少的成果。 对于爆破振动波形的预测，广泛采用S c a n l a n 和S a c h s 提出的利用三角级数叠加来模拟振动加速 度波形，其基本思路是用一组三角级数之和构造一 个近似的平稳高斯过程，然后乘以强度包络线从而 得到非平稳的地面加速度波形⋯，此种方法需要建 立在大量实测数据的基础之上，在爆破振动控制中 由于实测数据样本有限导致此方法应用受限；吴从 师等在1 9 9 0 年研究合理的微差时间时将单孔爆破 振动实际记录波形作为源函数，通过线性叠加实现 多孑L 爆破振动波形的模拟预测旧1 ；徐全军等在1 9 9 8 年采用单孑L 样本函数的线性叠加来模拟多段微差爆 破振动的波形，并通过混凝土实验模型对此进行了 验证旧1 ；另外，刘军HJ 、D o w d i n g 等也对爆破振动的 预测模型做了深入的研究口j 。 以上的各种研究都是建立在单段波形的线性叠 加基础之上，此类方法对已测数据的依赖性很强且 误差较大。由于在爆破信号传播过程中不同频率下 的振速衰减情况存在差异【6 ．7 J ，因此采用分段的方 法分别预测各频带的振速衰减状况较传统的整个振 动信号衰减规律采用单一的k 、a 值要精确。截至目 前，很少有关于此方面研究的报道。 1多频小波包系数爆破振动波形预测 模型 短时傅里叶变换对信号的频带划分是线性等间 隔的。多分辨分析可以对信号进行有效地时频分 解，但由于其尺度是按二进制变化的，所以在高频段 其频率分辨率较差，而在低频段其时间分辨率较差。 在对频带进行多层次划分时，小波包分析对小波分 析过程中没有细分的高频部分做了进一步的分解， 从而提高了爆破振动信号的时一频分辨率。因此，在 对未知结果的爆破振动预测中小波包分析具有更广 泛的应用价值“ J 。 此波形预测模型以有限元方法为基本思想，通 过小波包分析将爆破信号分解成多个频带，然后利 用修正的萨道夫斯基公式分别预测各个频带上的信 号衰减情况，这样弥补了原始的波形预测模型中因 各频带的衰减情况不同而造成误差的不足，基本思 路如下 1 利用小波包分解技术将实测数据样本进行 小波包分解，然后利用牛顿迭代法对修正萨道夫斯 基公式 见式1 进行求解，分别提取各频带上具有 场地特征的振动衰减系数k 、a 、口。 对于爆破振动速度预测，我国爆破安全规程 G B 6 7 2 2 2 0 0 3 中推荐的是萨道夫斯基经验公式， 然而大量的实践表明药量和爆心距的最佳指数比例 并不是1 ／3 的固定模式一。1 1 1 。为了提高预测精度， 本模型采用修正的萨道夫斯基公式2 l ，其表达式如 下 V k Q 8 舻 1 式中k 为经验系数；a ，口为与场地有关的衰减指 数；Q 、R 、V 分别表示爆破中的装药量 k g 、爆心距 m 和质点振动峰值速度 m ／s 。 根据最t J 、- - 乘法公式拟合原理，待定系数I | } ，a ， 口应使实测值与拟合值的误差最小，数值分析时也 就是保证B ∑ K J | } Q 硝 2 最小。由多元函数 的极值定理知，曰取最小值的必要条件是 蓑 蓑 差 0 2 a 后a aa 8 、7 即 ∑ K 一忌Q 砰 Q 砰 0 i 1 ∑ K 一| | } Q 硝 Q 磁I nQ i 0 i 一 ∑ v i k Q T R Q 群I nR j 0 i 1 3 然后利用M a t l a b 编程对场地特征系数J | l 、a 、j i B 进行回归拟合。 2 从实测数据样本中选取一组典型的振动信 号作为基函数，并对其进行小波包分解，得到各频带 的小波包系数后分别对其进行归一化处理。其中频 带小波包系数的计算公式为“ ] a i J 2J - - ms t f i j t d t 4 式中s t 为爆破振动的基函数；Z ， t 表示二进小 波函数Z J t 的共轭函数，i 为频带数，_ 『为采样点 数。 3 利用前两步中得到的参数，建立多频带小波 包系数预测模型，相应的计算表达式为 A u 口i i 反Q q 群‘ 5 式中a i ． ／为基函数多频带小波包系数口。峰值的归 一化；k i ，a ；，反分别表示场地特征参数。 对小波包系数进行归一化时采用t r a m n m x 函 数，即 。。 2 r i i 三皇乏亏亨‰一 6 万方数据 第2 8 卷第l 期燕永峰，陈士海，张安康，等爆破振动信号波形预测模型的优化研究 2 5 使得归一化后的数据口0 分布在[ 一l ，1 ] 区间内。 4 将A “重构得出爆破振动预测信号。 ∞ ∞ j s z ∑∑A i Z J ￡ 7 式中Z ．j t 2 “砜2 1 t - j 为二进小波函数，其中 2 ‘为尺度、2 7 为平移。 对于微差爆破，在确定微差时间的前提下，可将 多段波形归结为单段波形在不同微差时刻共同叠加 作用的结果，进而可用已预测的单段波形叠加整合 成多段爆破振动信号波形，其原理表达式可写成 s T ∑A t ‘s I f 一“ ’艿 瓦 8 式中J s r 为全局坐标系下监测点处总的质点运动 速度，m ／s ；T 为全局坐标时间，s ；A 。为第k 段炸药的 药量振动系数，当爆源至测点周围的地质条件基本 相同，各段药量和炸药品种均相差不大时取l ； 5 。 t t 。 为局部坐标系下第k 段炸药在测点产生的 质点运动速度，m ／s ；t 。为局部坐标系下第k 段炸药 从起爆到震动波传至监测点所耗时间，s ；n 为段数； 瓦为全局坐标下第k 段炸药的起爆时刻，s ；6 瓦 为分段函数，其取值如式 9 所示。 r 0 TZ 6 瓦 { ‘ 9 L l T ≥瓦 同实测信号一样，预测信号仍是速度关于时间 变化的信号，所以可直接计算得到预测信号的P P V 和总能量；另外，通过快速傅立叶变换可将预测信号 转化为功率谱密度关于频率变化的信号，进而可得 出预测信号的主频；同理，利用第3 步的计算结果可 得到各频带的主频、能量及所占的能量比例。因此， 此模型可方便、全面地预测爆破过程中产生的各种 参数，为实际工程提供详细的依据。 2 工程应用 2 ．1 工程概况 临沂引水隧洞竖井爆破工程总长约2k m ，途经 人口密集的村庄，周边建筑均为一般砖砌或砌块结 构。工程采用全断面掘进微差爆破，总装药4 8k g ， 使用l 一9 段雷管起爆，最大段药量为1 5k ，现场采 用I D T S 3 8 5 0 型三维爆破振动测量仪，每个测点布置 3 个传感器，分别对振动信号在水平径向 测点与爆 源之间连线的水平方向 、水平切向 与径向相垂直 的水平方向 和垂直方向 垂直于地面方向 3 个方 向的分量进行监测，测点布置见图1 。 民房 民房 测 测点1 回 回回 囫回回 回回回 睦圆尸 l 回亮 图1 测点布置示意图 F i g ．1L a y o u to ft h em e a s u r i n gp o i n t s 2 ．2 数据分析 为了减少工作量，只对所采集的1 1 个样本数据 的垂直方向振动信号进行分析3 | 。在小波包分析 时，选用d b 6 作为小波基函数，对各测点爆破振动信 号进行尺度为3 的小波包分解和重构，则在分解后的 8 个频带 0 1 2 5H z 、1 2 5 ～2 5 0H z 、2 5 0 ～3 7 5H z 、 3 7 5 ～5 0 0H z 、5 0 0 ～6 2 5H z 、6 2 5 7 5 0H z 、7 5 0 ～ 8 7 5H z 、8 7 5 10 0 0H z 上的信号仍是速度关于时 间变化的信号，最后经快速傅立叶变换转化为功率 谱密度关于频率变化的信号，分解后各频带上的信 号参数见表1 。 表1实测振动信号各频带振动参数汇总表 T a b l e1 S u m m a r ys h e e to fv i b r a t i o np a r a m e t e r si nv a r i o u sb a n d so fm e a s u r e dv i b r a t i o ns i g n a l s 万方数据 爆破 2 0 1 1 年3 月 结合样本信号，分别利用修正萨道夫斯基公式 简称修正公式 与传统萨道夫斯基公式 简称传统 公式 2 种模型对爆心距离为2 6 ．3i n ，最大段药量 为1 5k g 处测点的爆破振动峰值 P P V 进行了预 测，预测结果见表2 。 表2 单一衰减因子下各预测模型的预测结果对比表 T a b l e 2C o n t r a s tt a b l eo ft h ep r e d i c t e dr e s u l t so fd i f f e r e n t p r e d i c t i o nm o d e li nt h es i n g l ea t t e n u a t i o nf a c t o r 模型 c m P P ．V s ／．。 F 查型氅 误差／％ 实测值 1 ．3 9 1 修正公式 1 ．3 7 5 2 5 2 ．7 4 10 ．3 1 21 ．8 5 31 ．1 5 传统公式 1 ．6 0 61 9 3 ．7 1 3O ．6 7 5 2 ．0 2 51 5 ．4 6 可以看出，修正公式预测的P P V 误差仅1 ．5 ％． 其效果远远优越于传统公式的预测效果 误差为 1 5 ．4 6 ％ 。因此，在波形预测中采用修正公式才能 得到更高的预测精度；另外，利用修正公式得到的 n 月的比值并非通常采用的一1 ／3 ，即一1 ／3 并非所 有场地的最佳指数比。 为了克服单一衰减因子情况下各预测模型对振 动信号振动峰值外其它特性的不可预测性，表3 给 出了采用修正萨道夫斯基公式多频小波包系数波形 预测模型 简称修正萨氏多频带小波包系数预测模 型 得到的测点 爆心距离为2 6 ．3m ，最大段药量为 1 5k g 处振动信号各频带的峰值速度和主频。 表3 预测点实测信号与预测信号各频带参数汇总表 T a b l e 3S u m m a r ys h e e to ft h em e a s u r e ds i g n a la n dt h ep r e d i c t e ds i g n a li nv a r i o u sb a n d so fp r e d i c t i n gp o i n t s 结果显示，利用修正萨氏多频带小波包系数预 测模型得到的预测信号的主频与实测信号的主频的 相对误差不足6 ．5 0 ％，且各频带主频的相对误差均 保持在1 0 ％以内；幅值方面，实测值与预测值间的 相对误差仅2 ．1 3 ％，且主振频带 介于0 6 2 5H z 间 的预测精度比高频带高；此外，预测结果中虽有 个别频带的P P V 预测误差偏大 达1 7 ．1 0 ％ ，但整 合信号的可靠度较高 误差仅2 ．1 3 ％ ，这可避免单 一衰减因子情况下因一次的预测误差大而导致预测 结果失真的现象。因此，多频带小波包系数预测方 法较单一衰减因子分析方法的预测效果稳定。测点 振动信号的预测波形效果图见图2 。 综合表2 、表3 可以看出，利用修正萨氏多频带 小波包系数预测模型得到的P P V 与单一衰减因子 情况下由修正公式得到的P P V 几乎与真实值相等 误差均在2 ．1 5 ％以内 ，但前者得到了爆破振动信 号的整个波形，更有利于研究人员全面地把握爆破 振动信号各要素 P P V 、主频和振动持时 的综合作 用机制，因此前者较后者先进。 p l 帕 ● 暑 X t ／s 实测信号 【．1 L。【 。L 【Lk f7 盯’r 叮ff 旷 f ，B b 预测信号波形图 图2 实测信号和预测信号波形图 F 逅．2W a v e f o r mo ft h ea c t u a ls i g n a la n dt h e p r e d i c t e ds i g r l 且l 3结语 围绕修正萨道夫斯基公式和多频带小波包系 万方数据 第2 8 卷第1 期燕永峰，陈士海，张安康，等爆破振动信号波形预测模型的优化研究 2 7 数，利用M a t l a b 语言编程，提出了一种新型的爆破 振动信号波形预测模型。结合工程实例对比分析， 发现本波形预测模型的预测精度较高 幅值预测误 差仅2 ．1 3 ％，主频预测误差不足6 ．5 0 ％ ，预测信号 中个别频带的预测误差较大不一定导致最终预测信 号失真，这较单一衰减因子分析方法稳定，预测结果 可靠度高；此外，非线性分析求解的修正公式的预测 精度明显优越于传统公式，且由修正公式拟合得到的 装药量和爆心距之间的指数比并非常规采用的一1 ／3 ， 建议设计人员对特定场地进行爆破预测时，不要局限 于一1 ／3 的固定模式，应利用修正公式回归出适合本 爆破场地振动衰减规律的最佳指数比。 参考文献 R e f e r e n c e s [ 1 ]熊代余，顾毅成．岩石爆破理论与技术新进展[ M ] ．北 京冶金工业出版社，2 0 0 2 ． [ 2 ] 吴从师，吴其苏．爆破地震模拟初探[ J ] ．爆炸与冲击， 1 9 9 0 ，1 0 2 1 7 0 - 1 7 5 ． [ 3 ] 徐全军，张庆明，恽寿榕．爆破地震峰值的神经网络预 报模型[ J ] ．北京理工大学学报，1 9 9 8 ，1 8 1 4 4 7 2 4 7 5 ． [ 4 ]刘军，吴从师，高全臣．建筑结构对爆破震动的响应 预测[ J ] ．爆炸与冲击，2 0 0 2 ，2 0 4 3 3 3 - 3 3 7 ． [ 5 ] D O W D I N GCH ．M o n i t o r i n ga n dc o n t r o lo fb l a s t i n ge f f e c t s ． M i n i n gE n g i n e e r i n gH a n d b o o k [ M ] ．U p p e rS a d d l eR i v e r ， N e wJ e r s e y P r e n t i c eH a l l ，1 9 8 5 7 4 6 - 7 6 1 ． [ 6 ]陈士海，魏海霞，杜荣强．爆破震动信号的多分辨小波 分析[ J ] ．岩土力学，2 0 0 9 ，3 0 S 1 1 3 5 ．1 4 3 ． [ 7 ]毕明芽，李名山，刘朝红，等．爆破地震预测误差的因 素分析[ J ] ．爆破，2 0 0 9 ，2 6 2 9 6 - 9 8 ． [ 8 ] 高成．M a t l a b 小波分析与应用[ M ] ．北京国防工业 出版社，2 0 0 7 ． [ 9 ]毕卫国，石崇．爆破振动速度衰减公式的优化选择 [ J ] ．岩土力学，2 0 0 4 ，2 5 s 1 9 9 ．1 0 2 ． [ 1 0 ] 杜汉清．爆破振动衰减规律的现场试验研究[ J ] ．爆 破，2 0 0 7 ，2 4 2 1 7 2 0 ． [ 1 1 ] 林丽群，林从谋，蔡丽光，等．隧道爆破振动的B P 神 经网络预测及应用研究[ J ] ．爆破，2 0 0 9 ，2 6 4 4 ．7 ． [ 1 2 ]罗正，陈寿如．广义的爆破震动速度计算公式及其 应用[ J ] ．采矿技术，2 0 0 7 ，7 3 1 4 7 1 4 8 ． [ 1 3 ] 董陇军，李夕兵，赵国彦，等．爆炸冲击波信号处理方 法比较[ J ] ．岩石力学与工程学报，2 0 0 9 ，2 8 4 7 5 0 7 5 6 ． 上接第2 2 页 [ 2 ] 邓国强，杨秀敏．S P H 方法在爆炸冲击效应计算中的 应用[ J ] ．防护工程，2 0 0 4 ，2 6 6 4 6 - 4 9 ． [ 3 ] 赵海鸥．1 ．S D Y N A 动力分析指南[ M ] ．北京兵器工业 出版社，2 0 0 3 ． [ 4 ]邓国强，周早生，杨秀敏．爆炸冲击效应数值仿真中的 几项关键技术[ J ] ．系统仿真学报，2 0 0 5 ，1 7 5 1 0 5 9 1 0 6 2 ． [ 5 ] 杨冬梅，王晓鸣．混凝土中爆炸数值仿真算法研究 [ J ] ．爆炸与冲击，2 0 0 5 ，2 5 6 5 6 9 - 5 7 2 ． [ 6 ] 张雄，陆明万，王建军．任意拉格郎1 3 - 欧拉描述法研 究进展[ J ] ．计算力学学报，1 9 9 7 ，1 4 1 9 1 - 9 9 ． [ 7 ] 孙晓艳，王军．S P H 方法的理论及应用[ J ] ．水利水 电技术，2 0 0 7 ，3 8 3 4 4 4 6 ． [ 8 ] G e n e v i eT o u s s a i n t ，R O B E R TD u r o c h e r ．F i n i t ee l e m e n t s i m u l a t i o nu s i n gS P Hp a r t i c l e sf i t 8l o a d i n go nt y p i c a ll i g h t a n n o u r e dv e h i c l e s [ c ] ∥1 0I n t e r n a t i o n a lL S ．D Y N AU s e r s C o n f e r e n c e ，2 0 0 8 ． [ 9 ] L S ．D Y N Av e r s i o n9 7 1k e y w o r du s e r ’8m a n u a l [ M ] ．L i v e ． m o r eS o f t w a r eT e c h n o l o g yC o r p o r a t i o n ，2 0 0 7 ． [ 1 0 ] W A N GJ ．S i m u l a t i o no fl a n d m i n ee x p l o s i o nu s i n gL S - D Y N A 3 Ds o f t w a r e b e n c h m a r kw o r ko fs i m u l a t i o no fe x p l 0 6 i o ni n s o i la n da i r [ M ] ．F i s h e r m a n sB e n d ，V i c t o r i a ，A u s t F a l i a _ ． D S r I DA e r o n a u t i c a la n dM a r i t i m e R e s e a r c hI ．a 1 o r a t o r y ． 2 0 0 1 ． [ 1 1 ]白金泽．L 5 一D Y N A3 D 理论基础与实例分析[ M ] ．北 京科学出版社，2 0 0 5 ． [ 1 2 ] H O L M Q U I S T TJ ，J O H N S O NGR ，C O O KWH ．Ac o m - p u t a f i o n a lc o n s t i t u t i v em o d e lf o rc o n c r e t es u b j e c t e dt o l a r g es t r a i n s ，h i g hs w a i nr a t e sa n dh i g hp r e s s u r e s [ C ] ／／ 1 4 t hI n t e r n a t i o n a lS y m p o s i u mo nB a l l i s t i c s ，1 9 9 5 5 9 1 ．6 0 0 ． [ 1 2 ]龚自明，方秦．混凝土靶体尺寸及边界约束对侵彻 深度影响的数值研究[ J ] ．防护工程，2 0 0 4 ，2 6 6 5 0 - 5 7 ． [ 1 3 ]郑全平．钢筋混凝土板爆炸震塌试验研究与数值模 拟[ D ] ．南京中国人民解放军理工大学，2 0 0 3 ． 万方数据