CN103678788B - Spatial data interpolation and curved surface fitting method based on surface theory - Google Patents

Abstract

The present invention provides a kind of spatial data interpolation based on surface theory and curved surface fitting method, including: will treat that fitting surface spatial spreading turns to grid points form, calculate possessive case site the first Equations of The Second Kind fundamental quantity coefficient andI=1,2；J=1,2；K=1,2, based on constraint function, set up object function according to the three of surface theory fundamental equations, treat the interior zone grid points of fitted area, it is desirable to minimize the finite difference quadratic sum of these three equation.With Optimization Solution device, Optimal Control Problem is solved, obtain the Digitized surface of matching；Judging whether gained Digitized surface meets required precision, if be unsatisfactory for, then iteration said process, until drawing the curved surface meeting required precision.Method provided by the invention, breaks away from traditional curved surface fitting method based on surface theory and treats the dependence of fitted area boundary value, it is possible to accurately and reliably obtain high-precision simulation curved surface.

Description

Spatial data interpolation and curved surface fitting method based on surface theory

Technical field

The present invention relates to learn a kind ofly, the curved surface modeling method in computer aided design and manufacture field, be specifically related to a kind of spatial data interpolation based on surface theory and curved surface fitting method.

Background technology

HASM (HighAccuracySurfaceModeling, High Accuracy Surface Modeling) is the curved surface modeling new method based on Differential Geometry surface theory proposed at the beginning of Chinese scholar 21 century.The basic thought of surface theory is: a curved surface is except its locus, and its shape is determined by first kind basic form and Equations of The Second Kind basic form.The exact value of two class basic form coefficients of curved surface is obtained, thus the shape obtaining a curved surface can be simulated by iterative computation.

At present, HASM method is mainly used in the research being simulated curved surface by discrete point.Existing HASH method mainly comprises the following steps: (1) needs sampled data to be interpolated calculating to obtain boundary value, wherein, interpolation method includes the method such as IDW, Spline (cubic spline interpolation) and Kriging (Kriging regression method)；Concrete, IDW (InverseDistanceWeighted) is a kind of conventional and the spatial interpolation methods of simplicity, principle is: be weighted on average with the distance between interpolation point and sample point for weight, from interpolation point more close to sample point give weight more big.(2) then, boundary value and other sampled datas being inputted object function jointly, by object function is resolved, simulation obtains curved surface.That is, boundary value is a necessary requirement in existing HASH method simulation curved surface process, and, the HASH precision simulating curved surface is had material impact effect by boundary value, when other conditions are identical, different boundary values can cause that low precision is from bigger simulation curved surface.Visible, boundary value is relied on relatively big by existing HASM method, thus hindering the development of HASM technology.

Further, patent 200710064595.7 " method that based on surface theory philosophy, terrain surface is set up mathematical model ", patent 201110021504.8 " curved surface modeling method based on surface theory and Optimal Control Theory ", number of patent application 201310002472.6 " High Accuracy Surface Modeling Intelligentized method and device " and document [1-15] are all without explicitly pointing out the scope treating fitting surface grid points, also without the concrete preparation method of zoning boundary value needed for pointing out curved surface modeling, and regional edge dividing value is bigger on the impact of 201310002472.6 methods and resultses, make 201310002472.6 and 201110021504.8 patent user in implementing process, there is very big ambiguity and uncertainty, HASM unstable result is increased, this patent overcomes these ambiguities and unstability, explicitly point out the scope treating matching grid points, explicitly point out the processing method of boundary value and initial value, the fitting result of this patent method has the feature such as stability and uniqueness.

As follows, for the details of document [1]-[16] published in prior art, these 16 sections of documents can be cited in the follow-up introduction of the present invention.

Yue Tianxiang, Du Zhengping, Song Dunjiang. High Accuracy Surface Modeling (HASM4) [J], Journal of Image and Graphics, 2006,12 (2): 343-348.

Tian-XiangYue, Zheng-PingDu, Dun-JiangSong, YunGong.2007.Anewmethodofsurfacemdinganditsapplicationto DEMconstruction.Geomorphology91 (1-2): 161-172.

W.J.Shi, J.Y.Liu, Z.P.Du, Y.J.Song, C.F.Chen, andT.X.Yue, SurfacemdingofsoilpH.Geoderma, 2009,150 (1-2), pp.113-119.

Tian-XiangYue, Dun-JiangSong, Zheng-PingDu, WeiWang.2010.High-accuracysurfacemdlinganditsapplication toDEMgeneration.InternationalJournalofRemoteSensing31 (8): 2205-2226.

Tian-XiangYue, Shi-HaiWang.2010.AdjustmentcomputationofHASM:ahigh-accur acyandhigh-speedmethod.InternationalJournalofGeographica lInformationScience24 (11): 1725-1743.

Tian-XiangYue, Chuan-FaChen, Bai-LianLi.2010.Anadaptivemethodofhighaccuracysurfacemdi nganditsapplicationtosimulatingelevationsurfaces.Transac tionsinGIS14 (5): 615-630.

Chen, C.F., Yue, T.X.:AmethodofDEMconstructionandrelatederroranalysis.Com put.Geosci.36,717-725 (2010).

Tian-XiangYue.2011.SurfaceMdling:HighAccuracyandHighSpee dMethods.NewYork:CRCPress (Taylor&Francisgroup).

Tian-XiangYue, Chuan-FAChen&Bai-LianLi (2012): Ahigh-accuracymethodforfillingvoidsanditsverification, InternationalJournalofRemoteSensing, 33:9,2815-2830.

Chen, C.F., Yue, T.X., Li, Y.Y.:AhighspeedmethodofSMTS.Comput.Geosci.41,64-71 (2012).

Chen, Chuan-Fa, Yue, Tian-Xiang, Dai, Hong-Lei, Tian, Mao-Yi.2013.ThemoothnessofHASM.InternationalJournalofGeo graphicalInformationScience27 (8), 1651-1667.

Chen, Chuan-Fa, YanyanLi, TianxiangYue.2013.SurfacemdingofDEMsbasedonasequentialad justmentmethod.InternationalJournalofGeographicalInforma tionScience27 (7), 1272-1291.

ShengxiangTang&DaijunJiang, ApplicableAnalysis (2013): Twoefficientalgorithmsforsurfaceconstruction, ApplicableAnalysis:AnInternationalJournal, DOI:10.1080/00036811.2013.781157.

Yue, T.-X., Zhao, N., Yang, H., Song, Y.-J., Du, Z.-P., Fan, Z.-M.andSong, D.-J. (2013), AMulti-GridMethodofHighAccuracySurfaceMdingandItsValidat ion.TransactionsinGIS.doi:10.1111/tgis.12019.

NaZhao, TianxiangYue, MingweiZhao, Animprovedversionofahighaccuracysurfacemdingmethod [J], GEM-InternationalJournalonGeomathematics, DOI10.1007/s13137-013-0051-z, 2013 (inpress).

Toponogov, V.A., 2006.Differentialgeometryofcurvesandsurfaces.NewYork:Bir khauser.

Summary of the invention

For the defect that prior art exists, the present invention provides a kind of spatial data interpolation based on surface theory and curved surface fitting method, breaks away from the conventional curved-surface approximating method based on surface theory and treats the dependence of fitted area boundary value, it is possible to obtains high-precision simulation curved surface.

The technical solution used in the present invention is as follows:

The present invention provides a kind of spatial data interpolation based on surface theory and curved surface fitting method, comprises the following steps:

S1, sets initial value matrix F 0 and the sampled data for the treatment of fitting surface；Set up constraint function；

S2, will treat that fitting surface place original area spatial spreading turns to grid points form, and wherein, grid points includes two classes: border grid points and internal grids point；Each border grid points composition borderline region, each internal grids point composition interior zone；

Then according to S1 obtain initial value matrix F 0 and sampled data, calculate the first kind fundamental quantity coefficient E of possessive case site, F, G, Equations of The Second Kind fundamental quantity coefficient L, M, N andI=1,2；J=1,2；K=1,2；Wherein, the area of described first kind fundamental quantity coefficient angle and terrain surface territory for representing both direction in the arc length of terrain surface upper curve, terrain surface；Described Equations of The Second Kind fundamental quantity coefficient is for portraying the bendability in terrain surface space；

S3, utilizes Algebra modeling instrument, sets up object function according to the three of surface theory fundamental equations, and described object function, described constraint function and described initial value matrix F 0 are passed to Optimization Solution device；Wherein, described object function is based upon described interior zone；

S4, described object function is optimized and solves by described Optimization Solution device, obtains the Digitized surface of analog form；

S5, described Optimization Solution device judges whether obtained described Digitized surface meets required precision, if it is satisfied, then directly export described Digitized surface, and process ends；If be unsatisfactory for, then iteration S2-S5, until drawing the Digitized surface meeting required precision, and export the described Digitized surface meeting required precision.

Preferably, in S1, described initial value matrix F 0 includes following several situation:

(1) described initial value matrix F 0 is set as null matrix；

(2) meansigma methods that the element in described initial value matrix F 0 is all sampled data is set.

Preferably, in S1, described constraint function is the equality constraints functions set up according to sampled data, it may be assumed that the constraint function of SIST is f_{M, n}The equality constraint of=height form, wherein m, n represent at grid points (x_m, y_n) there is sampled point at place, the property value of this sampled point is height.

Preferably, S2 specifically includes following steps:

S21, by treating that fitting surface place original area spatial spreading turns to the grid points rectangular array of r row c row, includes r*c grid points altogether, and wherein, line number is followed successively by: 0,1,2 ... r-1；Column number is followed successively by: 0,1,2 ... c-1；

S22, if curved surface expression formula is that (x, y), for any one grid points (x for s=f_i, y_j), wherein, i ∈ (0,1,2...r-1), j ∈ (0,1,2 ... c-1), calculate the E of this grid points, F, G-value by equation (1)；L, M, N value of this grid points is calculated by equation (2)；This grid points is calculated by equation (4)-(9)I=1,2；J=1,2；K=1,2[see reference document 16]；

$\left\{\begin{array}{c}E=1+{f_x}^2\\ G=1+{f_y}^2\\ F=f_x\·f_y\end{array}\right.,---\left(19\right)$

Wherein, f_xRepresent f (x, y) first-order partial derivative in the x direction；F_yRepresent f (x, y) first-order partial derivative in y-direction；

$\left\{\begin{array}{c}L=\frac{f_{\mathrm{xx}}}{\sqrt{1+{f_x}^2+{f_y}^2}}\\ N=\frac{f_{\mathrm{yy}}}{\sqrt{1+{f_x}^2+{f_y}^2}}\\ M=\frac{f_{\mathrm{xy}}}{\sqrt{1+{f_x}^2+{f_y}^2}}\end{array}\right.,---\left(20\right)$

Wherein, f_xx, f_yy, f_xyRepresent respectively f (x, y) to the second dervative of x, to the second dervative of y and to each first derivative of x, y；

For any one grid points (x_i, y_j), Gauss equation group is expressed as formula (3) [see reference document 16, p128]:

$\left\{\begin{array}{c}{f}_{\mathrm{xx}}={\mathrm{\Γ}}_{11}^1\·f_x+{\mathrm{\Γ}}_{11}^2\·f_y+L\·{(E\·G-F^2)}^{-\frac{1}{2}}\\ f_{\mathrm{yy}}={\mathrm{\Γ}}_{22}^1\·f_x+{\mathrm{\Γ}}_{22}^2\·f_y+L\·{(E\·G-F^2)}^{-\frac{1}{2}}\\ f_{\mathrm{xy}}={\mathrm{\Γ}}_{12}^1\·f_x+{\mathrm{\Γ}}_{12}^2\·f_y+L\·{(E\·G-F^2)}^{-\frac{1}{2}}\end{array}\right.,---\left(21\right)$

Then:

${\mathrm{\Γ}}_{11}^1=\frac{1}{2}(G\·E_x-2F\·F_x+F\·E_y)\·{(E\·G-F^2)}^{-1}---\left(22\right)$

${\mathrm{\Γ}}_{12}^1=\frac{1}{2}(G\·E_y-F\·G_x)\·{(E\·G-F^2)}^{-1}---\left(23\right)$

${\mathrm{\Γ}}_{22}^1=\frac{1}{2}(2G\·F_y-G\·G_x+F\·G_y)\·{(E\·G-F^2)}^{-1}---\left(24\right)$

${\mathrm{\Γ}}_{11}^2=\frac{1}{2}(2E\·F_x-E\·E_y-F\·E_x)\·{(E\·G-F^2)}^{-1}---\left(25\right)$

${\mathrm{\Γ}}_{12}^2=\frac{1}{2}(E\·G_x-F\·E_y)\·{(E\·G-F^2)}^{-1}---\left(26\right)$

${\mathrm{\Γ}}_{22}^2=\frac{1}{2}(E\·G_y-2F\·F_y+F\·G_x)\·{(E\·G-F^2)}^{-1},---\left(27\right)$

Wherein, E_x、E_y、F_x、F_y、G_x、G_yRepresent that x, y are sought partial derivative by E, F, G respectively.

Preferably, in S3, utilize Algebra modeling instrument, set up object function according to the three of surface theory fundamental equations particularly as follows:

S31, by the equal sign left end centered finite difference formal expansion of formula (3), then is multiplied by h at equal sign two ends simultaneously², and keep primitive form constant as constant term formula (3) right-hand member, using formula (3) left end as variable, then formula (3) is re-expressed as:

$\left\{\begin{array}{c}{f}_{i+1,j}-2f_{i,j}+f_{i-1,j}=({\mathrm{\Γ}}_{11}^1\·f_x+{\mathrm{\Γ}}_{11}^2\·f_y+L\·{(E\·G-F^2)}^{-\frac{1}{2}})\·h^2\\ f_{i,j+1}-2f_{i,j}+f_{i,j-1}=({\mathrm{\Γ}}_{22}^1\·f_x+{\mathrm{\Γ}}_{22}^2\·f_y+N\·{(E\·G-F^2)}^{-\frac{1}{2}})\·h^2\\ f_{i+1,j+1}-2f_{i-1,j+1}+f_{i-1,j-1}-=f_{i+1,j-1}({\mathrm{\Γ}}_{12}^1\·f_x+{\mathrm{\Γ}}_{12}^2\·f_y+M\·{(E\·G-F^2)}^{-\frac{1}{2}})\·h^2\end{array}\right.,---\left(28\right)$

Wherein, h is material calculation, i.e. pixel resolution；F_{I, j}(x, y) at grid points (x to represent f_i, y_i) value at place；

S32, is abbreviated as RHS1, RHS2 and RHS3 respectively by the right-hand vector of formula (10), right-hand member moves to left end simultaneously, obtains below equation expression formula:

$\left\{\begin{array}{c}\mathrm{ex}1=f_{i+1,j}-2f_{i,j}+f_{i-1,j}-{\mathrm{RHS}1}_{i,j}\\ \mathrm{ex}2=f_{i,j+1}-2f_{i,j}+f_{i,j-1}-{\mathrm{RHS}2}_{i,j}\\ \mathrm{ex}3=f_{i+1,j+1}-2f_{i-1,j+1}+f_{i-1,j-1}-f_{i+1,j-1}-{\mathrm{RHS}3}_{i,j}\end{array}\right.,---\left(29\right)$

S33, in formula (11), the expression formula that minimizes of the quadratic sum of three expression formulas is the object function of SIST method, and wherein, the span of i is 1,2 ..., r-2；The span of j is 1,2 ..., c-2；Being write as the form of mathematical function, the object function of SIST method is expressed as:

$\min \left\{\underset{i=1}{\overset{r-2}{\mathrm{\Σ}}}{\underset{j=1}{\overset{c-2}{\mathrm{\Σ}}}}_{+{(f_{i+1,j+1}-f_{i-1,j+1}+f_{i-1,j-1}-f_{i+1,j-1}-{\mathrm{RHS}3}_{i,j})}^2}^{{(f_{i+1,j}-{2f}_{i,j}+f_{i-1,j}-{\mathrm{RHS}1}_{i,j})}^2+{(f_{i,j+1}-{2f}_{i,j}+f_{i,j-1}-{\mathrm{RHS}2}_{i,j})}^2}\right\}$ 。

Preferably, also include:

If z=is (f_0,0, f_0,1..., f_{0, c-1}, f_1,0, f_1,1..., f_{I, 0}..., f_{I, c-1}..., f_{R-1, c-1})^T, then equation (10) is expressed as following matrix form:

$\left\{\begin{array}{c}\mathrm{Az}=\mathrm{RHS}1\\ \mathrm{Bz}=\mathrm{RHS}2\\ \mathrm{Cz}=\mathrm{RHS}3\end{array}\right.,---\left(30\right)$

A, B and C are coefficient matrix, express in the following ways:

Definition Z_{M, n}Represent a null matrix, the line number of this null matrix and columns respectively m and n, then: coefficient matrices A is expressed as:

Wherein, A_{C, (r-2) * c}Represent the matrix of a line number and columns respectively c and (r-2) * c, it may be assumed that

Coefficient matrix B is expressed as:

Wherein, B_{C, (r-2) * c}Represent the matrix of a line number and columns respectively c and (r-2) * c, be expressed as:

Coefficient matrix C is:

Wherein, C_{C, (r-2) * c}Represent the matrix of a line number and columns respectively c and (r-2) * c, be expressed as

Preferably, in S4, described Optimization Solution device described object function is optimized solve particularly as follows:

Under the constraints of described constraint function, described object function is optimized and solves.

Preferably, in S5, described Optimization Solution device judges whether obtained described Digitized surface meets required precision particularly as follows: assume twice curve modeling result respectively two the column vector z in front and back_iAnd z_i-1, then when | | z_i-z_i+1||₂/||z_i||₂During≤delta, show that precision meets the conclusion of requirement；Wherein, delta span is: equal to or less than the numeral of 1e-10.

Spatial data interpolation based on surface theory provided by the invention and curved surface fitting method, have the advantage that

(1) SIST method relaxes the requirement of boundary value, decreases the dependence to other interpolation methods.Further, specify that the acquisition mode of initial value.

(2) SIST regards Fitting as the Optimal Control Problem of one standard, and three equations are only set up in the interior zone grid points treating fitted area, improve space interpolation and the reliability of surface fitting result.

Accompanying drawing explanation

Fig. 1 is the schematic flow sheet of the spatial data interpolation based on surface theory provided by the invention and curved surface fitting method；

Fig. 2 is the schematic diagram of the discrete grid points rectangular array turning to r row c row；

Fig. 3 is be original surface figure when 40 for contrasting the sampling interval of Rosenbrock mathematical surface；

Fig. 4 is the sampling interval uses when being 40 HASM4 method to simulate the Rosenbrock mathematical surface figure obtained；

Fig. 5 is the sampling interval uses when being 40 SIST method to simulate the Rosenbrock mathematical surface figure obtained.

Fig. 6 is be original surface figure when 20 for contrasting the sampling interval of Rosenbrock mathematical surface；

Fig. 7 is the sampling interval uses when being 20 HASM4 method to simulate the Rosenbrock mathematical surface figure obtained；

Fig. 8 is the sampling interval uses when being 20 SIST method to simulate the Rosenbrock mathematical surface figure obtained；

Fig. 9 is be original surface figure when 80 for contrasting the sampling interval of F1 mathematical surface；

Figure 10 is the sampling interval uses when being 80 HASM4 method to simulate the F1 mathematical surface figure obtained；

Figure 11 is the sampling interval uses when being 80 SIST method to simulate the F1 mathematical surface figure obtained；

Figure 12 is be original surface figure when 40 for contrasting the sampling interval of F1 mathematical surface；

Figure 13 is the sampling interval uses when being 40 HASM4 method to simulate the F1 mathematical surface figure obtained；

Figure 14 is the sampling interval uses when being 40 SIST method to simulate the F1 mathematical surface figure obtained；

Figure 15 uses SIST to be fitted the surface chart obtained after obtaining initial value by IDW Spatial Interpolation Method；

Figure 16 uses SIST to be fitted the surface chart obtained after obtaining initial value by Kriging Spatial Interpolation Method；

Figure 17 uses SIST to be fitted the surface chart obtained after obtaining initial value by Spline Spatial Interpolation Method；

Figure 18 is that initial value adopts null matrix to use SIST to be fitted the surface chart obtained.

Detailed description of the invention

Below in conjunction with accompanying drawing, the present invention is described in detail:

For convenience of description, by the spatial data interpolation based on surface theory provided by the invention and curved surface fitting method referred to as SIST method.The present invention is as without particularly pointing out, the alleged series methods that HASM method is the HASM difference variants such as HASM4, HASM5 is referred to as.Wherein, HASM4 only considers two equations of surface theory, and HASM5 considers 3 equations of surface theory.

As it is shown in figure 1, the present invention provides a kind of spatial data interpolation based on surface theory and curved surface fitting method, comprise the following steps:

S1, sets the initial value matrix F 0 treating fitting surface and obtains sampled data；Set up constraint function；

Initial value matrix F 0 includes following several situation:

(1) initial value matrix F 0 is set as null matrix；

(2) meansigma methods that the element in described initial value matrix F 0 is all sampled data is set.

Patent 200710064595.7 " method that based on surface theory philosophy, terrain surface is set up mathematical model ", patent 201110021504.8 " curved surface modeling method based on surface theory and Optimal Control Theory ", number of patent application 201310002472.6 " High Accuracy Surface Modeling Intelligentized method and device " and document [1-12] all do not explicitly point out curved surface modeling needed for the concrete preparation method of zoning initial value, and the fitting result of these methods is affected bigger by initial value, these articles existing and implementing of patented method user is made still to have very big difficulty, very big ambiguity and uncertainty occur.

And SIST method of the present invention is also required to initial value, further, specify that the preparation method of initial value, it may be assumed that if user does not provide initial value, SIST can utilize the meansigma methods of sampled data to carry out surface fitting as initial value automatically, or adopts null matrix as initial value.The results showed, when initial value adopts the meansigma methods that null matrix still adopts sampled data, the precision of surface fitting result is affected less.

Additionally, it is emphasized that the particular type solving the constraint function that object function adopts is not limiting as by the present invention, it is possible to for conventional any kind of constraint function, such as: can be the equality constraints functions set up according to sampled data, it may be assumed that the constraint function of SIST is f_{M, n}The equality constraint of=heigt form, wherein m, n represent at grid points (x_m, y_n) there is sampled point at place, the property value of this sampled point is height.

S2, will treat that fitting surface place original area spatial spreading turns to grid points form, and wherein, grid points includes two classes: border grid points and internal grids point；Each border grid points composition borderline region, each internal grids point composition interior zone；

Then according to S1 obtain initial value matrix F 0 and sampled data, calculate the first kind fundamental quantity coefficient E of possessive case site, F, G, Equations of The Second Kind fundamental quantity coefficient L, M, N andI=1,2；J=1,2；K=1,2；Wherein, the area of described first kind fundamental quantity coefficient angle and terrain surface territory for representing both direction in the arc length of terrain surface upper curve, terrain surface；Described Equations of The Second Kind fundamental quantity coefficient is for portraying the bendability in terrain surface space；

This step specifically includes following steps:

S21, by treating that fitting surface place original area spatial spreading turns to the grid points rectangular array of r row c row, includes (r*c) individual grid points altogether, and wherein, line number is followed successively by: 0,1,2 ... foretell 1；Column number is followed successively by: 0,1,2 ... c-1；As in figure 2 it is shown, be the schematic diagram of the discrete grid points rectangular array turning to r row c row.In fig. 2, the grid points drawing slash is border grid points, forms borderline region；Other do not draw the grid points composition interior zone of slash.

S22, if curved surface expression formula is that (x, y), for any one grid points (x for s=f_i, y_j), wherein, i ∈ (0,1,2 ... r-1), j ∈ (0,1,2 ... c-1), calculate the E of this grid points, F, G-value by equation (1)；L, M, N value of this grid points is calculated by equation (2)；This grid points is calculated by equation (4)-(9)I=1,2；J=1,2；K=1,2[see reference document 16]；

$\left\{\begin{array}{c}E=1+{f_x}^2\\ G=1+{f_y}^2\\ F=f_x\·f_y\end{array}\right.,---\left(37\right)$

Wherein, f_xRepresent f (x, y) first-order partial derivative in the x direction；F_yRepresent f (x, y) first-order partial derivative in y-direction；

$\left\{\begin{array}{c}L=\frac{f_{\mathrm{xx}}}{\sqrt{1+{f_x}^2+{f_y}^2}}\\ N=\frac{f_{\mathrm{yy}}}{\sqrt{1+{f_x}^2+{f_y}^2}}\\ M=\frac{f_{\mathrm{xy}}}{\sqrt{1+{f_x}^2+{f_y}^2}}\end{array}\right.,---\left(38\right)$

Wherein, f_xx, f_yy, f_xyRepresent respectively f (x, y) to the second dervative of x, to the second dervative of y and to cutter, each first derivative of y；

For any one grid points (x_i, y_i), Gauss equation group is expressed as formula (3) [see reference document 16, p128]:

$\left\{\begin{array}{c}{f}_{\mathrm{xx}}={\mathrm{\Γ}}_{11}^1\·f_x+{\mathrm{\Γ}}_{11}^2\·f_y+L\·{(E\·G-F^2)}^{-\frac{1}{2}}\\ f_{\mathrm{yy}}={\mathrm{\Γ}}_{22}^1\·f_x+{\mathrm{\Γ}}_{22}^2\·f_y+L\·{(E\·G-F^2)}^{-\frac{1}{2}}\\ f_{\mathrm{xy}}={\mathrm{\Γ}}_{12}^1\·f_x+{\mathrm{\Γ}}_{12}^2\·f_y+L\·{(E\·G-F^2)}^{-\frac{1}{2}}\end{array}\right.,---\left(39\right)$

Then:

${\mathrm{\Γ}}_{11}^1=\frac{1}{2}(G\·E_x-2F\·F_x+F\·E_y)\·{(E\·G-F^2)}^{-1}---\left(40\right)$

${\mathrm{\Γ}}_{12}^1=\frac{1}{2}(G\·E_y-F\·G_x)\·{(E\·G-F^2)}^{-1}---\left(41\right)$

${\mathrm{\Γ}}_{22}^1=\frac{1}{2}(2G\·F_y-G\·G_x+F\·G_y)\·{(E\·G-F^2)}^{-1}---\left(42\right)$

${\mathrm{\Γ}}_{11}^2=\frac{1}{2}(2E\·F_x-E\·E_y-F\·E_x)\·{(E\·G-F^2)}^{-1}---\left(43\right)$

${\mathrm{\Γ}}_{12}^2=\frac{1}{2}(E\·G_x-F\·E_y)\·{(E\·G-F^2)}^{-1}---\left(44\right)$

${\mathrm{\Γ}}_{22}^2=\frac{1}{2}(E\·G_y-2F\·F_y+F\·G_x)\·{(E\·G-F^2)}^{-1}---\left(45\right)$

Wherein, E_x、E_y、F_x、F_y、G_x、G_yRepresent that x, y are sought partial derivative by E, F, G respectively.

S3, utilizes Algebra modeling instrument, sets up object function according to the three of surface theory fundamental equations, and described object function, described constraint function and described initial value matrix F 0 are passed to Optimization Solution device；Wherein, object function is based upon interior zone, it may be assumed that in object function expression formula, and the spatial dimension of variable is the grid points of interior zone.

Algebra modeling instrument can adopt the industrial conventional modeling tool such as GAMS (TheGeneralAlgebraicModelingSystem), AMPL (AMathematicalProgrammingLanguage).

Solver can adopt CPLEX, KNITRO, IPOPT, MINOS, SNOPT etc. to have secondary or the solver of nonlinear problem Optimization Solution function.

This step particularly as follows:

S31, by the equal sign left end centered finite difference formal expansion of formula (3), then is multiplied by h at equal sign two ends simultaneously², and keep primitive form constant as constant term formula (3) right-hand member, using formula (3) left end as variable, then formula (3) is re-expressed as:

$\left\{\begin{array}{c}{f}_{i+1,j}-2f_{i,j}+f_{i-1,j}=({\mathrm{\Γ}}_{11}^1\·f_x+{\mathrm{\Γ}}_{11}^2\·f_y+L\·{(E\·G-F^2)}^{-\frac{1}{2}})\·h^2\\ f_{i,j+1}-2f_{i,j}+f_{i,j-1}=({\mathrm{\Γ}}_{22}^1\·f_x+{\mathrm{\Γ}}_{22}^2\·f_y+N\·{(E\·G-F^2)}^{-\frac{1}{2}})\·h^2\\ f_{i+1,j+1}-2f_{i-1,j+1}+f_{i-1,j-1}-=f_{i+1,j-1}({\mathrm{\Γ}}_{12}^1\·f_x+{\mathrm{\Γ}}_{12}^2\·f_y+M\·{(E\·G-F^2)}^{-\frac{1}{2}})\·h^2\end{array}\right.,---\left(46\right)$

Wherein, h is material calculation, i.e. pixel resolution；F_{I, j}(x, y) at grid points (x to represent f_i, y_j) value at place；

S32, is abbreviated as RHS1, RHS2 and RHS3 respectively by the right-hand vector of formula (10), right-hand member moves to left end simultaneously, obtains below equation expression formula:

$\left\{\begin{array}{c}\mathrm{ex}1=f_{i+1,j}-2f_{i,j}+f_{i-1,j}-{\mathrm{RHS}1}_{i,j}\\ \mathrm{ex}2=f_{i,j+1}-2f_{i,j}+f_{i,j-1}-{\mathrm{RHS}2}_{i,j}\\ \mathrm{ex}3=f_{i+1,j+1}-2f_{i-1,j+1}+f_{i-1,j-1}-f_{i+1,j-1}-{\mathrm{RHS}3}_{i,j}\end{array}\right.,---\left(47\right)$

S33, in formula (11), the expression formula that minimizes of the quadratic sum of three expression formulas is the object function of SIST method, and wherein, the span of i is 1,2 ..., r-2；The span of j is 1,2 ..., C-2.Being write as the form of mathematical function, the object function of SIST method is expressed as:

$\min \left\{\underset{i=1}{\overset{r-2}{\mathrm{\Σ}}}\underset{j=4}{\overset{c-2}{\mathrm{\Σ}}}\begin{array}{c}{(f_{i+1,j}-2j_{i,j}+j_{i-1,j}-{\mathrm{RHS}1}_{i,j})}^2\\ +{(f_{i,j+1}-2f_{i,j}+f_{i,j-1}-{\mathrm{RHS}2}_{i,j})}^2\\ +{(f_{i+1,j+1}-f_{i-1,j+1}+f_{i-1,j-1}-f_{i+1,j-1}-{\mathrm{RHS}3}_{i,j})}^2\end{array}\right\};$

The constraint function of SIST is f_{M, n}The equality constraint of=height form, wherein m, n represent at grid points (x_m, y_n) there is sampled point at place, the property value of this sampled point is height.

S4, described object function is optimized and solves by described Optimization Solution device, obtains the Digitized surface of analog form；Wherein, Optimization Solution device can adopt the industrial conventional Optimization Solution device such as CPLEX, KNITRO, IPOPT, MINOS, SNOPT.Wherein, Algebra modeling instrument is used to call these Optimization Solution devices, it is possible to complicated traffic issues shows as Mathematical Planning (MathematicProgramming) model of standard.Industrial conventional Optimization Solution device all achieves the function of parallel computation, in multi-CPU hardware environment, it is possible to accelerate space interpolation and the surface fitting process of SIST.

This step particularly as follows:

Under the constraints of described constraint function, described object function is optimized and solves.

S5, described Optimization Solution device judges whether obtained described Digitized surface meets required precision, if it is satisfied, then directly export described Digitized surface, and process ends；If be unsatisfactory for, then iteration S2-S5, until drawing the Digitized surface meeting required precision, and export the described Digitized surface meeting required precision.

Wherein it is possible to twice curve modeling result respectively two column vector z before and after assuming_iAnd z_i+1, then when | | z_i-z_i+1||₂/||z_i||₂During≤delta, show that precision meets the conclusion of requirement；Wherein, delta span is: equal to or less than the numeral of 1e-10.

The separation of process is just as professional division, and generally, professional division is more thin, is more conducive to the development of this specialty.The process that realizes of SIST adopts independent Algebra modeling instrument, and modeling process and Optimization Solution process are separated, and the model of foundation can export as .nl file, be conducive to upgrading and the maintenance of SIST.Assuming that z=(f_0,0, f_0,1..., f_{0, c-1}, f_1,0, f_1,1..., f_{I, 0}..., f_{I, c-1}..., f_{R-1, c-1})^T, then the dependent equation (10) of SIST method can be expressed as following matrix form:

$\left\{\begin{array}{c}\mathrm{Az}=\mathrm{RHS}1\\ \mathrm{Bz}=\mathrm{RHS}2\\ \mathrm{Cz}=\mathrm{RHS}3\end{array}\right.,---\left(48\right)$

Right-hand vector RHS1, RHS2, the RHS3 of formula (12) expresses identical object with above-mentioned RHS1, RHS2, RHS3.Difference between SIST method described in detail below and the dependent equation of HASM method:

(1) coefficient matrices A

First Z is defined_{M, n}Represent a null matrix, the line number of this null matrix and columns respectively m and n.Coefficient matrices A in the equation (10) of SIST can be expressed as:

Wherein, A_{C, (r-2) * c}Represent the matrix of a line number and columns respectively c and (r-2) * c, be represented by

The HASM's corresponding with the coefficient matrices A of SIST[deriving from document 8, P35-37, lower same] is:

It should be noted that I here_c-2It is the unit square formation on (c-2) rank, HASM'sIt is the square formation on (r-2) * (c-2) rank, uses hereRepresenting the coefficient matrix of HASM dependent equation, be to represent that HASM and SIST method dependent equation coefficient matrix has some similarities, but both structures are different, vary in size, the number of non-zero element, position also differ.

(2) coefficient matrix B

The coefficient matrix B of SIST method dependent equation (10) is

Wherein, B_{C, (r-2) * c}Represent the matrix of a line number and columns respectively c and (r-2) * c, be represented by:

HASM's[deriving from document 8, P35-37, lower same] is

Wherein, B_c-2It is the reflection matrix on (c-2) rank:

It is pointed out that HASM'sBe one by (r-2) individual B_c-2The square formation of composition, therefore be also the square formation on (r-2) * (c-2) rank, use hereRepresenting the coefficient matrix of HASM dependent equation, be to represent that HASM and SIST method dependent equation coefficient matrix has some similarities, but both matrix structures are different, vary in size, the number of non-zero element, position also differ.

(3) coefficient matrix C

The coefficient matrix C of the equation (10) of SIST method is

Wherein, C_{C, (r-2) * c}Represent a line number and columns respectively c and the matrix of (foretelling 2) * c, be represented by

HASM's[deriving from document 8, P35-37, lower same] is

Wherein, C_c-2It it is the reflection matrix on (c-2) rank

It is pointed out that HASM'sIt is the square formation on (r-2) * (c-2) rank, uses hereRepresenting the coefficient matrix of HASM dependent equation, be to represent that HASM and SIST method dependent equation coefficient matrix has some similarities, but both structures are different, vary in size, the number of non-zero element, position also differ.

It addition, in the equation group that obtains of HASM method, its coefficient matrix is symmetric positive definite sparse matrix, and in SIST method of the present invention, coefficient matrix is not symmetric positive definite matrix.Only when user provides accurate boundary value, SIST becomes identical with HASM, and the coefficient matrix that SIST equation is formed is only symmetric positive definite.

It can thus be seen that SIST method and HASM method difference mainly have 2 points: 1) z in SIST method is the column vector of a r*c row, and HASMColumn vector for (r-2) * (c-2) row；2) coefficient matrices A of SIST equation (10), coefficient matrix that B, C and HASM are correspondingDifference, the coefficient matrices A of SIST, B, C are the square formations on r*c rank, and the coefficient matrix of HASMIt it is all the square formation on (r-2) * (c-2) rank.

From the aforegoing it can be seen that the coefficient matrix of the coefficient matrices A of SIST, B, C and HASMThere is a lot of difference, for identical sampled data, not only matrix size is different, in matrix, the distributing position of nonzero element also has very big difference, due in the coefficient matrices A of SIST, B, C except containing 0 row vector, the number of all the other each row nonzero elements is just as, and the coefficient matrix more corresponding than HASM is more regular, thus SIST to realize process simpler than HASM.

The difference of SIST and HASM is also embodied in use procedure, SIST method does not need SIST user and provides the property value of borderline region, and HASM method needs HASM user first to provide the property value of borderline region, and in the article delivered and patent, HASM boundary value acquisition methods is not all discussed in detail, but analyze the coefficient matrix of HASM, it is known that boundary value is simply simply processed by really as given value, therefore, HASM boundary value acquisition methods has the feature of ambiguity, causes that HASM realizes the uncertainty of process.Similar with SIST, HASM is also based on the basic thought of surface theory, but HASM needs the boundary information that the offer of HASM user is extra.HASM method needs HASM4 user to provide the property value of borderline region, therefore the actual fitted area of HASM is inner area, SIST method provides the property value of borderline region without SIST user, that is, SIST need not obtain boundary value by the pre-interpolation method such as kriging, IDW, and historical facts or anecdotes border fitted area is whole region.HASM and SIST also differ in that: HASM4 only employs the first two equation in equation (9), although HASM5 employs whole three equations in equation (9), but assume F=0 by force, SIST employs whole three equations of equation (9), it does not have any assume by force.

In sum, the coefficient matrix of the coefficient matrices A of SIST, B, C and HASMAbove-mentioned difference, SIST method is different from HASM method the most significantly one of mark just, based on mathematical theory, as well known to those skilled in the art, the difference of coefficient matrix is different with process and cause by the approximating method of the two just, it is all very big difference that coefficient matrix differs one 0, or, it is all very big improvement that certain parameter adds numeral 1；Additionally, SIST method is compared with HASM method, although [surface theory is shown in document 16 to be all based on surface theory equally, P128] basic thought, but concrete principle has a lot of difference with the process of realization, and as all different in coefficient matrix size, coefficient matrix structure etc., the use procedure of SIST is also clearly distinguishable from HASM method, SIST need not provide boundary value by SIST user, and HASM needs HASM user first by other software (such as ARCGIS) interpolation to provide boundary value；Additionally, multiple contrast tests below also all demonstrate, SIST method surface fitting precision is apparently higher than HASM method, it is seen then that the new method that SIST method is a kind of surface fitting that the present invention proposes.

Contrast test 1

This experimental example is for comparing HASM4 and the SIST surface fitting precision for Rosenbrock mathematical surface.

Boundary value and initial value that HASM4 uses derive from IDW (InverseDistanceWeight) method in ArcGIS10.1, take IDW default parameter value.To even things up, the initial value of SIST also uses the initial value that IDW method obtains.

Use the precision of Rosenbrock function test HASM4 and SIST.Wherein, Rosenbrock function is a non-convex function for test optimization algorithm performance, and span is [0,2509].The mathematic(al) representation of Rosenbrock function is: s=(1-x)²+100(y-x²)², wherein ,-2≤x≤2 ,-1≤y≤3.

As it is shown on figure 3, be original surface when 40, as can be seen from Figure 3 sampled data for the sampling interval, wherein, original grid cell number is 401*401.Fig. 4 is the sampling interval uses when being 40 HASM4 method to simulate the Rosenbrock mathematical surface obtained；Fig. 5 is the sampling interval uses when being 40 SIST method to simulate the Rosenbrock mathematical surface obtained.

As shown in Figure 6, it is original surface when 20, as can be seen from Figure 6 sampled data for the sampling interval.Fig. 7 is the sampling interval uses when being 20 HASM4 method to simulate the Rosenbrock mathematical surface obtained；Fig. 8 is the sampling interval uses when being 20 SIST method to simulate the Rosenbrock mathematical surface obtained.

Fig. 4 and Fig. 5 is compared, Fig. 7 and Fig. 8 is compared, it can be seen that SIST analog result is closer to real mathematical surface, and HASM4 is very big in boundary deformation, and two kinds of methods are less in internal difference.

It addition, adopt following formula (12) to calculate the root-mean-square error RMSE of different sampled point number and SIST and the HASM4 under the sampling interval respectively, result is in Table 1.

$\mathrm{RMSE}=\sqrt{\frac{\mathrm{\Σ}{(F_{i,j}-\stackrel{\‾}{f_{i,j}})}^2}{r*c}},---\left(12\right)$

Wherein, f_{I, j}Represent real curved surface property value,Representing the result of surface fitting, r, c represent the line number and columns for the treatment of matching Discrete Surfaces respectively.

Table 1RMSE statistical table (original grid cell number 401*401)

Sequence number	Sampling interval	Sampled point number	HASM4	SIST
					1	40	121	381.36	51.75
2	20	441	288.07	20.65
					3	10	1681	216.54	7.90
4	5	6561	167.15	2.41

From the error statistics result shown in table 1 it can be seen that SIST fitting surface precision is better than HASM4.

Contrast test 2

This practicality example is used for comparing HASM4 and the SIST surface fitting precision for F1 mathematical surface [curved surface derives from document 2].

Boundary value and initial value that HASM4 uses derive from IDW (InverseDistanceWeight) method in ArcGIS10.1, take IDW default parameter value.The initial value of SIST also uses IDW.

As it is shown in figure 9, be original surface when 80, as can be seen from Figure 9 sampled data for the sampling interval, wherein, original grid cell number is 401*401.Figure 10 is the sampling interval uses when being 80 HASM4 method to simulate the F1 mathematical surface obtained；Figure 11 is the sampling interval uses when being 80 SIST method to simulate the F1 mathematical surface obtained.

As shown in figure 12, being original surface when 40, as can be seen from Figure 12 sampled data for the sampling interval, wherein, original grid cell number is 401*401.Figure 13 is the sampling interval uses when being 40 HASM4 method to simulate the F1 mathematical surface obtained；Figure 14 is the sampling interval uses when being 40 SIST method to simulate the F1 mathematical surface obtained.

Figure 10 and Figure 11 is compared, Figure 13 and Figure 14 is compared, it can be seen that SIST analog result is closer to real mathematical surface, and HASM4 is very big in boundary deformation, and two kinds of methods are less in internal difference.

It addition, the formula of employing (12) calculates the root-mean-square error RMSE of different sampled point number and SIST and the HASM4 under the sampling interval respectively, result is in Table 2.

Table 2RMSE statistical table (original grid cell number 401*401)

Sequence number	Sampling interval	Sampled point number	HASM4	SIST
					1	100	25	1.33	0.66
2	80	36	1.00	0.42
					3	50	81	0.72	0.31
4	40	121	0.65	0.24

From the error statistics result shown in table 2 it can be seen that SIST fitting surface precision is better than HASM4.

Contrast test 3

This practicality example is used for comparing HASM4 and the SIST surface fitting precision for F2-F8 mathematical surface [curved surface derives from document 2].

Employing formula (12) calculates the root-mean-square error RMSE of different sampled point number and SIST and the HASM4 under the sampling interval respectively, and result is in Table 3.

Table 3RMSE statistical table (original grid cell number 401*401)

As can be seen from the above table, in most cases, the surface fitting precision of SIST method is above the surface fitting precision of HASM4 method；Only in limited instances, SIST precision relatively HASM4 is low, it is possible to reason is: the boundary value relatively actual value that IDW method interpolation obtains.

Contrast test 4

This test is for the terrain surface of matching Jiangxi Province Ji'an City soil pH value, and this experiment is for illustrating the stability of SIST method, and the impact by initial value is only small.

Study area is typical steep hilly region in southern China, is positioned at the Ji'an City area under one's jurisdiction of Central Jiangxi Province, Jian County and Taihe County, gross area 6156.92km².Weather is middle subtropical zone monsoon climate, the average annual precipitation of Ji'an City area under one's jurisdiction, Jian County and Taihe County respectively 1458mm, 1438mm and 1381mm, average annual temperature respectively 18.1 DEG C, 18.4 DEG C and 19.0 DEG C.The elevation of study area tapers into from four Zhou Dao centers, ranges for 1204.5～42.0m.

Spatial distribution state according to the In The Soils of study area, land use pattern and parent rock, take into account the factors such as vegetation, landform and transportation condition, in October, 2007, in study area, gather 150,0～20em topsoil sample, record the longitude and latitude of sampled point, elevation, soil types, land use pattern, vegetation etc. simultaneously.Determine soil quick-effective phosphor AP, lithium Li, chromium Cr, pH, alkali-hydrolyzable nitrogen AN and full potassium TK.Wherein, available P adopts HCl-H₂SO₄Lixiviate-molybdenum antimony resistance colorimetric method measures；Li, Cr adopt HF-HNO₃-HClO₄Nitrated, ICP-OES measures；The soil ratio measuring employing 1:2.5 of soil pH, pH meter measures.This experiment mainly utilizes the space interpolation of the pH value of 150 soil sampling points to test.

As shown in figure 15, for using SIST to be fitted the surface chart obtained after obtaining initial value by IDW Spatial Interpolation Method；

As shown in figure 16, for using SIST to be fitted the surface chart obtained after obtaining initial value by Kriging Spatial Interpolation Method；

As shown in figure 17, for using SIST to be fitted the surface chart obtained after obtaining initial value by Spline Spatial Interpolation Method；

As shown in figure 18, null matrix is adopted to use SIST to be fitted the surface chart obtained for initial value.

Comparison diagram 15-Figure 18, it can be seen that the surface fitting precision of Figure 15-Figure 18 is closely, it is seen then that when adopting SIST method provided by the invention to carry out surface fitting, the dependence of initial value is less, has higher stability and reliability.

In sum, the spatial data interpolation based on surface theory provided by the invention and curved surface fitting method, Fitting is converted to an Optimal Control Problem, constraint function derives from sampled data, and object function is the quadratic sum minimizing three equations.In object function expression formula, what variable related to ranges for the inside in region undetermined.Have the advantage that

(1) SIST method relaxes the requirement of boundary value, decreases the dependence to other interpolation methods, but, if user is aware of when comparing exact boundry Value Data, it is possible to using boundary value as sampled data, in input SIST.

(2) processing method of initial value: the HASM4 of existing patent and document is required for initial value, but existing patent and document are this details is not carried, what the very difficult initial value learning HASM4 and correlation technique from these documents and patent is, the HASM4 method user in existing patent and document is made to be difficult to select or use correct initial value so that the result uncertainty of HASM4 surface fitting is very big.The fit procedure of SIST is also required to initial value, and initial value affects SIST and calculates the time that process needs, but the impact of result of calculation is only small.If SIST user does not provide initial value, SIST can utilize the meansigma methods of sampled data to carry out surface fitting as initial value automatically, if the user of service of this method provides comparatively ideal initial value, then SIST can use the initial value initial value as SIST of SIST user offer, reduces the operation time of SIST.

(3) with the addition of the 3rd equation in object function.With the addition of the cross term equation of surface theory equation, and suppose that F is not equal to 0, be different from the hypothesis of document [15].

(4) Fitting is regarded as the Optimal Control Problem in an operational research field by SIST, to treat that whole grid points of fitting surface are as variable undetermined, the equality constraint of sampled data is used as constraint function, equations turned for object function by the three of surface theory, treat the interior zone grid points of fitted area, setting up the finite difference form of these three equation respectively, the target of SIST is exactly the quadratic sum minimizing these three equation.This details is not carried or is not discussed in detail by the HASM4 method of existing patent and document, it is difficult to learn the equation of HASM4 and correlation technique for which grid points is set up from these documents and patent, the HASM4 method user in existing patent and document is made to be difficulty with correct HASM4 so that the result uncertainty of HASM4 surface fitting is very big.

The above is only the preferred embodiment of the present invention; it should be pointed out that, for those skilled in the art, under the premise without departing from the principles of the invention; can also making some improvements and modifications, these improvements and modifications also should look protection scope of the present invention.

Claims (8)

1. the spatial data interpolation based on surface theory and curved surface fitting method, it is characterised in that comprise the following steps:

S1, sets initial value matrix F 0 and the sampled data for the treatment of fitting surface；Set up constraint function；

S2, will treat that fitting surface place original area spatial spreading turns to grid points form, and wherein, grid points includes two classes: border grid points and internal grids point；Each border grid points composition borderline region, each internal grids point composition interior zone；

Then according to S1 obtain initial value matrix F 0 and sampled data, calculate the first kind fundamental quantity coefficient E of possessive case site, F, G, Equations of The Second Kind fundamental quantity coefficient L, M, N andI=1,2；J=1,2；K=1,2；Wherein, the area of described first kind fundamental quantity coefficient angle and terrain surface territory for representing both direction in the arc length of terrain surface upper curve, terrain surface；Described Equations of The Second Kind fundamental quantity coefficient is for portraying the bendability in terrain surface space；

S3, utilizes Algebra modeling instrument, sets up object function according to the three of surface theory fundamental equations, and described object function, described constraint function and described initial value matrix F 0 are passed to Optimization Solution device；Wherein, described object function is based upon described interior zone；

S4, described object function is optimized and solves by described Optimization Solution device, obtains the Digitized surface of analog form；

S5, described Optimization Solution device judges whether obtained described Digitized surface meets required precision, if it is satisfied, then directly export described Digitized surface, and process ends；If be unsatisfactory for, then iteration S2-S5, until drawing the Digitized surface meeting required precision, and export the described Digitized surface meeting required precision.

2. the spatial data interpolation based on surface theory according to claim 1 and curved surface fitting method, it is characterised in that in S1, described initial value matrix F 0 includes following several situation:

(1) described initial value matrix F 0 is set as null matrix；

(2) meansigma methods that the element in described initial value matrix F 0 is all sampled data is set.

3. the spatial data interpolation based on surface theory according to claim 1 and curved surface fitting method, it is characterised in that in S1, described constraint function is the equality constraints functions set up according to sampled data, it may be assumed that the constraint function of SIST is f_{M, n}The equality constraint of=height form, wherein m, n represent at grid points (x_m, y_n) there is sampled point at place, the property value of this sampled point is height.

4. the spatial data interpolation based on surface theory according to claim 1 and curved surface fitting method, it is characterised in that S2 specifically includes following steps:

S21, by treating that fitting surface place original area spatial spreading turns to the grid points rectangular array of r row c row, includes r*c grid points altogether, and wherein, line number is followed successively by: 0,1,2 ... r-1；0,1,2...c-1 column number is followed successively by:；

S22, if curved surface expression formula is that (x, y), for any one grid points (x for s=f_i, y_j), wherein, i ∈ (0,1,2...r-1), j ∈ (0,1,2...c-1), calculate the E of this grid points, F, G-value by equation (1)；L, M, N value of this grid points is calculated by equation (2)；This grid points is calculated by equation (4)-(9)I=1,2；J=1,2；K=1,2；

$\left\{\begin{array}{c}E=1+{f_x}^2\\ G=1+{f_y}^2\\ F=f_x\·f_y\end{array}\right.,---\left(1\right)$

Wherein, f_xRepresent f (x, y) first-order partial derivative in the x direction；F_yRepresent f (x, y) first-order partial derivative in y-direction；

$\left\{\begin{array}{c}L=\frac{f_{\mathrm{xx}}}{\sqrt{1+{f_x}^2+{f_y}^2}}\\ N=\frac{f_{\mathrm{yy}}}{\sqrt{1+{f_x}^2+{f_y}^2}}\\ M=\frac{f_{\mathrm{xy}}}{\sqrt{1+{f_x}^2+{f_y}^2}}\end{array}\right.,---\left(2\right)$

Wherein, f_xx, f_yy, f_xyRepresent respectively f (x, y) to the second dervative of x, to the second dervative of y and to each first derivative of x, y；

For any one grid points (x_i, y_i), Gauss equation group is expressed as formula (3):

$\left\{\begin{array}{c}{f}_{\mathrm{xx}}={\mathrm{\Γ}}_{11}^1\·f_x+{\mathrm{\Γ}}_{11}^2\·f_y+L\·{(E\·G-F^2)}^{-\frac{1}{2}}\\ f_{\mathrm{yy}}={\mathrm{\Γ}}_{22}^1\·f_x+{\mathrm{\Γ}}_{22}^2\·f_y+L\·{(E\·G-F^2)}^{-\frac{1}{2}}\\ f_{\mathrm{xy}}={\mathrm{\Γ}}_{12}^1\·f_x+{\mathrm{\Γ}}_{12}^2\·f_y+L\·{(E\·G-F^2)}^{-\frac{1}{2}}\end{array}\right.,---\left(3\right)$

Then:

${\mathrm{\Γ}}_{11}^1=\frac{1}{2}(G\·E_x-2F\·F_x+F\·E_y)\·{(E\·G-F^2)}^{-1}---\left(4\right)$

${\mathrm{\Γ}}_{12}^1=\frac{1}{2}(G\·E_y-F\·G_x)\·{(E\·G-F^2)}^{-1}---\left(5\right)$

${\mathrm{\Γ}}_{22}^1=\frac{1}{2}(2G\·F_y-G\·G_x+F\·G_y)\·{(E\·G-F^2)}^{-1}---\left(6\right)$

${\mathrm{\Γ}}_{11}^2=\frac{1}{2}(2E\·F_x-E\·E_y-F\·E_x)\·{(E\·G-F^2)}^{-1}---\left(7\right)$

${\mathrm{\Γ}}_{12}^2=\frac{1}{2}(E\·G_x-F\·E_y)\·{(E\·G-F^2)}^{-1}---\left(8\right)$

${\mathrm{\Γ}}_{22}^2=\frac{1}{2}(E\·G_y-2F\·F_y+F\·G_x)\·{(E\·G-F^2)}^{-1}---\left(9\right)$

Wherein, E_x、E_y、F_x、F_y、G_x、G_yRepresent that x, y are sought partial derivative by E, F, G respectively.

5. the spatial data interpolation based on surface theory according to claim 4 and curved surface fitting method, it is characterised in that in S3, utilize Algebra modeling instrument, set up object function according to the three of surface theory fundamental equations particularly as follows:

S31, by the equal sign left end centered finite difference formal expansion of formula (3), then is multiplied by h at equal sign two ends simultaneously², and keep primitive form constant as constant term formula (3) right-hand member, using formula (3) left end as variable, then formula (3) is re-expressed as:

$\left\{\begin{array}{c}{f}_{i+1,j}-2f_{i,j}+f_{i-1,j}=({\mathrm{\Γ}}_{11}^1\·f_x+{\mathrm{\Γ}}_{11}^2\·f_y+L\·{(E\·G-F^2)}^{-\frac{1}{2}})\·h^2\\ f_{i,j+1}-2f_{i,j}+f_{i,j-1}=({\mathrm{\Γ}}_{22}^1\·f_x+{\mathrm{\Γ}}_{22}^2\·f_y+N\·{(E\·G-F^2)}^{-\frac{1}{2}})\·h^2\\ f_{i+1,j+1}-2f_{i-1,j+1}+f_{i-1,j-1}-=f_{i+1,j-1}({\mathrm{\Γ}}_{12}^1\·f_x+{\mathrm{\Γ}}_{12}^2\·f_y+M\·{(E\·G-F^2)}^{-\frac{1}{2}})\·h^2\end{array}\right.;---\left(10\right)$

Wherein, h is material calculation, i.e. pixel resolution；F_{I, j}(x, y) at grid points (x to represent f_i, y_j) value at place；

S32, is abbreviated as RHS1, RHS2 and RHS3 respectively by the right-hand vector of formula (10), right-hand member moves to left end simultaneously, obtains below equation expression formula:

$\left\{\begin{array}{c}\mathrm{ex}1=f_{i+1,j}-2f_{i,j}+f_{i-1,j}-{\mathrm{RHS}1}_{i,j}\\ \mathrm{ex}2=f_{i,j+1}-2f_{i,j}+f_{i,j-1}-{\mathrm{RHS}2}_{i,j}\\ \mathrm{ex}3=f_{i+1,j+1}-2f_{i-1,j+1}+f_{i-1,j-1}-f_{i+1,j-1}-{\mathrm{RHS}3}_{i,j}\end{array}\right.;---\left(11\right)$

S33, in formula (11), the expression formula that minimizes of the quadratic sum of three expression formulas is the object function of SIST method, and wherein, the span of i is 1,2 ..., r-2；The span of j is 1,2 ..., c-2；Being write as the form of mathematical function, the object function of SIST method is expressed as:

$\min \left\{\underset{i=1}{\overset{r-2}{\mathrm{\Σ}}}{\underset{j=1}{\overset{c-2}{\mathrm{\Σ}}}}_{+{(f_{i+1,j+1}-f_{i-1,j+1}+f_{i-1,j-1}-f_{i+1,j-1}-{\mathrm{RHS}3}_{i,j})}^2}^{{(f_{i+1,j}-{2f}_{i,j}+f_{i-1,j}-{\mathrm{RHS}1}_{i,j})}^2+{(f_{i,j+1}-{2f}_{i,j}+f_{i,j-1}-{\mathrm{RHS}2}_{i,j})}^2}\right\}.$

6. the spatial data interpolation based on surface theory according to claim 5 and curved surface fitting method, it is characterised in that also include:

If z=is (f_0,0, f_0,1..., f_{0, c-1}, f_1,0, f_1,1..., f_{I, 0}..., f_{I, c-1}..., f_{R-1, c-1})^T, then equation (10) is expressed as following matrix form:

$\left\{\begin{array}{c}\mathrm{Az}=\mathrm{RHS}1\\ \mathrm{Bz}=\mathrm{RHS}2\\ \mathrm{Cz}=\mathrm{RHS}3\end{array}\right.,---\left(12\right)$

A, B and C are coefficient matrix, express in the following ways:

Definition Z_{M, n}Represent a null matrix, the line number of this null matrix and columns respectively m and n, then: coefficient matrices A is expressed as:

Wherein, A_{C, (r-2) * c}Represent the matrix of a line number and columns respectively c and (r-2) * c, it may be assumed that

Coefficient matrix B is expressed as:

Wherein, B_{C, (r-2) * c}Represent the matrix of a line number and columns respectively c and (r-2) * c, be expressed as:

Coefficient matrix C is:

Wherein, C_{C, (r-2) * c}Represent the matrix of a line number and columns respectively c and (r-2) * c, be expressed as

7. the spatial data interpolation based on surface theory according to claim 5 and curved surface fitting method, it is characterised in that in S4, described Optimization Solution device described object function is optimized solve particularly as follows:

Under the constraints of described constraint function, described object function is optimized and solves.

8. the spatial data interpolation based on surface theory according to claim 1 and curved surface fitting method, it is characterized in that, in S5, described Optimization Solution device judges whether obtained described Digitized surface meets required precision particularly as follows: assume twice curve modeling result respectively two the column vector z in front and back_iAnd z_i+1, then when | | z_i-z_i+1||₂/||z_i||₂During≤delta, show that precision meets the conclusion of requirement；Wherein, delta span is: equal to or less than the numeral of 1e-10.