本發(fā)明涉及gpgpu和三維地理信息可視化技術(shù)領(lǐng)域，特別是涉及一種針對天然氣泄漏擴散條件下的擴散濃度數(shù)據(jù)生成和基于球面的動態(tài)可視化一體化的方法。

背景技術(shù)：

天然氣的泄露擴散事故危及國家能源的運輸安全，往往帶來不可估量的經(jīng)濟損失和人員傷亡。及時準(zhǔn)確地預(yù)測天然氣泄漏擴散后的濃度分布范圍，對于事故后的應(yīng)急決策十分必要。

對于天然氣泄露擴散范圍的可視化表達(dá)，有關(guān)文獻提出的常用方法是結(jié)合gis和大氣擴散模型，張斌才,趙軍.大氣污染擴散的高斯煙羽模型及其gis集成研究[j].環(huán)境監(jiān)測管理與技術(shù),2008,05:17-19+55；歐陽坤.基于三維和時態(tài)gis的大氣污染擴散模擬系統(tǒng)研究與實現(xiàn)[d].清華大學(xué),2011；鄭茂輝,金敏,郭飛.gis支持下的街區(qū)毒氣擴散過程建模與模擬[j].武漢大學(xué)學(xué)報(信息科學(xué)版),2013,08:935-939；簡洪登,范湘濤,王金鑫.地鐵站點有害氣體擴散三維可視化研究[j].測繪科學(xué),2013,01:136-138+141。在結(jié)合gis和大氣擴散模型來表達(dá)天然氣泄漏擴散范圍時，許多研究側(cè)重于泄漏擴散范圍在二維地圖上的可視化；泄漏擴散范圍在三維空間的擴散模擬主要是基于marchingcubes算法繪制濃度等值面的三維可視化，存在的問題在于：(1)濃度等值面只能表示數(shù)據(jù)中部分等值面的結(jié)構(gòu)和空間關(guān)系，不能反映天然氣泄漏擴散過程中濃度分布的全貌和細(xì)節(jié)；(2)對于大規(guī)模的濃度數(shù)據(jù)，濃度等值點的確定和濃度等值面的構(gòu)建影響了渲染效率，因此常常是將濃度等值點數(shù)據(jù)的生成作為預(yù)處理過程，難以實現(xiàn)濃度數(shù)據(jù)生成和可視化的一體化。

體繪制技術(shù)是科學(xué)可視化技術(shù)的一種，能夠產(chǎn)生三維體數(shù)據(jù)的整體圖像，包括體內(nèi)的有用細(xì)節(jié)。因此，體繪制技術(shù)能夠解決上述第一個問題。有關(guān)文獻結(jié)合體繪制技術(shù)實現(xiàn)了三維數(shù)據(jù)場的可視化。侯美華,吳志紅,陳楷民.基于osg的實時煙霧模擬的研究與應(yīng)用[j].計算機工程與設(shè)計,2011,06:2088-2091；胡自和,劉坡,龔建華,王群.基于虛擬地球的臺風(fēng)多維動態(tài)可視化系統(tǒng)的設(shè)計與實現(xiàn)[j].武漢大學(xué)學(xué)報(信息科學(xué)版),2015,10:1299-1305。文獻中用于體繪制的三維數(shù)據(jù)場均為可以直接使用的規(guī)則體數(shù)據(jù)，不涉及體數(shù)據(jù)的生成過程。目前體繪制技術(shù)還沒有應(yīng)用到天然氣泄露擴散模擬方面。

gpgpu是利用gpu進行的通用計算技術(shù)，經(jīng)典gpgpu技術(shù)借助gpu圖形流水線的大規(guī)模并行計算能力，通過紋理映射功能實現(xiàn)通用算法的并行計算，極大地加速了通用計算過程。gpgpu技術(shù)和體繪制技術(shù)的結(jié)合能夠解決上述第二個問題，可以實現(xiàn)擴散濃度數(shù)據(jù)生成和動態(tài)可視化的一體化。卞燕山.三維數(shù)字化戰(zhàn)場中的核擴散表現(xiàn)技術(shù)[d].國防科學(xué)技術(shù)大學(xué),2010。文獻中結(jié)合gpgpu技術(shù)和體繪制技術(shù)實現(xiàn)了三維數(shù)字化戰(zhàn)場中核素的動態(tài)擴散過程。該方法的不足之處在于將gpgpu計算得到核素擴散體數(shù)據(jù)和基于核素擴散體數(shù)據(jù)實現(xiàn)體繪制分裂為了兩個過程，沒有實現(xiàn)體數(shù)據(jù)生成和動態(tài)可視化的一體化。在使用gpgpu計算得到核素擴散體數(shù)據(jù)時沒有考慮三維球面的影響，并且使用固定的三維擴散范圍，難以完整地展示動態(tài)過程中不同時刻擴散范圍的全貌。

相關(guān)術(shù)語：

gpgpu通用計算圖形處理器

gis地理信息系統(tǒng)

gpu圖形處理器

cpu中央處理器

opengl開放圖形庫，是行業(yè)領(lǐng)域中最為廣泛接納的2d/3d圖形api

pbo像素緩沖區(qū)對象

fbo幀緩沖區(qū)對象

技術(shù)實現(xiàn)要素：

本發(fā)明的目的是針對現(xiàn)有的對于天然氣泄漏擴散范圍的三維動態(tài)可視化的不足，提出一種針對天然氣泄漏擴散條件下的基于gpgpu和體繪制的擴散濃度數(shù)據(jù)生成和基于球面的動態(tài)可視化一體化的方法。

本發(fā)明的技術(shù)方案提供一種基于gpgpu的天然氣泄漏擴散的動態(tài)可視化方法，包括以下過程，

首先，初始化繪制環(huán)境，使用opengl創(chuàng)建兩個pbo和一個fbo，創(chuàng)建三維紋理對象和二維紋理對象用于存儲gpgpu的計算結(jié)果，定義擴散坐標(biāo)系和體素坐標(biāo)系，初始化視點、光源在wgs84坐標(biāo)系中的位置；

設(shè)當(dāng)前的擴散時刻為t，先基于cpu進行以下泄露擴散范圍的頂點陣列生成處理，

步驟a1，由cpu在擴散坐標(biāo)系中，計算擴散時刻t時天然氣分別在x軸、y軸、z軸的擴散距離rx、ry、rz和擴散范圍在x軸的最小值minx；

步驟a2，由cpu在擴散坐標(biāo)系中，將擴散范圍內(nèi)的采樣距離與三維場景的lod級別相關(guān)聯(lián)，在當(dāng)前三維場景的lod級別下以deltax、deltay、deltaz的距離分別在x軸、y軸、z軸方向采樣，生成t時刻泄露擴散范圍的頂點陣列，并重新調(diào)整rx、ry、rz，使得x軸、y軸、z軸方向的采樣點個數(shù)均為2的整次冪；將擴散坐標(biāo)系中泄漏擴散范圍的體包圍盒轉(zhuǎn)換到wgs84坐標(biāo)系，繼而轉(zhuǎn)換到世界坐標(biāo)系中，得到世界坐標(biāo)系中泄漏擴散范圍的體包圍盒；由兩個平行平面與世界坐標(biāo)系中泄漏擴散范圍的體包圍盒構(gòu)建擴大之后的新的體包圍盒，其中一個平面經(jīng)過((rx+minx)/2,0,rz)在世界坐標(biāo)系中的對應(yīng)點，另一個平面經(jīng)過世界坐標(biāo)系中泄漏擴散范圍體包圍盒底面的四個頂點；

步驟a3，由cpu將t時刻泄露擴散范圍的頂點陣列載入一個空閑的pbo，令t＝t+deltat，deltat為時間間隔，針對下一個當(dāng)前的擴散時刻返回步驟a1，執(zhí)行步驟a1‐a3，將新的擴散時刻t泄露擴散范圍的頂點陣列載入另一個pbo；

當(dāng)cpu將t時刻泄露擴散范圍的頂點陣列載入pbo后，針對下一擴散時刻執(zhí)行步驟a1‐a3時，由gpu將頂點陣列從pbo傳輸?shù)脚c頂點陣列大小相同的三維紋理對象，并完成以下繪制處理，

繪制兩個長方體，將濃度值體數(shù)據(jù)和頂點光照透射率體數(shù)據(jù)通過紋理映射的方式在片段著色器中實現(xiàn)gpgpu計算，把計算結(jié)果渲染并輸出到fbo，則濃度值體數(shù)據(jù)和頂點光照透射率體數(shù)據(jù)分別作為三維紋理被寫入紋理緩存；繪制一個矩形，將屏幕點顏色值數(shù)據(jù)通過紋理映射的方式，結(jié)合對濃度值體數(shù)據(jù)三維紋理和頂點光照透射率體數(shù)據(jù)三維紋理的采樣，在片段著色器中實現(xiàn)gpgpu計算，計算結(jié)果寫入二維紋理并直接渲染到屏幕，實現(xiàn)了擴散濃度數(shù)據(jù)生成和可視化的一體化。

而且，步驟1中，對于擴散坐標(biāo)系和體素坐標(biāo)系的定義如下，

擴散坐標(biāo)系定義為以泄漏源p0在地面的投影點為原點，以下風(fēng)向方向為x軸，以地面上垂直于x軸方向為y軸，以垂直于地面豎直向上方向為z軸；在擴散坐標(biāo)系中，擴散范圍由長方體包圍盒(minx,‐ry,0)、(rx,‐ry,0)、(rx,ry,0)、(minx,ry,0)、(minx,‐ry,rz)、(rx,‐ry,rz)、(rx,ry,rz)、(minx,ry,rz)確定，rx、ry、rz分別為x軸、y軸、z軸的最大擴散距離，minx為擴散范圍在x軸的最小值；

體素坐標(biāo)系定義為以擴散坐標(biāo)系中(minx,‐ry,0)為原點，x軸、y軸、z軸方向均與擴散坐標(biāo)系相同；在體素坐標(biāo)系中，擴散范圍由長方體包圍盒(0,0,0)、(m‐1,0,0)、(m‐1,n‐1,0)、(0,n‐1,0)、(0,0,l‐1)、(m‐1,0,l‐1)、(m‐1,n‐1,l‐1)、(0,n‐1,l‐1)確定，m‐1、n‐1、l‐1分別為x軸、y軸、z軸的最大擴散距離，m、n、l分別表示x軸、y軸、z軸的采樣點個數(shù)。

而且，步驟a1中，擴散坐標(biāo)系中x軸、y軸、z軸的擴散距離rx、ry、rz、minx的計算方法如下，

(1)計算minx的方法為，

當(dāng)t＜t時，令minx＝0，保存minx值；

當(dāng)t≥t時，

①初始化minx＝x0，x0為預(yù)設(shè)初始值，計算(minx,0,h)處的濃度值c；

②若c＜c0，則令minx＝minx+delta，delta為預(yù)設(shè)步長，重新計算(minx,0,h)處的濃度值c，直到c≥c0，進入③；若c≥c0，則直接進入③；

③若c≥c0，則令minx＝minx-1，重新計算(minx,0,h)處的濃度值c，直到c＜c0，保存minx值；

(2)計算rx的方法為，

①初始化rx＝minx+x0，x0為預(yù)設(shè)初始值，計算(rx,0,h)處的濃度值c；

②若c≥c0，則令rx＝rx+delta，delta為預(yù)設(shè)步長，重新計算(rx,0,h)處的濃度值c，直到c＜c0，進入③；若c＜c0，則直接進入③；

③若c＜c0，則令rx＝rx-1，重新計算(rx,0,h)處的濃度值c，直到c≥c0，保存rx值；(3)計算ry的方法為，

①初始化ry＝0，temp_x＝rx，center_x＝0；temp_x為臨時存儲x值的變量，center_x為ry處對應(yīng)的x值；

②計算(temp_x,ry,h)處的濃度值c，若c≥c0，則令ry＝ry+1，重新計算(temp_x,ry,h)處的濃度值c，直到c＜c0，進入③；若c＜c0，則直接進入③；

③令center_x＝temp_x，temp_x＝temp_x-1，計算(temp_x,ry,h)處的濃度值c，若c＜c0，則保存ry值；若c≥c0，則返回②，重復(fù)②—③；

(4)計算rz的方法為，

①初始化rz＝h，temp_x＝center_x，x_left＝temp_x，x_right＝temp_x；temp_x、x_left、x_right均為臨時存儲x值的變量，center_x為ry處對應(yīng)的x值；

②計算(temp_x,0,rz)處的濃度值c，若c≥c0，則令rz＝rz+1，重新計算(temp_x,0,rz)處的濃度值c，直到c＜c0，進入③；若c＜c0，則直接進入③；

③令x_left＝temp_x-1，x_right＝temp_x+1，判斷(x_left,0,rz)處濃度c1、(x_right,0,rz)處濃度c2和給定濃度值c0的大小，

若c1＜c0且c2＜c0，則保存rz值；

若c1≥c0且c2＜c0，則temp_x＝temp_x-1，返回②，重復(fù)②—③；

若c1＜c0且c2≥c0，則temp_x＝temp_x+1，返回②，重復(fù)②—③。

而且，步驟a2中，將泄露擴散范圍的可視化精度劃分為若干個級別，將擴散范圍內(nèi)的采樣距離與三維場景的lod級別相關(guān)聯(lián)，通過多尺度采樣的策略獲得最優(yōu)的m、n、l值，優(yōu)化紋理映射效率。

而且，設(shè)兩個pbo分別記為pbo0和pbo1，

設(shè)cpu綁定pbo0，進行泄露擴散范圍的頂點陣列生成處理，將擴散時刻t的頂點陣列載入pbo0，gpu將頂點陣列從pbo0傳輸?shù)饺S紋理對象并進行繪制；

gpu將頂點陣列從pbo0傳輸?shù)饺S紋理對象并進行繪制的同時，cpu綁定pbo1,進行泄露擴散范圍的頂點陣列生成處理，將下一擴散時刻t的頂點陣列載入pbo1；

gpu將頂點陣列從pbo1傳輸?shù)饺S紋理對象并進行繪制，同時cpu重新綁定pbo0，如此類推，pbo0、pbo1不斷交替?zhèn)鬏旐旤c陣列到三維紋理對象。

而且，gpu的繪制處理包括以下步驟，

步驟b1，首先繪制一個與三維紋理對象等大的長方體，通過紋理映射功能，將三維紋理對象映射到長方體，通過gpgpu技術(shù)，在頂點著色器中實現(xiàn)泄露擴散范圍的頂點陣列的坐標(biāo)轉(zhuǎn)換，在片段著色器中實現(xiàn)泄露擴散范圍的頂點陣列濃度值的通用計算；將計算結(jié)果渲染并輸出到fbo，則擴散空間的濃度值體數(shù)據(jù)作為新的三維紋理被寫入紋理緩存；

步驟b2，再次繪制一個與三維紋理對象等大的長方體，通過紋理映射功能，將三維紋理對象映射到長方體，通過gpgpu技術(shù)，在頂點著色器中實現(xiàn)泄露擴散范圍的頂點陣列的坐標(biāo)轉(zhuǎn)換，在片段著色器中，由頂點向光源發(fā)出的射線在世界坐標(biāo)系中擴大后的泄漏擴散范圍體包圍盒內(nèi)按步長進行空間點采樣，結(jié)合對步驟b2中三維紋理的采樣，實現(xiàn)泄漏擴散體數(shù)據(jù)頂點光線透射率的通用計算，將計算結(jié)果渲染并輸出到fbo，則泄漏擴散體數(shù)據(jù)的頂點光線透射率作為新的三維紋理被寫入紋理緩存；

步驟b3，最后繪制一個屏幕大小的矩形，通過紋理映射功能，將等大的二維紋理對象映射到矩形，通過gpgpu技術(shù)，在頂點著色器中實現(xiàn)矩形頂點的坐標(biāo)轉(zhuǎn)換，在片段著色器中由視點向屏幕像素點發(fā)出的射線在世界坐標(biāo)系中擴大后的泄漏擴散范圍體包圍盒內(nèi)按步長進行空間點采樣，結(jié)合對步驟b1和步驟b2中三維紋理的采樣，實現(xiàn)屏幕像素點顏色值的計算，最終將矩形輸出至屏幕，實現(xiàn)了t時刻天然氣泄漏擴散范圍的三維可視化。

而且，步驟b2中，在擴大之后的包圍盒內(nèi)采樣，需要判斷采樣點是否位于泄露擴散范圍內(nèi)，包括將采樣點從世界坐標(biāo)系轉(zhuǎn)換到wgs84坐標(biāo)系，若采樣點的高程小于rz對應(yīng)的高程并大于0，則采樣點位于泄漏擴散范圍內(nèi)，否則放棄該采樣點；若采樣點位于泄漏擴散范圍內(nèi)，則將采樣點轉(zhuǎn)換到擴散坐標(biāo)系，再從擴散坐標(biāo)系轉(zhuǎn)換到體素坐標(biāo)系，繼而映射得到該采樣點的三維紋理坐標(biāo)，在紋理坐標(biāo)處對三維紋理采樣。

而且，步驟b3中，在擴大之后的包圍盒內(nèi)采樣，需要判斷采樣點是否位于泄露擴散范圍內(nèi)，包括將采樣點從世界坐標(biāo)系轉(zhuǎn)換到wgs84球坐標(biāo)系，若采樣點的高程小于rz對應(yīng)的高程并大于0，則采樣點位于泄漏擴散范圍內(nèi)，否則放棄該采樣點；若采樣點位于泄漏擴散范圍內(nèi)，則將采樣點轉(zhuǎn)換到擴散坐標(biāo)系，再從擴散坐標(biāo)系轉(zhuǎn)換到體素坐標(biāo)系，繼而映射得到該采樣點的三維紋理坐標(biāo)，在紋理坐標(biāo)處對三維紋理采樣。

本發(fā)明創(chuàng)造性地將濃度值體數(shù)據(jù)和頂點光照透射率體數(shù)據(jù)通過gpgpu計算寫入三維紋理，繼而通過采樣三維紋理完成屏幕像素點顏色值的計算，通過三次連續(xù)的繪制過程實現(xiàn)了擴散濃度數(shù)據(jù)生成和可視化的一體化；將擴散空間的采樣距離與三維場景的lod級別相關(guān)聯(lián)，通過多尺度采樣的策略優(yōu)化了紋理映射效率；利用cpu和gpu異步處理的特性來實現(xiàn)不同時刻擴散范圍頂點陣列的異步生成和傳輸，提高了動態(tài)渲染性能；將擴散坐標(biāo)系中泄漏擴散范圍的體包圍盒轉(zhuǎn)換到了世界坐標(biāo)系中，并在世界坐標(biāo)系中擴大后的泄漏擴散范圍體包圍盒內(nèi)進行采樣，考慮了三維球面的影響，實現(xiàn)了基于球面的可視化。本發(fā)明能夠支持天然氣泄漏擴散條件下的擴散濃度數(shù)據(jù)生成和基于球面的動態(tài)可視化一體化。

附圖說明

圖1為本發(fā)明實施例的擴散坐標(biāo)系示意圖。

圖2為本發(fā)明實施例的體素坐標(biāo)系示意圖。

圖3為本發(fā)明實施例的計算rx和minx的流程圖。

圖4為本發(fā)明實施例的計算ry的流程圖。

圖5為本發(fā)明實施例的計算rz的流程圖。

圖6為本發(fā)明實施例的世界坐標(biāo)系中擴大后的泄漏擴散范圍體包圍盒示意圖。

圖7為本發(fā)明實施例的流程圖。

具體實施方式

以下結(jié)合附圖和實施例詳細(xì)說明本發(fā)明技術(shù)方案。

本發(fā)明中使用的泄露模型為小孔泄露模型，小孔泄漏強度取決于泄漏過程中燃?xì)獾牧魉?，其判斷依?jù)為臨界壓力比cpr：

當(dāng)時，小孔泄露的流速可達(dá)到音速，其泄漏強度為：

當(dāng)時，小孔泄露的流速為亞音速，其泄露強度為：

式中，q—小孔泄漏強度，kg/s；c0—孔口泄露系數(shù)；a—泄漏口面積，m²；p0—環(huán)境壓力，pa；p—泄漏中心處壓力，pa；t—泄漏中心處溫度；z—壓縮因子；r—氣體常數(shù)，j/(kmol*k)；m—氣體摩爾質(zhì)量，kg/mol；k—燃?xì)獾褥刂笖?shù)。

本發(fā)明中基于高斯煙羽擴散模型和高斯煙團擴散模型，使用瞬時點源的移動煙團積分模式來求解連續(xù)點源的天然氣濃度分布：

式中，c(x,y,z,t)—擴散坐標(biāo)系中某點(x,y,z)處泄漏氣體的平均濃度，kg/m³；q—泄漏率，kg/s；σx、σy、σz—下風(fēng)向、橫風(fēng)向和垂直方向上氣體的擴散參數(shù)，m；u—平均風(fēng)速，m/s；t—擴散時刻，s；t—天然氣泄露的總時間，s；h—有效源高，m。

實施例采用本發(fā)明的方法，采用opengl和glsl著色語言，基于osgearth三維渲染引擎實現(xiàn)了天然氣泄漏擴散條件下的基于gpgpu和體繪制的擴散濃度數(shù)據(jù)生成和基于球面的動態(tài)可視化一體化。

本發(fā)明實施例提供的一種基于gpgpu的天然氣泄漏擴散的動態(tài)可視化方法，包括以下過程，

首先，初始化繪制環(huán)境，使用opengl創(chuàng)建兩個pbo和一個fbo，創(chuàng)建三維紋理對象和二維紋理對象用于存儲gpgpu的計算結(jié)果，定義擴散坐標(biāo)系和體素坐標(biāo)系，初始化視點、光源在wgs84坐標(biāo)系中的位置；

設(shè)當(dāng)前的擴散時刻為t，先基于cpu進行以下泄露擴散范圍的頂點陣列生成處理，

步驟a1，由cpu在擴散坐標(biāo)系中，計算擴散時刻t時天然氣分別在x軸、y軸、z軸的擴散距離rx、ry、rz和擴散范圍在x軸的最小值minx；

步驟a2，由cpu在擴散坐標(biāo)系中，將擴散范圍內(nèi)的采樣距離與三維場景的lod級別相關(guān)聯(lián)，在當(dāng)前三維場景的lod級別下以deltax、deltay、deltaz的距離分別在x軸、y軸、z軸方向采樣，生成t時刻泄露擴散范圍的頂點陣列，并重新調(diào)整rx、ry、rz，使得x軸、y軸、z軸方向的采樣點個數(shù)均為2的整次冪；將擴散坐標(biāo)系中泄漏擴散范圍的體包圍盒轉(zhuǎn)換到wgs84坐標(biāo)系，繼而轉(zhuǎn)換到世界坐標(biāo)系中，得到世界坐標(biāo)系中泄漏擴散范圍的體包圍盒；由兩個平行平面與世界坐標(biāo)系中泄漏擴散范圍的體包圍盒構(gòu)建擴大之后的新的體包圍盒，其中一個平面經(jīng)過((rx+minx)/2,0,rz)在世界坐標(biāo)系中的對應(yīng)點，另一個平面經(jīng)過世界坐標(biāo)系中泄漏擴散范圍體包圍盒底面的四個頂點；

步驟a3，由cpu將t時刻泄露擴散范圍的頂點陣列載入一個空閑的pbo，令t＝t+deltat，deltat為時間間隔，針對下一個當(dāng)前的擴散時刻返回步驟a1，執(zhí)行步驟a1‐a3，將新的擴散時刻t泄露擴散范圍的頂點陣列載入另一個pbo；

當(dāng)cpu將t時刻泄露擴散范圍的頂點陣列載入pbo后，針對下一擴散時刻t執(zhí)行步驟a1‐a3時，由gpu將頂點陣列從pbo傳輸?shù)脚c頂點陣列大小相同的三維紋理對象，并完成以下繪制處理，

繪制兩個長方體，將濃度值體數(shù)據(jù)和頂點光照透射率體數(shù)據(jù)通過紋理映射的方式在片段著色器中實現(xiàn)gpgpu計算，把計算結(jié)果渲染并輸出到fbo，則濃度值體數(shù)據(jù)和頂點光照透射率體數(shù)據(jù)分別作為三維紋理被寫入紋理緩存；繪制一個矩形，將屏幕點顏色值數(shù)據(jù)通過紋理映射的方式，結(jié)合對濃度值體數(shù)據(jù)三維紋理和頂點光照透射率體數(shù)據(jù)三維紋理的采樣，在片段著色器中實現(xiàn)gpgpu計算，計算結(jié)果寫入二維紋理并直接渲染到屏幕，實現(xiàn)了擴散濃度數(shù)據(jù)生成和可視化的一體化。

具體實施時，可以采用計算機軟件方式實現(xiàn)自動流程運行，提供以上一種基于gpgpu的天然氣泄漏擴散的動態(tài)可視化方法。

為便于參考實施起見，提供實施例具體流程設(shè)計如下，參見圖7：

步驟1，初始化繪制環(huán)境，使用opengl的pbo、fbo技術(shù)創(chuàng)建兩個pbo和一個fbo，創(chuàng)建三維紋理對象和二維紋理對象用于存儲gpgpu的計算結(jié)果，定義擴散坐標(biāo)系和體素坐標(biāo)系，初始化視點、光源在wgs84坐標(biāo)系中的位置；

初始化步驟可由cpu執(zhí)行。

本發(fā)明定義擴散坐標(biāo)系(如圖1所示)和體素坐標(biāo)系(如圖2所示)：

擴散坐標(biāo)系定義為以泄漏源p0在地面的投影點為原點，以下風(fēng)向方向為x軸，以地面上垂直于x軸方向為y軸，以垂直于地面豎直向上方向為z軸，在擴散坐標(biāo)系中，擴散范圍由長方體包圍盒(minx,-ry,0)、(rx,-ry,0)、(rx,ry,0)、(minx,ry,0)、(minx,-ry,rz)、(rx,-ry,rz)、(rx,ry,rz)、(minx,ry,rz)確定，rx、ry、rz分別為x軸、y軸、z軸的最大擴散距離，minx為擴散范圍在x軸的最小值。

體素坐標(biāo)系定義為以擴散坐標(biāo)系中(minx,-ry,0)為原點，x軸、y軸、z軸方向均與擴散坐標(biāo)系相同，即以下風(fēng)向方向為x軸，以地面上垂直于x軸方向為y軸，以垂直于地面豎直向上方向為z軸。在體素坐標(biāo)系中，擴散范圍由長方體包圍盒(0,0,0)、(m-1,0,0)、(m-1,n-1,0)、(0,n-1,0)、(0,0,l-1)、(m-1,0,l-1)、(m-1,n-1,l-1)、(0,n-1,l-1)確定，m-1、n-1、l-1分別為x軸、y軸、z軸的最大擴散距離，m、n、l分別表示x軸、y軸、z軸的采樣點個數(shù)。

體素坐標(biāo)系中點(i,j,k)對應(yīng)擴散坐標(biāo)系中點(i×deltax+minx,j×deltay-ry,k×deltaz)，其中，i＝0，1，…，m‐1；j＝0，1，…，n‐1；k＝0，1，…，l‐1；m‐1、n‐1、l‐1分別為x軸、y軸、z軸的最大距離，m、n、l分別表示x軸、y軸、z軸的采樣點個數(shù)；deltax、deltay、deltaz為采樣距離。在體素坐標(biāo)系中，擴散范圍由長方體包圍盒(0,0,0)、(m‐1,0,0)、(m‐1,n‐1,0)、(0,n‐1,0)、(0,0,l‐1)、(m‐1,0,l‐1)、(m‐1,n‐1,l‐1)、(0,n‐1,l‐1)確定。

實施例中，初始化繪制環(huán)境，創(chuàng)建兩個pbo，分別為pbo0、pbo1。創(chuàng)建三個三維紋理對象texture0、texture1、texture2，分別用作存儲泄露擴散范圍的頂點陣列、泄露擴散濃度值體數(shù)據(jù)、光線透射率體數(shù)據(jù)。創(chuàng)建一個二維紋理對象texture3用作渲染屏幕。創(chuàng)建一個fbo，用作渲染到紋理。

初始化泄漏源的位置為psource(lon,lat,h)，視點的初始位置為peye(lon,lat,h)，光源的位置為plig(lon,lat,h)，psource、peye、plig均為wgs84坐標(biāo)。將視點位置peye和光源位置plig從wgs84坐標(biāo)系轉(zhuǎn)換到世界坐標(biāo)系，則可得到視點的世界坐標(biāo)eyeworld(xw,yw,zw)和光源的世界坐標(biāo)lightworld(xw,yw,zw)。

具體實施時，可以根據(jù)實際情況自行確定泄漏源、視點以及光源在wgs84坐標(biāo)系中的位置。

步驟2，在擴散坐標(biāo)系中，對于泄漏源psource、天然氣泄露的總時間t、擴散時刻t和低危濃度c0，計算擴散距離rx、ry、rz和minx；

本發(fā)明提出：

對于minx，若擴散時刻t小于天然氣泄露總時間t，則令minx等于0；否則為minx賦一個較小值(本領(lǐng)域技術(shù)人員可預(yù)設(shè)取值)，根據(jù)天然氣濃度計算公式c(x,y,z,t)計算濃度c，使t時刻minx處的濃度c大于指定濃度c0，然后逐漸減小minx，直到對應(yīng)的濃度c小于或等于指定濃度c0；

對于rx，首先為rx賦一個較大值(本領(lǐng)域技術(shù)人員可預(yù)設(shè)取值)，根據(jù)天然氣濃度計算公式c(x,y,z,t)計算濃度c，使t時刻rx處的濃度c小于指定濃度c0，然后逐漸減小rx，直到對應(yīng)的濃度c大于或等于指定濃度c0；

對于ry，從(rx,0,h)處開始，逐漸增大y值，根據(jù)天然氣濃度計算公式c(x,y,z,t)計算濃度c，直到濃度c小于指定濃度c0，保存此時的y為ry。若減小x后，(x,ry,h)處的濃度c小于指定濃度c0，則ry為最大值，否則繼續(xù)增大y值求解ry；

對于rz，求出ry對應(yīng)的x值后，從(x,0,h)處開始，逐漸增大z值，根據(jù)天然氣濃度計算公式c(x,y,z,t)計算濃度c，直到濃度c小于指定濃度c0，保存此時的z為rz。若減小x或增大x時，(x,0,rz)處的濃度c均小于指定濃度c0，則rz為最大值，否則，沿著濃度c大于指定濃度c0的一側(cè)重復(fù)上述操作繼續(xù)求解rz。

實施例中具體實現(xiàn)如下：

初始化天然氣泄露的總時間t、擴散時刻t、低危濃度c0、x0和delta；

如圖3，計算minx的方法為：

當(dāng)t＜t時：令minx＝0，保存minx值；

當(dāng)t≥t時：

①初始化minx＝x0，計算(minx,0,h)處的濃度值c,本領(lǐng)域技術(shù)人員可預(yù)設(shè)初始值x0取值，優(yōu)選地x0的取值范圍為0≤x0≤10；

②若c＜c0，則令minx＝minx+delta(本領(lǐng)域技術(shù)人員可預(yù)設(shè)步長delta取值，例如可取50或100或200)，重新計算(minx,0,h)處的濃度值c，直到c≥c0，進入③；若c≥c0，則直接進入③；

③若c≥c0，則令minx＝minx-1，重新計算(minx,0,h)處的濃度值c，直到c＜c0，保存minx值。

如圖3，計算rx的方法為：

①初始化rx＝minx+x0(本領(lǐng)域技術(shù)人員可預(yù)設(shè)初始值x0取值，x0優(yōu)選的取值范圍為0x0≤10)，計算(rx,0,h)處的濃度值c；

②若c≥c0，則令rx＝rx+delta(本領(lǐng)域技術(shù)人員可預(yù)設(shè)步長delta取值，例如可取50或100或200)，重新計算(rx,0,h)處的濃度值c，直到c＜c0，進入③；若c＜c0，則直接進入③；

③若c＜c0，則令rx＝rx-1，重新計算(rx,0,h)處的濃度值c，直到c≥c0，保存rx值。

本實施例取x0＝10，delta＝100。

如圖4，計算ry的方法為：

①初始化ry＝0，temp_x＝rx，center_x＝0；temp_x為臨時存儲x值的變量，center_x為ry處對應(yīng)的x值；

②計算(temp_x,ry,h)處的濃度值c，若c≥c0，則令ry＝ry+1，重新計算(temp_x,ry,h)處的濃度值c，直到c＜c0，進入③；若c＜c0，則直接進入③；

③令center_x＝temp_x，temp_x＝temp_x-1，計算(temp_x,ry,h)處的濃度值c，若c＜c0，則保存ry值；若c≥c0，則返回②，重復(fù)②—③；

如圖5，計算rz的方法為：

①初始化rz＝h，temp_x＝center_x，x_left＝temp_x，x_right＝temp_x；temp_x、x_left、x_right均為臨時存儲x值的變量，center_x為ry處對應(yīng)的x值；

②計算(temp_x,0,rz)處的濃度值c，若c≥c0，則令rz＝rz+1，重新計算(temp_x,0,rz)處的濃度值c，直到c＜c0，進入③；若c＜c0，則直接進入③；

③令x_left＝temp_x-1，x_right＝temp_x+1，判斷(x_left,0,rz)處濃度c1、(x_right,0,rz)處濃度c2和給定濃度值c0的大?。?/p>

若c1＜c0且c2＜c0，則保存rz值；

若c1≥c0且c2＜c0，則temp_x＝temp_x-1，返回②，重復(fù)②—③；

若c1＜c0且c2≥c0，則temp_x＝temp_x+1，返回②，重復(fù)②—③。

具體實施時，可以根據(jù)實際情況自行確定泄露事故持續(xù)時間t、擴散時刻t和低危濃度c0。

步驟3，在擴散坐標(biāo)系中，將擴散范圍內(nèi)的采樣距離與三維場景的lod級別相關(guān)聯(lián)，在當(dāng)前三維場景的lod級別下以deltax、deltay、deltaz的距離分別在x軸、y軸、z軸方向采樣，生成t時刻泄露擴散范圍的頂點陣列，并重新調(diào)整rx、ry、rz，使得x軸、y軸、z軸方向的采樣點個數(shù)均為2的整次冪。將擴散坐標(biāo)系中泄漏擴散范圍的體包圍盒轉(zhuǎn)換到wgs84坐標(biāo)系，繼而轉(zhuǎn)換到世界坐標(biāo)系中，得到世界坐標(biāo)系中泄漏擴散范圍的體包圍盒(如圖6實線所示)。此時的體包圍盒由兩個球面和四個平面包圍而成，由兩個平行平面與世界坐標(biāo)系中泄漏擴散范圍的體包圍盒可構(gòu)建擴大之后的新的體包圍盒，如圖6虛線所示。

本發(fā)明進一步提出，將泄露擴散范圍的可視化精度劃分為若干個級別，將擴散范圍內(nèi)的采樣距離與三維場景的lod級別相關(guān)聯(lián)，通過多尺度采樣的策略優(yōu)化了紋理映射效率：

當(dāng)前三維場景的lod級別是隨著視點的拉近或遠(yuǎn)離而動態(tài)變化的，使用當(dāng)前三維場景的lod級別對應(yīng)的采樣距離在泄露擴散范圍內(nèi)采樣，生成泄漏擴散范圍的頂點陣列。令指數(shù)k＝log2m(k取整數(shù))，即m＝2^k，同理對n、l重新賦值，使m、n、l均為2的整次冪，并根據(jù)m、n、l重新調(diào)整rx、ry、rz。多尺度采樣的策略在滿足泄漏擴散范圍三維可視化視覺精度的同時，能夠獲得最優(yōu)的m、n、l值，由于m、n、l共同決定了紋理映射時三維紋理的大小，因此通過多尺度采樣的策略可以優(yōu)化紋理映射效率。

步驟3中，在隨時間變化的、精確計算的三維擴散范圍內(nèi)采樣空間點，可以以最小的數(shù)據(jù)量完整地展示動態(tài)過程中不同時刻擴散范圍的全貌：

在擴散坐標(biāo)系中由x軸、y軸、z軸的最大擴散距離rx、ry、rz和擴散范圍在x軸的最小值minx確定的泄露擴散范圍內(nèi)采樣空間點，可以完整地展示某時刻擴散范圍的全貌。并且rx、ry、rz、minx隨著時間的變化也會相應(yīng)地改變，這樣每次都能夠以變化的最小的采樣點數(shù)據(jù)量來精確地表達(dá)泄漏擴散的范圍。

實施例中，將泄露擴散范圍的可視化精度劃分為5個級別，與三維場景的lod級別相關(guān)聯(lián)，用lod1～lod5來表示，即在擴散坐標(biāo)系中，在不同的lod級別下對擴散空間采樣時使用不同的采樣距離：

當(dāng)前三維場景的lod級別為lod5或大于lod5時，deltax＝deltay＝deltaz＝1m；

當(dāng)前三維場景的lod級別為lod4時，deltax＝deltay＝deltaz＝2m；

當(dāng)前三維場景的lod級別為lod3時，deltax＝deltay＝deltaz＝4m；

當(dāng)前三維場景的lod級別為lod2時，deltax＝deltay＝deltaz＝8m；

當(dāng)前三維場景的lod級別為lod1或小于lod1時，deltax＝deltay＝deltaz＝16m。以deltax、deltay、deltaz的距離分別在x、y、z方向采樣空間坐標(biāo)，生成頂點陣列。令m-1＝(rx-minx)/deltax、n-1＝2×ry/deltay、l-1＝rz/deltaz，取k＝log_2m(k取整數(shù))，即m＝2^k，同理對n、l重新賦值，令rx＝(m-1)×deltax+minx、ry＝(n-1)×deltay/2、rz＝(l-1)×deltaz。x軸、y軸、z軸方向的采樣點個數(shù)分別為m、n、l，頂點陣列為三維數(shù)組volumearray[m][n][l]，任意頂點的坐標(biāo)為擴散坐標(biāo)(i×deltax+minx,j×deltay-ry,k×deltaz)，其中，i＝0，1，…，m-1；j＝0，1，…，n-1；k＝0，1，…，l-1。初始化三維數(shù)組的元素值均為0。將擴散坐標(biāo)系中泄露擴散范圍的體包圍盒的頂點(minx,-ry,0)、(rx,-ry,0)、(rx,ry,0)、(minx,ry,0)、(minx,-ry,rz)、(rx,-ry,rz)、(rx,ry,rz)、(minx,ry,rz)分別轉(zhuǎn)換到wgs84坐標(biāo)系中，再轉(zhuǎn)換到世界坐標(biāo)系中，獲得世界坐標(biāo)系中泄漏擴散范圍的體包圍盒頂點p1、p2、p3、p4、p5、p6、p7、p8。將擴散坐標(biāo)系中的坐標(biāo)點((rx+minx)/2,0,rz)轉(zhuǎn)換到wgs84坐標(biāo)系中，再轉(zhuǎn)換到世界坐標(biāo)系中，獲得其在世界坐標(biāo)系中的空間位置pc。如圖6中所示由兩個平行平面(其中一個平面經(jīng)過pc,另一個平面經(jīng)過p1、p2、p3、p4)與世界坐標(biāo)系中泄漏擴散范圍的體包圍盒可構(gòu)建擴大之后的新的體包圍盒。具體實施時，可以根據(jù)實際情況自行確定泄漏擴散范圍的可視化精度級別以及不同精度級別下的采樣距離deltax、deltay、deltaz。

步驟4，針對擴散時刻t，gpu將頂點陣列從pbo傳輸?shù)脚c頂點陣列大小相同的三維紋理對象并執(zhí)行三次連續(xù)的繪制過程，同時cpu將會實現(xiàn)把擴散時刻t＝t+deltat的頂點陣列載入另一個pbo。利用cpu和gpu異步處理的特性，對于擴散時刻t和下一時刻的泄露擴散范圍的頂點陣列，使用兩個pbo交替執(zhí)行紋理源的更新和紋理的復(fù)制，實現(xiàn)不同時刻泄露擴散范圍的頂點陣列的異步生成和繪制。

本發(fā)明提出，使用兩個pbo來實現(xiàn)不同時刻泄漏擴散范圍的頂點陣列的異步生成和繪制：

cpu負(fù)責(zé)執(zhí)行計算擴散時刻t＝t+deltat的擴散范圍rx、ry、rz和minx、生成頂點陣列、將頂點陣列載入pbo的一系列過程,gpu負(fù)責(zé)將頂點陣列從pbo0傳輸?shù)饺S紋理對象以及執(zhí)行三次連續(xù)的繪制過程。

擴散時刻t，gpu將頂點陣列從pbo0傳輸?shù)饺S紋理對象，并執(zhí)行三次連續(xù)的繪制過程。同時，cpu綁定pbo1并且執(zhí)行計算擴散時刻t＝t+deltat的擴散范圍rx、ry、rz和minx、生成頂點陣列、將頂點陣列載入pbo1的一系列過程。當(dāng)擴散時刻t＝t+deltat到來時，gpu將已經(jīng)載入完畢的pbo1傳輸?shù)饺S紋理對象，同時cpu重新綁定pbo0，為pbo0傳送下一擴散時刻的頂點陣列。如此類推，pbo0、pbo1不斷交替?zhèn)鬏旐旤c陣列到三維紋理對象。充分利用cpu和gpu異步處理的特性，cpu將紋理源更新為pbo，gpu從其他pbo復(fù)制紋理并執(zhí)行繪制過程，由于異步dma傳輸，更新和繪制過程可以同時執(zhí)行，提高了動態(tài)渲染性能。

實施例中，設(shè)置時間間隔deltat＝10s。在擴散時刻t，gpu將已經(jīng)載入pbo0的頂點陣列volumearray[m][n][l]傳輸?shù)饺S紋理對象texture0(大小為m×n×l)并執(zhí)行三次連續(xù)的繪制過程。使用pbo可以充分利用cpu和gpu異步處理的特性，opengl可以在pbo0和紋理對象texture0之間執(zhí)行異步dma傳輸，所以在gpu將頂點陣列從pbo0傳輸?shù)饺S紋理對象texture0并執(zhí)行三次連續(xù)的繪制過程的同時，cpu綁定pbo1并且執(zhí)行計算擴散時刻t＝t+deltat的擴散范圍rx、ry、rz和minx、生成頂點陣列的三維數(shù)組volumearray[m][n][l]、將頂點陣列volumearray[m][n][l]載入pbo1的一系列過程。當(dāng)擴散時刻t＝t+deltat到來時，gpu將已經(jīng)載入完畢的pbo1傳輸?shù)饺S紋理對象texture0，同時cpu重新綁定pbo0，為pbo0傳送下一擴散時刻的頂點陣列。如此類推，pbo0、pbo1不斷交替?zhèn)鬏旐旤c陣列到三維紋理對象texture0，可以實現(xiàn)不同時刻泄漏擴散范圍的頂點陣列的異步生成和繪制。充分利用不同處理器之間的異步特性，cpu將紋理源更新為pbo，gpu從其他pbo復(fù)制紋理并執(zhí)行繪制過程，由于異步dma傳輸，更新和繪制過程可以同時執(zhí)行，提高了動態(tài)渲染性能。

具體實施時，可以根據(jù)實際情況自行設(shè)置時間間隔deltat。

步驟5，針對擴散時刻t，首先繪制一個與步驟4中三維紋理對象等大的長方體，通過紋理映射功能，將三維紋理對象映射到長方體，通過gpgpu技術(shù)，在頂點著色器中實現(xiàn)泄露擴散范圍的頂點陣列的坐標(biāo)轉(zhuǎn)換，在片段著色器中實現(xiàn)泄露擴散范圍的頂點陣列濃度值的通用計算。將計算結(jié)果渲染并輸出到fbo，則擴散空間的濃度值體數(shù)據(jù)作為新的三維紋理被寫入紋理緩存；

實施例中，針對擴散時刻t，首先繪制一個m×n×l大小的長方體作為紋理映射的對象，長方體在體素坐標(biāo)系中的頂點坐標(biāo)為(0,0,0)、(m-1,0,0)、(m-1,n-1,0)、(0,n-1,0)、(0,0,l-1)、(m-1,0,l-1)、(m-1,n-1,l-1)、(0,n-1,l-1)，對應(yīng)于三維紋理坐標(biāo)(0,0,0)、(1,0,0)、(1,1,0)、(0,1,0)、(0,0,1)、(1,0,1)、(1,1,1)、(0,1,1)。將三維紋理texture0映射到長方體，利用glsl著色語言編寫針對一個紋理元的計算泄漏擴散濃度的通用算法。

在頂點著色器中，輸入變量為頂點坐標(biāo)gl_vertex、頂點對應(yīng)的三維紋理坐標(biāo)gl_texcoord。對gl_vertex執(zhí)行模型視圖變換、投影變換、視口變換后，轉(zhuǎn)換成了屏幕坐標(biāo)gl_position。進行投影變換時以xoy平面為投影平面，使用正對投影平面的正交投影，進行視口變換時定義了一個m×n大小的視口，以保證投影過程中像素的大小和位置保持不變且像素數(shù)據(jù)不會丟失。頂點著色器的輸出變量為屏幕坐標(biāo)gl_position和紋理坐標(biāo)gl_texcoord。頂點變換后的長方體經(jīng)光柵化后得到m×n×l個片元，每個片元對應(yīng)一個三維紋理坐標(biāo)gl_texcoord。

在片段著色器中，輸入變量為紋理坐標(biāo)gl_texcoord，將紋理坐標(biāo)gl_texcoord(u,v,w)映射到擴散坐標(biāo)系中的空間坐標(biāo)(u×(m-1)×deltax+minx,v×(n-1)×deltay-ry,w×(l-1)×deltaz)，根據(jù)天然氣濃度計算公式c(x,y,z,t)得到該片元對應(yīng)的濃度值c，并將濃度值c存儲為gl_fragcolor的r通道值。片段著色器的輸出變量為gl_fragcolor。

將texture1關(guān)聯(lián)到fbo，把計算結(jié)果渲染并輸出到fbo，則擴散空間的濃度值體數(shù)據(jù)作為三維紋理texture1被寫入紋理緩存，然后將texture1與fbo分離。

步驟6，針對擴散時刻t，再次繪制一個與步驟5中三維紋理對象等大的長方體，通過紋理映射功能，將三維紋理對象映射到長方體，通過gpgpu技術(shù)，在頂點著色器中實現(xiàn)泄露擴散范圍的頂點陣列的坐標(biāo)轉(zhuǎn)換，在片段著色器中，由頂點向光源發(fā)出的射線在世界坐標(biāo)系中擴大后的泄漏擴散范圍體包圍盒內(nèi)以一定的步長進行空間點采樣，結(jié)合對步驟5中三維紋理的采樣，實現(xiàn)泄漏擴散體數(shù)據(jù)頂點光線透射率的通用計算。將計算結(jié)果渲染并輸出到fbo，則泄漏擴散體數(shù)據(jù)的頂點光線透射率作為新的三維紋理被寫入紋理緩存。

本發(fā)明進一步提出，為了實現(xiàn)基于三維球面的可視化，在片段著色器中，頂點向光源發(fā)出的射線是在世界坐標(biāo)系中擴大后的泄漏擴散范圍體包圍盒內(nèi)以一定的步長進行空間點采樣。由于是在擴大之后的包圍盒內(nèi)采樣，因此需要判斷采樣點是否位于泄露擴散范圍內(nèi)。將采樣點從世界坐標(biāo)系轉(zhuǎn)換到wgs84坐標(biāo)系：若采樣點的高程小于rz對應(yīng)的高程并大于0，則采樣點位于泄漏擴散范圍內(nèi)，否則，放棄該采樣點。若采樣點位于泄漏擴散范圍內(nèi)，則將采樣點轉(zhuǎn)換到擴散坐標(biāo)系，再從擴散坐標(biāo)系轉(zhuǎn)換到體素坐標(biāo)系，繼而映射得到該采樣點的三維紋理坐標(biāo)，在紋理坐標(biāo)處對三維紋理采樣。

實施例的步驟6中，再次繪制一個m×n×l大小的長方體作為紋理映射的對象，長方體在體素坐標(biāo)系中的頂點坐標(biāo)為(0,0,0)、(m-1,0,0)、(m-1,n-1,0)、(0,n-1,0)、(0,0,l-1)、(m-1,0,l-1)、(m-1,n-1,l-1)、(0,n-1,l-1)，對應(yīng)于三維紋理坐標(biāo)(0,0,0)、(1,0,0)、(1,1,0)、(0,1,0)、(0,0,1)、(1,0,1)、(1,1,1)、(0,1,1)。將三維紋理texture1映射到長方體。利用glsl著色語言編寫針對一個紋理元的計算頂點光照透射率的通用算法。

在頂點著色器中，輸入變量為頂點坐標(biāo)gl_vertex、頂點對應(yīng)的三維紋理坐標(biāo)gl_texcoord、光源的世界坐標(biāo)lightworld和世界坐標(biāo)系中擴大后的泄漏擴散范圍體包圍盒。對gl_vertex執(zhí)行模型視圖變換、投影變換、視口變換后，轉(zhuǎn)換成了屏幕坐標(biāo)gl_position。進行投影變換時以xoy平面為投影平面，使用正對投影平面的正交投影，進行視口變換時定義了一個m×n大小的視口，以保證投影過程中像素的大小和位置保持不變且像素數(shù)據(jù)不會丟失。頂點著色器的輸出變量為屏幕坐標(biāo)gl_position、紋理坐標(biāo)gl_texcoord、光源的世界坐標(biāo)lightworld和世界坐標(biāo)系中擴大后的泄漏擴散范圍體包圍盒。頂點變換后的長方體經(jīng)光柵化后得到m×n×l個片元，每個片元對應(yīng)一個三維紋理坐標(biāo)gl_texcoord。

在片段著色器中，輸入變量為紋理坐標(biāo)gl_texcoord、光源的世界坐標(biāo)lightworld和世界坐標(biāo)系中泄漏擴散范圍的體包圍盒，將紋理坐標(biāo)gl_texcoord(u,v,w)映射到擴散坐標(biāo)系中的空間坐標(biāo)(u×(m-1)×deltax+minx,v×(n-1)×deltay-ry,w×(l-1)×deltaz)，再從擴散坐標(biāo)系轉(zhuǎn)換到wgs84坐標(biāo)系、從wgs84坐標(biāo)系轉(zhuǎn)換到世界坐標(biāo)系中，得到該紋理元的世界坐標(biāo)vertexworld，由vertexworld和lightworld可得到光線的方向為：

dir＝lightworld-vertexworld

求出光線與世界坐標(biāo)系中擴大后的泄漏擴散范圍體包圍盒相交的入點start和出點end。本步驟中光線的入點start即為頂點vertexworld。從vertexworld開始，在世界坐標(biāo)系中以間距step＝0.1沿著光線方向dir采樣，每次得到的采樣點為：

pos＝vertexworld+dir×(pos+step)

將采樣點pos依次存儲到數(shù)組points中，直到出點end。遍歷數(shù)組points來累積計算頂點的透射率。由于世界坐標(biāo)系中泄漏擴散范圍的體包圍盒是擴大之后的包圍盒，因此需要判斷采樣點是否位于泄露擴散范圍內(nèi)。將采樣點從世界坐標(biāo)系轉(zhuǎn)換到wgs84球坐標(biāo)系：若采樣點的高程小于rz對應(yīng)的高程并且大于0，則采樣點位于泄漏擴散范圍內(nèi)，將其映射為三維紋理坐標(biāo)(u,v,w)；否則，放棄該采樣點。在紋理坐標(biāo)(u,v,w)處對三維紋理texture1采樣得到r通道值(即濃度c)，利用光線透射率計算公式：

transmittance＝transmittance×exp(-ε×step×c)(ε為常數(shù))

累積得到頂點的透射率transmittance，并將transmittance存儲為gl_fragcolor的r通道值。片段著色器的輸出變量為gl_fragcolor。

將texture2關(guān)聯(lián)到fbo，把計算結(jié)果渲染并輸出到fbo，則光線透射率體數(shù)據(jù)作為三維紋理texture2被寫入紋理緩存，然后將texture2與fbo分離；

具體實施時，可以根據(jù)實際情況自行設(shè)置采樣間距step。

步驟7，針對擴散時刻t，最后繪制一個屏幕大小的矩形，通過紋理映射功能，將等大的二維紋理對象映射到矩形，通過gpgpu技術(shù)，在頂點著色器中實現(xiàn)矩形頂點的坐標(biāo)轉(zhuǎn)換，在片段著色器中由視點向屏幕像素點發(fā)出的射線在世界坐標(biāo)系中擴大后的泄漏擴散范圍體包圍盒內(nèi)以一定的步長進行空間點采樣，結(jié)合對步驟5和步驟6中三維紋理的采樣，實現(xiàn)屏幕像素點顏色值的計算，最終將矩形輸出至屏幕，實現(xiàn)了t時刻天然氣泄漏擴散范圍的三維可視化；

本發(fā)明進一步提出：為了實現(xiàn)基于三維球面的可視化，在片段著色器中，由視點向屏幕像素點發(fā)出的射線是在世界坐標(biāo)系中擴大后的泄漏擴散范圍體包圍盒內(nèi)以一定的步長進行空間點采樣，由于是在擴大之后的包圍盒內(nèi)采樣，因此需要判斷采樣點是否位于泄露擴散范圍內(nèi)。將采樣點從世界坐標(biāo)系轉(zhuǎn)換到wgs84球坐標(biāo)系：若采樣點的高程小于rz對應(yīng)的高程并大于0，則采樣點位于泄漏擴散范圍內(nèi)，否則，放棄該采樣點。若采樣點位于泄漏擴散范圍內(nèi)，則將采樣點轉(zhuǎn)換到擴散坐標(biāo)系，再從擴散坐標(biāo)系轉(zhuǎn)換到體素坐標(biāo)系，繼而映射得到該采樣點的三維紋理坐標(biāo)，在紋理坐標(biāo)處對三維紋理采樣。

實施例的步驟7中，最后繪制一個resolx×resoly大小的矩形(resolx×resoly為屏幕分辨率)作為紋理映射的對象，矩形的四個頂點為(0,0)、(resolx-1,0)、(resolx-1,resoly-1)、(0,resoly)，分別對應(yīng)二維紋理坐標(biāo)(0,0)、(1,0)、(1,1)、(0,1)。將二維紋理texture3映射到矩形。利用glsl著色語言編寫針對一個紋理元的計算屏幕像素顏色值的通用算法。

在頂點著色器中，輸入變量為矩形頂點坐標(biāo)gl_vertex、頂點對應(yīng)的二維紋理坐標(biāo)gl_texcoord、視點的世界坐標(biāo)eyeworld和世界坐標(biāo)系中擴大后的泄漏擴散范圍體包圍盒。gl_vertex經(jīng)過模型視圖變換、投影變換、視口變換后，轉(zhuǎn)換成了屏幕坐標(biāo)gl_position。進行投影變換時以xoy平面為投影平面，使用正對投影平面的正交投影，進行視口變換時定義了一個resolx×resoly大小的視口，以保證投影過程中像素的大小和位置保持不變且像素數(shù)據(jù)不會丟失。頂點著色器的輸出變量為屏幕坐標(biāo)gl_position、紋理坐標(biāo)gl_texcoord、視點的世界坐標(biāo)eyeworld和世界坐標(biāo)系中擴大后的泄漏擴散范圍體包圍盒。頂點變換后的矩形經(jīng)光柵化后得到resolx×resoly個片元，每個片元對應(yīng)一個屏幕坐標(biāo)gl_position和紋理坐標(biāo)gl_texcoord。

在片段著色器中，輸入變量為屏幕坐標(biāo)gl_position、紋理坐標(biāo)gl_texcoord、視點的世界坐標(biāo)eyeworld和世界坐標(biāo)系中擴大后的泄漏擴散范圍體包圍盒。由視點的世界坐標(biāo)eyeworld可以得到屏幕中心點的世界坐標(biāo)，繼而得到屏幕坐標(biāo)gl_position處的世界坐標(biāo)screenworld。通過eyeworld和screenworld可確定光線的方向為：

dir＝screenworld‐eyeworld

如果光線與世界坐標(biāo)系中擴大后的泄漏擴散范圍體包圍盒不相交，則放棄這條光線，屏幕像素的顏色值為背景顏色值bcolor；如果相交則求出兩個交點start、end，從start開始，以間距step＝0.1沿著光線方向dir采樣，每次得到的采樣點為：

pos＝start+dir×(pos+step)

將采樣點pos依次存儲到數(shù)組points中，直到end。遍歷數(shù)組points來累積計算屏幕像素點的顏色值。初始化透明度opacity＝1，顏色值rgb的計算公式為：

rgb＝pcolor×lcolor

pcolor為采樣點顏色值，初始化為(1,1,1)，lcolor為光照顏色值，初始化為(1,1,1)。由于世界坐標(biāo)系中泄漏擴散范圍的體包圍盒是擴大之后的包圍盒，因此需要判斷采樣點是否位于泄露擴散范圍內(nèi)。將采樣點從世界坐標(biāo)系轉(zhuǎn)換到wgs84球坐標(biāo)系：若采樣點的高程小于rz對應(yīng)的高程并且大于0，則采樣點位于泄漏擴散范圍內(nèi)，將其映射為三維紋理坐標(biāo)(u,v,w)；否則，放棄該采樣點。在紋理坐標(biāo)(u,v,w)處分別對三維紋理texture1采樣得到r通道值(即濃度c)、對三維紋理texture2采樣得到r通道值(即光照透射率transmittance)，利用：

t＝exp(-ε×step×c)(ε為常數(shù))

opacity＝opacity×t

累積得到屏幕像素發(fā)出的光線的透明度opacity，直到opacity接近0或者到達(dá)end。利用：

累積得到光線上的顏色值color，利用：

pixel＝color+opacity×bcolor

得到屏幕像素的顏色值pixel。停止使用fbo，則最終texture3渲染到幀緩存，實現(xiàn)體繪制。

具體實施時，可以根據(jù)實際情況自行設(shè)置采樣點顏色值plocor和光照顏色值lcolor。

步驟5、步驟6和步驟7中，通過連續(xù)的三次繪制過程實現(xiàn)了擴散濃度數(shù)據(jù)生成和可視化的一體化。

實施例以時間間隔deltat重復(fù)上述過程，實現(xiàn)天然氣泄漏擴散范圍的動態(tài)可視化。當(dāng)天然氣泄漏擴散動態(tài)可視化模擬結(jié)束時，刪除pbo0、pbo1和fbo，刪除紋理對象texture0、texture1、texture2、texture3。

通過以上具體實施可以看出，通過gpgpu技術(shù)來實現(xiàn)濃度值體數(shù)據(jù)和頂點光照透射率體數(shù)據(jù)的通用計算，然后通過采樣三維紋理完成屏幕像素點顏色值的計算，雖然在渲染時需要自定義著色過程而且算法復(fù)雜度上有一定的增加，但是使用了pbo、多尺度采樣等優(yōu)化策略和fbo技術(shù)，并且考慮了三維球面的影響，能夠?qū)崿F(xiàn)天然氣泄漏擴散條件下擴散濃度數(shù)據(jù)生成和基于球面的動態(tài)可視化的一體化并且展示了天然氣泄漏擴散過程中濃度分布的全貌和細(xì)節(jié)。