李普山 李偉波 馮智莉

(武漢工程大學(xué)計(jì)算機(jī)科學(xué)與工程學(xué)院 湖北 武漢 430000)

0 引 言

羅布泊鹽湖是我國(guó)超大型的含鉀鹵水礦床，鹽湖由東向西寬90千米，從南到北長(zhǎng)115千米，是世界上最大的干鹽湖之一。其鹵水屬硫酸鎂亞型鹵水，賦存有鉀、鎂、鈉等成分，其品位超過工業(yè)品位，其中鉀鹽近期儲(chǔ)量為1億噸，遠(yuǎn)景儲(chǔ)量為2.5億噸[1]。科學(xué)規(guī)劃鹽田年度的采收方案、優(yōu)化采鹽船的開采路徑、保證硫酸鉀廠吃礦的含鉀品位，對(duì)合理開發(fā)、高效利用不可再生鹽田資源，提高企業(yè)經(jīng)濟(jì)效益，實(shí)現(xiàn)企業(yè)可持續(xù)發(fā)展具有重要意義。而采收方案的基礎(chǔ)就是礦床模型，因此，如何準(zhǔn)確高效地構(gòu)建礦床模型是科學(xué)規(guī)劃采收方案的關(guān)鍵。

為此，本文提出了一種基于Delaunay三角網(wǎng)和并行計(jì)算技術(shù)的礦床三維可視化的分層模型。通過對(duì)整個(gè)鹽池進(jìn)行分層處理來準(zhǔn)確反映不同深度礦床的內(nèi)部情況，同時(shí)運(yùn)用并行技術(shù)對(duì)數(shù)據(jù)計(jì)算進(jìn)行并行處理，使得模型生成速度顯著提高。

1 問題分析

鹽礦開采，實(shí)際上是對(duì)鹵水析出的鹽層進(jìn)行開采，必須知道鹽池中的礦量以及屬性結(jié)構(gòu)。因此需要在相對(duì)固定的時(shí)間周期內(nèi)用測(cè)深儀等進(jìn)行空間數(shù)據(jù)和屬性數(shù)據(jù)的采集，再選擇規(guī)則網(wǎng)絡(luò)結(jié)構(gòu)或不規(guī)則三角網(wǎng)TIN(Triangular Irregular Network)建立鹽礦的數(shù)字高程模型DEM(Digital Elevation Model)。由于鹵水受環(huán)境條件的因素影響，不同深度鹽礦形成的品位和厚度也是變化的，單一的DEM顯然無法反映不同深度的真實(shí)情況。因此，必須建立多層的DEM，而且每一層是可以揭開以反映鹽礦內(nèi)部情況。另外，即使單一的DEM建模都有大量的數(shù)據(jù)計(jì)算，多層DEM構(gòu)建過程中如何提高建模效率也是必須重視的問題。

構(gòu)建可獨(dú)立分層的DEM模型，可以采用化整為零的思想，根據(jù)測(cè)深時(shí)間的不同，將礦床分成不同的鹽層，每一層鹽層利用TIN來擬合曲面，再將各層TIN組合起來就構(gòu)成了整個(gè)鹽礦[2-6]。這樣既可以對(duì)礦床內(nèi)部的結(jié)構(gòu)有所了解，又可以對(duì)整個(gè)礦床不同位置的品位進(jìn)行查詢，使礦量計(jì)算更加準(zhǔn)確。

單層DEM的構(gòu)建采用Delaunay三角剖分算法。經(jīng)過大量學(xué)者證明，在平面內(nèi)任意給定點(diǎn)集的Delaunay三角剖分具有整體最優(yōu)的性質(zhì)，即對(duì)于任意給定的平面點(diǎn)集，Delaunay三角劃分能夠得到整體最優(yōu)的三角形網(wǎng)格，能夠盡可能地避免病態(tài)三角形的出現(xiàn)[7]。為了使三角剖分盡可能優(yōu)化，一般在進(jìn)行三角剖分時(shí)遵循兩個(gè)準(zhǔn)則：(1) 空?qǐng)A特性，在三角網(wǎng)格中任一三角形的外接圓內(nèi)都不存在其他的點(diǎn)；(2) 最大化最小角特征，相鄰三角形構(gòu)成的凸多邊形的對(duì)角線，在交換后，六個(gè)內(nèi)角中最小的內(nèi)角最大[8]。根據(jù)這兩個(gè)準(zhǔn)則形成的三角網(wǎng)格是唯一的，這也是Delaunay三角劃分的劃分標(biāo)準(zhǔn)。

當(dāng)前，構(gòu)建三角網(wǎng)格的算法大致可以分為三類：區(qū)域生長(zhǎng)法、合并分割法(分治法)和逐點(diǎn)插入法。其中：區(qū)域生長(zhǎng)法中主要以三角網(wǎng)生長(zhǎng)為主，三角網(wǎng)生長(zhǎng)法的效率相對(duì)較低，雖然改進(jìn)空間較小，但是實(shí)用性很強(qiáng)；合并分割法的效率最高，但是其運(yùn)用了遞歸的思想，運(yùn)行過程會(huì)占用大量?jī)?nèi)存，不能處理大量數(shù)據(jù)；逐點(diǎn)插入算法比較簡(jiǎn)單，需要內(nèi)存空間不大，但是它的時(shí)間復(fù)雜度相對(duì)較高。

每一次采集的空間數(shù)據(jù)都比較多，同時(shí)多層構(gòu)建DEM運(yùn)算量急劇加大，在盡量保證系統(tǒng)效率的前提下，本文采用了多線程的并行處理辦法，讓不同層的數(shù)據(jù)并行處理，盡可能提高系統(tǒng)效率。

2 方案設(shè)計(jì)

提出本方案的最終目的是為了能夠準(zhǔn)確地計(jì)算礦量，礦量計(jì)算公式為：礦量=體積×品位，故品位和體積是計(jì)算礦量的關(guān)鍵。而前文已經(jīng)提到，受環(huán)境條件的因素影響，不同深度鹽礦形成的品位不一樣，因此需要將整個(gè)鹽礦分層以便對(duì)不同位置的鹽礦品位進(jìn)行處理。

對(duì)于鹽礦體積計(jì)算，這涉及到鹽層間的聯(lián)系(盡管各鹽層構(gòu)網(wǎng)過程相互獨(dú)立，但在進(jìn)行體積計(jì)算時(shí)需要知道相鄰鹽層對(duì)應(yīng)采樣點(diǎn)的過程)，計(jì)算出單位體積內(nèi)的礦量后再進(jìn)行積分求和，從而計(jì)算出整個(gè)礦床的礦量。因此，要求在每一次數(shù)據(jù)采集時(shí)，采樣點(diǎn)的水平面坐標(biāo)、采樣點(diǎn)編號(hào)相同，這樣相鄰兩層三角網(wǎng)格將礦體分割成無數(shù)小三棱柱，計(jì)算每個(gè)小三棱柱體積再積分求和得到每個(gè)礦層體積，而在進(jìn)行查找相鄰鹽層的對(duì)應(yīng)點(diǎn)高程時(shí)也會(huì)簡(jiǎn)單很多，從而降低算法難度。

為使三維仿真模型盡可能逼真，每次獲取的數(shù)據(jù)均比較多。其中存在大量數(shù)據(jù)冗余，而且受測(cè)深航線影響，數(shù)據(jù)分布不夠均勻，深度變化小，需要進(jìn)行合適的處理。為此，本文采用了一種設(shè)定影響圓半徑的方法剔除部分相似數(shù)據(jù)，采用誤差識(shí)別方法剔除異常數(shù)據(jù)。將采樣點(diǎn)集處理之后，對(duì)這些散亂點(diǎn)集進(jìn)行Delaunay三角剖分，即生成Delaunay三角網(wǎng)格。三角網(wǎng)是三維地質(zhì)建模中的一種常用方法，是按一定規(guī)則將離散點(diǎn)拓?fù)溥B接成覆蓋整個(gè)區(qū)域且互不重疊的三角形網(wǎng)格。當(dāng)三角形數(shù)量足夠多時(shí)，可以無限擬合曲面的起伏情況，從而建立了離散點(diǎn)間的空間曲面關(guān)系。三角剖分可以根據(jù)鹵水下鹽層的起伏情況來調(diào)整三角形的大小和數(shù)量，便于瀏覽三維鹽層的礦床模型。

由于每一鹽層的采樣數(shù)據(jù)較多，如逐一對(duì)每一鹽層進(jìn)行三角剖分會(huì)嚴(yán)重影響系統(tǒng)效率，故本文采用了多線程的并行處理的方式，來對(duì)各鹽層并行處理，以此來提高系統(tǒng)效率。對(duì)于多線程技術(shù)，核心問題是任務(wù)調(diào)度，需要從系統(tǒng)的性能評(píng)估入手，即從并行線程的規(guī)模數(shù)、數(shù)據(jù)競(jìng)爭(zhēng)、鎖競(jìng)爭(zhēng)、線程安全、存儲(chǔ)和數(shù)據(jù)組織等方面對(duì)系統(tǒng)進(jìn)行分析和優(yōu)化[9]。一般對(duì)共享數(shù)據(jù)的訪問會(huì)引發(fā)數(shù)據(jù)競(jìng)爭(zhēng)，從而導(dǎo)致結(jié)果的不確定性。為此，需要將各鹽層的采樣點(diǎn)劃分好，使其相互獨(dú)立，從而消除數(shù)據(jù)共享可能引發(fā)的數(shù)據(jù)競(jìng)爭(zhēng)、鎖競(jìng)爭(zhēng)以及存儲(chǔ)問題。并且在程序構(gòu)造上加強(qiáng)結(jié)構(gòu)化和模塊化，從而大大降低了多線程編程的復(fù)雜程度。

3 算法實(shí)現(xiàn)

根據(jù)用戶需求，方案的總體工作流程如圖1所示。

圖1 系統(tǒng)總體功能流程圖

本方案需要對(duì)鹽池采樣點(diǎn)進(jìn)行分析，獲取采樣點(diǎn)的位置信息以及采樣點(diǎn)附近鹽礦樣品的屬性(組份、含量、時(shí)間)，分別建立鹽池采樣點(diǎn)的空間數(shù)據(jù)庫和屬性數(shù)據(jù)庫。根據(jù)空間數(shù)據(jù)庫里的三維坐標(biāo)信息建立鹵水下礦床的三維模型并計(jì)算出鹽層體積，然后利用鹽層體積結(jié)合屬性數(shù)據(jù)庫里的鹽礦品位信息計(jì)算礦量。

3.1 獨(dú)立分層數(shù)據(jù)的組織

根據(jù)系統(tǒng)設(shè)計(jì)方案，本系統(tǒng)數(shù)據(jù)文件包含兩類：圖形數(shù)據(jù)文件和屬性數(shù)據(jù)文件。圖形數(shù)據(jù)文件是可視化的圖形文件，本文選用WRL文件作為系統(tǒng)的空間圖形數(shù)據(jù)文件。同時(shí)該文件也可用于與外部軟件進(jìn)行交換的接口。屬性數(shù)據(jù)文件組織，包含項(xiàng)目數(shù)據(jù)配置文件、測(cè)深儀數(shù)據(jù)文件以及層樣品分析數(shù)據(jù)文件。其中：項(xiàng)目數(shù)據(jù)配置文件用來定義初始數(shù)據(jù)、鏈接采樣數(shù)據(jù)文件和樣品分析數(shù)據(jù)文件；測(cè)深儀數(shù)據(jù)文件記錄的是各采樣點(diǎn)的坐標(biāo)信息以及其編號(hào)；層樣品分析數(shù)據(jù)文件記錄了各采樣點(diǎn)的鹽礦品位信息、坐標(biāo)信息以及采樣時(shí)間。

為了提高檢索效率，方便存儲(chǔ)和維護(hù)鹽池信息，需要建立相關(guān)數(shù)據(jù)庫。表1、表2是根據(jù)分層模型數(shù)據(jù)組織建立的數(shù)據(jù)庫部分情況。表1是按鹽層建立的數(shù)據(jù)庫鹽層表；表2是數(shù)據(jù)庫采樣點(diǎn)表，存儲(chǔ)的信息是鹽層里面的采樣點(diǎn)信息。

表1 數(shù)據(jù)庫鹽層表

表2 數(shù)據(jù)庫采樣點(diǎn)表

3.2 基于三角網(wǎng)生長(zhǎng)法的Delaunay三角剖分

本文采用的是三角網(wǎng)生長(zhǎng)法，相對(duì)于逐點(diǎn)插入法，不需要進(jìn)行大量定位三角形的操作; 相對(duì)于分割合并算法，生長(zhǎng)法也沒有大量的搜尋子集邊集、并網(wǎng)的操作[10-11]。其基本思想為：以一個(gè)三角形為初始三角網(wǎng)，逐漸向外尋找離該三角網(wǎng)各條邊最近的點(diǎn)，與對(duì)應(yīng)邊組成新的三角形，加入到三角網(wǎng)中形成新的三角網(wǎng)。在生長(zhǎng)過程中，要隨時(shí)對(duì)已經(jīng)生長(zhǎng)好的三角網(wǎng)進(jìn)行拓?fù)浣Y(jié)構(gòu)檢測(cè)，確定三角網(wǎng)中各邊之間的關(guān)系，同時(shí)判斷生長(zhǎng)過程是否終止。其具體算法如下：

(1) 在平面點(diǎn)集中選取任意一點(diǎn)(盡量選擇平面點(diǎn)集中間區(qū)域的點(diǎn))，尋找離該點(diǎn)最近的點(diǎn)，形成線段，并以此線段作為定向基線。

(2) 根據(jù)Delaunay三角形的判別標(biāo)準(zhǔn)，將離基線兩端點(diǎn)距離之和最小的點(diǎn)作為新三角形的第三點(diǎn)，并以此確定新三角形。

(3) 把新三角形的兩條邊(新三角形中除去原基線的兩條邊)作新的基線。

(4) 重復(fù)(2)、(3)，直到所有基線都被用過為止。

其具體構(gòu)網(wǎng)流程如圖2所示。

圖2 Delaunay三角剖分三角形生長(zhǎng)法構(gòu)網(wǎng)流程圖

Delaunay三角剖分構(gòu)網(wǎng)過程復(fù)雜，算法代價(jià)較大。因此，本方案利用多線程的并行化思想，在同一坐標(biāo)系下，讓各個(gè)鹽層的離散點(diǎn)集并行構(gòu)建三角網(wǎng)格，完成鹵水下礦層的三維曲面仿真。

3.3 基于多線程的并行化構(gòu)網(wǎng)

并行線程并不是規(guī)模越大，性能越好。運(yùn)用多線程是為了更好地利用系統(tǒng)資源，而線程本身也會(huì)占用系統(tǒng)資源，如果線程數(shù)過多會(huì)嚴(yán)重影響程序本身的性能，過少則不能使系統(tǒng)資源得到充分的利用。另外，線程之間切換時(shí)新線程所需的數(shù)據(jù)需要恢復(fù)，被切換的線程數(shù)據(jù)需要保存，這也能直接導(dǎo)致系統(tǒng)性能降低，故需要在線程數(shù)與系統(tǒng)性能之間找一個(gè)平衡點(diǎn)。

并行處理算法由MFC結(jié)合OpenGL以及OpenMP實(shí)現(xiàn)。OpenMP提供了對(duì)并行算法的高層抽象描述，在程序功能模塊中加入專用的pragma來表明意圖，由編譯器來自動(dòng)將程序并行化，并在可能發(fā)生數(shù)據(jù)沖突之處加入同步互斥及通信[12-14]。其運(yùn)用形式示例如下：

#pragma omp parallel for schedule(dynamic)

//并行調(diào)度分配迭代

for(int k=0;k

{

…

tin.createTin();//構(gòu)網(wǎng)

OutputWrl(sWrlFileName);

//生成WRL模型文件

…

#pragma omp critical//沖突檢查

{

…

WrltoDatabase();

//將層數(shù)據(jù)寫入數(shù)據(jù)庫

}

}

得到的結(jié)果為整個(gè)礦床內(nèi)部各鹽層的三維結(jié)構(gòu)，如圖3所示。

圖3 各鹽層構(gòu)網(wǎng)結(jié)果圖

從圖3可以很明顯看出礦床的各鹽層的三維結(jié)構(gòu)情況，對(duì)構(gòu)網(wǎng)結(jié)果進(jìn)行渲染后的鹽層三維模型如圖4所示。

圖4 礦床分層模型三維顯示圖

系統(tǒng)有以下三個(gè)特性：

(1) 一致性，對(duì)同一鹽層進(jìn)行多次運(yùn)行得到的構(gòu)網(wǎng)結(jié)果基本相同，即組成Delaunay三角形的頂點(diǎn)編號(hào)一致。

(2) 高效性，利用了多線程的并行化思想，讓各個(gè)鹽層并行構(gòu)網(wǎng)，大大提高了系統(tǒng)效率。

(3) 真實(shí)性，能夠較真實(shí)地反映鹵水下礦床的內(nèi)部空間結(jié)構(gòu)以及成礦過程。

4 應(yīng)用結(jié)果分析

本方案可以根據(jù)從勘探現(xiàn)場(chǎng)獲得的測(cè)深數(shù)據(jù)生成鹽池礦床分層的三維仿真模型，反映鹵水下礦床的內(nèi)部屬性結(jié)構(gòu)，方便用戶在進(jìn)行開礦生產(chǎn)時(shí)做出相關(guān)決策。本文選擇C++作為開發(fā)語言，利用面向?qū)ο蟮脑O(shè)計(jì)思想，以MFC作為框架，把點(diǎn)、邊、三角形分別作為一個(gè)新的數(shù)據(jù)結(jié)構(gòu)以對(duì)象的形式進(jìn)行封裝，結(jié)合OpenGL的三維圖形API，較為真實(shí)地實(shí)現(xiàn)了對(duì)鹽層的三維仿真。同時(shí)，為了提高系統(tǒng)效率，利用多線程并行處理的思想，結(jié)合OpenMP編程，讓各鹽層并行構(gòu)網(wǎng)，大大提高了系統(tǒng)效率。表3列出線程數(shù)與系統(tǒng)運(yùn)行時(shí)間的對(duì)比。

表3 不同線程數(shù)與系統(tǒng)運(yùn)行時(shí)間對(duì)比

由此可知，利用多線程優(yōu)化后的程序執(zhí)行效率比不使用多線程時(shí)有顯著提高，并且線程數(shù)對(duì)系統(tǒng)效率影響很大，并不是越多效率越高，當(dāng)線程數(shù)達(dá)到一定數(shù)目時(shí)系統(tǒng)效率反而降低，不能簡(jiǎn)單通過提高線程數(shù)來進(jìn)一步提高效率。

以該項(xiàng)目中5#光鹵石池為例，在相同時(shí)間間隔內(nèi)，前后8次獲取了其鹵水下地形數(shù)據(jù)，取部分?jǐn)?shù)據(jù)進(jìn)行實(shí)驗(yàn)，得到的結(jié)果如圖4所示，共8個(gè)鹽層。通過填充之后能夠真實(shí)地反映其鹵水下的礦床情況，具有真實(shí)性。將其中第一個(gè)鹽層轉(zhuǎn)化為WRL文件格式，其內(nèi)部結(jié)構(gòu)如下，選取17個(gè)點(diǎn)，根據(jù)坐標(biāo)x、y、z排列：

132.777 110 0.625

80.574 99.337 1.875

80.331 40.435 0.875

123.92 46.646 0.875

159.676 40.435 0.375

80.396 10 0

100.69 16.625 0.625

140.102 16.625 0.875

159.782 10.414 1.875

139.891 52.545 0.625

107.395 110 0.875

160 100.373 0.625

80 70.455 0.625

100.69 52.959 0

120.581 85.155 0.375

159.782 70.455 0

80 70.455 0.75

經(jīng)三角剖分后，得到的三角網(wǎng)格如下：其中正數(shù)為三角形的三個(gè)頂點(diǎn)編號(hào)，-1為標(biāo)志號(hào)。

0, 10, 1, -1

0, 14, 1, -1

0, 11, 14, -1

1, 16, 14, -1

14, 15, 11, -1

1, 12, 16, -1

14, 13, 16, -1

14, 9, 15, -1

11, 8, 15, -1

16, 15, 12, -1

14, 3, 13, -1

16, 2, 13, -1

14, 3, 9, -1

15, 4, 9, -1

15, 4, 8, -1

13, 6, 3, -1

13, 6, 2, -1

9, 7, 3, -1

9, 7, 4, -1

8, 7, 4, -1

3, 7, 6, -1

2, 5, 6, -1

8, 5, 7, -1

6, 5, 7, -1

經(jīng)多次實(shí)驗(yàn)，三角剖分結(jié)果基本相同，系統(tǒng)具有一致性。另外，將系統(tǒng)以8個(gè)線程同時(shí)對(duì)8個(gè)鹽層進(jìn)行構(gòu)網(wǎng)，系統(tǒng)效率大大提高。

5 結(jié) 語

本文針對(duì)鹽池鹵水下礦床的基本特征提出了一種礦床分層模型的三維可視化方案。詳細(xì)說明了分層數(shù)據(jù)結(jié)構(gòu)和基于三角網(wǎng)生長(zhǎng)法的Delaunay三角剖分算法構(gòu)建分層DEM方法，并運(yùn)用OpenMP的多線程并行處理技術(shù)和沖突共享機(jī)制，實(shí)現(xiàn)各鹽層快速生成Delaunay三角網(wǎng)后再組合構(gòu)建多層DEM。應(yīng)用結(jié)果表明，這是一種分層快速建模的方法，能夠比較高效準(zhǔn)確地反映礦床內(nèi)部情況，也為復(fù)雜DEM建模提供了一種新的研究手段和方法。