CN107093207A

CN107093207A - A kind of dynamic and visual method of the natural gas leaking diffusion based on GPGPU

Info

Publication number: CN107093207A
Application number: CN201710237761.2A
Authority: CN
Inventors: 程若桢; 陈静
Original assignee: Wuhan University WHU
Current assignee: Wuhan University WHU
Priority date: 2017-04-12
Filing date: 2017-04-12
Publication date: 2017-08-25
Anticipated expiration: 2037-04-12
Also published as: CN107093207B

Abstract

The present invention provides a kind of dynamic and visual method of the natural gas leaking diffusion based on GPGPU, the array of vertices generation processing of leakage range of scatter is carried out by CPU, after the array of vertices that t reveals range of scatter is loaded into PBO by CPU, when performing array of vertices generation processing for next spreading moment and being loaded into another PBO, array of vertices is transferred to three-D grain object and drawing modification from PBO by GPU, concentration value volume data and per-vertex lighting transmissivity volume data are calculated into write-in three-D grain by GPGPU, screen pixels point color value is completed by sampled three-dimensional texture to calculate, realize the generation of diffusion concentration data and the dynamic and visual integration based on sphere under natural gas leaking diffusion conditions.The present invention uses multi-scale sampling policy optimization texture mapping efficiency, and improving dynamic using the characteristic of CPU and GPU asynchronous process renders performance, illustrates the overall picture and details of concentration distribution in natural gas leaking diffusion process.

Description

A kind of dynamic and visual method of the natural gas leaking diffusion based on GPGPU

Technical field

It is more particularly to a kind of to be directed to natural gas the present invention relates to GPGPU and three-dimensional geographic information visualization technique field Diffusion concentration data generation under the conditions of Release and dispersion and the method for the dynamic and visual integration based on sphere.

Background technology

The leakage diffusion accident of natural gas jeopardizes the transportation safety of national energy, often brings immeasurable economic loss And casualties.The timely and accurately concentration distribution scope after Natural Gas Prediction Release and dispersion, for the Emergency decision after accident It is very necessary.

For the Visualization of natural gas leakage range of scatter, the common method that relevant document is proposed be combine GIS and Model of atmospheric diffusion, Zhang Bincai, Gaussian plume model and its GIS Integrated research [J] the environment prison of Zhao Jun air pollution diffusions Test tube is managed and technology, and 2008,05:17-19+55；Air pollution diffusion simulation systems of the Ou Yangkun based on three peacekeeping Temporal GISs Research is with realizing [D] Tsing-Hua University, 2011；Block toxic gas diffusion process model building under Zheng Maohui, Jin Min, Guo Fei .GIS supports With simulation [J] Wuhan University Journals (information science version), 2013,08:935-939；Simple flood is stepped on, Fan Xiangtao, the prosperous subways of Wang Jin Website harmful gases diffusion three-dimensional visualization studies [J] Surveying and mappings, 2013,01:136-138+141.In combination GIS and greatly Gas diffusion model is expressed during natural gas leaking range of scatter, and many researchs lay particular emphasis on Release and dispersion scope on two-dimensional map Visualization；Release and dispersion scope is mainly based upon Marching Cubes algorithms in the diffusion simulations of three dimensions and draws concentration etc. The three-dimensional visualization in value face, has problems in that：(1) concentration contour surface can only represent in data the structure of part contour surface and Spatial relationship, it is impossible to reflect the overall picture and details of concentration distribution in natural gas leaking diffusion process；(2) for large-scale concentration Data, the determination of concentration equivalent point and the structure of concentration contour surface have impact on rendering efficiency, therefore often by concentration equivalent point The generation of data is used as preprocessing process, it is difficult to realize concentration data generation and visual integration.

Volume Rendering Techniques are one kind of scientific visualization technology, can produce the general image of three-dimensional data, including body Interior useful details.Therefore, Volume Rendering Techniques can solve the problem that above-mentioned first problem.Relevant document combination Volume Rendering Techniques are realized The visualization of 3 d data field.The research and application of Real-time Smoke simulations of Hou Meihua, Wu Zhihong, the Chen Kai people based on OSG [J] computer engineering and design, 2011,06:2088-2091；Recklessly from and, Liu Po, Gong Jianhua, Wang Qun are based on virtual earth The design of typhoon Dynamic and Multi dimensional visualization system and realization [J] Wuhan University Journals (information science version), 2015,10:1299- 1305.3 d data field in document for volume drawing is the Regular datasets that can directly use, and is not related to volume data Generating process.In terms of current Volume Rendering Techniques are also not applied to natural gas leakage diffusion simulations.

GPGPU be using GPU carry out general-purpose computations technology, classical GPGPU technologies by GPU graphics pipelines big rule Mould computation capability, the parallel computation of general-purpose algorithm is realized by texture mapping function, has significantly speeded up general-purpose computations mistake Journey.The combination of GPGPU technologies and Volume Rendering Techniques can solve the problem that above-mentioned Second Problem, it is possible to achieve diffusion concentration data are generated With the integration of dynamic and visual.Nuclear proliferation performance technology [D] science and techniques of defence in Bian Yanshan three-dimensional digitals battlefield University, 2010.GPGPU technologies and Volume Rendering Techniques are combined in document realizes the dynamic diffusion of three-dimensional digital battlefield Radionuclide Process.This method, which is disadvantageous in that to calculate GPGPU, obtains nucleic diffusion volume data and real based on nucleic diffusion volume data Existing volume drawing division is not carried out volume data generation and the integration of dynamic and visual for two processes.Counted using GPGPU Calculate and obtain not accounting for the influence of three-dimensional sphere during nucleic diffusion volume data, and use fixed three-dimensional diffusion scope, it is difficult to The intactly not overall picture of range of scatter in the same time in displaying dynamic process.

Relational language：

GPGPU general-purpose computations graphics processors

GIS GIS-Geographic Information System

GPU graphics processors

CPU central processing units

OpenGL open graphic libraries, are the 2D/3D figures API received the most extensively in industry field

PBO pixel buffers object

FBO frame buffer zones object

The content of the invention

The purpose of the present invention is the deficiency for the existing three dimensional dynamic FEM for natural gas leaking range of scatter, Propose a kind of generation of the diffusion concentration data based on GPGPU and volume drawing being directed under natural gas leaking diffusion conditions and based on ball The method of the dynamic and visual integration in face.

Technical scheme provides a kind of dynamic and visual method of the natural gas leaking diffusion based on GPGPU, bag Include procedure below,

First, initialization draw environment, using OpenGL create two PBO and FBO, create three-D grain object and 2 d texture object is used for the result of calculation for storing GPGPU, definition diffusion coordinate system and voxel coordinate system, initialization viewpoint, light Position of the source in WGS84 coordinate systems；

If current spreading moment is t, the array of vertices generation processing of following leakage range of scatter is first carried out based on CPU,

Step a1, by CPU diffusion coordinate system in, calculate spreading moment t when natural gas respectively x-axis, y-axis, z-axis expansion Dissipate the minimum M inx in x-axis apart from Rx, Ry, Rz and range of scatter；

Step a2, by CPU in diffusion coordinate system, by the LOD ranks of the sampled distance in range of scatter and three-dimensional scenic It is associated, with deltaX, deltaY, deltaZ distance respectively in x-axis, y-axis, z-axis under the LOD ranks of current three-dimensional scenic Direction is sampled, and generation t reveals the array of vertices of range of scatter, and readjusts Rx, Ry, Rz so that x-axis, y-axis, z-axis side To sampled point number be 2 whole power；The body bounding box for spreading Release and dispersion scope in coordinate system is transformed into WGS84 to sit Mark system, is then transformed into world coordinate system, obtains the body bounding box of Release and dispersion scope in world coordinate system；It is parallel by two The body bounding box of plane and Release and dispersion scope in world coordinate system builds the new body bounding box after expanding, and one of them is put down Corresponding points of the face by ((Rx+Minx)/2,0, Rz) in world coordinate system, another plane, which is passed through in world coordinate system, to be leaked Four summits of range of scatter body bounding box bottom surface；

Step a3, is loaded into an idle PBO by the array of vertices that t reveals range of scatter by CPU, makes t=t+ DeltaT, deltaT are time interval, for next current spreading moment return to step a1, perform step a1-a3, will be new Spreading moment t leakage range of scatter array of vertices be loaded into another PBO；

After the array of vertices that t reveals range of scatter is loaded into PBO by CPU, step is performed for next spreading moment During a1-a3, by GPU by array of vertices from PBO be transferred to array of vertices size identical three-D grain object, and complete following Drawing modification,

Two cuboids are drawn, by concentration value volume data and per-vertex lighting transmissivity volume data by way of texture mapping Realize that GPGPU is calculated in fragment shader, result of calculation is rendered and FBO is output to, then concentration value volume data and summit light According to transmissivity volume data texture cache is written into respectively as three-D grain；A rectangle is drawn, by screen point color value data By way of texture mapping, with reference to concentration value volume data three-D grain and per-vertex lighting transmissivity volume data three-D grain Sampling, realizes that GPGPU is calculated in fragment shader, and result of calculation write-in 2 d texture is simultaneously directly rendered into screen, realizes Diffusion concentration data are generated and visual integration.

Moreover, in step 1, be defined as follows for diffusion coordinate system and voxel coordinate system,

Diffusion coordinate system is defined as with source of leaks p₀Subpoint on ground is origin, and following wind direction is x-axis, with ground On face perpendicular to x-axis direction be y-axis, using perpendicular to ground straight up direction as z-axis；In diffusion coordinate system, range of scatter By cuboid bounding box (Minx ,-Ry, 0), (Rx ,-Ry, 0), (Rx, Ry, 0), (Minx, Ry, 0), (Minx ,-Ry, Rz), (Rx ,-Ry, Rz), (Rx, Ry, Rz), (Minx, Ry, Rz) determine, Rx, Ry, Rz be respectively x-axis, y-axis, z-axis maximally diffuse away from From Minx is minimum value of the range of scatter in x-axis；

Voxel coordinate system is defined as to spread in coordinate system (Minx ,-Ry, 0) for origin, x-axis, y-axis, z-axis direction with Spread coordinate system identical；In voxel coordinate system, range of scatter is by cuboid bounding box (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) are determined, M-1, N-1, L- 1 is respectively that x-axis, y-axis, z-axis maximally diffuse distance, and M, N, L represent the sampled point number of x-axis, y-axis, z-axis respectively.

Moreover, in step a1, x-axis, y-axis, diffusion length Rx, Ry, Rz, Minx of z-axis calculating side in diffusion coordinate system Method is as follows,

(1) calculate Minx method be,

As t ＜ T, Minx=0 is made, Minx values are preserved；

As t >=T,

1. Minx=x0 is initialized, x0 is preset initial value, calculates the concentration value c at (Minx, 0, H) place；

If 2. c ＜ c0, make Minx=Minx+delta, delta recalculate (Minx, 0, H) place to preset step-length Concentration value c, until c >=c0, into 3.；If c >=c0, it is directly entered 3.；

If 3. c >=c0, makes Minx=Minx-1, the concentration value c at (Minx, 0, H) place is recalculated, until c ＜ c0, is protected Deposit Minx values；

(2) calculate Rx method be,

1. Rx=Minx+x0 is initialized, x0 is preset initial value, calculates the concentration value c at (Rx, 0, H) place；

If 2. c >=c0, makes Rx=Rx+delta, delta recalculate the concentration value at (Rx, 0, H) place to preset step-length C, until c ＜ c0, into 3.；If 3. c ＜ c0, be directly entered；

If 3. c ＜ c0, make Rx=Rx-1, the concentration value c at (Rx, 0, H) place is recalculated, until c >=c0, Rx is preserved Value；(3) calculate Ry method be,

1. Ry=0, Temp_x=Rx, center_x=0 are initialized；Temp_x is the variable of interim storage x values, Center_x is corresponding x values at Ry；

2. the concentration value c at (Temp_x, Ry, H) place is calculated, if c >=c0, Ry=Ry+1 is made, recalculate (Temp_x, Ry, H) place concentration value c, until c ＜ c0, into 3.；If 3. c ＜ c0, be directly entered；

3. center_x=Temp_x, Temp_x=Temp_x-1 are made, the concentration value c at (Temp_x, Ry, H) place is calculated, if C ＜ c0, then preserve Ry values；If c >=c0, return 2., repeat 2. -3.；

(4) calculate Rz method be,

1. Rz=H, Temp_x=center_x, x_left=Temp_x, x_right=Temp_x are initialized；Temp_x、 X_left, x_right are the variable of interim storage x values, and center_x is corresponding x values at Ry；

2. the concentration value c at (Temp_x, 0, Rz) place is calculated, if c >=c0, Rz=Rz+1 is made, recalculate (Temp_x, 0, Rz) the concentration value c at place, until c ＜ c0, into 3.；If 3. c ＜ c0, be directly entered；

3. x_left=Temp_x-1, x_right=Temp_x+1 are made, (x_left, 0, Rz) place's concentration c 1, (x_ is judged Right, 0, Rz) place's concentration c 2 and given concentration value c₀Size,

If c1 ＜ c0 and c2 ＜ c0, Rz values are preserved；

If c1 >=c0 and c2 ＜ c0, Temp_x=Temp_x-1, return 2., repeat 2. -3.；

If c1 ＜ c0 and c2 >=c0, Temp_x=Temp_x+1, return 2., repeat 2. -3..

Moreover, in step a2, the visualization precision for revealing range of scatter is divided into several ranks, by range of scatter Sampled distance it is associated with the LOD ranks of three-dimensional scenic, optimal M, N, L value is obtained by the strategy of multi-scale sampling, it is excellent Change texture mapping efficiency.

Moreover, set two PBO is designated as PBO0 and PBO1 respectively,

If CPU binds PBO0, the array of vertices generation processing of leakage range of scatter is carried out, by spreading moment t summit battle array List and array of vertices is transferred to three-D grain object from PBO0 and drawn into PBO0, GPU；

While array of vertices is transferred to three-D grain object from PBO0 and drawn by GPU, CPU binding PBO1 enter The array of vertices generation processing of row leakage range of scatter, PBO1 is loaded into by next spreading moment t array of vertices；

Array of vertices is transferred to three-D grain object from PBO1 and drawn by GPU, while CPU binds PBO0 again, And so on, the continuous alternate transport array of vertices of PBO0, PBO1 to three-D grain object.

Moreover, GPU drawing modification comprises the following steps,

Step b1, draws a cuboid big with three-D grain object etc., by texture mapping function, by three-dimensional first Texture object is mapped to cuboid, and by GPGPU technologies, the array of vertices of leakage range of scatter is realized in vertex shader Coordinate Conversion, realizes the general-purpose computations of the array of vertices concentration value of leakage range of scatter in fragment shader；By result of calculation Render and be output to FBO, then the concentration value volume data of diffusion space is written into texture cache as new three-D grain；

Step b2, draws a cuboid big with three-D grain object etc., by texture mapping function, by three-dimensional again Texture object is mapped to cuboid, and by GPGPU technologies, the array of vertices of leakage range of scatter is realized in vertex shader Coordinate Conversion, in fragment shader, the ray sent from summit to light source expand in world coordinate system after Release and dispersion Space point sampling is carried out by step-length in scope body bounding box, with reference to the sampling to three-D grain in step b2, Release and dispersion is realized The general-purpose computations of volume data summit light transmittance, result of calculation is rendered and FBO is output to, then the top of Release and dispersion volume data Point light transmittance is written into texture cache as new three-D grain；

Step b3, finally draws the rectangle of a screen size, by texture mapping function, by etc. big 2 d texture pair As being mapped to rectangle, by GPGPU technologies, the Coordinate Conversion on rectangle summit is realized in vertex shader, in fragment shader In expanded from the ray that viewpoint is sent to screen pixels point in world coordinate system after Release and dispersion scope body bounding box in by Step-length carries out space point sampling, with reference to the sampling to three-D grain in step b1 and step b2, realizes screen pixels point color value Calculating, most rectangle is exported to screen at last, realizes the three-dimensional visualization of t natural gas leaking range of scatter.

Moreover, in step b2, the sampling in the bounding box after expansion is, it is necessary to judge whether sampled point is located at leakage diffusion In the range of, including sampled point from world coordinate system is transformed into WGS84 coordinate systems, if the elevation of sampled point is less than the corresponding height of Rz Journey is simultaneously more than 0, then sample abandons the sampled point in the range of Release and dispersion, otherwise；If sample is in Release and dispersion model In enclosing, then sampled point is transformed into diffusion coordinate system, then voxel coordinate system is transformed into from diffusion coordinate system, then mapping is somebody's turn to do The three-D grain coordinate of sampled point, samples at texture coordinate to three-D grain.

Moreover, in step b3, the sampling in the bounding box after expansion is, it is necessary to judge whether sampled point is located at leakage diffusion In the range of, including sampled point from world coordinate system is transformed into WGS84 spherical coordinate systems, if to be less than Rz corresponding for the elevation of sampled point Elevation is simultaneously more than 0, then sample abandons the sampled point in the range of Release and dispersion, otherwise；If sample is in Release and dispersion In the range of, then sampled point is transformed into diffusion coordinate system, then voxel coordinate system is transformed into from diffusion coordinate system, then mapping is obtained The three-D grain coordinate of the sampled point, samples at texture coordinate to three-D grain.

Concentration value volume data and per-vertex lighting transmissivity volume data are creatively calculated write-in three by the present invention by GPGPU Texture is tieed up, the calculating of screen pixels point color value is then completed by sampled three-dimensional texture, passes through three continuous drawing process Realize the generation of diffusion concentration data and visual integration；By the LOD levels of the sampled distance of diffusion space and three-dimensional scenic Not Xiang Guanlian, the texture mapping efficiency by the policy optimization of multi-scale sampling；Using the characteristic of CPU and GPU asynchronous process come The not asynchronous generation and transmission of range of scatter array of vertices in the same time is realized, dynamic is improved and renders performance；Coordinate system will be spread The body bounding box of middle Release and dispersion scope has been transformed into world coordinate system, and the Release and dispersion in world coordinate system after expansion Sampled in scope body bounding box, it is contemplated that the influence of three-dimensional sphere, realize the visualization based on sphere.The present invention can Support the diffusion concentration data generation under natural gas leaking diffusion conditions and the dynamic and visual integration based on sphere.

Brief description of the drawings

Fig. 1 is the diffusion coordinate system schematic diagram of the embodiment of the present invention.

Fig. 2 is the voxel coordinate system schematic diagram of the embodiment of the present invention.

Fig. 3 is the calculating Rx and Minx of embodiment of the present invention flow chart.

Fig. 4 is the calculating Ry of embodiment of the present invention flow chart.

Fig. 5 is the calculating Rz of embodiment of the present invention flow chart.

Fig. 6 is the Release and dispersion scope body bounding box schematic diagram after expanding in the world coordinate system of the embodiment of the present invention.

Fig. 7 is the flow chart of the embodiment of the present invention.

Embodiment

Technical solution of the present invention is described in detail below in conjunction with drawings and examples.

The leakage model used in the present invention is that aperture reveals model, and small hole leaking intensity depends on combustion gas in leakage process Flow velocity, its basis for estimation be critical pressure ratio CPR：

WhenWhen, the flow velocity of aperture leakage can reach velocity of sound, and its leakage intensity is：

WhenWhen, the flow velocity of aperture leakage is subsonic speed, and its leakage intensity is：

In formula, q-small hole leaking intensity, kg/s；C₀Reveal coefficient in-aperture；A-leakage open area, m²；P₀- environment pressure Power, Pa；P-leakage center pressure, Pa；T-leakage center temperature；Z-compressibility factor；R-gas constant, J/ (kmol* K)；M-gas molar quality, kg/mol；K-combustion gas isentropic index.

Gauss plume dispersion model and Gauss puff diffusion model are based in the present invention, the mobile tobacco group of instantaneous point source is used Integral mode is distributed to solve the concentration of natural gas of continuous point source：

In formula, the mean concentration of certain point (x, y, z) place gas leakage, kg/m in C (x, y, z, t)-diffusion coordinate system³； Q-slip, kg/s；σ_x、σ_y、σ_z- lower wind direction, beam wind are to the diffusion parameter with gas in vertical direction, m；U-average wind Speed, m/s；T-spreading moment, s；The total time of T-natural gas leakage, s；H-effective source height, m.

Embodiment is using method of the invention, using OpenGL and GLSL shading languages, based on osgEarth three-dimensional renderings The generation of the diffusion concentration data based on GPGPU and volume drawing under engine implementation natural gas leaking diffusion conditions and based on sphere Dynamic and visual integration.

A kind of dynamic and visual method of natural gas leaking diffusion based on GPGPU provided in an embodiment of the present invention, including Procedure below,

After the array of vertices that t reveals range of scatter is loaded into PBO by CPU, step is performed for next spreading moment t During a1-a3, by GPU by array of vertices from PBO be transferred to array of vertices size identical three-D grain object, and complete following Drawing modification,

When it is implemented, can realize that automatic flow operation is based on there is provided one of the above using computer software mode The dynamic and visual method of GPGPU natural gas leaking diffusion.

For ease of designing as follows there is provided embodiment idiographic flow for the sake of reference implementation, referring to Fig. 7：

Step 1, environment is drawn in initialization, is created two PBO and FBO using OpenGL PBO, FBO technology, is created Three-D grain object and 2 d texture object are used for the result of calculation for storing GPGPU, definition diffusion coordinate system and voxel coordinate system, Initialize the position of viewpoint, light source in WGS84 coordinate systems；

Initialization step can be performed by CPU.

Present invention definition diffusion coordinate system (as shown in Figure 1) and voxel coordinate system (as shown in Figure 2)：

Diffusion coordinate system is defined as with source of leaks p₀Subpoint on ground is origin, and following wind direction is x-axis, with ground On face perpendicular to x-axis direction be y-axis, using perpendicular to ground straight up direction as z-axis, diffusion coordinate system in, range of scatter By cuboid bounding box (Minx ,-Ry, 0), (Rx ,-Ry, 0), (Rx, Ry, 0), (Minx, Ry, 0), (Minx ,-Ry, Rz), (Rx ,-Ry, Rz), (Rx, Ry, Rz), (Minx, Ry, Rz) determine, Rx, Ry, Rz be respectively x-axis, y-axis, z-axis maximally diffuse away from From Minx is minimum value of the range of scatter in x-axis.

Voxel coordinate system is defined as to spread in coordinate system (Minx ,-Ry, 0) for origin, x-axis, y-axis, z-axis direction with Coordinate system is spread identical, i.e., following wind direction is x-axis, using on ground perpendicular to x-axis direction as y-axis, with perpendicular perpendicular to ground Straight upward direction is z-axis.In voxel coordinate system, range of scatter by cuboid bounding box (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) determine, M-1, N-1, L-1 is respectively that x-axis, y-axis, z-axis maximally diffuse distance, and M, N, L represent the sampled point number of x-axis, y-axis, z-axis respectively.

Voxel coordinate system midpoint (i, j, k) correspondence diffusion coordinate system midpoint (i × deltaX+Minx, j × deltaY-Ry, k × deltaZ), wherein, i=0,1 ..., M-1；J=0,1 ..., N-1；K=0,1 ..., L-1；M-1, N-1, L-1 are respectively x Axle, y-axis, the ultimate range of z-axis, M, N, L represent the sampled point number of x-axis, y-axis, z-axis respectively；deltaX、deltaY、 DeltaZ is sampled distance.In voxel coordinate system, range of scatter by cuboid bounding box (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) are determined.

In embodiment, environment is drawn in initialization, creates two PBO, respectively PBO0, PBO1.Create three three-D grains Object texture0, texture1, texture2, array of vertices, the leakage diffusion for being used separately as storage leakage range of scatter are dense Angle value volume data, light transmittance volume data.A 2 d texture object texture3 is created to be used as rendering screen.Create one FBO, as render to texture.

The position for initializing source of leaks is p_source(lon, lat, h), the initial position of viewpoint is p_eye(lon, lat, h), The position of light source is p_lig(lon, lat, h), p_source、p_eye、p_ligIt is WGS84 coordinates.By viewpoint position p_eyeAnd light source position p_ligWorld coordinate system is transformed into from WGS84 coordinate systems, then can obtain the world coordinates Eye of viewpoint_world(x_w,y_w,z_w) and light source World coordinates Light_world(x_w,y_w,z_w)。

When it is implemented, can voluntarily be determined source of leaks, viewpoint and light source in WGS84 coordinate systems according to actual conditions Position.

Step 2, in diffusion coordinate system, for source of leaks p_source, total time T of natural gas leakage, spreading moment t and Low danger concentration c₀, calculate diffusion length Rx, Ry, Rz and Minx；

The present invention is proposed：

For Minx, if spreading moment t is less than natural gas leakage total time T, Minx is made to be equal to 0；Otherwise assigned for Minx One smaller value (the predeterminable value of those skilled in the art), concentration is calculated according to concentration of natural gas calculation formula C (x, y, z, t) C, makes the concentration c at t Minx be more than prescribed concentration c0, Minx is then gradually reduced, until corresponding concentration c is less than or waits In prescribed concentration c0；

For Rx, a higher value (the predeterminable value of those skilled in the art) is assigned for Rx first, according to concentration of natural gas Calculation formula C (x, y, z, t) calculates concentration c, the concentration c at t Rx is less than prescribed concentration c0, is then gradually reduced Rx, Until corresponding concentration c is more than or equal to prescribed concentration c0；

For Ry, start from (Rx, 0, H) place, gradually increase y values, according to concentration of natural gas calculation formula C (x, y, z, t) Concentration c is calculated, until concentration c is less than prescribed concentration c0, the y preserved now is Ry.If reducing after x, the concentration c at (x, Ry, H) place Less than prescribed concentration c0, then Ry is maximum, otherwise continues to increase y values solution Ry；

For Rz, obtain after the corresponding x values of Ry, start from (x, 0, H) place, gradually increase z values, according to concentration of natural gas meter Calculate formula C (x, y, z, t) and calculate concentration c, until concentration c is less than prescribed concentration c0, the z preserved now is Rz.If reducing x or increasing During big x, the concentration c at (x, 0, Rz) place is respectively less than prescribed concentration c0, then Rz is maximum, otherwise, is more than along concentration c and specifies dense The side for spending c0 repeats aforesaid operations continuation solution Rz.

It is implemented as follows in embodiment：

Initialize total time T, spreading moment t, low danger concentration c 0, x0 and the delta of natural gas leakage；

Such as Fig. 3, the method for calculating Minx is：

As t ＜ T：Minx=0 is made, Minx values are preserved；

As t >=T：

1. Minx=x0 is initialized, the concentration value c at (Minx, 0, H) place, the predeterminable initial value of those skilled in the art is calculated The span of x0 values, preferably x0 is 0≤x0≤10；

If 2. c ＜ c0, make Minx=Minx+delta (the predeterminable step-length delta values of those skilled in the art, for example Desirable 50 or 100 or 200), the concentration value c at (Minx, 0, H) place is recalculated, until c >=c0, into 3.；If c >=c0, directly Tap into 3.；

If 3. c >=c0, makes Minx=Minx-1, the concentration value c at (Minx, 0, H) place is recalculated, until c ＜ c0, is protected Deposit Minx values.

Such as Fig. 3, the method for calculating Rx is：

1. Rx=Minx+x0 (the predeterminable initial value x0 values of those skilled in the art, span preferred x0 are initialized For 0x0≤10), the concentration value c at (Rx, 0, H) place is calculated；

If 2. c >=c0, making Rx=Rx+delta, (the predeterminable step-length delta values of those skilled in the art for example can use 50 or 100 or 200), the concentration value c at (Rx, 0, H) place is recalculated, until c ＜ c0, into 3.；If c ＜ c0, are directly entered ③；

If 3. c ＜ c0, make Rx=Rx-1, the concentration value c at (Rx, 0, H) place is recalculated, until c >=c0, Rx is preserved Value.

The present embodiment takes x0=10, delta=100.

Such as Fig. 4, the method for calculating Ry is：

Such as Fig. 5, the method for calculating Rz is：

3. x_left=Temp_x-1, x_right=Temp_x+1 are made, (x_left, 0, Rz) place's concentration c 1, (x_ is judged Right, 0, Rz) place's concentration c 2 and given concentration value c₀Size：

If c1 ＜ c0 and c2 ＜ c0, Rz values are preserved；

If c1 >=c0 and c2 ＜ c0, Temp_x=Temp_x-1, return 2., repeat 2. -3.；

If c1 ＜ c0 and c2 >=c0, Temp_x=Temp_x+1, return 2., repeat 2. -3..

When it is implemented, leakage duration of fault T, spreading moment t and low danger can voluntarily be determined according to actual conditions Concentration c₀。

It is step 3, in diffusion coordinate system, the sampled distance in range of scatter is associated with the LOD ranks of three-dimensional scenic, Adopted respectively in x-axis, y-axis, z-axis direction with deltaX, deltaY, deltaZ distance under the LOD ranks of current three-dimensional scenic Sample, generation t reveals the array of vertices of range of scatter, and readjust Rx, Ry, Rz so that x-axis, y-axis, z-axis direction are adopted Sampling point number is 2 whole power.The body bounding box for spreading Release and dispersion scope in coordinate system is transformed into WGS84 coordinate systems, Then be transformed into world coordinate system, obtain the body bounding box of Release and dispersion scope in world coordinate system (such as Fig. 6 is shown in solid). Body bounding box now is surrounded and formed by two spheres and four planes, is leaked and is expanded in two parallel planes and world coordinate system The new body bounding box after expanding can be built by dissipating the body bounding box of scope, as shown in Fig. 6 dotted lines.

Present invention further propose that, the visualization precision for revealing range of scatter is divided into several ranks, model will be spread Sampled distance in enclosing is associated with the LOD ranks of three-dimensional scenic, the texture mapping effect by the policy optimization of multi-scale sampling Rate：

The LOD ranks of current three-dimensional scenic be with viewpoint drawing close or remote from and dynamic change, use current three-dimensional The corresponding sampled distance of LOD ranks of scene is sampled in leakage range of scatter, generates the array of vertices of Release and dispersion scope.Order Index K=log₂M (K round numbers), i.e. M=2^K, similarly to N, L again assignment, make the whole power that M, N, L are 2, and according to M, N, L readjust Rx, Ry, Rz.The strategy of multi-scale sampling is meeting the same of Release and dispersion scope three-dimensional visualization visual acuity When, optimal M, N, L value is resulted in, the size of three-D grain when having together decided on texture mapping due to M, N, L, therefore pass through The strategy of multi-scale sampling can optimize texture mapping efficiency.

In step 3, the sample space point in the range of three-dimensional diffusion changing over time, accurately calculating can be with minimum Data volume intactly show in dynamic process not the overall picture of range of scatter in the same time：

Diffusion coordinate system in by x-axis, y-axis, z-axis maximally diffuse apart from Rx, Ry, Rz and range of scatter in x-axis most Sample space point in the leakage range of scatter that small value Minx is determined, can intactly show the overall picture of certain moment range of scatter.And , so every time can be with the minimum sampling of change and the change of Rx, Ry, Rz, Minx over time also can correspondingly change Point data amount accurately expresses the scope of Release and dispersion.

In embodiment, the visualization precision for revealing range of scatter is divided into 5 ranks, the LOD ranks with three-dimensional scenic It is associated, represented with LOD1~LOD5, i.e., in diffusion coordinate system, when being sampled under different LOD ranks to diffusion space Use different sampled distances：

When the LOD ranks of current three-dimensional scenic are for LOD5 or more than LOD5, deltaX=deltaY=deltaZ=1m；

When the LOD ranks of current three-dimensional scenic are LOD4, deltaX=deltaY=deltaZ=2m；

When the LOD ranks of current three-dimensional scenic are LOD3, deltaX=deltaY=deltaZ=4m；

When the LOD ranks of current three-dimensional scenic are LOD2, deltaX=deltaY=deltaZ=8m；

When the LOD ranks of current three-dimensional scenic are for LOD1 or less than LOD1, deltaX=deltaY=deltaZ=16m. With deltaX, deltaY, deltaZ distance respectively in x, y, z direction sample space coordinate, array of vertices is generated.Make M-1= (Rx-Minx)/deltaX, N-1=2 × Ry/deltaY, L-1=Rz/deltaZ, takes K=log_2M (K round numbers), i.e. M=2 ^K, similarly to N, L again assignment, make Rx=(M-1) × deltaX+Minx, Ry=(N-1) × deltaY/2, Rz=(L-1) × deltaZ.X-axis, y-axis, the sampled point number in z-axis direction are respectively M, N, L, and array of vertices is three-dimensional array VolumeArray [M] [N] [L], the coordinate on any summit is diffusion coordinate (i × deltaX+Minx, j × deltaY-Ry, k × deltaZ), its In, i=0,1 ..., M-1；J=0,1 ..., N-1；K=0,1 ..., L-1.The element value for initializing three-dimensional array is 0.It will expand Dissipate the summit (Minx ,-Ry, 0) of body bounding box of leakage range of scatter in coordinate system, (Rx ,-Ry, 0), (Rx, Ry, 0), (Minx, Ry, 0), (Minx ,-Ry, Rz), (Rx ,-Ry, Rz), (Rx, Ry, Rz), (Minx, Ry, Rz) are transformed into WGS84 respectively In coordinate system, reconvert into world coordinate system, obtain world coordinate system in Release and dispersion scope body bounding box summit p1, p2、p3、p4、p5、p6、p7、p8.The coordinate points ((Rx+Minx)/2,0, Rz) spread in coordinate system are transformed into WGS84 coordinates In system, reconvert obtains its locus pc in world coordinate system into world coordinate system.As shown in Figure 6 by two Parallel plane (one of plane passes through pc, and another plane passes through p1, p2, p3, p4) and Release and dispersion in world coordinate system The body bounding box of scope can build the new body bounding box after expanding.When it is implemented, can be voluntarily true according to actual conditions Determine the visualization precision grade and sampled distance deltaX under different accuracy rank of Release and dispersion scope, deltaY, deltaZ。

Step 4, for spreading moment t, array of vertices is transferred to three-dimensional with array of vertices size identical by GPU from PBO Texture object simultaneously performs three continuous drawing process, while CPU will realize the summit battle array spreading moment t=t+deltaT List into another PBO.Using the characteristic of CPU and GPU asynchronous process, the leakage for spreading moment t and subsequent time is spread The array of vertices of scope, the renewal in texture source and the duplication of texture are alternately performed using two PBO, realize that leakage is not expanded in the same time Dissipate the asynchronous generation and drafting of the array of vertices of scope.

The present invention is proposed, the asynchronous generation of the array of vertices of Release and dispersion scope in the same time is realized not using two PBO And drafting：

CPU is responsible for performing range of scatter Rx, Ry, Rz and Minx for calculating spreading moment t=t+deltaT, generation summit battle array Row, by array of vertices be loaded into PBO a series of processes, GPU be responsible for by array of vertices from PBO0 be transferred to three-D grain object with And perform three continuous drawing process.

Array of vertices is transferred to three-D grain object by spreading moment t, GPU from PBO0, and performs three continuous draftings Process.Meanwhile, CPU binding PBO1 and performing calculate spreading moment t=t+deltaT range of scatter Rx, Ry, Rz and Minx, Generate array of vertices, array of vertices is loaded into a series of PBO1 processes.When spreading moment t=t+deltaT arrives, GPU will It is already loaded into the PBO1 finished and is transferred to three-D grain object, while CPU binds PBO0 again, when being that PBO0 transmits next diffusion The array of vertices at quarter.And so on, the continuous alternate transport array of vertices of PBO0, PBO1 to three-D grain object.Make full use of CPU With the characteristic of GPU asynchronous process, texture source is updated to PBO by CPU, and GPU is from other PBO copy textures and performs drawing process, Due to Asynchronous DMA transmission, renewal and drawing process can be performed simultaneously, improved dynamic and rendered performance.

In embodiment, set interval deltaT=10s.PBO0 summit will be already loaded into spreading moment t, GPU Array VolumeArray [M] [N] [L] is transferred to three-D grain object texture0 (size is M × N × L) and performs three companies Continuous drawing process.The characteristic of CPU and GPU asynchronous process can be made full use of using PBO, OpenGL can be in PBO0 and texture Asynchronous DMA transmission is performed between object texture0, so array of vertices is transferred into three-D grain object from PBO0 in GPU While three continuous drawing process of texture0 and execution, CPU binds PBO1 and performs calculating spreading moment t=t+ DeltaT range of scatter Rx, Ry, Rz and Minx, the three-dimensional array VolumeArray [M] [N] [L] for generating array of vertices, general Array of vertices VolumeArray [M] [N] [L] is loaded into a series of PBO1 processes.When spreading moment t=t+deltaT arrives, GPU will be already loaded into the PBO1 that finishes and be transferred to three-D grain object texture0, while CPU binds PBO0 again, it is PBO0 Transmit the array of vertices of next spreading moment.And so on, the continuous alternate transport array of vertices of PBO0, PBO1 to three-D grain pair As texture0, it is possible to achieve the not asynchronous generation and drafting of the array of vertices of Release and dispersion scope in the same time.Make full use of not With the asynchronous behavior between processor, texture source is updated to PBO by CPU, and GPU was drawn from other PBO copy textures and performing Journey, due to Asynchronous DMA transmission, renewal and drawing process can be performed simultaneously, improved dynamic and rendered performance.

When it is implemented, the deltaT that can voluntarily be set interval according to actual conditions.

Step 5, for spreading moment t, a cuboid big with three-D grain object etc. in step 4 is drawn first, is led to Texture mapping function is crossed, by three-D grain object map to cuboid, by GPGPU technologies, realizes and lets out in vertex shader Reveal the Coordinate Conversion of the array of vertices of range of scatter, the array of vertices concentration value of leakage range of scatter is realized in fragment shader General-purpose computations.Result of calculation is rendered and FBO is output to, then the concentration value volume data of diffusion space is used as new three-D grain It is written into texture cache；

In embodiment, for spreading moment t, the cuboid that a M × N × L size is drawn first is used as texture mapping Object, apex coordinate of the cuboid in voxel coordinate system be (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), corresponding to three-D grain coordinate (0,0,0), (1, 0,0)、(1,1,0)、(0,1,0)、(0,0,1)、(1,0,1)、(1,1,1)、(0,1,1).Three-D grain texture0 is mapped to Cuboid, the general-purpose algorithm of the calculating Release and dispersion concentration for a texture primitive is write using GLSL shading languages.

In vertex shader, input variable is apex coordinate gl_Vertex, the three-D grain coordinate gl_ of vertex correspondence TexCoord.Gl_Vertex is performed after model view transform, projective transformation, the viewport transform, screen coordinate gl_ is converted to position.Using XOY plane as projection plane when carrying out projective transformation, using the rectangular projection just to projection plane, depending on The viewport of M × N size is defined during mouth conversion, constant and picture is kept with the size and location for ensureing pixel in projection process Prime number evidence will not lose.The output variable of vertex shader is screen coordinate gl_position and texture coordinate gl_ TexCoord.Obtain the piece of M × N × L member after cuboid after the conversion of summit is rasterized, the corresponding three-dimensional pattern of each piece member Manage coordinate gl_TexCoord.

In fragment shader, input variable be texture coordinate gl_TexCoord, by texture coordinate gl_TexCoord (u, V, w) it is mapped to space coordinate (u × (M-1) × deltaX+Minx, v × (N-1) × deltaY-Ry, w in diffusion coordinate system × (L-1) × deltaZ), the corresponding concentration value c of piece member is obtained according to concentration of natural gas calculation formula C (x, y, z, t), and will Concentration value c is stored as gl_FragColor R channel values.The output variable of fragment shader is gl_FragColor.

Texture1 is associated with FBO, result of calculation is rendered and FBO is output to, then the concentration value body number of diffusion space Texture cache is written into according to as three-D grain texture1, then separates texture1 with FBO.

Step 6, for spreading moment t, a cuboid big with three-D grain object etc. in step 5 is drawn again, is led to Texture mapping function is crossed, by three-D grain object map to cuboid, by GPGPU technologies, realizes and lets out in vertex shader Reveal the Coordinate Conversion of the array of vertices of range of scatter, in fragment shader, the ray sent from summit to light source is sat in the world Space point sampling is carried out with certain step-length in Release and dispersion scope body bounding box after expanding in mark system, with reference to in step 5 The sampling of three-D grain, realizes the general-purpose computations of Release and dispersion volume data summit light transmittance.Result of calculation is rendered and defeated Go out to FBO, then the summit light transmittance of Release and dispersion volume data is written into texture cache as new three-D grain.

Present invention further propose that, in order to realize the visualization based on three-dimensional sphere, in fragment shader, summit is to light The ray that source is sent is to carry out sky in the Release and dispersion scope body bounding box after expanding in world coordinate system with certain step-length Between point sampling.Due to being sampled in the bounding box after expansion, it is therefore desirable to judge whether sampled point is located at leakage diffusion model In enclosing.Sampled point is transformed into WGS84 coordinate systems from world coordinate system：If the elevation of sampled point is less than the corresponding elevations of Rz and big In 0, then sample is in the range of Release and dispersion, otherwise, abandons the sampled point.If sample is in the range of Release and dispersion, Sampled point is then transformed into diffusion coordinate system, then voxel coordinate system is transformed into from diffusion coordinate system, then mapping obtains the sampling The three-D grain coordinate of point, samples at texture coordinate to three-D grain.

In the step 6 of embodiment, the cuboid of a M × N × L size is drawn again as the object of texture mapping, it is long Apex coordinate of the cube in voxel coordinate system is (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), corresponding to three-D grain coordinate (0,0,0), (1,0,0), (1,1,0)、(0,1,0)、(0,0,1)、(1,0,1)、(1,1,1)、(0,1,1).Three-D grain texture1 is mapped to rectangular Body.The general-purpose algorithm of the calculating per-vertex lighting transmissivity for a texture primitive is write using GLSL shading languages.

In vertex shader, input variable is apex coordinate gl_Vertex, the three-D grain coordinate gl_ of vertex correspondence TexCoord, light source world coordinates Light_worldWith the Release and dispersion scope body bounding box after expansion in world coordinate system.It is right Gl_Vertex is performed after model view transform, projective transformation, the viewport transform, is converted to screen coordinate gl_position.Enter Using XOY plane as projection plane during row projective transformation, using the rectangular projection just to projection plane, the definition during viewport transform is carried out The viewport of one M × N size, keeps constant and pixel data not lose with the size and location for ensureing pixel in projection process Lose.The output variable of vertex shader is world's seat of screen coordinate gl_position, texture coordinate gl_TexCoord, light source Mark Light_worldWith the Release and dispersion scope body bounding box after expansion in world coordinate system.Cuboid after the conversion of summit is through grating M × N × L piece member, one three-D grain coordinate gl_TexCoord of each piece member correspondence are obtained after change.

In fragment shader, input variable is texture coordinate gl_TexCoord, the world coordinates Light of light source_worldWith The body bounding box of Release and dispersion scope in world coordinate system, diffusion coordinate is mapped to by texture coordinate gl_TexCoord (u, v, w) Space coordinate (u × (M-1) × deltaX+Minx, v × (N-1) × deltaY-Ry, w × (L-1) × deltaZ) in system, then WGS84 coordinate systems are transformed into from diffusion coordinate system, from WGS84 coordinate systems are transformed into world coordinate system, obtain the texture primitive World coordinates Vertex_world, by Vertex_worldAnd Light_worldThe direction that can obtain light is：

Dir=Light_world-Vertex_world

Access point start and go out that Release and dispersion scope body bounding box after obtaining light and expanding in world coordinate system intersects Point end.The access point start of light is summit Vertex in this step_world.From Vertex_worldStart, in world coordinate system In sampled with distance s tep=0.1 along radiation direction dir, the sampled point obtained every time is：

Pos=Vertex_world+dir×(pos+step)

Sampled point pos is stored in array points successively, until going out point end.Array points is traveled through to accumulate meter Calculate the transmissivity on summit.Because the body bounding box of Release and dispersion scope in world coordinate system is the bounding box after expanding, therefore Need to judge whether sampled point is located in leakage range of scatter.Sampled point is transformed into WGS84 spherical coordinate systems from world coordinate system： If the elevation of sampled point is less than the corresponding elevations of Rz and more than 0, sample is mapped as in the range of Release and dispersion Three-D grain coordinate (u, v, w)；Otherwise, the sampled point is abandoned.Three-D grain texture1 is adopted at texture coordinate (u, v, w) place Sample obtains R channel values (i.e. concentration c), utilizes light transmittance calculation formula：

Transmittance=transmittance × exp (- ε × step × c) (ε is constant)

Accumulation obtains the transmissivity transmittance on summit, and transmittance is stored as into gl_FragColor R channel values.The output variable of fragment shader is gl_FragColor.

Texture2 is associated with FBO, result of calculation is rendered and FBO is output to, then the conduct of light transmittance volume data Three-D grain texture2 is written into texture cache, then separates texture2 with FBO；

When it is implemented, sampling interval step can voluntarily be set according to actual conditions.

Step 7, for spreading moment t, the rectangle of a screen size is finally drawn, by texture mapping function, will etc. Big 2 d texture object map is to rectangle, by GPGPU technologies, realizes that the coordinate on rectangle summit turns in vertex shader Change, the Release and dispersion model after being expanded in fragment shader from the ray that viewpoint is sent to screen pixels point in world coordinate system Space point sampling is carried out with certain step-length in containment body bounding box, it is real with reference to the sampling to three-D grain in step 5 and step 6 The calculating of existing screen pixels point color value, most rectangle is exported to screen at last, realizes t natural gas leaking range of scatter Three-dimensional visualization；

Present invention further propose that：In order to realize the visualization based on three-dimensional sphere, in fragment shader, from viewpoint to The ray that screen pixels point is sent be in world coordinate system expand after Release and dispersion scope body bounding box in certain step It is long to carry out spatial point sampling, due to being sampled in the bounding box after expansion, it is therefore desirable to judge whether sampled point is located at and let out Reveal in range of scatter.Sampled point is transformed into WGS84 spherical coordinate systems from world coordinate system：If the elevation of sampled point is less than Rz correspondences Elevation and more than 0, then sample is in the range of Release and dispersion, otherwise, abandons the sampled point.If sample expands in leakage Dissipate in scope, then sampled point is transformed into diffusion coordinate system, then voxel coordinate system is transformed into from diffusion coordinate system, then map To the three-D grain coordinate of the sampled point, three-D grain is sampled at texture coordinate.

In the step 7 of embodiment, finally drawing the rectangle of resolX × resolY size, (resolX × resolY is Screen resolution) as the object of texture mapping, four summits of rectangle are (0,0), (resolX-1,0), (resolX-1, ResolY-1), (0, resolY), corresponds to 2 d texture coordinate (0,0), (1,0), (1,1), (0,1) respectively.By 2 d texture Texture3 is mapped to rectangle.Using GLSL shading languages write for texture primitive calculating screen pixels color value it is logical Use algorithm.

In vertex shader, input variable is rectangle apex coordinate gl_Vertex, the 2 d texture coordinate of vertex correspondence Gl_TexCoord, viewpoint world coordinates Eye_worldWith the Release and dispersion scope body bounding box after expansion in world coordinate system. Gl_Vertex is converted to screen coordinate gl_position after model view transform, projective transformation, the viewport transform.Enter Using XOY plane as projection plane during row projective transformation, using the rectangular projection just to projection plane, the definition during viewport transform is carried out The viewport of one resolX × resolY size, constant and pixel is kept with the size and location for ensureing pixel in projection process Data will not lose.The output variable of vertex shader is screen coordinate gl_position, texture coordinate gl_TexCoord, regarded The world coordinates Eye of point_worldWith the Release and dispersion scope body bounding box after expansion in world coordinate system.Rectangle after the conversion of summit ResolX × resolY piece member, each piece member one screen coordinate gl_position of correspondence and texture are obtained after rasterized Coordinate gl_TexCoord.

In fragment shader, input variable is screen coordinate gl_position, texture coordinate gl_TexCoord, viewpoint World coordinates Eye_worldWith the Release and dispersion scope body bounding box after expansion in world coordinate system.By the world coordinates of viewpoint Eye_worldThe world coordinates of center's point of screen can be obtained, the world coordinates at screen coordinate gl_position is then obtained Screen_world.Pass through Eye_worldAnd Screen_worldThe direction that can determine that light is：

Dir=Screen_world‐Eye_world

If light with world coordinate system expand after Release and dispersion scope body bounding box it is non-intersect, abandon this light Line, the color value of screen pixels is background color value bcolor；Two intersection points start, end are obtained if intersecting, from start Start, sampled with distance s tep=0.1 along radiation direction dir, the sampled point obtained every time is：

Pos=start+dir × (pos+step)

Sampled point pos is stored in array points successively, until end.Traversal array points calculates screen to accumulate The color value of curtain pixel.Transparency opacity=1 is initialized, color value rgb calculation formula is：

Rgb=pcolor × lcolor

Pcolor is sampled point color value, is initialized as (1,1,1), and lcolor is lighting color value, be initialized as (1,1, 1).Because the body bounding box of Release and dispersion scope in world coordinate system is the bounding box after expanding, it is therefore desirable to judge sampling Whether point is located in leakage range of scatter.Sampled point is transformed into WGS84 spherical coordinate systems from world coordinate system：If the height of sampled point Journey is less than the corresponding elevations of Rz and more than 0, then sample is mapped as three-D grain coordinate in the range of Release and dispersion (u,v,w)；Otherwise, the sampled point is abandoned.R is obtained to three-D grain texture1 samplings respectively at texture coordinate (u, v, w) place Channel value (i.e. concentration c), to three-D grain texture2 sampling obtain R channel values (i.e. illumination transmissivity transmittance), Utilize：

T=exp (- ε × step × c) (ε is constant)

Opacity=opacity × t

Accumulation obtains the transparency opacity for the light that screen pixels are sent, until opacity is close to 0 or reaches end.Utilize：

Accumulation obtains the color value color on light, utilizes：

Pixel=color+opacity × bcolor

Obtain the color value pixel of screen pixels.FBO is stopped using, then final texture3 is rendered into frame buffer, realized Volume drawing.

When it is implemented, sampled point color value plocor and lighting color value can voluntarily be set according to actual conditions lcolor。

In step 5, step 6 and step 7, by continuous three drawing process realize diffusion concentration data generation and can Depending on the integration of change.

Embodiment repeats said process with time interval deltaT, realizes the dynamic visual of natural gas leaking range of scatter Change.At the end of natural gas leaking dispersion ability visual Simulation, PBO0, PBO1 and FBO are deleted, texture object is deleted texture0、texture1、texture2、texture3。

Specific implementation is as can be seen that realize that concentration value volume data and per-vertex lighting are saturating by GPGPU technologies more than The general-purpose computations of rate volume data are penetrated, the calculating of screen pixels point color value are then completed by sampled three-dimensional texture, although in wash with watercolours Need have certain increase on self-defined coloring process and algorithm complex during dye, but used PBO, multi-scale sampling etc. Optimisation strategy and FBO technologies, and consider the influence of three-dimensional sphere, can realize spread under natural gas leaking diffusion conditions it is dense Degrees of data generate and dynamic and visual based on sphere integration and illustrate concentration point in natural gas leaking diffusion process The overall picture and details of cloth.

Claims

1. a kind of dynamic and visual method of the natural gas leaking diffusion based on GPGPU, it is characterised in that：Including procedure below,

First, environment is drawn in initialization, and two PBO and FBO are created using OpenGL, creates three-D grain object and two dimension Texture object is used for the result of calculation for storing GPGPU, definition diffusion coordinate system and voxel coordinate system, and initialization viewpoint, light source exist Position in WGS84 coordinate systems；

Step a1, by CPU diffusion coordinate system in, calculate spreading moment t when natural gas respectively x-axis, y-axis, z-axis diffusion away from From Rx, Ry, Rz and range of scatter x-axis minimum M inx；

Step a2 is related to the LOD ranks of three-dimensional scenic by the sampled distance in range of scatter by CPU in diffusion coordinate system Connection, with deltaX, deltaY, deltaZ distance respectively in x-axis, y-axis, z-axis direction under the LOD ranks of current three-dimensional scenic Sampling, generation t reveals the array of vertices of range of scatter, and readjust Rx, Ry, Rz so that x-axis, y-axis, z-axis direction Sampled point number is 2 whole power；The body bounding box for spreading Release and dispersion scope in coordinate system is transformed into WGS84 coordinates System, is then transformed into world coordinate system, obtains the body bounding box of Release and dispersion scope in world coordinate system；It is parallel flat by two The body bounding box of face and Release and dispersion scope in world coordinate system builds the new body bounding box after expanding, one of plane By the corresponding points of ((Rx+Minx)/2,0, Rz) in world coordinate system, another plane, which passes through to leak in world coordinate system, to expand Dissipate four summits of scope body bounding box bottom surface；

Step a3, is loaded into an idle PBO by the array of vertices that t reveals range of scatter by CPU, makes t=t+deltaT, DeltaT is time interval, for next current spreading moment return to step a1, step a1-a3 is performed, by new diffusion The array of vertices of moment t leakage range of scatters is loaded into another PBO；

After the array of vertices that t reveals range of scatter is loaded into PBO by CPU, step a1-a3 is performed for next spreading moment When, array of vertices is transferred to being drawn below array of vertices size identical three-D grain object, and completion from PBO by GPU Processing, draw two cuboids, by concentration value volume data and per-vertex lighting transmissivity volume data by way of texture mapping Realize that GPGPU is calculated in fragment shader, result of calculation is rendered and FBO is output to, then concentration value volume data and per-vertex lighting Transmissivity volume data is written into texture cache respectively as three-D grain；A rectangle is drawn, screen point color value data is led to The mode of texture mapping is crossed, is adopted with reference to concentration value volume data three-D grain and per-vertex lighting transmissivity volume data three-D grain Sample, realizes that GPGPU is calculated in fragment shader, and result of calculation write-in 2 d texture is simultaneously directly rendered into screen, realizes expansion Dissipate concentration data generation and visual integration.

2. a kind of dynamic and visual method that natural gas leaking based on GPGPU spreads as claimed in claim 1, its feature exists In：In step 1, it is defined as follows for diffusion coordinate system and voxel coordinate system,

Diffusion coordinate system is defined as with source of leaks p₀Subpoint on ground is origin, and following wind direction is x-axis, with ground Perpendicular to x-axis direction be y-axis, using perpendicular to ground straight up direction as z-axis；In diffusion coordinate system, range of scatter is by growing Cube bounding box (Minx ,-Ry, 0), (Rx ,-Ry, 0), (Rx, Ry, 0), (Minx, Ry, 0), (Minx ,-Ry, Rz), (Rx ,- Ry, Rz), (Rx, Ry, Rz), (Minx, Ry, Rz) determine that Rx, Ry, Rz are respectively that x-axis, y-axis, z-axis maximally diffuse distance, Minx is minimum value of the range of scatter in x-axis；

Voxel coordinate system is defined as to spread in coordinate system (Minx ,-Ry, 0) for origin, x-axis, y-axis, z-axis direction with diffusion Coordinate system is identical；In voxel coordinate system, range of scatter by cuboid bounding box (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) are determined, M-1, N-1, L-1 Respectively x-axis, y-axis, z-axis maximally diffuse distance, and M, N, L represent the sampled point number of x-axis, y-axis, z-axis respectively.

3. a kind of dynamic and visual method that natural gas leaking based on GPGPU spreads as claimed in claim 2, its feature exists In：In step a1, x-axis, y-axis, diffusion length Rx, Ry, Rz, Minx of z-axis computational methods are as follows in diffusion coordinate system,

(1) calculate Minx method be,

As t ＜ T, Minx=0 is made, Minx values are preserved；

As t >=T,

If 2. c ＜ c0, make Minx=Minx+delta, delta recalculate the concentration at (Minx, 0, H) place to preset step-length Value c, until c >=c0, into 3.；If c >=c0, it is directly entered 3.；

If 3. c >=c0, makes Minx=Minx-1, the concentration value c at (Minx, 0, H) place is recalculated, until c ＜ c0, is preserved Minx values；

(2) calculate Rx method be,

If 2. c >=c0, makes Rx=Rx+delta, delta recalculate the concentration value c at (Rx, 0, H) place, directly to preset step-length To c ＜ c0, into 3.；If 3. c ＜ c0, be directly entered；

If 3. c ＜ c0, make Rx=Rx-1, the concentration value c at (Rx, 0, H) place is recalculated, until c >=c0, Rx values are preserved；

(3) calculate Ry method be,

1. Ry=0, Temp_x=Rx, center_x=0 are initialized；Temp_x is the variable of interim storage x values, and center_x is Corresponding x values at Ry；

2. the concentration value c at (Temp_x, Ry, H) place is calculated, if c >=c0, Ry=Ry+1 is made, recalculates (Temp_x, Ry, H) The concentration value c at place, until c ＜ c0, into 3.；If 3. c ＜ c0, be directly entered；

(4) calculate Rz method be,

1. Rz=H, Temp_x=center_x, x_left=Temp_x, x_right=Temp_x are initialized；Temp_x、x_ Left, x_right are the variable of interim storage x values, and center_x is corresponding x values at Ry；

2. the concentration value c at (Temp_x, 0, Rz) place is calculated, if c >=c0, Rz=Rz+1 is made, recalculates (Temp_x, 0, Rz) The concentration value c at place, until c ＜ c0, into 3.；If 3. c ＜ c0, be directly entered；

If c1 ＜ c0 and c2 ＜ c0, Rz values are preserved；

If c1 >=c0 and c2 ＜ c0, Temp_x=Temp_x-1, return 2., repeat 2. -3.；

If c1 ＜ c0 and c2 >=c0, Temp_x=Temp_x+1, return 2., repeat 2. -3..

4. a kind of dynamic and visual method of the natural gas leaking diffusion based on GPGPU, its feature as described in Claims 2 or 3 It is：In step a2, by reveal range of scatter visualization precision be divided into several ranks, by the sampling in range of scatter away from It is associated from the LOD ranks of three-dimensional scenic, optimal M, N, L value is obtained by the strategy of multi-scale sampling, optimization texture reflects Penetrate efficiency.

5. a kind of dynamic and visual method of the natural gas leaking diffusion based on GPGPU, its feature as described in Claims 2 or 3 It is：If two PBO are designated as PBO0 and PBO1 respectively,

If CPU binds PBO0, the array of vertices generation processing of leakage range of scatter is carried out, spreading moment t array of vertices is carried Enter PBO0, array of vertices is transferred to three-D grain object from PBO0 and drawn by GPU；

While array of vertices is transferred to three-D grain object from PBO0 and drawn by GPU, CPU binding PBO1 are let out Reveal the array of vertices generation processing of range of scatter, next spreading moment t array of vertices is loaded into PBO1；

Array of vertices is transferred to three-D grain object from PBO1 and drawn by GPU, while CPU binds PBO0 again, so Analogize, the continuous alternate transport array of vertices of PBO0, PBO1 to three-D grain object.

6. a kind of dynamic and visual method of the natural gas leaking diffusion based on GPGPU, its feature as described in Claims 2 or 3 It is：GPU drawing modification comprises the following steps,

Step b1, draws a cuboid big with three-D grain object etc., by texture mapping function, by three-D grain first Object map, by GPGPU technologies, realizes the coordinate of the array of vertices of leakage range of scatter to cuboid in vertex shader Conversion, realizes the general-purpose computations of the array of vertices concentration value of leakage range of scatter in fragment shader；Result of calculation is rendered And being output to FBO, then the concentration value volume data of diffusion space is written into texture cache as new three-D grain；

Step b2, draws a cuboid big with three-D grain object etc., by texture mapping function, by three-D grain again Object map, by GPGPU technologies, realizes the coordinate of the array of vertices of leakage range of scatter to cuboid in vertex shader Conversion, in fragment shader, the ray sent from summit to light source expand in world coordinate system after Release and dispersion scope Space point sampling is carried out by step-length in body bounding box, with reference to the sampling to three-D grain in step b2, Release and dispersion body number is realized According to the general-purpose computations of summit light transmittance, result of calculation is rendered and FBO is output to, then the summit light of Release and dispersion volume data Line transmissivity is written into texture cache as new three-D grain；

Step b3, finally draws the rectangle of a screen size, by texture mapping function, by etc. big 2 d texture object reflect Be mapped to rectangle, by GPGPU technologies, the Coordinate Conversion on rectangle summit realized in vertex shader, in fragment shader by The ray that viewpoint is sent to screen pixels point expand in world coordinate system after Release and dispersion scope body bounding box in by step-length Space point sampling is carried out, with reference to the sampling to three-D grain in step b1 and step b2, the meter of screen pixels point color value is realized Calculate, most rectangle is exported to screen at last, realizes the three-dimensional visualization of t natural gas leaking range of scatter.

7. a kind of dynamic and visual method that natural gas leaking based on GPGPU spreads as claimed in claim 6, its feature exists In：In step b2, the sampling in the bounding box after expansion is revealed in range of scatter, it is necessary to judge whether sampled point is located at, bag Include sampled point and to be transformed into WGS84 coordinate systems from world coordinate system, if the elevation of sampled point is less than the corresponding elevations of Rz and is more than 0, then sample the sampled point is abandoned in the range of Release and dispersion, otherwise；, will if sample is in the range of Release and dispersion Sampled point is transformed into diffusion coordinate system, then is transformed into voxel coordinate system from diffusion coordinate system, and then mapping obtains the sampled point Three-D grain coordinate, samples at texture coordinate to three-D grain.

8. a kind of dynamic and visual method that natural gas leaking based on GPGPU spreads as claimed in claim 6, its feature exists In：In step b3, the sampling in the bounding box after expansion is revealed in range of scatter, it is necessary to judge whether sampled point is located at, bag Include sampled point and to be transformed into WGS84 spherical coordinate systems from world coordinate system, if the elevation of sampled point is less than the corresponding elevations of Rz and big In 0, then sample abandons the sampled point in the range of Release and dispersion, otherwise；If sample is in the range of Release and dispersion, Sampled point is transformed into diffusion coordinate system, then voxel coordinate system is transformed into from diffusion coordinate system, then mapping obtains the sampled point Three-D grain coordinate, at texture coordinate to three-D grain sample.