CN102324106A - SFS (Shape From Shading) three-dimensional reconstruction sparse-DEM (Digital Elevation Model) encrypting method considering surface spectral information - Google Patents

Abstract

The invention relates to an SFS (Shape From Shading) three-dimensional reconstruction sparse-DEM (Digital Elevation Model) encrypting method considering surface spectral information. The invention creatively provides a coarse-resolution DEM-grid interpolation-encrypting method by combining multi-spectral remote-sensing images with SFS three-dimensional reconstruction; ground-object spectral information is utilized to estimate the reflectivity of different ground-object types, and the influences of surface-vegetation coverage types on the precision of SFS reconstruction are removed; in the invention, multi-spectral images are utilized to replace full-color images, and the precision of three-dimensional reconstruction is improved by applying the assistance of the ground-object spectral information; the influences of the surface-vegetation coverage types and ground-object shadows on establishing a surface three-dimensional model and image spectral-brightness transforming relation can be effectively inhibited; and relative to a commonly-used DEM interpolation algorithm, the precision is markedly improved, and examples validate that the method can be used for improving the DEM spatial resolution to be two times of original data and finishing large-range rapid DEM establishment based on sparse control points.

Description

A kind of SFS three-dimensional reconstruction of taking face of land spectral information into account is encrypted sparse DEM method

Technical field

The present invention relates to a kind of SFS three-dimensional reconstruction and encrypt sparse DEM method, especially relate to a kind of SFS three-dimensional reconstruction of taking face of land spectral information into account and encrypt sparse DEM method.

Background technology

Digital elevation model (DEM) has widespread use at aspects such as mapping remote sensing, military affairs, and directly receives the influence of DEM precision and resolution as the landform basic data.Geologic hazard emergency response rescue also needs high-resolution face of land three-dimensional visualization model in addition, can assist the rescue work of burst accident and disaster accurately to carry out, and therefore setting up high-resolution DEM has very important meaning.

It is according to same sensor areal to be obtained stereogram that traditional DEM sets up, and produces based on photogrammetric principle.Yet in the time also can't obtaining high-resolution stereogram fast and effectively in certain areas, high-resolution DEM sets up the very difficulty that seems.On the other hand, classical interpolation method is the strategy commonly used that obtains high resolving power DEM, through inserting the continuity of face in keeping, on the basis of original coarse resolution data, obtains more high-resolution DEM graticule mesh.At present more effective dem data interpolating method like methods such as golden interpolation in the gram all based on known altitude figures; Can't increase new detailed information; In case ground control point or sampled point have been lost certain key character elevation; Interpolation method just can't recover this regional topographic relief conversion, and therefore traditional interpolation method has very big limitation to encrypting DEM.

It is a kind of classical way of three-dimensional information extraction that light and shade is recovered shape (SFS); Also be that 3D shape is recovered the key of problem technology in the computer vision; Its objective is parameter values such as the relative height utilizing body surface optical density light and shade in the single image to change to recover its surperficial each point or surperficial method direction vector, lay the foundation for further object being carried out three-dimensionalreconstruction.The foreign scholar launches the theoretical research work of SFS successively since the seventies; Horn has proposed SFS thought in 1979; Follow-up other scholars constantly propose all kinds of solutions of SFS problem; Zhang (1999) summarizes to classic algorithm, based on this Jean-Denis Durou (2008) SFS is resolved system and has carried out further perfect.

Thinking according to the SFS algorithm; The remote sensing image spectrum change that this patent utilizes topography and geomorphology to cause is used the principle that light and shade is recovered shape, in conjunction with original coarse resolution dem data and sparse reference mark altitude figures; Obtain encrypting the three-dimensional elevation of DEM according to the mathematical model inverting; Can satisfy stereogram and can't obtain the requirement that the area makes up DEM, take the inner face of land of DEM graticule mesh real change trend again into account, fine-resolution meshes DEM is significant for setting up fast and effectively.Utilize the spectral information of remote sensing image atural object simultaneously, repair the image shadow region, recover the true brightness value of each point on the image, different atural objects are discerned, accurately obtain the reflectance value of different atural object classifications under the anisotropy face of land; Carry out three-dimensional reconstruction in the homogeneous region based on each atural object inside, improve the precision that the single width image three-dimensional is rebuild, the inverting face of land truly rises and falls, and reduces the error that SFS introduces.

Rajabi (2002) at first introduces the SFS algorithm in the DEM graticule mesh and inserts; Recover face of land 3 D surface shape more accurately based on single image; Cooperate sparse ground control point, can overcome the shortcoming of stereogram and traditional interpolation method, generate the high resolving power DEM that objectively responds landform.Yet Rajabi (2003,2004 years) all is based on gray level image, under the known luminaire condition, resolves SFS reflectogram equation, and with the bilinear interpolation result as SFS algorithm initial value, resultant error is bigger, the reflection landform structural information be easy to generate aliasing ambiguity.

Summary of the invention

The present invention mainly is that single width gray level image light source evaluation method precision is low in the existing SFS technology of solution, does not consider that also image shade and vegetation cover the problem in science such as error that cause for single width remote sensing image three-dimensional reconstruction.This patent expands to the multispectral image field with the light source evaluation method, has proposed a kind of multispectral image that utilizes and has replaced full-colour image, and auxiliary raising of utilization object spectrum information rebuild precision; Cause and set up surface relief and cause the influence with image spectral brightness transformation relation through suppressing surface vegetation cover type and atural object shade simultaneously; Take the SFS three-dimensional reconstruction of face of land spectral information into account and encrypt sparse DEM method.

Another purpose of the present invention is that solution existing in prior technology DEM interpolation algorithm precision commonly used is low, and the DEM on a large scale that can not accomplish based on sparse reference mark sets up technical matters fast; Provide a kind of a kind of SFS three-dimensional reconstruction of taking face of land spectral information into account of setting up fast based on the DEM on a large scale at sparse reference mark to encrypt sparse DEM method.

Above-mentioned technical matters of the present invention mainly is able to solve through following technical proposals:

A kind of SFS three-dimensional reconstruction of taking face of land spectral information into account is encrypted sparse DEM method, it is characterized in that, may further comprise the steps:

Step 1; Be given in gradient and this coordinate points (x of the point set height value on the three-dimensional surface on x direction, y direction under the complicated earth surface coverage condition; Y) locate the relation of spectral brightness value, make up the reflectogram model under the different assumed conditions, and then set up reflectogram equation towards different atural object cover types; Said data comprise: coarse resolution DEM graticule mesh, areal high-resolution remote sensing image;

Step 2, the reflectogram model that step 1 is set up carries out pre-service;

Step 3 to the reflectogram model inversion face of land normal vector of completing steps 3, to different atural object objects, is utilized the face of land normal vector in the object-based Method for minimization inverting atural object;

Step 4 is utilized coarse resolution DEM grid points, recovers the space convergent-divergent yardstick of image and real surface, and utilization slip least square method is found the solution partial points;

Step 5; The space convergent-divergent yardstick that the face of land normal vector that integrating step 3 is obtained and step 4 are recovered image and real surface carries out under the multi-constraint condition inserting in the dem data: in the optimization interpolation process; Introduce the slip least-squares algorithm and carry out interior iterative computation of inserting; And, produce high resolving power DEM graticule mesh with inserting in the edge fit between the different atural object objects of the ICBI algorithm process that keeps structural texture.

A kind of SFS three-dimensional reconstruction of taking face of land spectral information into account above-mentioned is encrypted sparse DEM method, and described reflectogram equation is based on formula:

I _(x，y)＝R(p，q)

And hypothesis:

1) topographical surface is lambert's phantom type;

2) solar source can be regarded infinite distant place pointolite parallel projection as;

3) imaging geometry is rectangular projection;

Wherein, original optics image data is designated as I _{(x, y)}, according to light source parameters (Sx, Sy Sz) are designated as R with the reflectogram that reflectivity albedo recovers to obtain, and p, q are surface direction, be z at x, the derivative of y direction:

Lambert's body reflection model R to embody formula following:

R(p，q)＝albedo*(Sz-p*Sx-q*Sy)/sqrt(p ²+q ²+1)

If surface point (x, y, z _{(x, y)}) normal located is N=(cos α sin β, sin α sin β, cos β), light source direction is L=(cos τ sin γ, sin τ sin γ, cos γ), and τ is the light source drift angle, and γ is the light source inclination angle, and α is the normal drift angle, and β is the normal inclination angle, and following formula is:

$R(p,q)=\mathrm{\η}\frac{\cos \mathrm{\γ}-p\cos \mathrm{\τ}\sin \mathrm{\γ}-q\sin \mathrm{\τ}\sin \mathrm{\γ}}{\sqrt{1+p^2+q^2}}$

$=\mathrm{\η}(\cos (\mathrm{\α}-\mathrm{\τ})\sin \mathrm{\β}\sin \mathrm{\γ}+\cos \mathrm{\β}\cos \mathrm{\γ}).$

A kind of SFS three-dimensional reconstruction of taking face of land spectral information into account above-mentioned is encrypted sparse DEM method, and in the described step 2, the reflectogram model that step 1 is set up carries out pre-service and may further comprise the steps:

Step 3.1 is carried out cutting and registration to multispectral remote sensing image and dem data, makes the corresponding location matches of shadow information;

Whether step 3.2 has the spectral brightness of existence on the inspection multispectral remote sensing image, and promptly the DN value is zero zone, because the graticule mesh of zero luminance lacks second-order properties, is regarded as shade, carries out this area image reparation in step 3.4.

Step 3.3, the atural object segmentation and classification that the remote sensing image of completing steps 3.2 is carried out based on image spectrum and texture; For cutting apart back image, the surface reflectivity of estimation light source direction and each object;

Step 3.4 utilizes the constant theoretical model of color to carry out image repair, the filtering shadow region to the remote sensing image of completing steps 3.3;

Step 3.5, the remote sensing image of completing steps 3.4 is carried out under the anisotropy condition light and shade recover shape: add the shade of gray constraint in the object, the slickness constraint retrains with integrability, carry out the variational method and fast the method for advancing resolve, obtain face of land relative altitude.

The proposition of the invention property utilizes the multispectral remote sensing image to combine the coarse resolution DEM graticule mesh interpolating method of SFS three-dimensional reconstruction result; Proposition utilizes the reflectivity of the different types of ground objects of object spectrum data-evaluation, removes the surface vegetation cover type is rebuild precision to SFS influence; Proposition reduces the influence of image shade to the atural object spectral brightness based on the single width multispectral remote sensing image three-dimensional method for reconstructing of shadow removal; Set up the light source estimation of anisotropy complicated earth surface, introduce multispectral information, add the ground additional control point, improve the light source direction estimation precision; Based on the result of three-dimensional reconstruction, utilize the slip least square to carry out the elevation value interpolation, cooperate the atural object edge fit method of inserting in the ICBI, avoided the level and smooth phenomenon of mistake of DEM interpolating method, improve DEM and encrypt precision.The method that patent proposes is mainly used in the DEM graticule mesh in the deficient area of high-resolution three-dimension modeling, depopulated zone or data of the quick foundation of 3-dimensional digital terrain model (DSM), complex region under the case of emergency and sets up.

Encrypt sparse DEM method at above-mentioned a kind of SFS three-dimensional reconstruction of taking face of land spectral information into account, in the described step 3.1, the concrete steps of remote sensing image and dem data being carried out cutting and registration are following:

1) reads dem data and remote sensing image data file header, obtain the coverage of its data;

2) set up the affined transformation model

$\left(\begin{array}{c}{x}^{\′}\\ y^{\′}\end{array}\right)=\left(\begin{array}{cc}{a}_{11}& a_{12}\\ a_{21}& a_{\mathrm{xx}}\end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right)+\left(\begin{array}{c}{t}_x\\ t_y\end{array}\right);$

3) according to image four angle point value ranges up and down, bring the affined transformation model into, obtain eight equations, find the solution 6 transformation model parameters;

4) geographic coordinate after the conversion of calculating every bit.

A kind of SFS three-dimensional reconstruction of taking face of land spectral information into account above-mentioned is encrypted sparse DEM method, and in the described step 3.3, the concrete steps of atural object segmentation and classification are following:

1) combines the atural object texture features according to spectral information, definition classification number c, cost function threshold epsilon, maximum iteration time t; The instance (see figure 2) that this patent provides in this step, definition classification number is 3, and the cost function threshold value is 0.1, and maximum iteration time is 50;

2) be subordinate to matrix U with the random number initialization of value between 0-1, make it satisfy the constraint condition in the following formula:

$\underset{i=1}{\overset{c}{\mathrm{\Σ}}}{u}_{\mathrm{ij}}=1,\∀j=1,...,n$

3) with computes c cluster centre ci, i=1 ..., c:

$c_i=\frac{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{u}_{\mathrm{ij}}^m{x}_j}{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{u}_{\mathrm{ij}}^m}$

Wherein, x _i(i=1,2 ..., n) be n sample point spectrum vector;

4) according to the computes cost function: if cost function less than the threshold epsilon of confirming, or the change amount of cost function relative last time of cost function value is less than threshold value, then algorithm stops:

$J(U,{\mathrm{<figure-callout\; id="1"\; label="c"\; filenames="HDA0000065913090000022.png,HDA0000065913090000032.png"\; state="\{\{state\}\}">c</figure-callout>}}_1,...,c_c)=\underset{i=1}{\overset{c}{\mathrm{\&Sigma;}}}{J}_i=\underset{i=1}{\overset{c}{\mathrm{\&Sigma;}}}\underset{j}{\overset{n}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ij}}^m{d}_{\mathrm{ij}}^2$

Here u _IjBetween 0 to 1; c _iBe the cluster centre of ambiguity group I, d _Ij=|| c _i-x _j|| be the Euclidean distance between i cluster centre and j data points; And m ∈ [1, ∞) be a weighted index, the m value is 2 here;

If cost function is greater than threshold value, then with the new U matrix of computes; Return step 2);

$u_{\mathrm{ij}}=\frac{1}{\underset{k=1}{\overset{c}{\mathrm{\&Sigma;}}}{\left(\frac{d_{\mathrm{ij}}}{d_{\mathrm{kj}}}\right)}^{2/(m-1)}};$

5) if iterations surpasses threshold value t, then calculating stops, and finally obtains the terrain classification figure of multispectral remote sensing image.

A kind of SFS three-dimensional reconstruction of taking face of land spectral information into account above-mentioned is encrypted sparse DEM method, and in the described step 3.3, to cutting apart the back image, the concrete steps of estimation light source direction are following:

1) utilize the linear function conversion that the brightness of image value is carried out normalization:

I _i′＝I _i/I _max

I wherein _iBe a certain wave band original brightness value at i point place, I _MaxBe this wave band overall situation image brilliance maximal value, I _i' be the brightness value after the normalization;

2) get the original DEM graticule mesh height value of same classification atural object and the brightness value on the correspondence remote sensing image of position thereof,, set up a plurality of reflectogram equations as the reference mark:

I _i＝albedo*(Sz-p _i*Sx-q _i*Sy)/sqrt(p _i ²+q _i ²+1)；

Wherein i ∈ N is a DEM grid points number, I _iBe the brightness value of this point, (p q) is the x of this point, the gradient of y direction, and for same type of atural object, reflectivity albedo is a constant; Promptly utilize N equation (N＞3) resolve 3 light source parameters (Sx, Sy, Sz); Utilize least square method to obtain estimates of parameters.

A kind of SFS three-dimensional reconstruction of taking face of land spectral information into account above-mentioned is encrypted sparse DEM method, and in the described step 3.3, the concrete steps of surface reflectivity of obtaining each object are following: wherein, each object is meant different atural object classifications:

$\mathrm{\&tau;}=\arctan \left(\frac{E\left\{\frac{y_L}{\sqrt{x_L^2+y_L^2}}\right\}}{E\left\{\frac{x_L}{\sqrt{x_L^2+y_L^2}}\right\}}\right);$

Wherein E is a population mean, $\left(\begin{array}{c}{x}_L\\ y_L\end{array}\right)={\left(B^{\&prime;}B\right)}^{-1}{B}^{\&prime;}\mathrm{DI},\mathrm{DI}=\left(\begin{array}{c}\mathrm{\&delta;}{I}_1\\ \mathrm{\&delta;}{I}_2\\ M\\ \mathrm{\&delta;}{I}_N\end{array}\right),B=\left(\begin{array}{cc}\mathrm{\&delta;}{x}_1& \mathrm{\&delta;}{\mathrm{<figure-callout\; id="1"\; label="y"\; filenames="HDA0000065913090000022.png,HDA0000065913090000032.png"\; state="\{\{state\}\}">y</figure-callout>}}_1\\ \mathrm{\&delta;}{x}_2& \mathrm{\&delta;}{y}_2\\ M& M\\ \mathrm{\&delta;}{x}_N& \mathrm{\&delta;}{y}_N\end{array}\right)$

$B^{\&prime;}=\left(\begin{array}{cccccccc}1& \frac{\sqrt{2}}{2}& 0& -\frac{\sqrt{2}}{2}& -1& -\frac{\sqrt{2}}{2}& 0& \frac{\sqrt{2}}{2}\\ 0& \frac{\sqrt{2}}{2}& 1& \frac{\sqrt{2}}{2}& 0& -\frac{\sqrt{2}}{2}& -1& -\frac{\sqrt{2}}{2}\end{array}\right);$

Order:

$f_1\left(\mathrm{\&gamma;}\right)\&ap;\underset{i=0}{\overset{7}{\mathrm{\&Sigma;}}}{a}_i{\cos }^i\mathrm{\&gamma;}$

$f_2\left(\mathrm{\&gamma;}\right)\&ap;\underset{i=0}{\overset{7}{\mathrm{\&Sigma;}}}{b}_i{\cos }^i\mathrm{\&gamma;},$

$f_3\left(\mathrm{\&gamma;}\right)\&ap;\underset{i=0}{\overset{7}{\mathrm{\&Sigma;}}}{c}_i{\cos }^i\mathrm{\&gamma;}$

If $\frac{E\left\{I\right\}}{\sqrt{E\left\{I^2\right\}}}>f_3\left(0\right)=0.9691,$ γ=0 then; Otherwise γ does $\frac{E\left\{I\right\}}{\sqrt{E\left\{I^2\right\}}}=f_3\left(\mathrm{\&gamma;}\right)$ Separate;

For reflectivity: $\mathrm{\&eta;}=\frac{1}{{\mathrm{<figure-callout\; id="1"\; label="f"\; filenames="HDA0000065913090000022.png,HDA0000065913090000032.png"\; state="\{\{state\}\}">f</figure-callout>}}_1^2\left(\mathrm{\&gamma;}\right)+f_2\left(\mathrm{\&gamma;}\right)}\&CenterDot;\left\{E\right\{I\}\&CenterDot;f_1\left(\mathrm{\&gamma;}\right)+\sqrt{E\left\{I^2\right\}\&CenterDot;f_2\left(\mathrm{\&gamma;}\right)}\},$ Wherein, I is the normalization brightness value of same type atural object on the image.

A kind of SFS three-dimensional reconstruction of taking face of land spectral information into account above-mentioned is encrypted sparse DEM method, in the described step 3.4, carries out image repair, and the concrete steps of filtering shadow region are following:

1) carrying out the shadow region according to each wave band DN value of multispectral remote sensing image and structural information detects: think that it is the shadow region that reflectivity on all wave bands all approaches 0 zone;

2) read in original image I and detected shadow region template, according to template with the original image separated into two parts: the shadow region image I _sWith the non-hatched area image I _Ns

3) respectively to the shadow region image I _sWith the non-hatched area image I _NsCarry out color constancy with Shade of Gray algorithm and calculate, obtain shadow region illumination color and estimate [s_r, s_g, s_b] and nonshaded area illumination color estimation [n_r, n_g, n_b], computing formula is:

$e=k{\left(\frac{\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\left(f{(i,j)}^p\right)}{M*N}\right)}^{1/p}$

In the formula, the light source colour in

expression present image, f={R; G; The DN value of each wave band of B} presentation video is assumed to be red, green here; Blue three wave bands, p=6;

The correcting gain of each wave band is:

4) shadow region transform to non-hatched area the ratio constant vector is

5) ratio constant vector k is transformed into number space

6) the shadow region image I _sBe transformed into the logarithm Space L _s=log (I _s);

7) obtain removing image behind the shade, I=exp (L _s+ c)+I _Ns

A kind of SFS three-dimensional reconstruction of taking face of land spectral information into account above-mentioned is encrypted sparse DEM method, and the concrete steps of carrying out light and shade recovery shape are following:

1) SFS reflectogram equation is rewritten as the function of face of land gradient, and adds the brightness constraint, conditions such as slickness constraint are come solving equation; With the global error structure cost functional of image and reflectogram, and use the variational method to find the solution the three-dimensional surface function;

Possibly there is error between the determined body surface brightness of image of known image gray scale and reflection model, luminance equation is become the form of error function, be i.e. brightness constraint:

$I_1=\underset{\mathrm{\&Omega;}}{\&Integral;\&Integral;}{(E(x,y)-R(p(x,y),q(x,y)))}^2\mathrm{dxdy};$

2) introduce shade of gray constraint and integrability constraint, set up error equation:

∫∫F(p，q，z)dxdy＝∫∫((I-R) ²+(R _x-I _x) ²+(R _y-I _y) ²+μ((z _x-p) ²+(z _y-q) ²))dxdy；

3) utilize Eulerian equation to ask the functional minimal value, through successive iteration, (p is q) with elevation z to obtain the true normal vector in the face of land.

A kind of SFS three-dimensional reconstruction of taking face of land spectral information into account above-mentioned is encrypted sparse DEM method, in the described step 5, carries out following two steps:

Step 5.1 for same atural object object, has promptly kept the face of land changes in the object continuity and slickness, utilizes the result of remote sensing image three-dimensional reconstruction, obtains surface slope aspect figure according to formula:

$\mathrm{slope}=\arctan \sqrt{{f_x}^2+{f_y}^2}*180/\mathrm{\&pi;}$

$\mathrm{aspect}=\arctan \left(\frac{f_x}{f_y}\right)$

Wherein, (f _x, f _y) be the relative height of three-dimensional reconstruction at x, the Grad on the y direction;

Set up cost function, insert the difference minimum of the surface slope aspect figure of back elevation graticule mesh surface slope aspect and three-dimensional reconstruction in making, with moving least squares method iterative calculation DEM graticule mesh;

Step 5.2 is if adjacent twice iteration result's difference less than threshold value 0.1, then with inserting in the edge fit between the different atural object objects of ICBI algorithm process, produces high resolving power DEM graticule mesh, that is:

If the even number site in the graticule mesh is the elevation known point, for the position be (2i+1,2j+1) point at first calculate the value on the both direction:

I ₁₁(2i+1，2j+1)＝I(2i-2，2j+2)+I(2i，2j)+I(2i+2，2j-2)

-3I(2i，2j+2)-3I(2i+2，2j)+I(2i，2j+4)+I(2i+2，2j+2)+I(2i+4，2j)；

I ₂₂(2i+1，2j+1)＝I(2i，2j-2)+I(2i+2，2j)+I(2i+4，2j+2)

-3I(2i，2j)-3I(2i+2，2j+2)+I(2i-2，2j)+I(2i，2j+2)+I(2i+2，2j+4)；

Final its interpolate value is pressed following formula and is confirmed:

Therefore, the present invention has following advantage:

1. utilize multispectral image to replace panchromatic image, the auxiliary accuracy of reconstruction that improves of utilization object spectrum information; Can effectively suppress surface vegetation cover type, atural object shade for the influence of setting up face of land three-dimensional model and image spectral brightness transformation relation;

2. be significantly increased with respect to DEM interpolation algorithm interpolation precision commonly used, case verification can improve 2 times with coarse grid DEM spatial resolution, and the DEM on a large scale that can accomplish based on sparse control point sets up fast.

Description of drawings

Fig. 1 workflow diagram of the present invention.

Fig. 2 a1 is the pilot region multi-spectral remote sensing image;

Fig. 2 a2 is the original DEM image of pilot region (2.5 meters);

Fig. 2 a3 is a classification results, is respectively rock, vegetation, water body;

Fig. 2 b1 is a this patent SFS arithmetic result;

Fig. 2 b2 is traditional SFS arithmetic result;

Fig. 2 c1 is a DEM encrypted result of the present invention;

Fig. 2 c2 is the 3-D display of figure Fig. 2 c1;

Fig. 2 d1 encrypts graticule mesh for the DEM that inserts DEM result in the weighted least-squares surface fitting of the present invention;

Fig. 2 d2 is the 3-D display of Fig. 2 d1;

Fig. 2 e1 is down that the distance weighted interior DEM result's of inserting DEM encrypts graticule mesh;

Fig. 2 e2 is the 3-D display of figure Fig. 2 e1;

Fig. 2 f1 encrypts for golden interpolation DEM in the gram;

Fig. 2 f2 is the 3-D display DEM graticule mesh of figure Fig. 2 f1;

Fig. 3 a generates the DEM elevation relatively for each interpolation algorithm;

Fig. 3 b1 is the level line of algorithm of the present invention;

Fig. 3 b2 is golden algorithm level line in the gram;

Embodiment

Pass through embodiment below, and combine accompanying drawing, do further bright specifically technical scheme of the present invention.

Embodiment:

As shown in Figure 1, a kind of SFS three-dimensional reconstruction of taking face of land spectral information into account is encrypted sparse DEM method, may further comprise the steps:

Step 1; Be given in gradient and the coordinate points (x of the point set height value on the three-dimensional surface on x direction, y direction under the complicated earth surface coverage condition; Y) locate the relation of spectral brightness value, make up the reflectogram model under the different assumed conditions, and then set up reflectogram equation towards different atural object cover types; Said data comprise: coarse resolution DEM graticule mesh, areal high-resolution remote sensing image; The reflectogram equation is based on formula:

I _(x，y)＝R(p，q)

And hypothesis:

1) topographical surface is lambert's phantom type;

2) solar source can be regarded infinite distant place pointolite parallel projection as;

3) imaging geometry is rectangular projection;

Wherein, original optics image data is designated as I _{(x, y)}, according to light source parameters (Sx, Sy Sz) are designated as R with the reflectogram that reflectivity albedo recovers to obtain, and p, q are surface direction, be z at x, the derivative of y direction:

Lambert's body reflection model R to embody formula following:

R(p，q)＝albedo*(Sz-p*Sx-q*Sy)/sqrt(p ²+q ²+1)

If surface point (x, y, z _{(x, y)}) normal located is N=(cos α sin β, sin α sin β, cos β), light source direction is L=(cos τ sin γ, sin τ sin γ, cos γ), and τ is the light source drift angle, and γ is the light source inclination angle, and α is the normal drift angle, and β is the normal inclination angle, and following formula is:

$R(p,q)=\mathrm{\&eta;}\frac{\cos \mathrm{\&gamma;}-p\cos \mathrm{\&tau;}\sin \mathrm{\&gamma;}-q\sin \mathrm{\&tau;}\sin \mathrm{\&gamma;}}{\sqrt{1+p^2+q^2}}$

$=\mathrm{\&eta;}(\cos (\mathrm{\&alpha;}-\mathrm{\&tau;})\sin \mathrm{\&beta;}\sin \mathrm{\&gamma;}+\cos \mathrm{\&beta;}\cos \mathrm{\&gamma;})$

Step 2, the reflectogram model that step 1 is set up carries out pre-service; Pre-service may further comprise the steps:

Step 3.1 is carried out cutting and registration to multispectral remote sensing image and dem data, makes the corresponding location matches of shadow information; Concrete steps are following:

1) reads dem data and remote sensing image data file header, obtain the coverage of its data;

2) set up the affined transformation model

$\left(\begin{array}{c}{x}^{\&prime;}\\ y^{\&prime;}\end{array}\right)=\left(\begin{array}{cc}{a}_{11}& a_{12}\\ a_{21}& a_{\mathrm{xx}}\end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right)+\left(\begin{array}{c}{t}_x\\ t_y\end{array}\right);$

3) according to image four angle point value ranges up and down, bring the affined transformation model into, obtain eight equations, find the solution 6 transformation model parameters;

4) geographic coordinate after the conversion of calculating every bit.

Step 3.2, whether having spectral brightness (being the DN value) on the inspection multispectral remote sensing image is zero zone, because the graticule mesh of zero luminance lacks second-order properties, is regarded as shade, carries out this area image reparation in step 3.4;

Step 3.3, the atural object segmentation and classification that the remote sensing image of completing steps 3.2 is carried out based on image spectrum and texture; For cutting apart back image, the surface reflectivity of estimation light source direction and each object; The concrete steps of atural object segmentation and classification are following:

1) according to spectral information and atural object texture features definition classification number c, cost function threshold epsilon, maximum iteration time t; The example definition classification number that this patent provides in this step is 3, and the cost function threshold value is 0.1, and maximum iteration time is 50;

2) be subordinate to matrix U with the random number initialization that is worth at 0,1, make it satisfy the constraint condition in the following formula:

$\underset{i=1}{\overset{c}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ij}}=1,\&ForAll;j=1,...,n$

3) with computes c cluster centre ci, i=1 ..., c:

$c_i=\frac{\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ij}}^m{x}_j}{\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ij}}^m}$

Wherein, x _i(i=1,2 ..., n) be n sample point spectrum vector;

4) according to the computes cost function: if cost function less than the threshold epsilon of confirming, or the change amount of cost function relative last time of cost function value is less than threshold value, then algorithm stops:

$J(U,{\mathrm{<figure-callout\; id="1"\; label="c"\; filenames="HDA0000065913090000022.png,HDA0000065913090000032.png"\; state="\{\{state\}\}">c</figure-callout>}}_1,...,c_c)=\underset{i=1}{\overset{c}{\mathrm{\&Sigma;}}}{J}_i=\underset{i=1}{\overset{c}{\mathrm{\&Sigma;}}}\underset{j}{\overset{n}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ij}}^m{d}_{\mathrm{ij}}^2$

Here u _IjBetween 0 to 1; c _iBe the cluster centre of ambiguity group I, d _Ij=|| c _i-x _j|| be the Euclidean distance between i cluster centre and j data points; And m ∈ [1, ∞) be a weighted index, the m value is 2 here;

If cost function is greater than threshold value, then with the new U matrix of computes; Return step 2);

$u_{\mathrm{ij}}=\frac{1}{\underset{k=1}{\overset{c}{\mathrm{\&Sigma;}}}{\left(\frac{d_{\mathrm{ij}}}{d_{\mathrm{kj}}}\right)}^{2/(m-1)}};$

5) if iterations surpasses threshold value t, then calculating stops, and finally obtains the terrain classification figure of multispectral remote sensing image.

Then, to cutting apart the back image, the concrete steps of estimation light source direction are following:

1) utilize the linear function conversion that the brightness of image value is carried out normalization:

I _i′＝I _i/I _max

I wherein _iBe a certain wave band original brightness value at i point place, I _MaxBe this wave band overall situation image brilliance maximal value, I _i' be the brightness value after the normalization;

2) get the original DEM graticule mesh height value of same classification atural object and the brightness value on the correspondence remote sensing image of position thereof,, set up a plurality of reflectogram equations as the reference mark:

I _i＝albedo*(Sz-p _i*Sx-q _i*Sy)/sqrt(p _i ²+q _i ²+1)；

Wherein i ∈ N is a DEM grid points number, I _iBe the brightness value of this point, (p q) is the x of this point, the gradient of y direction, and for same type of atural object, reflectivity albedo is a constant; Promptly utilize N equation (N＞3) resolve 3 light source parameters (Sx, Sy, Sz); Utilize least square method to obtain estimates of parameters.

At last, it is following to obtain the concrete steps of surface reflectivity of each object: wherein, each object is meant different atural object classifications:

$\mathrm{\&tau;}=\arctan \left(\frac{E\left\{\frac{y_L}{\sqrt{x_L^2+y_L^2}}\right\}}{E\left\{\frac{x_L}{\sqrt{x_L^2+y_L^2}}\right\}}\right);$

Wherein E is a population mean, $\left(\begin{array}{c}{x}_L\\ y_L\end{array}\right)={\left(B^{\&prime;}B\right)}^{-1}{B}^{\&prime;}\mathrm{DI},\mathrm{DI}=\left(\begin{array}{c}\mathrm{\&delta;}{I}_1\\ \mathrm{\&delta;}{I}_2\\ M\\ \mathrm{\&delta;}{I}_N\end{array}\right),B=\left(\begin{array}{cc}\mathrm{\&delta;}{x}_1& \mathrm{\&delta;}{\mathrm{<figure-callout\; id="1"\; label="y"\; filenames="HDA0000065913090000022.png,HDA0000065913090000032.png"\; state="\{\{state\}\}">y</figure-callout>}}_1\\ \mathrm{\&delta;}{x}_2& \mathrm{\&delta;}{y}_2\\ M& M\\ \mathrm{\&delta;}{x}_N& \mathrm{\&delta;}{y}_N\end{array}\right)$

$B^{\&prime;}=\left(\begin{array}{cccccccc}1& \frac{\sqrt{2}}{2}& 0& -\frac{\sqrt{2}}{2}& -1& -\frac{\sqrt{2}}{2}& 0& \frac{\sqrt{2}}{2}\\ 0& \frac{\sqrt{2}}{2}& 1& \frac{\sqrt{2}}{2}& 0& -\frac{\sqrt{2}}{2}& -1& -\frac{\sqrt{2}}{2}\end{array}\right);$

Order:

$f_1\left(\mathrm{\&gamma;}\right)\&ap;\underset{i=0}{\overset{7}{\mathrm{\&Sigma;}}}{a}_i{\cos }^i\mathrm{\&gamma;}$

$f_2\left(\mathrm{\&gamma;}\right)\&ap;\underset{i=0}{\overset{7}{\mathrm{\&Sigma;}}}{b}_i{\cos }^i\mathrm{\&gamma;},$

$f_3\left(\mathrm{\&gamma;}\right)\&ap;\underset{i=0}{\overset{7}{\mathrm{\&Sigma;}}}{c}_i{\cos }^i\mathrm{\&gamma;}$

If $\frac{E\left\{I\right\}}{\sqrt{E\left\{I^2\right\}}}>f_3\left(0\right)=0.9691,$ γ=0 then; Otherwise γ does $\frac{E\left\{I\right\}}{\sqrt{E\left\{I^2\right\}}}=f_3\left(\mathrm{\&gamma;}\right)$ Separate;

For reflectivity: $\mathrm{\&eta;}=\frac{1}{{\mathrm{<figure-callout\; id="1"\; label="f"\; filenames="HDA0000065913090000022.png,HDA0000065913090000032.png"\; state="\{\{state\}\}">f</figure-callout>}}_1^2\left(\mathrm{\&gamma;}\right)+f_2\left(\mathrm{\&gamma;}\right)}\&CenterDot;\left\{E\right\{I\}\&CenterDot;f_1\left(\mathrm{\&gamma;}\right)+\sqrt{E\left\{I^2\right\}\&CenterDot;f_2\left(\mathrm{\&gamma;}\right)}\},$ Wherein, I is the normalization brightness value of same type atural object on the image.

Step 3.4 utilizes the constant theoretical model of color to carry out image repair, the filtering shadow region to the remote sensing image of completing steps 3.3; Concrete steps are following:

1) carrying out the shadow region according to each wave band DN value of multispectral remote sensing image and structural information detects: think that it is the shadow region that reflectivity on all wave bands all approaches 0 zone;

2) read in original image I and detected shadow region template, according to template with the original image separated into two parts: the shadow region image I _sWith the non-hatched area image I _Ns

3) respectively to the shadow region image I _sWith the non-hatched area image I _NsCarry out color constancy with Shade of Gray algorithm and calculate, obtain shadow region illumination color and estimate [s_r, s_g, s_b] and nonshaded area illumination color estimation [n_r, n_g, n_b], computing formula is:

$e=k{\left(\frac{\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\left(f{(i,j)}^p\right)}{M*N}\right)}^{1/p}$

In the formula, the light source colour in

expression present image, f={R; G; The spectral brightness DN value of each wave band of B} presentation video is assumed to be red, green here; Blue three wave bands, p=6;

The correcting gain of each wave band is:

4) shadow region transform to non-hatched area the ratio constant vector is

5) ratio constant vector k is transformed into number space

6) the shadow region image I _sBe transformed into the logarithm Space L _s=log (I _s);

7) obtain removing image behind the shade, I=exp (L _s+ c)+I _Ns

Step 3.5; The remote sensing image of completing steps 3.4 is carried out under the anisotropy condition light and shade recover shape: add the shade of gray constraint in the object; Slickness constraint and integrability constraint; Carry out the variational method and fast the method for advancing resolve, obtain face of land relative altitude, it is following that light and shade is recovered the shape concrete steps:

1) SFS reflectogram equation is rewritten as the function of face of land gradient, and adds the brightness constraint, conditions such as slickness constraint are come solving equation; With the global error structure cost functional of image and reflectogram, and use the variational method to find the solution the three-dimensional surface function;

Possibly there is error between the determined body surface brightness of image of known image gray scale and reflection model, luminance equation is become the form of error function, be i.e. brightness constraint:

$I_1=\underset{\mathrm{\&Omega;}}{\&Integral;\&Integral;}{(E(x,y)-R(p(x,y),q(x,y)))}^2\mathrm{dxdy};$

2) introduce shade of gray constraint and integrability constraint, set up error equation:

∫∫F(p，q，z)dxdy＝∫∫((I-R) ²+(R _x-I _x) ²+(R _y-I _y) ²+μ((z _x-p) ²+(z _y-q) ²))dxdy；

3) utilize Eulerian equation to ask the functional minimal value, through successive iteration, (p is q) with elevation z to obtain the true normal vector in the face of land.

Step 3 to the reflectogram model inversion face of land normal vector of completing steps 3, to different atural object objects, is utilized the face of land normal vector in the object-based Method for minimization inverting atural object;

Step 4 is utilized coarse resolution DEM grid points, recovers the space convergent-divergent yardstick of image and real surface, and utilization slip least square method is found the solution partial points;

Step 5; The space convergent-divergent yardstick that the face of land normal vector that integrating step 3 is obtained and step 4 are recovered image and real surface carries out under the multi-constraint condition inserting in the dem data: in the optimization interpolation process; Introduce the slip least-squares algorithm and carry out interior iterative computation of inserting; And with inserting in the edge fit between the ICBI algorithm process difference atural object objects that keep structural texture, produce high resolving power DEM graticule mesh, following two steps of the concrete execution of this step:

Step 5.1 for same atural object object, has promptly kept the face of land changes in the object continuity and slickness, utilizes the result of remote sensing image three-dimensional reconstruction, obtains surface slope aspect figure according to formula:

$\mathrm{slope}=\arctan \sqrt{{f_x}^2+{f_y}^2}*180/\mathrm{\&pi;}$

$\mathrm{aspect}=\arctan \left(\frac{f_x}{f_y}\right)$

Wherein, (f _x, f _y) be the relative height of three-dimensional reconstruction at x, the Grad on the y direction;

Set up cost function, insert the difference minimum of the surface slope aspect figure of back elevation graticule mesh surface slope aspect and three-dimensional reconstruction in making, with moving least squares method iterative calculation DEM graticule mesh;

Step 5.2 is if adjacent twice iteration result's difference less than threshold value 0.1, then with inserting in the edge fit between the different atural object objects of ICBI algorithm process, produces high resolving power DEM graticule mesh, that is:

If the even number site in the graticule mesh is the elevation known point, for the position be (2i+1,2j+1) point at first calculate the value on the both direction:

I ₁₁(2i+1，2j+1)＝I(2i-2，2j+2)+I(2i，2j)+I(2i+2，2j-2)

-3I(2i，2j+2)-3I(2i+2，2j)+I(2i，2j+4)+I(2i+2，2j+2)+I(2i+4，2j)；

I ₂₂(2i+1，2j+1)＝I(2i，2j-2)+I(2i+2，2j)+I(2i+4，2j+2)

-3I(2i，2j)-3I(2i+2，2j+2)+I(2i-2，2j)+I(2i，2j+2)+I(2i+2，2j+4)；

Final its interpolate value is pressed following formula and is confirmed:

As shown in Figure 2, adopt the SPOT remote sensing images and the DEM graticule mesh image of registration.The remote sensing images spatial resolution is 5m, and original DEM graticule mesh resolution is 2.5m.Being the verification algorithm precision, is 10m with DEM graticule mesh down-sampling to spatial resolution, and returns to 5m resolution, comparative analysis reconstructed results, the validity of verification method with the method for remote sensing image three-dimensional reconstruction and the method for DEM interpolation commonly used.At first to the remote sensing image classification that exercises supervision, the face of land covered be divided into vegetation, water body, rock and zone, farmland.For different face of land cover types, respectively with statistical information estimation clutter reflections rate between class and light source position estimation.For different cover types zone, carry out OO SFS and recover inverting face of land relative relief state.According to surface slope aspect information, utilize least square to combine to insert in the ICBI method high resolving power DEM graticule mesh again.And with least square fitting, the distance weighted interpolation of falling, the result of spline interpolation compares analysis, and following table is single width image three-dimensional method for reconstructing of the present invention and the contrast of traditional interpolation algorithm precision:

	Equal error	Maximum error	Least error	Absolute all errors
					SFS rebuilds	-0.1399	13	-38	0.4673
Least square fitting	0.2574	66	-69	1.7914
					Distance weighted inverse ratio interpolation	-0.2638	12	-68	0.7487
Golden interpolation in the gram	0.6373	22	-21	1.3490

Left side figure is a SFS algorithm image restored relative height of the present invention among Fig. 2 b, and right figure is the relative height that traditional SFS algorithm recovers.At first for three types of atural objects that exist on the remote sensing image among the figure (a): rock, vegetation, water body, this paper are estimated three's earth surface reflection rate value respectively, the gradation of image difference of in SFS resolves, bringing with face of land reflection parameters filtering atural object cover type difference.Therefore SFS algorithm of the present invention can obviously separate vegetation and rock place, calculates vegetation area and has relatively low height value for the rock zone on image top.And traditional SFS algorithm causes vegetation area and rock zone sliding linking to each other of height flat among the result owing to do not consider the influence of atural object cover type, catches less than the elevation that changes between classification.

At Fig. 2 c is the DEM graticule mesh of rebuilding based on algorithm of the present invention, can see that from 3-D display the DEM details after the reconstruction is outstanding, and the topographic structure in the original dem data comprises ridge, the mountain valley, and the water body pond can both clearly be seen, architectural feature keeps better.

The result of Fig. 2 d and e shows that there was level and smooth problem in distance weighted inverse ratio interpolation, lost a lot of surface relief details.Then there is vibration in least square fitting, causes fake information.Simultaneously find out that by errors table-1 both compare with algorithm of the present invention, have bigger minimax error, differ greatly with true high resolving power DEM graticule mesh.

Golden interpolation result is shown in Fig. 2 f in the gram; Can see that the DEM morphosis is more clear; It is lighter to cross level and smooth phenomenon than least square fitting and distance weighted inverse ratio interpolation; But still the result not as the inventive method is abundant in conjunction with the detailed information that can find out the result that golden DEM encrypts in the gram among Fig. 3 a and the 3b, and this is because this patent algorithm has increased the high-resolution remote sensing image spectral information.

Can find out that from Fig. 3 a least square fitting exists vibration and fake information, especially when noise exists, disturb very serious.There was level and smooth problem in distance weighted inverse ratio interpolation, lost a lot of surface relief details, and table-1 finds out that both compare with algorithm of the present invention simultaneously, has bigger minimum and maximum error, differs greatly with true high resolving power DEM graticule mesh.In conjunction with golden interpolation result in the gram, (shown in the f, can see that the DEM morphosis is more clear, it is lighter than least square fitting and distance weighted inverse ratio interpolation to cross level and smooth phenomenon, but also has level and smooth phenomenon to a certain degree like figure.Can find out from Fig. 3, Fig. 2 and table-1; The present invention has increased single width remote sensing image and the constraint of face of land spectral information; Utilize the SFS three-dimensional reconstruction to encrypt sparse DEM method, can recover more rich details information, significant advantage is relatively arranged with method such as golden interpolation algorithm in the gram.

The present invention can set up the complicated earth surface three-dimensional model fast; Be mainly used in high resolving power DEM graticule mesh and obtain fast, the three-dimension terrian model is set up fast under case of emergency, the application of satellite remote-sensing image Target Recognition and remote sensing image measurement aspect; The high-resolution three-dimension modeling of complex region; The DEM graticule mesh in depopulated zone is set up, the emergency landing of aircraft under MODEL OVER COMPLEX TOPOGRAPHY, and three-dimensional modeling data compression aspect, the face of land is significant.

The present invention is that a kind of SFS three-dimensional reconstruction of taking face of land spectral information into account is encrypted sparse DEM method; Specifically comprise:, utilize heterogeneous face of land spectral characteristic to remove the influence of the grey scale change that the surface vegetation cover type causes surface relief based on the single width remote sensing image three-dimensional reconstruction of anisotropy complicated earth surface.For different atural object type of topical estimation earth surface reflection rate values, realized the light source position evaluation method of multispectral image.Remove atural object and block shade, avoid the image spectrum brightness value can not truly reflect the situation of surface relief.Clearly set up complicated earth surface gradation of image and light source direction; Relation between face of land normal vector and the earth surface reflection rate; Utilize the slip least square method to find the solution partial points and ICBI interpolation method; Utilize coarse resolution DEM grid points, recover the convergent-divergent yardstick of image and real surface, resolve high-density gridding DEM elevation.

The present invention is used under the case of emergency three-dimension terrian model and sets up fast; Aspects such as identification of satellite remote-sensing image military target and remote sensing image measurement; High-resolution three-dimension modeling in the complex-terrain zone; The three-dimensional integratedization demonstration of the foundation of the graticule mesh of DEM on a large scale, the data compression of DEM graticule mesh, DEM and the image of data starvation areas, the aspects such as emergency landing of aircraft under complex-terrain also has wide practical use.

The described specific embodiment of this patent only is that inventive concept is illustrated.The technical staff of the technical field of the invention can make various modifications or additional to described specific embodiment, or adopts similar mode to substitute, but can't depart from spirit of the present invention or surmount the defined scope of appended claims.

Claims (10)

1. a SFS three-dimensional reconstruction of taking face of land spectral information into account is encrypted sparse DEM method, it is characterized in that, may further comprise the steps:

Step 1; Be given in gradient and the coordinate points (x of the point set height value on the three-dimensional surface on x direction, y direction under the complicated earth surface coverage condition; Y) relation of spectral brightness value makes up the reflectogram model under the different assumed conditions, and then sets up the reflectogram equation towards different atural object cover types; Said data comprise: coarse resolution DEM graticule mesh, areal high-resolution remote sensing image;

Step 2, the reflectogram model that step 1 is set up carries out pre-service;

Step 3 to the reflectogram model inversion face of land normal vector of completing steps 2, to different atural object objects, is utilized the face of land normal vector in the object-based Method for minimization inverting atural object;

Step 4 is utilized coarse resolution DEM grid points, recovers the space convergent-divergent yardstick of image and real surface, and utilization slip least square method is found the solution partial points;

Step 5; The space convergent-divergent yardstick that the face of land normal vector that integrating step 3 is obtained and step 4 are recovered image and real surface carries out under the multi-constraint condition inserting in the dem data: in the optimization interpolation process; Introduce the slip least-squares algorithm and carry out interior iterative computation of inserting; And, produce high resolving power DEM graticule mesh with inserting in the edge fit between the different atural object objects of the ICBI algorithm process that keeps structural texture.

2. a kind of SFS three-dimensional reconstruction of taking face of land spectral information into account according to claim 1 is encrypted sparse DEM method, it is characterized in that, described reflectogram equation is based on formula:

I _(x，y)＝R(p，q)

And hypothesis:

1) topographical surface is lambert's phantom type;

2) solar source can be regarded infinite distant place pointolite parallel projection as;

3) imaging geometry is rectangular projection;

Wherein, original optics image data is designated as I _{(x, y)}, according to light source parameters (Sx, Sy Sz) are designated as R with the reflectogram that reflectivity albedo recovers to obtain, and p, q are surface direction, be z at x, the derivative of y direction:

Lambert's body reflection model R to embody formula following:

R(p，q)＝albedo*(Sz-p*Sx-q*Sy)/sqrt(p ²+q ²+1)

If surface point (x, y, z _{(x, y)}) normal located is N=(cos α sin β, sin α sin β, cos β), light source direction is L=(cos τ sin γ, sin τ sin γ, cos γ), and τ is the light source drift angle, and γ is the light source inclination angle, and α is the normal drift angle, and β is the normal inclination angle, and following formula is:

$R(p,q)=\mathrm{\&eta;}\frac{\cos \mathrm{\&gamma;}-p\cos \mathrm{\&tau;}\sin \mathrm{\&gamma;}-q\sin \mathrm{\&tau;}\sin \mathrm{\&gamma;}}{\sqrt{1+p^2+q^2}}$

$=\mathrm{\&eta;}(\cos (\mathrm{\&alpha;}-\mathrm{\&tau;})\sin \mathrm{\&beta;}\sin \mathrm{\&gamma;}+\cos \mathrm{\&beta;}\cos \mathrm{\&gamma;}).$

3. a kind of SFS three-dimensional reconstruction of taking face of land spectral information into account according to claim 1 is encrypted sparse DEM method, it is characterized in that, in the described step 2, the reflectogram model that step 1 is set up carries out pre-service and may further comprise the steps:

Step 3.1 is carried out cutting and registration to multispectral remote sensing image and dem data, makes the corresponding location matches of shadow information;

Step 3.2, inspection multispectral remote sensing image spectral brightness DN value is zero zone, because the graticule mesh of zero luminance lacks second-order properties, is regarded as shade, carries out this area image reparation in step 3.4;

Step 3.3, the atural object segmentation and classification that the remote sensing image of completing steps 3.2 is carried out based on image spectrum and texture; For cutting apart back image, the surface reflectivity of estimation light source direction and each object;

Step 3.4 utilizes the constant theoretical model of color to carry out image repair, the filtering shadow region to the remote sensing image of completing steps 3.3;

Step 3.5, the remote sensing image of completing steps 3.4 is carried out under the anisotropy condition light and shade recover shape: add the shade of gray constraint in the object, the slickness constraint retrains with integrability, carry out the variational method and fast the method for advancing resolve, obtain face of land relative altitude.

4. a kind of SFS three-dimensional reconstruction of taking face of land spectral information into account according to claim 1 is encrypted sparse DEM method, it is characterized in that, in the described step 3.1, the concrete steps of remote sensing image and dem data being carried out cutting and registration are following:

1) reads dem data and remote sensing image data file header, obtain the coverage of its data;

2) set up the affined transformation model

$\left(\begin{array}{c}{x}^{\&prime;}\\ y^{\&prime;}\end{array}\right)=\left(\begin{array}{cc}{a}_{11}& a_{12}\\ a_{21}& a_{\mathrm{xx}}\end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right)+\left(\begin{array}{c}{t}_x\\ t_y\end{array}\right);$

3) according to the image value range of four angle points up and down, bring the affined transformation model into, obtain eight equations, find the solution 6 transformation model parameters;

4) geographic coordinate after the conversion of calculating every bit.

5. a kind of SFS three-dimensional reconstruction of taking face of land spectral information into account according to claim 1 is encrypted sparse DEM method, it is characterized in that in the described step 3.3, the concrete steps of atural object segmentation and classification are following:

1) combines the atural object texture features according to spectral information, definition classification number c, cost function threshold epsilon, maximum iteration time t; Usually there are 6 big classification atural objects in the remote sensing image, comprise water body, vegetation, the farmland, buildings, road, bare area, the example definition classification number that this patent provides in this step is 3, and the cost function threshold value is 0.1, and maximum iteration time is 50;

2) be subordinate to matrix U with the random number initialization of value between 0-1, make it satisfy the constraint condition in the following formula:

$\underset{i=1}{\overset{c}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ij}}=1,\&ForAll;j=1,...,n$

3) with computes c cluster centre ci, i=1 ..., c:

$c_i=\frac{\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ij}}^m{x}_j}{\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ij}}^m}$

Wherein, x _i(i=1,2 ..., n) be n sample point spectrum vector;

4) according to the computes cost function: if cost function less than the threshold epsilon of confirming, or the change amount of cost function relative last time of cost function value is less than threshold value, then algorithm stops;

$J(U,c_1,...,c_c)=\underset{i=1}{\overset{c}{\mathrm{\&Sigma;}}}{J}_i=\underset{i=1}{\overset{c}{\mathrm{\&Sigma;}}}\underset{j}{\overset{n}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ij}}^m{d}_{\mathrm{ij}}^2$

Here u _IjBetween 0 to 1; c _iBe the cluster centre of ambiguity group I, d _Ij=|| c _i-x _j|| be the Euclidean distance between i cluster centre and j data points; And m ∈ [1, be a weighted index ∞), the example definition m value that this patent provides in this step is 2;

If cost function is greater than threshold value, then with the new U matrix of computes; Return step 2);

$u_{\mathrm{ij}}=\frac{1}{\underset{k=1}{\overset{c}{\mathrm{\&Sigma;}}}{\left(\frac{d_{\mathrm{ij}}}{d_{\mathrm{kj}}}\right)}^{2/(m-1)}};$

5) if iterations surpasses threshold value t, then calculating stops, and finally obtains the terrain classification figure of multispectral remote sensing image.

6. a kind of SFS three-dimensional reconstruction of taking face of land spectral information into account according to claim 1 is encrypted sparse DEM method, it is characterized in that, in the described step 3.3, to cutting apart the back image, the concrete steps of estimation light source direction are following:

1) utilize the linear function conversion that the brightness of image value is carried out normalization:

I _i′＝I _i/I _max

I wherein _iBe a certain spectral band original brightness value at i point place, I _MaxBe this wave band overall situation image brilliance maximal value, I _i' be the brightness value after the normalization;

2) get the original DEM graticule mesh height value of same classification atural object and the brightness value on the correspondence remote sensing image of position thereof,, set up a plurality of reflectogram equations as the reference mark:

I _i＝albedo*(Sz-p _i*Sx-q _i*Sy)/sqrt(p _i ²+q _i ²+1)；

Wherein i ∈ N is a DEM grid points number, I _iBe the brightness value of this point, (p q) is the x of this point, the gradient of y direction, and for same type of atural object, reflectivity albedo is a constant; Promptly utilize N equation (N＞3) resolve 3 light source parameters (Sx, Sy, Sz); Utilize least square method to obtain estimates of parameters.

7. a kind of SFS three-dimensional reconstruction of taking face of land spectral information into account according to claim 1 is encrypted sparse DEM method; It is characterized in that; In the described step 3.3, the concrete steps of surface reflectivity of obtaining each object are following: wherein, each object is meant different atural object classifications:

$\mathrm{\&tau;}=\arctan \left(\frac{E\left\{\frac{y_L}{\sqrt{x_L^2+y_L^2}}\right\}}{E\left\{\frac{x_L}{\sqrt{x_L^2+y_L^2}}\right\}}\right);$

Wherein E is a population mean, $\left(\begin{array}{c}{x}_L\\ y_L\end{array}\right)={\left(B^{\&prime;}B\right)}^{-1}{B}^{\&prime;}\mathrm{DI},\mathrm{DI}=\left(\begin{array}{c}\mathrm{\&delta;}{I}_1\\ \mathrm{\&delta;}{I}_2\\ M\\ \mathrm{\&delta;}{I}_N\end{array}\right),B=\left(\begin{array}{cc}\mathrm{\&delta;}{x}_1& \mathrm{\&delta;}{y}_1\\ \mathrm{\&delta;}{x}_2& \mathrm{\&delta;}{y}_2\\ M& M\\ \mathrm{\&delta;}{x}_N& \mathrm{\&delta;}{y}_N\end{array}\right)$

$B^{\&prime;}=\left(\begin{array}{cccccccc}1& \frac{\sqrt{2}}{2}& 0& -\frac{\sqrt{2}}{2}& -1& -\frac{\sqrt{2}}{2}& 0& \frac{\sqrt{2}}{2}\\ 0& \frac{\sqrt{2}}{2}& 1& \frac{\sqrt{2}}{2}& 0& -\frac{\sqrt{2}}{2}& -1& -\frac{\sqrt{2}}{2}\end{array}\right);$

Order:

$f_1\left(\mathrm{\&gamma;}\right)\&ap;\underset{i=0}{\overset{7}{\mathrm{\&Sigma;}}}{a}_i{\cos }^i\mathrm{\&gamma;}$

$f_2\left(\mathrm{\&gamma;}\right)\&ap;\underset{i=0}{\overset{7}{\mathrm{\&Sigma;}}}{b}_i{\cos }^i\mathrm{\&gamma;},$

$f_3\left(\mathrm{\&gamma;}\right)\&ap;\underset{i=0}{\overset{7}{\mathrm{\&Sigma;}}}{c}_i{\cos }^i\mathrm{\&gamma;}$

If $\frac{E\left\{I\right\}}{\sqrt{E\left\{I^2\right\}}}>f_3\left(0\right)=0.9691,$ γ=0 then; Otherwise γ does $\frac{E\left\{I\right\}}{\sqrt{E\left\{I^2\right\}}}=f_3\left(\mathrm{\&gamma;}\right)$ Separate;

For reflectivity: $\mathrm{\&eta;}=\frac{1}{f_1^2\left(\mathrm{\&gamma;}\right)+f_2\left(\mathrm{\&gamma;}\right)}\&CenterDot;\left\{E\right\{I\}\&CenterDot;f_1\left(\mathrm{\&gamma;}\right)+\sqrt{E\left\{I^2\right\}\&CenterDot;f_2\left(\mathrm{\&gamma;}\right)}\},$ Wherein, I is the normalization brightness value of same type atural object on the image.

8. kind according to claim 1 is taken the SFS three-dimensional reconstruction of face of land spectral information into account and is encrypted sparse DEM method, it is characterized in that, in the described step 3.4, carry out image repair, the concrete steps of filtering shadow region are following:

1) carrying out the shadow region according to each wave band brightness DN value of multispectral remote sensing image and structural information detects: think that it is the shadow region that reflectivity on all wave bands all approaches 0 zone;

2) read in original image I and detected shadow region template, according to template with the original image separated into two parts: the shadow region image I _sWith the non-hatched area image I _Ns

3) respectively to the shadow region image I _sWith the non-hatched area image I _NsCarry out color constancy with Shade of Gray algorithm and calculate, obtain shadow region illumination color and estimate [s_r, s_g, s_b] and nonshaded area illumination color estimation [n_r, n_g, n_b], computing formula is:

$e=k{\left(\frac{\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\left(f{(i,j)}^p\right)}{M*N}\right)}^{1/p}$

In the formula, the light source colour in

expression present image, f={R; G; The brightness DN value of each wave band of B} presentation video, hypothesis has red, green here; Blue three wave bands, p=6;

The correcting gain of each wave band is:

4) shadow region transform to non-hatched area the ratio constant vector is

5) ratio constant vector k is transformed into number space

6) the shadow region image I _sBe transformed into the logarithm Space L _s=log (I _s);

7) obtain removing image behind the shade, I=exp (L _s+ c)+I _Ns

9. kind according to claim 1 is taken the SFS three-dimensional reconstruction of face of land spectral information into account and is encrypted sparse DEM method, it is characterized in that, the concrete steps of carrying out light and shade recovery shape are following:

1) SFS reflectogram equation is rewritten as the function of face of land gradient, and adds the brightness constraint, conditions such as slickness constraint are come solving equation; With the global error structure cost functional of image and reflectogram, and use the variational method to find the solution the three-dimensional surface function;

Possibly there is error between the determined body surface brightness of image of known image gray scale and reflection model, luminance equation is become the form of error function, be i.e. brightness constraint:

$I_1=\underset{\mathrm{\&Omega;}}{\&Integral;\&Integral;}{(E(x,y)-R(p(x,y),q(x,y)))}^2\mathrm{dxdy};$

2) introduce shade of gray constraint and integrability constraint, set up error equation:

∫∫F(p，q，z)dxdy＝∫∫((I-R) ²+(R _x-I _x) ²+(R _y-I _y) ²+μ((z _x-p) ²+(z _y-q) ²))dxdy；

3) utilize Eulerian equation to ask the functional minimal value, through successive iteration, (p is q) with elevation z to obtain the true normal vector in the face of land.

10. kind according to claim 1 is taken the SFS three-dimensional reconstruction of face of land spectral information into account and is encrypted sparse DEM method, it is characterized in that, in the described step 5, carries out following two steps:

Step 5.1 for same atural object object, has promptly kept the face of land changes in the object continuity and slickness, utilizes the result of remote sensing image three-dimensional reconstruction, obtains surface slope aspect figure according to formula:

$\mathrm{slope}=\arctan \sqrt{{f_x}^2+{f_y}^2}*180/\mathrm{\&pi;}$

$\mathrm{aspect}=\arctan \left(\frac{f_x}{f_y}\right)$

Wherein, (f _x, f _y) be the relative height of three-dimensional reconstruction at x, the Grad on the y direction;

Set up cost function, insert the difference minimum of the surface slope aspect figure of back elevation graticule mesh surface slope aspect and three-dimensional reconstruction in making, with moving least squares method iterative calculation DEM graticule mesh;

Step 5.2 is if adjacent twice iteration result's difference less than threshold value 0.1, then with inserting in the edge fit between the different atural object objects of ICBI algorithm process, produces high resolving power DEM graticule mesh, that is:

If the even number site in the graticule mesh is the elevation known point, for the position be (2i+1,2j+1) point at first calculate the value on the both direction:

I ₁₁(2i+1，2j+1)＝I(2i-2，2j+2)+I(2i，2j)+I(2i+2，2j-2)

-3I(2i，2j+2)-3I(2i+2，2j)+I(2i，2j+4)+I(2i+2，2j+2)+I(2i+4，2j)；

I ₂₂(2i+1，2j+1)＝I(2i，2j-2)+I(2i+2，2j)+I(2i+4，2j+2)

-3I(2i，2j)-3I(2i+2，2j+2)+I(2i-2，2j)+I(2i，2j+2)+I(2i+2，2j+4)；

Final its interpolate value is pressed following formula and is confirmed:

Images