CN104634761B - Suspended solution radiation characteristic parameter inversion method based on GPU parallel acceleration - Google Patents

Abstract

The invention discloses a suspended solution radiation characteristic parameter inversion method based on GPU parallel acceleration, and relates to a suspended solution radiation characteristic parameter method based on GPU parallel acceleration. The suspended solution radiation characteristic parameter inversion method aims to solve the problems that the accuracy rate of measurement results of suspended solution radiation characteristic parameters is low and a great number of computers are needed in the prior art. The suspended solution radiation characteristic parameter inversion method is specifically performed by virtue of the following steps: step 1, preparing a to-be-measured suspended solution, and feeding the to-be-measured suspended solution into a sample container which is made from quartz glass; step 2, measuring BSDF of the sample container filled with the suspended solution to acquire a group of BSDF experimental measurement data BSDFexp in different scattering directions; and step 3, performing inversion on radiation characteristic parameters of the sample container suspended solution based on the combination of a GPU acceleration algorithm and an optimization algorithm, wherein the radiation characteristic parameters of the suspended solution comprise extinction coefficients beta, scattering coefficients sigma s and asymmetry factors g. The suspended solution radiation characteristic parameter inversion method disclosed by the invention can be applied to the technical field of optical characteristic measurement of the suspended solutions.

Description

The method of the aaerosol solution Radiation Characteristics Parameters inverting for being accelerated based on GPU parallel

Technical field

The present invention relates to the method for the parallel aaerosol solution Radiation Characteristics Parameters for accelerating of GPU.

Background technology

The conventional single thread MC simulation based on CPU is calculated and needs substantial amounts of computer, tall and handsome up to company from 2006 (NVIDIA) isomery (CUDA) parallel computation framework is proposed, makes to obtain using the programming language of figure video card (GPU) parallel computation Great simplification, at the same time the performance of video card also constantly improving, make MC methods solve the simulation of equation of radiative transfer Time is greatly reduced, and this is to combine optimized algorithm using MC methods (such as particle swarm optimization algorithm PSO, Genetic Algorithms) Carry out inverting to lay the foundation.

In general for the method that suspended particles solution measures its Radiation Characteristics Parameters has following two：

Method measured directly and aaerosol solution radiation characteristic is indirectly obtained with reference to inverse problem model by scattered signal Method；

Method measured directly requires that aaerosol solution has scattering properties, it is also contemplated that the container of filling solution is to scattering letter Number impact, cause aaerosol solution Radiation Characteristics Parameters measurement result accuracy rate low.

Expend substantial amounts of by scattered signal with reference to the method that inverse problem model indirectly obtains aaerosol solution radiation characteristic Computer.

The content of the invention

The present invention is low in order to solve aaerosol solution Radiation Characteristics Parameters measurement result accuracy rate in prior art, and expends The problem of substantial amounts of computer, and the method for proposing the aaerosol solution Radiation Characteristics Parameters inverting for accelerating based on GPU parallel.

Above-mentioned goal of the invention is achieved through the following technical solutions：

Step one, preparation aaerosol solution to be measured, aaerosol solution to be measured are mounted in the shuttle that quartz glass does, are tested Aaerosol solution to be measured is well mixed in measurement process；

The BSDF of step 2, shuttle of the measurement equipped with aaerosol solution to be measured：

LASER Light Source is incided into shuttle left surface along the direction vertical with sample holder surface, using being arranged in The scattering in the different scattering directions in the range of right flank of the rotatable detector measurement shuttle left surface outside shuttle Signal, shuttle left surface to right flank scope are 0 degree～180 degree, that is, obtain one group of BSDF experiment in different scattering directions Measurement data BSDF_exp；

BSDF is two-way dispersion distribution function, and which is defined as follows；

In formula, θ_iFor incidence zenith angle,For incident orientation angle, θ_sFor scattering zenith angle,For scattering azimuth, sr is Three-dimensional angular unit, L_iIt is the intensity of incident radiation of unit solid angle, dL_sScattering radiation intensity is represented, d represents differentiating operator, ω X represents incident solid angles；

Step 3, the Radiation Characteristics Parameters for carrying out shuttle aaerosol solution based on GPU accelerating algorithms with reference to optimized algorithm Inverting, wherein, the Radiation Characteristics Parameters of the aaerosol solution are extinction coefficient β, scattering coefficient σ_sWith dissymmetry factor g；

1) initial value of the Radiation Characteristics Parameters of aaerosol solution, the initial value β ', scattering coefficient σ including extinction coefficient β are set_s's Initial value σ_s' and dissymmetry factor g initial value g ', then use MC algorithms based on GPU in the shuttle equipped with suspension Radiation transmittance process solved, obtain one group of BSDF analogue data BSDF on shuttle light scattering different directions_sim；

Wherein, the β ' and σ_s' it is to randomly select numerical value；

G ' is to randomly select numerical value in the range of [- 1,1]；

Derivation algorithms of the MC for direct problem equation of radiative transfer；

GPU is graphic process unit；

2) initial value β ', scattering coefficient σ of the optimized algorithm to extinction coefficient β are used_sInitial value σ_s' first with dissymmetry factor g Value extinction coefficient β after value g ', scattering coefficient σ_sIt is optimized with dissymmetry factor g, the value of object function F (x) is constantly subtracted It is little；

When the value of object function F (x) is less than the accuracy value of setting or reaches the iterative steps of setting, then stop optimization, Extinction coefficient β, the scattering coefficient σ obtained by inverting_sWith the value of dissymmetry factor g as shuttle Radiation Characteristics Parameters.

Invention effect

Using the method for the aaerosol solution Radiation Characteristics Parameters inverting for being accelerated based on GPU parallel of the present invention,

(1) present invention is surveyed by setting up the solution of inverse problems model that aaerosol solution shuttle scattered signal is measured, experiment Amount obtains one group of BSDF measured data of experiment BSDF in different scattering directions_exp；Can be to radiation using the MC methods based on GPU Intensity field carries out three-dimensional solution, obtains one group of BSDF analogue data BSDF on shuttle light scattering different directions_sim；Solve Equation of radiative transfer combines optimized algorithm, the initial value, scattering coefficient σ to extinction coefficient β_sInitial value and dissymmetry factor g just Value after value is optimized, and the value of object function F (x) is constantly reduced；When the value of object function F (x) is less than the precision for setting Value or when reaching the iterative steps of setting, then stop optimization, extinction coefficient β, the scattering coefficient σ obtained by inverting_sWith it is asymmetric Radiation Characteristics Parameters of the value of factor g as shuttle, the Radiation biodosimetry for completing the aaerosol solution based on GPU are solved, The method can realize load aaerosol solution shuttle on the scattered signal of aaerosol solution without impact, solve aaerosol solution The low problem of Radiation Characteristics Parameters measurement result accuracy rate, improves aaerosol solution Radiation Characteristics Parameters measurement result accuracy rate 22%.

(2) present invention obtains sample with reference to the powerful computation capability of GPU using LASER Light Source irradiating sample container The scattered signal of container, combines the Radiation Characteristics Parameters that solution of inverse problems algorithm obtains shuttle indirectly based on scattered signal, Radiation Characteristics Parameters after the initial value of the Radiation Characteristics Parameters for optimizing aaerosol solution using particle swarm optimization algorithm, the Algorithm for Solving Have the advantages that during optimization problem that little by initial value affecting and sensitivity is high, it is not necessary to expend substantial amounts of computer, a computer The radiation characteristic inversion method of the aaerosol solution based on GPU parallel computations can just be completed.

The method of the aaerosol solution Radiation Characteristics Parameters inverting for being accelerated based on GPU parallel proposed in the present invention is in time The speed-up ratio of more than 2 orders of magnitude is realized compared with CPU is calculated, and GPU proposed by the present invention is thick in optics compared to existing CPU Degree take 0.1, albedo take 0.1 and dissymmetry factor take under conditions of 0.9, it is possible to achieve 163.5 times of speed-up ratio.

Description of the drawings

Fig. 1 is flow chart of the present invention；

Fig. 2 is scattered signal measurement when being subject to laser to irradiate on the left of aaerosol solution sample cell described in specific embodiment one Schematic diagram；Three layer model in figure positioned at middle body is aaerosol solution sample cell, and a week of periphery is detector, and A is laser Light source, B are aaerosol solution, and C is glass, and D is detector；X is reference axis X-axis, and y is reference axis y-axis, and z is reference axis z-axis, and o is Reference axis o axle；

Fig. 3 is BSDF Data Comparison figures, and abscissa is angle of scattering, and for spending, ordinate is two-way dispersion distribution function to unit, Unit is sr^-1；

Fig. 4 is object function convergence graph, and abscissa is iteration step, and ordinate is object function；

Fig. 5 is optical thickness convergence graph, and abscissa is iteration step, and ordinate is optical thickness；

Fig. 6 is albedo convergence graph, and abscissa is iteration step, and ordinate is albedo；

Fig. 7 is dissymmetry factor convergence graph, and abscissa is iteration step, and ordinate is dissymmetry factor.

Specific embodiment

Specific embodiment one：Present embodiment is illustrated with reference to Fig. 1, the aaerosol solution radiation for accelerating based on GPU parallel is special What the method for property parametric inversion was specifically followed the steps below：

Step one, preparation aaerosol solution to be measured, aaerosol solution to be measured are mounted in the shuttle that quartz glass does, are tested Aaerosol solution to be measured is well mixed in measurement process；

Step 2, measurement equipped with aaerosol solution to be measured shuttle BSDF, by LASER Light Source along with shuttle Shuttle left surface is incided in the vertical direction in surface, using the rotatable detector measurement sample being arranged in outside shuttle Product container left surface difference in the range of right flank scatters the scattered signals in directions, and shuttle left surface to right flank scope is 0 degree～180 degree, that is, obtain one group of BSDF measured data of experiment BSDF in different scattering directions_exp；Such as Fig. 2；

BSDF is two-way dispersion distribution function, and which is defined as follows；

In formula, θ_iFor incidence zenith angle,For incident orientation angle, θ_sFor scattering zenith angle,For scattering azimuth, sr is Three-dimensional angular unit, L_iIt is the intensity of incident radiation of unit solid angle, dL_sScattering radiation intensity is represented, d represents differentiating operator, ω_iRepresent incident solid angles；

Step 3, the Radiation Characteristics Parameters for carrying out shuttle aaerosol solution based on GPU accelerating algorithms with reference to optimized algorithm Inverting, wherein, the Radiation Characteristics Parameters of the aaerosol solution are extinction coefficient β, scattering coefficient σ_sWith dissymmetry factor g；

1) initial value of the Radiation Characteristics Parameters of aaerosol solution, the initial value β ', scattering coefficient σ including extinction coefficient β are set_s's Initial value σ_s' and dissymmetry factor g initial value g ', then use MC algorithms based on GPU in the shuttle equipped with suspension Radiation transmittance process solved, obtain one group of BSDF analogue data BSDF on shuttle light scattering different directions_sim；

Wherein, the β ' and σ_s' it is to randomly select numerical value；

G ' is to randomly select numerical value in the range of [- 1,1]；

Derivation algorithms of the MC for direct problem equation of radiative transfer；

GPU is graphic process unit；

2) initial value β ', scattering coefficient σ of the optimized algorithm to extinction coefficient β are used_sInitial value σ_s' first with dissymmetry factor g Value extinction coefficient β after value g ', scattering coefficient σ_sIt is optimized with dissymmetry factor g, the value of object function F (x) is constantly subtracted It is little；

When the value of object function F (x) is less than the accuracy value of setting or reaches the iterative steps of setting, then stop optimization, Extinction coefficient β, the scattering coefficient σ obtained by inverting_sWith the value of dissymmetry factor g as shuttle Radiation Characteristics Parameters.

Optimized algorithm is particle swarm optimization algorithm PSO or Genetic Algorithms.

Specific embodiment two：Present embodiment from unlike specific embodiment one：Assume in the step 3 outstanding The initial value of the Radiation Characteristics Parameters of floating solution, including extinction coefficient β, scattering coefficient σ_sWith the initial value of dissymmetry factor g, then make The radiation transmittance process in the shuttle equipped with suspension is solved with the MC algorithms based on GPU, obtain shuttle One group of BSDF analogue data BSDF on light scattering different directions_sim；Detailed process is：

Radiation transmittance process in shuttle is solved using radiation transfer equation：

In formula, I (r, s) is radiation intensity, and r is radiation field position vector, and s is direction vector, and β ' is the first of extinction coefficient β Value, Φ (s ' → s) be Scattering Phase Function, Ω ' be solid angle, σ_s' it is scattering coefficient σ_sInitial value, Φ is Scattering Phase Function, and Θ is Angle of scattering,For Hamiltonian operator, initial values of the g ' for dissymmetry factor g.

Other steps and parameter are identical with specific embodiment one.

Specific embodiment three：Present embodiment from unlike specific embodiment one or two：Mesh in the step 3 Scalar functions F (x) are：

With reference to below equation：

In formula, τ is optical thickness, and ω is albedo, and g is dissymmetry factor, and g spans are [- 1,1], N_dIt is use Scattering angle quantity.

Other steps and parameter are identical with specific embodiment one or two.

Beneficial effects of the present invention are verified using following examples：

Embodiment 1：

The method of the aaerosol solution Radiation Characteristics Parameters inverting for being accelerated based on GPU parallel is specifically followed the steps below 's：

Step one, preparation aaerosol solution to be measured, aaerosol solution to be measured are mounted in the shuttle that quartz glass does, are tested Aaerosol solution to be measured is well mixed in measurement process；

Step 2, measurement equipped with aaerosol solution shuttle BSDF, by LASER Light Source along with sample holder surface Shuttle left surface is incided in vertical direction, is held using the rotatable detector measurement sample being arranged in outside shuttle The scattered signal in device left surface different scattering directions in the range of right flank, shuttle left surface to right flank scope is 0 degree ～180 degree, that is, obtain one group of BSDF measured data of experiment BSDF in different scattering directions_exp；

BSDF is two-way dispersion distribution function, and which is defined as follows；

In formula, θ_iFor incidence zenith angle,For incident orientation angle, θ_sFor scattering zenith angle,For scattering azimuth, sr is Three-dimensional angular unit, L_iIt is the intensity of incident radiation of unit solid angle, dL_sScattering radiation intensity is represented, d represents differentiating operator, ω_iRepresent incident solid angles；

Step 3, the Radiation Characteristics Parameters for carrying out shuttle aaerosol solution based on GPU accelerating algorithms with reference to optimized algorithm Inverting, wherein, the Radiation Characteristics Parameters of the aaerosol solution are extinction coefficient β, scattering coefficient σ_sWith dissymmetry factor g；

1) initial value of the Radiation Characteristics Parameters of aaerosol solution, the initial value β ', scattering coefficient σ including extinction coefficient β are set_s's Initial value σ_s' and dissymmetry factor g initial value g ', then use MC algorithms based on GPU in the shuttle equipped with suspension Radiation transmittance process solved, obtain one group of BSDF analogue data BSDF on shuttle light scattering different directions_sim；

Wherein, the β ' and σ_s' it is to randomly select numerical value；

G ' is to randomly select numerical value in the range of [- 1,1]；

Derivation algorithms of the MC for direct problem equation of radiative transfer；

GPU is graphic process unit；

2) initial value β ', scattering coefficient σ of the optimized algorithm to extinction coefficient β are used_sInitial value σ_s' first with dissymmetry factor g Value extinction coefficient β after value g ', scattering coefficient σ_sIt is optimized with dissymmetry factor g, the value of object function F (x) is constantly subtracted It is little；

When the value of object function F (x) is less than the accuracy value (such as 10 for setting^-2) or reach the iterative steps (such as 30 steps) of setting When, then stop optimization, extinction coefficient β, the scattering coefficient σ obtained by inverting_sWith the value of dissymmetry factor g as shuttle Radiation Characteristics Parameters.

Optimized algorithm is particle swarm optimization algorithm PSO or Genetic Algorithms.

Fig. 3 be τ=1.0, ω=0.5, under g=0.9 Parameter Conditions a inversion result, in order to verify inversion result Stability, has carried out three invertings to this group of parameter, can be seen that inversion result stability is fine from the convergence graph of inversion result；

Fig. 4 be τ=1.0, ω=0.5, under g=0.9 Parameter Conditions a inversion result, in order to verify inversion result Stability, has carried out three invertings to this group of parameter, can be seen that inversion result stability is fine from the convergence graph of inversion result；

Fig. 5 be τ=1.0, ω=0.5, under g=0.9 Parameter Conditions a inversion result, in order to verify inversion result Stability, has carried out three invertings to this group of parameter, can be seen that inversion result stability is fine from the convergence graph of inversion result；

Fig. 6 be τ=1.0, ω=0.5, under g=0.9 Parameter Conditions a inversion result, in order to verify inversion result Stability, has carried out three invertings to this group of parameter, can be seen that inversion result stability is fine from the convergence graph of inversion result；

Fig. 7 be τ=1.0, ω=0.5, under g=0.9 Parameter Conditions a inversion result, in order to verify inversion result Stability, has carried out three invertings to this group of parameter, can be seen that inversion result stability is fine from the convergence graph of inversion result；

The present invention is combined the method for the aaerosol solution Radiation Characteristics Parameters inverting for being accelerated based on GPU parallel and is compared in time The speed-up ratio being capable of achieving more than 2 orders of magnitude is calculated in existing CPU, and the GPU (GTX 660) that the present invention is used is compared to existing CPU (Intel i53570K 3.4GHz) takes 0.1 albedo in optical thickness and dissymmetry factor is taken under conditions of 0.9, can To realize 163.5 times of speed-up ratio, for the hardware used in the present invention can be reduced to the Inversion Calculation time 30 minutes with It is interior, if can be Inversion Calculation time restriction within 10 minutes using GPU expections best in the market.

Claims (1)

1. the method for the aaerosol solution Radiation Characteristics Parameters inverting for being accelerated based on GPU parallel, it is characterised in that：It is parallel based on GPU What the method for the aaerosol solution Radiation Characteristics Parameters inverting of acceleration was specifically followed the steps below：

Step one, preparation aaerosol solution to be measured, aaerosol solution to be measured is mounted in the shuttle that quartz glass does；

The BSDF of step 2, shuttle of the measurement equipped with aaerosol solution to be measured：

LASER Light Source is incided into shuttle left surface along the direction vertical with sample holder surface, using being arranged in sample The scattered signal in the rotatable detector measurement shuttle left surface of external container different scattering directions in the range of right flank, Shuttle left surface is 0 degree～180 degree to right flank scope, that is, obtain one group of BSDF experiment measurement number in different scattering directions According to BSDF_exp；

BSDF is two-way dispersion distribution function, and which is defined as follows；

In formula, θ_iFor incidence zenith angle,For incident orientation angle, θ_sFor scattering zenith angle,For scattering azimuth, sr is solid Angular unit, L_iIt is the intensity of incident radiation of unit solid angle, dL_sScattering radiation intensity is represented, d represents differentiating operator, ω_iTable Show incident solid angles；

Step 3, carried out with reference to optimized algorithm based on GPU accelerating algorithms shuttle aaerosol solution Radiation Characteristics Parameters it is anti- Drill, wherein, the Radiation Characteristics Parameters of the aaerosol solution are extinction coefficient β, scattering coefficient σ_sWith dissymmetry factor g；

1) initial value of the Radiation Characteristics Parameters of aaerosol solution, the initial value β ', scattering coefficient σ including extinction coefficient β are set_sInitial value σ_s' and dissymmetry factor g initial value g ', then use MC algorithms based on GPU to the spoke in the shuttle equipped with suspension Penetrate transmittance process to be solved, obtain one group of BSDF analogue data BSDF on shuttle light scattering different directions_sim；

Wherein, the β ' and σ_s' it is to randomly select numerical value；

G ' is to randomly select numerical value in the range of [- 1,1]；

Derivation algorithms of the MC for direct problem equation of radiative transfer；

GPU is graphic process unit；

2) using optimized algorithm to extinction coefficient β, scattering coefficient σ_sIt is optimized with the value of dissymmetry factor g, makes object function F X the value of () constantly reduces；

When the value of object function F (x) is less than the accuracy value of setting or reaches the iterative steps of setting, then stop optimization, anti- Drill extinction coefficient β, the scattering coefficient σ for obtaining_sWith the value of dissymmetry factor g as shuttle Radiation Characteristics Parameters；

The initial value of the Radiation Characteristics Parameters of aaerosol solution is assumed in the step 3, including extinction coefficient β, scattering coefficient σ_sNo The initial value of symmetrical factor g, then uses the MC algorithms based on GPU to the radiation transmittance process in the shuttle equipped with suspension Solved, obtained one group of BSDF analogue data BSDF on shuttle light scattering different directions_sim；Detailed process is：

Radiation transmittance process in shuttle is solved using radiation transfer equation：

$s\·\▿I(r,s)+\mathrm{\β}I(r,s)=\frac{{\mathrm{\σ}}_s}{4\mathrm{\π}}{\∫}_{{\mathrm{\Ω}}^{\′}=4\mathrm{\π}}I(r,s)\mathrm{\Φ}(s^{\′}\→s){\mathrm{d\Ω}}^{\′}---\left(1\right)$

$\mathrm{\Φ}(s^{\′}\→s)=\mathrm{\Φ}\left(cos\mathrm{\Θ}\right)=\frac{1-g^2}{{(1+g^2-2gcos\mathrm{\Θ})}^{3/2}}---\left(2\right)$

In formula, I (r, s) is radiation intensity, and r is radiation field position vector, and s is direction vector, and β is extinction coefficient, Φ (s ' → s) For Scattering Phase Function, Ω ' is solid angle, σ_sFor scattering coefficient, Φ is Scattering Phase Function, and Θ is angle of scattering,Calculate for Hamilton Symbol, g is dissymmetry factor；

In the step 3, object function F (x) is：

With reference to below equation：

$F\left(x\right)=F(\mathrm{\τ},\mathrm{\ω},g)=\underset{i=1}{\overset{N_d}{\Σ}}{v}_i^2{({\mathrm{BSDF}}_{i,sim}-{\mathrm{BSDF}}_{i,\exp })}^2---\left(3\right)$

$v_i=\frac{1}{{\mathrm{BSDF}}_{\exp ,i}},i=1,...,N_d---\left(4\right)$

In formula, τ is optical thickness, and ω is albedo, and g is dissymmetry factor, and g spans are [- 1,1], N_dIt is the scattering for using Number of angles.