CN106683185A - High accuracy surface modeling method based on big data - Google Patents

Abstract

The invention relates to a high accuracy surface modeling method based on big data. The method comprises the following steps: creating geographic coordinate information and to-be-tested variable sampling values of all sampling points; turning a to-be-tested regional space into a grid point form through discretization, and establishing a sampling equation; conducting high order difference discretization on a partial differential equation set of a surface, obtaining a corresponding algebraic system, combining the algebraic system and the sampling equation to form an equality constrain least squares problem, converting the problem into a solve cross cutting object function minimum value problem, and then converting the minimum value problem into a solve symmetry indeterminate equation set; randomly selecting iteration initial values; blocking a coefficient matrix, and storing the block matrixes after decomposition; solving an HASM equation set through a block line projection iterative method, and determining whether a solving result is convergent; determining whether a solution meets demands of a Gauss and C.odazzi equation set; and outputting a high accuracy simulation surface model. Large-scale problems can be solved, the needed storage space is small, and defects of HASM for solving large-scale problems are prevented.

Description

A kind of High Accuracy Surface Modeling method based on big data

Technical field

The present invention relates to be related to a kind of method that space curved surface for extensive problem is modeled, can be used for big data background The fields such as generation, meteorology, Biomass, the spatial distribution simulation of soil of lower digital elevation model, it is also possible to be considered as curved surface grid A kind of method approached, can be used for the surface modeling of the aspects such as large-scale physics, chemistry, medical science.

Background technology

In order to perplex the error problem and Issues On Multi-scales of curved surface modeling since solving half a century, we are several with differential What principle and optimal control opinion are theoretical basiss, establish one with approximate data of overall importance (including remotely-sensed data and global mould Type coarse resolution analog data) be for driving field, with local high accuracy data (including monitoring network data and investigation sampled data) High Accuracy Surface Modeling (High Accuracy Surface Modelling, the be abbreviated as HASM) method of optimal control condition, And widely applied in more than 20 years on Research foundation, refinement defines epigeosphere modeling philosophy.

HASM can finally be converted into an equality constraint least square problem (Equality by ground sampling constraints constraint least squares problem:LSE), in order to sample ensureing that the sample point analogue value is equal to Under conditions of value, keep overall simulation error minimum.Make full use of sample information, be also ensure iteration level off to best simulation effect The effective means of fruit.Although HASM methods have very big improvement than traditional interpolation method on simulation precision, table is studied Bright, current HASM can only process small-scale problem, and above-mentioned solution procedure is that LSE problems are converted into into normal equation system to solve.One Aspect, the equation group condition number of coefficient matrix after conversion is very big, causes the sensitivity and pathosis of problem, and this causes to solve The alternative manner convergence rate of equation group is very slow, or does not restrain, and thereby results in and can cannot get the approximate of former Solve problems Solution；On the other hand, currently employed preconditioning conjugate gradient solves HASM equation group to be needed equation group coefficient matrix whole Disposable input internal memory, with the increase of Solve problems scale, its memory space can rise appreciably, for extensive problem work as Front solution strategies are not used.Current HASM often processes plan for the simulation of national 1km resolution datas using burst Slightly, the error at partition boundaries is which increased, so as to reduce the overall precision of analog result.

The content of the invention

Defect of the present invention for existing HASM models in extensive problem is solved, there is provided a kind of improved HASM, required memory space is little, can solve extensive problem in the case where precision is ensured.

The technical scheme that the present invention solves above-mentioned technical problem is as follows：A kind of High Accuracy Surface Modeling side based on big data Method, comprises the following steps：

Step 1：The geographic coordinate information and variable sampled value to be measured of each sampled point are created, is wrapped in the geographic coordinate information Include the longitude information of sampled point and the latitude information of sampled point；

Step 2：Mesh point form is turned to by regional space to be measured is discrete, mesh point centrifugal pump is obtained, according to geographical coordinate Information and variable sampled value to be measured set up sampling equation, and the sampling equation is used to whether judge sampled point in mesh point, its In, if sampled point is on the mesh point, the value of the mesh point is variable sampled value to be measured；If sampled point is in grid Approximation sample value that is interior, then will being obtained using Taylor expansion on the nearest mesh point of the sampled point on the mesh point；

Step 3：First kind fundamental quantity E, F, G and second of each mesh point of region to be measured is calculated according to mesh point centrifugal pump Class fundamental quantity L, M, N, wherein first fundamental quantity be used for represent simulation Curves on surfaces length, simulation curved surface area and The curvature of simulation Curves on surfaces, the second fundamental quantity is used to represent the local buckling intensity of variation of simulation curved surface, will be with the It is discrete that the partial differential equations of the curved surface that one fundamental quantity and second fundamental quantity are represented carry out higher difference, obtains discrete equation group pair The algebra system answered, equality constraint least square problem is combined into by the algebra system with the sampling equation；

Step 4：Equality constraint least square is converted into into solution for topic and blocks object function minimum problem；

Step 5：Object function minimum problem is blocked in solution and is converted into the symmetrical Indeterminate Equation Group of solution, it is described it is symmetrical not Equation group is determined for High Accuracy Surface Modeling HASM equation group；

Specifically,

Step 6：Randomly select the iterative initial value of HASM equation group；

Step 7：Coefficient matrix in HASM equation group is entered into every trade piecemeal, and the block matrix after decomposing to coefficient matrix enters Row storage；

Step 8：Kuai Hang projection iterative methods are adopted to iterative initial value, HASM equation group is solved, and whether judges solving result Convergence；

Step 9：When solving result is not restrained, Kuai Hang projection iterative methods are adopted again to solving result, solve HASM side Journey group, and judge whether solving result restrains, if convergence, execution step 10 otherwise, re-executes step 9；

Step 10：When HASM solution of equations is restrained, determine whether whether HASM solution of equations meets Gauss section Up to neat equation group, if being unsatisfactory for, execution step 6；If meeting, exported with regard to variable to be measured according to HASM solution of equations High-precision analog surface model.

Present invention seek to address that the big data problem in space curved surface modeling field, has the advantage that：

1st, improve and developed existing model, when the model after improvement has relatively low storage demand and less calculating Between；

2nd, truncated system is introduced, gives the optimization method for solving constrained least-squares problem, it is as a result more stable；

3rd, give large-scale equation group coefficient matrix specific partitioned mode, it is ensured that the Line independent between matrix rows Property, it is to avoid information redundancy, so as to give the accelerating algorithm Ji Kuaihang projection iterative method of row projection algorithm, improves method Convergence rate；

4th, the method after improving can solve extensive problem, and the calculating time is linear with calculating grid number, required to deposit Storage space is far smaller than nonzero element number in HASM equation group coefficient matrixes；In the case where precision is ensured, HASM is solved Defect in extensive problem is solved.

On the basis of above-mentioned technical proposal, the present invention can also do following improvement.

Further, the partial differential equations of the curved surface are：

Wherein, E=1+f_x ², F=f_x·f_y,

X for spatial discretization mesh point abscissa, y for spatial discretization mesh point vertical coordinate, f_xFor function f First-order partial derivative to x, f_xxFor second-order partial differential coefficients of the function f to x, f_yFor first-order partial derivatives of the function f to y, f_yyFor function f pair The second-order partial differential coefficient E of y_x、F_x、G_xThe respectively first-order partial derivative of E, F, G to x, E_y、F_y、G_yRespectively E, F, G is inclined to the single order of y Derivative, For Equations of The Second Kind Cristoffel symbols.

It is using the beneficial effect of above-mentioned further scheme so that HASM is set up in complete differential geometry theoretical basiss On, it is ensured that the stability of HASM, improve the simulation precision of HASM.

Further, the step 3 is specifically included：

Step 3.1：Calculated according to the finite difference approximation of first kind fundamental quantity E, F, G and Equations of The Second Kind fundamental quantity L, M, N and treated Survey the first kind fundamental quantity of each mesh point of region and the value of Equations of The Second Kind fundamental quantity；

Step 3.2：Finite difference approximation, the finite difference approximation of Equations of The Second Kind fundamental quantity with reference to first kind fundamental quantity and The finite difference of two class Cristoffel symbols obtains discrete equation group；

Step 3.3：The discrete equation group is converted into into the algebra system of matrix form；

Step 3.4：The algebra system is combined into into equality constraint least square problem with sampling equation.

It is that HASM is finally translated into an equation by ground sampling constraints using the beneficial effect of above-mentioned further scheme Constrained least-squares problem, in order under conditions of ensureing that the sample point analogue value is equal to sampled value, keep non-detachable mold Intend error minimum.Sample information is made full use of, is also to ensure that iteration levels off to the effective means of best simulation effect.

Further, it is by the concrete grammar that the coefficient matrix in HASM equation group enters every trade piecemeal in the step 7：If being Each sub-block A of matrix number A_iBe up to u rows, and A_i∈R^u×n, it is independent that u ＜ ＜ n are full rank, the i.e. row per block, it is assumed that Generating sub-block A_i, the sub-block is comprising the j rows (j ＜ u) in matrix A, and known A_i ^TQR be decomposed into A_i ^T=Q_jR_j, Block A_i ^TRow dependency estimationWhereinr_hhIt is positioned at R_j(h, H) element of position, a₁For a line in matrix A, orderAccording to A_i ^TQR decompose, using QR decompose update method To obtainQR decompose, then calculateRow correlation estimation, if the estimation is still less than κ, a₁Add block A_i, its In, κ is a positive number set in advance.

Further, the QR decomposes using Francis QR decomposition strategies.

It is, for the QR being related to decomposes, to be decomposed using Francis QR using the beneficial effect of above-mentioned further scheme Strategy, can be calculated cost by O (u³) it is reduced to O (u²), u ＜ ＜ n, so as to improve computational efficiency.The partitioning of matrix can be accelerated The convergence of row projection iterative method, reduces HASM solving equations and calculates the time.

Further, the mode for being stored to the block matrix after coefficient matrix decomposition in the step 7 is row compression storage Mode CSR.

It is that CSR forms are most common and flexible compressed format using the beneficial effect of above-mentioned further scheme, it is by square The nonzero element of battle array carries out being stored by row compression, and the original position of nonzero element is recorded with special array, and compression efficiency is high, Compression process readily appreciates, thereby dramatically reduces memory requirements.

Further, the step 8 and step 9 are specifically included：Assume that current iteration point is to judge x^kWhether restrain, if not Convergence, then find out from x^kOrthogonal chooses x apart from farthest block^kProjection on described piece is used as next iteration point.

It is to accelerate convergence rate using the beneficial effect of above-mentioned further scheme.

Description of the drawings

Fig. 1 is a kind of High Accuracy Surface Modeling method and step flow chart based on big data provided in an embodiment of the present invention；

Fig. 2 is standard mathematical curved surface and sampling point distributions figure；

Fig. 3 is the Error Graph that HASM, FRK, IDW and FRS simulate curved surface and true curved surface.

Specific embodiment

The principle and feature of the present invention are described below in conjunction with accompanying drawing, example is served only for explaining the present invention, and It is non-for limiting the scope of the present invention.

As shown in figure 1, a kind of High Accuracy Surface Modeling method based on big data provided in an embodiment of the present invention, including with Lower step：

Step 1：The geographic coordinate information and variable sampled value to be measured of each sampled point are created, is wrapped in the geographic coordinate information Include the longitude information of sampled point and the latitude information of sampled point；

Step 2：Mesh point form is turned to by regional space to be measured is discrete, mesh point centrifugal pump is obtained, according to geographical coordinate Information and variable sampled value to be measured set up sampling equation, and the sampling equation is used to whether judge sampled point in mesh point, its In, if sampled point is on the mesh point, the value of the mesh point is variable sampled value to be measured；If sampled point is in grid Approximation sample value that is interior, then will being obtained using Taylor expansion on the nearest mesh point of the sampled point on the mesh point；

Step 3：First kind fundamental quantity E, F, G and second of each mesh point of region to be measured is calculated according to mesh point centrifugal pump Class fundamental quantity L, M, N, wherein first fundamental quantity be used for represent simulation Curves on surfaces length, simulation curved surface area and The curvature of simulation Curves on surfaces, the second fundamental quantity is used to represent the local buckling intensity of variation of simulation curved surface, will be with the It is discrete that the partial differential equations of the curved surface that one fundamental quantity and second fundamental quantity are represented carry out higher difference, obtains discrete equation group pair The algebra system answered, equality constraint least square problem is combined into by the algebra system with the sampling equation；

The step 3 is specifically included：

Step 3.1：Calculated according to the finite difference approximation of first kind fundamental quantity E, F, G and Equations of The Second Kind fundamental quantity L, M, N and treated Survey the first kind fundamental quantity of each mesh point of region and the value of Equations of The Second Kind fundamental quantity；

Step 3.2：Finite difference approximation, the finite difference approximation of Equations of The Second Kind fundamental quantity with reference to first kind fundamental quantity and The finite difference of two class Cristoffel symbols obtains discrete equation group；

Step 3.3：The discrete equation group is converted into into the algebra system of matrix form；

Step 3.4：The algebra system is combined into into equality constraint least square problem with sampling equation.

Specifically, according to surface theory fundamental theorem：If the first kind of curved surface and Equations of The Second Kind fundamental quantity E, F, G, L, M and N are full Sufficient symmetry, E, F, G positive definite, E, F, G, L, M and N meet the partial differential equations of curved surface, i.e. Gauss equation group, then total differential Equation group is in f (x, y)=f (x₀,y₀) (x=x₀, y=y₀) under initial condition, there is unique solution z=f (x, y).

The expression formula of Gauss equation group is,

X for spatial discretization mesh point abscissa, y for spatial discretization mesh point vertical coordinate, f_xFor function f First-order partial derivative to x, f_xxFor second-order partial differential coefficients of the function f to x, f_yFor first-order partial derivatives of the function f to y, f_yyFor function f pair The second-order partial differential coefficient E of y_x、F_x、G_xThe respectively first-order partial derivative of E, F, G to x, E_y、F_y、G_yRespectively E, F, G is inclined to the single order of y Derivative, For Equations of The Second Kind Cristoffel symbols.

If { (x_i,y_j) be the orthogonal subdivision of computational fields Ω, [0, Lx] × [0, Ly] be dimensionless standardization computational fields, max Numerical value is rounded after { Lx, Ly }=1, the normalization that Lx is survey region horizontal length.Ly is survey region vertical-direction length Normalization after round numerical value.For material calculation, { (x_i,y_j) | 0≤i≤I+1,0≤j≤J+1 } it is mark The grid of standardization computational fields, then the finite difference approximation of first kind fundamental quantity is,

Wherein, E_i,j, F_i,j, G_i,j, be respectively the value of E, F, G on mesh point (i, j), f_i,jTo simulate curved surface in mesh point Value on (i, j)；

The finite difference approximation of Equations of The Second Kind fundamental quantity is,

Wherein, L_i,j, M_i,j, N_i,jThe respectively value of L, M, N on mesh point (i, j)；

The finite difference of Equations of The Second Kind Cristoffel symbols is expressed as,

Wherein,For Equations of The Second Kind Christoffer Value of that symbol on mesh point (i, j)；

The finite difference form of Gauss equation group is,

(1.1.2) matrix form can be written as：

Wherein,

Wherein,

For value of above-mentioned (1.1.2) the formula right-hand vector using finite difference scheme after discrete.

Control with reference to the operative constraint of sample information, above-mentioned constrained least-squares problem (1.1.3) can be expressed as HASM institutes The least square problem of the equality constraint of solution,

Sx=k is the sampling equation that sampled point meets, and wherein S and k is respectively sampling matrix and vector of samples；If It is z=f (x, y) in pth sampled point (x_i,y_j) value, then s_p,(i-1)×J+j=1,

HASM is finally translated into an equality constraint least square problem (Equality by ground sampling constraints constraint least squares problem:LSE)

Step 4：Equality constraint least square is converted into into solution for topic and blocks object function minimum problem；

Specifically, by studying discovery, a parameter lambda is introduced, LSE problems can first pass through Lagrange multiplier methods by its turn To turn to ask as follows and block object function minimum problem：

Wherein

Step 5：Object function minimum problem is blocked in solution and is converted into the symmetrical Indeterminate Equation Group of solution, it is described it is symmetrical not Equation group is determined for High Accuracy Surface Modeling HASM equation group；

Specifically, it is rightDerivation, another its derivative is 0, can be obtainedWith reference toAnd Sx=k can obtain symmetrical indefinite system as follows System：

Wherein, I is all 1 for diagonal entry, and other elements are 0 unit matrix.For matrixTransposed matrix.S^T For the transposed matrix of S.λ is the weight matrix that sampled point weight is constituted, and value is between 1~2 in HASM.

Step 6：Randomly select the iterative initial value of HASM equation group；

Step 7：Coefficient matrix in HASM equation group is entered into every trade piecemeal, and the block matrix after decomposing to coefficient matrix enters Row storage；

Specifically, for the partitioned mode of coefficient matrix, this research gives a block algorithm of A, using the algorithm per block Be up to u rows, and every piece of A_i∈R^u×n, it is independent that u ＜ ＜ n are full rank, the i.e. row per block.Hypothesis is generating a block A_i, the block is comprising the j rows (j ＜ u) in matrix A, and known A_i ^TQR be decomposed into A_i ^T=Q_jR_j, make r_hhIt is positioned at R_j(h, H) element of position.Block A_i ^TRow dependency estimationWherein For a line a in matrix A₁Whether block A is added_iCriterion be a₁It is not allocated to any other piece, and a₁After addition A_iRow correlation estimation β_iLess than κ, κ is a positive number set in advance.In order to determine a new row a₁Whether block is added A_i, orderAnd according to existing A_i ^TQR decompose obtainingQR decompose, to this using QR decompose renewal side Method is realizing.Then calculateRow correlation estimation, if the estimation is still less than κ, a₁Add block A_i。

The block decomposition algorithm is described as follows, and the algorithm is realized based on C++.

Γ=1,2 are made ..., m, Γ_s={ l:A in A₁Row have distributed to certain block }

Initialization：Make i=1, Γ_s=φ,

whileΓ_s≠Γ,do

Make Γ_c=Γ Γ_s, j=1, Γ_i=φ

Choose l ∈ Γ_c, make A_i=[a₁],R₁=[| | a₁||₂,0]^T,

CalculateQ₁=I_n-2vv^T

Make d_max=| | a₁||₂, d_min=| | a₁||₂

Γ_c←Γ_c{ l }, Γ_s←Γ_s∪ { l }, Γ_i←Γ_i∪{l}

While (j ＜ u and Γ_c≠φ)do

Choose l ∈ Γ_c

OrderV=[| | a | |₂,0]^T∈R^n-j

Calculate

t_max←max{d_max,||a||₂},t_min←min{d_min,||a||₂},

if t_min≠ 0, andThen

Update A_i←[A_i,a₁]

Γ_c←Γ_c{ l }, Γ_s←Γ_s∪ { l }, Γ_i←Γ_i∪{l}

J=j+1

d_max←t_max,d_min←t_min

else

Γ_c←Γ_c\{l}

end

end

i←i+1

end

In order to improve computational efficiency, for the QR being related to decomposes, we employ Francis QR decomposition strategies, and will It calculates and spends by O (u³) it is reduced to O (u²), u ＜ ＜ n.

In above-mentioned algorithm, it is related to sparse matrix block A_iBetween multiplication and sparse matrix and the computing taken advantage of of vector, in vector Product and vector addition computing.Sparse matrix is typically stored as compressed format.With application scenarios and the difference of calculating platform, matrix Generally it is compressed to different storage formats.The selection of compressed format will consider the sparse feature and calculating platform of matrix. To matrix storage using row compression storage mode (CSR), CSR forms are most common and flexible compressed format for this research, and it will Nonzeros carries out being stored by row compression, and the original position of nonzero element, compression efficiency are recorded with special array Height, compression process is readily appreciated.CSR forms and in being transplanted in different platform.Many new compressed formats are also all based on CSR Form modifying.During this method storage n rank matrix As, it is assumed that l nonzero element is had in A, then need with a l tie up to Amount x is deposited successively the nonzero element in A by the order of Row Column, with a l dimensional vector x^JA is deposited successively by same order In these nonzero element row numbers, while also need to introduce n+1 tieing up integer vector x^(R),Indicate i-th in A First nonzero element is stored in the position in x in row,The coefficient matrix of Indeterminate Equation Group (4) in this research For symmetrical matrix, therefore triangular portions nonzero element need to be only stored in practice.It is corresponding to employ block row pressure for Huge types Contracting storage mode.

The matrix realized based on above-mentioned form can be written as with vector multiplication y=Av：

y(x^(R)(i))=x (i) × v (x^J(i))

end

Storage mode is compressed based on CSR rows, can analyze the calculating of above-mentioned algorithm spends and is O (n) to the maximum, wherein n is meter Calculate grid number.

Step 8：Kuai Hang projection iterative methods are adopted to iterative initial value, HASM equation group is solved, and whether judges solving result Convergence；

Step 9：When solving result is not restrained, Kuai Hang projection iterative methods are adopted again to solving result, solve HASM side Journey group, and judge whether solving result restrains, if convergence, execution step 10 otherwise, re-executes step 9；

Specifically, for the method for solving of above-mentioned symmetrical uncertain system is generally direct method and iterative method.Direct method is related to Decomposition to matrix etc. is not appropriate for extensive problem.Conventional iterative method is Krylov subspace iterative method, is located more in advance Reason conjugate gradient method PCG, the QR decomposition method LSQR based on least square, and minimum residual error method GMRES of broad sense etc..Due to Such iterative method needs that equation group coefficient matrix is all disposably input into internal memory, with the increase of calculation scale, such method Do not use.Row projection iterative method can avoid problems, row projection iterative method from being widely used, and algorithm is simple, classical row Projection iterative method such as Cimmio algorithms, Kaczmar algorithms etc..

Define H_iFor set：RⁿThe space that real number vector is constituted, i=1 ..., p are tieed up for n. Therefore it is anyThat is x^*For p sub-spaces H_iCommon factor, be linear system A^TThe solution of x=b.If required (1.1.5) In symmetrical uncertain system be Ax=b, is deducted marks according to the grid of zoning and matrix A is decomposed into block row form：

And define P_k(A^iT) it is A^iTProjection operator in codomain.General thought is first by x_kProject to P hyperplane

Algorithm:Block Cimmino row sciagraphies:

Choose x⁽⁰⁾, k=0 is made,

Steps be repeated alternatively until convergence,

begin

Do in parallel i=1 ... p

δ^i(k)=Aⁱ⁺bⁱ-P_A(A^iT)x^(k)=Aⁱ⁺(bⁱ-Aⁱx^(k))

end parallel

Set k=k+1

end

The method can be avoided without fully entering matrix element every time, and one piece of row of input matrix is only needed every time, protected The solution probability of extensive problem is demonstrate,proved.But have the disadvantage that convergence rate is slow, therefore less application in practice.

Kaczmar algorithms are based on for this this research, a new block-iterative algorithm with row projection is given.By choose from work as Front iteration point carries out orthogonal projection apart from farthest block, and will project as next iteration point, so as to accelerate convergence speed Degree.Assume that current iteration point is x^k, x is calculated first^kIn all H_iOn orthogonal projection P_i(x^k), for i=1 ..., then p compares It is more all of | | P_i(x^k)-x^k||₂And obtain wherein the maximum, that is, find out from x^kOrthogonal is apart from farthest block H_i, it is assumed that | | P_j (x^k)-x^k||₂Maximum, i.e.,Choose x^kIn H_jOn projection P_j(x^k) as next Individual iteration point, i.e. x^k+1=P_j(x^k), the process for continuing the above obtains a sequence { x^k, continuous iteration is until convergence.The algorithm is retouched State as follows：

Initialization：Piecemeal is carried out to matrix A, initial point x is given⁰∈Rⁿ, 0≤ε ＜ 1 make k=0

Step 10：When HASM solution of equations is restrained, determine whether whether HASM solution of equations meets Gauss section Up to neat equation group, if being unsatisfactory for, execution step 6；If meeting, exported with regard to variable to be measured according to HASM solution of equations High-precision analog surface model..

The checking of the embodiment of the present invention is described below,

1st, numerical experimentation

First by taking standard mathematical curved surface as an example, the effectiveness of algorithm is studied.Compare the HASM methods after improving and be suitable for The Kriging methods (FRK) of extensive Solve problems, the Spline methods (FRS) for being suitable for extensive Solve problems and it is anti-away from From the simulated performance difference of method of weighting IDW.

Domain of definition is [0,1] × [0,1].

As shown in figure 1, randomly selecting 1/100 sampled point on standard mathematical curved surface, all experiments are transported on 64 machines OK, it is Gauss-Codazii equation group that criterion is shut down in HASM.Zoning is deducted marks as 1000 × 1000.u=100, κ= 5. the measurement of pair calculation error, using mean absolute error RMSE, its computing formula is as follows：

Wherein f_i,For the i-th th actual value and the analogue value, N is sampled point number, is counted Calculate result as shown in table 1：

The distinct methods computational efficiency com-parison and analysis of table 1

Method	HASM	G-IDW	FRK	FRS
					RMSE	1.3076×10-4	1.66×10-2	2.20×10-3	5.78×10-4
Calculating time (s)	54.06s	20.11s	40.00s	24.79s

As a result show, HASM computational accuracy highests, next to that FRS methods.HASM computational accuracies are respectively IDW, FRK and 127 times of FRS, 17 times and 4 times.But the calculating speed of HASM is less than additive method.This is because needing in HASM calculating process It is discrete to partial differential equations finite difference, need group of equations right-hand vector, and solving equation group etc., calculating process is more multiple It is miscellaneous.

(a) HASM according to Fig. 2, (b) FRK, (c) G-IDW, and (d) FRS simulate the Error Graph of curved surface and true curved surface, Therein, it can be seen that IDW simulation errors are maximum, maximum error is at border and peak value.Maximum deviation is 0.2139.FRK Also larger error is generated, maximum error occurs in boundary for -0.0123. errors most for the maximum deviation of 0.02737.FRS Little is HASM methods, it can be seen that error distribution is relatively uniform.

Although the calculating time of HASM is more than additive method, it is linearly proportional with calculating grid that it calculates the time. Table 2 gives the calculating time of HASM under different calculation scales.The calculating time can be expressed as with the relational expression for calculating grid：

T=1.1099 × 10^-5n+26.7447,R²=0.91

The calculating time of HASM under the different calculation scales of table 2.

Calculate grid number	1million	4million	9million	16million	25million
						Calculating time (s)	54.06	88.21	107.54	155.19	339.15

2nd, real case

With China to study area, the distribution of many mean annual precipitations of past 50 years 1961-2010 is calculated using said method, Spatial resolution is 1km, and it is 19,606,916 that correspondence calculates grid number.

All sides are can be seen that from many mean annual precipitation distribution results of Chinese scope of HASM, FRK, IDW and FRS simulation Method can preferably reflect the distribution trend of precipitation.The curved surface of HASM and FRK methods simulation is with respect to additive method more more It is smooth.IDW and FRS methods occur in that many buphthalmos phenomenons, and FRS methods generate vibration by a relatively large margin in boundary.

Random selection 15% as check post verify distinct methods.This process is repeated 10 times, and calculates this 10 times Average RMSE.As a result it is as shown in table 3.It can be seen that HASM computational accuracies are highests, its simulation precision distinguishes my FRK, IDW and 1.4,1.5 and 1.6 times of FRS.The simulation error of FRS is maximum.IDW calculates the time at least, next to that HASM.

The comparison of distinct methods in the real case of table 3.

Method	HASM	FRK	IDW	FRS
					RMSE(mm)	90.87	123.51	133.44	147.10
Calculating time (s)	166	237	63	172

Numerical experimentation and real case, show the HASM methods after improving while precision is ensured so that during its calculating Between with calculate grid number it is linear, in storage, due to solve constrained least-squares problem when, employ block row projection Algorithm, only needs one piece of row of storage matrix, and it is stored using row compression storage mode, it is ensured that the rule of Solve problems Mould.

The foregoing is only presently preferred embodiments of the present invention, not to limit the present invention, all spirit in the present invention and Within principle, any modification, equivalent substitution and improvements made etc. should be included within the scope of the present invention.

Claims (7)

1. a kind of High Accuracy Surface Modeling method based on big data, comprises the following steps：

Step 1：The geographic coordinate information and variable sampled value to be measured of each sampled point are created, the geographic coordinate information includes adopting The longitude information of sampling point and the latitude information of sampled point；

Step 2：Mesh point form is turned to by regional space to be measured is discrete, mesh point centrifugal pump is obtained, according to geographic coordinate information Sampling equation is set up with variable sampled value to be measured, the sampling equation is used to whether judge sampled point in mesh point, wherein, such as On the mesh point, then the value of the mesh point is variable sampled value to be measured to fruit sampled point；If sampled point is in grid, The approximation sample value that will be obtained using Taylor expansion on the nearest mesh point of the sampled point on the mesh point；

Step 3：First kind fundamental quantity E, F, G and Equations of The Second Kind base of each mesh point of region to be measured are calculated according to mesh point centrifugal pump This amount L, M, N, wherein first fundamental quantity is used to represent the length of simulation Curves on surfaces, the area of simulation curved surface and simulation The curvature of Curves on surfaces, the second fundamental quantity is used for the local buckling intensity of variation of expression simulation curved surface, will use the first base It is discrete that the partial differential equations of the curved surface that this amount and second fundamental quantity are represented carry out higher difference, obtains discrete equation group corresponding Algebra system, equality constraint least square problem is combined into by the algebra system with the sampling equation；

Step 4：Equality constraint least square is converted into into solution for topic and blocks object function minimum problem；

Step 5：Object function minimum problem is blocked in solution and is converted into the symmetrical Indeterminate Equation Group of solution, the symmetrical indefinite side Journey group is High Accuracy Surface Modeling HASM equation group；

Step 6：Randomly select the iterative initial value of HASM equation group；

Step 7：Coefficient matrix in HASM equation group is entered into every trade piecemeal, and the block matrix after decomposing to coefficient matrix is deposited Storage；

Step 8：Kuai Hang projection iterative methods are adopted to iterative initial value, HASM equation group is solved, and judges whether solving result restrains；

Step 9：When solving result is not restrained, Kuai Hang projection iterative methods are adopted again to solving result, solve HASM equation group, And judge whether solving result restrains, if convergence, execution step 10 otherwise, re-executes step 9；

Step 10：When HASM solution of equations is restrained, determine whether whether HASM solution of equations meets Gauss Kodak neat Equation group, if being unsatisfactory for, execution step 6；If meeting, exported with regard to the high-precision of variable to be measured according to HASM solution of equations Degree simulation surface model.

2. method according to claim 1, it is characterised in that the partial differential equations of the curved surface are：

$\left\{\begin{array}{c}{f}_{xx}={\mathrm{\Γ}}_{11}^1{f}_x+{\mathrm{\Γ}}_{11}^2{f}_y+\frac{L}{\sqrt{E+G-1}}\\ f_{yy}=T_{22}^1{f}_x+{\mathrm{\Γ}}_{22}^2{f}_y+\frac{N}{\sqrt{E+G-1}}\\ f_{xy}={\mathrm{\Γ}}_{12}^1{f}_x+{\mathrm{\Γ}}_{12}^2{f}_y+\frac{M}{\sqrt{E+G-1}}\end{array}\right.$

Wherein,F=f_x·f_y,

$L=\frac{f_{xx}}{\sqrt{1+f_x^2+f_y^2}},M=\frac{f_{xy}}{\sqrt{1+f_x^2+f_y^2}},N=\frac{f_{yy}}{\sqrt{1+f_x^2+f_y^2}},$

${\mathrm{\Γ}}_{11}^1=\frac{1}{2}({\mathrm{GE}}_x-2{\mathrm{FF}}_x+{\mathrm{FE}}_y){(EG-F^2)}^{-1},{\mathrm{\Γ}}_{11}^2=\frac{1}{2}(2{\mathrm{EF}}_x-{\mathrm{EE}}_y-{\mathrm{FE}}_x){(EG-F^2)}^{-1},$

${\mathrm{\Γ}}_{22}^1=\frac{1}{2}(2{\mathrm{GF}}_y-{\mathrm{GG}}_x-{\mathrm{FG}}_y){(EG-F^2)}^{-1},{\mathrm{\Γ}}_{22}^2=\frac{1}{2}({\mathrm{EG}}_y-2{\mathrm{FF}}_y+{\mathrm{FG}}_x){(EG-F^2)}^{-1},$

${\mathrm{\Γ}}_{12}^1=\frac{1}{2}({\mathrm{GE}}_y-{\mathrm{FG}}_x){(EG-F^2)}^{-1},{\mathrm{\Γ}}_{12}^2=\frac{1}{2}({\mathrm{EG}}_x-{\mathrm{FE}}_y){(EG-F^2)}^{-1},$

X for spatial discretization mesh point abscissa, y for spatial discretization mesh point vertical coordinate, f_xIt is function f to x's First-order partial derivative, f_xxFor second-order partial differential coefficients of the function f to x, f_yFor first-order partial derivatives of the function f to y, f_yyBe function f to y two Rank partial derivative E_x、F_x、G_xThe respectively first-order partial derivative of E, F, G to x, E_y、F_y、G_yThe respectively first-order partial derivative of E, F, G to y, For Equations of The Second Kind Cristoffel symbols.

3. method according to claim 2, it is characterised in that the step 3 is specifically included：

Step 3.1：Area to be measured is calculated according to the finite difference approximation of first kind fundamental quantity E, F, G and Equations of The Second Kind fundamental quantity L, M, N The first kind fundamental quantity of each mesh point of domain and the value of Equations of The Second Kind fundamental quantity；

Step 3.2：Finite difference approximation, the finite difference approximation of Equations of The Second Kind fundamental quantity and Equations of The Second Kind with reference to first kind fundamental quantity The finite difference of Cristoffel symbols obtains discrete equation group；

Step 3.3：The discrete equation group is converted into into the algebra system of matrix form；

Step 3.4：The algebra system is combined into into equality constraint least square problem with sampling equation.

4. method according to claim 3, it is characterised in that by the coefficient matrix in HASM equation group in the step 7 The concrete grammar for entering every trade piecemeal is：If each sub-block A of coefficient matrices A_iBe up to u rows, and A_i∈R^u×n, u ＜ ＜ n are Full rank, the i.e. row per block are independent, it is assumed that generating sub-block A_i, the sub-block is comprising j rows (the j ＜ in matrix A And known A u),_i ^TQR be decomposed into A_i ^T=Q_jR_j, block A_i ^TRow dependency estimationWherein r_hhIt is positioned at R_j(h, h) position element, a₁For a line in matrix A, orderAccording to A_i ^TQR decompose, update method is decomposed obtaining using QRQR decompose, then calculateRow Correlation estimation, if the estimation is still less than κ, a₁Add block A_i, wherein, κ is a positive number set in advance.

5. method according to claim 4, it is characterised in that the QR decomposes and adopts FrancisQR decomposition strategies.

6. the method according to claim 4 or 5, it is characterised in that to the block square after coefficient matrix decomposition in the step 7 The mode that battle array is stored is row compression storage mode CSR.

7. method according to claim 6, it is characterised in that the step 8 and step 9 are specifically included：Assume current changing Generation point is to judge x^kWhether restrain, if not restraining, find out from x^kOrthogonal chooses x apart from farthest block^kThrowing on described piece Shadow is used as next iteration point.