空洞探测(一等奖)

摘要
本文以介质均匀的平板为研究对象，对平板内部的空洞位置进行探究。

建模过程中，先基于对数据的统计分析，确定出测量误差，用小的计算量找出了空洞的位置及洞径，另外，我们在确定空洞位置的基础上，对如何减少波源和接收器进行了分析，主要内容如下：
针对问题一，运用几何排除法，确定空洞在网格中的位置，首先根据洞径长度求解洞径的个数，画几何分布图，连接无空洞的直线，得到结论：洞在网格中的具体位置为：2223,253233344453,,,,,,a a a a a a a a ，并在此基础上，计算洞径，画出洞径在小方格内的大致位置，空洞呈现椭圆形。

针对问题二，分析出了只由出i P 的弹性波到达j Q 时间(),=1,,7ij t i j ⋅⋅⋅，不能确定空洞位置的两点原因：
1.平行于PQ 连线方向的若干个空洞无法确定其位置。

2.发现只存在i P 面的波源与j Q 面的接收器存在一定的“死区”无法测量
针对问题三，在问题二的基础上通过增加的i R 面的波源与j S 面的接收器来避免“死区”，同理增加的i P 面的波源与j Q 面的接收器个数。

综合分析增加的波源的传播路径，最终得出最优方案：每边减少为5对波源与接收器。

关键词：空洞位置 几何排除法 探测“死区”
一、问题重述
1.1问题背景与条件
本题目来源于吉林丰满水电站水库大坝的检测与维修，这种类型的大坝每隔一定时间就需进行较全面的维修。

为此, 首先就要对坝的内部结在维修利用山体、隧洞及坝体时，人们必须考虑可能由于微生物、地质等因素引起山体等的内部结构变化而产生的不良后果，例如空洞的产生。

因而，在维修之前，必须运用一定的技术对山体进行空洞的探测，了解内部的缺陷情况，而对大坝的检测运用的就是弹性波的断层技术，根据弹性波在不同介质中的传播速度不同来确定空洞的位置。

通过接收器接收的数据来进行分析，但由于数据量大，必须由计算机来完成，由此我们需要对这类空洞探测问题建立合乎实际的模型及算法，所以本题以介质均匀的平板为研究对象，探究空洞的位置。

1.2需要解决的问题
本题将山体、隧洞及坝体等大型物体简化为一块均匀介质构成的矩形平板进行研究，在平板内有着一些充满空气的空洞，在平板的两个邻边分别等距地设置若干波源，在它们的对边对等地安放同样多的接收器，记录弹性波由每个波源到达对边上每个接收器的时间，根据弹性波在介质中和在空气中不同的传播速度，来确定板内空洞的位置。

现考察如下的具体问题：
一块240（米）⨯240（米）的平板（如图），在 AB 边等距地设置7个波源，
()=1
,,7i P i ⋅⋅⋅，CD 边对等地安放7个接收器()=1,,7j Q j ⋅⋅⋅，记录由i P 发出的弹性波到达j Q 的时间ij t （秒); 在AD 边等距地设置7个波源()=1,,7i R i ⋅⋅⋅，BC 边对等地安放7个接收器()=1,,7j S j ⋅⋅⋅，录由i R 发出的弹性波到达j S 的时间ij r （秒），弹性波在介质和在空中的传播速度分别为2880（米/秒）和320（米/秒），且弹性波沿板边缘的传播速度与在
介质中的传播速度相同。

1）确定该平板内空洞的位置。

2）只根据由出i P 的弹性波到达j Q 时间(),=1,,7ij t i j ⋅⋅⋅，确定空洞的位置吗；讨论在同样能够确定空洞位置的前提下，减少波源和接受器的方法。

二、模型假设
结合本题的实际，为了确保模型求解的准确性和合理性，我们排除了一些位置因素的干扰，提出以下几点假设：
1、弹性波在介质和空气中的传播是沿直线传播的，在传播过程中互不干扰且不发生折射与衍射。

2、题中所给的对于弹性波传播时间的数据真实可信，但存在一定的误差。

3、只求空洞的大概位置，形状近似认为椭圆形。

4、所有的探测均在同一平面内进行，所谓的空洞即平面上的一个区域。

5、误差允许范围内的空洞忽略不计。

三、符号说明及名词解释
为了便于问题的求解，我们给出以下符号说明：（其他未说明的符号在文中第一次出
四、问题分析
利用弹性波进行空洞位置的探测，是根据弹性波在介质及空气中的传播速度不同而确定的。

对于240（米）⨯240（米）的平板（如图），在 AB 边等距地设置7个波源，
()=1
,,7i P i ⋅⋅⋅，CD 边对等地安放7个接收器()=1,,7j Q j ⋅⋅⋅，记录由i P 发出的弹性波到达j Q 的时间ij t （秒); 在AD 边等距地设置7个波源()=1,,7i R i ⋅⋅⋅，本题的目的就是根据接收器测得的时间，求出空洞的具体位置及探讨减少波源和接收器的方法。

对于问题一，分析题目中所给的要求和条件，发现要想直接根据题目中的时间测得空洞的位置是十分困难的，鉴于本题中已将平板分为66⨯网格形式，对此我们将充分利用。

根据题中的条件 ，“弹性波沿板边缘的传播速度与在介质中的传播速度相同”，即弹性波在沿板边缘传播时的实际值应该与理论值相同，都为2402880/=0.0833m m s s ÷，但题中表格提供的实际测量值从=0.0583ij T s 到=0.1024ij t s 不等，这说明了在测量时存在误差，最大误差为：0.0833-0.0583=0.025s s s ，在用最大误差时间乘以弹性波在空气中的传播速度，得出结论：由于测量误差的限制，弹性波只能测量由传播方向上洞径大于8m 的空洞。

然后根据弹性波在纯介质中的传播时间和弹性波的传播距离，利用公式
100
2
212
-==-ij i j i j ij v t PQ PQ t v v v v ⨯⨯ 求出空洞的长度，并将洞径小于8m 的忽略不计，根据表三中数据，于是得出结论：空洞位于36个网格的内部，在利用空洞长度和弹性波的总长度的百分比计算每条弹性波穿过空洞的个数，运用几何排除法，连接没有空洞的波线，余下即为空洞在网格中的具体
位置。

由于空洞只是位于网格内，但并不一定占满网格，为了更准确，本文通过对洞径求解算数平均值来确定其在以40m 为边长的正方形中的大致分布。

对于问题二，根据问题一的求解过程可以知道想确定空洞的位置必须有横纵两种波线，另外，通过波线传播的路径图，直观的看出在平板的周边区域，还是有很大的空隙，并没有被路径覆盖，不能探测出在周边是否有空洞存在，即“死区”无法测量。

并且对“死区”面积进行计算，并添加少量的波源使“死区”内的空洞得以被探测。

对于问题三，在问题二的基础上进行探讨，必须使弹性波能探测死区面积在RS 上添加相应必须的波源与接收器弥补问题二的缺陷，按照同样的方法可以确定在PQ 平面上必须存在的波源与接收器的个数与位置，达到用最少的波源与接收器基本覆盖所有方格且没有死区的目的。

五、模型的建立与求解
经过以上的分析和准备，我们将逐步建立以下数学模型，进一步阐述模型的实际建立过程。

5.1问题一的模型建立与求解 5.1.1空洞位置的确定
依据波源和接收器的位置，已经将平板划分为36个方格区域，而弹性波总有98种传播路径，为求空洞的位置，可以通过弹性波的总长度和空洞的总长度求出每种传播路径上空洞的数量，运用向上取整的方法确定占用的网格数，画出几何分布图，连接没有经过空洞的波线，余下即为空洞在网格中的位置。

1step ：空洞总长度的计算
本题将一块240（米）⨯240（米）的平板，在AB 和AD 边等距地设置7个波源，()=1
,,7i P i ⋅⋅⋅，分别在CD 和BC 边对等地安放7个接收器()=1,,7j Q j ⋅⋅⋅，正好将平板划分为36个小正方形，每个正方形的边长为40m ，由此可以求出弹性波在平板中的传播路径的长度为：
=R i j i j PQ S
运用Mathematica 编程（见附录一）得到弹性波从i P 到j Q 的传播距离见表1。

表1 弹性波从i P 到j Q 的传播距离
i j 先假设该平板介质均匀，无空洞，则弹性波在平板中的传播速度是固定值，为
2880/m s ，在结合上表1，根据公式：
1
=
i j
ij PQ t v
求出弹性波在纯介质中传播时间的理论值，运用Mathematica 编程（见附录一）得到结果见下表2。

表2 从波源i P 的弹性波到达j Q 的理论时间
i j 观察数据，在结合题目中给出的条件，“弹性波沿板边缘的传播速度与在介质中的传播速度相同”，即弹性波在沿板边缘传播时的实际值应该与理论值相同，都为2402880/=0.0833m m s s ÷，观察表三和题目中的表格，发现并不相同，实际测量值从=0.0583ij T s 到=0.1024ij t s 不等，这说明了在测量时存在误差，最大误差为：0.0833-0.0583=0.025s s s ，
在用最大误差时间乘以弹性波在空气中的传播速度，得出结论：由于测量误差的限制，弹性波只能测量由传播方向上洞径大于8m 的空洞。

空洞总长度002=i j ij PQ t v ⨯，00
2=i j ij R S T v ⨯，题目中已经给出了空气的传播速度，只需求出空洞传播的时间，对于空洞存在的传播路径中，其总路程可以划分为，介质传播路径和空气的传播路径，两者相加即为总路径的长度，公式为
()0021=t +-i j ij ij ij PQ v t t v ，()00
21=+i j ij ij ij R S T v T T v -
即10
12
-=
-ij i j
ij v t PQ t v v ， 1012-=
-ij i j
ij v T R S T v v
又结合题中已知的2=320/v m s ，所以
1002212-==
-ij i j i j ij v t PQ PQ t v v v v ，100
2212
-==-ij i j i j ij v T R S R S T v v v v 运用Mathematica 编程（见附录二）得到P 到Q 空洞的长度见下表3。

表3 P 到Q 空洞的长度
注：R 到S 传播路径上的空洞长度见附录三。

观察表中的数值，由于误差测量的限制，可以将沿传播方向上洞径小于8m 的路径直接忽略，结合网格图可知：i j PQ 和i j R S 的传播路径上均无空洞，即空洞的位置位于36个网格内。

2step ：求解空洞数量
在66⨯的网格中，要求解空洞的数量，即空洞在网格中所占用的格数，而空洞的总长度和弹性波的总长度已知，可以用空洞总长度与弹性波的总长度的比值乘以6即为空洞占网格的个数，但由于本题要求的解空洞的数量必须为整数，所以要对求解的数值向上取整，具体公式如下：
()0
=6
+0.5i j
i j i j PQ N PQ PQ ⎡⎤⎢⎥⎢⎥⎣⎦，()0=6+0.5i j i j i j R S N R S R S ⎡⎤⎢⎥⎢⎥⎣⎦
根据公式运用Matlab 编程（见附录四）求出空洞数量见下表4
i j 3step ：确定空洞的具体位置
观察上表，很多弹性波经过的空洞数量都大于1，不易直接观察空洞的具体位置，所以采用几何排除法，先出几何分布图，连接不穿过空洞的波线，余下即为空洞在网格中所处的位置，因此将表格中个数为0的线连接起来，剩下的区域即为空洞区域，见图1。

图1 空洞的位置
所以空洞在网格中的具体位置为：2223,253233344453,,,,,,a a a a a a a a
5.1.2 空洞直径的估计
根据图2，可以粗略的看出洞径在网格的位置，但是洞径并非占满整个网格，为了是结果更加精确，我们计算洞径的大小，大致画出洞径的示意图，使结果更加直观。

对于求解空洞56a 的洞径：
可以求出弹性波56R S 经过的空洞距离和65R S 经过的空洞距离，取两者的平均值，减少误差。

计算如下：
056R S =0
15656
562212-=
=37.05-v T R S T v v v v
065R S =0
16565652212
-==43.67-v T R S T v v v v
所以 00
566525+=40.362
R S R S d =
即56a 的洞径为40.36m 同样对于存在空洞的每行、每列上的波线长度运用取平均值的方法，可以得到8个方程：
0000
5665566525530000
544545544444340000
23323443
223232333400233223335325222++==40.3622==41.3162
22++==36.2942+==86.5502
22
++==78.6112++=117.06
22
++=119.76++2R S R S P Q P Q d d P Q P Q R S R S d d d R S R S PQ P Q d d d d d R S S R d d d d d d +=+=00
2332
3=123.378
2
P Q PQ += 求解方程组得到如下结果见表5
图2 空洞的具体位置
5.2问题二的模型建立与求解
有如下2点理由说明只根据由出i P 的弹性波到达j Q 时间(),=1,,7ij t i j ⋅⋅⋅，无法确定空洞的位置。

1、问题一的求解过程，可以看出要求得空洞的位置，需要根据空洞的总长度与弹性波的总长度求每条弹性波经过的空洞数量，但是要得出空洞在网格中的具体位置，仅根据i j PQ 上的弹性波经过的空洞个数是无法确定的，因为如果存在平行于PQ 连线方向的若干个空洞, 其总长接近240米, 且空洞的洞径在沿垂直于PQ 的方向上接近相等，则因这些洞引起的由各i P 发出的弹性波到达所对应接受器的时间延迟都近似相同，此时便不能确定空洞存在的上下位置。

2、由对称性分析, 发现“死区”主要有2部分，第一“死区”主要集中在探测器的周围，若在其周围存在直径大于8m 的空洞，且弹性波无法涉及到，则产生了一部分的“死区”；其次，在直线11PQ 和77P Q 中点附近的20m 处，也是测量仪的死区，在平板的周边区域，还是有很大的空隙，并没有被路径覆盖，不能探测出在周边是否有空洞存在。

基于以上两点理由可以确由如果只在PQ 方向上有波源和接收器，是不能具体求出的空洞区域。

下面计算“死区”面积：
（1）探测器周围的“死区”面积：
过点234,,P P P 分别作172737,,PQ P Q PQ 的垂线，观察下图3
图3 探测器周围“死区”面积
得出总共有12个“死区”，首先计算12P P 间的“死区”面积即求解121S PP h ∆，
140sin(/4)28.3h m π==
所以2211
400()2
S h m =≈
对于其他“死区”，与第一个“死区”相比较，是等底不等高的三角形，所以先求出其余“死区”三角形的高
272340sin h Q P P =∠ 373440sin h Q P P =∠
因为7347234
Q P P Q P P π
∠>∠>
，所以32128(h h h >>>米)。

又由问题一求解过程中的
误差分析可知，洞径大于8m 的均能被探测的，但在问题二中由于只有i P 到j Q 的弹性波
导致了存在盲区，每一个探测器周围的死区的大概面积为：2400m （2）平板周边区域的“死区”面积
观察弹性波从i P 到j Q 的传播路径，只做出平板周边区域的传播路径，见下图4：11PQ 和77P Q 的中点直线附近大约20m 处，存在“死区”，而平板两边的死区位置及面积是相对称的，“死区”总面积26s=240020m m =（）,其中h
图4 平板周边的“死区”
5.3.探究减少波源和接收器的方法
结合上一问的结论，仅根据i P 到j Q 的传播路径是无法确定空洞位置的，忽略了盲区，在这一问中，要求确定空洞位置，并减少波源和接收器，也就是尽可能少地增加波源和接收器使“死区”能够被探测到。

（1）在仅有i P 到j Q 方向的探测器的情况下，分别在()11R Q 、12R R 中点、4R 、67R R 中
点和()71R P 处安装五个探测器，同时在其对称处77()S Q 、67S S 中点、4S 和()17S P 处安装
五个接收器，其波线的传播路径见下图5。

图5 增加波源与接收器 图6 增加波源与接收器效果图
观察图6，发现增加这五对波线，探测器周围的“死区”都将能被探测到，由于在77
Q S 处存在两个接收器，在17P R 处也同时存在两个波源，为了使所用的波源和探测器尽可能少，将两组中多余的一个去除，所以最终在()11R Q 、12R R 中点、4R 、67R R 中点处安装四个探测器，同时在其对称处77S Q 、67S S 中点、4S 处安装四个接收器。

（2）同理在仅有i R 到j S 方向的探测器的情况下,与PQ 方向的波线一样，也会导致两部分的死区，为了消除“死区”，需要适当地在PQ 方向增加波源和探测器的数量，分别在()17P R 、12P P 中点、4P 、67P P 中点和()71P S 处安装五个探测器，同时在其对称处
11()Q R 、12Q Q 中点、4Q 、67Q Q 和77Q S 处安装五个接收器,同样将()17P R ，()77Q S 处重
复安装的波源和探测器去除，所以最终在12P P 中点、4P 、67P P 中点和()71P S 处各安装一个探测器，同时在其对称处11()Q R 、12Q Q 中点、4Q 、67Q Q 和处安装四个接收器。

（3）综合以上增加的波源和探测器，做出新波源和接收器传播路径效果图见下图7。

图7 改良后弹性波的传播路径
由图7中波线的传播路径可知，该平板内，“死区”均被探测到，所以可以将波源与接收器的组数减少为5组，根据问题一的空洞探测模型同样能够确定空洞位置。

六、模型的评价与改进
6.1优点
1．问题一采用几何排除法，将弹性波穿过空洞个数为0的路径划出，余下即为空洞的位置，此方法比较直观，易于计算。

2．对于问题一，本文不仅求出空洞在网格中的位置，还详细求出洞径，画出空洞的分布图，使结果更加明确，直观。

3.本文多处将文字与图表相结合，使结果直观，便于研究。

4.本文使用AutoCAD对图形进行处理，图形准确且清晰，使模型更加精确，可信。

6.2模型的缺点
1.该模型的建立依赖于网格的分布，在实际生活中，山体或隧道的内部空洞形状不规则，数量较多，所以该模型不具有普遍性。

2.只能求解出空洞的大体位置和洞径，精确度较低，不适用于密集、复杂的空洞探测
但两方案均依赖于平板的分区, 若空洞的分布跨越多个方格区域, 则应寻求更为普遍的
模型与算法.答卷是用排除法(或称剔除法) 求解的. 针对此问题的特殊情况, 此方法
是有效的, 也是较简单的; 但有一定的局限性, 而不具备较普遍的适用性.
6.3模型的改进
本模型提供了确定空洞存在、数量和位置易行的方法，相关管理人员可以通过适时观测发现隐患，及时排除。

但由于文字数据全是人工拟定的，不可避免地与实际情形有一定的差距和误差，在实际操作时应尽量提高精确度。

本模型是在允许偏差的情况下，根据7个接收器测得的时间建立模型，从而确定空洞位置，在实际测量中，对于不同的接收器来说，其灵敏度、大小及系统本身的误差是不同的，所以为了是结果更准确，应该对每个接收器做单独的误差分析。

七、模型推广
(1)该模型可以适用于精度较低的空洞探测，在精度不高的空洞探测问题中，能够很好的降低成本，达到探测出空洞的目的。

合理安排波源和接收器的位置，不仅能够使平板被完全覆盖，同时也能够减少机器的使用，降低能耗，降低成本，减少探测使用的时间，达到最好的效果。

但此模型不适用于精度较高的空洞探测，因为考虑到存在误差的问题，以及在连线处也可能存在空洞的问题，应该寻找更精细，也更合理的模型进行求解。

(2)将平板按小格划分时，可以将小格数目进一步扩大，用逐步求精的方法细化空洞的具体位置，使空洞位置更清晰。

在实际生活中我们也仅仅只能估算出空洞的大概位置, 然后再用其它方法确定其比较准确位置。

八、参考文献
[1] 邹建光，刘丹娟，曾德刚，空洞探测，数学的实践与知识，2001年1月；
[2] 杜秀丽，朱群生，王建锋，探测空洞的方法，数学的实践与认识，2001年3月；
[3]姜启源等编，数学模型（第三版），高等教育出版社，2003年08月；
[4]张德丰等编，MATLAB高级语言编程，机械工业出版社，2010年01月
九、 附录
附录一：弹性波从
i
P 到
j
Q 的传播距离
Export["F:\\p 到q 的距离和r 到S 的距离.xls", Table[40*Sqrt[6^2 + (j - i)^2], {i, 7}, {j, 7}] // N] Export["F:\\无洞情况下的时间.xls", Table[40*Sqrt[6^2 + (j - i)^2]/2880, {i, 7}, {j, 7}] // N]
附录二：P 到Q 和R 到S 传播路径上的空洞长度
spq = Table[40*Sqrt[6^2 + (j - i)^2], {i, 7}, {j, 7}] // N; tpq = Flatten[Import["F:\\p 到q 的时间.xls"], 1]; spq0 = (tpq*2880 - spq)/(2880 - 320)*320;
Export["F:\\p 到q 穿越空洞的时间.xls", (tpq*2880 - spq)/(2880 - 320)] Export["F:\\p 到q 穿越空洞总长度.xls", spq0]
srs = Table[40*Sqrt[6^2 + (j - i)^2], {i, 7}, {j, 7}] // N; trs = Flatten[Import["F:\\r 到s 的时间.xls"], 1]; srs0 = (trs*2880 - srs)/(2880 - 320)*320;
Export["F:\\r 到s 穿越空洞的时间.xls", (trs*2880 - srs)/(2880 - 320)]
Export["F:\\r 到s 穿越空洞总长度.xls", srs0]Export["F:\\r s .xls",srs0]
附录四：P 到Q 和R 到S 传播路径上空洞数量
clc
clear all
load pqkdjuli.txt %导入数据 load rskdjuli.txt %导入数据 load pqrsjuli.txt %导入数据 a=rskdjuli() b=pqkdjuli() c=pqrsjuli()
x=fix(6*b./c+0.5) %取整 y=fix(6*a./c+0.5)
xlswrite('pq 穿越空洞的平均个数', x) %导出数据 xlswrite('rs 穿越空洞的平均个数', y) %导出数据
附录五：R 到S 空洞空洞的数量
附录六：
Matlab画图：
clc
clear all
close all
a=linspace(0,2*pi,100)
x1=(40+45.176/2)+45.176/2*cos(a); y1=(40+45.176/2)+45.176/2*sin(a); plot(x1,y1)
hold on %使下一张图不被覆盖x2=(40+37.8/2)+37.8/2*cos(a);
y2=(80+37.8/2)+37.8/2*sin(a);
plot(x2,y2)
hold on
x3=(40+40.3/2)+40.3/2*cos(a);
y3=(160+40.3/2)+40.3/2*sin(a); plot(x3,y3)
hold on
x4=(80+33.4/2)+33.4/2*cos(a);
y4=(40+33.4/2)+33.4/2*sin(a);
plot(x4,y4)
hold on
x5=(80+33.3/2)+33.3/2*cos(a);
y5=(80+33.3/2)+33.3/2*sin(a);
plot(x5,y5)
hold on
x6=(80+50.2/2)+50.2/2*cos(a);
y6=(120+50.2/2)+50.2/2*sin(a); plot(x6,y6)
hold on
x7=(120+36.29/2)+36.29/2*cos(a); y7=(120+36.29/2)+36.29/2*sin(a);
plot(x7,y7)
hold on
x8=(160+41.3/2)+41.3/2*cos(a); y8=(80+41.3/2)+41.3/2*sin(a); plot(x8,y8)
axis([1 240 1 240]) %范围xlabel('p波源')
text(-30,80,'r波源')
text(108,245,'q接收器') legend('空洞位置')
text(240,80,'s接受器')。