CN106683185B

CN106683185B - High-precision curved surface modeling method based on big data

Info

Publication number: CN106683185B
Application number: CN201710013712.0A
Authority: CN
Inventors: 赵娜; 岳天祥
Original assignee: Institute of Geographic Sciences and Natural Resources of CAS
Current assignee: Institute of Geographic Sciences and Natural Resources of CAS
Priority date: 2017-01-09
Filing date: 2017-01-09
Publication date: 2020-04-03
Anticipated expiration: 2037-01-09
Also published as: CN106683185A

Abstract

The invention relates to a high-precision curved surface modeling method based on big data, which comprises the following steps: establishing geographic coordinate information and a variable sampling value to be measured of each sampling point; discretizing the space of the area to be measured into a grid point form, and establishing a sampling equation; carrying out high-order difference dispersion on a partial differential equation set of the curved surface to obtain a corresponding algebraic system, combining the algebraic system and a sampling equation into an equality constraint least square problem, then converting the equality constraint least square problem into a problem of solving a minimal value of a truncated objective function, and then converting the problem into a problem of solving a symmetrical indeterminate equation set; randomly selecting an iteration initial value; partitioning the coefficient matrix into blocks, and storing the decomposed block matrix; solving the HASM equation set by adopting a block line projection iteration method, and judging whether the solution result is converged; judging whether the solution meets a Gauss Kjeldahl equation set or not; and outputting the high-precision simulation curved surface model. The invention can solve the large-scale problem, has small required storage space, and solves the defects of HASM in solving the large-scale problem under the condition of ensuring the precision.

Description

High-precision curved surface modeling method based on big data

Technical Field

The invention relates to a space curved surface modeling method aiming at large-scale problems, which can be used in the fields of generation of a digital elevation model under a big data background, space distribution simulation of meteorology, biomass and soil and the like, can also be regarded as a method for curved surface grid approximation, and can be used for large-scale curved surface modeling in the aspects of physics, chemistry, medicine and the like.

Background

In order to solve the error problem and the multi-scale problem which puzzle the curved Surface modeling for half a century, a High Accuracy Surface Modeling (HASM) method which takes global approximate data (including remote sensing data and global model coarse resolution simulation data) as a driving field and local High Accuracy data (including monitoring network data and survey sampling data) as an optimization control condition is established by taking a differential geometry principle and an optimization control theory as theoretical bases, and the basic law of earth Surface modeling is refined and formed on the basis of a large number of application researches for more than 20 years.

HASM can finally be transformed into an equality constrained least squares problem (LSE) constrained by ground sampling in order to keep the overall simulation error to a minimum under the condition that the simulation value at the sampling point is guaranteed to be equal to the sampling value. And the method fully utilizes sampling information and is an effective means for ensuring that iteration approaches to the optimal simulation effect. Although the HASM method is greatly improved in simulation accuracy compared with the traditional interpolation method, research shows that the current HASM can only process small-scale problems, and the solution process is to convert the LSE problem into a normal equation system for solution. On one hand, the condition number of the transformed equation set coefficient matrix is very large, so that the sensitivity and the ill-conditioned performance of the problem are caused, the convergence speed of an iterative method for solving the equation set is very slow or is not converged, and therefore an approximate solution of the original solution problem cannot be obtained; on the other hand, the currently adopted preprocessing conjugate gradient method for solving the HASM equation set needs to input all the coefficient matrixes of the equation set into a memory at one time, and along with the increase of the problem solving scale, the storage space of the HASM equation set is obviously increased, so that the current solving strategy for large-scale problems is not used. The current HASM is used for simulating national 1km resolution data, a fragment processing strategy is often adopted, errors at partition boundaries are increased, and therefore the overall precision of a simulation result is reduced.

Disclosure of Invention

Aiming at the defects of the existing HASM model in solving the large-scale problem, the invention provides the improved HASM, which has small required storage space and can solve the large-scale problem under the condition of ensuring the precision.

The technical scheme for solving the technical problems is as follows: a high-precision curved surface modeling method based on big data comprises the following steps:

step 1: establishing geographic coordinate information and a variable sampling value to be measured of each sampling point, wherein the geographic coordinate information comprises longitude information and latitude information of the sampling points;

step 2: the method comprises the steps of discretizing a region to be measured into a grid point form to obtain grid point discrete values, establishing a sampling equation according to geographic coordinate information and a variable sampling value to be measured, wherein the sampling equation is used for judging whether a sampling point is on a grid point or not, and if the sampling point is on the grid point, the value of the grid point is the variable sampling value to be measured; if the sampling point is in the grid, obtaining an approximate sampling value on the grid point by using Taylor expansion on the grid point closest to the sampling point;

and step 3: calculating a first basic quantity E, F, G and a second basic quantity L, M, N of each grid point of the region to be measured according to the discrete values of the grid points, wherein the first basic quantity is used for expressing the length of a curve on a simulated curved surface, the area of the simulated curved surface and the curvature of the curve on the simulated curved surface, the second basic quantity is used for expressing the local bending change degree of the simulated curved surface, performing high-order differential discretization on partial differential equations of the curved surface expressed by the first basic quantity and the second basic quantity to obtain an algebraic system corresponding to a discrete equation set, and combining the algebraic system and the sampling equation into an equation constrained least square problem;

and 4, step 4: converting an equality constraint least square problem into a solution of a cut-off objective function minimum value problem;

and 5: converting the solution of the minimal value problem of the truncated objective function into a solution of a symmetrical and uncertain equation set, wherein the symmetrical and uncertain equation set is a high-precision curved surface modeling HASM equation set;

in particular, the method comprises the following steps of,

step 6: randomly selecting an iteration initial value of the HASM equation set;

and 7: performing row blocking on a coefficient matrix in the HASM equation set, and storing a block matrix after coefficient matrix decomposition;

and 8: solving the HASM equation set by adopting a block line projection iteration method for the initial iteration value, and judging whether the solution result is converged;

and step 9: when the solution result is not converged, solving the HASM equation set by adopting the block-line projection iterative method again for the solution result, judging whether the solution result is converged, if so, executing the step 10, otherwise, executing the step 9 again;

step 10: when the solution of the HASM equation set is converged, further judging whether the solution of the HASM equation set meets the Gauss-Kyori equation set, if not, executing the step 6; and if so, outputting a high-precision simulation curved surface model related to the variable to be measured according to the solution of the HASM equation set.

The invention aims to solve the big data problem in the field of space curved surface modeling, and has the following advantages:

1. the existing model is perfected and developed, and the improved model has lower storage requirement and shorter calculation time;

2. a truncation system is introduced, an optimization solving method for solving the constraint least square problem is provided, and the result is more stable;

3. a specific blocking mode of a large equation system coefficient matrix is provided, the linear independence between each row of the matrix is ensured, and the information redundancy is avoided, so that an acceleration algorithm of a row projection algorithm, namely a block row projection iterative method is provided, and the convergence speed of the method is improved;

4. the improved method can solve a large-scale problem, the calculation time and the calculation grid number are in a linear relation, and the required storage space is far smaller than the number of non-zero elements in a coefficient matrix of the HASM equation set; under the condition of ensuring the precision, the defects of the HASM in solving a large-scale problem are overcome.

On the basis of the technical scheme, the invention can be further improved as follows.

Further, the partial differential equation system of the curved surface is as follows:

wherein E is 1+ f_x ²，F＝f_x·f_y，

x is the abscissa of the spatially discretized grid point, y is the ordinate of the spatially discretized grid point, f_xIs the first partial derivative of the function f on x, f_xxIs the second partial derivative of the function f to x, f_yAs the first partial derivative of the function f with respect to y, f_yySecond partial derivative E of function f to y_x、F_x、G_xE, F, G first partial derivatives of x, E_y、F_y、G_yRespectively E, F, G for the first partial derivative of y,

a second class of cristokes symbols.

The method has the advantages that the HASM is built on the basis of a complete differential geometry theory, the stability of the HASM is guaranteed, and the simulation accuracy of the HASM is improved.

Further, the step 3 specifically includes:

step 3.1: calculating the values of the first type basic quantity and the second type basic quantity of each grid point of the area to be measured according to the finite difference approximation of the first type basic quantity E, F, G and the second type basic quantity L, M, N;

step 3.2: combining the finite difference approximation of the first type basic quantity, the finite difference approximation of the second type basic quantity and the finite difference of the second type Crisefer symbol to obtain a discrete equation set;

step 3.3: converting the discrete equation set into an algebraic system in a matrix form;

step 3.4: and combining the algebraic system and the sampling equation into an equality constraint least square problem.

The HASM is finally converted into an equality constraint least square problem constrained by ground sampling, and the purpose is to keep the minimum integral simulation error under the condition of ensuring that the simulation value at the sampling point is equal to the sampling value. And the method fully utilizes sampling information and is an effective means for ensuring that iteration approaches to the optimal simulation effect.

Further, the specific method for partitioning the coefficient matrix in the HASM equation set in step 7 is as follows: setting each sub-block A of the coefficient matrix A_iAt most u rows, and A_i∈R^u×nU < n are all full ranks, i.e. the rows of each block are independent, assuming that one sub-block A is being generated_iThe sub-block already contains j rows (j < u) in matrix A, and knowing A_i ^TDecomposition of QR into A_i ^T＝Q_jR_jBlock A_i ^TEstimation of line correlation

Wherein

r_hhIs located at R_jThe element at the (h, h) position of (a)₁For a row in matrix A, order

According to A_i ^TBy using QR decomposition update method

Is decomposed and then calculated

If the estimate is still less than k, then a₁Adding block A_iWherein κ is a preset positive number.

Further, the QR decomposition adopts a Francis QR decomposition strategy.

The beneficial effect of adopting the further scheme is that for related QR decomposition, the calculation cost of the Francis QR decomposition strategy is O (u) by adopting the Francis QR decomposition strategy³) Reduced to O (u)²) U < n, thereby improving computational efficiency. The convergence of the line projection iterative method can be accelerated by partitioning the matrix, and the calculation time for solving the HASM equation set is reduced.

Further, the method for storing the block matrix after the coefficient matrix decomposition in step 7 is a row compression storage method CSR.

The CSR format is the most common and flexible compression format, the non-zero elements of the matrix are compressed and stored according to rows, the original positions of the non-zero elements are recorded by a special array, the compression efficiency is high, the compression process is convenient to understand, and therefore the memory requirement is obviously reduced.

Further, the steps 8 and 9 specifically include: assuming the current iteration point is, judge x^kIf it is not, find out x^kSelecting x from the block with the farthest orthogonal distance^kThe projection on the block serves as the next iteration point.

The beneficial effect of adopting the further scheme is that the convergence speed is accelerated.

Drawings

FIG. 1 is a flowchart of steps of a big data-based high-precision curved surface modeling method according to an embodiment of the present invention;

FIG. 2 is a diagram of a standard mathematical surface and a distribution of sampling points;

fig. 3 is an error diagram of the simulated curved surface and the real curved surface of the HASM, FRK, IDW and FRS.

Detailed Description

The principles and features of this invention are described below in conjunction with the following drawings, which are set forth by way of illustration only and are not intended to limit the scope of the invention.

As shown in fig. 1, a high-precision curved surface modeling method based on big data provided by an embodiment of the present invention includes the following steps:

the step 3 specifically includes:

Specifically, according to the basic theorem of the surface theory: assuming that the first and second basis quantities E, F, G, L, M and N of the curved surface satisfy symmetry, E, F, G is positive, E, F, G, L, M and N satisfy the partial differential equation set of the curved surface, i.e., Gauss equation set, then the full differential equation set is set at f (x, y) ═ f (x, y)₀,y₀)(x＝x₀,y＝y₀) Under initial conditions, there is a unique solution z ═ f (x, y).

The expression of the Gauss system of equations is,

a second class of cristokes symbols.

If { (x)_i,y_j) Is an orthogonal subdivision of the computational domain omega, [0, Lx }]×[0,Ly]For the dimensionless normalized calculation domain, max { Lx, Ly } ═ 1, Lx is the normalized integer value of the horizontal length of the study region. Ly is the normalized integer value of the length of the study region in the vertical direction.

To calculate the step size, { (x)_i,y_j) I is more than or equal to |0 and less than or equal to I +1, and J is more than or equal to 0 and less than or equal to J +1 is a grid of a standardized calculation domain, the finite difference approximation of the first type basic quantity is,

wherein E is_i,j，F_i,j，G_i,jE, F, G at grid points (i, j), respectively, f_i,jIs the value of the simulated surface at grid point (i, j);

a second class of finite difference approximations of the basis quantities is,

wherein L is_i,j，M_i,j，N_i,jL, M, N at grid points (i, j), respectively;

the finite difference of the second class of cristokes-fler symbols is expressed as,

wherein the content of the first and second substances,

the values at lattice points (i, j) for the second class of Crisefer symbols;

the finite difference form of the system of gaussian equations is,

the matrix form of (1.1.2) can be written as:

wherein the content of the first and second substances,

wherein the content of the first and second substances,

discrete values in finite difference format are used for the right term of the above equation (1.1.2).

In conjunction with the efficient constrained control of the sampled information, the constrained least squares problem (1.1.3) described above can be expressed as an equality-constrained least squares problem solved by HASM,

sx is a sampling equation which is satisfied by sampling points, wherein S and k are a sampling matrix and a sampling vector respectively; if it is not

Is that z is f (x, y) at the p-th sampling point (x)_i,y_j) A value of (1), then s_p,(i-1)×J+j＝1，

HASM finally translates into an equality constrained least squares problem (LSE) constrained by ground sampling

specifically, it is found through research that, by introducing a parameter λ, the LSE problem can be first converted into the following minimal value problem of the truncated objective function by Lagrange multiplier method:

wherein

specifically, a pair

Taking the derivative of 0 to obtain

Bonding of

And Sx ═ k, the following symmetric undefined system can be obtained:

wherein, I is an identity matrix in which all diagonal elements are 1 and other elements are 0.

Is a matrix

The transposed matrix of (2). S^TIs the transposed matrix of S. And lambda is a weight matrix formed by the weights of the sampling points, and the value of the lambda in the HASM is between 1 and 2.

Step 6: randomly selecting an iteration initial value of the HASM equation set;

specifically, for the block division mode of the coefficient matrix, the research provides a block division algorithm of A, each block has at most u rows by applying the algorithm, and each block A has_i∈R^u×nU < n is full rank, i.e. the rows of each block are independent. Suppose that a block A is being generated_iThe block already containsJ rows in matrix A (j < u), and knowing A_i ^TDecomposition of QR into A_i ^T＝Q_jR_jLet r_hhIs located at R_jThe (h, h) position of (c). Block A_i ^TEstimation of line correlation

Wherein

For a row a in matrix A₁Whether or not to add the block A_iHas a criterion of₁Is not assigned to any other block, and a₁After addition of A_iLine correlation estimation β_iAnd k is not more than k, and is a preset positive number. To determine a new column a₁Whether or not to add block A_iLet us order

And according to the existing A_i ^TBy QR decomposition to obtain

For which the QR decomposition update method is adopted. Then calculate

If the estimate is still less than k, then a₁Adding block A_i。

The block decomposition algorithm is described below, and is based on a C + + implementation.

Let Γ equal 1,2, …, m, Γ_sA in ═ l: A₁Column has been assigned to a block }

Initialization: let i equal to 1, Γ_s＝φ，

whileΓ_s≠Γ,do

Let F be_c＝Γ\Γ_s，j＝1,Γ_i＝φ

Selecting l epsilon gamma_cLet A_i＝[a₁],R₁＝[||a₁||₂,0]^T,

Computing

Q₁＝I_n-2vv^T

Let d_max＝||a₁||₂，d_min＝||a₁||₂

Γ_c←Γ_c\{l}，Γ_s←Γ_s∪{l}，Γ_i←Γ_i∪{l}

while (j < u and Γ)_c≠φ)do

Selecting l epsilon gamma_c

Order to

v＝[||a||₂,0]^T∈R^n-j

Computing

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

if t_minNot equal to 0, and

then

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

i←i+1

end

To improve the computational efficiency, for the involved QR decomposition, we adopt the Francis QR decomposition strategy and set the computation cost of the strategy to O (u)³) Reduced to O (u)²),u＜＜n.

The algorithm involves sparse matrix block A_iMultiplication between them, operation of sparse matrix and vector multiplication, vector inner product and vector addition operation. Sparse matrices are typically stored in compressed format. The matrices are typically compressed into different storage formats depending on the application scenario and the computing platform. The selection of the compression format comprehensively considers the sparse characteristic of the matrix and a computing platform. The research adopts a row compression storage mode (CSR) for matrix storage, the CSR format is the most common and flexible compression format, non-zero elements of a matrix are compressed and stored according to rows, a special array is used for recording the original positions of the non-zero elements, the compression efficiency is high, and the compression process is convenient to understand. The CSR format is then ported to a different platform. Many new compression formats are also modified based on the CSR format. When the method stores the n-order matrix A, assuming that the A has l non-zero elements in total, the non-zero elements in the A need to be stored in sequence by using one l-dimensional vector x according to the sequence of the first row and the later row, and one l-dimensional vector x is used^JStoring the non-zero element column numbers in A in sequence according to the same sequence, and introducing an n +1 dimensional integer vector x^(R)，

Indicates where in x the first non-zero element in row i in a is stored,

the coefficient matrix of the indefinite equation set (4) in the research is a symmetric matrix, so that only the non-zero elements of the upper triangular part need to be stored in practice. For the Huge type, a block line compression storage mode is correspondingly adopted.

The matrix and vector multiplication y ═ Av implemented based on the above format can be written as:

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

end

based on the CSR line compression storage mode, the maximum calculation cost of the algorithm can be analyzed and obtained to be O (n), wherein n is the number of calculation grids.

specifically, the solution method for the above-mentioned asymmetric system is generally a direct method and an iterative method. The direct method involving matrix decomposition and the like is not suitable for large-scale problems. The commonly used iteration method is a Krylov subspace iteration method, a comparison pretreatment conjugate gradient method PCG, a QR decomposition method LSQR based on least square, a generalized minimum residual error method GMRES and the like. Since such iterative methods require that the coefficient matrix of the equation set is input into the memory all at once, such methods are no longer used as the scale of the calculation increases. The line projection iterative method can avoid the problems, the line projection iterative method is wide in application and simple in algorithm, and classical line projection iterative methods such as a Cimmiio algorithm and a Kaczmar algorithm are adopted.

Definition H_iAs a set:

Rⁿi is 1, …, p, which is the space formed by the n-dimensional real number vector. Therefore any

I.e. x^*Is p subspaces H_iIs a linear system A^Tx is the solution of b. Set a desiredThe asymmetric system in (1.1.5) is Ax ═ b, and the matrix a is decomposed into block-row form according to the mesh planing of the computation region:

and define P_k(A^iT) Is A^iTProjection operator over a range of values. The general idea is to firstly x_kProjected onto p hyperplanes

Algorithm block Cimmino line projection:

selecting x⁽⁰⁾When k is equal to 0, the value of k is,

the above-mentioned steps are repeated 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 avoid that all matrix elements are not required to be input every time, only one row of the matrix is required to be input every time, and the solving possibility of large-scale problems is guaranteed. But has the disadvantage of slow convergence and is therefore less practical.

Therefore, the research is based on the Kaczmar algorithm, and a new line projection block iterative algorithm is provided. The block which is farthest away from the current iteration point is selected to perform orthogonal projection, and the projection is used as the next iteration point, so that the convergence speed is accelerated. Assume that the current iteration point is x^kFirst, x is calculated^kIn all H_iOrthogonal projection P of_i(x^k) For i | | | P, 1, …, P, all P's are then compared_i(x^k)-x^k||₂Finding the maximum of them, i.e. finding the distance x^kBlock H with farthest orthogonal distance_iSuppose | | | P_j(x^k)-x^k||₂At a maximum, i.e.

Selecting x^kAt H_jProjection P of_j(x^k) As the next iteration point, i.e. x^k+1＝P_j(x^k) Continuing the above process, a sequence { x is obtained^kAnd continuously iterating until convergence. The algorithm is described as follows:

initialization: partitioning the matrix A, giving an initial point x⁰∈R ⁿ0 ≤ epsilon < 1, and let k equal to 0

Step 10: when the solution of the HASM equation set is converged, further judging whether the solution of the HASM equation set meets the Gauss-Kyori equation set, if not, executing the step 6; and if so, outputting a high-precision simulation curved surface model related to the variable to be measured according to the solution of the HASM equation set. .

The following describes the verification of an embodiment of the present invention,

1. numerical test

Firstly, taking a standard mathematical surface as an example, the effectiveness of the algorithm is researched. Comparing the simulation performance difference of the improved HASM method with a Kriging method (FRK) which can adapt to large-scale problem solving, a Spline method (FRS) which can adapt to large-scale problem solving and an IDW (inverse distance weighted method).

The domain is defined as [0,1] × [0,1 ].

As shown in FIG. 1, 1/100 sampling points are randomly selected on a standard mathematical surface, all experiments are carried out on a 64-bit computer, and the halt criterion in HASM is Gauss-Codazii equation set. The calculation area is divided into 1000 × 1000.u is 100, k is 5. the measurement of the calculation error is measured by the average absolute error RMSE, which is calculated as follows:

wherein f is_i，

The ith real value and the ith analog value are obtained, N is the number of sampling points, and the calculation result is shown in Table 1:

TABLE 1 comparison of computational efficiency analysis by different methods

Method of producing a composite material	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

The results show that HASM has the highest calculation accuracy, followed by FRS method. HASM calculation precision is 127 times, 17 times and 4 times of IDW, FRK and FRS respectively. But the computation speed of the HASM is lower than other methods. This is because finite difference dispersion of partial differential equation sets is required in the hash calculation process, right-end terms of the equation sets need to be calculated, the equation sets need to be solved, and the calculation process is more complicated.

Fig. 2 is a graph of the error of a simulated surface from a real surface according to (a) HASM, (b) FRK, (c) G-IDW, and (d) FRS, wherein it can be seen that the IDW simulation error is the largest, and the largest error occurs at the boundary and peak. The maximum deviation of 0.2139.FRK also produced a large error, the maximum deviation of 0.02737.FRS appeared at the boundary of-0.0123. the smallest error is the HASM method, it can be seen that the error distribution is relatively uniform.

Although the computation time of the HASM is more than that of other methods, the computation time is linearly proportional to the computation grid. Table 2 shows the calculated time of the HASM at different calculation scales. The relationship of computation time to computation grid can be expressed as:

t＝1.1099×10^-5n+26.7447,R²＝0.91

TABLE 2 calculation time of HASM at different calculation scales

Calculating the number of grids	1million	4million	9million	16million	25million
						Calculating time(s)	54.06	88.21	107.54	155.19	339.15

2. Practical case

With China as a research area, the method is adopted to calculate the average precipitation distribution in the past 50 years of 1961-2010, the spatial resolution is 1km, and the corresponding calculation grid number is 19,606,916.

As can be seen from the average precipitation distribution results of HASM, FRK, IDW and FRS in China for years, all the methods can better reflect the precipitation distribution trend. The curved surface simulated by the HASM and FRK method is smoother than that simulated by other methods. The IDW and FRS methods exhibit a number of bulls-eye phenomena, and the FRS methods produce large amplitude oscillations at the boundaries.

Randomly select 15% as verification points to verify different methods. This process was repeated 10 times and the average RMSE was calculated for these 10 times. The results are shown in Table 3. The HASM calculation precision is the highest, and the simulation precision is 1.4,1.5 and 1.6 times of those of FRK, IDW and FRS respectively. The simulation error of FRS is the largest. IDW is the least time to compute, followed by the HASM.

TABLE 3 comparison of the different methods in the actual case

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

The numerical test and the practical case show that the improved HASM method ensures the accuracy and simultaneously leads the calculation time to be in linear relation with the calculation grid number, and in the aspect of storage, because a block line projection algorithm is adopted when solving the constraint least square problem, only one block line of the matrix needs to be stored, and a line compression storage mode is adopted to store the block line, so that the scale of the problem solving is ensured.

The above description is only for the purpose of illustrating the preferred embodiments of the present invention and is not to be construed as limiting the invention, and any modifications, equivalents, improvements and the like that fall within the spirit and principle of the present invention are intended to be included therein.

Claims

1. A high-precision curved surface modeling method based on big data comprises the following steps:

and step 3: calculating a first type basic quantity E, F, G and a second type basic quantity L, M, N of each grid point of the region to be measured according to the discrete values of the grid points, wherein the first type basic quantity is used for expressing the length of a curve on a simulated curved surface, the area of the simulated curved surface and the curvature of the curve on the simulated curved surface, the second type basic quantity is used for expressing the local bending change degree of the simulated curved surface, performing high-order differential discretization on a partial differential equation set of the curved surface expressed by the first type basic quantity and the second type basic quantity to obtain an algebraic system corresponding to the discrete equation set, and combining the algebraic system and the sampling equation into a constraint least square problem;

and 4, step 4: converting an equality constraint least square problem into a problem of solving a cut-off objective function minimum value;

step 6: randomly selecting an iteration initial value of the HASM equation set;

the specific method for partitioning the coefficient matrix in the HASM equation set in the step 7 is as follows: setting each sub-block A of the coefficient matrix A_iAt most u rows, and A_i∈R^u×nU < n are all full ranks, i.e. the rows of each block are independent, assuming that one sub-is being generatedBlock A_iThe sub-block already contains j rows (j < u) in matrix A, and knowing A_i ^TDecomposition of QR into A_i ^T＝Q_jR_jBlock A_i ^TEstimation of line correlation

Wherein

According to A_i ^TBy using QR decomposition update method

Is decomposed and then calculated

If the estimate is still less than k, then a₁Adding block A_iWherein, k is a preset positive number;

2. The method of claim 1, wherein the system of partial differential equations for the curved surface is:

wherein the content of the first and second substances,

F＝f_x·f_y，

a second class of cristokes symbols.

3. The method according to claim 2, wherein the step 3 specifically comprises:

4. The method of claim 1, wherein the QR decomposition employs a Francis QR decomposition strategy.

5. The method according to claim 4, wherein the block matrix after coefficient matrix decomposition in step 7 is stored in a row compression storage mode CSR.

6. The method according to claim 5, wherein the steps 8 and 9 specifically comprise: assume that the current iteration point is x^kJudgment of x^kIf it is not, find out x^kSelecting x from the block with the farthest orthogonal distance^kThe projection on the block serves as the next iteration point.