CN102521882A - Method for obtaining seabed terrain data based on discrete elevation and adaptive mixed weighting - Google Patents

Abstract

The invention provides a method for processing seabed terrain data based on discrete elevation and adaptive mixed weighting. The method comprises the steps of: dividing the horizontal plane of seabed terrain into uniform grids; taking the top points of the grids as control points, wherein the height from each control point to the horizontal plane is elevation data of a read-in digital chart, and the heights of these control points are stored in a three-dimensional array; according to elevation points, firstly calculating basic sub-point reference data by using methods of weighting interpolation and weighting extrapolation, then dynamically calculating new sub-points by using a mixed weighting dynamic correction method, finally performing dynamic scheduling according to an adaptive LOD (Level Of Detail) interpolation correlation table model; establishing a curved surface passing through all the points of all the LOD layers to generate the seabed terrain; and finally attaching textures onto the curved surface to obtain vivid seabed terrain.

Description

Obtain the method for sea bed terrain data based on disperse elevation and ADAPTIVE MIXED weighting

Technical field

The present invention relates to a kind ofly obtain the method for sea bed terrain data, belong to the virtual reality technology field based on discrete elevation and ADAPTIVE MIXED weighting.

Background technology

It is important content in the Virtual Battlefield technical research that the large-scale terrain 3 d image data generates, and the method that generates terrain data at present mainly contains:

1. create landform based on fractal algorithm, be characterized in that polygon quantity is controlled fully, do not have topological mistake, but the iterative computation number of times is more, computing cost was considerable when landform was on a grand scale, and was difficult to generate real landform;

2. generate the method for dimensional topography according to the contour map of BMP form, its advantage is to generate real terrain, but the triangle generating algorithm of net point is complicated between each elevation;

3. utilize the electronic chart data to generate submarine topography; Be characterized in that the gained terrain data is general more coarse; Since the sea bed numerical map that precision height and data point are intensive to obtain difficulty bigger; Need generate the further processing of fine granularity sea bed landform to coarseness sea bed landform, but single weighting extrapolation algorithm can cause the landform distortion.

During popular Vega 3-D view develops software both at home and abroad; Except basic module; Developed several expansion modules, comprised modules such as sensor, cloud, ocean, lacked the hydrospace module that the following space environment in ocean surface, ocean is developed as instrument; One of task that will solve in the exploitation of module under water is a sea bed landform generation technique, promptly generates the sea bed landform that rises and falls according to given electronic chart elevation discrete data.Because the limitation of underwater acoustic measurement method, it is big that general sea chart exists the adjacent data spacing distance, and the curved surface that is linked to be in view of the above is unsmooth poor with the image fidelity that generates, and existing other landform method for building up exists complex algorithm, defective that computing cost is big.

Summary of the invention

For fear of the weak point of prior art, the present invention proposes a kind of method that obtains the sea bed terrain data based on disperse elevation and ADAPTIVE MIXED weighting.

A kind ofly obtain the method for sea bed terrain data, it is characterized in that step is following based on discrete elevation and ADAPTIVE MIXED weighting:

Step 1: according to the digital chart data, the sea bed landform zone that desire is generated evenly is divided into grid by the data break size, and grid vertex wherein is reference mark P;

Step 2: grid is pressed the terrain graph data demand, is sub-point with certain grid vertex, divides three processes to ask for the altitude figures of this son point;

Step a: the weighting interpolation algorithm of putting the altitude figures at all reference mark of peripheral ground floor through son generates certain the sub altitude information z that puts in the zone that the x that makes dot grid is capable and the y row surround _Xy:

$z_{\mathrm{Xy}}=\frac{\underset{i=x}{\overset{x+2}{\mathrm{\Σ}}}\underset{j=y}{\overset{y+1}{\mathrm{\Σ}}}{k}_{(i,j)}{z}_{(i,j)}}{\underset{(i,j)=(x,y)}{\overset{(x+2,y+1)}{\mathrm{\Σ}}}{k}_{(i,j)}}+\mathrm{\ϵ},$ Wherein: ε is the random number between-1～1

Utilize linear relationship obtain son point apart from each reference mark apart from r _{(i, j)}, get 1/r _{(i, j)}Be weighting coefficient k _{(i, j)}(i, j ∈ Z ⁺).

Step b: capable certain the interior son point altitude information z ' of zone that surrounds with the y row of x that generates the reference mark grid through son point with the weighting extrapolation algorithm of the altitude figures at adjacent two reference mark in two-layer on longitude or the latitude _Xy:

$z_{\mathrm{xy}}^{\&prime;}=(\frac{z_{i+2,j+2}}{{(l+r_{\mathrm{xy}})}^2}+\frac{z_{i+2,j+1}}{{{\mathrm{<figure-callout\; id="2"\; label="r\; xy"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">r</figure-callout>}}_{\mathrm{xy}}}^2}+\frac{z_v^{\&prime;}}{r_{\mathrm{xy}}})/2(\frac{1}{{(l+r_{\mathrm{xy}})}^2}+\frac{1}{{{\mathrm{<figure-callout\; id="2"\; label="r\; xy"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">r</figure-callout>}}_{\mathrm{xy}}}^2}+\frac{1}{r_{\mathrm{xy}}})$

$+(\frac{z_{i+2,j-1}}{{(2l-r_{\mathrm{xy}})}^2}+\frac{z_{i+2,j}}{{(l-r_{\mathrm{xy}})}^2}+\frac{z_v^{\&prime;\&prime;}}{l-r_{\mathrm{xy}}})/2(\frac{1}{{(2l-r_{\mathrm{xy}})}^2}+\frac{1}{{(l-r_{\mathrm{xy}})}^2}+\frac{1}{l-r_{\mathrm{xy}}})$

Wherein:

z _XyCapable certain interior son point height of zone that surrounds with the y row of x for the reference mark grid;

X, y are respectively row, the row sequence number of reference mark grid;

n is the total line number in reference mark, and m is the total columns in reference mark;

r _XyExpression point is at the subpoint Q of surface level _XyPut the subpoint Q of adjacent right reference mark with son at surface level _{I, j}Between distance;

Said $z_v^{\&prime;}=z_{i+2,j+1}+\frac{r_{\mathrm{Xy}}}{l}(z_{i+2,j+1}-z_{i+2,j+2})$

Said $z_v^{\&prime;\&prime;}=z_{i+2,j}+\frac{l-r_{\mathrm{Xy}}}{l}(z_{i+2,j}-z_{i+2,j-1})$

I, j is natural number;

Step c: certain the son point that obtains the son point through comprehensive step a and two kinds of methods of step b highly carries out the altitude figures z that the mixed weighting dynamic correcting method generates new son point _x, mixed weighting dynamic calibration formula is:

$z_x=\frac{z_{\mathrm{xy}}+z_{\mathrm{xy}}^{\&prime;}}{2}+\frac{|z_{\mathrm{xy}}-z_{\mathrm{xy}}^{\&prime;}|}{2}\mathrm{\&zeta;}$

Wherein: ζ is a random number in-1～1;

Step 3: the new son point that generates with step 2 belongs to same new LOD level with the reference mark, with this LOD level repetitive cycling step 2, obtains the LOD level that next level point forms next level; Cycle criterion is smaller or equal to 17.

After step 3 is accomplished; Count and the LOD hierarchical relationship is set up 17 grades adaptive multistage LOD interpolation association table lattice model based on sighting distance, interpolation; The altitude data that each son that utilizes OpenGL that step 3 is obtained is put is connected with original control point successively; Generate triangle surface, obtain the sea bed terrain surface of forming by each triangle;

Said adaptive multistage LOD interpolation association table lattice model is:

Supermacro is seen sighting distance	The macroscopic view sighting distance	The microcosmic sighting distance	The ultramicroscopic view sighting distance
				10000	1000	60	3
5000	500	20	1
				4000	200	10	0.5
2000	100	5	0.3
				1～4 grade of LOD	5～8 grades of LOD	9～12 grades of LOD	0.1

Wherein: the ultramicroscopic view sighting distance is divided into 5 grades, the special one deck LOD that increases of refined model, and the corresponding interpolation of number of levels is counted,

The conversion step-length is 5000m, 2000 meters, 1000 meters, 500 meters, 300 meters, 100 meters, 40 meters, 10 meters, 5 meters, 2 meters, 1 meter, 0.5 meter, 0.2 meter.

In the described step c of step 2, the x of reference mark grid is capable, and certain the interior son point of zone that surrounds with the y row highly adopts following formula to generate:

$z_x=\frac{z_{\mathrm{Xy}}+z_{\mathrm{Xy}}^{\&prime;}}{2}+\frac{|z_{\mathrm{Xy}}-z_{\mathrm{Xy}}^{\&prime;}|}{2}\mathrm{\&zeta;},$ Wherein: ζ is a random number in-1～1,

$z_{\mathrm{Xy}}=\frac{\underset{i=x}{\overset{x+1}{\mathrm{\&Sigma;}}}\underset{j=y}{\overset{y+2}{\mathrm{\&Sigma;}}}{k}_{(i,j)}{z}_{(i,j)}}{\underset{(i,j)=(x,y)}{\overset{(x+1,y+2)}{\mathrm{\&Sigma;}}}{k}_{(i,j)}}+\mathrm{\&epsiv;},$ ε is the random number between-1～1

$z_{\mathrm{xy}}^{\&prime;}=(\frac{z_{i+2,j+2}}{{(l+r_{\mathrm{xy}})}^2}+\frac{z_{i+2,j+1}}{{{\mathrm{<figure-callout\; id="2"\; label="r\; xy"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">r</figure-callout>}}_{\mathrm{xy}}}^2}+\frac{z_v^{\&prime;}}{r_{\mathrm{xy}}})/2(\frac{1}{{(l+r_{\mathrm{xy}})}^2}+\frac{1}{{{\mathrm{<figure-callout\; id="2"\; label="r\; xy"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">r</figure-callout>}}_{\mathrm{xy}}}^2}+\frac{1}{r_{\mathrm{xy}}})$

$+(\frac{z_{i+2,j-1}}{{(2l-r_{\mathrm{xy}})}^2}+\frac{z_{i+2,j}}{{(l-r_{\mathrm{xy}})}^2}+\frac{z_v^{\&prime;\&prime;}}{l-r_{\mathrm{xy}}})/2(\frac{1}{{(2l-r_{\mathrm{xy}})}^2}+\frac{1}{{(l-r_{\mathrm{xy}})}^2}+\frac{1}{l-r_{\mathrm{xy}}}),$

$z_v^{\&prime;}=z_{i+2,j+1}+\frac{r_{\mathrm{xy}}}{l}(z_{i+2,j+1}-z_{i+2,j+2}),$ $z_v^{\&prime;\&prime;}=z_{i+2,j}+\frac{l-r_{\mathrm{xy}}}{l}(z_{i+2,j}-z_{i+2,j-1}).$

On the sea bed terrain surface, add texture, generate sea bed landform 3 d image data.

A kind of sea bed terrain data disposal route that the present invention proposes based on discrete elevation and ADAPTIVE MIXED weighting.The surface level of sea bed landform is divided into uniform grid, and grid vertex is the reference mark, and the reference mark is the altitude figures of reading in digital chart apart from the height of surface level; The height at these reference mark is kept in the three-dimensional array; According to spot elevation, at first use the method for inserting in the weighting to calculate basic son point reference data with the weighting extrapolation, use the mixed weighting dynamic correcting method dynamically to calculate the sub-point that makes new advances then; Carry out dynamic dispatching according to self-adaptation LOD interpolation association table lattice model at last; Set up the curved surface that all these LOD layers of a process are had a few, generate the sea bed landform, on curved surface, stick texture at last and just obtain sea bed landform true to nature.

The inventive method; Can obtain each altitude figures that meets sea bed 3-D view rendering request; But and the fast requirement of real time of computing velocity; Thereby solved the very few and discontinuous difficult problem that is difficult to generate real terrain of altitude figures in the electronic chart, can generate the sea bed landform that rises and falls according to given electronic chart elevation discrete data.Having overcome general chart data, to exist the adjacent data spacing distance big, the unsmooth problem with the image fidelity difference that generates of the curved surface that is linked to be in view of the above.The sea bed landform 3-D view that utilizes this method to generate is true to nature, practical, and algorithm is simple, computing cost is little.

Description of drawings

Fig. 1: sea bed terrain data process flow figure

Fig. 2: the interpolation synoptic diagram of one dimension

Fig. 3: the inner interpolation synoptic diagram of two dimensional surface

Fig. 4: one dimension weighting Extrapolating model is analyzed

Fig. 5: the outside interpolation synoptic diagram of calculating of two dimensional surface

Table 1: adaptive multistage LOD interpolation degree of association form

Fig. 6: the uniform grid reference mark altitude figures at 4 * 4 reference mark

Fig. 7: the uniform grid reference mark altitude figures at 16 * 16 newborn reference mark behind the mixed weighting

Embodiment

Combine embodiment, accompanying drawing that the present invention is further described at present:

The surface level of sea bed landform is divided into uniform grid, and grid vertex is the reference mark, and the reference mark is the altitude figures of reading in digital chart apart from the height of surface level; The height at these reference mark is kept in the three-dimensional array; According to spot elevation, use the method for mixing inside and outside weighted correction and self-adaptation LOD interpolation, set up one and pass through the curved surface of being had a few; Generate the sea bed landform, on curved surface, stick texture at last and just obtain sea bed landform true to nature.Sea bed terrain data process flow figure sees Fig. 1.

Natural quality by sea bed can know that its face of land is a fluctuations.It is milder to suppose to rise and fall, and the height change between the then adjacent reference mark is continuous.Promptly the height of son point is determined by the height at reference mark around it, the model of following surface analysis foundation and scheduling sublayer point height.

The ultimate principle of method: in FOV, megarelief has two kinds by the method that the coarseness model generates the fine granularity model: inner interpolation is calculated with outside.Each method can both generate the sub-point of landform; But it does not have self-correcting positron point ability; The existing child point that two kinds of methods are produced carries out mixed weighting and comes the syndrome point, generates the best sub-point that meets the multistage gradient field, and realizes the adaptive hierarchical scheduling according to sighting distance through multilayer LOD.

The physical significance of the dynamic projectional technique of mixed weighting: the variation tendency of the original True Data in two reference mark is reflected in the sub-point on its line; Also be reflected in the sub-point beyond the line; And it is far away more that the son point departs from the reference mark line; It is more little influenced by its variation tendency, and the son point is big more with the direction and the graded angular separation of the most contiguous RP line, receives its variation tendency to influence more little.

(1) inner interpolation method:

General relief block all is equally spaced network, sees from the side, and the interpolation of two three-dimensional coordinate data points is just developed into the geometric relationship of one dimension, its interpolation synoptic diagram such as Fig. 2.Here newly-generated son point P _xLatitude and longitude information can obtain from interpolation number and adjacent 2, its elevation information can draw through weighting.Random number can fluctuate.

(wherein, get k respectively ₁, k ₂For

ε is the random number between-1～1) (1)

Though more than this method correct in theory, distortion is too big, main cause is to consider that not adjacent other point is to its influence.From the two-dimensional plane analysis, its interpolation synoptic diagram such as Fig. 3.Here newly-generated son point Q _xLatitude and longitude information can obtain from interpolation number and adjacent 2, its elevation information can draw through weighting.Z _{(i, j)}The height value of representing the capable j row of i.Then the general calculation formula is following:

$z_{\mathrm{Xy}}=\frac{\underset{i=x}{\overset{x+2}{\mathrm{\&Sigma;}}}\underset{j=y}{\overset{y+1}{\mathrm{\&Sigma;}}}{k}_{(i,j)}{z}_{(i,j)}}{\underset{(i,j)=(x,y)}{\overset{(x+2,y+1)}{\mathrm{\&Sigma;}}}{k}_{(i,j)}}+\mathrm{\&epsiv;},$ Wherein: ε is positioned at the sub-point on the x axle,

Utilize linear relationship obtain son point apart from each reference mark apart from r _{(i, j)}, get 1/r _{(i, j)}Be weighting coefficient k _{(i, j)}(i, j ∈ Z ⁺), ε is the random number between-1～1.

(2) outside projectional technique:

If any 2 P in space ₁And P ₂, at the Q that is projected as of reference plane (being parallel to surface level) ₁And Q ₂, be respectively Z with respect to the height of reference field ₁And Z ₂, as shown in Figure 4.Calculate Q now ₁And Q ₂Arbitrary son point Q outside 2 intervals _xHeight, it is apart from Q ₁, Q ₂Point is respectively r ₁, r ₂At first consider the most simply method of linear extrapolation, then obtain the P among the figure _VPoint, but the bottom contour line that forms like this is a straight line does not meet the characteristics that fluctuate with curved surface in true seabed.Because the height value of son point fluctuates near the height value at its adjacent reference mark; And the reference mark that distance is near more, the influence of antithetical phrase point is big more, so can adopt near the reference mark height of antithetical phrase point to come the height of estimator point as weighted-average method; Distance is far away more, and weight is more little.At first selected distance square reciprocal value be weighting coefficient; After the reference mark carried out weighted mean, suppose to obtain spot elevation and be for

its height calculation formula:

$z_x^{\&prime;}=\frac{\frac{z_1}{{\mathrm{<figure-callout\; id="1"\; label="r"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">r</figure-callout>}}_1^2}+\frac{z_2}{{\mathrm{<figure-callout\; id="2"\; label="r"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">r</figure-callout>}}_2^2}}{\frac{1}{{\mathrm{<figure-callout\; id="1"\; label="r"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">r</figure-callout>}}_1^2}+\frac{1}{{\mathrm{<figure-callout\; id="2"\; label="r"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">r</figure-callout>}}_2^2}}---\left(3\right)$

The following formula arrangement can get:

$z_x^{\&prime;}=\frac{{\mathrm{<figure-callout\; id="1"\; label="z"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">z</figure-callout>}}_1\&CenterDot;{\mathrm{<figure-callout\; id="2"\; label="r"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">r</figure-callout>}}_2^2+{\mathrm{<figure-callout\; id="2"\; label="z"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">z</figure-callout>}}_2\&CenterDot;{\mathrm{<figure-callout\; id="1"\; label="r"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">r</figure-callout>}}_1^2}{{\mathrm{<figure-callout\; id="1"\; label="r"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">r</figure-callout>}}_1^2+{\mathrm{<figure-callout\; id="2"\; label="r"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">r</figure-callout>}}_2^2}=\frac{{\mathrm{<figure-callout\; id="1"\; label="z"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">z</figure-callout>}}_1\&CenterDot;{\mathrm{<figure-callout\; id="2"\; label="r"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">r</figure-callout>}}_2^2+({\mathrm{<figure-callout\; id="2"\; label="z"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">z</figure-callout>}}_2\&CenterDot;{\mathrm{<figure-callout\; id="1"\; label="r"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">r</figure-callout>}}_1^2+{\mathrm{<figure-callout\; id="2"\; label="z"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">z</figure-callout>}}_2\&CenterDot;{\mathrm{<figure-callout\; id="2"\; label="r"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">r</figure-callout>}}_2^2)-{\mathrm{<figure-callout\; id="2"\; label="z"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">z</figure-callout>}}_2\&CenterDot;{\mathrm{<figure-callout\; id="2"\; label="r"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">r</figure-callout>}}_2^2}{{\mathrm{<figure-callout\; id="1"\; label="r"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">r</figure-callout>}}_1^2+{\mathrm{<figure-callout\; id="2"\; label="r"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">r</figure-callout>}}_2^2}=\frac{z_1-z_2}{{\mathrm{<figure-callout\; id="1"\; label="r"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">r</figure-callout>}}_1^2+{\mathrm{<figure-callout\; id="2"\; label="r"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">r</figure-callout>}}_2^2}{\&CenterDot;\mathrm{<figure-callout\; id="2"\; label="r"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">r</figure-callout>}}_2^2+z_2---\left(4\right)$

Can know by P ₁Point is to P ₂The height field graded of point reduces, but can know the following formula analysis,

Point always exceeds P ₂Point is worked as r ₂When increasing gradually, z ' _xWith respect to z ₂The analysis of trend of incremental portion following:

$\underset{r_{}}{\mathrm{<figure-callout\; id="2"\; label="lim\; r"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">lim</figure-callout>}}\frac{z_1-z_2}{{\mathrm{<figure-callout\; id="1"\; label="r"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">r</figure-callout>}}_1^2+{\mathrm{<figure-callout\; id="2"\; label="r"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">r</figure-callout>}}_2^2}\&CenterDot;{\mathrm{<figure-callout\; id="2"\; label="r"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">r</figure-callout>}}_2^2=(z_1-z_2)\&CenterDot;\underset{r_{}}{\mathrm{<figure-callout\; id="2"\; label="lim\; r"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">lim</figure-callout>}}\frac{1}{{\left(\frac{r_1}{r_2}\right)}^2+1}=(z_1-z_2)\&CenterDot;\underset{r_{}}{\mathrm{<figure-callout\; id="2"\; label="lim\; r"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">lim</figure-callout>}}\frac{1}{{\left(\frac{{\mathrm{<figure-callout\; id="2"\; label="r"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">r</figure-callout>}}_2+\stackrel{\&OverBar;}{Q_1{\mathrm{<figure-callout\; id="2"\; label="Q"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">Q</figure-callout>}}_2}}{r_2}\right)}^2+1}$

$=\frac{z_1-{\mathrm{<figure-callout\; id="2"\; label="z"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">z</figure-callout>}}_2}{2}---\left(5\right)$

In the formula

Be Q ₁And Q ₂The distance of point-to-point transmission,

The change curve of point is shown in segment of curve among Fig. 41.

Obviously, prediction equation (3) does not meet the graded trend at reference mark, and r ₂Big more, error is big more.Consider the linear P that obtains that calculates _VThe point height is less than P ₂So the point height is can be with P _VPoint also weighting is incorporated in the top prediction equation, wherein weighting coefficient elect as son point and nearer reference mark apart from r ₂Square inverse, proof, the son point P that obtains of formula thus below " _xLess than P ₂, coincidence control point variation tendency.Expectation point change curve is shown in segment of curve among Fig. 42.Revised prediction equation is:

$z_x^{\&prime;\&prime;}=\frac{\frac{z_1}{{\mathrm{<figure-callout\; id="1"\; label="r"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">r</figure-callout>}}_1^2}+\frac{z_2}{{\mathrm{<figure-callout\; id="2"\; label="r"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">r</figure-callout>}}_2^2}+\frac{{\mathrm{<figure-callout\; id="2"\; label="z\; v\; r"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">z</figure-callout>}}_v}{r_{}^{}}}{\frac{1}{{\mathrm{<figure-callout\; id="1"\; label="r"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">r</figure-callout>}}_1^2}+\frac{1}{{\mathrm{<figure-callout\; id="2"\; label="r"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">r</figure-callout>}}_2^2}+\frac{1}{{\mathrm{<figure-callout\; id="2"\; label="r"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">r</figure-callout>}}_2^2}}---\left(6\right)$

In like manner, work as Z ₁＜Z ₂The time, and along P ₂And P ₁When point is calculated, can release above conclusion similarly left, correction formula is always rational.

By

The change curve and the expectation curve of point can be known, along with r ₂Increase, its error strengthens, this variation tendency still exists in revised (6) formula, therefore can strengthen linear dead-reckoning point P _VWeight, the inverse of only getting apart from first power is a weight coefficient, the weighted mean formula shown in (7) formula, the P that tries to achieve _xThe height z of point _xCan compare z ₂Littler.

$z_x=\frac{\frac{z_1}{{\mathrm{<figure-callout\; id="1"\; label="r"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">r</figure-callout>}}_1^2}+\frac{z_2}{{\mathrm{<figure-callout\; id="2"\; label="r"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">r</figure-callout>}}_2^2}+\frac{{\mathrm{<figure-callout\; id="2"\; label="z\; v\; r"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">z</figure-callout>}}_v}{r_{}}}{\frac{1}{{\mathrm{<figure-callout\; id="1"\; label="r"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">r</figure-callout>}}_1^2}+\frac{1}{{\mathrm{<figure-callout\; id="2"\; label="r"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">r</figure-callout>}}_2^2}+\frac{1}{r_2}}---\left(7\right)$

Need to prove, for generating son point Q _{X, y}Do the weighting extrapolation and need do six calculating, the extrapolation RP is respectively: Q _{I+1, j+1}Q _{I+2, j+1}, Q _{I+1, j}Q _{I+2, j}, Q _{I+3, j}Q _{I+2, j}, Q _{I+3, j+1}Q _{I+2, j+1}, Q _{I+2, j-1}Q _{I+2, j}, Q _{I+2, j+2}Q _{I+2, j+1}But the influence of preceding four weightings extrapolation antithetical phrase point is little; Because also being the graded direction, two RPs extrapolation direction can only explain that it generates the variation tendency of son point on the x direction of principal axis; Can not show its situation of change on the y axle, and newborn sub-point must be along y direction of principal axis two adjacent RP Q _{I+2, j}Q _{I+2, j+1}In.So can only be from Q _{I+2, j-1}To Q _{I+2, j}, from Q _{I+2, j+2}To Q _{I+2, j+1}Doing twice weighting extrapolation calculates.Then average is done in twice extrapolation and can be obtained Q _{X, y}The outside reckoning interpolation synoptic diagram of two dimensional surface weighting is seen Fig. 5.

$z_{\mathrm{xy}}^{\&prime;}=(\frac{z_{i+2,j+2}}{{(l+r_{\mathrm{xy}})}^2}+\frac{z_{i+2,j+1}}{{{\mathrm{<figure-callout\; id="2"\; label="r\; xy"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">r</figure-callout>}}_{\mathrm{xy}}}^2}+\frac{z_v^{\&prime;}}{r_{\mathrm{xy}}})/2(\frac{1}{{(l+r_{\mathrm{xy}})}^2}+\frac{1}{{{\mathrm{<figure-callout\; id="2"\; label="r\; xy"\; filenames="HDA0000116018740000021.png,HDA0000116018740000022.png"\; state="\{\{state\}\}">r</figure-callout>}}_{\mathrm{xy}}}^2}+\frac{1}{r_{\mathrm{xy}}})$ (8)

$+(\frac{z_{i+2,j-1}}{{(2l-r_{\mathrm{xy}})}^2}+\frac{z_{i+2,j}}{{(l-r_{\mathrm{xy}})}^2}+\frac{z_v^{\&prime;\&prime;}}{l-r_{\mathrm{xy}}})/2(\frac{1}{{(2l-r_{\mathrm{xy}})}^2}+\frac{1}{{(l-r_{\mathrm{xy}})}^2}+\frac{1}{l-r_{\mathrm{xy}}})$

Wherein:

z _XyBe capable certain interior son point height of zone that surrounds with the y row of x of reference mark grid,

X; Y is respectively the row of reference mark grid; The row sequence number;

n is the total line number in reference mark; M is the total columns in reference mark

r _XyExpression Q _Xy(the son point is at the subpoint of surface level) and Q _{I+2, j+1}Distance between (the adjacent right reference mark of son point is at the subpoint of surface level);

$z_v^{\&prime;}=z_{i+2,j+1}+\frac{r_{\mathrm{Xy}}}{l}(z_{i+2,j+1}-z_{i+2,j+2})$ With $z_v^{\&prime;\&prime;}=z_{i+2,j}+\frac{l-r_{\mathrm{Xy}}}{l}(z_{i+2,j}-z_{i+2,j-1})$ (i, j are natural numbers)

(3) mixed weighting dynamic correcting method:

Inner interpolation method mainly draws the elevation of son point according to its adjacent RP, its advantage is that gained point has certain confidence level, and shortcoming is not consider the influence of outside RP to it.Outside projectional technique then is to calculate to the graded of adjacent RP according to the outside, and the graded of its son point meets the landform Changing Pattern, but has lost the influence of other adjacent RP antithetical phrase point.Comprehensive two kinds of methods; Can consider this two reasons, to the influence of the graded of RP and adjacent high weight RP antithetical phrase point method, carry out quadratic fit to the outside through weighting; Two son points are proofreaied and correct each other, draw up-to-date son point and satisfy the real terrain data characteristics.It is poor that the child point that two kinds of methods obtain is done, and the value of establishing is a, establishes the scope that [1,1] is new stochastic variable ζ then.Get the average of two sub-point heights, then can obtain the correction position of new son point.

$z_x=\frac{z_{\mathrm{Xy}}+z_{\mathrm{Xy}}^{\&prime;}}{2}+\frac{|z_{\mathrm{Xy}}-z_{\mathrm{Xy}}^{\&prime;}|}{2}\mathrm{\&zeta;}$ (ζ is a random number in-1～1) (9)

(4) adaptive hierarchical scheduling strategy:

In order to realize adaptation function, it is related to combine multilayer LOD scheduling to carry out according to the visual effect needs to sighting distance and terrain mesh interpolation points number, conventional LOD divided method be every layer according to 2 ⁿ-1 point carries out interpolation, and its saltus step width of clothization is bigger.Can use in order step by step according to natural number and insert the point-to-points layer of a son LOD and carry out pre-service this moment, and it is corresponding one by one that LOD number of levels and interpolation are counted, and sets 17 grades of LOD altogether, can use cloud computing to handle above more interpolation of 17 grades.Adaptive multistage LOD interpolation association table lattice model such as table 1.

Sighting distance just is meant that viewpoint arrives the bee-line of relief block; Effect of visualization is very poor when sighting distance surpasses 10000; The relief block of this moment can use a plane to replace; From 10000 to 0 meters sighting distance can be divided into supermacro sight, macroscopic view, microcosmic and four kinds of sighting distances of ultramicroscopic view; Each sighting distance is all set the conversion step-length of multistage LOD layer respectively with the form of successively decreasing, decrement step size is respectively 5000m, 2000 meters, 1000 meters, 500 meters, 300 meters, 100 meters, 40 meters, 10 meters, 5 meters, 2 meters, 1 meter, 0.5 meter, 0.2 meter etc.

Table 1 adaptive multistage LOD interpolation degree of association form

Relief block thus, we have obtained thinner terrain mesh and have divided, and application formula (10) is very easy to programming and calculates.As long as utilize the function of the drafting triangular piece of OpenGL that adjacent net point is linked to each other successively, just obtained the sea bed landform that rises and falls.

For example from digital chart, read the uniform grid reference mark altitude figures at 4 * 4 reference mark.And be stored among the array AN, its sea bed landform is as shown in Figure 6.Re-use the ADAPTIVE MIXED weighting and dynamically calculate terrain generation algorithm; If sighting distance is 2000 meters, call the degree of association model this moment automatically, dynamic dispatching and insertion point become 4 grades of LOD landform; It also is sub-point of each grade insertion in each cell; Insert four sub-points altogether, array becomes 16 * 16 reference mark, is kept at A _SIn.Its sea bed landform is as shown in Figure 7.

$A_N=\left(\begin{array}{cccc}2.911& 2.707& 2.461& 2.240\\ 2.473& 3.333& 3.320& 3.100\\ 2.735& 2.761& 2.961& 2.387\\ 2.864& 2.651& 2.372& 2.456\end{array}\right),$ $A_S=\left(\begin{array}{cccc}{A}_{11}& A_{12}& A_{13}& A_{14}\\ A_{21}& A_{22}& A_{23}& A_{24}\\ A_{31}& A_{32}& A_{33}& A_{34}\\ A_{41}& A_{42}& A_{43}& A_{44}\end{array}\right)$

$A_{11}=\left(\begin{array}{cccc}2.911& 2.599& 2.547& 2.549\\ 2.980& 3.185& 3.030& 3.239\\ 3.190& 3.037& 3.284& 3.024\\ 3.258& 3.019& 2.819& 3.097\end{array}\right),$ $A_{12}=\left(\begin{array}{cccc}2.537& 2.707& 2.973& 2.554\\ 2.833& 2.853& 3.217& 2.807\\ 3.109& 3.263& 3.325& 3.059\\ 3.185& 2.903& 3.169& 3.309\end{array}\right)$

$A_{13}=\left(8\begin{array}{cccc}2.834& 2.698& 2.461& 2.269\\ 3.136& 3.272& 2.841& 2.708\\ 3.180& 3.037& 3.014& 2.624\\ 3.089& 3.325& 2.918& 2.499\end{array}\right),$ $A_{14}=\left(\begin{array}{cccc}2.367& 2.319& 2.198& 2.240\\ 2.354& 2.324& 2.466& 2.761\\ 2.554& 2.795& 2.771& 2.992\\ 2.961& 3.376& 3.111& 3.214\end{array}\right)$

$A_{21}=\left(\begin{array}{cccc}2.903& 2.814& 2.803& 3.083\\ 2.473& 2.527& 2.748& 2.758\\ 3.008& 2.899& 2.926& 3.188\\ 3.199& 2.979& 2.844& 2.934\end{array}\right),$ $A_{22}=\left(\begin{array}{cccc}2.800& 3.025& 3.374& 3.100\\ 2.975& 3.333& 3.510& 3.258\\ 3.181& 3.346& 3.172& 3.103\\ 3.306& 3.293& 3.418& 3.353\end{array}\right)$

$A_{23}=\left(\begin{array}{cccc}2.853& 2.686& 2.672& 3.061\\ 3.234& 3.195& 3.320& 3.150\\ 3.102& 3.121& 3.063& 3.326\\ 2.999& 2.576& 2.516& 2.830\end{array}\right),$ $A_{24}=\left(\begin{array}{cccc}3.241& 3.367& 3.174& 3.231\\ 3.081& 3.096& 3.028& 3.100\\ 3.105& 2.740& 2.578& 2.821\\ 2.833& 2.942& 2.642& 2.504\end{array}\right)$

$A_{31}=\left(\begin{array}{cccc}2.907& 2.644& 2.873& 3.219\\ 2.560& 2.663& 2.827& 2.871\\ 2.735& 2.717& 2.959& 3.296\\ 3.114& 3.060& 3.203& 2.983\end{array}\right),$ $A_{32}=\left(\begin{array}{cccc}3.442& 2.959& 3.128& 3.483\\ 3.169& 2.962& 2.810& 3.146\\ 3.105& 2.761& 2.683& 2.653\\ 3.264& 3.198& 3.209& 2.873\end{array}\right)$

$A_{33}=\left(\begin{array}{cccc}3.344& 2.933& 2.666& 2.396\\ 3.388& 2.967& 2.880& 2.602\\ 2.860& 3.044& 2.961& 2.569\\ 2.574& 2.757& 2.815& 2.737\end{array}\right),$ $A_{34}=\left(\begin{array}{cccc}2.536& 2.743& 2.748& 2.403\\ 2.391& 2.754& 2.625& 2.426\\ 2.435& 2.236& 2.572& 2.387\\ 2.449& 2.286& 2.329& 2.395\end{array}\right)$

$A_{41}=\left(\begin{array}{cccc}3.184& 3.805& 2.967& 2.738\\ 2.995& 2.776& 2.944& 2.934\\ 2.976& 2.682& 2.607& 2.614\\ 2.864& 2.752& 2.418& 2.326\end{array}\right),$ $A_{42}=\left(\begin{array}{cccc}2.796& 2.759& 3.098& 3.116\\ 2.579& 2.723& 3.060& 3.073\\ 2.487& 2.747& 3.806& 3.162\\ 2.651& 2.737& 2.895& 2.856\end{array}\right)$

$A_{43}=\left(\begin{array}{cccc}2.684& 2.400& 2.561& 2.839\\ 2.954& 2.608& 2.311& 2.712\\ 3.056& 2.720& 2.292& 2.477\\ 2.758& 2.770& 2.372& 2.227\end{array}\right),$ $A_{44}=\left(\begin{array}{cccc}2.647& 2.394& 2.240& 2.125\\ 2.564& 2.721& 2.618& 2.583\\ 2.325& 2.502& 2.484& 2.513\\ 2.158& 2.444& 2.406& 2.456\end{array}\right)$

Visible by present embodiment: the inventive method can be obtained each altitude figures that meets sea bed 3-D view rendering request; But and the fast requirement of real time of computing velocity; Thereby solved the very few and discontinuous difficult problem that is difficult to generate real terrain of altitude figures in the electronic chart, can generate the sea bed landform that rises and falls according to given electronic chart elevation discrete data.Having overcome general chart data, to exist the adjacent data spacing distance big, the unsmooth problem with the image fidelity difference that generates of the curved surface that is linked to be in view of the above.The sea bed landform 3-D view that utilizes this method to generate is true to nature, practical, and algorithm is simple, computing cost is little.

Claims (4)

1. one kind obtains the method for sea bed terrain data based on discrete elevation and ADAPTIVE MIXED weighting, it is characterized in that step is following:

Step 1: according to the digital chart data, the sea bed landform zone that desire is generated evenly is divided into grid by the data break size, and grid vertex wherein is reference mark P;

Step 2: grid is pressed the terrain graph data demand, is sub-point with certain grid vertex, divides three processes to ask for the altitude figures of this son point;

Step a: the weighting interpolation algorithm of putting the altitude figures at all reference mark of peripheral ground floor through son generates certain the sub altitude information z that puts in the zone that the x that makes dot grid is capable and the y row surround _Xy:

$z_{\mathrm{Xy}}=\frac{\underset{i=x}{\overset{x+2}{\mathrm{\&Sigma;}}}\underset{j=y}{\overset{y+1}{\mathrm{\&Sigma;}}}{k}_{(i,j)}{z}_{(i,j)}}{\underset{(i,j)=(x,y)}{\overset{(x+2,y+1)}{\mathrm{\&Sigma;}}}{k}_{(i,j)}}+\mathrm{\&epsiv;},$ Wherein: ε is the random number between-1～1

Utilize linear relationship obtain son point apart from each reference mark apart from r _{(i, j)}, get 1/r _{(i, j)}Be weighting coefficient k _{(i, j)}(i, j ∈ Z ⁺).

Step b: capable certain the interior son point altitude information z ' of zone that surrounds with the y row of x that generates the reference mark grid through son point with the weighting extrapolation algorithm of the altitude figures at adjacent two reference mark in two-layer on longitude or the latitude _Xy:

$z_{\mathrm{xy}}^{\&prime;}=(\frac{z_{i+2,j+2}}{{(l+r_{\mathrm{xy}})}^2}+\frac{z_{i+2,j+1}}{{r_{\mathrm{xy}}}^2}+\frac{z_v^{\&prime;}}{r_{\mathrm{xy}}})/2(\frac{1}{{(l+r_{\mathrm{xy}})}^2}+\frac{1}{{r_{\mathrm{xy}}}^2}+\frac{1}{r_{\mathrm{xy}}})$

$+(\frac{z_{i+2,j-1}}{{(2l-r_{\mathrm{xy}})}^2}+\frac{z_{i+2,j}}{{(l-r_{\mathrm{xy}})}^2}+\frac{z_v^{\&prime;\&prime;}}{l-r_{\mathrm{xy}}})/2(\frac{1}{{(2l-r_{\mathrm{xy}})}^2}+\frac{1}{{(l-r_{\mathrm{xy}})}^2}+\frac{1}{l-r_{\mathrm{xy}}})$

Wherein:

z _XyCapable certain interior son point height of zone that surrounds with the y row of x for the reference mark grid;

X, y are respectively row, the row sequence number of reference mark grid;

n is the total line number in reference mark, and m is the total columns in reference mark;

r _XyExpression point is at the subpoint Q of surface level _XyPut the subpoint Q of adjacent right reference mark with son at surface level _{I, j}Between distance;

Said $z_v^{\&prime;}=z_{i+2,j+1}+\frac{r_{\mathrm{Xy}}}{l}(z_{i+2,j+1}-z_{i+2,j+2})$

Said $z_v^{\&prime;\&prime;}=z_{i+2,j}+\frac{l-r_{\mathrm{Xy}}}{l}(z_{i+2,j}-z_{i+2,j-1})$

I, j is natural number;

Step c: certain the son point that obtains the son point through comprehensive step a and two kinds of methods of step b highly carries out the altitude figures z that the mixed weighting dynamic correcting method generates new son point _x, mixed weighting dynamic calibration formula is:

$z_x=\frac{z_{\mathrm{xy}}+z_{\mathrm{xy}}^{\&prime;}}{2}+\frac{|z_{\mathrm{xy}}-z_{\mathrm{xy}}^{\&prime;}|}{2}\mathrm{\&zeta;}$

Wherein: ζ is a random number in-1～1;

Step 3: the new son point that generates with step 2 belongs to same new LOD level with the reference mark, with this LOD level repetitive cycling step 2, obtains the LOD level that next level point forms next level; Cycle criterion is smaller or equal to 17.

2. based on the described sea bed terrain data processing method of claim 1 based on discrete elevation and ADAPTIVE MIXED weighting; It is characterized in that: after step 3 is accomplished; Count and the LOD hierarchical relationship is set up 17 grades adaptive multistage LOD interpolation association table lattice model based on sighting distance, interpolation; The altitude data that each son that utilizes OpenGL that step 3 is obtained is put is connected with original control point successively, generates triangle surface, obtains the sea bed terrain surface of being made up of each triangle;

Said adaptive multistage LOD interpolation association table lattice model is:

Supermacro is seen sighting distance The macroscopic view sighting distance The microcosmic sighting distance The ultramicroscopic view sighting distance 10000 1000 60 3 5000 500 20 1 4000 200 10 0.5 2000 100 5 0.3 1～4 grade of LOD 5～8 grades of LOD 9～12 grades of LOD 0.1

Wherein: the ultramicroscopic view sighting distance is divided into 5 grades, the special one deck LOD that increases of refined model, and the corresponding interpolation of number of levels is counted,

The conversion step-length is 5000m, 2000 meters, 1000 meters, 500 meters, 300 meters, 100 meters, 40 meters, 10 meters, 5 meters, 2 meters, 1 meter, 0.5 meter, 0.2 meter.

3. the sea bed terrain data disposal route based on discrete elevation and ADAPTIVE MIXED weighting according to claim 1; It is characterized in that: in the described step c of step 2, the x of reference mark grid is capable, and certain the interior son point of zone that surrounds with the y row highly adopts following formula to generate:

$z_x=\frac{z_{\mathrm{Xy}}+z_{\mathrm{Xy}}^{\&prime;}}{2}+\frac{|z_{\mathrm{Xy}}-z_{\mathrm{Xy}}^{\&prime;}|}{2}\mathrm{\&zeta;},$ Wherein: ζ is a random number in-1～1,

$z_{\mathrm{Xy}}=\frac{\underset{i=x}{\overset{x+1}{\mathrm{\&Sigma;}}}\underset{j=y}{\overset{y+2}{\mathrm{\&Sigma;}}}{k}_{(i,j)}{z}_{(i,j)}}{\underset{(i,j)=(x,y)}{\overset{(x+1,y+2)}{\mathrm{\&Sigma;}}}{k}_{(i,j)}}+\mathrm{\&epsiv;},$ ε is the random number between-1～1

$z_{\mathrm{xy}}^{\&prime;}=(\frac{z_{i+2,j+2}}{{(l+r_{\mathrm{xy}})}^2}+\frac{z_{i+2,j+1}}{{r_{\mathrm{xy}}}^2}+\frac{z_v^{\&prime;}}{r_{\mathrm{xy}}})/2(\frac{1}{{(l+r_{\mathrm{xy}})}^2}+\frac{1}{{r_{\mathrm{xy}}}^2}+\frac{1}{r_{\mathrm{xy}}})$

$+(\frac{z_{i+2,j-1}}{{(2l-r_{\mathrm{xy}})}^2}+\frac{z_{i+2,j}}{{(l-r_{\mathrm{xy}})}^2}+\frac{z_v^{\&prime;\&prime;}}{l-r_{\mathrm{xy}}})/2(\frac{1}{{(2l-r_{\mathrm{xy}})}^2}+\frac{1}{{(l-r_{\mathrm{xy}})}^2}+\frac{1}{l-r_{\mathrm{xy}}}),$

$z_v^{\&prime;}=z_{i+2,j+1}+\frac{r_{\mathrm{xy}}}{l}(z_{i+2,j+1}-z_{i+2,j+2}),$ $z_v^{\&prime;\&prime;}=z_{i+2,j}+\frac{l-r_{\mathrm{xy}}}{l}(z_{i+2,j}-z_{i+2,j-1}).$

4. the sea bed terrain data disposal route based on discrete elevation and ADAPTIVE MIXED weighting according to claim 2 is characterized in that: on the sea bed terrain surface, add texture, generate sea bed landform 3 d image data.

Images