CN107093207B

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

Info

Publication number: CN107093207B
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: 2019-07-09
Anticipated expiration: 2037-04-12
Also published as: CN107093207A

Abstract

The present invention provides a kind of dynamic and visual method of 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 of t moment leakage range of scatter is loaded into PBO by CPU, when executing array of vertices generation processing for next spreading moment and be 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 spherical surface under natural gas leaking diffusion conditions.The present invention improves dynamic rendering performance using multi-scale sampling policy optimization texture mapping efficiency, using the characteristic of CPU and GPU asynchronous process, 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

The present invention relates to GPGPU and three-dimensional geographic information visualization technique field, are directed to natural gas more particularly to one kind Diffusion concentration data under the conditions of Release and dispersion generate and the integrated method of dynamic and visual based on spherical surface.

Background technique

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 range 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 related document proposes 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 diffusion Test tube reason and technology, 2008,05:17-19+55；Air pollution diffusion simulation system of the Ou Yangkun based on three peacekeeping Temporal GISs Research and realization [D] Tsinghua University, 2011；Block toxic gas diffusion process model building under Zheng Maohui, Jin Min, Guo Fei .GIS support With simulation [J] Wuhan University Journal (information science version), 2013,08:935-939；Simple flood is stepped on, Fan Xiangtao, Wang Jinxin subway Website harmful gases diffusion three-dimensional visualization studies [J] Surveying and mapping, 2013,01:136-138+141.In combination GIS and greatly Gas diffusion model is come when expressing natural gas leaking range of scatter, many researchs lay particular emphasis on Release and dispersion range on two-dimensional map Visualization；Release and dispersion range is mainly based upon Marching Cubes algorithm in the diffusion simulations of three-dimensional space and draws concentration etc. The three-dimensional visualization in value face, has problems in that: (1) concentration contour surface can only indicate in data the structure of part contour surface and Spatial relationship cannot 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 building of concentration contour surface affect rendering efficiency, therefore often by concentration equivalent point The generation of data is as preprocessing process, it is difficult to realize that concentration data generates and visual integration.

Volume Rendering Techniques are one kind of scientific visualization technology, can generate the general image of three-dimensional data, including body Interior useful details.Therefore, Volume Rendering Techniques are able to solve above-mentioned first problem.Related document combination Volume Rendering Techniques are realized The visualization of 3 d data field.The research and application that Hou Meihua, Wu Zhihong, Chen Kaimin are simulated based on the Real-time Smoke of OSG [J] computer engineering and design, 2011,06:2088-2091；Recklessly from and, Liu Po, Gong Jianhua, Wang Qun is based on virtual earth The design and realization [J] Wuhan University Journal (information science version) of typhoon Dynamic and Multi dimensional visualization system, 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 Volume Rendering Techniques are applied to natural gas leakage diffusion simulations not yet at present.

GPGPU be using GPU carry out general-purpose computations technology, classical GPGPU technology by GPU graphics pipeline big rule Mould computation capability is realized the parallel computation of general-purpose algorithm by texture mapping function, has significantly speeded up general-purpose computations mistake Journey.The combination of GPGPU technology and Volume Rendering Techniques is able to solve above-mentioned Second Problem, and the generation of diffusion concentration data may be implemented With the integration of dynamic and visual.Nuclear proliferation in Bian Yanshan three-dimensional digital battlefield shows technology [D] science and techniques of defence University, 2010.GPGPU technology and Volume Rendering Techniques is combined to realize the dynamic diffusion of three-dimensional digital battlefield Radionuclide in document Process.The disadvantage of this method is that nucleic diffusion volume data is calculated in GPGPU and based on nucleic diffusion volume data reality Existing volume drawing division is not carried out the integration of volume data generation and dynamic and visual for two processes.It is counted using GPGPU It calculates and obtains not accounting for the influence of three-dimensional sphere when nucleic diffusion volume data, and use fixed three-dimensional diffusion range, it is difficult to Completely show the overall picture of different moments range of scatter in dynamic process.

Relational language:

GPGPU general-purpose computations graphics processor

GIS GIS-Geographic Information System

GPU graphics processor

CPU central processing unit

OpenGL open graphic library is the 2D/3D figure API received the most extensively in industry field

The pixel buffer PBO object

The frame buffer zone FBO object

Summary of the invention

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

Technical solution of the present invention provides a kind of dynamic and visual method of natural gas leaking diffusion based on GPGPU, packet Following procedure is included,

Firstly, initialization draw environment, using OpenGL create two PBO and FBO, creation three-D grain object and 2 d texture object is used to store the calculated result of GPGPU, definition diffusion coordinate system and voxel coordinate system, initializes viewpoint, light Position of the source in WGS84 coordinate system；

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 in diffusion coordinate system, natural gas is respectively in the expansion of x-axis, y-axis, z-axis when calculating spreading moment t Distance Rx, Ry, Rz and range of scatter are dissipated in the minimum M inx of x-axis；

Step a2, by CPU in diffusion coordinate system, by the LOD rank of sampled distance and three-dimensional scenic in range of scatter It is associated, with the distance of deltaX, deltaY, deltaZ respectively in x-axis, y-axis, z-axis under the LOD rank of current three-dimensional scenic Direction sampling, generates the array of vertices of t moment leakage range of scatter, and readjusts Rx, Ry, Rz, so that x-axis, y-axis, z-axis side To number of sampling points be 2 whole power；The body bounding box for spreading Release and dispersion range 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 range in world coordinate system；It is parallel by two New body bounding box after the body bounding box building of Release and dispersion range expands in plane and world coordinate system, one of them is flat Corresponding points of the face by ((Rx+Minx)/2,0, Rz) in world coordinate system, another plane in world coordinate system by leaking Four vertex of range of scatter body bounding box bottom surface；

The array of vertices of t moment leakage range of scatter is loaded into an idle PBO by CPU, enables t=t+ by step a3 DeltaT, deltaT are time interval, for next current spreading moment return step a1, execute 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 of t moment leakage range of scatter is loaded into PBO by CPU, step is executed for next spreading moment When a1-a3, array of vertices is transferred to three-D grain object identical with array of vertices size from PBO by GPU, 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, and renders and be output to FBO calculated result, then concentration value volume data and vertex light in fragment shader Texture cache is written into respectively as three-D grain according to transmissivity volume data；A rectangle is drawn, by screen point color value data By way of texture mapping, in conjunction with 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 calculated result write-in 2 d texture is simultaneously directly rendered into screen, realizes Diffusion concentration data generate and visual integration.

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

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 be defined as with spread in coordinate system (Minx ,-Ry, 0) be origin, x-axis, y-axis, z-axis direction with It is identical to spread coordinate system；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 respectively indicate the number of sampling points of x-axis, y-axis, z-axis.

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

(1) calculate Minx method be,

As t < T, Minx=0 is enabled, saves Minx value；

As t >=T,

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

2. enabling Minx=Minx+delta, delta is preset step-length if c < c0, recalculate at (Minx, 0, H) Concentration value c, until c >=c0, into 3.；If 3. c >=c0 is directly entered；

3. enabling Minx=Minx-1 if c >=c0, the concentration value c at (Minx, 0, H) is recalculated, until c < c0, is protected Deposit Minx value；

(2) calculate Rx method be,

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

2. enabling Rx=Rx+delta, delta is preset step-length if c >=c0, the concentration value at (Rx, 0, H) is recalculated C, until c < c0, into 3.；If 3. c < c0, is directly entered；

3. recalculating the concentration value c at (Rx, 0, H) if c < c0, enables Rx=Rx-1, until c >=c0, Rx is saved Value；(3) calculate Ry method be,

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

2. calculating the concentration value c at (Temp_x, Ry, H) enables Ry=Ry+1 if c >=c0, recalculate (Temp_x, Ry, H) at concentration value c, until c < c0, into 3.；If 3. c < c0, is directly entered；

3. enabling center_x=Temp_x, Temp_x=Temp_x-1, the concentration value c at (Temp_x, Ry, H) is calculated, if C < c0, then save Ry value；If c >=c0, return 2., repeat 2. -3.；

(4) calculate Rz method be,

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

2. calculating the concentration value c at (Temp_x, 0, Rz) enables Rz=Rz+1 if c >=c0, recalculate (Temp_x, 0, Rz) the concentration value c at, until c < c0, into 3.；If 3. c < c0, is directly entered；

3. enabling x_left=Temp_x-1, x_right=Temp_x+1, concentration c 1, (x_ at judgement (x_left, 0, Rz) Right, 0, Rz) concentration c 2 and given concentration value c at₀Size,

If c1 < c0 and c2 < c0, save Rz value；

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, the visualization precision for revealing range of scatter is divided into several ranks in step a2, it will be in range of scatter Sampled distance it is associated with the LOD rank 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 denoted as PBO0 and PBO1 respectively,

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

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

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.

Moreover, the drawing modification of GPU includes the following steps,

Step b1, the big cuboid of one and three-D grain object etc. of drafting first will be three-dimensional by texture mapping function Texture object is mapped to cuboid, and by GPGPU technology, the array of vertices of leakage range of scatter is realized in vertex shader The general-purpose computations of the array of vertices concentration value of leakage range of scatter are realized in coordinate conversion in fragment shader；By calculated result FBO is rendered and is output to, then the concentration value volume data of diffusion space is written into texture cache as new three-D grain；

Step b2, drawing the big cuboid such as one and three-D grain object again will be three-dimensional by texture mapping function Texture object is mapped to cuboid, and by GPGPU technology, the array of vertices of leakage range of scatter is realized in vertex shader Coordinate conversion, in fragment shader, the ray issued from vertex to light source expand in world coordinate system after Release and dispersion Space point sampling is carried out by step-length in range body bounding box, in conjunction with the sampling to three-D grain in step b2, realizes Release and dispersion Calculated result, is rendered and is output to FBO by the general-purpose computations of volume data vertex light transmittance, 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 technology, the coordinate conversion on rectangle vertex is realized in vertex shader, in fragment shader In expand in world coordinate system from viewpoint to the ray that screen pixels point issues after Release and dispersion range body bounding box in press Step-length carries out space point sampling, in conjunction with the sampling to three-D grain in step b1 and step b2, realizes screen pixels point color value Calculating, finally rectangle is exported to screen, realizes the three-dimensional visualization of t moment natural gas leaking range of scatter.

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

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

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

Detailed description of the invention

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 flow chart of the calculating Rx and Minx of the embodiment of the present invention.

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

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

Fig. 6 is the Release and dispersion range 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.

Specific embodiment

Below in conjunction with drawings and examples the present invention will be described in detail technical solution.

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, judgment basis be critical pressure ratio CPR:

At that time, the flow velocity of aperture leakage can reach velocity of sound, leak intensity are as follows:

WhenWhen, the flow velocity of aperture leakage is subsonic speed, reveals intensity are as follows:

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

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

In formula, C (x, y, z, t)-spreads the mean concentration of gas leakage at certain point (x, y, z) in coordinate system, kg/m³； 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；T-natural gas leakage total time, s；H-effective source height, m.

Embodiment uses method of the invention, using OpenGL and GLSL shading language, is based on osgEarth three-dimensional rendering Spherical surface is generated and is based on based on the diffusion concentration data of GPGPU and volume drawing under engine implementation natural gas leaking diffusion conditions 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 Following procedure,

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

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

For the sake of the implementation that is convenient for reference, provide embodiment detailed process design it is as follows, referring to Fig. 7:

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

Initialization step can be executed by CPU.

The present invention defines 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 be defined as with spread in coordinate system (Minx ,-Ry, 0) be origin, x-axis, y-axis, z-axis direction with It is identical to spread coordinate system, 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) are determined, M-1, N-1, L-1 is respectively that x-axis, y-axis, z-axis maximally diffuse distance, and M, N, L respectively indicate the number of sampling points of x-axis, y-axis, z-axis.

Corresponding diffusion coordinate system midpoint (i × deltaX+Minx, j × deltaY-Ry, k of voxel coordinate system midpoint (i, j, 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 Axis, y-axis, the maximum distance of z-axis, M, N, L respectively indicate the number of sampling points of x-axis, y-axis, z-axis；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.It creates a 2 d texture object texture3 and is used as rendering screen.Creation one FBO, as being rendered into texture.

The position for initializing source of leaks is p_source(lon, lat, h), the initial position of viewpoint are p_eye(lon, lat, h), The position of light source is p_lig(lon, lat, h), p_source、p_eye、p_ligIt is WGS84 coordinate.By viewpoint position p_eyeAnd light source position p_ligIt is transformed into world coordinate system from WGS84 coordinate system, then the world coordinates Eye of viewpoint can be obtained_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 determine source of leaks, viewpoint and light source in WGS84 coordinate system according to the actual situation 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 proposes:

For Minx, if spreading moment t is less than natural gas leakage total time T, Minx is enabled to be equal to 0；Otherwise it is assigned for Minx One smaller value (those skilled in the art can preset value) calculates concentration according to concentration of natural gas calculation formula C (x, y, z, t) C makes the concentration c at t moment Minx be greater 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 larger value (those skilled in the art can preset value) is assigned for Rx first, according to concentration of natural gas Calculation formula C (x, y, z, t) calculates concentration c, so that the concentration c at t moment Rx is less than prescribed concentration c0, is then gradually reduced Rx, Until corresponding concentration c is greater than or equal to prescribed concentration c0；

For Ry, start at (Rx, 0, H), be gradually increased y value, 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 saved at this time is Ry.If after reducing x, the concentration c at (x, Ry, H) Less than prescribed concentration c0, then Ry is maximum value, otherwise continues to increase y value solution Ry；

Rz is started at (x, 0, H) after finding out the corresponding x value of Ry, is gradually increased z value, according to concentration of natural gas meter It calculates formula C (x, y, z, t) and calculates concentration c, until concentration c is less than prescribed concentration c0, the z saved at this time is Rz.If reducing x or increasing When big x, the concentration c at (x, 0, Rz) is respectively less than prescribed concentration c0, then Rz is maximum value, otherwise, dense greater than specifying along concentration c The side of degree c0 repeats aforesaid operations and continues to solve 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 are as follows:

As t < T: enabling Minx=0, save Minx value；

As t >=T:

1. initializing Minx=x0, the concentration value c at (Minx, 0, H) is calculated, those skilled in the art can preset initial value X0 value, the preferably value range of x0 are 0≤x0≤10；

2. if c < c0, enable Minx=Minx+delta (those skilled in the art can preset step-length delta value, such as Desirable 50 or 100 or 200), the concentration value c at (Minx, 0, H) is recalculated, until c >=c0, into 3.；If c >=c0, directly It taps into 3.；

3. enabling Minx=Minx-1 if c >=c0, the concentration value c at (Minx, 0, H) is recalculated, until c < c0, is protected Deposit Minx value.

Such as Fig. 3, the method for calculating Rx are as follows:

1. initialize Rx=Minx+x0 (those skilled in the art can preset initial value x0 value, the preferred value range of x0 For 0x0≤10), calculate the concentration value c at (Rx, 0, H)；

2. if c >=c0, enable Rx=Rx+delta (those skilled in the art can preset step-length delta value, such as it is desirable 50 or 100 or 200), the concentration value c at (Rx, 0, H) is recalculated, until c < c0, into 3.；If c < c0, is directly entered ③；

3. recalculating the concentration value c at (Rx, 0, H) if c < c0, enables Rx=Rx-1, until c >=c0, Rx is saved Value.

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

Such as Fig. 4, the method for calculating Ry are as follows:

Such as Fig. 5, the method for calculating Rz are as follows:

3. enabling x_left=Temp_x-1, x_right=Temp_x+1, concentration c 1, (x_ at judgement (x_left, 0, Rz) Right, 0, Rz) concentration c 2 and given concentration value c at₀Size:

If c1 < c0 and c2 < c0, save Rz value；

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 be determined voluntarily according to the actual situation Concentration c₀。

It is step 3, in diffusion coordinate system, the sampled distance in range of scatter is associated with the LOD rank of three-dimensional scenic, It is adopted respectively in x-axis, y-axis, z-axis direction under the LOD rank of current three-dimensional scenic with the distance of deltaX, deltaY, deltaZ Sample generates the array of vertices of t moment leakage range of scatter, and readjusts 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 range in coordinate system is transformed into WGS84 coordinate system, Then it is transformed into world coordinate system, obtains the body bounding box of Release and dispersion range in world coordinate system (such as Fig. 6 is shown in solid). Body bounding box at this time is surrounded by two spherical surfaces and four planes, is expanded by leaking in two parallel planes and world coordinate system The body bounding box of scattered range can construct the new body bounding box after expanding, as shown in Fig. 6 dotted line.

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 rank of three-dimensional scenic, passes through the policy optimization of multi-scale sampling texture mapping effect Rate:

The LOD rank 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 rank of scene samples in leakage range of scatter, generates the array of vertices of Release and dispersion range.It enables Index K=log₂M (K round numbers), i.e. M=2^K, similarly to N, L again assignment, making M, N, L is 2 whole power, and according to M, N, L readjusts Rx, Ry, Rz.The strategy of multi-scale sampling is meeting the same of Release and dispersion range three-dimensional visualization visual acuity When, optimal M, N, L value can be obtained, due to the size of three-D grain when M, N, L have codetermined texture mapping, is passed through The strategy of multi-scale sampling can optimize texture mapping efficiency.

In step 3, the sample space point within the scope of three-dimensional diffusion changing over time, accurately calculating can be with minimum Data volume completely show the overall picture of different moments range of scatter in dynamic process:

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

In embodiment, the visualization precision for revealing range of scatter is divided into 5 ranks, the LOD rank with three-dimensional scenic It is associated, it is indicated 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 rank of current three-dimensional scenic is LOD5 or is greater than LOD5, deltaX=deltaY=deltaZ=1m；

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

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

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

When the LOD rank of current three-dimensional scenic is LOD1 or is less than LOD1, deltaX=deltaY=deltaZ=16m. With the distance of deltaX, deltaY, deltaZ respectively in x, y, z direction sample space coordinate, array of vertices is generated.Enable 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, enable Rx=(M-1) × deltaX+Minx, Ry=(N-1) × deltaY/2, Rz=(L-1) × deltaZ.X-axis, y-axis, the number of sampling points 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 vertex are diffusion coordinate (i × deltaX+Minx, j × deltaY-Ry, k × deltaZ), 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 vertex (Minx ,-Ry, 0) of body bounding box that range of scatter is revealed 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 range body bounding box vertex p1, p2,p3,p4,p5,p6,p7,p8.The coordinate points ((Rx+Minx)/2,0, Rz) spread in coordinate system are transformed into WGS84 coordinate In system, reconvert obtains its spatial position pc in world coordinate system into world coordinate system.As shown in Figure 6 by two Parallel plane (one of plane passes through pc, another plane passes through p1, p2, p3, p4) and Release and dispersion in world coordinate system The body bounding box of range can construct the new body bounding box after expanding.When it is implemented, can be voluntarily true according to the actual situation Determine Release and dispersion range visualization precision grade and sampled distance deltaX, deltaY under different accuracy rank, deltaZ。

Step 4, array of vertices is transferred to three-dimensional identical with array of vertices size from PBO for spreading moment t, GPU Texture object simultaneously executes continuous drawing process three times, while CPU will realize the vertex battle array spreading moment t=t+deltaT It lists into another PBO.Using the characteristic of CPU and GPU asynchronous process, the leakage of spreading moment t and subsequent time are spread The array of vertices of range is alternately performed the update in texture source and the duplication of texture using two PBO, realizes that different moments leakage is expanded Dissipate the asynchronous generation and drafting of the array of vertices of range.

The present invention proposes, the asynchronous generation of the array of vertices of different moments Release and dispersion range is realized using two PBO And drafting:

CPU is responsible for executing range of scatter Rx, Ry, Rz and Minx for calculating spreading moment t=t+deltaT, generates vertex battle array Column, 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 execute continuous drawing process three times.

Array of vertices is transferred to three-D grain object from PBO0 by spreading moment t, GPU, and executes continuous three times draw Process.Meanwhile CPU binding PBO1 and execute calculate spreading moment t=t+deltaT range of scatter Rx, Ry, Rz and Minx, It generates array of vertices, array of vertices is loaded into a series of processes of PBO1.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 transmitting next diffusion for PBO0 The array of vertices at quarter.So analogize, 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 executes drawing process, It since Asynchronous DMA transmits, updates and drawing process may be performed simultaneously, improve dynamic rendering performance.

In embodiment, set interval deltaT=10s.The vertex of PBO0 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 executes to be connected three times 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 executed between object texture0, so array of vertices is transferred to three-D grain object from PBO0 in GPU Texture0 and execute three times continuous drawing process while, CPU binding PBO1 and execute calculate spreading moment t=t+ Range of scatter Rx, Ry, Rz and Minx of deltaT, it the three-dimensional array VolumeArray [M] [N] [L] for generating array of vertices, incites somebody to action A series of processes of array of vertices VolumeArray [M] [N] [L] loading PBO1.When spreading moment t=t+deltaT arrives, GPU will be already loaded into the PBO1 finished and be transferred to three-D grain object texture0, while CPU binds PBO0 again, be PBO0 Transmit the array of vertices of next spreading moment.So analogize, the continuous alternate transport array of vertices of PBO0, PBO1 to three-D grain pair As texture0, the asynchronous generation and drafting of the array of vertices of different moments Release and dispersion range may be implemented.It makes 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 executing Journey updates since Asynchronous DMA transmits and drawing process may be performed simultaneously, and improves dynamic rendering performance.

When it is implemented, can self-setting time interval deltaT according to the actual situation.

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, three-D grain object is mapped to cuboid, by GPGPU technology, realizes and lets out in vertex shader The coordinate conversion for revealing the array of vertices of range of scatter, realizes the array of vertices concentration value of leakage range of scatter in fragment shader General-purpose computations.Calculated result is rendered to and is output to FBO, then the concentration value volume data of diffusion space is as new three-D grain It is written into texture cache；

In embodiment, for spreading moment t, drafting one M × N × L size cuboid first is 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), correspond 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 writes the general-purpose algorithm of the calculating Release and dispersion concentration for a texture primitive using GLSL shading language.

In vertex shader, input variable is the three-D grain coordinate gl_ of apex coordinate gl_Vertex, vertex correspondence TexCoord.After executing model view transform, projective transformation, the viewport transform to gl_Vertex, it is converted to screen coordinate gl_ position.Using XOY plane as projection plane when carrying out projective transformation, using the rectangular projection of face projection plane, depending on The viewport of M × N size is defined when mouthful transformation, to guarantee that the size and location of pixel in projection process remains unchanged and picture Prime number evidence will not lose.The output variable of vertex shader is screen coordinate gl_position and texture coordinate gl_ TexCoord.M × N × L piece member, the corresponding three-dimensional pattern of each member are obtained after the transformed cuboid in vertex is rasterized 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 the R channel value of gl_FragColor.The output variable of fragment shader is gl_FragColor.

Texture1 is associated with FBO, calculated result is rendered and be output to FBO, then the concentration value body number of diffusion space It is written into texture cache 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, three-D grain object is mapped to cuboid, by GPGPU technology, realizes and lets out in vertex shader The coordinate conversion for revealing the array of vertices of range of scatter, in fragment shader, the ray issued from vertex to light source is sat in the world Space point sampling is carried out with certain step-length in Release and dispersion range body bounding box after expanding in mark system, in conjunction in step 5 The general-purpose computations of Release and dispersion volume data vertex light transmittance are realized in the sampling of three-D grain.By calculated result rendering and it is defeated FBO is arrived out, then the vertex 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, vertex is to light The ray that source issues be carried out in Release and dispersion range body bounding box after expanding in world coordinate system with certain step-length it is empty 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 system from world coordinate system: if the elevation of sampled point is less than the corresponding elevation of Rz and big In 0, then otherwise sample, abandons the sampled point within the scope of Release and dispersion.If sample within the scope of Release and dispersion, Sampled point is then transformed into diffusion coordinate system, then is transformed into voxel coordinate system from diffusion coordinate system, then mapping obtains the sampling The three-D grain coordinate of point samples three-D grain at texture coordinate.

In the step 6 of embodiment, object of a M × N × L size cuboid as texture mapping is drawn again, 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), correspond 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 language.

In vertex shader, input variable is the three-D grain coordinate gl_ of apex coordinate gl_Vertex, vertex correspondence The world coordinates Light of TexCoord, light source_worldWith the Release and dispersion range body bounding box after expansion in world coordinate system.It is right After gl_Vertex executes model view transform, projective transformation, the viewport transform, it is converted to screen coordinate gl_position.Into Using XOY plane as projection plane when row projective transformation, using the rectangular projection of face projection plane, the definition when viewport transform is carried out The viewport of one M × N size, to guarantee that the size and location of pixel in projection process remains unchanged and pixel data will not be lost It loses.The output variable of vertex shader is screen coordinate gl_position, the world of texture coordinate gl_TexCoord, light source is sat Mark Light_worldWith the Release and dispersion range body bounding box after expansion in world coordinate system.The transformed cuboid in vertex is through grating M × N × L piece member, the corresponding three-D grain coordinate gl_TexCoord of each member are obtained after change.

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

Dir=Light_world-Vertex_world

Find out the access point start and go out that light intersects with the Release and dispersion range body bounding box after expansion in world coordinate system Point end.The access point start of light is vertex Vertex in this step_world.From Vertex_worldStart, in world coordinate system In with distance s tep=0.1 along radiation direction dir sample, the sampled point obtained every time are as follows:

Pos=Vertex_world+dir×(pos+step)

Sampled point pos is successively stored in array points, until going out point end.Traversal array points counts to accumulate Calculate the transmissivity on vertex.Since the body bounding box of Release and dispersion range in world coordinate system is the bounding box after expanding, Need to judge whether sampled point is located in leakage range of scatter.Sampled point is transformed into WGS84 spherical coordinate system from world coordinate system: If the elevation of sampled point is less than the corresponding elevation of Rz and greater than 0, sample is mapped as within the scope 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) Sample obtains R channel value (i.e. concentration c), utilizes light transmittance calculation formula:

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

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

Texture2 is associated with FBO, calculated result is rendered and be output to FBO, 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, can self-setting sampling interval step according to the actual situation.

Step 7, for spreading moment t, the rectangle for finally drawing a screen size will be waited by texture mapping function Big 2 d texture object is mapped to rectangle, by GPGPU technology, realizes that the coordinate on rectangle vertex turns in vertex shader It changes, the Release and dispersion model after expanding in world coordinate system in fragment shader from viewpoint to the ray that screen pixels point issues Space point sampling is carried out with certain step-length in containment body bounding box, it is real in conjunction with the sampling to three-D grain in step 5 and step 6 The calculating of existing screen pixels point color value, rectangle is finally exported to screen, realize t moment 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 Screen pixels point issue ray be in world coordinate system expand after Release and dispersion range 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 system from world coordinate system: if the elevation of sampled point is corresponding less than Rz Elevation and be greater than 0, then otherwise sample, abandons the sampled point within the scope of Release and dispersion.If sample expands in leakage It dissipates in range, then sampled point is transformed into diffusion coordinate system, then be transformed into voxel coordinate system from diffusion coordinate system, then mapped 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) object as texture mapping, four vertex of rectangle are (0,0), (resolX-1,0), (resolX-1, ResolY-1), (0, resolY) respectively corresponds 2 d texture coordinate (0,0), (1,0), (1,1), (0,1).By 2 d texture Texture3 is mapped to rectangle.Using GLSL shading language write for texture primitive calculating screen pixels color value it is logical Use algorithm.

In vertex shader, input variable is the 2 d texture coordinate of rectangle apex coordinate gl_Vertex, vertex correspondence The world coordinates Eye of gl_TexCoord, viewpoint_worldWith the Release and dispersion range 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.Into Using XOY plane as projection plane when row projective transformation, using the rectangular projection of face projection plane, the definition when viewport transform is carried out The viewport of one resolX × resolY size, to guarantee that the size and location of pixel in projection process remains unchanged and pixel Data will not lose.The output variable of vertex shader is screen coordinate gl_position, texture coordinate gl_TexCoord, regards The world coordinates Eye of point_worldWith the Release and dispersion range body bounding box after expansion in world coordinate system.The transformed rectangle in vertex ResolX × resolY piece member, the corresponding screen coordinate gl_position of each member 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 range body bounding box after expansion in world coordinate system.By the world coordinates of viewpoint Eye_worldThe world coordinates of available center's point of screen then obtains the world coordinates at screen coordinate gl_position Screen_world.Pass through Eye_worldAnd Screen_worldIt can determine the direction of light are as follows:

Dir=Screen_world‐Eye_world

If light and the Release and dispersion range body bounding box after expanding in world coordinate system are non-intersecting, this light is abandoned Line, the color value of screen pixels are background color value bcolor；Two intersection points start, end are found out if intersection, from start Start, sampled with distance s tep=0.1 along radiation direction dir, the sampled point obtained every time are as follows:

Pos=start+dir × (pos+step)

Sampled point pos is successively stored in array points, until end.Traversal array points calculates screen to accumulate The color value of curtain pixel.Initialize transparency opacity=1, the calculation formula of color value rgb are as follows:

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).Since the body bounding box of Release and dispersion range in world coordinate system is the bounding box after expanding, it is therefore desirable to judgement sampling Whether point is located in leakage range of scatter.Sampled point is transformed into WGS84 spherical coordinate system from world coordinate system: if the height of sampled point Journey is less than the corresponding elevation of Rz and is greater than 0, then sample is within the scope of Release and dispersion, being mapped as three-D grain coordinate (u,v,w)；Otherwise, the sampled point is abandoned.Respectively three-D grain texture1 is sampled to obtain R at texture coordinate (u, v, w) Channel value (i.e. concentration c) samples three-D grain texture2 to obtain R channel value (i.e. illumination transmissivity transmittance), It utilizes:

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

Opacity=opacity × t

Accumulation obtains the transparency opacity of the light of screen pixels sending, until opacity close to 0 or is reached end.It utilizes:

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, realizes Volume drawing.

When it is implemented, can self-setting sampled point color value plocor and lighting color value according to the actual situation lcolor。

In step 5, step 6 and step 7, by the way that continuously drawing process realizes the generation of diffusion concentration data and can three times Depending on the integration changed.

Embodiment is repeated the above 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, deletes texture object texture0、texture1、texture2、texture3。

By being embodied above as can be seen that realizing that concentration value volume data and per-vertex lighting are saturating by GPGPU technology 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 The increase for needing to have certain on customized coloring process and algorithm complexity when dye, but used PBO, multi-scale sampling etc. Optimisation strategy and FBO technology, and consider the influence of three-dimensional sphere, can be realized spread under natural gas leaking diffusion conditions it is dense Degree is according to generation and the integration of the dynamic and visual based on spherical surface and illustrates 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 following procedure,

Firstly, environment is drawn in initialization, two PBO and FBO are created using OpenGL, create three-D grain object and two dimension Texture object is used to store the calculated result of GPGPU, definition diffusion coordinate system and voxel coordinate system, and initialization viewpoint, light source exist Position in WGS84 coordinate system；

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 rank of three-dimensional scenic by the sampled distance in range of scatter by CPU in diffusion coordinate system Connection, with the distance of deltaX, deltaY, deltaZ respectively in x-axis, y-axis, z-axis direction under the LOD rank of current three-dimensional scenic Sampling generates the array of vertices of t moment leakage range of scatter, and readjusts Rx, Ry, Rz, so that x-axis, y-axis, z-axis direction Number of sampling points is 2 whole power；The body bounding box for spreading Release and dispersion range in coordinate system is transformed into WGS84 coordinate System, is then transformed into world coordinate system, obtains the body bounding box of Release and dispersion range in world coordinate system；It is parallel flat by two New body bounding box after the body bounding box building of Release and dispersion range expands in face and world coordinate system, one of plane By the corresponding points of ((Rx+Minx)/2,0, Rz) in world coordinate system, another plane expands by leaking in world coordinate system Dissipate four vertex of range body bounding box bottom surface；

The array of vertices of t moment leakage range of scatter is loaded into an idle PBO by CPU, enables t=t+deltaT by step a3, DeltaT is time interval, for next current spreading moment return step a1, step a1-a3 is executed, by new diffusion The array of vertices of moment t leakage range of scatter is loaded into another PBO；

After the array of vertices of t moment leakage range of scatter is loaded into PBO by CPU, step a1-a3 is executed for next spreading moment When, array of vertices is transferred to three-D grain object identical with array of vertices size from PBO by GPU, and complete following draw Processing,

Two cuboids are drawn, by concentration value volume data and per-vertex lighting transmissivity volume data in piece by way of texture mapping It realizes that GPGPU is calculated in section tinter, calculated result is rendered and be output to FBO, then concentration value volume data and per-vertex lighting are saturating It penetrates rate volume data and is written into texture cache respectively as three-D grain；A rectangle is drawn, screen point color value data is passed through The mode of texture mapping is adopted in conjunction with 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 calculated result write-in 2 d texture is simultaneously directly rendered into screen, realizes expansion Concentration data is dissipated to generate and visual integration.

2. a kind of dynamic and visual method of the natural gas leaking diffusion based on GPGPU, feature exist as described in claim 1 In:

In step 1, diffusion coordinate system and voxel coordinate system are defined as follows,

Diffusion coordinate system is defined as with source of leaks p₀Subpoint on ground is origin, and following wind direction is x-axis, on 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 be defined as with spread in coordinate system (Minx ,-Ry, 0) be 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 respectively indicate the number of sampling points of x-axis, y-axis, z-axis.

3. a kind of dynamic and visual method of the natural gas leaking diffusion based on GPGPU, feature exist as claimed in claim 2 In: in step a1, the calculation method for spreading diffusion length Rx, Ry, Rz, Minx of x-axis, y-axis, z-axis in coordinate system is as follows,

(1) calculate Minx method be,

As t < T, Minx=0 is enabled, saves Minx value；

As t >=T,

2. enabling Minx=Minx+delta, delta is preset step-length if c < c0, the concentration at (Minx, 0, H) is recalculated Value c, until c >=c0, into 3.；If 3. c >=c0 is directly entered；

3. enabling Minx=Minx-1 if c >=c0, the concentration value c at (Minx, 0, H) is recalculated, until c < c0, is saved Minx value；

Wherein, T is natural gas leakage total time, and c0 is low danger concentration, and H is effective source height；

(2) calculate Rx method be,

2. enabling Rx=Rx+delta, delta is preset step-length if c >=c0, the concentration value c at (Rx, 0, H) is recalculated, directly To c < c0, into 3.；If 3. c < c0, is directly entered；

3. recalculating the concentration value c at (Rx, 0, H) if c < c0, enables Rx=Rx-1, until c >=c0, Rx value is saved；

(3) calculate Ry method be,

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

2. calculating the concentration value c at (Temp_x, Ry, H) enables Ry=Ry+1 if c >=c0, (Temp_x, Ry, H) is recalculated The concentration value c at place, until c < c0, into 3.；If 3. c < c0, is directly entered；

3. enabling center_x=Temp_x, Temp_x=Temp_x-1, the concentration value c at (Temp_x, Ry, H) is calculated, if c < C0 then saves Ry value；If c >=c0, return 2., repeat 2. -3.；

(4) calculate Rz method be,

1. initializing Rz=H, Temp_x=center_x, x_leff=Temp_x, x_right=Temp_x；Temp_x,x_ Left, x_right are the variable of interim storage x value, and center_x is corresponding x value at Ry；

2. calculating the concentration value c at (Temp_x, 0, Rz) enables Rz=Rz+1 if c >=c0, (Temp_x, 0, Rz) is recalculated The concentration value c at place, until c < c0, into 3.；If 3. c < c0, is directly entered；

3. enabling x_leff=Temp_x-1, x_right=Temp_x+1, concentration c 1, (x_ at judgement (x_left, 0, Rz) Right, 0, Rz) concentration c 2 and given concentration value c at₀Size,

If c1 < c0 and c2 < c0, save Rz value；

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, feature as described in Claims 2 or 3 Be: in step a2, the visualization precision for revealing range of scatter being divided into several ranks, by the sampling in range of scatter away from It is associated from the LOD rank 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, feature as described in Claims 2 or 3 It is: sets two PBO and be denoted as PBO0 and PBO1 respectively,

If CPU binds PBO0, the array of vertices generation processing of leakage range of scatter is carried out, the array of vertices of spreading moment t 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 and is drawn by GPU from PBO0, CPU binds PBO1, is let out The array of vertices generation processing for revealing range of scatter, is loaded into PBO1 for the array of vertices of next spreading moment t；

6. a kind of dynamic and visual method of the natural gas leaking diffusion based on GPGPU, feature as described in Claims 2 or 3 Be: the drawing modification of GPU includes the following steps,

Step b1, the first big cuboid of one and three-D grain object etc. of drafting, by texture mapping function, by three-D grain Object is mapped to cuboid, and by GPGPU technology, the coordinate of the array of vertices of leakage range of scatter is realized in vertex shader The general-purpose computations of the array of vertices concentration value of leakage range of scatter are realized in conversion in fragment shader；Calculated result 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 the big cuboid such as one and three-D grain object, by texture mapping function, by three-D grain again Object is mapped to cuboid, and by GPGPU technology, the coordinate of the array of vertices of leakage range of scatter is realized in vertex shader Conversion, in fragment shader, the ray issued from vertex to light source expand in world coordinate system after Release and dispersion range Space point sampling is carried out by step-length in body bounding box, in conjunction with the sampling to three-D grain in step b2, realizes Release and dispersion body number According to the general-purpose computations of vertex light transmittance, calculated result is rendered to and is output to FBO, then the vertex 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 technology, the coordinate conversion on rectangle vertex realized in vertex shader, in fragment shader by Viewpoint expand in world coordinate system to the ray that screen pixels point issues after Release and dispersion range body bounding box in by step-length Space point sampling is carried out, in conjunction with the sampling to three-D grain in step b1 and step b2, realizes the meter of screen pixels point color value It calculates, finally rectangle is exported to screen, realizes the three-dimensional visualization of t moment natural gas leaking range of scatter.

7. a kind of dynamic and visual method of the natural gas leaking diffusion based on GPGPU, feature exist as claimed in claim 6 In: it in step b2, is sampled in the bounding box after expansion, needs to judge whether sampled point is located in leakage range of scatter, packet It includes and sampled point is transformed into WGS84 coordinate system from world coordinate system, if the elevation of sampled point is less than the corresponding elevation of Rz and is greater than 0, then otherwise sample abandons the sampled point within the scope of Release and dispersion；If sample is incited somebody to action within the scope 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 three-D grain at texture coordinate.

8. a kind of dynamic and visual method of the natural gas leaking diffusion based on GPGPU, feature exist as claimed in claim 6 In: it in step b3, is sampled in the bounding box after expansion, needs to judge whether sampled point is located in leakage range of scatter, packet It includes and sampled point is transformed into WGS84 spherical coordinate system from world coordinate system, if the elevation of sampled point is less than the corresponding elevation of Rz and big In 0, then otherwise sample abandons the sampled point within the scope of Release and dispersion；If sample within the scope of Release and dispersion, Sampled point is transformed into diffusion coordinate system, then is transformed into voxel coordinate system from diffusion coordinate system, then mapping obtains the sampled point Three-D grain coordinate, at texture coordinate to three-D grain sample.