CN104036277A - Method and equipment for extracting road characteristics - Google Patents

Abstract

The embodiment of the invention discloses a method and equipment for extracting road characteristics. The method can comprise the steps that: an original SAR (Synthetic Aperture Radar) image is preprocessed through nonlinear quantization and an image filtering method based on an MAR (Multi-scale Auto Regressive) model, and an image after the SAR image preprocessing is obtained; the preprocessed image is subjected to road characteristic extraction by a road extraction algorithm based on the ratio and the direction fusion, and the road characteristics of the original SAR image are obtained; and thus the road characteristics in the SAR image can be efficiently and accurately extracted.

Description

A kind of method and apparatus that extracts roadway characteristic

Technical field

The present invention relates to image processing techniques, relate in particular to a kind of method and apparatus that extracts roadway characteristic.

Background technology

Current, road extraction in Technologies Against Synthetic Aperture Radar (SAR, Synthetic Aperture Radar) image can comprise three levels conventionally: the road feature extraction of low level, middle-level road primitives extract and high-level road network is set up.And the effect of road feature extraction in road extraction is most important, determining performance and the efficiency of whole road extraction algorithm.

Summary of the invention

The embodiment of the present invention is expected to provide a kind of method and apparatus that extracts roadway characteristic, can efficiently and exactly extract the roadway characteristic in SAR image.

For achieving the above object, technical scheme of the present invention is achieved in that

First aspect, the embodiment of the present invention provides a kind of method of extracting roadway characteristic, and described method comprises:

By nonlinear quantization and the image filtering method based on Multiscale Autoregressive (MAR, Multiscale Auto Regressive) model, original SAR image is carried out to pre-service, obtain the pretreated image of described SAR image;

By the road extraction algorithm based on ratio and direction fusion, described pretreated image is carried out to road feature extraction, obtain the roadway characteristic of described original SAR image.

The implementation possible according to the first, in conjunction with first aspect, describedly carries out pre-service by nonlinear quantization and the image filtering method based on MAR model to original SAR image, comprising:

Remove a small amount of supersaturation point of described original SAR image by the grey level histogram of described original SAR image, and obtain the dynamic range of the gray scale of described original SAR image;

By Nonlinear extension algorithm by the dynamic range expansion to 0 of the gray scale of described original SAR image between 255, obtain the SAR image after nonlinear quantization.

The implementation possible according to the second, in conjunction with the possible implementation of the first, describedly carries out pre-service by nonlinear quantization and the image filtering method based on MAR model to original SAR image, also comprises:

First, the SAR image after described nonlinear quantization is carried out to pyramid decomposition, obtain the SAR image sequence Φ that resolution reduces successively _p=[x ₀x ₁... x _p];

Wherein, p+1 is described SAR image sequence Φ _pin SAR image number, x ₀for described SAR image sequence Φ _pthe SAR image that intermediate-resolution is the highest, x _pfor described SAR image sequence Φ _pthe SAR image that intermediate-resolution is minimum; x _lfor x _l-1obtain according to formula (1):

$x_l(m,n)=\underset{i=2m}{\overset{2m+1}{\mathrm{\Σ}}}\underset{j=2n}{\overset{2n+1}{\mathrm{\Σ}}}{x}_{l-1}(i,j)---\left(1\right)$

Wherein, x _i(m, n) presentation video x _ithe pixel value of the capable n row of m, x _l-1(i, j) presentation video x _l-1the pixel value of the capable j of i row, l is greater than zero and be less than the natural number of p;

Then, according to described SAR image sequence Φ _pbuild MAR model; Wherein, described MAR model is suc as formula shown in (2):

$x_0=\underset{i=1}{\overset{p}{\mathrm{\Σ}}}{a}_i{x}_i+\mathrm{\ϵ}---\left(2\right)$

Wherein, a=[a ₁a ₂... a _p] be the autoregressive coefficient of described MAR model, ε is error term, p is the order of described MAR model;

Finally, according to default estimation criterion, by least square method, the autoregressive coefficient a of described MAR model is estimated, obtain described SAR image sequence Φ _pthe SAR image x that intermediate-resolution is the highest ₀predicted picture wherein, predicted picture for the pretreated image of described SAR image, shown in (3):

${\hat{x}}_0=\underset{i=1}{\overset{p}{\mathrm{\Σ}}}{a}_i{x}_i---\left(3\right).$

The implementation possible according to the third, in conjunction with first aspect, describedly carries out road feature extraction by the road extraction algorithm merging based on ratio and direction to described pretreated image, comprising:

Described pretreated image is carried out to piecemeal, obtain at least one image block;

By the road extraction algorithm based on ratio and direction fusion, the roadway characteristic of each image block is extracted;

Roadway characteristic after the extraction of each image block is spliced, obtain the roadway characteristic of described original SAR image.

According to the 4th kind of possible implementation, in conjunction with the third possible implementation, describedly by the road extraction algorithm merging based on ratio and direction, the roadway characteristic of each image block is extracted, comprising:

The detection of sliding on described image block by rotatable ROA operator, obtains rate information R and the directional information θ of described image block; Wherein, described rotatable ROA operator comprise a main window W1 and two centrosymmetric from window W2 and W3 with main window;

Convert according to the rate information R of described image block, directional information θ and La Dong Radon the road principal direction of obtaining described image block and according to the road principal direction of described image block the best of intersecting roads in described image block is decomposed, obtain two exploded views of described image block;

The Road Detection operator extraction that two of described image block exploded views are all used to ROA operator and simple crosscorrelation Operator Fusion road information separately, and the road information of two exploded views of described image block is merged, thereby obtain the roadway characteristic of described image block.

According to the 5th kind of possible implementation, in conjunction with the 4th kind of possible implementation, described by image block is slided on the described image block detection of rotatable ROA operator, obtain rate information R and the directional information θ of described image block, comprising:

Calculate described rotatable ROA operator current location (x ₀, y ₀) maximum rate value R (x ₀, y ₀) and the maximum rate value R (x of described current location ₀, y ₀) corresponding direction θ (x ₀, y ₀);

Described rotatable ROA operator is slided into next position (x ₁, y ₁) calculate the maximum rate value R (x of next position described in described rotatable ROA operator ₁, y ₁) and the maximum rate value R (x of described next position ₁, y ₁) corresponding direction θ (x ₁, y ₁), until described rotatable ROA operator slides complete by described image block;

Wherein, calculating the direction θ that the maximum rate value R of described rotatable ROA operator current location and the maximum rate value R of described current location are corresponding is obtained by formula (4) to formula (6):

$r(x,y,{\mathrm{\θ}}_k)=\min (1-\frac{M_1(x,y,{\mathrm{\θ}}_k)}{M_2(x,y,{\mathrm{\θ}}_k)},1-\frac{M_1(x,y,{\mathrm{\θ}}_k)}{M_3(x,y,{\mathrm{\θ}}_k)})---\left(4\right)$

$R(x,y)=\underset{k=1:N}{\max }\left(r(x,y,{\mathrm{\θ}}_k)\right)---\left(5\right)$

$\mathrm{\θ}(x,y)=\underset{{\mathrm{\θ}}_k,k=1:N}{\arg }\max \left\{r(x,y,{\mathrm{\θ}}_k)\right\}---\left(6\right)$

Wherein, r (x, y, θ _k) be that the described rotatable ROA operator of current location is at direction θ _kratio, M ₁(x, y, θ _k), M ₂(x, y, θ _k) and M ₃(x, y, θ _k) represent that respectively window W1, the W2 of described rotatable ROA operator and W3 are at direction θ _kpixel value average; N represents the total degree that described rotatable ROA operator is rotated; represent that current described rotatable ROA operator is in directive ratio maximal value; $\underset{{\mathrm{\θ}}_k,k=1:N}{\arg }\max \left\{r(x,y,{\mathrm{\θ}}_k)\right\}$ Represent $\underset{k=1:N}{\max }\left(r(x,y,{\mathrm{\θ}}_k)\right)$ Corresponding direction.

According to the 6th kind of possible implementation, in conjunction with the 4th kind of possible implementation, convert according to the rate information R of described image block, directional information θ and Radon the road principal direction of obtaining described image block and according to the road principal direction of described image block the best of intersecting roads in described image block is decomposed, obtains two exploded views of described image block, comprising:

The rate information R of the each pixel of described image block and directional information θ are transformed to the direction vector (p under Descartes's rectangular coordinate that described each pixel is corresponding according to formula (7) _x, p _y);

$\left\{\begin{array}{c}{p}_x=R\cos \mathrm{\θ}\\ p_y=R\sin \mathrm{\θ}\end{array}\right.---\left(7\right)$

Obtain the local peaking of formula (8) in described image block according to the rate information R of described image block;

Wherein, D is integral domain, represents described image block plane in the present embodiment; R (x, y) is the value that the rate information R of described image block locates at (x, y); δ is Dirac function; ρ be the initial point of Descartes's right-angle plane to the distance of straight line, for the initial point of described Descartes's right-angle plane is to the angle of the vertical line of straight line and the x axle of described Descartes's right-angle plane;

Described in inciting somebody to action what direction was corresponding owns local peaking superpose and obtain described image block and exist there is the probability level of straight line in direction

Searching probability index maximal value, and by described probability level the corresponding angle of maximal value as the road principal direction of described image block, wherein, described probability level the corresponding angle of maximal value through type (9) represents:

By the direction vector (p under corresponding Descartes's rectangular coordinate of the each pixel of described image block _x, p _y) the extremely road principal direction of described image block of through type (10) rotation the best that obtains the intersecting roads to described image block is decomposed, and wherein, formula (10) is:

Wherein, (p _x, p _y) be the direction vector under Descartes's rectangular coordinate before the each pixel rotation of described image block, (p _x', p _y') be the direction vector under Descartes's rectangular coordinate after the each pixel rotation of described image block.

According to the 7th kind of possible implementation, in conjunction with the 4th kind of possible implementation, describedly all use ROA operator and simple crosscorrelation operator to detect road information separately to two of described image block exploded views, and the road information of two exploded views of described image block is merged to the roadway characteristic that obtains described image block, comprising:

Two of described image block exploded views are detected by ROA operator respectively, obtain the rate information R ' of described two exploded views;

Two of described image block exploded views are detected by simple crosscorrelation operator respectively, obtain the cross correlation value ξ corresponding to two exploded views of described image block;

The cross correlation value ξ that the rate information R ' of two exploded views of described image block is corresponding with two exploded views of described image block merges, and by the result passing threshold processing after merging, obtains two roadway characteristics that exploded view is corresponding of described image block;

Two roadway characteristics corresponding to exploded view of described image block are merged to the roadway characteristic that obtains described image block.

Second aspect, the embodiment of the present invention provides a kind of equipment that extracts roadway characteristic, and described equipment comprises: pretreatment unit and road feature extraction unit, wherein,

Described pretreatment unit, for original synthetic-aperture radar SAR image being carried out to pre-service by nonlinear quantization and the image filtering method based on Multiscale Autoregressive MAR model, obtains the pretreated image of described SAR image;

Described road feature extraction unit, for by the road extraction algorithm based on ratio and direction fusion, described pretreated image being carried out to road feature extraction, obtains the roadway characteristic of described original SAR image.

The implementation possible according to the first, in conjunction with second aspect, described pretreatment unit comprises nonlinear quantization subelement, remove a small amount of supersaturation point of described original SAR image for the grey level histogram by described original SAR image, and obtain the dynamic range of the gray scale of described original SAR image;

And by Nonlinear extension algorithm by the dynamic range expansion of the gray scale of described original SAR image between 0-255, obtain the SAR image after nonlinear quantization.

The implementation possible according to the second, in conjunction with the possible implementation of the first, described pretreatment unit also comprises filtering subelement, carries out pyramid decomposition for the SAR image to after described nonlinear quantization, obtains the SAR image sequence Φ that resolution reduces successively _p=[x ₀x ₁... x _p];

Wherein, p+1 is described SAR image sequence Φ _pin SAR image number, x ₀for described SAR image sequence Φ _pthe SAR image that intermediate-resolution is the highest, x _pfor described SAR image sequence Φ _pthe SAR image that intermediate-resolution is minimum; x _lfor x _l-1obtain according to formula (11):

$x_l(m,n)=\underset{i=2m}{\overset{2m+1}{\mathrm{\Σ}}}\underset{j=2n}{\overset{2n+1}{\mathrm{\Σ}}}{x}_{l-1}(i,j)---\left(11\right)$

Wherein, x _i(m, n) presentation video x _ithe pixel value of the capable n row of m, x _l-1(i, j) presentation video x _l-1the pixel value of the capable j of i row, l is greater than zero and be less than the natural number of p;

And according to described SAR image sequence Φ _pbuild MAR model; Wherein, described MAR model is suc as formula shown in (12):

$x_0=\underset{i=1}{\overset{p}{\mathrm{\Σ}}}{a}_i{x}_i+\mathrm{\ϵ}---\left(12\right)$

Wherein, a=[a ₁a ₂... a _p] be the autoregressive coefficient of described MAR model, ε is error term, p is the order of described MAR model;

And according to default estimation criterion, by least square method, the autoregressive coefficient a of described MAR model is estimated, obtain described SAR image sequence Φ _pthe SAR image x that intermediate-resolution is the highest ₀predicted picture described predicted picture wherein, predicted picture for the pretreated image of described SAR image, shown in (13):

${\hat{x}}_0=\underset{i=1}{\overset{p}{\mathrm{\Σ}}}{a}_i{x}_i---\left(13\right).$

The implementation possible according to the third, in conjunction with second aspect, described road feature extraction unit comprises that piecemeal subelement, image block characteristics extract subelement and splicing subelement, wherein,

Described piecemeal subelement, for described pretreated image is carried out to piecemeal, obtains at least one image block;

Described image block characteristics extracts subelement, for the roadway characteristic of each image block being extracted by the road extraction algorithm based on ratio and direction fusion;

Described splicing subelement, for the roadway characteristic after the extraction of each image block is spliced, obtains the roadway characteristic of described original SAR image.

According to the 4th kind of possible implementation, in conjunction with the third possible implementation, described image block characteristics extracts subelement and is used for:

The detection of sliding on described image block by rotatable ROA operator, obtains rate information R and the directional information θ of described image block; Wherein, described rotatable ROA operator comprise a main window W1 and two centrosymmetric from window W2 and W3 with main window;

Convert according to the rate information R of described image block, directional information θ and La Dong Radon the road principal direction of obtaining described image block and according to the road principal direction of described image block the best of intersecting roads in described image block is decomposed, obtain two exploded views of described image block;

All use ROA operator and simple crosscorrelation operator to detect road information separately to two of described image block exploded views, and the road information of two exploded views of described image block is merged to the roadway characteristic that obtains described image block.

According to the 5th kind of possible implementation, in conjunction with the 4th kind of possible implementation, described image block characteristics extracts subelement and is used for:

Calculate described rotatable ROA operator current location (x ₀, y ₀) maximum rate value R (x ₀, y ₀) and the maximum rate value R (x of described current location ₀, y ₀) corresponding direction θ (x ₀, y ₀);

Described rotatable ROA operator is slided into next position (x ₁, y ₁) calculate the maximum rate value R (x of next position described in described rotatable ROA operator ₁, y ₁) and the maximum rate value R (x of described next position ₁, y ₁) corresponding direction θ (x ₁, y ₁), until described rotatable ROA operator slides complete by described image block;

Wherein, calculating the direction θ that the maximum rate value R of described rotatable ROA operator current location and the maximum rate value R of described current location are corresponding is obtained by formula (14) to formula (16):

$r(x,y,{\mathrm{\θ}}_k)=\min (1-\frac{M_1(x,y,{\mathrm{\θ}}_k)}{M_2(x,y,{\mathrm{\θ}}_k)},1-\frac{M_1(x,y,{\mathrm{\θ}}_k)}{M_3(x,y,{\mathrm{\θ}}_k)})---\left(14\right)$

$R(x,y)=\underset{k=1:N}{\max }\left(r(x,y,{\mathrm{\θ}}_k)\right)---\left(15\right)$

$\mathrm{\θ}(x,y)=\underset{{\mathrm{\θ}}_k,k=1:N}{\arg }\max \left\{r(x,y,{\mathrm{\θ}}_k)\right\}---\left(16\right)$

Wherein, r (x, y, θ _k) be that the described rotatable ROA operator of current location is at direction θ _kratio, M ₁(x, y, θ _k), M ₂(x, y, θ _k) and M ₃(x, y, θ _k) represent that respectively window W1, the W2 of described rotatable ROA operator and W3 are at direction θ _kpixel value average; N represents the total degree that described rotatable ROA operator is rotated; represent that current described rotatable ROA operator is in directive ratio maximal value; $\underset{{\mathrm{\θ}}_k,k=1:N}{\arg }\max \left\{r(x,y,{\mathrm{\θ}}_k)\right\}$ Represent $\underset{k=1:N}{\max }\left(r(x,y,{\mathrm{\θ}}_k)\right)$ Corresponding direction.

According to the 6th kind of possible implementation, in conjunction with the 4th kind of possible implementation, described image block characteristics extracts subelement and is used for:

The rate information R of the each pixel of described image block and directional information θ are transformed to the direction vector (p under Descartes's rectangular coordinate that described each pixel is corresponding according to formula (17) _x, p _y);

$\left\{\begin{array}{c}{p}_x=R\cos \mathrm{\θ}\\ p_y=R\sin \mathrm{\θ}\end{array}\right.---\left(17\right)$

And obtain the local peaking of formula (18) in described image block according to the rate information R of described image block;

Wherein, D is integral domain, represents described image block plane in the present embodiment; R (x, y) is the value that the rate information R of described image block locates at (x, y); δ is Dirac function; ρ be the initial point of Descartes's right-angle plane to the distance of straight line, for the initial point of described Descartes's right-angle plane is to the angle of the vertical line of straight line and the x axle of described Descartes's right-angle plane;

And described in inciting somebody to action what direction was corresponding owns local peaking superpose and obtain described image block and exist there is the probability level of straight line in direction

And searching probability index maximal value, and by described probability level the corresponding angle of maximal value as the road principal direction of described image block, wherein, described probability level the corresponding angle of maximal value through type (19) represents:

And by the direction vector (p under corresponding Descartes's rectangular coordinate of the each pixel of described image block _x, p _y) the extremely road principal direction of described image block of through type (20) rotation the best that obtains the intersecting roads to described image block is decomposed, and wherein, formula (10) is:

Wherein, (p _x, p _y) be the direction vector under the Descartes's rectangular coordinate before the each pixel rotation of described image block, (p _x', p _y') be the direction vector under the postrotational Descartes's rectangular coordinate of the each pixel of described image block.

According to the 7th kind of possible implementation, in conjunction with the 4th kind of possible implementation, described image block characteristics extracts subelement and is used for:

Two of described image block exploded views are detected by ROA operator respectively, obtain the rate information R ' of described two exploded views;

Two of described image block exploded views are detected by simple crosscorrelation operator respectively, obtain the cross correlation value ξ corresponding to two exploded views of described image block;

And the rate information R ' of two exploded views of the described image block cross correlation value ξ corresponding with two exploded views of described image block merged, and by the result passing threshold processing after merging, obtain two roadway characteristics that exploded view is corresponding of described image block;

And two roadway characteristics corresponding to exploded view of described image block are merged to the roadway characteristic that obtains described image block.

The embodiment of the present invention provides a kind of method and apparatus that extracts roadway characteristic, by nonlinear quantization and the image filtering method based on Multiscale Autoregressive model, original SAR image is carried out pre-service and by the road extraction algorithm based on ratio and direction fusion, described pretreated image carried out to road feature extraction, can efficiently and exactly extract the roadway characteristic in SAR image.

Brief description of the drawings

The original SAR image that certain radar satellite that Fig. 1 provides for the embodiment of the present invention acquires B region, A city;

A kind of method flow schematic diagram that extracts roadway characteristic that Fig. 2 provides for the embodiment of the present invention;

The original SAR image to shown in Fig. 1 that Fig. 3 provides for the embodiment of the present invention carries out the design sketch after nonlinear quantization;

Fig. 4 carries out filtered design sketch for the image obtaining after to the nonlinear quantization shown in Fig. 3 based on MAR model that the embodiment of the present invention provides;

Fig. 5 is the enlarged drawing of the image block in the white box shown in Fig. 4;

A kind of realization flow schematic diagram roadway characteristic of image block being extracted by the road extraction algorithm based on ratio and direction fusion that Fig. 6 provides for the embodiment of the present invention;

The schematic diagram of a kind of rotatable ROA operator that Fig. 7 provides for the embodiment of the present invention;

Fig. 8 a is the ratio chart of image block shown in Fig. 5;

Fig. 8 b is the directional diagram of image block shown in Fig. 5;

The conversion according to the rate information R of image block, directional information θ and La Dong Radon that Fig. 9 provides for the embodiment of the present invention obtains the realization flow schematic diagram of two exploded views of described image block;

Figure 10 a is the result schematic diagram after changing through Radon of image block shown in Fig. 5;

The probability level that Figure 10 b provides for the embodiment of the present invention and the corresponding relation figure of direction;

Figure 11 a is that image block shown in Fig. 5 is at the exploded view of directions X;

Figure 11 b is that image block shown in Fig. 5 is at the exploded view of Y-direction;

Figure 12 a is the roadway characteristic schematic diagram of directions X exploded view shown in Figure 11 a;

Figure 12 b is the roadway characteristic schematic diagram of Y-direction exploded view shown in Figure 11 b;

Figure 13 is the road feature extraction result schematic diagram of image block shown in Fig. 5;

Figure 14 is the road feature extraction result schematic diagram of the original SAR image shown in Fig. 1;

Figure 15 is the result that shown in Fig. 5, image block carries out manual extraction roadway characteristic;

Figure 16 a only carries out, after autoregression filtering, image block shown in Fig. 5 is carried out to road extraction result schematic diagram to original SAR image shown in Fig. 1 for what the embodiment of the present invention provided;

Figure 16 b for the embodiment of the present invention provide for original SAR image shown in Fig. 1 only being carried out image block shown in Fig. 5 being carried out after nonlinear quantization the result schematic diagram of road extraction;

Figure 16 c does not carry out any pre-service and directly image block shown in Fig. 5 is carried out the result schematic diagram of road extraction for what the embodiment of the present invention provided to original SAR image shown in Fig. 1;

The structural representation of a kind of equipment that extracts roadway characteristic that Figure 17 provides for the embodiment of the present invention;

The another kind that Figure 18 provides for the embodiment of the present invention extracts the structural representation of the equipment of roadway characteristic.

Embodiment

Below in conjunction with the accompanying drawing in the embodiment of the present invention, the technical scheme in the embodiment of the present invention is clearly and completely described.

For the technical scheme of the embodiment of the present invention can be clearly described, the SAR image that the embodiment of the present invention is all acquired B region, A city using certain radar satellite describes as the image object of processing.As shown in Figure 1, the original SAR image that its certain radar satellite that shows the embodiment of the present invention and provide acquires B region, A city.Those skilled in the art are understandable, do not represent that the embodiment of the present invention is only for processing this SAR image.

Referring to Fig. 2, it shows a kind of method flow that extracts roadway characteristic that the embodiment of the present invention provides, and as shown in Figure 2, the method can comprise:

S201: original SAR image is carried out to pre-service by nonlinear quantization and the image filtering method based on Multiscale Autoregressive model;

Particularly, nonlinear quantization process can comprise: first remove a small amount of supersaturation point of described original SAR image by the grey level histogram of original SAR image, and obtain the dynamic range of the gray scale of described original SAR image; Then by Nonlinear extension algorithm, described original SAR image gray levels is expanded between 0-255, thereby by the upwards lifting of gray-scale value of described original SAR dark picture areas; Preferably, described Nonlinear extension algorithm can be selected gamma (Gamma) correcting algorithm, and in the time that described original SAR integral image is partially dark, Gamma value should preferably arrange and be less than 1.

As shown in Figure 3, it shows the original SAR image shown in Fig. 1 is carried out to nonlinear quantization design sketch afterwards, from figure, can find out significantly, compared with the original SAR image shown in Fig. 1, the effect of visualization that original SAR image is carried out to the design sketch after nonlinear quantization shown in Fig. 3 obviously promotes, and the lines that link table reveals are also more obvious than original image.But owing to passing through nonlinear quantization afterwards by integral raising compared with the pixel value of dark areas in original SAR image, therefore, the image coherent noise in the design sketch shown in Fig. 3 can expand thereupon.

The problem expanding in order to solve image coherent noise that nonlinear quantization causes, after original SAR image is carried out to nonlinear quantization, also needs to continue the image obtaining after nonlinear quantization to carry out filtering;

Preferably, the image that the embodiment of the present invention obtains after adopting image filtering method based on Multiscale Autoregressive (MAR, Multiscale Auto Regressive) model to described nonlinear quantization carries out filtering, and concrete process is as follows:

First, the image obtaining after described nonlinear quantization is carried out to pyramid decomposition, obtain the SAR image sequence Φ that resolution reduces successively _p=[x ₀x ₁... x _p], it should be noted that, image is carried out to the common technology means that pyramid decomposition is those skilled in the art, do not repeat them here; Wherein, p+1 is SAR image sequence Φ _pin SAR image number, x ₀for described SAR image sequence Φ _pthe SAR image that intermediate-resolution is the highest, is generally the SAR image that carries out pyramid decomposition, in the present embodiment, and x ₀for the image obtaining after described nonlinear quantization; x _pfor described SAR image sequence Φ _pthe SAR image that intermediate-resolution is minimum; x _lfor x _l-1obtain according to formula (1):

$x_l(m,n)=\underset{i=2m}{\overset{2m+1}{\mathrm{\Σ}}}\underset{j=2n}{\overset{2n+1}{\mathrm{\Σ}}}{x}_{l-1}(i,j)---\left(1\right)$

Wherein, x _i(m, n) presentation video x _ithe pixel value of the capable n row of m, x _l-1(i, j) presentation video x _l-1the pixel value of the capable j of i row, l is greater than zero and be less than the natural number of p.

Then, according to described SAR image sequence Φ _pbuild MAR model, described MAR model is suc as formula shown in (2):

$x_0=\underset{i=1}{\overset{p}{\mathrm{\Σ}}}{a}_i{x}_i+\mathrm{\ϵ}---\left(2\right)$

Wherein, a=[a ₁a ₂... a _p] be the autoregressive coefficient of described MAR model, ε is error term, and p is the order of described MAR model, and in the present embodiment, the p in formula (2) equals the p in formula (1).

Finally, according to default estimation criterion, by least square method, the autoregressive coefficient a of described MAR model is estimated, obtain described SAR image sequence Φ _pthe SAR image x that intermediate-resolution is the highest ₀predicted picture shown in (3), predicted picture for:

${\hat{x}}_0=\underset{i=1}{\overset{p}{\mathrm{\Σ}}}{a}_i{x}_i---\left(3\right)$

Preferably, in the present embodiment, described default estimation criterion can be minimum mean square error criterion, meets a be the autoregressive coefficient of described MAR model, wherein, min (*) represents the minimum value of *.

Obtaining predicted picture afterwards, predicted picture the image obtaining after nonlinear quantization is carried out to the image obtaining after filtering; As shown in Figure 4, it shows the image obtaining after to the nonlinear quantization shown in Fig. 3 based on MAR model and carries out filtered design sketch, and Fig. 4 and Fig. 3 are contrasted and can be found out: the flatness of the SAR image obtaining afterwards after filtering improves, road area is more obvious.

Its reason is: first, and the SAR image sequence Φ obtaining by pyramid decomposition _pin SAR image there is good flatness, therefore, obtain described SAR image sequence Φ by MAR model _pthe SAR image x that intermediate-resolution is the highest ₀predicted picture equally also there is good flatness; Secondly, because MAR model is to be based upon SAR image sequence Φ _pon correlativity basis between the pixel of middle two adjacent images, therefore, predicted picture the hold facility of details better.

So far, completed the preprocessing process to original SAR image.

S202: described pretreated image is carried out to road feature extraction by the road extraction algorithm based on ratio and direction fusion;

It should be noted that, because the trend of road in SAR image is normally inconsistent.Therefore, preferably, can first pretreated image be carried out to piecemeal, obtain at least one image block; Then by the road extraction algorithm based on ratio and direction fusion, the roadway characteristic of each image block is extracted; Then the roadway characteristic after the extraction of each image block is spliced, thereby obtain the roadway characteristic of entire image.

For the technical scheme of the embodiment of the present invention can be clearly described, the present embodiment describes as an example of the image block in the white box shown in Fig. 4 example, referring to Fig. 5, it shows the enlarged image of image block in the white box shown in Fig. 4, understandable, other image blocks that those skilled in the art can be applied to S202 in image carry out road feature extraction.

Exemplary, referring to Fig. 6, it shows a kind of realization flow roadway characteristic of image block being extracted by the road extraction algorithm based on ratio and direction fusion that the embodiment of the present invention provides, and as shown in Figure 6, this realization flow can comprise:

S2021: by the detection of sliding of rotatable ROA (Ratio of Average) operator, obtain rate information R and the directional information θ of described image block on described image block;

Understandably, because the coherent noise in SAR image is generally multiplicative noise, cause being applicable to the detection operator of additive noise and being not suitable for the processing to SAR image, and ROA operator can reduce the influence of fluctuations that multiplicative noise causes single pixel, thereby the present embodiment is preferably ROA operator.

Concrete, rotatable ROA operator as shown in Figure 7, can comprise a main window W1 and two centrosymmetric from window W2 and W3 with main window, wherein, main window W1 fills and represents with grid, W2 fills and represents with horizontal stripe, W3 represents with nicking filling; It should be noted that, the size of these three windows is provided with certain criterion, and in simple terms, these criterions can comprise: when the width of the window of ROA operator arranges when excessive, there will be details fuzzy and not obvious to thin road response; When the window width of ROA operator arranges when too small, there will be false edge and the repeatedly response to a road; When the length of window of ROA operator arranges when too short, can cause the road continuity extracted poor; When the length of window of ROA operator arranges when long, can cause the response on short-track road not obvious.In the present embodiment, two width from window W2 and W3 are all set to 8 pixels, and two length from window W2 and W3 are all set to and 51 pixels, and the width of main window W1 is set to 1 pixel, and the length of main window W1 is set to 31 pixels.Also it should be noted that, rotatable ROA operator shown in Fig. 7 can also be rotated, the angle of each rotation can be determined according to concrete needs or these default pivot rules of empirical data, in the present embodiment, described default pivot rule can be set to described rotatable ROA operator and turn clockwise 5 ° at every turn.

Window size and pivot rule for ROA operator can arrange in advance, need in the process that realizes S2021, just not arrange.

Concrete, the described process detecting of image block being slided on described image block by rotatable ROA operator can be: calculate described rotatable ROA operator current location (x ₀, y ₀) maximum rate value R (x ₀, y ₀) and the maximum rate value R (x of described current location ₀, y ₀) corresponding direction θ (x ₀, y ₀); Described rotatable ROA operator is slided into next position (x ₁, y ₁) calculate the maximum rate value R (x of next position described in described rotatable ROA operator ₁, y ₁) and the maximum rate value R (x of described next position ₁, y ₁) corresponding direction θ (x ₁, y ₁), until described rotatable ROA operator slides complete by described image block.

More detailed, calculate the maximum rate value R of described rotatable ROA operator current location ₁and the maximum rate value R of described current location ₁corresponding direction θ ₁can obtain according to formula (4) to formula (6):

$r(x,y,{\mathrm{\θ}}_k)=\min (1-\frac{M_1(x,y,{\mathrm{\θ}}_k)}{M_2(x,y,{\mathrm{\θ}}_k)},1-\frac{M_1(x,y,{\mathrm{\θ}}_k)}{M_3(x,y,{\mathrm{\θ}}_k)})---\left(4\right)$

$R(x,y)=\underset{k=1:N}{\max }\left(r(x,y,{\mathrm{\θ}}_k)\right)---\left(5\right)$

$\mathrm{\θ}(x,y)=\underset{{\mathrm{\θ}}_k,k=1:N}{\arg }\max \left\{r(x,y,{\mathrm{\θ}}_k)\right\}---\left(6\right)$

Wherein, r (x, y, θ _k) be that current described rotatable ROA operator is at direction θ _kratio, M ₁(x, y, θ _k), M ₂(x, y, θ _k) and M ₃(x, y, θ _k) represent that respectively window W1, the W2 of described rotatable ROA operator and W3 are at direction θ _kpixel value average; N represents the total degree that described rotatable ROA operator is rotated; represent that current described rotatable ROA operator is in directive ratio maximal value; $\underset{{\mathrm{\θ}}_k,k=1:N}{\arg }\max \left\{r(x,y,{\mathrm{\θ}}_k)\right\}$ Represent $\underset{k=1:N}{\max }\left(r(x,y,{\mathrm{\θ}}_k)\right)$ Corresponding direction;

Understandable, in the time that rotatable ROA operator slides into a position, all can obtain maximum rate value and direction corresponding to this maximum rate value of rotatable ROA operator in this position according to formula (4) to formula (6); In the present embodiment, as shown in Fig. 8 a and Fig. 8 b, the result that can obtain sliding on the image block shown in Fig. 5 through rotatable ROA operator and detect, wherein, Fig. 8 a is the ratio chart of image block shown in Fig. 5, for the rate information of image block shown in presentation graphs 5; Fig. 8 b is the directional diagram of image block shown in Fig. 5, for the directional information of image block shown in presentation graphs 5.

Through type (4) not only can obtain the rate information R of image block to formula (6), can also obtain directional information θ, and θ has certain robustness.This is because with respect to R, and θ has certain trend of road instruction ability, such as, the θ of non-roadway area is disorderly and unsystematic, and the θ of roadway area has coherence and continuity; And in some road fragment, for example, because the interference (: greenbelt, building, water body etc.) of roadway area can cause R less, but the deviation of θ is not too large.Therefore θ can be used as the important supplement of R, thereby improves the operational performance of ROA operator.

S2022: convert the road principal direction of obtaining described image block according to the rate information R of described image block, directional information θ and La Dong (Radon) and according to the road principal direction of described image block the best of intersecting roads in described image block is decomposed, obtain two exploded views of described image block;

Concrete, as shown in Figure 9, the specific implementation flow process of S2022 can comprise:

S20221: the rate information R of the each pixel of described image block and directional information θ are transformed to the direction vector (p under Descartes's rectangular coordinate that described each pixel is corresponding according to formula (7) _x, p _y);

$\left\{\begin{array}{c}{p}_x=R\cos \mathrm{\θ}\\ p_y=R\sin \mathrm{\θ}\end{array}\right.---\left(7\right)$

It should be noted that, through the variation of formula (6), can make the coordinate of directions X under Descartes's rectangular coordinate and Y-direction be the function of R and θ, so not only 0 ° and 90 ° of cross one another straight lines can be separated in different planes, thereby make the lineal energy on both direction more concentrated; But also remove vertical information interference.

S20222: obtain the local peaking of formula (8) in described image block according to the rate information R of described image block;

Wherein, D is integral domain, represents whole described image block plane in the present embodiment; R (x, y) is the value that the rate information R of described image block locates at (x, y); δ is Dirac function; ρ be the initial point of Descartes's right-angle plane to the distance of straight line, for the initial point of described Descartes's right-angle plane is to the angle of the vertical line of straight line and the x axle of described Descartes's right-angle plane;

It should be noted that, formula (8) is the formula of Radon conversion, and in the present embodiment, the result after described image block changes through Radon can be as shown in Figure 10 a, can find in Figure 10 a the local peaking of angle corresponding to white box region.

S20223: described in inciting somebody to action what direction was corresponding owns local peaking superpose and obtain described image block and exist there is the probability level of straight line in direction

S20224: searching probability index maximal value, and by described probability level the corresponding angle of maximal value as the road principal direction of described image block;

At length, S20224 can represent by through type (9):

Concrete, in the present embodiment, as shown in Figure 10 b, for there is the probability level of straight line in direction, can find out probability level from Figure 10 b maximal value appears at 67 ° and locates, and shows that the road principal direction of described image block and y axle clamp angle are 67 °.

S20225: by the direction vector (p under corresponding Descartes's rectangular coordinate of the each pixel of described image block _x, p _y) the extremely road principal direction of described image block of rotation decompose with the best that realizes the intersecting roads to described image block;

Detailed rotary course can obtain by through type (10):

Wherein, (p _x, p _y) be the direction vector under the Descartes's rectangular coordinate before rotation, (p _x', p _y') be the direction vector under postrotational Descartes's rectangular coordinate.

After the intersecting roads of described image block being decomposed by S20221 to S20225, can obtain two exploded views of described image block.These two exploded views are respectively as shown in Figure 11 a and Figure 11 b, and wherein, Figure 11 a is the exploded view of described image block at directions X, and Figure 11 b is the exploded view of described image block in Y-direction.From Figure 11 a with Figure 11 b, can find out convert by Radon after the road vectors of principal direction of road obtained separating to the full extent with y directional dependency at x, be conducive to the testing process for the second time described in ensuing S2023.

S2023: all use ROA operator and simple crosscorrelation operator to detect road information separately to two of described image block exploded views, and the road information of two exploded views of described image block is merged to the roadway characteristic that can obtain described image block.

Concrete, first two of described image block exploded views are detected by ROA operator respectively, obtain the rate information R ' of two exploded views of described image block; Detailed testing process, as described in S2021, does not repeat them here, and it should be noted that, the bright dark-part of exploded view is contrary with the bright dark-part of former figure, therefore need to change the min in formula (4) into max;

Then, two of described image block exploded views are detected by simple crosscorrelation operator respectively, to the cross correlation value ξ corresponding to two exploded views of described image block; Detailed testing process can be: in two exploded views of described image block, slide respectively by described ROA operator, and calculate described ROA operator at the cross correlation value ξ=min of the main window of current sliding position (ξ in the process of sliding _{w1, W2}, ξ _{w1, W3})

Wherein, ξ _{w1, W2}represent main window W1 and the cross correlation value from window W2, ξ _{w1, W3}represent main window W1 and the cross correlation value from window W3, and

${\mathrm{\ξ}}_{W1,W2}=\frac{1}{1+(n_{W1}+n_{W2})\frac{n_{W1}{{\mathrm{\γ}}_{W1}}^2{c_{W1,W2}}^2+n_{W3}{{\mathrm{\γ}}_{W2}}^2}{n_{W1}{n}_{W2}{(c_{W1,W2}-1)}^2}}$

${\mathrm{\ξ}}_{W1,W3}=\frac{1}{1+(n_{W1}+n_{W3})\frac{n_{W1}{{\mathrm{\γ}}_{W1}}^2{c_{W1,W3}}^2+n_{W3}{{\mathrm{\γ}}_{W3}}^2}{n_{W1}{n}_{W3}{(c_{W1,W3}-1)}^2}}$

Wherein, n _w1for the pixel number of main window W1, n _w2for the pixel number from window W2, n _w3for the pixel number from window W3, for main window W1 with from window W2 regional average value ratio, for main window W1 with from window W3 regional average value ratio, γ _w1for the coefficient of variation of main window W1, γ _w2for the coefficient of variation from window W2, γ _w3for the coefficient of variation from window W3.It should be noted that, by the detection of simple crosscorrelation operator, can will reduce the false response being caused by isolated point, and then the extraction accuracy that improves roadway characteristic.

Then the rate information R ' of two exploded views of the described image block cross correlation value ξ corresponding with two exploded views of described image block passed through r ', ξ ∈ [0,1] merges, and the result after merging is carried out to simple thresholding, obtains described two roadway characteristics that exploded view is corresponding; In the present embodiment, as shown in Figure 12 a and 12b, be respectively the roadway characteristic of directions X exploded view and the roadway characteristic of Y-direction exploded view.

Finally two roadway characteristics corresponding to exploded view of described image block are merged to the roadway characteristic that obtains described image block; In the present embodiment, as shown in figure 13, be the road feature extraction result of described image block.

Understandable, other image blocks for described pretreated image also can pass through the above-mentioned corresponding roadway characteristic of other image blocks of procedure extraction, finally the roadway characteristic of each image block is merged to the roadway characteristic that can obtain described pretreated image, the routine techniques means that concrete mode is those skilled in the art, do not repeat them here.The roadway characteristic result that last original SAR image obtains by the technical scheme of the embodiment of the present invention, as shown in figure 14.

For the validity of S201 to S202 is described, the present embodiment also only carries out after autoregression filtering described image block to carry out road extraction result, original SAR image is only carried out after nonlinear quantization described image block carry out the result of road extraction and original SAR image do not carried out to the result that any pre-service directly carries out road extraction to described image block and analyze and contrast to original SAR image.The result of analyzing and contrast is as follows:

As shown in figure 15, for the described image block of original SAR image being carried out to the result of manual extraction roadway characteristic, the roadway characteristic shown in this figure is completely consistent with the roadway characteristic of original SAR image, as the reference of analyzing and contrasting;

As shown in Figure 16 a, for original SAR image is only carried out, after autoregression filtering, described image block is carried out to road extraction result; As shown in Figure 16 b, for original SAR image only being carried out to described image block is carried out after nonlinear quantization the result of road extraction; As Figure 16 c be depicted as to original SAR image do not carry out any pre-service directly to as described in image block carry out the result of road extraction.

Figure 15 is contrasted with Figure 13, Figure 16 a, Figure 16 b, Figure 16 c respectively, can find out significantly, in Figure 13, the extraction ratio of road is all best with the roadway characteristic quality extracting.

By extraction ratio, accuracy and 3 parameters of quality factor are assessed Figure 13, Figure 16 a, Figure 16 b, Figure 16 c, can obtain result as shown in table 1,

?	Extraction ratio (%)	Accuracy (%)	Quality factor (%)
				Figure 13	91.9	93.1	86.1
Figure 16 (a)	83.4	96.6	81.0
				Figure 16 (b)	90.9	93.3	85.3
Figure 16 (c)	85.5	97.7	83.8

Table 1

Also can find out from table 1, the extraction ratio of the roadway characteristic of Figure 13 and quality factor are apparently higher than the extraction effect of other 3 kinds of modes.

The present embodiment provides a kind of method of extracting roadway characteristic, by nonlinear quantization and the image filtering method based on Multiscale Autoregressive model, original SAR image is carried out pre-service and by the road extraction algorithm based on ratio and direction fusion, described pretreated image carried out to road feature extraction, can efficiently and exactly extract the roadway characteristic in SAR image.

Referring to Figure 17, the structure of a kind of equipment 170 that extracts roadway characteristic providing for the embodiment of the present invention, this equipment 170 can comprise: pretreatment unit 1701 and road feature extraction unit 1702, wherein,

Described pretreatment unit 1701, for original synthetic-aperture radar SAR image being carried out to pre-service by nonlinear quantization and the image filtering method based on Multiscale Autoregressive MAR model, obtains the pretreated image of described SAR image;

Described road feature extraction unit 1702, for by the road extraction algorithm based on ratio and direction fusion, described pretreated image being carried out to road feature extraction, obtains the roadway characteristic of described original SAR image.

Exemplary, referring to Figure 18, described pretreatment unit 1701 can comprise nonlinear quantization subelement 17011 and filtering subelement 17012, wherein,

Described nonlinear quantization subelement 17011 can be for, removes a small amount of supersaturation point of described original SAR image by the grey level histogram of described original SAR image, and obtain the dynamic range of the gray scale of described original SAR image;

And by Nonlinear extension algorithm by the dynamic range expansion of the gray scale of described original SAR image between 0-255, obtain the SAR image after nonlinear quantization.

Described filtering subelement 17012 can carry out pyramid decomposition for the SAR image to after described nonlinear quantization, obtains the SAR image sequence Φ that resolution reduces successively _p=[x ₀x ₁... x _p];

Wherein, p+1 is described SAR image sequence Φ _pin SAR image number, x ₀for described SAR image sequence Φ _pthe SAR image that intermediate-resolution is the highest, x _pfor described SAR image sequence Φ _pthe SAR image that intermediate-resolution is minimum; x _lfor x _l-1obtain according to formula (11):

$x_l(m,n)=\underset{i=2m}{\overset{2m+1}{\mathrm{\Σ}}}\underset{j=2n}{\overset{2n+1}{\mathrm{\Σ}}}{x}_{l-1}(i,j)---\left(11\right)$

Wherein, x _i(m, n) presentation video x _ithe pixel value of the capable n row of m, x _l-1(i, j) presentation video x _l-1the pixel value of the capable j of i row, l is greater than zero and be less than the natural number of p;

And according to described SAR image sequence Φ _pbuild MAR model; Wherein, described MAR model is suc as formula shown in (12):

$x_0=\underset{i=1}{\overset{p}{\mathrm{\Σ}}}{a}_i{x}_i+\mathrm{\ϵ}---\left(12\right)$

Wherein, a=[a ₁a ₂... a _p] be the autoregressive coefficient of described MAR model, ε is error term, p is the order of described MAR model;

And according to default estimation criterion, by least square method, the autoregressive coefficient a of described MAR model is estimated, obtain described SAR image sequence Φ _pthe SAR image x that intermediate-resolution is the highest ₀predicted picture described predicted picture wherein, predicted picture for the pretreated image of described SAR image, shown in (13):

${\hat{x}}_0=\underset{i=1}{\overset{p}{\mathrm{\Σ}}}{a}_i{x}_i---\left(13\right).$

Exemplary, referring to Figure 18, described road feature extraction unit 1702 comprises that piecemeal subelement 17021, image block characteristics extract subelement 17022 and splicing subelement 17023, wherein,

Described piecemeal subelement 17021, for described pretreated image is carried out to piecemeal, obtains at least one image block;

Described image block characteristics extracts subelement 17022, for the roadway characteristic of each image block being extracted by the road extraction algorithm based on ratio and direction fusion;

Described splicing subelement 17023, for the roadway characteristic after the extraction of each image block is spliced, obtains the roadway characteristic of described original SAR image.

Further, described image block characteristics extract subelement 17022 for:

The detection of sliding on described image block by rotatable ROA operator, obtains rate information R and the directional information θ of described image block; Wherein, described rotatable ROA operator comprise a main window W1 and two centrosymmetric from window W2 and W3 with main window;

Convert according to the rate information R of described image block, directional information θ and La Dong Radon the road principal direction of obtaining described image block and according to the road principal direction of described image block the best of intersecting roads in described image block is decomposed, obtain two exploded views of described image block;

All use ROA operator and simple crosscorrelation operator to detect road information separately to two of described image block exploded views, and the road information of two exploded views of described image block is merged to the roadway characteristic that obtains described image block.

Preferably, described image block characteristics extract subelement 17022 for: calculate described rotatable ROA operator current location (x ₀, y ₀) maximum rate value R (x ₀, y ₀) and the maximum rate value R (x of described current location ₀, y ₀) corresponding direction θ (x ₀, y ₀);

Described rotatable ROA operator is slided into next position (x ₁, y ₁) calculate the maximum rate value R (x of next position described in described rotatable ROA operator ₁, y ₁) and the maximum rate value R (x of described next position ₁, y ₁) corresponding direction θ (x ₁, y ₁), until described rotatable ROA operator slides complete by described image block;

Wherein, calculating the direction θ that the maximum rate value R of described rotatable ROA operator current location and the maximum rate value R of described current location are corresponding is obtained by formula (14) to formula (16):

$r(x,y,{\mathrm{\θ}}_k)=\min (1-\frac{M_1(x,y,{\mathrm{\θ}}_k)}{M_2(x,y,{\mathrm{\θ}}_k)},1-\frac{M_1(x,y,{\mathrm{\θ}}_k)}{M_3(x,y,{\mathrm{\θ}}_k)})---\left(14\right)$

$R(x,y)=\underset{k=1:N}{\max }\left(r(x,y,{\mathrm{\θ}}_k)\right)---\left(15\right)$

$\mathrm{\θ}(x,y)=\underset{{\mathrm{\θ}}_k,k=1:N}{\arg }\max \left\{r(x,y,{\mathrm{\θ}}_k)\right\}---\left(16\right)$

Wherein, r (x, y, θ _k) be that the described rotatable ROA operator of current location is at direction θ _kratio, M ₁(x, y, θ _k), M ₂(x, y, θ _k) and M ₃(x, y, θ _k) represent that respectively window W1, the W2 of described rotatable ROA operator and W3 are at direction θ _kpixel value average; N represents the total degree that described rotatable ROA operator is rotated; represent that current described rotatable ROA operator is in directive ratio maximal value; $\underset{{\mathrm{\θ}}_k,k=1:N}{\arg }\max \left\{r(x,y,{\mathrm{\θ}}_k)\right\}$ Represent $\underset{k=1:N}{\max }\left(r(x,y,{\mathrm{\θ}}_k)\right)$ Corresponding direction.

Preferably, described image block characteristics extract subelement 17022 for:

The rate information R of the each pixel of described image block and directional information θ are transformed to the direction vector (p under Descartes's rectangular coordinate that described each pixel is corresponding according to formula (17) _x, p _y);

$\left\{\begin{array}{c}{p}_x=R\cos \mathrm{\θ}\\ p_y=R\sin \mathrm{\θ}\end{array}\right.---\left(17\right)$

And obtain the local peaking of formula (18) in described image block according to the rate information R of described image block;

Wherein, D is integral domain, represents described image block plane in the present embodiment; R (x, y) is the value that the rate information R of described image block locates at (x, y); δ is Dirac function; ρ be the initial point of Descartes's right-angle plane to the distance of straight line, for the initial point of described Descartes's right-angle plane is to the angle of the vertical line of straight line and the x axle of described Descartes's right-angle plane;

And described in inciting somebody to action what direction was corresponding owns local peaking superpose and obtain described image block and exist there is the probability level of straight line in direction

And searching probability index maximal value, and by described probability level the corresponding angle of maximal value as the road principal direction of described image block, wherein, described probability level the corresponding angle of maximal value through type (19) represents:

And by the direction vector (p under corresponding Descartes's rectangular coordinate of the each pixel of described image block _x, p _y) the extremely road principal direction of described image block of through type (20) rotation the best that obtains the intersecting roads to described image block is decomposed, and wherein, formula (10) is:

Wherein, (p _x, p _y) be the direction vector under the Descartes's rectangular coordinate before the each pixel rotation of described image block, (p _x', p _y') be the direction vector under postrotational Descartes's rectangular coordinate of the each pixel of described image block.

Preferably, described image block characteristics extract subelement 17022 for:

Two of described image block exploded views are detected by ROA operator respectively, obtain the rate information R ' of described two exploded views;

Two of described image block exploded views are detected by simple crosscorrelation operator respectively, obtain cross correlation value ξ=min (ξ corresponding to two exploded views of described image block _{w1, W2}, ξ _{w1, W3}),

Wherein, ξ _{w1, W2}represent main window W1 and the cross correlation value from window W2, ξ _{w1, W3}represent main window W1 and the cross correlation value from window W3, and

${\mathrm{\ξ}}_{W1,W2}=\frac{1}{1+(n_{W1}+n_{W2})\frac{n_{W1}{{\mathrm{\γ}}_{W1}}^2{c_{W1,W2}}^2+n_{W3}{{\mathrm{\γ}}_{W2}}^2}{n_{W1}{n}_{W2}{(c_{W1,W2}-1)}^2}}$

${\mathrm{\ξ}}_{W1,W3}=\frac{1}{1+(n_{W1}+n_{W3})\frac{n_{W1}{{\mathrm{\γ}}_{W1}}^2{c_{W1,W3}}^2+n_{W3}{{\mathrm{\γ}}_{W3}}^2}{n_{W1}{n}_{W3}{(c_{W1,W3}-1)}^2}}$

Wherein, n _w1for the pixel number of main window W1, n _w2for the pixel number from window W2, n _w3for the pixel number from window W3, for main window W1 with from window W2 regional average value ratio, for main window W1 with from window W3 regional average value ratio, γ _w1for the coefficient of variation of main window W1, γ _w2for the coefficient of variation from window W2, γ _w3for the coefficient of variation from window W3.It should be noted that, by the detection of simple crosscorrelation operator, the isolated point in exploded view can be removed, is road because these isolated points can be mistaken as conventionally, thereby affects the extraction accuracy of roadway characteristic.；

And the rate information R' of two exploded views of the described image block cross correlation value ξ corresponding with two exploded views of described image block passed through r', ξ ∈ [0,1] merges, and by the result passing threshold processing after merging, obtains two roadway characteristics that exploded view is corresponding of described image block;

And two roadway characteristics corresponding to exploded view of described image block are merged to the roadway characteristic that obtains described image block.

The present embodiment provides a kind of equipment that extracts roadway characteristic, by nonlinear quantization and the image filtering method based on Multiscale Autoregressive model, original SAR image is carried out pre-service and by the road extraction algorithm based on ratio and direction fusion, described pretreated image carried out to road feature extraction, can efficiently and exactly extract the roadway characteristic in SAR image.

Those skilled in the art should understand, embodiments of the invention can be provided as method, system or computer program.Therefore, the present invention can adopt hardware implementation example, implement software example or the form in conjunction with the embodiment of software and hardware aspect.And the present invention can adopt the form at one or more upper computer programs of implementing of computer-usable storage medium (including but not limited to magnetic disk memory and optical memory etc.) that wherein include computer usable program code.

The present invention is with reference to describing according to process flow diagram and/or the block scheme of the method for the embodiment of the present invention, equipment (system) and computer program.Should understand can be by the flow process in each flow process in computer program instructions realization flow figure and/or block scheme and/or square frame and process flow diagram and/or block scheme and/or the combination of square frame.Can provide these computer program instructions to the processor of multi-purpose computer, special purpose computer, Embedded Processor or other programmable data processing device to produce a machine, the instruction that makes to carry out by the processor of computing machine or other programmable data processing device produces the device for realizing the function of specifying at flow process of process flow diagram or multiple flow process and/or square frame of block scheme or multiple square frame.

These computer program instructions also can be stored in energy vectoring computer or the computer-readable memory of other programmable data processing device with ad hoc fashion work, the instruction that makes to be stored in this computer-readable memory produces the manufacture that comprises command device, and this command device is realized the function of specifying in flow process of process flow diagram or multiple flow process and/or square frame of block scheme or multiple square frame.

These computer program instructions also can be loaded in computing machine or other programmable data processing device, make to carry out sequence of operations step to produce computer implemented processing on computing machine or other programmable devices, thereby the instruction of carrying out is provided for realizing the step of the function of specifying in flow process of process flow diagram or multiple flow process and/or square frame of block scheme or multiple square frame on computing machine or other programmable devices.

The above, be only preferred embodiment of the present invention, is not intended to limit protection scope of the present invention.

Claims (16)

1. a method of extracting roadway characteristic, is characterized in that, described method comprises:

By nonlinear quantization and the image filtering method based on Multiscale Autoregressive MAR model, original synthetic-aperture radar SAR image is carried out to pre-service, obtain the pretreated image of described SAR image;

By the road extraction algorithm based on ratio and direction fusion, described pretreated image is carried out to road feature extraction, obtain the roadway characteristic of described original SAR image.

2. method according to claim 1, is characterized in that, describedly by nonlinear quantization and the image filtering method based on MAR model, original SAR image is carried out to pre-service, comprising:

Remove a small amount of supersaturation point of described original SAR image by the grey level histogram of described original SAR image, and obtain the dynamic range of the gray scale of described original SAR image;

By Nonlinear extension algorithm by the dynamic range expansion of the gray scale of described original SAR image between 0-255, obtain the SAR image after nonlinear quantization.

3. method according to claim 2, is characterized in that, describedly by nonlinear quantization and the image filtering method based on MAR model, original SAR image is carried out to pre-service, also comprises:

SAR image after described nonlinear quantization is carried out to pyramid decomposition, obtain the SAR image sequence Φ that resolution reduces successively _p=[x ₀x ₁... x _p];

Wherein, p+1 is described SAR image sequence Φ _pin SAR image number, x ₀for described SAR image sequence Φ _pthe SAR image that intermediate-resolution is the highest, x _pfor described SAR image sequence Φ _pthe SAR image that intermediate-resolution is minimum; x _lfor x _l-1obtain according to formula (1):

$x_l(m,n)=\underset{i=2m}{\overset{2m+1}{\mathrm{\Σ}}}\underset{j=2n}{\overset{2n+1}{\mathrm{\Σ}}}{x}_{l-1}(i,j)---\left(1\right)$

Wherein, x _i(m, n) presentation video x _ithe pixel value of the capable n row of m, x _l-1(i, j) presentation video x _l-1the pixel value of the capable j of i row, l is greater than zero and be less than the natural number of p;

According to described SAR image sequence Φ _pbuild MAR model; Wherein, described MAR model is suc as formula shown in (2):

$x_0=\underset{i=1}{\overset{p}{\mathrm{\Σ}}}{a}_i{x}_i+\mathrm{\ϵ}---\left(2\right)$

Wherein, a=[a ₁a ₂... a _p] be the autoregressive coefficient of described MAR model, ε is error term, p is the order of described MAR model;

According to default estimation criterion, by least square method, the autoregressive coefficient a of described MAR model is estimated, obtain described SAR image sequence Φ _pthe SAR image x that intermediate-resolution is the highest ₀predicted picture wherein, predicted picture for the pretreated image of described SAR image, shown in (3):

${\hat{x}}_0=\underset{i=1}{\overset{p}{\mathrm{\Σ}}}{a}_i{x}_i---\left(3\right).$

4. method according to claim 1, is characterized in that, describedly by the road extraction algorithm merging based on ratio and direction, described pretreated image is carried out to road feature extraction, comprising:

Described pretreated image is carried out to piecemeal, obtain at least one image block;

By the road extraction algorithm based on ratio and direction fusion, the roadway characteristic of each image block is extracted;

Roadway characteristic after the extraction of each image block is spliced, obtain the roadway characteristic of described original SAR image.

5. method according to claim 4, is characterized in that, describedly by the road extraction algorithm merging based on ratio and direction, the roadway characteristic of each image block is extracted, and comprising:

The detection of sliding on described image block by rotatable ROA operator, obtains rate information R and the directional information θ of described image block; Wherein, described rotatable ROA operator comprise a main window W1 and two centrosymmetric from window W2 and W3 with main window;

Convert according to the rate information R of described image block, directional information θ and La Dong Radon the road principal direction of obtaining described image block and according to the road principal direction of described image block the best of intersecting roads in described image block is decomposed, obtain two exploded views of described image block;

The Road Detection operator extraction that two of described image block exploded views are all used to ROA operator and simple crosscorrelation Operator Fusion road information separately, and the road information of two exploded views of described image block is merged, thereby obtain the roadway characteristic of described image block.

6. method according to claim 5, is characterized in that, described by image block is slided on the described image block detection of rotatable ROA operator, obtains rate information R and the directional information θ of described image block, comprising:

Calculate described rotatable ROA operator current location (x ₀, y ₀) maximum rate value R (x ₀, y ₀) and the maximum rate value R (x of described current location ₀, y ₀) corresponding direction θ (x ₀, y ₀);

Described rotatable ROA operator is slided into next position (x ₁, y ₁) calculate the maximum rate value R (x of next position described in described rotatable ROA operator ₁, y ₁) and the maximum rate value R (x of described next position ₁, y ₁) corresponding direction θ (x ₁, y ₁), until described rotatable ROA operator slides complete by described image block;

Wherein, calculating the direction θ that the maximum rate value R of described rotatable ROA operator current location and the maximum rate value R of described current location are corresponding is obtained by formula (4) to formula (6):

$r(x,y,{\mathrm{\θ}}_k)=\min (1-\frac{M_1(x,y,{\mathrm{\θ}}_k)}{M_2(x,y,{\mathrm{\θ}}_k)},1-\frac{M_1(x,y,{\mathrm{\θ}}_k)}{M_3(x,y,{\mathrm{\θ}}_k)})---\left(4\right)$

$R(x,y)=\underset{k=1:N}{\max }\left(r(x,y,{\mathrm{\θ}}_k)\right)---\left(5\right)$

$\mathrm{\θ}(x,y)=\underset{{\mathrm{\θ}}_k,k=1:N}{\arg }\max \left\{r(x,y,{\mathrm{\θ}}_k)\right\}---\left(6\right)$

Wherein, r (x, y, θ _k) be that the described rotatable ROA operator of current location is at direction θ _kratio, M ₁(x, y, θ _k), M ₂(x, y, θ _k) and M ₃(x, y, θ _k) represent that respectively window W1, the W2 of described rotatable ROA operator and W3 are at direction θ _kpixel value average; N represents the total degree that described rotatable ROA operator is rotated; represent that current described rotatable ROA operator is in directive ratio maximal value; $\underset{{\mathrm{\θ}}_k,k=1:N}{\arg }\max \left\{r(x,y,{\mathrm{\θ}}_k)\right\}$ Represent $\underset{k=1:N}{\max }\left(r(x,y,{\mathrm{\θ}}_k)\right)$ Corresponding direction.

7. method according to claim 5, is characterized in that, converts according to the rate information R of described image block, directional information θ and Radon the road principal direction of obtaining described image block and according to the road principal direction of described image block the best of intersecting roads in described image block is decomposed, obtains two exploded views of described image block, comprising:

The rate information R of the each pixel of described image block and directional information θ are transformed to the direction vector (p under Descartes's rectangular coordinate that described each pixel is corresponding according to formula (7) _x, p _y);

$\left\{\begin{array}{c}{p}_x=R\cos \mathrm{\θ}\\ p_y=R\sin \mathrm{\θ}\end{array}\right.---\left(7\right)$

Obtain the local peaking of formula (8) in described image block according to the rate information R of described image block;

Wherein, D is integral domain, represents described image block plane in the present embodiment; R (x, y) is the value that the rate information R of described image block locates at (x, y); δ is Dirac function; ρ be the initial point of Descartes's right-angle plane to the distance of straight line, for the initial point of described Descartes's right-angle plane is to the angle of the vertical line of straight line and the x axle of described Descartes's right-angle plane;

Described in inciting somebody to action what direction was corresponding owns local peaking superpose and obtain described image block and exist there is the probability level of straight line in direction

Searching probability index maximal value, and by described probability level the corresponding angle of maximal value as the road principal direction of described image block, wherein, described probability level the corresponding angle of maximal value through type (9) represents:

By the direction vector (p under corresponding Descartes's rectangular coordinate of the each pixel of described image block _x, p _y) the extremely road principal direction of described image block of through type (10) rotation the best that obtains the intersecting roads to described image block is decomposed, and wherein, formula (10) is:

Wherein, (p _x, p _y) be the direction vector under Descartes's rectangular coordinate before the each pixel rotation of described image block, (p _x', p _y') be the direction vector under Descartes's rectangular coordinate after the each pixel rotation of described image block.

8. method according to claim 5, it is characterized in that, describedly all use ROA operator and simple crosscorrelation operator to detect road information separately to two of described image block exploded views, and the road information of two exploded views of described image block is merged to the roadway characteristic that obtains described image block, comprising:

Two of described image block exploded views are detected by ROA operator respectively, obtain the rate information R ' of described two exploded views;

Two of described image block exploded views are detected by simple crosscorrelation operator respectively, obtain the cross correlation value ξ corresponding to two exploded views of described image block;

The cross correlation value ξ that the rate information R ' of two exploded views of described image block is corresponding with two exploded views of described image block merges, and by the result passing threshold processing after merging, obtains two roadway characteristics that exploded view is corresponding of described image block;

Two roadway characteristics corresponding to exploded view of described image block are merged to the roadway characteristic that obtains described image block.

9. an equipment that extracts roadway characteristic, is characterized in that, described equipment comprises: pretreatment unit and road feature extraction unit, wherein,

Described pretreatment unit, for original synthetic-aperture radar SAR image being carried out to pre-service by nonlinear quantization and the image filtering method based on Multiscale Autoregressive MAR model, obtains the pretreated image of described SAR image;

Described road feature extraction unit, for by the road extraction algorithm based on ratio and direction fusion, described pretreated image being carried out to road feature extraction, obtains the roadway characteristic of described original SAR image.

10. equipment according to claim 9, it is characterized in that, described pretreatment unit comprises nonlinear quantization subelement, remove a small amount of supersaturation point of described original SAR image for the grey level histogram by described original SAR image, and obtain the dynamic range of the gray scale of described original SAR image;

And by Nonlinear extension algorithm by the dynamic range expansion of the gray scale of described original SAR image between 0-255, obtain the SAR image after nonlinear quantization.

11. equipment according to claim 10, is characterized in that, described pretreatment unit also comprises filtering subelement, carry out pyramid decomposition for the SAR image to after described nonlinear quantization, obtain the SAR image sequence Φ that resolution reduces successively _p=[x ₀x ₁... x _p];

Wherein, p+1 is described SAR image sequence Φ _pin SAR image number, x ₀for described SAR image sequence Φ _pthe SAR image that intermediate-resolution is the highest, x _pfor described SAR image sequence Φ _pthe SAR image that intermediate-resolution is minimum; x _lfor x _l-1obtain according to formula (11):

$x_l(m,n)=\underset{i=2m}{\overset{2m+1}{\mathrm{\Σ}}}\underset{j=2n}{\overset{2n+1}{\mathrm{\Σ}}}{x}_{l-1}(i,j)---\left(11\right)$

Wherein, x _i(m, n) presentation video x _ithe pixel value of the capable n row of m, x _l-1(i, j) presentation video x _l-1the pixel value of the capable j of i row, l is greater than zero and be less than the natural number of p;

And according to described SAR image sequence Φ _pbuild MAR model; Wherein, described MAR model is suc as formula shown in (12):

$x_0=\underset{i=1}{\overset{p}{\mathrm{\Σ}}}{a}_i{x}_i+\mathrm{\ϵ}---\left(12\right)$

Wherein, a=[a ₁a ₂... a _p] be the autoregressive coefficient of described MAR model, ε is error term, p is the order of described MAR model;

And according to default estimation criterion, by least square method, the autoregressive coefficient a of described MAR model is estimated, obtain described SAR image sequence Φ _pthe SAR image x that intermediate-resolution is the highest ₀predicted picture described predicted picture wherein, predicted picture for the pretreated image of described SAR image, shown in (13):

${\hat{x}}_0=\underset{i=1}{\overset{p}{\mathrm{\Σ}}}{a}_i{x}_i---\left(13\right).$

12. equipment according to claim 9, is characterized in that, described road feature extraction unit comprises that piecemeal subelement, image block characteristics extract subelement and splicing subelement, wherein,

Described piecemeal subelement, for described pretreated image is carried out to piecemeal, obtains at least one image block;

Described image block characteristics extracts subelement, for the roadway characteristic of each image block being extracted by the road extraction algorithm based on ratio and direction fusion;

Described splicing subelement, for the roadway characteristic after the extraction of each image block is spliced, obtains the roadway characteristic of described original SAR image.

13. equipment according to claim 12, is characterized in that, described image block characteristics extracts subelement and is used for:

The detection of sliding on described image block by rotatable ROA operator, obtains rate information R and the directional information θ of described image block; Wherein, described rotatable ROA operator comprise a main window W1 and two centrosymmetric from window W2 and W3 with main window;

Convert according to the rate information R of described image block, directional information θ and La Dong Radon the road principal direction of obtaining described image block and according to the road principal direction of described image block the best of intersecting roads in described image block is decomposed, obtain two exploded views of described image block;

All use ROA operator and simple crosscorrelation operator to detect road information separately to two of described image block exploded views, and the road information of two exploded views of described image block is merged to the roadway characteristic that obtains described image block.

14. equipment according to claim 13, is characterized in that, described image block characteristics extracts subelement and is used for:

Calculate described rotatable ROA operator current location (x ₀, y ₀) maximum rate value R (x ₀, y ₀) and the maximum rate value R (x of described current location ₀, y ₀) corresponding direction θ (x ₀, y ₀);

Described rotatable ROA operator is slided into next position (x ₁, y ₁) calculate the maximum rate value R (x of next position described in described rotatable ROA operator ₁, y ₁) and the maximum rate value R (x of described next position ₁, y ₁) corresponding direction θ (x ₁, y ₁), until described rotatable ROA operator slides complete by described image block;

Wherein, calculating the direction θ that the maximum rate value R of described rotatable ROA operator current location and the maximum rate value R of described current location are corresponding is obtained by formula (14) to formula (16):

$r(x,y,{\mathrm{\θ}}_k)=\min (1-\frac{M_1(x,y,{\mathrm{\θ}}_k)}{M_2(x,y,{\mathrm{\θ}}_k)},1-\frac{M_1(x,y,{\mathrm{\θ}}_k)}{M_3(x,y,{\mathrm{\θ}}_k)})---\left(14\right)$

$R(x,y)=\underset{k=1:N}{\max }\left(r(x,y,{\mathrm{\θ}}_k)\right)---\left(15\right)$

$\mathrm{\θ}(x,y)=\underset{{\mathrm{\θ}}_k,k=1:N}{\arg }\max \left\{r(x,y,{\mathrm{\θ}}_k)\right\}---\left(16\right)$

Wherein, r (x, y, θ _k) be that the described rotatable ROA operator of current location is at direction θ _kratio, M ₁(x, y, θ _k), M ₂(x, y, θ _k) and M ₃(x, y, θ _k) represent that respectively window W1, the W2 of described rotatable ROA operator and W3 are at direction θ _kpixel value average; N represents the total degree that described rotatable ROA operator is rotated; represent that current described rotatable ROA operator is in directive ratio maximal value; $\underset{{\mathrm{\θ}}_k,k=1:N}{\arg }\max \left\{r(x,y,{\mathrm{\θ}}_k)\right\}$ Represent $\underset{k=1:N}{\max }\left(r(x,y,{\mathrm{\θ}}_k)\right)$ Corresponding direction.

15. equipment according to claim 13, is characterized in that, described image block characteristics extracts subelement and is used for:

The rate information R of the each pixel of described image block and directional information θ are transformed to the direction vector (p under Descartes's rectangular coordinate that described each pixel is corresponding according to formula (17) _x, p _y);

$\left\{\begin{array}{c}{p}_x=R\cos \mathrm{\θ}\\ p_y=R\sin \mathrm{\θ}\end{array}\right.---\left(17\right)$

And obtain the local peaking of formula (18) in described image block according to the rate information R of described image block;

Wherein, D is integral domain, represents described image block plane in the present embodiment; R (x, y) is the value that the rate information R of described image block locates at (x, y); δ is Dirac function; ρ be the initial point of Descartes's right-angle plane to the distance of straight line, for the initial point of described Descartes's right-angle plane is to the angle of the vertical line of straight line and the x axle of described Descartes's right-angle plane;

And described in inciting somebody to action what direction was corresponding owns local peaking superpose and obtain described image block and exist there is the probability level of straight line in direction

And searching probability index maximal value, and by described probability level the corresponding angle of maximal value as the road principal direction of described image block, wherein, described probability level the corresponding angle of maximal value through type (19) represents:

And by the direction vector (p under corresponding Descartes's rectangular coordinate of the each pixel of described image block _x, p _y) the extremely road principal direction of described image block of through type (20) rotation the best that obtains the intersecting roads to described image block is decomposed, and wherein, formula (10) is:

Wherein, (p _x, p _y) be the direction vector under the Descartes's rectangular coordinate before the each pixel rotation of described image block, (p _x', p _y') be the direction vector under the postrotational Descartes's rectangular coordinate of the each pixel of described image block.

16. extraction equipments according to claim 13, is characterized in that, described image block characteristics extracts subelement and is used for:

Two of described image block exploded views are detected by ROA operator respectively, obtain the rate information R ' of described two exploded views;

Two of described image block exploded views are detected by simple crosscorrelation operator respectively, obtain the cross correlation value ξ corresponding to two exploded views of described image block;

And the rate information R ' of two exploded views of the described image block cross correlation value ξ corresponding with two exploded views of described image block merged, and by the result passing threshold processing after merging, obtain two roadway characteristics that exploded view is corresponding of described image block;

And two roadway characteristics corresponding to exploded view of described image block are merged to the roadway characteristic that obtains described image block.