基于DEM数据体的山地地形地类划分方法及在物探经济技术评价中的应用

2023-02-14杨连刚江民潘树林罗新生宋金鹏罗姜民

石油地球物理勘探 2023年1期

杨连刚，江民，潘树林，罗新生，宋金鹏，罗姜民

(1.中国石油塔里木油田公司，新疆库尔勒 841000； 2.西南石油大学地球科学与技术学院，四川成都 610500)

0 引言

塔里木盆地库车和塔西南山前受构造挤压、沉积风化等因素影响，地表发育山体、沟壑等多种地形，施工条件复杂[1-5]。在地震采集项目部署、技术方案论证、工程造价测算等环节均需对山地地形地类进行详细划分。传统方法通常采用多轮次联合野外踏勘方式，在踏勘效率、成本、安全性与全面性等方面均存在较大不足，主要体现在：①室内地形地类划分以遥感信息和人工目视解译为主，多为定性评价，对地形、地貌信息的定量表征较少；②野外踏勘能够到达的实际区域受限，存在“以偏概全、以点概面”的问题；③室内地形地类划分需收集大量资料，同时需要人工识别和统计，耗时较长；④野外踏勘需要大量人力和物力，同时难以进入尚未修路的工区，短期内难以实现精细分区。

近年来，随着中秋1、博孜9等油气藏的相继发现，复杂山前构造带勘探力度持续加大，地震采集呈现出工作量大、周期短、节奏快、复杂程度日趋增高等特点[6-9]，仅仅通过野外踏勘划分地形地类已无法满足油田生产快节奏、科学化的发展需求，亟需探索一种科学合理、客观、经济高效、可重复性强的地形地类划分方法，以便准确地掌握复杂山地工区的炮、检点地形地类构成，从而更好地满足地震采集项目成本测算的需求。

高精度遥感信息广泛运用于国土资源调查、区域地貌研究、地质灾害预防、地震采集方案优选等方面。封志明等[10]利用遥感信息研究地形起伏度与人口分布的相关性；朗玲玲等[11]提出多尺度数字高程模型(Digital Elevation Model，DEM)提取起伏度划分低山丘陵区地形、地貌；牛全富等[12]从坡度、起伏度等地形特征要素中优选敏感因子实现了黄土区滑坡灾害分级预测；叶勇等[13]为改善高、陡山地极低信噪比地区测线的采集条件和可实施性，通过人工加权高程、坡度、相对湿度、岩性等信息计算出评价因子优选地震测线部署方案；陈学强等[14]、张友焱等[15]提出利用高程、坡度等逐点优选检波点改善地震资料品质；赵邦六等[16]、陈德武等[17]介绍了遥感技术在地震采集方面发挥的重要作用，并指出遥感信息在地震勘探中的应用潜力需进一步挖掘。山地地形地类划分通常根据海拔高程并结合坡度、起伏度，由于地区特性和研究需要不同，划分方式和标准也大不相同。前人的研究由于计算地理尺度较大，同时地震采集施工难度表征具有一定特殊性，无法直接应用，但为本次研究提供了“利用DEM数据体来划分山地地形地类”的技术思路。

结合多年生产实践，本文提出一种基于高精度DEM数据体炮、检点地形地类量化评价方法，通过提取遥感信息中的坡度、起伏度、纵深距离等地形特征要素定量表征野外施工难度，并在当前塔里木盆地地震采集工程造价水平下，依据不同要素特点，合理建立地类划分区间，实现地形地类的准确划分。

1 方法原理

目前影响塔里木盆地山地地震采集项目施工难度的主要因素有山地高大程度、山体风化破碎程度、山地纵深距离、山地地表岩性等。其中地表岩性主要影响钻井施工，而该部分在工程造价中已有明确的补偿系数，因此在地形地类中暂不考虑。在近年山地地形地类划分实践过程中，将地形特征要素中的起伏度、坡度和纵深距离作为划分地貌类型的3个重要指标。

1.1 地形特征要素基本概念

1.1.1 起伏度

起伏度是指在当前位置一定空间范围内，最高点与最低点海拔高度的差值[18]，可作为地震地形地类划分过程中描述山体区域地形特征的宏观指标，主要反映山体高大程度(图1)。它表征地震施工纵向上采集设备搬迁等施工难度。地形起伏度R的计算公式为

图1 地形起伏度计算示意图

R=Hmax-Hmin

(1)

式中Hmax、Hmin分别表示计算范围内DEM高程栅格数据体的最大值、最小值。

1.1.2 坡度

图2 高程网格(左)及地形坡度计算示意图(右)

(2)

利用三阶反距离平方权差分代替微分，则离散形式为

(3)

(4)

式中：h11、h12、h13、…、h32、h33为对应高程栅格数据的值；dx为栅格间距，即DEM数据体的精度；Hx、Hy分别为x、y方向差分算子。

1.1.3 纵深距离

为完善起伏度、坡度两个地形特征要素在表征山地纵、横向上施工搬迁难度的不足，引入纵深距离的概念，其为从当前位置沿山体边界法线方向向内延伸的距离，它也是描述山体区域地形特征的一个宏观指标，在起伏度和坡度划分结果的基础上调整地类，从而更准确的表征野外施工难度。

为计算方便，利用角平分线法将山体按照不同距离建立若干个与山体边界等距离区域，即多个山体边界纵深距离分析区。纵深距离分析区内任意位置沿其法线方向距离山体边界相等。

图3 山体边界纵深距离分析示意图

(6)

(7)

(8)

(9)

重复计算式(6)～式(9)即可求得纵深距离为L的分析区边界点坐标。

1.2 邻域大小计算和地类划分区间确定方法

上述起伏度、坡度、纵深距离等地形特征要素能否准确表征山体施工难度，关键在于3个方面：DEM数据体的分辨率大小；地形特征要素中起伏度、坡度的最佳邻域大小；各地形特征要素的地类划分区间。

由于地形地类划分主要应用于地震采集技术经济技术一体化评价，侧重于地震采集工程造价，不需采集方案设计中炮、检点优化选点布设那么高要求如(大于或等于5 m的DEM数据体分辨率)，同时结合数据体获取难度和近年生产实践，认为DEM数据体分辨率大于或等于10 m即可。

为适应塔里木盆地山地复杂多样，这里分别采用数理统计中的两种方法：一是均值变点法，用于求取地形特征要素的计算邻域大小；二是自然间断分级法，用于求取地形特征要素的地类划分区间。

1.2.1 均值变点法

(10)

(11)

Ti=S-Si

(12)

当Ti取最大值时对应位置为变点位置，即最佳的地形特征要素计算邻域。

1.2.2 自然间断分级法

自然间断法是地理信息中的一种分类区间确定方法[20]，该方法能依据数据的分布特征，将数据按照区间内相似性最高、区间之间差异性最大划分为多个区间。

数据序列平均值的偏差平方和SDAM记作

(13)

SDCM1=SDCMi1+SDCMi2+…+SDCMik

(14)

最后利用方差拟合优度GVF验证

(15)

GVF的取值范围为0～1，其值越大表示组间差异越大。在地形地类划分中，一般认为GVF≥0.7时，自然间断法计算的分类结果较好。

1.3 综合划分地形地类

依照中国石油天然气集团公司现行定额标准，塔里木盆地地形地类有山地、戈壁、平原农田、沼泽水网和沙漠共5种地表类型。山地地震项目一般包括山地地表类型与戈壁、农田、沼泽水网3种非山地地表类型。山地地形划分标准如表1所示，当前划分标准多为定性描述，受主观因素影响大。

表1 塔里木盆地山地地形地类划分现行定额标准

为使山地地类划分工作在现行定额下相对客观合理、划分过程利于量化统一，提出了一种划分方法(图4)，利用10 m精度DEM数据体计算坡度、起伏度、纵深距离，并依据施工难度将山地地类由小到大划分为Ⅰ、Ⅱ、Ⅲ、Ⅳ、Ⅴ及Ⅵ类；非山地由于区域连续且范围较小，可直接通过高清卫片结合现场踏勘情况划分地表类型及难度类别。

图4 塔里木油田地震项目地形地类划分示意图

首先利用起伏度、坡度划分地形地类，将任意DEM数据体栅格位置由起伏度地类划分结果记作IR(x,y)，坡度地类划分结果记作Iθ(x,y)，取值为0、1、2、3、4、5、6，分别对应非山地、Ⅰ类山地、Ⅱ类山地、Ⅲ类山地、Ⅳ类山地、Ⅴ类山地、Ⅵ类山地等7种地类，计算起伏度、坡度两者划分结果最大值，得到的地类记作YZ0(x,y)，即

YZ0(x,y)=max[IR(x,y),Iθ(x,y)]

(16)

式中YZ0(x,y)表征山地施工难度。实际生产表明，Ⅲ类及Ⅲ类以上山体包围的区域(如山间相对平坦的区域或更加陡峭的区域)施工难度会较该区域直接计算结果大，为进一步刻画此类区域搬迁距离对施工难度的影响，利用纵深距离对YZ0(x,y)做一定调整。沿Ⅲ类山体边界建立3个不同纵深距离分析区，不同纵深距离地类调整因子记作D0(x,y)，即

(17)

(18)

图5为山地地形地类自动划分的计算流程，在利用上述方法计算部署地震采集项目综合地形地类之后，根据炮点、检波点、微测井点坐标提取、统计每个点对应地类，结合地表岩性即可快速得到地震采集工程计价清单，从而实现高效测算项目工程造价，在采集设计阶段完成经济技术一体化评价，从而优选采集设计方案。

图5 塔里木油田山地地形地类自动划分计算流程

2 地形地类划分关键参数

2.1 坡度、起伏度的计算邻域的确定

为确定地形地类划分中坡度、起伏度的最佳计算邻域，选取秋里塔格构造带典型山地地表西秋—却勒—佳木段中的佳木地区作为分析数据(图6a蓝色区域)，并利用前文介绍的均值变点法分别计算(采用航拍的1 m精度DEM数据体)100 m×100 m、200 m×200 m、…、1000 m×1000 m范围下的起伏度平均值和3.0 m×3.0 m、7.5 m×7.5 m、15.0 m×15.0 m、…、90.0×90.0 m范围下坡度平均值作为均值变点分析序列。

如图6所示，计算结果表明起伏度平均值随邻域大小增加而增加，变点位置对应邻域大小为500 m，坡度平均值随邻域大小增加而减小，变点位置计算邻域大小为30 m，即起伏度、坡度对应的最佳计算邻域大小分别为500 m和30 m。

图6 坡度、起伏度对应的计算邻域

根据近年塔里木油田地震采集现场地形地类划分经验，将塔里木盆地地震采集地形地类划分中山地定义为： 500 m计算范围内起伏度大于30 m或30 m计算范围内坡度大于5°(以DEM数据体分辨率≥10 m计算)的地表或第四系西域组以下老地层出露区。

依据上述定义计算出的西秋—却勒—佳木段山地区域(图6f黑色阴影部分)与DEM数据体吻合较好，同时秋里塔格构造带和克拉苏构造带划分出的山地分布区域与现场实际地形高度吻合，表明此类山地与非山地划分方法可行、有效。

2.2 起伏度、坡度地类划分区间确定

为得到适用于整个塔里木盆地山地地类坡度和起伏度的合理划分区间，选取近三年新采集且具有代表性的6个地震采集项目用于分析，包括库车山前的佳木3、却勒4、克深1三维，柯坪地区的柯坪4、沙井子二维以及塔西南山前的乌泊尔二维。将山地地形地类划分为非山地和山地共7种类型，利用自然间断法分别计算坡度和起伏度计算划分区间，并计算方差拟合优度，计算结果见表2、表3。计算结果表明，采用自然间断法分类计算出的各项目起伏度、坡度方差拟合优度为0.96～0.98，表明该方法在山地地类划分中是有效的。各项目之间的划分结果虽然存在一定差异，但整体比较接近。由于二维项目工区跨度大于三维项目，其区间划分范围较三维项目稍大一些。

为统一地类划分的标准，同时与当前塔里木油田地震采集工程造价水平和现有地类划分认识相适应，需对坡度和起伏度地类划分区间做出一定调整，调整后分类方案如表2、表3所示。从表可见，调整后所有项目方差拟合优度均大于0.91，表明调整后地类划分区间的聚类效果仍然合理。

表2 自然间断法起伏度分析结果

表3 自然间断法坡度分析结果

3 应用实例——以库北1三维为例

库北1三维构造位置位于库车坳陷北部构造带库北单斜带，南部紧邻克拉苏构造带，位于阿克苏地区拜城县境内，部署工作量偏移处理前的满覆盖面积为194.72 km2，其中东西长19.8 km、南北宽9.4 km。

工区地形地貌以山地地表(图7)为主，中部主体为一般山地，北部为落差较大的高难山体，非山体地表为戈壁。高密度采集中炮点范围海拔高程为1550～2750 m，检波点范围海拔高程为1450～2900 m，属于库车山前典型的复杂山地地表类型。以往利用野外踏勘与卫片结合的方式划分地形地类，难以高效准确地表征项目施工难度。

图7 库北1三维工区地理信息

工区选取5处典型地表位置(图8a中的A～E处)现场验证标定，并利用起伏度和坡度划分地形地类(图8b)，仔细对比可以看出，划分结果与实际地表情况基本相符，但工区北部和中部Ⅲ类及Ⅲ类以上包围的高难山体区域局部地类偏低(图8b中的C和E处)，沿Ⅲ类山体边界建立不同纵深距离分析区(图8c)来调整地类，并依据纵深距离大小对地形地类进行调整，调整后的地类(图8d)与实际地表更加吻合，划分结果与现场标定结果一致。

图8 库北1三维工区地形地类划分

在地震采集部署设计阶段，一般会结合地质任务和工区特点提出2或3种采集方案，库北1三维工区共有2种设计了2种采集方案(表4)。

分别根据2种方案炮点、检波点坐标位置提取对应的地类，并统计得到炮点点、接收点不同地类占比(图9)。两种方案不同地表类型炮点、接收点占比较为接近，非山地约占8%，Ⅰ类～Ⅲ类山体约占72%，Ⅳ类～Ⅵ类山体占20%。

图9 库北1三维工区激发点、接收点不同地类占比柱状图

在观测系统参数基础上结合激发参数、接收参数和近地表岩性及厚度情况，即可快速获得量化后的用于项目工程造价测算清单，为经济技术一体化评价优选出最佳采集方案提供决策依据。

依据塔里木油田地震项目工程计价标准分别计算表4和表5测算清单工程造价，方案1较方案2在炮道密度提升36.64万道/km2的同时，每平方公里地震采集项目施工费用单价降低0.99%，节约成本约115.94万元。综合考虑成本和采集效果，方案1明显优于方案2，这种“以道换炮”提高炮道密度、降低勘探成本的策略在塔里木探区山地地震采集经常用到。

表4 库北1三维观测系统设计方案

表5 库北1三维项目采集方案1测算清单

目前基于起伏度、坡度、纵深距离三地形特征要素的复杂山地炮、检点地形地类高效划分方法已经应用到塔里木油田近3年全部山地地震采集项目经营测算和经济技术一体化评价中，并以此研发了基于GIS平台的复杂山地地震智能工程造价系统。该方法大大提高了地形地类划分的客观性和科学性，野外踏勘周期从1周以上缩短至1天，同时踏勘范围也从以往全工区踏勘转变成安全有效的局部针对性踏勘，整个采集项目地类划分周期由以往的1周以上缩减至2天以内。