CN103353923B - Adaptive space interpolation method and system thereof based on space characteristics analysis - Google Patents

Abstract

The invention discloses a kind of adaptive space interpolation method based on spatial distribution characteristic, including: extract pollutant monitoring value；Interpolation model is set up based on variation function；Incoming direction feature is to interpolation model and variation function；Judge that the spatial distribution of pollutant is as isotropism or anisotropy；Based on isotropism or anisotropy setting search strategy to determine that ginseng is estimated a little；The pollutant levels property value a little determining interpolation point is estimated according to the variation function expression formula in all directions and ginseng.Compared with prior art, interpolation method of the present invention is to generate after sample point attribute data is made space characteristics analysis, under considering the premise of spatial distribution characteristic of sample point attribute data, ensure that amount of calculation is little, and improve interpolation precision, thus the attribute character value of predictive study region unknown point more exactly.The present invention discloses a kind of adaptive space interplotation system based on spatial distribution characteristic.

Description

Adaptive space interpolation method and system thereof based on space characteristics analysis

Technical field

The present invention relates to environmental monitoring technology field, relate more specifically to a kind of adaptive space interpolation method based on space characteristics analysis and system thereof.

Background technology

On current automatic inspection criterion, points all in space cannot be observed by the points distributing method of regulation, but, based on the monitoring network laid, a number of space sample can be obtained, these samples reflect all or part of feature of spatial distribution, according to these known sample, adopt suitable spatial interpolation methods, it is possible to the accurately feature of the unknown geographical space of prediction.

Common spatial interpolation methods has anti-distance weighting interpolation method (IDW), Kriging regression method (Kriging).Above two spatial interpolation methods is respectively present following defect: IDW method only considers the distance of space point-to-point transmission, does not consider spatial distribution characteristic, and therefore interpolation precision is generally relatively low；Kriging method is although it is contemplated that spatial distribution characteristic, but it is computationally intensive to relate to solving equations, the problem that interpolation speed is slow.Atmosphere environment supervision website is generally discrete, uneven, limited amount, adopt above-mentioned spatial interpolation algorithm can there is interpolation precision not high, or can not accurately reflecting the feature of unknown point attribute character, this is also that current regional air mass space is distributed the important problem faced.

Therefore, be necessary to provide a kind of sample point attribute data is made further space characteristics analysis, generation considers the novel interpolation method of spatial distribution characteristic and system thereof, and this interpolation method and system thereof can guarantee that amount of calculation is little, and interpolation precision is higher, thus the attribute character value of predictive study region unknown point more exactly.

Summary of the invention

It is an object of the invention to provide a kind of adaptive space interpolation method based on spatial distribution characteristic, this interpolation method considers the spatial distribution characteristic of sample point attribute data, amount of calculation is little, and can improve interpolation precision, thus the attribute character value of predictive study region unknown point more exactly.

It is a further object of the present invention to provide a kind of adaptive space interplotation system based on spatial distribution characteristic, this interplotation system considers the spatial distribution characteristic of sample point attribute data, amount of calculation is little, and can improve interpolation precision, thus the attribute character value of predictive study region unknown point more exactly.

To achieve these goals, the invention provides a kind of adaptive space interpolation method based on spatial distribution characteristic, including:

S1: extract geographical location information and pollutant monitoring value；

S2: based on Geostatistical principle, the spatial distribution characteristic of analyzed area ambient air quality；

S3: set up adaptive space interpolation model based on variation function, the expression formula of described adaptive space interpolation model is:

$Z\left(x_0\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\λ}}_{i}Z\left(x_i\right)=\frac{{\mathrm{\Σ}}_{i=1}^n\frac{1}{\mathrm{\γ}\left(d_i\right)}Z\left(x_i\right)}{{\mathrm{\Σ}}_{i=1}^n\frac{1}{\mathrm{\γ}\left(d_i\right)}}---\left(1\right)$

Wherein, λ_iFor known sample point x_iWeight coefficient；Z (x_o) for the pollutant levels property value of interpolation point, i is natural number, and n is the sample point sum affecting interpolation point, d_iFor interpolation point x_oAn x is estimated with ginseng_iBetween distance, Z (x_i) estimate pollutant levels property value a little, γ (d for ginseng_i) it is with d_iVariation function for variable；

S4: incoming direction feature to described adaptive space interpolation model obtains formula (2):

$Z\left(x_o\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\λ}}_{i}Z\left(x_i\right)=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{1}{\mathrm{\γ}(d_i,\mathrm{\θ})}Z\left(x_i\right)}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{1}{\mathrm{\γ}(d_i,\mathrm{\θ})}}---\left(2\right)$

Wherein, λ_iFor known sample point x_iWeight coefficient；Z (x_o) for the pollutant levels property value of interpolation point, i is natural number, and n is the sample point sum affecting interpolation point, d_iFor interpolation point x_oAn x is estimated with ginseng_iBetween distance, θ is interpolation point x_oAn x is estimated with ginseng_iBetween position angle, γ (d_i, θ) and it is with d_i, θ be the variation function of variable；

S5: to described variation function incoming direction feature to obtain the variation function expression formula in all directions；

S6: judge that the spatial distribution of described pollutant is as isotropism or anisotropy according to the variation function expression formula in described all directions；

S7: based on isotropism or anisotropy setting search strategy to determine that described ginseng estimates an x_i；

S8: estimate the pollutant levels property value Z (x a little determining interpolation point in formula (2) according to the variation function expression formula in described all directions and described ginseng_o)。

Compared with prior art, the present invention first extracts geographical location information and pollutant monitoring value based on the adaptive space interpolation method of spatial distribution characteristic, again based on the spatial distribution characteristic of Geostatistical principle analysis regional environment air quality, it is then based on variation function and sets up adaptive space interpolation model, incoming direction feature is to adaptive space interpolation model and variation function afterwards, judge that the spatial distribution of pollutant is as isotropism or anisotropy further according to the variation function of incoming direction feature, and based on isotropism and anisotropy setting search strategy to determine that ginseng is estimated a little, final variation function and ginseng according to all directions estimates the pollution concentration property value a little determining interpolation point；This interpolation method is the novel interpolation method considering spatial distribution characteristic generated after sample point attribute data is made further space characteristics analysis, under considering the premise of spatial distribution characteristic of sample point attribute data data, amount of calculation is little, and improve interpolation precision, thus the attribute character value of predictive study region unknown point more exactly.

Correspondingly, present invention also offers a kind of adaptive space interplotation system based on spatial distribution characteristic, including:

Extraction module, is used for extracting geographical location information and pollutant monitoring value；

Analysis module, for based on Geostatistical principle, the spatial distribution characteristic of analyzed area ambient air quality；

Model building module, for setting up adaptive space interpolation model based on variation function, the expression formula of described adaptive space interpolation model is:

$Z\left(x_0\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\λ}}_{i}Z\left(x_i\right)=\frac{{\mathrm{\Σ}}_{i=1}^n\frac{1}{\mathrm{\γ}\left(d_i\right)}Z\left(x_i\right)}{{\mathrm{\Σ}}_{i=1}^n\frac{1}{\mathrm{\γ}\left(d_i\right)}}---\left(1\right)$

Wherein, λ_iFor known sample point x_iWeight coefficient；Z (x_o) for the pollutant levels property value of interpolation point, i is natural number, and n is the sample point sum affecting interpolation point, d_iFor interpolation point x_oAn x is estimated with ginseng_iBetween distance, Z (x_i) estimate pollutant levels property value a little, γ (d for ginseng_i) it is with d_iVariation function for variable；

First introduces module, obtains formula (2) for incoming direction feature to described adaptive space interpolation model:

$Z\left(x_o\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\λ}}_{i}Z\left(x_i\right)=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{1}{\mathrm{\γ}(d_i,\mathrm{\θ})}Z\left(x_i\right)}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{1}{\mathrm{\γ}(d_i,\mathrm{\θ})}}---\left(2\right)$

Wherein, λ_iFor known sample point x_iWeight coefficient；Z (x_o) for the pollutant levels property value of interpolation point, i is natural number, and n is the sample point sum affecting interpolation point, d_iFor interpolation point x_oAn x is estimated with ginseng_iBetween distance, θ is interpolation point x_oAn x is estimated with ginseng_iBetween position angle, γ (d_i, θ) and it is with d_i, θ be the variation function of variable；

Second introduce module, for described variation function incoming direction feature to obtain the variation function expression formula in all directions；

Determination module, for judging that according to the variation function expression formula in described all directions the spatial distribution of described pollutant is as isotropism or anisotropy；

Search module, for based on isotropism or anisotropy setting search strategy to determine that described ginseng estimates an x_i；

Determine module, for estimating the pollutant levels property value Z (x a little determining interpolation point in formula (2) according to the variation function expression formula in described all directions and described ginseng_o)。

By description below and in conjunction with accompanying drawing, the present invention will become more fully apparent, and these accompanying drawings are used for explaining embodiments of the invention.

Accompanying drawing explanation

Fig. 1 is the present invention flow chart based on adaptive space interpolation method one embodiment of space characteristics analysis.

The variation function scatterplot that Fig. 2 draws when being and variation function is fitted.

The fitting result chart of three kinds of models when Fig. 3 is that variation function is fitted.

Fig. 4 is the sub-process figure of S5 in Fig. 1.

Fig. 5 is the variation function scatterplot on different directions on eight directions in step S54 shown in Fig. 4, and the model that wherein the figure shows is spherical model, and pollutant are SO₂。

Fig. 6 is the sub-process figure of S6 in Fig. 1.

Fig. 7 is the rose diagram of S61 in Fig. 6.

Fig. 8 is the sub-process figure of S7 in Fig. 1.

Fig. 9 is the present invention structured flowchart based on adaptive space interpolation method one embodiment of space characteristics analysis.

Figure 10 is the structured flowchart of fitting module shown in Fig. 9.

Figure 11 is the structured flowchart of the second introducing module shown in Fig. 9.

Figure 12 is the structured flowchart of determination module shown in Fig. 9.

Figure 13 is the structured flowchart of search module shown in Fig. 9.

Detailed description of the invention

With reference now to accompanying drawing, describing embodiments of the invention, element numbers similar in accompanying drawing represents similar element.

Refer to Fig. 1, the present invention comprises the steps: based on the adaptive space interpolation method of space characteristics analysis

S1: according to the atmospheric environment automatic monitoring data in single timeslice, extracts geographical location information and pollutant monitoring value, and wherein, geographical location information and pollutant monitoring value need one_to_one corresponding；Specifically, with 62, region, Pearl River Delta on the 24th August in 2012 website SO₂Monitoring Data is example, first extracts geographical location information and the SO of 62 websites₂Concentration information, adopts MATLAB software, is respectively adopted x, y, z coordinate and represents the longitude of each website, latitude and pollutant levels data respectively, specific as follows:

113.624623.65028

113.234723.14226

113.259723.13314

113.276523.15448

113.261223.105032

113.322023.13196

S2: based on Geostatistical principle, the spatial distribution characteristic of analyzed area ambient air quality；

S3: set up adaptive space interpolation model based on variation function, the expression formula of adaptive space interpolation model is:

$Z\left(x_0\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\λ}}_{i}Z\left(x_i\right)=\frac{{\mathrm{\Σ}}_{i=1}^n\frac{1}{\mathrm{\γ}\left(d_i\right)}Z\left(x_i\right)}{{\mathrm{\Σ}}_{i=1}^n\frac{1}{\mathrm{\γ}\left(d_i\right)}}---\left(1\right)$

Wherein, λ_iFor known sample point x_iWeight coefficient；Z (x_o) for the pollutant levels property value of interpolation point, i is natural number, and n is the sample point sum affecting interpolation point, d_iFor interpolation point x_oAn x is estimated with ginseng_iBetween distance, Z (x_i) estimate pollutant levels property value a little, γ (d for ginseng_i) it is with d_iVariation function for variable；

Specifically, the expression formula of variation function is:

Nondirectional ideal expression

$\mathrm{\γ}\left(h\right)=\mathrm{\γ}\left(\right|x_i-x_j\left|\right)=\frac{1}{2}E{[Z\left(x_i\right)-Z\left(x_j\right)]}^2---\left(8\right)$

Nondirectional discrete expression

${\mathrm{\γ}}^{*}\left(h\right)={\mathrm{\γ}}^{*}\left(\right|x_i-x_j\left|\right)=\frac{1}{2m}{\mathrm{\Σ}}_{k=1}^m{[Z\left(x_i\right)-Z\left(x_j\right)]}^2---\left(9\right)$

Wherein γ (h) is desirable variation function, γ^*H () is discrete variation function, x_i,x_jFor the sample point pair of distance h, Z (x_i)、Z(x_j) for sample point x_i,x_jProperty value, h is sample point x_i,x_jBetween distance, (x_i',y_i') for known sample point x_iLatitude and longitude coordinates, (x_j',y_j') for known sample point x_jLatitude and longitude coordinates, m represents the number of the sample point pair that distance is h；

Distance h adopts following formula to be calculated:

$h=\sqrt{{({x_i}^{\′}-{x_j}^{\′})}^2+{({y_i}^{\′}-{y_j}^{\′})}^2}---\left(7\right)$

Wherein, (x_i',y_i') for known sample point x_iLatitude and longitude coordinates, (x_j',y_j') for known sample point x_jLatitude and longitude coordinates.

S4: incoming direction feature to described adaptive space interpolation model obtains formula (2):

$Z\left(x_o\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\λ}}_{i}Z\left(x_i\right)=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{1}{\mathrm{\γ}(d_i,\mathrm{\θ})}Z\left(x_i\right)}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{1}{\mathrm{\γ}(d_i,\mathrm{\θ})}}---\left(2\right)$

Wherein, λ_iFor known sample point x_iWeight coefficient；Z (x_o) for the pollutant levels property value of interpolation point, i is natural number, and n is the sample point sum affecting interpolation point, d_iFor interpolation point x_oAn x is estimated with ginseng_iBetween distance, θ is interpolation point x_oAn x is estimated with ginseng_iBetween position angle, γ (d_i, θ) and it is with d_i, θ be the variation function of variable；

S5: to described variation function incoming direction feature to obtain the variation function expression formula in all directions；

S6: judge that the spatial distribution of described pollutant is as isotropism or anisotropy according to the variation function expression formula in described all directions；

S7: based on isotropism or anisotropy setting search strategy to determine that described ginseng estimates an x_i；

S8: estimate the pollutant levels property value Z (x a little determining interpolation point in formula (2) according to the variation function expression formula in described all directions and described ginseng_o)；

S9: Research on partition region, generates multiple mesh point, and obtains the center point coordinate of each described grid；

S10: using the center point coordinate of each described grid as described interpolation point, adopts S8 that the central point of each described grid is interpolated calculating, to obtain the interpolation result of each described mesh point；

S11: adopt cross validation method that described interpolation method is verified, with mean error and root-mean-square error for index, it is thus achieved that best model；Specifically, the thinking of cross validation is: 62 websites in region, known Pearl River Delta, whole websites are begun stepping through, as taken i-th website (i=1,2 ... 62), this point is interpolated by the adaptive space interpolation algorithm using remaining 61 website property values and step S1 to S8 about, it is achieved with 62 interpolation points estimated value Z* (xi), actual property value Z (xi) in conjunction with 62 websites, it is averaging the evaluation index such as error (ME) and root-mean-square error (RMSIE), in order to verify interpolation model interpolation precision at different conditions.It should be noted that in cross validation, the pollutant levels property value Z (x of the interpolation point that interpolation point estimated value Z* (xi) is in S8_o).In short, namely suppose that a certain sample point numerical value is unknown, it is interpolated estimation by the key element value of sample point around, utilizes the error of the measured value and estimated value that calculate sample to evaluate, finally determine the best model of interpolation precision and parameter.Wherein, the calculation of mean error (ME) and root-mean-square error (RMSIE) is as follows respectively:

$\mathrm{ME}=\underset{j=1}{\overset{m}{\mathrm{\Σ}}}(C_{a,j}-C_{e,j})/m---\left(13\right)$

$\mathrm{RMSIE}=\sqrt{\frac{\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{(C_{a,j}-C_{e,j})}^2}{m}}---\left(14\right)$

Wherein, C_a,jFor the actual monitoring concentration value of jth website, i.e. actual property value Z (xi)；C_e,jFor the predictive value of jth website, i.e. interpolation point estimated value Z* (xi)；M is that authentication station is counted；

Statistics obtains SO₂Interpolation precision results contrast is as shown in table 2:

The interpolation result cross validation analysis contrast of table 2 two kinds of interpolation algorithms of different directions number

S12: the interpolation result according to each the described mesh point obtained in S10, gives different colors to different numerical value sections, and the space based on GIS platform drawing area air quality renders figure, it is achieved Visualization.

Specifically, the concrete calculation procedure of position angle θ includes:

S311: set reference vector H(100,0)；

S312: calculate two known sample point x_i,x_jVector M (the x formed_i'-x_j',y_i'-y_j'), it is designated as M (x_M',y_M'), then vector M and reference vector H institute angle degree P are then two described sample point x_i,x_jBetween position angle θ；

S313: if y_M' > 0, then have

$P=\arccos \left(\frac{H*M}{\left|H\right|*\left|M\right|}\right)*\frac{180}{\mathrm{\π}}=\arccos \left(\frac{x_M^{\′}}{\sqrt{{\left(x_M^{\′}\right)}^2+{\left(y_M^{\′}\right)}^2}}\right)*\frac{180}{\mathrm{\π}}---\left(3\right)$

If y_M' < 0, then have

$P=\arccos \left(\frac{H*M}{\left|H\right|*\left|M\right|}\right)*\frac{180}{\mathrm{\&pi;}}+180=\arccos \left(\frac{x_M^{\&prime;}}{\sqrt{{\left(x_M^{\&prime;}\right)}^2+{\left(y_M^{\&prime;}\right)}^2}}\right)*\frac{180}{\mathrm{\&pi;}}+180---\left(4\right)$

Wherein, π is pi, (x_i',y_i') for known sample point x_iLatitude and longitude coordinates, (x_j',y_j') for known sample point x_jLatitude and longitude coordinates, arccos is inverse cosine function.

It should be noted that position angle between any two sample point and interpolation point x in formula (2)_oAn x is estimated with ginseng_iBetween position angle θ step S311 to S313 all can be adopted to be calculated.

Preferably, needing variation function is fitted before carrying out step S5, its concrete fit procedure is:

1) the variation function value under different distance is calculated according to formula (7) to formula (9)；

$\mathrm{\&gamma;}\left(h\right)=\mathrm{\&gamma;}\left(\right|x_i-x_j\left|\right)=\frac{1}{2}E{[Z\left(x_i\right)-Z\left(x_j\right)]}^2---\left(8\right)$

${\mathrm{\&gamma;}}^{*}\left(h\right)-{\mathrm{\&gamma;}}^{*}\left(\right|x_i-x_j\left|\right)-\frac{1}{2m}\underset{k=1}{\overset{m}{\mathrm{\&Sigma;}}}{[Z\left(x_i\right)-Z\left(x_j\right)]}^2---\left(9\right)$

Wherein γ (h) is desirable variation function, γ^*H () is discrete variation function, x_i,x_jFor the sample point pair of distance h, Z (x_i)、Z(x_j) for sample point x_i,x_jProperty value, h is sample point x_i,x_jBetween distance, (x_i',y_i') for known sample point x_iLatitude and longitude coordinates, (x_j',y_j') for known sample point x_jLatitude and longitude coordinates, m represents the number of the sample point pair that distance is h, and E represents expectation；

2) with distance h for abscissa, variation function γ^*H the value of () is vertical coordinate, draw variation function scatterplot, as shown in Figure 2；

3) after above-mentioned variation function 3 being normalized, choose spherical model, exponential model or Gauss model to variation function scatterplot matching to obtain the highest model of degree of fitting and function expression γ (h) thereof, in the present embodiment, spherical model is the model that degree of fitting is the highest；The fitting effect of three kinds of models is as it is shown on figure 3, the point in figure represents variation function, and 1 represents exponential model, and 2 represent spherical model, and 3 represent Gauss model；Wherein, the function expression of spherical model, exponential model and Gauss model is respectively as shown in formula (10), (11), (12):

$\mathrm{\&gamma;}\left(h\right)=\left\{\begin{array}{cc}0& \\ C_o+C(\frac{3h}{2a}-\frac{h^3}{{2a}^3})& (0<h\&le;a)\\ C_o+C& (h>a)\end{array}\right.---\left(10\right)$

$\mathrm{\&gamma;}\left(h\right)=C_0+C(1-e^{\frac{h}{a}})---\left(11\right)$

$\mathrm{\&gamma;}\left(h\right)=C_0+C(1-e\frac{h^2}{a^2})---\left(12\right)$

Wherein, a, C_o, C, e be constant.

Specifically, as shown in Figure 4, step S5 specifically includes:

S51: selected directions number, each direction is a direction group, and orientation angle corresponding to each direction isWherein p=1 ..., N；

S52: determine tolerance angle, and calculate the angular range of each direction group；About tolerance angle Δ θ, it is set asSpecifically, choose 8 directions, i.e. N=8, then there are 8 direction groups, orientation angle corresponding to eight directions respectively 0, 45, 90, 135, 180, 225, 270, 315 degree, and it is numbered 1, 2, 3, 4, 5, 6, 7, 8 direction groups, determine that tolerance angle is 22.5 degree, then think that x ± 22.5 degree broadly fall into x direction, as 45 ± 22.5 degree broadly fall into 45 degree of directions, now 0, 45, 90, 135, 180, 225, 270, the angular range of 315 degree of 8 direction groups is respectively as follows: 0 ± 22.5 degree, 0-22.5+360 degree, 45 ± 22.5 degree, 90 ± 22.5 degree, 135 ± 22.5 degree, 180 ± 22.5 degree, 225 ± 22.5 degree, 270 ± 22.5 degree, 315 ± 22.5 degree；

S53: all of sample point is numbered, selects any two sample point to form sample point pair, calculates the distance of described sample point pair, it is determined that tolerance spacing, and computed range scope；The distance of described sample point pair is calculated according to formula (7), and tolerates that the setting of separation delta h is according to direction-adaptive, is set to the half of the minimum range of sample point pair corresponding to each direction；

S54: calculate the position angle between two sample points according to step S311 to S313, and according to position angle by sample point to being divided in the group of corresponding direction；Specifically, by all sample points according to locus, from left to right arrange along X-axis and be numbered 1～62, first calculate sample point 1 and sample point 2,3,4,5 ... 62 position angles being, according to sample point 1 and sample point 2,3,4,5 ... 62 position angles being by this point to (sample point 1 and sample point 2, sample point 1 and sample point 3 ... sample point 1 and sample point 62) be divided to corresponding direction group, travel through all sample points successively, finally give the sample point pair of all directions group；

S55: according to the described sample point in direction group each described to and formula (3) to formula (7) calculate the variation function expression formula in each described direction group, wherein, formula (5) to (7) is respectively as follows:

$\mathrm{\&gamma;}(h,\mathrm{\&theta;})-\mathrm{\&gamma;}\left(\right|x_i-x_j|,\mathrm{\&theta;})-\frac{1}{2}E{[Z\left(x_i\right)-Z\left(x_j\right)|\mathrm{\&theta;}]}^2---\left(5\right)$

${\mathrm{\&gamma;}}^{*}(h,\mathrm{\&theta;})={\mathrm{\&gamma;}}^{*}\left(\right|x_i-x_j|,\mathrm{\&theta;})=\frac{1}{2m}\underset{k=1}{\overset{m}{\mathrm{\&Sigma;}}}{[Z\left(x_i\right)-Z\left(x_j\right)|\mathrm{\&theta;}]}^2---\left(6\right)$

$h=\sqrt{{({x_i}^{\&prime;}-{x_j}^{\&prime;})}^2+{({y_i}^{\&prime;}-{y_j}^{\&prime;})}^2}---\left(7\right)$

Wherein γ (h, θ) is desirable variation function, γ^*(h, θ) is discrete variation function, x_i,x_jFor the sample point pair of distance h, Z (x_i)、Z(x_j) for sample point x_i,x_jProperty value, h is that sample point is to x_i,x_jBetween distance, (x_i',y_i') for known sample point x_iLatitude and longitude coordinates, (x_j',y_j') for known sample point x_jLatitude and longitude coordinates, θ is described sample point x_i,x_jBetween position angle, m represents the number of the sample point pair that distance is h；Specifically, SO is obtained using spherical model as the theoretical fitting model of variation function₂Variation function expression formula (spherical model) upwards is as follows respectively from all directions:

1) 0 degree/180 degree

$\mathrm{\&gamma;}\left(h\right)=\left\{\begin{array}{cc}0& h=0\\ 3.1+0.339h-5.065*{10}^{-5}{h}^3& 0<h\&le;47.3\\ 13.8& h>47.3\end{array}\right.$

2) 45 degree/225 degree

$\mathrm{\&gamma;}\left(h\right)=\left\{\begin{array}{cc}0& h=0\\ 3.51+0.247h-3.685*{10}^{-5}{h}^3& 0<h\&le;47.3\\ 11.4& h>47.3\end{array}\right.$

3) 90 degree/270 degree

$\mathrm{\&gamma;}\left(h\right)=\left\{\begin{array}{cc}0& h=0\\ 2.67+0.492h-8.794*{10}^{-5}{h}^3& 0<h\&le;43.2\\ 16.85& h>43.2\end{array}\right.$

4) 135/315 degree

$\mathrm{\&gamma;}\left(h\right)=\left\{\begin{array}{cc}0& h=0\\ 5.72+0.366h-6.143*{10}^{-5}{h}^3& 0<h\&le;44.6\\ 16.62& h>44.6\end{array}\right.$

With distance h for abscissa, the value of variation function γ (h) is vertical coordinate, SO₂Variation function scatterplot upwards is as shown in Figure 5 from all directions.

Specifically, refer to Fig. 6, step S6 specifically includes:

S61: calculate range value according to the variation function in described all directions, and draw rose diagram according to described range value；Specifically, according to the variation function expression formula on eight directions of step S55, it may be determined that go out the range in all directions in spherical model, as shown in table 1:

Direction	0°/180°	45°/225°	90°/270°	135°/315°
					Range/km	47.3	47.3	43.2	44.6

Range value list (the SO that region, table 1 Pearl River Delta is downward from all directions₂)

The rose diagram drawn is as shown in Figure 7；

S62: judge that the spatial distribution of described pollutant is as isotropism or anisotropy according to described rose diagram；Specifically, as it is shown in fig. 7, left figure is SO₂Rose diagram, in the characteristic of circle, therefore, it is determined that SO₂Spatial distribution there is isotropism, and the PM of right figure_2.5Rose diagram in oval feature, therefore, it is determined that PM_2.5Spatial distribution there is anisotropy.

Refer to Fig. 8 again, step S7 specifically includes:

S71: set minimum described ginseng and estimate a number；

S72: when the spatial distribution of described pollutant is isotropism, with the average range value of the variation function in all directions be radius, described interpolation point draw one circular for the center of circle, an x is estimated for searching for described ginseng in described circular inner zone_iRegion of search；

S73: determine whether that searching described ginseng estimates an x_i, if so, then terminate, otherwise, then perform S74；

S74: estimate an x when not searching described ginseng in S72_iTime, minimum range and described average range value ratio according to whole sample points pair expand described radius by the mode changed again, until searching out minimum described ginseng set in S71 to estimate a number；

S75: when the spatial distribution of described pollutant is anisotropy, it is one oval that the major semiaxis a simulated according to the range value of the variation function in all directions, semi-minor axis b, described interpolation point are that axle center is drawn, and described oval interior zone estimates an x for searching for described ginseng_iRegion of search；

S76: judge that whether searching described ginseng in S75 estimates an x_i, if so, then terminate, otherwise, then perform S77；

S77: estimate an x when not searching described ginseng in S75_iTime, according to whole sample points pair minimum range withRatio expand described major semiaxis and semi-minor axis by the mode changed again, until searching out in S71 set minimum described ginseng to estimate a number.

As can be seen from the above description, the present invention first extracts geographical location information and pollutant monitoring value based on the adaptive space interpolation method of spatial distribution characteristic, again based on the spatial distribution characteristic of Geostatistical principle analysis regional environment air quality, it is then based on variation function and sets up adaptive space interpolation model, incoming direction feature is to adaptive space interpolation model and variation function afterwards, judge that the spatial distribution of pollutant is as isotropism or anisotropy further according to the variation function of incoming direction feature, and based on isotropism and anisotropy setting search strategy to determine that ginseng is estimated a little, final variation function and ginseng according to all directions estimates the pollution concentration property value a little determining interpolation point；This interpolation method is the novel interpolation method considering spatial distribution characteristic generated after sample point attribute data is made further space characteristics analysis, under considering the premise of spatial distribution characteristic of sample point attribute data data, ensure that amount of calculation is little, and improve interpolation precision, thus the attribute character value of predictive study region unknown point more exactly.

It should be noted that the interpolation method of the present invention is for Kriging method, amount of calculation is to reduce, because Kriging method relates to solving an equation, often calculates a weight coefficient it is necessary to solve multiple equation, and is not related to solving equation group herein, therefore greatly reduce workload.

Correspondingly, present invention also offers a kind of adaptive space interplotation system based on spatial distribution characteristic to include:

Extraction module 10, is used for extracting geographical location information and pollutant monitoring value；

Analysis module 11, for based on Geostatistical principle, the spatial distribution characteristic of analyzed area ambient air quality；

Model building module 12, for setting up adaptive space interpolation model based on variation function, the expression formula of described adaptive space interpolation model is:

$Z\left(x_0\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_{i}Z\left(x_i\right)=\frac{{\mathrm{\&Sigma;}}_{i=1}^n\frac{1}{\mathrm{\&gamma;}\left(d_i\right)}Z\left(x_i\right)}{{\mathrm{\&Sigma;}}_{i=1}^n\frac{1}{\mathrm{\&gamma;}\left(d_i\right)}}---\left(1\right)$

Wherein, λ_iFor known sample point x_iWeight coefficient；Z (x_o) for the pollutant levels property value of interpolation point, i is natural number, and n is the sample point sum affecting interpolation point, d_iFor interpolation point x_oAn x is estimated with ginseng_iBetween distance, Z (x_i) estimate pollutant levels property value a little, γ (d for ginseng_i) it is with d_iVariation function for variable；

Specifically, the expression formula of variation function is:

Nondirectional ideal expression

$\mathrm{\&gamma;}\left(h\right)=\mathrm{\&gamma;}\left(\right|x_i-x_j\left|\right)=\frac{1}{2}E{[Z\left(x_i\right)-Z\left(x_j\right)]}^2---\left(8\right)$

Nondirectional discrete expression

${\mathrm{\&gamma;}}^{*}\left(h\right)={\mathrm{\&gamma;}}^{*}\left(\right|x_i-x_j\left|\right)=\frac{1}{2m}\underset{k=1}{\overset{m}{\mathrm{\&Sigma;}}}{[Z\left(x_i\right)-Z\left(x_j\right)]}^2---\left(9\right)$

Wherein γ (h) is desirable variation function, γ^*H () is discrete variation function, x_i,x_jFor the sample point pair of distance h, Z (x_i)、Z(x_j) for sample point x_i,x_jProperty value, h is sample point x_i,x_jBetween distance, (x_i',y_i') for known sample point x_iLatitude and longitude coordinates, (x_j',y_j') for known sample point x_jLatitude and longitude coordinates, m represents the number of the sample point pair that distance is h；

Distance h adopts following formula to be calculated:

$h=\sqrt{{({x_i}^{\&prime;}-{x_j}^{\&prime;})}^2+{({y_i}^{\&prime;}-{y_j}^{\&prime;})}^2}---\left(7\right)$

Wherein, (x_i',y_i') for known sample point x_iLatitude and longitude coordinates, (x_j',y_j') for known sample point x_jLatitude and longitude coordinates.

First introduces module 13, obtains formula (2) for incoming direction feature to described adaptive space interpolation model:

$Z\left(x_o\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_{i}Z\left(x_i\right)=\frac{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\frac{1}{\mathrm{\&gamma;}(d_i,\mathrm{\&theta;})}Z\left(x_i\right)}{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\frac{1}{\mathrm{\&gamma;}(d_i,\mathrm{\&theta;})}}---\left(2\right)$

Wherein, λ_iFor known sample point x_iWeight coefficient；Z (x_o) for the pollutant levels property value of interpolation point, i is natural number, and n is the sample point sum affecting interpolation point, d_iFor interpolation point x_oAn x is estimated with ginseng_iBetween distance, θ is interpolation point x_oAn x is estimated with ginseng_iBetween position angle, γ (d_i, θ) and it is with d_i, θ be the variation function of variable；

Fitting module 14, for being fitted variation function；

Second introduce module 15, for described variation function incoming direction feature to obtain the variation function expression formula in all directions；

Determination module 16, for judging that according to the variation function expression formula in described all directions the spatial distribution of described pollutant is as isotropism or anisotropy；

Search module 17, for based on isotropism or anisotropy setting search strategy to determine that described ginseng estimates an x_i；

Determine module 18, for estimating the pollutant levels property value Z (x a little determining interpolation point in formula (2) according to the variation function expression formula in described all directions and described ginseng_o)；

Divide acquisition module 19, for Research on partition region, generate multiple mesh point and obtain the center point coordinate of each described grid；

Interpolating module 20, for using the center point coordinate of each described grid as described interpolation point, adopt and described determine that the central point of each described grid is interpolated by module and calculate interpolation result to obtain each described mesh point；

Authentication module 21, is used for adopting cross validation method that described interpolation method is verified, obtaining best model with mean error and root-mean-square error for index；Specifically, the thinking of cross validation is: 62 websites in region, known Pearl River Delta, whole websites are begun stepping through, as taken i-th sample point (i=1,2 ... 62), this point is interpolated by the adaptive space interpolation algorithm using remaining 61 website property values and step S1 to S8 about, it is achieved with 62 interpolation points estimated value Z* (xi), actual property value Z (xi) in conjunction with 62 websites, it is averaging the evaluation index such as error (ME) and root-mean-square error (RMSIE), in order to verify interpolation model interpolation precision at different conditions.It should be noted that in cross validation, the pollutant levels property value Z (x of the interpolation point that interpolation point estimated value Z* (xi) is in S_o).In short, namely suppose that a certain sample point numerical value is unknown, it is interpolated estimation by the key element value of sample point around, utilizes the error of the measured value and estimated value that calculate sample to evaluate, finally determine the best model of interpolation precision and parameter.Wherein, the calculation of mean error (ME) and root-mean-square error (RMSIE) is as follows respectively:

$\mathrm{ME}=\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}(C_{a,j}-C_{e,j})/m---\left(13\right)$

$\mathrm{RMSIE}=\sqrt{\frac{\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}{(C_{a,j}-C_{e,j})}^2}{m}}---\left(14\right)$

Wherein, C_a,jFor the actual monitoring concentration value of jth website, i.e. actual property value Z (xi)；C_e,jFor the predictive value of jth website, i.e. interpolation point estimated value Z* (xi)；M is that authentication station is counted, and is m=62 in this；

Statistics obtains SO₂Interpolation precision results contrast is as shown in table 2:

The interpolation result cross validation analysis contrast of table 2 two kinds of interpolation algorithms of different directions number

Rendering module 22, the interpolation result of each the described mesh point for obtaining according to described interpolating module, different numerical value sections is given different color, space based on GIS platform drawing area air quality renders figure, it is achieved that Visualization.

Specifically, as shown in Figure 10, described fitting module 14 specifically includes:

Computing unit 141, is used for calculating the variation function value under different distance according to formula (7) with to formula (9)；

$\mathrm{\&gamma;}\left(h\right)=\mathrm{\&gamma;}\left(\right|x_i-x_j\left|\right)=\frac{1}{2}E{[Z\left(x_i\right)-Z\left(x_j\right)]}^2---\left(8\right)$

${\mathrm{\&gamma;}}^{*}\left(h\right)={\mathrm{\&gamma;}}^{*}\left(\right|x_i-x_j\left|\right)=\frac{1}{2m}\underset{k=1}{\overset{m}{\mathrm{\&Sigma;}}}{[Z\left(x_i\right)-Z\left(x_j\right)]}^2---\left(9\right)$

Wherein γ (h) is desirable variation function, γ^*H () is discrete variation function, x_i,x_jFor the sample point pair of distance h, Z (x_i)、Z(x_j) for sample point x_i,x_jProperty value, h is sample point x_i,x_jBetween distance, (x_i',y_i') for known sample point x_iLatitude and longitude coordinates, (x_j',y_j') for known sample point x_jLatitude and longitude coordinates, k is natural number, and m represents the number of the sample point pair that distance is h；

Drawing unit 142, is used for distance h for abscissa, variation function value γ^*H () is vertical coordinate, draw variation function scatterplot；

Fitting unit 143, for choosing spherical model, exponential model or Gauss model to variation function scatterplot matching to obtain the highest model of degree of fitting and function expression γ (h) thereof.

Specifically, as shown in figure 11, described second introducing module 15 specifically includes:

First module 151, for selected directions number N, each direction is a direction group, and orientation angle corresponding to each direction isWherein p=1 ..., N；

Second unit 152, for determining tolerance angle and calculating the angular range of each direction group；

3rd unit 153, for all of sample point is numbered, selects sample point described in any two with is formed sample point to, calculate the distance of described sample point pair, determine and tolerate spacing computed range scope；

4th unit 154, for calculate position angle between two described sample points and according to described position angle by described sample point to being divided in the group of corresponding described direction；

5th unit 155, for according to the described sample point in direction group each described to and formula (3) to formula (7) calculate the variation function in each described direction group；Wherein, formula (3) is respectively as follows: to formula (7)

If y_M' > 0, then have

$P=\arccos \left(\frac{H*M}{\left|H\right|*\left|M\right|}\right)*\frac{180}{\mathrm{\&pi;}}=\arccos \left(\frac{x_M^{\&prime;}}{\sqrt{{\left(x_M^{\&prime;}\right)}^2+{\left(y_M^{\&prime;}\right)}^2}}\right)*\frac{180}{\mathrm{\&pi;}}---\left(3\right)$

If y_M' < 0, then have

$P=\arccos \left(\frac{H*M}{\left|H\right|*\left|M\right|}\right)*\frac{180}{\mathrm{\&pi;}}+180=\arccos \left(\frac{x_M^{\&prime;}}{\sqrt{{\left(x_M^{\&prime;}\right)}^2+{\left(y_M^{\&prime;}\right)}^2}}\right)*\frac{180}{\mathrm{\&pi;}}+180---\left(4\right)$

$\mathrm{\&gamma;}(h,\mathrm{\&theta;})=\mathrm{\&gamma;}\left(\right|x_i-x_j|,\mathrm{\&theta;})=\frac{1}{2}E{[Z\left(x_i\right)-Z\left(x_j\right)|\mathrm{\&theta;}]}^2---\left(5\right)$

${\mathrm{\&gamma;}}^{*}(h,\mathrm{\&theta;})={\mathrm{\&gamma;}}^{*}\left(\right|x_i-x_j|,\mathrm{\&theta;})=\frac{1}{2m}\underset{k=1}{\overset{m}{\mathrm{\&Sigma;}}}{[Z\left(x_i\right)-Z\left(x_j\right)|\mathrm{\&theta;}]}^2---\left(6\right)$

$h=\sqrt{{({x_i}^{\&prime;}-{x_j}^{\&prime;})}^2+{({y_i}^{\&prime;}-{y_j}^{\&prime;})}^2}---\left(7\right)$

Wherein, H(100,0) for reference vector, M (x_M',y_M') it is two described sample point x_i,x_jVector M (the x formed_i'-x_j',y_i'-y_j'), P is vector M and reference vector H institute angle degree, and P is sample point x described in two_i,x_jBetween position angle θ, π be pi, γ (h, θ) is desirable variation function, γ^*(h, θ) is discrete variation function, x_i,x_jFor the sample point pair of distance h, Z (x_i)、Z(x_j) for sample point x_i,x_jProperty value, h is that sample point is to x_i,x_jBetween distance, (x_i',y_i') for known sample point x_iLatitude and longitude coordinates, (x_j',y_j') for known sample point x_jLatitude and longitude coordinates, θ is described sample point x_i,x_jBetween position angle, k is natural number, and m represents the number of the sample point pair that distance is h.

Specifically, as shown in figure 12, described determination module 16 specifically includes:

First identifying unit 161, for calculating range value according to the variation function in described all directions, and draws rose diagram according to described range value；

According to described rose diagram, second identifying unit 162, for judging that the spatial distribution of described pollutant is as isotropism or anisotropy.Specifically, when rose diagram is in the characteristic of circle, it is determined that the spatial distribution of pollutant has an isotropism, and rose diagram in oval feature time, then judge that the spatial distribution of pollutant has anisotropy.

Specifically, as shown in figure 13, described search module 17 specifically includes:

Setup unit 171, estimates a number for setting minimum described ginseng；

First search unit 172, estimates a little for searching for described ginseng when the spatial distribution of described pollutant is isotropism；Specifically, when the spatial distribution of described pollutant is isotropism, with the average range value of the variation function in all directions be radius, described interpolation point draw one circular for the center of circle, an x is estimated for searching for described ginseng in described circular inner zone_iRegion of search；Determine whether that searching described ginseng estimates an x_i；An x is estimated when not searching described ginseng_iTime, minimum range and described average range value ratio according to whole sample points pair expand described radius by the mode changed again, until searching out minimum described ginseng set in setup unit 171 to estimate a number

Second search unit 173, estimates a little for searching for described ginseng when the spatial distribution of described pollutant is anisotropy.Specifically, when the spatial distribution of described pollutant is anisotropy, simulating major semiaxis a according to the range value of the variation function in all directions, that semi-minor axis b, described interpolation point are that axle center is drawn is one oval, described oval interior zone estimates an x for searching for described ginseng_iRegion of search；Determine whether that searching described ginseng estimates an x_i；An x is estimated when not searching described ginseng_iTime, according to whole sample points pair minimum range withRatio expand described major semiaxis and semi-minor axis by the mode changed again, until searching out in setup unit 171 set minimum described ginseng to estimate a number.

As can be seen from the above description, the present invention is based on the adaptive space interplotation system of spatial distribution characteristic, first pass through extraction module 10 and extract geographical location information and pollutant monitoring value, again through the analysis module 11 spatial distribution characteristic based on Geostatistical principle analysis regional environment air quality, then pass through model building module 12 and set up adaptive space interpolation model based on variation function, afterwards by the first introducing module 13, second introduces module 15 incoming direction feature to adaptive space interpolation model and variation function, judge that the spatial distribution of pollutant is as isotropism or anisotropy again through determination module 16 according to the variation function of incoming direction feature, and by search module 17 based on isotropism and anisotropy setting search strategy to determine that ginseng is estimated a little, eventually through determining that module 18 estimates, according to variation function and the ginseng of all directions, the pollution concentration property value a little determining interpolation point；This interplotation system is the novel interpolation method considering spatial distribution characteristic generated after sample point attribute data is made further space characteristics analysis, under considering the premise of spatial distribution characteristic of sample point attribute data data, ensure that amount of calculation is little, and improve interpolation precision, thus the attribute character value of predictive study region unknown point more exactly.

It should be noted that the interpolation method of the present invention is for Kriging method, amount of calculation is to reduce, because Kriging method relates to solving an equation, often calculates a weight coefficient it is necessary to solve multiple equation, is not related to solving equation group herein, greatly reduces workload.

Above in association with most preferred embodiment, invention has been described, but the invention is not limited in embodiment disclosed above, and should contain amendment, the equivalent combinations that the various essence according to the present invention carries out.

Claims (17)

1. the adaptive space interpolation method based on space characteristics analysis, it is characterised in that include step:

S1: extract geographical location information and pollutant monitoring value；

S2: based on Geostatistical principle, the spatial distribution characteristic of analyzed area ambient air quality；

S3: set up adaptive space interpolation model based on variation function, the expression formula of described adaptive space interpolation model is:

$Z\left(x_0\right)=\underset{i=1}{\overset{n}{\&Sigma;}}{\mathrm{\&lambda;}}_{i}Z\left(x_i\right)=\frac{{\mathrm{\&Sigma;}}_{i=1}^n\frac{1}{\mathrm{\&gamma;}\left(d_i\right)}Z\left(x_i\right)}{{\mathrm{\&Sigma;}}_{i=1}^n\frac{1}{\mathrm{\&gamma;}\left(d_i\right)}}---\left(1\right)$

Wherein, λ_iFor known sample point x_iWeight coefficient；Z (x_o) for the pollutant levels property value of interpolation point, i is natural number, and n is the sample point sum affecting interpolation point, d_iFor interpolation point x_oAn x is estimated with ginseng_iBetween distance, Z (x_i) estimate pollutant levels property value a little, γ (d for ginseng_i) it is with d_iVariation function for variable；

S4: incoming direction feature to described adaptive space interpolation model obtains formula (2):

$Z\left(x_o\right)=\underset{i=1}{\overset{n}{\&Sigma;}}{\mathrm{\&lambda;}}_{i}Z\left(x_i\right)=\frac{\underset{i=1}{\overset{n}{\&Sigma;}}\frac{1}{\mathrm{\&gamma;}(d_i,\mathrm{\&theta;})}Z\left(x_i\right)}{\underset{i=1}{\overset{n}{\&Sigma;}}\frac{1}{\mathrm{\&gamma;}(d_i,\mathrm{\&theta;})}}---\left(2\right)$

Wherein, λ_iFor known sample point x_iWeight coefficient；Z (x_o) for the pollutant levels property value of interpolation point, i is natural number, and n is the sample point sum affecting interpolation point, d_iFor interpolation point x_oAn x is estimated with ginseng_iBetween distance, θ is interpolation point x_oAn x is estimated with ginseng_iBetween position angle, γ (d_i, θ) and it is with d_i, θ be the variation function of variable；

S5: to described variation function incoming direction feature to obtain the variation function expression formula in all directions；

S6: judge that the spatial distribution of described pollutant is as isotropism or anisotropy according to the variation function expression formula in described all directions；

S7: based on isotropism or anisotropy setting search strategy to determine that described ginseng estimates an x_i；

S8: estimate the pollutant levels property value Z (x a little determining interpolation point in formula (2) according to the variation function expression formula in described all directions and described ginseng_o)。

2. the adaptive space interpolation method based on space characteristics analysis as claimed in claim 1, it is characterised in that the concrete calculation procedure of position angle θ includes:

S311: set reference vector H (100,0)；

S312: calculate two known sample point x_i,x_jVector M (the x formed_i'-x_j',y_i'-y_j'), it is designated as M (x_M',y_M'), then vector M and reference vector H institute angle degree P are then two described sample point x_i,x_jBetween position angle θ；

S313: if y_M' > 0, then have

$P=\arccos \left(\frac{H*M}{\left|H\right|*\left|M\right|}\right)*\frac{180}{\mathrm{\&pi;}}=\arccos \left(\frac{x_M^{\&prime;}}{\sqrt{{\left(x_M^{\&prime;}\right)}^2+{\left(y_M^{\&prime;}\right)}^2}}\right)*\frac{180}{\mathrm{\&pi;}}---\left(3\right)$

If y_M' < 0, then have

$P=\arccos \left(\frac{H*M}{\left|H\right|*\left|M\right|}\right)*\frac{180}{\mathrm{\&pi;}}+180=\arccos \left(\frac{x_M^{\&prime;}}{\sqrt{{\left(x_M^{\&prime;}\right)}^2+{\left(y_M^{\&prime;}\right)}^2}}\right)*\frac{180}{\mathrm{\&pi;}}+180---\left(4\right)$

Wherein, π is pi, (x_i',y_i') for known sample point x_iLatitude and longitude coordinates, (x_j',y_j') for known sample point x_jLatitude and longitude coordinates.

3. the adaptive space interpolation method based on space characteristics analysis as claimed in claim 2, it is characterised in that step S5 specifically includes:

S51: selected directions number N, each direction is a direction group, and orientation angle corresponding to each direction isWherein p=1 ..., N；

S52: determine tolerance angle, and calculate the angular range of each described direction group；

S53: all of sample point is numbered, selects sample point described in any two to form sample point pair, to calculate the distance of described sample point pair, it is determined that tolerance spacing, and computed range scope；

S54: calculate the position angle between two described sample points according to step S311 to S313, and according to described position angle by described sample point to being divided in the group of corresponding described direction；

S55: according to the described sample point in direction group each described to and formula (3) to formula (7) calculate the variation function expression formula in each described direction group, wherein, formula (5) to (7) is respectively as follows:

$\mathrm{\&gamma;}(h,\mathrm{\&theta;})=\mathrm{\&gamma;}\left(\right|x_i-x_j|,\mathrm{\&theta;})=\frac{1}{2}E{\&lsqb;Z\left(x_i\right)-Z\left(x_j\right)|\mathrm{\&theta;}\&rsqb;}^2---\left(5\right)$

${\mathrm{\&gamma;}}^{*}\left(h,\mathrm{\&theta;}\right)={\mathrm{\&gamma;}}^{*}\left(|x_i-x_j|,\mathrm{\&theta;}\right)=\frac{1}{2m}{\mathrm{\&Sigma;}}_{k=1}^m{\&lsqb;Z\left(x_i\right)-Z\left(x_j\right)|\mathrm{\&theta;}\&rsqb;}^2---\left(6\right)$

$h=\sqrt{{({x_i}^{\&prime;}-{x_j}^{\&prime;})}^2+{({y_i}^{\&prime;}-{y_j}^{\&prime;})}^2}---\left(7\right)$

Wherein γ (h, θ) is desirable variation function, and γ * (h, θ) is discrete variation function, x_i,x_jThe sample point pair being h for distance, Z (x_i)、Z(x_j) for sample point x_i,x_jProperty value, h is that sample point is to x_i,x_jBetween distance, (x_i',y_i') for known sample point x_iLatitude and longitude coordinates, (x_j',y_j') for known sample point x_jLatitude and longitude coordinates, θ is described sample point x_i,x_jBetween position angle, m represents the number of the sample point pair that distance is h.

4. the adaptive space interpolation method based on space characteristics analysis as claimed in claim 3, it is characterised in that also include before step S5:

1) the variation function value under different distance is calculated according to formula (7) to formula (9)；

$\mathrm{\&gamma;}\left(h\right)=\mathrm{\&gamma;}\left(\right|x_i-x_j\left|\right)=\frac{1}{2}E{\&lsqb;Z\left(x_i\right)-Z\left(x_j\right)\&rsqb;}^2---\left(8\right)$

${\mathrm{\&gamma;}}^{*}\left(h\right)={\mathrm{\&gamma;}}^{*}\left(\right|x_i-x_j\left|\right)=\frac{1}{2m}{\mathrm{\&Sigma;}}_{k=1}^m{\&lsqb;Z\left(x_i\right)-Z\left(x_j\right)\&rsqb;}^2---\left(9\right)$

Wherein γ (h) is desirable variation function, γ^*H () is discrete variation function, x_i,x_jThe sample point pair being h for distance, Z (x_i)、Z(x_j) for sample point x_i,x_jProperty value, h is sample point x_i,x_jBetween distance, (x_i',y_i') for known sample point x_iLatitude and longitude coordinates, (x_j',y_j') for known sample point x_jLatitude and longitude coordinates, m represents the number of the sample point pair that distance is h；

2) with distance h for abscissa, variation function γ^*H the value of () is vertical coordinate, draw variation function scatterplot；

3) spherical model, exponential model or Gauss model are chosen to variation function scatterplot matching to obtain the highest model of degree of fitting and function expression γ (h) thereof.

5. the adaptive space interpolation method based on space characteristics analysis as claimed in claim 4, it is characterised in that step S6 specifically includes:

S61: calculate range value according to the variation function in described all directions, and draw rose diagram according to described range value；

S62: judge that the spatial distribution of described pollutant is as isotropism or anisotropy according to described rose diagram.

6. the adaptive space interpolation method based on space characteristics analysis as claimed in claim 5, it is characterised in that step S7 specifically includes:

S71: set minimum described ginseng and estimate a number；

S72: when the spatial distribution of described pollutant is isotropism, with the average range value of the variation function in all directions be radius, described interpolation point draw one circular for the center of circle, described circular inner zone is the region of search that described ginseng estimates a little；

S73: judge that whether searching described ginseng in S72 estimates a little；

A S74: when not searching described ginseng in S72 and estimating, minimum range and described average range value ratio according to whole sample points pair expand described radius by the mode changed again, until searching out minimum described ginseng set in S71 to estimate a number；

S75: when the spatial distribution of described pollutant is anisotropy, it is one oval that major semiaxis a that range value according to the variation function in all directions simulates, semi-minor axis b, described interpolation point are that axle center is drawn, and described oval interior zone is search for described ginseng to estimate region of search a little；

S76: judge that whether searching described ginseng in S75 estimates a little；

A S77: when not searching described ginseng in S75 and estimating, according to whole sample points pair minimum range withRatio expand described major semiaxis and semi-minor axis by the mode changed again, until searching out in S71 set minimum described ginseng to estimate a number.

7. the adaptive space interpolation method based on space characteristics analysis as claimed in claim 6, it is characterised in that also include:

S9: Research on partition region, generates multiple mesh point, and obtains the center point coordinate of each described grid；

S10: using the central point of each described grid as described interpolation point, adopts S8 that the central point of each described grid is interpolated calculating, to obtain the interpolation result of each described mesh point.

8. the adaptive space interpolation method based on space characteristics analysis as claimed in claim 7, it is characterised in that also include:

S11: adopt cross validation method that described interpolation method is verified.

9. as claimed in claim 7 or 8 based on the adaptive space interpolation method of space characteristics analysis, it is characterised in that also include:

S12: the interpolation result according to each the described mesh point obtained in S10, gives different colors to different numerical value sections, and the space based on GIS platform drawing area air quality renders figure.

10. the adaptive space interplotation system based on space characteristics analysis, it is characterised in that including:

Extraction module, is used for extracting geographical location information and pollutant monitoring value；

Analysis module, for based on Geostatistical principle, the spatial distribution characteristic of analyzed area ambient air quality；

Model building module, for setting up adaptive space interpolation model based on variation function, the expression formula of described adaptive space interpolation model is:

$Z\left(x_0\right)=\underset{i=1}{\overset{n}{\&Sigma;}}{\mathrm{\&lambda;}}_{i}Z\left(x_i\right)=\frac{{\mathrm{\&Sigma;}}_{i=1}^n\frac{1}{\mathrm{\&gamma;}\left(d_i\right)}Z\left(x_i\right)}{{\mathrm{\&Sigma;}}_{i=1}^n\frac{1}{\mathrm{\&gamma;}\left(d_i\right)}}---\left(1\right)$

Wherein, λ_iFor known sample point x_iWeight coefficient；Z (x_o) for the pollutant levels property value of interpolation point, i is natural number, and n is the sample point sum affecting interpolation point, d_iFor interpolation point x_oAn x is estimated with ginseng_iBetween distance, Z (x_i) estimate pollutant levels property value a little, γ (d for ginseng_i) it is with d_iVariation function for variable；

First introduces module, obtains formula (2) for incoming direction feature to described adaptive space interpolation model:

$Z\left(x_o\right)=\underset{i=1}{\overset{n}{\&Sigma;}}{\mathrm{\&lambda;}}_{i}Z\left(x_i\right)=\frac{\underset{i=1}{\overset{n}{\&Sigma;}}\frac{1}{\mathrm{\&gamma;}(d_i,\mathrm{\&theta;})}Z\left(x_i\right)}{\underset{i=1}{\overset{n}{\&Sigma;}}\frac{1}{\mathrm{\&gamma;}(d_i,\mathrm{\&theta;})}}---\left(2\right)$

Wherein, λ_iFor known sample point x_iWeight coefficient；Z (x_o) for the pollutant levels property value of interpolation point, i is natural number, and n is the sample point sum affecting interpolation point, d_iFor interpolation point x_oAn x is estimated with ginseng_iBetween distance, θ is interpolation point x_oAn x is estimated with ginseng_iBetween position angle, γ (d_i, θ) and it is with d_i, θ be the variation function of variable；

Second introduce module, for described variation function incoming direction feature to obtain the variation function expression formula in all directions；

Determination module, for judging that according to the variation function expression formula in described all directions the spatial distribution of described pollutant is as isotropism or anisotropy；

Search module, for based on isotropism or anisotropy setting search strategy to determine that described ginseng estimates an x_i；

Determine module, for estimating the pollutant levels property value Z (x a little determining interpolation point in formula (2) according to the variation function expression formula in described all directions and described ginseng_o)。

11. the adaptive space interplotation system based on space characteristics analysis as claimed in claim 10, it is characterised in that described second introduces module specifically includes:

First module, for selected directions number N, each direction is a direction group, and orientation angle corresponding to each direction isWherein p=1 ..., N；

Second unit, for determining tolerance angle and calculating the angular range of each direction group；

Unit the 3rd, for all of sample point is numbered, selects sample point described in any two with is formed sample point to, calculate the distance of described sample point pair, determine and tolerate spacing computed range scope；

Unit the 4th, for calculate position angle between two described sample points and according to described position angle by described sample point to being divided in the group of corresponding described direction；

Unit the 5th, for according to the described sample point in direction group each described to and formula (3) to formula (7) calculate the variation function in each described direction group；Wherein, formula (3) is respectively as follows: to formula (7)

If y_M' > 0, then have

$P=\arccos \left(\frac{H*M}{\left|H\right|*\left|M\right|}\right)*\frac{180}{\mathrm{\&pi;}}=\arccos \left(\frac{x_M^{\&prime;}}{\sqrt{{\left(x_M^{\&prime;}\right)}^2+{\left(y_M^{\&prime;}\right)}^2}}\right)*\frac{180}{\mathrm{\&pi;}}---\left(3\right)$

If y_M' < 0, then have

$P=\arccos \left(\frac{H*M}{\left|H\right|*\left|M\right|}\right)*\frac{180}{\mathrm{\&pi;}}+180=\arccos \left(\frac{x_M^{\&prime;}}{\sqrt{{\left(x_M^{\&prime;}\right)}^2+{\left(y_M^{\&prime;}\right)}^2}}\right)*\frac{180}{\mathrm{\&pi;}}+180---\left(4\right)$

$\mathrm{\&gamma;}(h,\mathrm{\&theta;})=\mathrm{\&gamma;}\left(\right|x_i-x_j|,\mathrm{\&theta;})=\frac{1}{2}E{\&lsqb;Z\left(x_i\right)-Z\left(x_j\right)|\mathrm{\&theta;}\&rsqb;}^2---\left(5\right)$

${\mathrm{\&gamma;}}^{*}\left(h,\mathrm{\&theta;}\right)={\mathrm{\&gamma;}}^{*}\left(|x_i-x_j|,\mathrm{\&theta;}\right)=\frac{1}{2m}{\mathrm{\&Sigma;}}_{k=1}^m{\&lsqb;Z\left(x_i\right)-Z\left(x_j\right)|\mathrm{\&theta;}\&rsqb;}^2---\left(6\right)$

$h=\sqrt{{({x_i}^{\&prime;}-{x_j}^{\&prime;})}^2+{({y_i}^{\&prime;}-{y_j}^{\&prime;})}^2}---\left(7\right)$

Wherein, H (100,0) is reference vector, M (x_M',y_M') it is two described sample point x_i,x_jVector M (the x formed_i'-x_j',y_i'-y_j'), P is vector M and reference vector H institute angle degree, and P is sample point x described in two_i,x_jBetween position angle θ, π be pi, γ (h, θ) is desirable variation function, and γ * (h, θ) is discrete variation function, x_i,x_jThe sample point pair being h for distance, Z (x_i)、Z(x_j) for sample point x_i,x_jProperty value, h is that sample point is to x_i,x_jBetween distance, (x_i',y_i') for known sample point x_iLatitude and longitude coordinates, (x_j',y_j') for known sample point x_jLatitude and longitude coordinates, θ is described sample point x_i,x_jBetween position angle, m represents the number of the sample point pair that distance is h.

12. the adaptive space interplotation system based on space characteristics analysis as claimed in claim 11, it is characterised in that also including fitting module, described fitting module specifically includes:

Computing unit, is used for calculating the variation function value under different distance according to formula (7) with to formula (9)；

$\mathrm{\&gamma;}\left(h\right)=\mathrm{\&gamma;}\left(\right|x_i-x_j\left|\right)=\frac{1}{2}E{\&lsqb;Z\left(x_i\right)-Z\left(x_j\right)\&rsqb;}^2---\left(8\right)$

${\mathrm{\&gamma;}}^{*}\left(h\right)={\mathrm{\&gamma;}}^{*}\left(\right|x_i-x_j\left|\right)=\frac{1}{2m}{\mathrm{\&Sigma;}}_{k=1}^m{\&lsqb;Z\left(x_i\right)-Z\left(x_j\right)\&rsqb;}^2---\left(9\right)$

Wherein γ (h) is desirable variation function, γ^*H () is discrete variation function, x_i,x_jThe sample point pair being h for distance, Z (x_i)、Z(x_j) for sample point x_i,x_jProperty value, h is sample point x_i,x_jBetween distance, (x_i',y_i') for known sample point x_iLatitude and longitude coordinates, (x_j',y_j') for known sample point x_jLatitude and longitude coordinates, m represents the number of the sample point pair that distance is h；

Drawing unit, is used for distance h for abscissa, variation function value γ^*H () is vertical coordinate, draw variation function scatterplot；

Fitting unit, for choosing spherical model, exponential model or Gauss model to variation function scatterplot matching to obtain the highest model of degree of fitting and function expression γ (h) thereof.

13. the adaptive space interplotation system based on space characteristics analysis as claimed in claim 12, it is characterised in that described determination module specifically includes:

First identifying unit, for calculating range value according to the variation function in described all directions, and draws rose diagram according to described range value；

According to described rose diagram, second identifying unit, for judging that the spatial distribution of described pollutant is as isotropism or anisotropy.

14. the adaptive space interplotation system based on space characteristics analysis as claimed in claim 13, it is characterised in that described search module specifically includes:

Setup unit, estimates a number for setting minimum described ginseng；

First search unit, estimates a little for searching for described ginseng when the spatial distribution of described pollutant is isotropism；

Second search unit, estimates a little for searching for described ginseng when the spatial distribution of described pollutant is anisotropy.

15. the adaptive space interplotation system based on space characteristics analysis as claimed in claim 14, it is characterised in that also include:

Divide acquisition module, for Research on partition region, generate multiple mesh point and obtain the center point coordinate of each described grid；

Interpolating module, for using the central point of each described grid as described interpolation point, adopt and described determine that the central point of each described grid is interpolated by module and calculate interpolation result to obtain each described mesh point.

16. the adaptive space interplotation system based on space characteristics analysis as claimed in claim 15, it is characterised in that also include:

Authentication module, is used for adopting cross validation method that described interpolation method is verified.

17. the adaptive space interplotation system based on space characteristics analysis as described in claim 15 or 16, it is characterised in that also include:

Rendering module, the interpolation result of each the described mesh point for obtaining according to described interpolating module, different numerical value sections is given different color, space based on GIS platform drawing area air quality renders figure.