CN104020508A - 一种用于地质雷达层析探测的直射线追踪算法 - Google Patents
一种用于地质雷达层析探测的直射线追踪算法 Download PDFInfo
- Publication number
- CN104020508A CN104020508A CN201410265945.6A CN201410265945A CN104020508A CN 104020508 A CN104020508 A CN 104020508A CN 201410265945 A CN201410265945 A CN 201410265945A CN 104020508 A CN104020508 A CN 104020508A
- Authority
- CN
- China
- Prior art keywords
- ray
- grid
- length
- coordinate
- section
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000004587 chromatography analysis Methods 0.000 title claims abstract description 16
- 238000000034 method Methods 0.000 claims abstract description 35
- 238000004364 calculation method Methods 0.000 claims description 14
- 238000001514 detection method Methods 0.000 claims description 12
- 239000011159 matrix material Substances 0.000 claims description 12
- 238000005516 engineering process Methods 0.000 description 4
- 238000010586 diagram Methods 0.000 description 2
- 238000011835 investigation Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000005065 mining Methods 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
Landscapes
- Image Generation (AREA)
Abstract
本发明涉及一种用于地质雷达层析探测的直射线追踪算法，属于地球物理反演算法领域，该方法包括如下步骤：建立平面直角相对坐标系；对剖面进行网格划分与编号；确定各网格的左右标点坐标；确定发射点与接收点的坐标；计算射线参数；计算射线在各网格内的长度；检验计算结果；将计算结果进行稀疏存储。本发明根据射线的左右截距与斜率三个参数追踪射线走向，只需计算射线经过的网格，相比传统算法遍历所有网格，大大降低了对内存和时间的消耗，提高了算法的效率与精度，最后将计算结果以稀疏格式存储，可节省存储空间并直观显示射线的分布状况。
Description
技术领域
本发明属于地球物理反演算法领域，涉及一种用于地质雷达层析探测的直射线追踪算法。
背景技术
地质雷达技术是一种快速、无损、高效的地球物理探测技术，近年来在矿井的应用逐渐增多，如煤矿含水层、陷落柱、裂隙带和瓦斯突出等灾害隐患的探测。目前地质雷达主要采用反射回波技术，探测深度局限于30~50米。由于现代化矿井回采区巷道之间的距离往往大于200米，使得地质雷达在探测采煤工作面内灾害源时无法达到所需要的探测深度。
层析探测技术是地质雷达领域一个新的方向，通过采集直达波的走时或能量进行反演计算，得到探测剖面的地质结构信息，探测距离可达到反射回波技术的2倍以上。反演计算的实质是求解一个矩阵方程组，其系数矩阵为射线在剖面离散化后的网格模型中各个网格内的长度。因此，射线追踪结果的精度与反演计算的精度密切相关。传统的射线追踪算法是针对每一条射线，均遍历整个网格模型，由于网格数目众多，并且每条射线经过的网格个数占整个模型的网格个数的比例很小，射线在大多数网格内的长度均为0，导致算法不仅存储量大、计算耗时长，并且效率很低。
发明内容
本发明所要解决的技术问题是针对现有射线追踪算法存储量大、计算耗时长，并且效率很低的缺点而提供一种用于地质雷达层析探测的直射线追踪算法。
本发明解决其技术问题采用的技术方案是：一种用于地质雷达层析探测的直射线追踪算法，具体包括：
S1、 建立平面直角相对坐标系：将地质雷达探测到的剖面全部置于坐标系的第一象限，其中层析探测方向平行于x轴，测线方向平行于y轴；坐标轴原点位于剖面顶点；层析探测垂直距离为L x ，剖面的测线长度为L y ，则剖面的四个端点的坐标分别为(0, 0)，(L x , 0)，(0, L y )，(L x , L y )；
S2、 对剖面进行网格划分与编号：在剖面内绘制边长为D的正方形网格，以坐标轴原点处的网格为起始，编号为1；进一步的以x轴方向为行，y轴方向为列，逐行依次编号，从而将剖面划分为N个网格；
S3、 确定各网格的左右标点坐标：将网格4个顶点中(x min , y min )的顶点确定为左标点，(x min , y max )的顶点确定为右标点，第i个网格的左标点的坐标记为(lx i , ly i )，第i个网格的右标点的坐标记为(rx i , ry i )，i=1,2,…,N；
S4、 确定发射点与接收点的坐标：根据发射点、接收点与探测剖面的相对位置，确定发射点、接收点在坐标系中的位置，发射点坐标记为（Tr x , Tr y ），接收点坐标记为(Re x , Re y )；进一步的定位射线的发射点所在的网格编号n 1 与接收点所在的网格编号n 2 ；
S5、 计算射线参数：k 1 = (Re y – Tr y ) / (Re x – Tr x )，k 2 = ，k 3 = k 2 /k 1 ；若k 1 = 0，则无需计算k 3 ；
S6、 计算射线在各网格内的长度：记为l i ，i=1,2,…,N；
S7、 检验计算结果：根据公式：计算出射线的实际长度L 1 ；根据公式：计算出射线的计算长度L 2 ；根据公式计算出射线追踪的误差?；如果?> 1.0E-04，检查之前步骤；
S8、 将计算结果进行稀疏存储：将每一条射线的计算结果存储为矩阵；矩阵行数为射线经过的网格个数，矩阵的第一列为射线经过的网格的编号，第二列为射线在相应网格内的长度。
进一步的，所述步骤S1中的剖面为二维矩形剖面。
进一步的，所述步骤S4中的发射点的横坐标均为0，接收点的横坐标均为L x 。
进一步的，所述步骤S4中定位发射点所在的网格编号n 1 的方法为：
,
进一步的，所述步骤S4中定位接收点所在的网格编号n 2 的方法为：
,
进一步的，所述步骤S6中计算k 1 = 0的射线在各网格内的长度的方法为：若n 1 ≤ n 2 ，射线在网格编号在[n 1 , n 2 ]范围的网格内的长度均为D；若n 2 < n 1 ，射线在网格编号在[n 2 , n 1 ]范围的网格内的长度均为D；其他网格内的长度均为0。
进一步的，所述步骤S6中计算k 1 ≠ 0的射线在各网格内的长度的方法为：
S61、 计算射线与网格n 1 左边界所在直线的交点的纵坐标y 1 , 根据公式：计算左截距d 1 ；
S62、 计算射线与网格n 1 右边界所在直线的交点的纵坐标y 2 , 根据公式：计算右截距d 2 ；
S63、 根据如下公式计算射线在网格n 1 内的长度：
若k 1 >0，
若k 1 <0，
S64、 根据如下公式判断射线经过的下一个网格的编号c：
若k 1 >0，
若k 1 <0，
S65、 重复以上过程，直到计算出射线在网格n 2 内的长度为止；
S66、 记射线在其他网格内的长度为0。
本发明的有益效果是：给出了一种用于地质雷达层析探测的直射线追踪算法，只需计算射线在经过的网格内的长度，算法步骤简便，大大降低了计算时间和内存需求，提高了算法的效率与精度；将计算结果以稀疏格式存储，既可以节省大量存储空间，又可直观的显示射线在探测剖面的分布状况。
附图说明
图1为本发明的算法流程图；
图2为本发明实施例的坐标系、网格模型及射线举例的示意图。
具体实施方式
下面将结合附图对本发明作进一步说明。
图1为本发明实施例的地质雷达层析探测的直射线追踪算法的流程图，其具体包括如下步骤：
S1、 建立平面直角相对坐标系：如图2所示，将地质雷达探测到的剖面全部置于坐标系的第一象限，其中层析探测方向平行于x轴，测线方向平行于y轴；坐标轴原点位于剖面顶点；层析探测垂直距离为200米，测线长度为250米，则剖面的四个端点的坐标分别为(0, 0)，(200, 0)，(0, 250)，(200, 250)；
S2、 对剖面进行网格划分与编号：在剖面内绘制边长为10米的正方形网格，以坐标轴原点处的网格为起始，编号为1；进一步的以x轴方向为行，y轴方向为列，逐行依次编号，从而将剖面划分为500个网格；
S3、 确定各网格的左右标点坐标：将网格4个顶点中(x min , y min )的顶点确定为左标点，(x min , y max )的顶点确定为右标点，第i个网格的左标点坐标记为(lx i , ly i )，第i个网格的右标点坐标记为(rx i , ry i )，i=1,2,…,500；
S4、 确定发射点与接收点的坐标：根据发射点、接收点与探测剖面的相对位置，确定发射点、接收点在坐标系中的位置，发射点坐标记为（Tr x , Tr y ），接收点坐标记为(Re x , Re y )；进一步的定位射线的发射点所在的网格编号n 1 与接收点所在的网格编号n 2 ；
所述步骤S4中定位发射点所在的网格编号n 1 的方法为：
,
所述步骤S4中定位接收点所在的网格编号n2的方法为：
,
举例说明：如图2所示，射线A：发射点坐标(0,125)，n 1 = 241，接收点坐标(200,125)，n 2 = 260；射线B：发射点坐标(0,125)，n 1 = 241，接收点坐标(200,15)，n 2 = 40；射线C：发射点坐标(0,125)，n 1 = 241，接收点坐标(200,245)，n 2 =500；
S5、 计算射线参数：k 1 = (Re y – Tr y ) / (Re x – Tr x )，k 2 = ，k 3 = k 2 /k 1 ；若k 1 = 0，则无需计算k 3 ；
举例说明： 射线A：k 1 = 0, k 2 = 1；射线B：k 1 = -0.55, k 2 = 1.14127, k 3 = -2.075；射线C：k 1 = 0.6 k 2 = 1.16619, k 3 = 1.94365；
S6、 计算射线在各网格内的长度：记为l i ，i = 1,2,…,N；
所述步骤S6中计算k 1 = 0的射线在各网格内的长度的方法为：若n 1 ≤ n 2 ，射线在网格编号在[n 1 , n 2 ]范围的网格内的长度均为D；若n 2 < n 1 ，射线在网格编号在[n 2 , n 1 ]范围的网格内的长度均为D；其他网格内的长度均为0。
举例说明：射线A，k 1 = 0, 因此在网格编号[241,260]范围的网格内的长度均为10，其他网格内的长度均为0；
所述步骤S6中计算k1≠0的射线在各网格内的长度的方法为：
S61、 计算射线与网格n1左边界所在直线的交点的纵坐标y1, 根据公式： 计算左截距d1；
S62、 计算射线与网格n1右边界所在直线的交点的纵坐标y2, 根据公式： 计算右截距d2；
S63、 根据如下公式计算射线在网格n1内的长度：
若k1>0，
若k1<0，
S64、 根据如下公式判断射线经过的下一个网格的编号c：
若k1>0，
若k1<0，
举例说明：射线B：d 1 = 5，d 2 = -0.5，= 10.375，c = 221；射线C：d 1 = 5，d 2 = 11，=9.718，c = 261；
S65、 重复以上过程，直到计算出射线在网格n2内的长度为止；
S66、 记射线在其他网格内的长度为0；
S7、 检验计算结果：根据公式：计算出射线的实际长度L1；根据公式：计算出射线的计算长度L2；根据公式计算出射线追踪的误差?；如果? > 1.0E-04，检查之前步骤；
举例说明：射线B：L 1 = 228.2542442，L 2 = 228.2542454，? = 5.03379E-09，?< 1.0E-04，满足精度要求；射线C：L 1 = 233.2380758，L 2 = 233.2380782，? = 1.05185E-08，? < 1.0E-04，满足精度要求；
S8、 将计算结果进行稀疏存储：将每一条射线的计算结果存储为矩阵；矩阵行数为射线经过的网格个数，矩阵的第一列为射线经过的网格的编号，第二列为射线在相应网格内的长度。
举例说明：射线A、B、C的计算结果分别存储为矩阵M A 、M B 、M C ，如下所示，
。
Claims (7)
1.一种用于地质雷达层析探测的直射线追踪算法，其特征在于，所述方法包含如下步骤：
建立平面直角相对坐标系：将地质雷达探测到的剖面全部置于坐标系的第一象限，其中层析探测方向平行于x轴，测线方向平行于y轴；坐标轴原点位于剖面顶点；层析探测垂直距离为L x ，剖面的测线长度为L y ，则剖面的四个端点的坐标分别为(0, 0)，(L x , 0)，(0, L y )，(L x , L y )；
S1、对剖面进行网格划分与编号：在剖面内绘制边长为D的正方形网格，以坐标轴原点处的网格为起始，编号为1；进一步的以x轴方向为行，y轴方向为列，逐行依次编号，从而将剖面划分为N个网格；
S2、确定各网格的左右标点坐标：将网格4个顶点中(x min , y min )的顶点确定为左标点，(x min , y max )的顶点确定为右标点，第i个网格的左标点的坐标记为(lx i , ly i )，第i个网格的右标点的坐标记为(rx i , ry i )，i=1,2,…,N；
S3、确定发射点与接收点的坐标：根据发射点、接收点与探测剖面的相对位置，确定发射点、接收点在坐标系中的位置，发射点坐标记为（Tr x , Tr y ），接收点坐标记为(Re x , Re y )；进一步的定位射线的发射点所在的网格编号n 1 与接收点所在的网格编号n 2 ；
S4、计算射线参数：k 1 = (Re y – Tr y ) / (Re x – Tr x )，k 2 = ，k 3 = k 2 /k 1 ；若k 1 = 0，则无需计算k 3 ；
S5、计算射线在各网格内的长度：记为l i ，i= 1,2,…,N；
S6、检验计算结果：根据公式：计算出射线的实际长度L 1 ；根据公式：计算出射线的计算长度L 2 ；根据公式计算出射线追踪的误差?；如果? > 1.0E-04，检查之前步骤；
S7、将计算结果进行稀疏存储：将每一条射线的计算结果存储为矩阵；矩阵行数为射线经过S8、的网格个数，矩阵的第一列为射线经过的网格的编号，第二列为射线在相应网格内的长度。
2.如权利要求1所述的方法，其特征在于，所述步骤S1中的剖面为二维矩形剖面。
3.如权利要求1所述的方法，其特征在于，所述步骤S4中的发射点的横坐标均为0，接收点的横坐标均为L x 。
4.如权利要求1所述的方法，其特征在于，所述步骤S4中定位发射点所在的网格编号n 1 的方法为：
, 。
5.如权利要求1所述的方法，其特征在于，所述步骤S4中定位接收点所在的网格编号n 2 的方法为：
, 。
6.如权利要求1所述的方法，其特征在于，所述步骤S6中计算k 1 = 0的射线在各网格内的长度的方法为：若n 1 ≤ n 2 ，射线在网格编号在[n 1 , n 2 ]范围的网格内的长度均为D；若n 2 < n 1 ，射线在网格编号在[n 2 , n 1 ]范围的网格内的长度均为D；其他网格内的长度均为0。
7.如权利要求1所述的方法，其特征在于，所述步骤S6中计算k 1 ≠ 0的射线在各网格内的长度的方法为：
S61 、计算射线与网格n 1 左边界所在直线的交点的纵坐标y 1 , 根据公式：计算左截距d 1 ；
S62 、计算射线与网格n 1 右边界所在直线的交点的纵坐标y 2 , 根据公式：计算右截距d 2 ；
S63 、根据如下公式计算射线在网格n 1 内的长度：
若k 1 > 0，
若k 1 < 0，
S64 、根据如下公式判断射线经过的下一个网格的编号c：
若k 1 > 0，
若k 1 <0，
S65 、重复以上过程，直到计算出射线在网格n 2 内的长度为止；
S66 、记射线在其他网格内的长度为0。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201410265945.6A CN104020508A (zh) | 2014-06-16 | 2014-06-16 | 一种用于地质雷达层析探测的直射线追踪算法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201410265945.6A CN104020508A (zh) | 2014-06-16 | 2014-06-16 | 一种用于地质雷达层析探测的直射线追踪算法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN104020508A true CN104020508A (zh) | 2014-09-03 |
Family
ID=51437363
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201410265945.6A Pending CN104020508A (zh) | 2014-06-16 | 2014-06-16 | 一种用于地质雷达层析探测的直射线追踪算法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN104020508A (zh) |
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109102553A (zh) * | 2018-06-27 | 2018-12-28 | 中国人民解放军战略支援部队航天工程大学 | 二维重构算法中极坐标系统矩阵计算方法和装置 |
| CN110045369A (zh) * | 2019-05-07 | 2019-07-23 | 中国矿业大学(北京) | 一种探地雷达层析探测曲线追踪方法 |
| CN110929375A (zh) * | 2019-10-17 | 2020-03-27 | 中国科学院电子学研究所 | 基于二维矩量法和射线追迹法的透镜高效仿真、优化方法 |
| CN112068128A (zh) * | 2020-09-19 | 2020-12-11 | 重庆大学 | 一种直道场景线段型雷达数据处理及位姿获取方法 |
| CN113156495A (zh) * | 2020-01-07 | 2021-07-23 | 中国石油天然气集团有限公司 | 网格层析反演反射点确定方法及装置 |
-
2014
- 2014-06-16 CN CN201410265945.6A patent/CN104020508A/zh active Pending
Cited By (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109102553A (zh) * | 2018-06-27 | 2018-12-28 | 中国人民解放军战略支援部队航天工程大学 | 二维重构算法中极坐标系统矩阵计算方法和装置 |
| CN110045369A (zh) * | 2019-05-07 | 2019-07-23 | 中国矿业大学(北京) | 一种探地雷达层析探测曲线追踪方法 |
| CN110045369B (zh) * | 2019-05-07 | 2021-04-27 | 中国矿业大学(北京) | 一种探地雷达层析探测曲线追踪方法 |
| CN110929375A (zh) * | 2019-10-17 | 2020-03-27 | 中国科学院电子学研究所 | 基于二维矩量法和射线追迹法的透镜高效仿真、优化方法 |
| CN113156495A (zh) * | 2020-01-07 | 2021-07-23 | 中国石油天然气集团有限公司 | 网格层析反演反射点确定方法及装置 |
| CN112068128A (zh) * | 2020-09-19 | 2020-12-11 | 重庆大学 | 一种直道场景线段型雷达数据处理及位姿获取方法 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Feng et al. | Sectional velocity model for microseismic source location in tunnels | |
| US10795053B2 (en) | Systems and methods of multi-scale meshing for geologic time modeling | |
| CN104020508A (zh) | 一种用于地质雷达层析探测的直射线追踪算法 | |
| CN103775071B (zh) | 采动煤岩体裂隙演化的测量方法 | |
| CN103245977B (zh) | 一种矿井回采区灾害源的地质雷达层析探测方法 | |
| EA201070615A8 (ru) | Формирование геологической модели | |
| CN105093319B (zh) | 基于三维地震数据的地面微地震静校正方法 | |
| CN102057368B (zh) | 利用最大连续场在三维体积模型中分布性质 | |
| CN105334542A (zh) | 任意密度分布复杂地质体重力场快速、高精度正演方法 | |
| CN105549077B (zh) | 基于多级多尺度网格相似性系数计算的微震震源定位方法 | |
| CN105005081B (zh) | 煤机采动激励下综采面近场煤岩动态层析成像系统及方法 | |
| CN108414983B (zh) | 一种基于逆时射线追踪方法的微地震定位技术 | |
| CN103529486B (zh) | 一种地球化学异常圈定方法 | |
| CN108845358B (zh) | 断层及构造异常体识别方法及装置 | |
| CN107132578A (zh) | 一种微地震地面监测速度模型校正算法 | |
| CN104360396B (zh) | 一种海上井间tti介质三种初至波走时层析成像方法 | |
| CN105092023B (zh) | 基于白噪声统计特征的锚杆振动信号矫正方法 | |
| CN110516650A (zh) | 一种基于震动传感器阵列的浅层盲空间震源定位系统 | |
| CN102692644A (zh) | 生成深度域成像道集的方法 | |
| CN103389489A (zh) | 基于大斜度井的微地震监测定位方法 | |
| CN103778298A (zh) | 改进的模拟多孔介质中二维水流运动的多尺度有限元方法 | |
| CN107045141B (zh) | 基于反时到时差数据库的微地震/地震震源快速定位方法 | |
| CN102877828A (zh) | 一种三维多井联合井地ct成像方法 | |
| CN104345336A (zh) | 一种基于目标区域照明度的观测系统优化方法 | |
| CN104570070B (zh) | 一种建立二维近地表地质模型的方法及设备 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| WD01 | Invention patent application deemed withdrawn after publication | ||
| WD01 | Invention patent application deemed withdrawn after publication |
Application publication date: 20140903 |