空間巖層屬性的預(yù)估計(jì)算方法

趙小平，王 偉，代佳和，李 捷

（1．西安石油大學(xué) 經(jīng)濟(jì)管理學(xué)院，陜西 西安 710065；2．西安石油大學(xué) 理學(xué)院，陜西 西安 710065）

引 言

空間巖層屬性的準(zhǔn)確預(yù)估對(duì)于研究地質(zhì)巖層結(jié)構(gòu)，開展空間地質(zhì)探索，以及合理安全地開采石油等資源具有十分重大的意義。近年來，人們?cè)诶脝我坏卣饏?shù)進(jìn)行空間巖層屬性的研究和預(yù)測(cè)巖性的分布情況方面取得了一定成果，但用這種方法不能客觀預(yù)測(cè)復(fù)雜、細(xì)微的巖石屬性變化，因而很難準(zhǔn)確描述巖性特征。目前，三維建模是解決空間巖層屬性預(yù)估的有效途徑之一，三維空間的屬性建?？梢杂糜谇蠼饪臻g內(nèi)部物理、化學(xué)等屬性參數(shù)［1-3］，然而，現(xiàn)階段的三維建模研究主要針對(duì)空間中實(shí)體進(jìn)行幾何形態(tài)建模，而對(duì)實(shí)體內(nèi)部屬性的建模研究仍存在一些不足。

三維空間的屬性建模中求解屬性值的較普遍的方法是對(duì)三維空間中的特定屬性進(jìn)行插值處理。在三維空間中選取部分?jǐn)?shù)據(jù)點(diǎn)進(jìn)行測(cè)量得到一定數(shù)量且離散的數(shù)據(jù)，這些數(shù)據(jù)可以部分或者全部反映出空間特征，利用這些已知數(shù)據(jù)進(jìn)行插值，可以預(yù)測(cè)未知數(shù)據(jù)的屬性特征。常見空間插值的方法有：反距離加權(quán)插值算法、樣條函數(shù)插值法和克里金插值法等［4-17］。這些方法普遍受觀察點(diǎn)和數(shù)據(jù)體規(guī)模影響較大，而導(dǎo)致計(jì)算插值速度較慢，插值效率降低，效果不理想。不僅如此，當(dāng)采樣點(diǎn)分布不均勻時(shí)得到的插值結(jié)果誤差較大，且存在誤差迭代積累等諸多因素影響屬性值的準(zhǔn)確性。

本文主要是通過對(duì)已有的空間巖層屬性的預(yù)估方法進(jìn)行研究對(duì)比，對(duì)現(xiàn)有的研究預(yù)測(cè)的數(shù)學(xué)計(jì)算方法加以改進(jìn)，提出一種基于多步迭代屬性插值算法，適用于處理三維空間巖層屬性插值預(yù)估問題和利用計(jì)算機(jī)構(gòu)建三維空間巖層屬性模型。

1 常見預(yù)估計(jì)算方法

1．1 反距離加權(quán)插值算法

地理事物或?qū)傩栽诳臻g分布上互為相關(guān)，存在聚集、隨機(jī)、規(guī)則分布的特征，即兩個(gè)物體間距離越近越相似?；诖?，陳鵬等［1］提出一種反距離加權(quán)插值算法。

反距離加權(quán)法是利用距離對(duì)空間離散數(shù)據(jù)點(diǎn)進(jìn)行確定屬性的一種插值方法。距離預(yù)測(cè)位置較近的測(cè)量值與較遠(yuǎn)的測(cè)量值相比，對(duì)預(yù)測(cè)值的影響更大。該方法為距離預(yù)測(cè)位置較近的點(diǎn)賦的權(quán)重較大，權(quán)重隨距離的增加而下降。

反距離加權(quán)法基于彼此距離較近的事物要比彼此距離較遠(yuǎn)的事物更相似，待估測(cè)點(diǎn)的屬性值為各個(gè)已知點(diǎn)屬性值的加權(quán)平均值的思想，構(gòu)造公式為

式中：fi為已知數(shù)據(jù)點(diǎn)的函數(shù)值；（x，y，z）為插值點(diǎn) 坐 標(biāo)； （） 為 離 散 點(diǎn) 坐 標(biāo)；di（x，y，z） ＝離散點(diǎn)的距離。μ為任意正實(shí)數(shù)，通常取值為2，定義的冪值越高，相鄰數(shù)據(jù)受影響越大，插值數(shù)據(jù)就越接近最近采樣點(diǎn)的值。

利用待估點(diǎn)的屬性值為各個(gè)已知點(diǎn)屬性值的加權(quán)平均值的基本思想，反距離加權(quán)插值算法是一種常用的空間插值方法。

1．2 基于B樣條的插值算法

基于B樣條的插值算法是一種大規(guī)模插值方法，在數(shù)據(jù)量較大時(shí)，B樣條插值算法插值效率較高且插值效果不易受外界因素的影響，插值效果比較理想，其插值效率快且插值效果能很好地逼近散亂數(shù)據(jù)點(diǎn)。該方法無法保證插值曲面能夠完全通過所有插值點(diǎn)，生成的是在一定誤差范圍內(nèi)的擬合曲面，具有一定的適用性。

設(shè)在三維空間中，有一離散點(diǎn)集合P ＝｛（x，y，z，p）｝，坐標(biāo)系 XOYOZ中有一數(shù)據(jù)體 Ω，（）是其中的離散點(diǎn)?？梢詷?gòu)造m×n×k個(gè)雙三次B樣條曲面片集合來逼近表示離散點(diǎn)集P，m表示沿X方向?qū)ⅵ阜殖蒻等分，同理沿Y方向?qū)ⅵ阜殖蒼等分、沿Z方向?qū)ⅵ阜殖蒶等分。由控制點(diǎn)網(wǎng)格體Φ來定義這個(gè)m×n×k個(gè)雙三次B樣條曲面片，Φ均勻地覆蓋在數(shù)據(jù)體Ω中，為（m＋3）×（n＋3）×（k＋3）的控制點(diǎn)網(wǎng)格體。

由這些控制點(diǎn)定義的B樣條函數(shù)為

B樣條函數(shù)插值算法則通過求解覆蓋一定數(shù)據(jù)體的控制點(diǎn)網(wǎng)格體，可對(duì)大規(guī)模且分布不均勻的三維數(shù)據(jù)點(diǎn)進(jìn)行插值［3］。

1．3 滑動(dòng)鄰域克里金插值算法

克里金算法是一種流行的空間插值方法，是國際上公認(rèn)的空間插值方法，也是地質(zhì)統(tǒng)計(jì)學(xué)的主要方法。該方法根據(jù)已知樣本的空間位置和相關(guān)程度，求出未知區(qū)域線性無偏、估計(jì)誤差最小的數(shù)值估計(jì)?？死锝鹚惴ㄔ从诤唵我仔械钠胀死锝?，傳統(tǒng)的克里金法選擇所有的已知點(diǎn)進(jìn)行插值，計(jì)算時(shí)間較長，尤其大規(guī)模采集數(shù)據(jù)的情況下，甚至無能為力；并且在實(shí)際中，真實(shí)的協(xié)方差是估計(jì)出來的，而數(shù)據(jù)測(cè)量時(shí)又不可避免產(chǎn)生誤差，這就有可能導(dǎo)致在大的鄰域上插值產(chǎn)生的均方差比在相對(duì)較小的鄰域上插值時(shí)的均方差要大。故在計(jì)算前對(duì)每個(gè)待插點(diǎn)分別選擇其鄰域內(nèi)的一部分已知點(diǎn)作為原始估值計(jì)算數(shù)據(jù)。因此，如何選取鄰域點(diǎn)成了克里金法亟需解決的問題?；诖?，在普通克里金的基礎(chǔ)上產(chǎn)生了一種高效的滑動(dòng)鄰域克里金方法［7］。以變程步長進(jìn)行矩形網(wǎng)格鄰域劃分，再通過以下步驟實(shí)現(xiàn)：

（1）原始點(diǎn)集合中加入采樣點(diǎn)數(shù)據(jù)。

（2）根據(jù)所有已知點(diǎn)進(jìn)行全局范圍求變程α。

（3）依托上述求得的變程α，采用鄰域劃分方法進(jìn)行正方形網(wǎng)格單元?jiǎng)澐帧?/p>

（4）將已知原始點(diǎn)分配到其各自所屬的網(wǎng)格單元。確定各點(diǎn)所在的鄰域位置。

（5）對(duì)所有待插點(diǎn)進(jìn)行遍歷確定其所在鄰域，運(yùn)用全局變程α建立局部范圍內(nèi)的克里金方程組，通過求解方程組得到待插點(diǎn)結(jié)果。

基于變程的滑動(dòng)鄰域克里金插值方法是一種全局求變程，并根據(jù)已知變程進(jìn)行局部克里金方程組計(jì)算的插值方法，可以在三維空間的曲面高程數(shù)據(jù)插值中得到應(yīng)用。

1．4 基于距離場(chǎng)的屬性插值算法

地理學(xué)第一定律指出“任何事物與其他事物之間都是相關(guān)的，距離越近相關(guān)性越強(qiáng)”。距離是地學(xué)研究中一個(gè)非常重要的空間屬性，常用于構(gòu)建三維地學(xué)變量，并在此基礎(chǔ)上進(jìn)行各類三維空間分析與地學(xué)模擬。根據(jù)表達(dá)側(cè)重點(diǎn)的不同，三維地質(zhì)模型分為實(shí)體模型和場(chǎng)模型兩類。實(shí)體模型旨在顯式構(gòu)建三維空間的地質(zhì)體結(jié)構(gòu)（如斷層、不整合面、礦體等），又稱為幾何模型。場(chǎng)模型是利用規(guī)則的或者不規(guī)則的體元對(duì)三維地學(xué)空間進(jìn)行劃分，用連續(xù)變化的體元屬性值來表示三維空間中的物理和化學(xué)屬性（如孔隙度、品位、電阻率等），又稱為屬性模型。

根據(jù)一系列屬性等值線剖面／斷面數(shù)據(jù)進(jìn)行三維空間插值問題，黃岸爍等［15］提出一種基于距離場(chǎng)的三維空間屬性插值方法。首先，根據(jù)對(duì)應(yīng)屬性等值線輪廓生成三維等值面；然后，構(gòu)建屬性等值面的三維距離場(chǎng)；最后，對(duì)于空間內(nèi)任意點(diǎn)屬性，根據(jù)其兩側(cè)等值面的距離場(chǎng)數(shù)據(jù)進(jìn)行線性插值求取，從而生成整個(gè)三維空間內(nèi)的屬性場(chǎng)模型。該方法可用于揭示地下空間物性分布和預(yù)測(cè)隱伏礦體等。

2 一種空間巖層屬性預(yù)估插值的新算法

2．1 方法

（1）選擇插值區(qū)域

選擇待預(yù)估節(jié)點(diǎn)的平面區(qū)域作為插值區(qū)域。區(qū)域中應(yīng)包含（n＋1）2個(gè)同類型離散插值采樣點(diǎn)及其相關(guān)屬性。

（2）生成插值結(jié)點(diǎn)

對(duì)于區(qū)域中的離散插值采樣點(diǎn) Bi，j＝（xij，yij），i，j＝0，1，2，…，n。其中 xij，yij分別為已知點(diǎn)橫、縱坐標(biāo)，tij為點(diǎn)某一已知巖層屬性值。構(gòu)造相應(yīng)的插值采樣結(jié)點(diǎn) Pi，j＝（xij，yij，tij），i，j＝0，1，2，…，n。

（3）多步迭代屬性插值

對(duì)于區(qū)域中任一點(diǎn)（x0，y0）的巖層屬性值，按如下步驟計(jì)算：

對(duì)于r＝1，2，…，n及i，j＝0，1，…，n-r，令P0，0i，j其中u，v為變量，u，v∈ ［0，1］。

（4）計(jì)算預(yù)估點(diǎn)屬性值

由Pn，n0，0＝（f（u，v），g（u，v），t（u，v）），令f（u，v）＝x0，g（u，v）＝y(tǒng)0，解方程求出u0，v0。代入t（u0，v0）＝t0，t0即為任一點(diǎn)（x0，y0）的巖層屬性值。

2．2 研究實(shí)例

選擇某一待預(yù)估節(jié)點(diǎn)的2 000m×2 000 m平面矩形區(qū)域作為插值區(qū)域。區(qū)域中應(yīng)包含100個(gè)同類型離散插值采樣點(diǎn)及其相關(guān)屬性。對(duì)于區(qū)域中100個(gè)離散插值采樣點(diǎn) Bi，j＝（xij，yij），i，j＝1，2，3，…，9，結(jié)合相應(yīng)的已知巖層屬性值tij，構(gòu)造相應(yīng)的插值采樣結(jié)點(diǎn) Pi，j＝（xij，yij，tij），i，j＝0，1，2，…，9。利用多步迭代屬性插值公式，可解方程求出任一點(diǎn)（x0，y0）的巖層屬性值t0，并作出預(yù)估曲面模擬結(jié)果，如圖1所示。

圖1 多步迭代插值曲面Fig．1 M ulti-step iterative interpolation surface

2．3 本文方法與常見算法的對(duì)比分析

文獻(xiàn)［1］在三維空間中利用反距離加權(quán)插值算法和基于B樣條插值算法進(jìn)行屬性插值建模實(shí)現(xiàn)；文獻(xiàn)［8］說明滑動(dòng)克里金算法無法在完成四維屬性數(shù)據(jù)體插值。文獻(xiàn)［15］將基于距離場(chǎng)的三維空間屬性插值方法應(yīng)用于某銅礦區(qū)的可控源音頻大地電磁法剖面數(shù)據(jù)，插值得到三維空間的視電阻率場(chǎng)。結(jié)合已有算法的解決方案，將各方案進(jìn)行對(duì)比，結(jié)果見表1。

表1 5種預(yù)估方案效果對(duì)比Tab．1 Com parison of five estimation schemes

3 結(jié) 論

通過分析對(duì)比已有空間巖層屬性的預(yù)估計(jì)算方法，本文提出了一種基于多步迭代屬性插值算法，該方法適用于三維空間巖層屬性插值預(yù)估問題。首先，通過選擇插值區(qū)域，生成插值結(jié)點(diǎn)；其次應(yīng)用多步迭代屬性插值，最終達(dá)到計(jì)算預(yù)估點(diǎn)屬性值的目的。基于此可以計(jì)算機(jī)構(gòu)建三維空間巖層屬性模型。本文算法方案對(duì)預(yù)估計(jì)算方法適于計(jì)算機(jī)實(shí)現(xiàn)，可提高預(yù)測(cè)的準(zhǔn)確度，降低計(jì)算復(fù)雜度，是地層勘探與巖層屬性研究的可行的有效途徑。論文研究為進(jìn)一步三維空間屬性建模提供了新的思路和方法。