基于ArcGIS與MATLAB的資源量計(jì)算方法
——以松遼盆地油頁(yè)巖為例

張靜平, 唐書恒

(1.中國(guó)石油勘探開發(fā)研究院廊坊分院,河北廊坊065007;2.中國(guó)地質(zhì)大學(xué)(北京)能源學(xué)院,北京100083)

基于ArcGIS與MATLAB的資源量計(jì)算方法
——以松遼盆地油頁(yè)巖為例

張靜平1, 唐書恒2

(1.中國(guó)石油勘探開發(fā)研究院廊坊分院,河北廊坊065007;2.中國(guó)地質(zhì)大學(xué)(北京)能源學(xué)院,北京100083)

基于鉆孔資料,采用體積法對(duì)松遼盆地青一段、嫩一段、嫩二段油頁(yè)巖進(jìn)行了資源量計(jì)算,依據(jù)分布面積最大的嫩二段建立基準(zhǔn)坐標(biāo)體系,確定了計(jì)算范圍,采用地理信息系統(tǒng)ArcGIS對(duì)3個(gè)層段油頁(yè)巖面積進(jìn)行了子塊劃分與子塊面積計(jì)算,青一段劃分為147個(gè)子塊,嫩一段劃分為197個(gè)子塊,嫩二段劃分為246個(gè)子塊,子塊單位實(shí)際面積為537 km2。采用趨勢(shì)面法,在嫩二段坐標(biāo)體系下以鉆孔地理坐標(biāo)為自變量,其對(duì)應(yīng)的厚度趨勢(shì)值為因變量,采用算法軟件MATLAB求得3層油頁(yè)巖段厚度趨勢(shì)值方程,并對(duì)油頁(yè)巖厚度空間展布進(jìn)行擬合,獲得了較為理想的擬合度。在確定子塊面積與厚度的基礎(chǔ)上,對(duì)每個(gè)子塊加權(quán)求和分別得出青一段資源量為12 724億t、嫩一段為8 873.1億t、嫩二段為9 872.4億t。

體積法; 趨勢(shì)面法; ArcGIS; MATLAB; 油頁(yè)巖; 松遼盆地

松遼盆地為大型坳陷盆地,油頁(yè)巖呈層狀大面積平緩分布,品位和厚度變化都比較穩(wěn)定,多位學(xué)者采用體積法進(jìn)行了油頁(yè)巖資源量的估算,但對(duì)計(jì)算中涉及的指標(biāo)取值方法并不統(tǒng)一。

張海龍[1]對(duì)松遼盆地三套油頁(yè)巖進(jìn)行了資源量計(jì)算,預(yù)測(cè)區(qū)邊界由巖礦厚度等值線圖圈定,厚度和面積由油頁(yè)巖厚度等值線圖確定,計(jì)算時(shí)以兩條厚度等值線圈定一個(gè)塊段(內(nèi)插0.7 m厚度等值線),再將各塊段資源累加。王慎余等[2]進(jìn)行油頁(yè)巖預(yù)測(cè)時(shí),含礦面積主要采用湖相(湖沼相)沉積面積,含礦層厚度直接采用或間接引出。2006年由國(guó)土資源部、財(cái)政部、國(guó)家發(fā)改委所啟動(dòng)的新一輪全國(guó)范圍的油頁(yè)巖資源評(píng)價(jià),主要采用體積法進(jìn)行,面積采用勘查區(qū)面積,厚度按油頁(yè)巖厚度等值線劃分若干塊段進(jìn)行,若無厚度等值線,利用勘查區(qū)油頁(yè)巖厚度用資源儲(chǔ)量加權(quán)[3]。

在采用體積法對(duì)松遼盆地油頁(yè)巖資源量進(jìn)行預(yù)測(cè)工作中,面積的確定有的采用湖相面積,有的采用盆地邊界,有的采用勘查區(qū)邊界,并不統(tǒng)一。厚度的確定大多以厚度等值線來的圈定,由于厚度等值線由差值生成,因此計(jì)算中厚度的準(zhǔn)確性值得商榷。

本文探討了體積法運(yùn)用中面積與厚度的確定方法,旨在提供一種較為合理的固體礦產(chǎn)資源量計(jì)算的方法,為其他地區(qū)的相關(guān)工作提供參考。

1 油頁(yè)巖資源量計(jì)算方法

1.1 體積法概述

油頁(yè)巖的資源量與其體積成正比,因此采用體積法估算資源量是可行的。體積法的一般公式:

式中,M為資源量,ˉC為某種礦產(chǎn)的平均含量,V為預(yù)測(cè)對(duì)象的體積。

將油頁(yè)巖礦的礦石量視為M,即認(rèn)為油頁(yè)巖礦石量為油頁(yè)巖資源量,對(duì)礦層劃分成子塊進(jìn)行體積的計(jì)算,所以將公式(1)改寫成如下形式:

式中,M為油頁(yè)巖礦資源量(礦石量),D為油頁(yè)巖礦平均體重,Si為子塊(見下文定義)面積,Hi為子塊內(nèi)可采油頁(yè)巖的平均厚度,i為子塊順序號(hào),SiHi是某子塊油頁(yè)巖礦的體積,SiHi是預(yù)測(cè)區(qū)油頁(yè)巖礦總體積。D為常數(shù),沿用以往D=2t/m3。在這個(gè)計(jì)算過程中,確定了子塊的面積與平均厚度,便可求得任一子塊的體積,進(jìn)而求和獲得各層油頁(yè)巖的資源量[4-6]。

1.2 面積計(jì)算

使用Arc GIS軟件完成面積計(jì)算、坐標(biāo)校正工作,它是一個(gè)完整的地理信息集成系統(tǒng),具有獨(dú)立模塊功能,提供數(shù)據(jù)管理,分析統(tǒng)計(jì)等功能,應(yīng)用領(lǐng)域廣泛[7-9]。

1.2.1 確定資源量計(jì)算區(qū)域 以分布最廣的嫩二段為基準(zhǔn)建立坐標(biāo)體系,以其平面分布圖最西邊做切線建立Y軸,其大致為垂直走向,以平面分布最南邊做切線建立X軸,其大致為水平方向(見圖1)。計(jì)算區(qū)域面積設(shè)定為185 504 km2,青一段、嫩一段完全涵蓋于此計(jì)算區(qū)域內(nèi)。

圖1 松遼盆地油頁(yè)巖資源量計(jì)算區(qū)域圖Fig.1 The caculation area of reserves about oil shale in songliao basin

1.2.2 油頁(yè)巖子塊面積劃分方案與計(jì)算 計(jì)算區(qū)域確定后,以相等間隔在ArcGIS軟件里生成所需格網(wǎng),然后通過手工方式刪除油頁(yè)巖分布范圍外的格網(wǎng)(見圖2)。

圖2 ArcGIS面積子塊劃分方案Fig.2 The procedure about subblock area division in ArcGIS

每個(gè)子塊單元實(shí)際邊長(zhǎng)為23 km,面積為527 km2,但子塊在邊界上或與剝蝕區(qū)相交的區(qū)域由于不規(guī)則均小于527 km2。子塊劃分過細(xì)會(huì)加大工作量,對(duì)于平面分布平緩的油頁(yè)巖層也沒有必要。本次計(jì)算嫩二段劃分為246子塊1 195個(gè)節(jié)點(diǎn)進(jìn)行計(jì)算,青一段劃分147子塊1 007個(gè)節(jié)點(diǎn),嫩一段為197子塊975個(gè)節(jié)點(diǎn)。

1.3 厚度的計(jì)算方法

面積確定后,計(jì)算資源量的最關(guān)鍵因素就是厚度的確定。每一個(gè)正方形子塊的厚度為4個(gè)節(jié)點(diǎn)厚度平均值,但邊界上不規(guī)則子塊的厚度由4～8個(gè)節(jié)點(diǎn)的平均值來控制,以減小誤差。厚度的確定采用多項(xiàng)式趨勢(shì)面分析法,并借助MATLAB算法軟件獲得厚度的趨勢(shì)值方程,可方便地計(jì)算任意節(jié)點(diǎn)的厚度值,進(jìn)而求取任一子塊的平均厚度進(jìn)行體積計(jì)算。

厚度確定原理:多元統(tǒng)計(jì)的趨勢(shì)面分析已經(jīng)被廣泛應(yīng)用于地質(zhì)、地理等多學(xué)科領(lǐng)域,在趨勢(shì)值與剩余值的確定中被證實(shí)是可行的有效的方法。其原理是以一個(gè)多項(xiàng)式曲面來模擬空間點(diǎn)的分布規(guī)律[10-11](見圖3)。

假設(shè)現(xiàn)有n個(gè)觀測(cè)點(diǎn),為得到觀測(cè)點(diǎn)空間分布的規(guī)律,假定空間存在一曲面來擬合原始值的變化規(guī)律,曲面方程為多項(xiàng)式,表達(dá)式為:^zi=a0+a1xi+a2yi+…。

圖3 多項(xiàng)式趨勢(shì)面分析法原理圖Fig.3 The principle map of the polynomial trend surface analysis

令A(yù)=XTX,B=(a1a2…aL)T,C=XT·Z, Z=(Z1Z2…Zn)T,則A·B=C,其中L為該多項(xiàng)式的系數(shù)個(gè)數(shù),其由多項(xiàng)式次數(shù)決定,運(yùn)行MATLAB軟件,解此正規(guī)方程組,得出多項(xiàng)式系數(shù),即可確定厚度趨勢(shì)面的方程。

趨勢(shì)面結(jié)果可經(jīng)過擬合度檢驗(yàn),偏差平方和與回歸平方和的比例一般認(rèn)為達(dá)到40%為佳,擬合指數(shù)公式為:

2 松遼盆地油頁(yè)巖資源量計(jì)算

每一個(gè)子塊資源量求和便得到該層油頁(yè)巖礦的資源量,每一個(gè)子塊的實(shí)際面積在上文面積劃分時(shí)已經(jīng)求得,子塊的厚度計(jì)算較為復(fù)雜,依靠趨勢(shì)面法擬合的多項(xiàng)式曲面方程來計(jì)算。

2.1 油頁(yè)巖厚度計(jì)算結(jié)果

(1)嫩二段厚度計(jì)算結(jié)果

嫩二段分別進(jìn)行了3次與4次多項(xiàng)式的計(jì)算對(duì)比。多項(xiàng)式次數(shù)以滿足計(jì)算要求為原則。厚度3次多項(xiàng)式趨勢(shì)面求解過程共遴選了68口鉆孔來進(jìn)行,鉆孔分布較為均勻,鉆孔坐標(biāo)在ArcGIS環(huán)境下可方便的由大地坐標(biāo)校正到本文設(shè)置的嫩二段直角坐標(biāo)系下,轉(zhuǎn)化后坐標(biāo)數(shù)值較小,方便MATLAB矩陣計(jì)算。根據(jù)趨勢(shì)面原理,3次趨勢(shì)面正規(guī)方程組為68行10列矩陣,68為觀測(cè)點(diǎn)個(gè)數(shù),本文為鉆孔數(shù)量,L為10,即3次多項(xiàng)式的項(xiàng)數(shù)。

在MATLAB命令窗口讀入該矩陣,獲得3次趨勢(shì)面多項(xiàng)式如下:

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。

基于此方程可計(jì)算出任意節(jié)點(diǎn)的厚度趨勢(shì)值,并進(jìn)一步獲得子塊的平均厚度值。調(diào)入以下命令生成該多項(xiàng)式在空間的曲面,即厚度的趨勢(shì)面[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)

依此方法對(duì)嫩二段進(jìn)行4次趨勢(shì)面分析,對(duì)比后,3次多項(xiàng)式趨勢(shì)面更能反映實(shí)際厚度變化趨勢(shì)(見圖4,圖5)。

圖4 嫩二段3次趨勢(shì)面厚度圖Fig.4 Three power polynomial about thickness trend of the 2nd member of Nenjiang group

圖5 嫩二段4次趨勢(shì)面厚度圖Fig.5 Four power polynomial about thickness trend of the 2nd member of Nenjiang group

經(jīng)檢驗(yàn),3次趨勢(shì)面實(shí)際擬合度為R2=52%,較為理想。趨勢(shì)值較好的反應(yīng)了實(shí)際厚度的變化趨勢(shì)(見表1),其高值與低值趨勢(shì)均與實(shí)際對(duì)應(yīng)。

表1 嫩二段厚度趨勢(shì)值與實(shí)際值對(duì)比表Table 1 The comparison of predicted and real thickness of the 2nd member in Nenjiang group

續(xù)表1

(2)青一段厚度計(jì)算結(jié)果

青一段分別進(jìn)行了2、3次趨勢(shì)面多項(xiàng)式擬合,共計(jì)54口鉆孔參與趨勢(shì)面分析。

優(yōu)選3次趨勢(shì)面方程為:

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

由于青山口厚度變化較大,造成數(shù)據(jù)點(diǎn)起伏較大,使得趨勢(shì)面在局部點(diǎn)擬合度不高,對(duì)比具體鉆孔數(shù)據(jù),其趨勢(shì)值整體上能夠反應(yīng)厚度平均值的變化趨勢(shì)(見表2),可以滿足子塊厚度計(jì)算要求。

表2 青一段厚度趨勢(shì)值與實(shí)際值對(duì)應(yīng)表Table 2 The comparison of predicted and real thickness of the 1st member in Qingshankou group

(3)嫩一段厚度計(jì)算結(jié)果

同理,嫩一段進(jìn)行了3次、4次多項(xiàng)式趨勢(shì)面比較,4次趨勢(shì)面不能反映局部厚度的變化,且出現(xiàn)大范圍畸變,最終確定3次趨勢(shì)面方程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,其擬合度達(dá)到35%。

2.2 松遼盆地資源量計(jì)算結(jié)果

以嫩二段資源量計(jì)算為例,計(jì)算面積一共劃分為246個(gè)子塊,將子塊進(jìn)行編號(hào),各個(gè)子塊的面積可以從ArcGIS軟件中讀取。子塊的平均厚度由對(duì)應(yīng)的趨勢(shì)面多項(xiàng)式方程求得。子塊的面積和厚度確定后,可以求取子塊的體積與資源量,將所有246個(gè)子塊資源量求和便獲得嫩二段油頁(yè)巖的資源量(見表3)。

經(jīng)過同樣的計(jì)算過程,可計(jì)算青一段、嫩一段的資源量。嫩二段246個(gè)子塊資源量為9 872.4億t、青一段147個(gè)子塊12 724億t、嫩一段197個(gè)子塊8 873.1億t。

3 結(jié)論

本文提供了一種體積法計(jì)算中面積與厚度的確定方法。子塊面積采用ArcGIS軟件進(jìn)行計(jì)算;子塊厚度計(jì)算依據(jù)趨勢(shì)面法的原理,以鉆孔坐標(biāo)為自變量,鉆孔所揭示的油頁(yè)巖厚度為因變量,采用MATLAB軟件進(jìn)行趨勢(shì)面多項(xiàng)式方程擬合,利用方程計(jì)算厚度值。

松遼盆地資源總量為31 469.5億t,其中青山口組一段油頁(yè)巖資源量規(guī)模最大,12 724億t,占盆地整體40.43%,嫩江組二段次之,9 872.4億t,占盆地整體資源量3 1.3 7%,嫩江組一段資源量為8 873.1億t,占盆地整體資源量28.20%。

表3 嫩二段油頁(yè)巖資源量計(jì)算參數(shù)表Table 3 The reserve calculation parameters of the 2nd member in Nenjiang group

[1] 張海龍.東北北部區(qū)油頁(yè)巖資源評(píng)價(jià)及評(píng)價(jià)方法研究[D].長(zhǎng)春:吉林大學(xué),2008.

[2] 王慎余,許家朋,王振海.我國(guó)油頁(yè)巖資源開發(fā)利用狀況及發(fā)展對(duì)策[J].中國(guó)地質(zhì)經(jīng)濟(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] 劉招君,柳蓉.中國(guó)油頁(yè)巖特征及開發(fā)利用前景[J].地學(xué)前緣,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] 孟慶濤.油頁(yè)巖資源評(píng)價(jià)方法研究—以松遼盆地南部農(nóng)安地區(qū)為例[D].長(zhǎng)春:吉林大學(xué),2007.

[5] 董清水,王立賢,于文斌,等.油頁(yè)巖資源評(píng)價(jià)關(guān)鍵參數(shù)及其求取方法[J].吉林大學(xué)學(xué)報(bào)(地球科學(xué)版),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] 汪云甲,黃宗文.礦產(chǎn)資源評(píng)價(jià)及其應(yīng)用研究[M].北京:中國(guó)礦業(yè)大學(xué)出版社,1998.

[7] 范建偉.基于ArcGIS評(píng)價(jià)油頁(yè)巖資源評(píng)價(jià)—以Piceance盆地綠河組油頁(yè)巖為例[D].北京:中國(guó)地質(zhì)大學(xué)(北京),2010.

[8] 黃波,華韻子.ArcGIS在氣象信息系統(tǒng)中的應(yīng)用[J].大氣科學(xué)研究與應(yīng)用,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] 高鵬,徐志剛,顏遠(yuǎn)欽.基于ArcGIS的縣級(jí)土地利用現(xiàn)狀調(diào)查研究[J].北京測(cè)繪,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] 李小明,劉德民,張品剛,等.趨勢(shì)面分析在東歡坨井田構(gòu)造研究中的應(yīng)用[J].中國(guó)煤炭地質(zhì),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] 孟祥霞,姚清洲,劉建新,等.趨勢(shì)面分析技術(shù)在塔北英買2井區(qū)斷裂識(shí)別中的應(yīng)用[J].新疆石油地質(zhì),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] 路鵬,周榮福,王玫.基于趨勢(shì)面分析的煤層空間分布特征研究[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] 徐建華.計(jì)量地理學(xué)[M].北京:高等教育出版社,2006:28-31.

[14] 張瑞林,李東印,李小軍.趨勢(shì)面擬合法在煤礦瓦斯地質(zhì)變量預(yù)測(cè)中的應(yīng)用研究[J].中國(guó)安全科學(xué)學(xué)報(bào),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] 李哲,張淑英.基于多項(xiàng)式趨勢(shì)面分析理論的天然氣需求預(yù)測(cè)[J].資源與產(chǎn)業(yè),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

長(zhǎng)江學(xué)者和創(chuàng)新團(tuán)隊(duì)發(fā)展計(jì)劃(IRT0864)。

張靜平(1981-),女,博士,工程師,從事非常規(guī)油氣地質(zhì)研究;E-mail:zhangjp69@petrochina.com.cn。