基于ArcGIS与MATLAB的资源量计算方法——以松辽盆地油页岩为例

张静平, 唐书恒

(1.中国石油勘探开发研究院廊坊分院,河北廊坊065007;2.中国地质大学(北京)能源学院,北京100083)

基于ArcGIS与MATLAB的资源量计算方法
——以松辽盆地油页岩为例

张静平1, 唐书恒2

(1.中国石油勘探开发研究院廊坊分院,河北廊坊065007;2.中国地质大学(北京)能源学院,北京100083)

基于钻孔资料,采用体积法对松辽盆地青一段、嫩一段、嫩二段油页岩进行了资源量计算,依据分布面积最大的嫩二段建立基准坐标体系,确定了计算范围,采用地理信息系统ArcGIS对3个层段油页岩面积进行了子块划分与子块面积计算,青一段划分为147个子块,嫩一段划分为197个子块,嫩二段划分为246个子块,子块单位实际面积为537 km2。采用趋势面法,在嫩二段坐标体系下以钻孔地理坐标为自变量,其对应的厚度趋势值为因变量,采用算法软件MATLAB求得3层油页岩段厚度趋势值方程,并对油页岩厚度空间展布进行拟合,获得了较为理想的拟合度。在确定子块面积与厚度的基础上,对每个子块加权求和分别得出青一段资源量为12 724亿t、嫩一段为8 873.1亿t、嫩二段为9 872.4亿t。

体积法; 趋势面法; ArcGIS; MATLAB; 油页岩; 松辽盆地

松辽盆地为大型坳陷盆地,油页岩呈层状大面积平缓分布,品位和厚度变化都比较稳定,多位学者采用体积法进行了油页岩资源量的估算,但对计算中涉及的指标取值方法并不统一。

张海龙[1]对松辽盆地三套油页岩进行了资源量计算,预测区边界由岩矿厚度等值线图圈定,厚度和面积由油页岩厚度等值线图确定,计算时以两条厚度等值线圈定一个块段(内插0.7 m厚度等值线),再将各块段资源累加。王慎余等[2]进行油页岩预测时,含矿面积主要采用湖相(湖沼相)沉积面积,含矿层厚度直接采用或间接引出。2006年由国土资源部、财政部、国家发改委所启动的新一轮全国范围的油页岩资源评价,主要采用体积法进行,面积采用勘查区面积,厚度按油页岩厚度等值线划分若干块段进行,若无厚度等值线,利用勘查区油页岩厚度用资源储量加权[3]。

在采用体积法对松辽盆地油页岩资源量进行预测工作中,面积的确定有的采用湖相面积,有的采用盆地边界,有的采用勘查区边界,并不统一。厚度的确定大多以厚度等值线来的圈定,由于厚度等值线由差值生成,因此计算中厚度的准确性值得商榷。

本文探讨了体积法运用中面积与厚度的确定方法,旨在提供一种较为合理的固体矿产资源量计算的方法,为其他地区的相关工作提供参考。

1 油页岩资源量计算方法

1.1 体积法概述

油页岩的资源量与其体积成正比,因此采用体积法估算资源量是可行的。体积法的一般公式:

式中,M为资源量,ˉC为某种矿产的平均含量,V为预测对象的体积。

将油页岩矿的矿石量视为M,即认为油页岩矿石量为油页岩资源量,对矿层划分成子块进行体积的计算,所以将公式(1)改写成如下形式:

式中,M为油页岩矿资源量(矿石量),D为油页岩矿平均体重,Si为子块(见下文定义)面积,Hi为子块内可采油页岩的平均厚度,i为子块顺序号,SiHi是某子块油页岩矿的体积,SiHi是预测区油页岩矿总体积。D为常数,沿用以往D=2t/m3。在这个计算过程中,确定了子块的面积与平均厚度,便可求得任一子块的体积,进而求和获得各层油页岩的资源量[4-6]。

1.2 面积计算

使用Arc GIS软件完成面积计算、坐标校正工作,它是一个完整的地理信息集成系统,具有独立模块功能,提供数据管理,分析统计等功能,应用领域广泛[7-9]。

1.2.1 确定资源量计算区域 以分布最广的嫩二段为基准建立坐标体系,以其平面分布图最西边做切线建立Y轴,其大致为垂直走向,以平面分布最南边做切线建立X轴,其大致为水平方向(见图1)。计算区域面积设定为185 504 km2,青一段、嫩一段完全涵盖于此计算区域内。

图1 松辽盆地油页岩资源量计算区域图Fig.1 The caculation area of reserves about oil shale in songliao basin

1.2.2 油页岩子块面积划分方案与计算 计算区域确定后,以相等间隔在ArcGIS软件里生成所需格网,然后通过手工方式删除油页岩分布范围外的格网(见图2)。

图2 ArcGIS面积子块划分方案Fig.2 The procedure about subblock area division in ArcGIS

每个子块单元实际边长为23 km,面积为527 km2,但子块在边界上或与剥蚀区相交的区域由于不规则均小于527 km2。子块划分过细会加大工作量,对于平面分布平缓的油页岩层也没有必要。本次计算嫩二段划分为246子块1 195个节点进行计算,青一段划分147子块1 007个节点,嫩一段为197子块975个节点。

1.3 厚度的计算方法

面积确定后,计算资源量的最关键因素就是厚度的确定。每一个正方形子块的厚度为4个节点厚度平均值,但边界上不规则子块的厚度由4～8个节点的平均值来控制,以减小误差。厚度的确定采用多项式趋势面分析法,并借助MATLAB算法软件获得厚度的趋势值方程,可方便地计算任意节点的厚度值,进而求取任一子块的平均厚度进行体积计算。

厚度确定原理:多元统计的趋势面分析已经被广泛应用于地质、地理等多学科领域,在趋势值与剩余值的确定中被证实是可行的有效的方法。其原理是以一个多项式曲面来模拟空间点的分布规律[10-11](见图3)。

假设现有n个观测点,为得到观测点空间分布的规律,假定空间存在一曲面来拟合原始值的变化规律,曲面方程为多项式,表达式为:^zi=a0+a1xi+a2yi+…。

图3 多项式趋势面分析法原理图Fig.3 The principle map of the polynomial trend surface analysis

令A=XTX,B=(a1a2…aL)T,C=XT·Z, Z=(Z1Z2…Zn)T,则A·B=C,其中L为该多项式的系数个数,其由多项式次数决定,运行MATLAB软件,解此正规方程组,得出多项式系数,即可确定厚度趋势面的方程。

趋势面结果可经过拟合度检验,偏差平方和与回归平方和的比例一般认为达到40%为佳,拟合指数公式为:

2 松辽盆地油页岩资源量计算

每一个子块资源量求和便得到该层油页岩矿的资源量,每一个子块的实际面积在上文面积划分时已经求得,子块的厚度计算较为复杂,依靠趋势面法拟合的多项式曲面方程来计算。

2.1 油页岩厚度计算结果

(1)嫩二段厚度计算结果

嫩二段分别进行了3次与4次多项式的计算对比。多项式次数以满足计算要求为原则。厚度3次多项式趋势面求解过程共遴选了68口钻孔来进行,钻孔分布较为均匀,钻孔坐标在ArcGIS环境下可方便的由大地坐标校正到本文设置的嫩二段直角坐标系下,转化后坐标数值较小,方便MATLAB矩阵计算。根据趋势面原理,3次趋势面正规方程组为68行10列矩阵,68为观测点个数,本文为钻孔数量,L为10,即3次多项式的项数。

在MATLAB命令窗口读入该矩阵,获得3次趋势面多项式如下:

Z=-1.2 8 4 0-0.4 0 4 0 x+4.3 8 0 2y+ 1.085 9x2-1.126 9xy-0.559 9y2-0.223 9x3+ 0.268 1x2y-0.071 5xy2+0.046 6y3。

基于此方程可计算出任意节点的厚度趋势值,并进一步获得子块的平均厚度值。调入以下命令生成该多项式在空间的曲面,即厚度的趋势面[14-15]。

[x,y]=meshgrid(0:0.25:15.5,0:0.25:21.8); Z=-1.284 0-0.404 0*x+4.380 2*y+ 1.085 9*x.^2-1.126 9*x.*y-0.559 9*y.^2-0.223 9*x.^3+0.268 1*x.^2.*y-0.071 5*x.*y.^2+0.046 6*y.^3; surface(x,y,z)

依此方法对嫩二段进行4次趋势面分析,对比后,3次多项式趋势面更能反映实际厚度变化趋势(见图4,图5)。

图4 嫩二段3次趋势面厚度图Fig.4 Three power polynomial about thickness trend of the 2nd member of Nenjiang group

图5 嫩二段4次趋势面厚度图Fig.5 Four power polynomial about thickness trend of the 2nd member of Nenjiang group

经检验,3次趋势面实际拟合度为R2=52%,较为理想。趋势值较好的反应了实际厚度的变化趋势(见表1),其高值与低值趋势均与实际对应。

表1 嫩二段厚度趋势值与实际值对比表Table 1 The comparison of predicted and real thickness of the 2nd member in Nenjiang group

续表1

(2)青一段厚度计算结果

青一段分别进行了2、3次趋势面多项式拟合,共计54口钻孔参与趋势面分析。

优选3次趋势面方程为:

Z=-44.749 0+9.498 8x+28.122 7y-3.284 7x2+0.311 1xy-4.736 5y2+0.336 7x3-0.047 1x2y-0.036 7xy2+0.241 0y3

由于青山口厚度变化较大,造成数据点起伏较大,使得趋势面在局部点拟合度不高,对比具体钻孔数据,其趋势值整体上能够反应厚度平均值的变化趋势(见表2),可以满足子块厚度计算要求。

表2 青一段厚度趋势值与实际值对应表Table 2 The comparison of predicted and real thickness of the 1st member in Qingshankou group

(3)嫩一段厚度计算结果

同理,嫩一段进行了3次、4次多项式趋势面比较,4次趋势面不能反映局部厚度的变化,且出现大范围畸变,最终确定3次趋势面方程Z=38.826 0-23.853 7x-7.095 6y+6.167 1x2+1.784 4xy+ 0.600 6y2-0.498 6x3-0.332 3x2y+0.134 4xy2-0.051 0y3,其拟合度达到35%。

2.2 松辽盆地资源量计算结果

以嫩二段资源量计算为例,计算面积一共划分为246个子块,将子块进行编号,各个子块的面积可以从ArcGIS软件中读取。子块的平均厚度由对应的趋势面多项式方程求得。子块的面积和厚度确定后,可以求取子块的体积与资源量,将所有246个子块资源量求和便获得嫩二段油页岩的资源量(见表3)。

经过同样的计算过程,可计算青一段、嫩一段的资源量。嫩二段246个子块资源量为9 872.4亿t、青一段147个子块12 724亿t、嫩一段197个子块8 873.1亿t。

3 结论

本文提供了一种体积法计算中面积与厚度的确定方法。子块面积采用ArcGIS软件进行计算;子块厚度计算依据趋势面法的原理,以钻孔坐标为自变量,钻孔所揭示的油页岩厚度为因变量,采用MATLAB软件进行趋势面多项式方程拟合,利用方程计算厚度值。

松辽盆地资源总量为31 469.5亿t,其中青山口组一段油页岩资源量规模最大,12 724亿t,占盆地整体40.43%,嫩江组二段次之,9 872.4亿t,占盆地整体资源量3 1.3 7%,嫩江组一段资源量为8 873.1亿t,占盆地整体资源量28.20%。

表3 嫩二段油页岩资源量计算参数表Table 3 The reserve calculation parameters of the 2nd member in Nenjiang group

[1] 张海龙.东北北部区油页岩资源评价及评价方法研究[D].长春:吉林大学,2008.

[2] 王慎余,许家朋,王振海.我国油页岩资源开发利用状况及发展对策[J].中国地质经济,1992,5:16-19.

Wang Shenyu,Xu Jiaming,Wang Zhenhai.Geology and mineral resources bureau of Guangdong province[J].Geological Economy of China,1992,5:16-19.

[3] 刘招君,柳蓉.中国油页岩特征及开发利用前景[J].地学前缘,2005,12(3):315-323.

Liu Zhaojun,Liu Rong.The characteristics and use of oil shale in China[J].Earth Science Frontiers,2005,12(3):315-323.

[4] 孟庆涛.油页岩资源评价方法研究—以松辽盆地南部农安地区为例[D].长春:吉林大学,2007.

[5] 董清水,王立贤,于文斌,等.油页岩资源评价关键参数及其求取方法[J].吉林大学学报(地球科学版),2006,36(6):898-903.

Dong Qingshui,Wang Lixian,Yu Wenbin,et al.The key parameters of oil shale resource appraisement and its evaluating methods[J].Journal of Jilin University(Earth Science Edition),2006,36(6):898-903.

[6] 汪云甲,黄宗文.矿产资源评价及其应用研究[M].北京:中国矿业大学出版社,1998.

[7] 范建伟.基于ArcGIS评价油页岩资源评价—以Piceance盆地绿河组油页岩为例[D].北京:中国地质大学(北京),2010.

[8] 黄波,华韵子.ArcGIS在气象信息系统中的应用[J].大气科学研究与应用,2010,2:54-59.

Huang Bo,Hua Yunzi.Application of ArcGIS to the meteorological information system[J].Atmospheric Science Research and Application,2010,2:54-59.

[9] 高鹏,徐志刚,颜远钦.基于ArcGIS的县级土地利用现状调查研究[J].北京测绘,2011(2):17-20.

Gao Peng,Xu Zhigang,Yan Yuanqin.Research of county-level land-use investigation based on ArcGIS[J].Beijing Surveying and Mapping,2011(2):17-20.

[10] 李小明,刘德民,张品刚,等.趋势面分析在东欢坨井田构造研究中的应用[J].中国煤炭地质,2010,22(4):6-9.

Li Xiaoming,Liu Demin,Zhang Pingang,et al.Application of trend surface analysis in Donghuantuo mine field structural study[J].Coal Geology of China,2010,22(4):6-9.

[11] 孟祥霞,姚清洲,刘建新,等.趋势面分析技术在塔北英买2井区断裂识别中的应用[J].新疆石油地质,2010,31(5):543-547.

Meng Xiangxia,Yao Qingzhou,Liu Jianxin,et al.Application of trend-surface analysis to fault recognition in Yingmai-2 well area in Tabei area,Tarim basin[J].Xinjiang Petroleum Geology,2010,31(5):543-547.

[12] 路鹏,周荣福,王玫.基于趋势面分析的煤层空间分布特征研究[J].中州煤炭,2009,5:9-11.

Lu Peng,Zhou Rongfu,Wang Mei.Research on coal seam s spatialdistribution features based on trend surface analysis [J].Zhongzhou Coal,2009,5:9-11.

[13] 徐建华.计量地理学[M].北京:高等教育出版社,2006:28-31.

[14] 张瑞林,李东印,李小军.趋势面拟合法在煤矿瓦斯地质变量预测中的应用研究[J].中国安全科学学报,2006,16(12): 25-29.

Zhang Ruilin,Li Dongyin,Li Xiaojun.Research on the application of trend surface fitting method to the forecast of gas and geology variables in coal mine[J].China Safety Science Journal,2006,16(12):25-29.

[15] 李哲,张淑英.基于多项式趋势面分析理论的天然气需求预测[J].资源与产业,2008,10(2):105-107.

Li Zhe,Zhang Shuying.Prediction of natural gas demands based on multi-trend analysis[J].Resources&Industries, 2008,10(2):105-107.

(编辑 王亚新)

Evaluating Method about Mineral Resource Based on ArcGIS and MATLAB: Take Oil Shale of Songliao Basin as an Example

Zhang Jingping1,Tang Shuheng2
(1.Research Institute of Petroleum Ex ploration and Development-Langfang,Langfang Hebei 065007,China; 2.School of Energy Resources,China University of Geosciences,Beijing 100083,China)

Based on bore hole data,the volumetric method to calculate the resources of oil shale in the 1st member of Qingshankou formation and the 1st-2nd member of Nenjiang formation in Songliao Basin were adopted.A reasonable area dividing scheme was established taking the 2nd member of Nenjiang formation as criterion when the area distribution of the three sets of oil shale was determined.Sub-block dividing and sub-block area calculation were proceeded on the three researching formations using geographic information system ArcGIS and the actual area of each sub-block was 537 km2.1st member of Qingshankou formation,1st member of Nenjiang formation and 2nd member of Nenjiang formation were divided into 147 sub-blocks,197 sub-blocks and 246 sub-blocks respectively.According to the theory of trend surface method,oil shale thickness tendency value equation was obtained using algorithm software MATLAB where the spatial distribution coordinates were considered as independent variable and the relative thickness tendency values as dependent variable.Further more,this paper did fitting analysis on oil shale thickness spatial distribution,and the fitting results were satisfied which ensured the accuracy of calculation.The oil shale resources of the three formations were obtained after the determination of the area and thickness of those sub-blocks.The results are 1 272 400 million tons,887 310 million tons and 987 240 million tons for 1st member of Qingshankou formation,1st member of Nenjiang formation and 2nd member of Nenjiang formation respectively.

Volumetric method;Trend surface method;ArcGIS;MATLAB;Oil shale;Songliao basin

TE155

A

10.3969/j.issn.1006-396X.2014.03.012

1006-396X(2014)03-0052-05

2012-12-22

2013-06-25

长江学者和创新团队发展计划(IRT0864)。

张静平(1981-),女,博士,工程师,从事非常规油气地质研究;E-mail:zhangjp69@petrochina.com.cn。