CN104700427A - 一种基于sar图像的道路损毁信息提取方法 - Google Patents

Abstract

本发明公开了一种基于SAR图像的道路损毁信息提取方法，通过SAR卫星或航拍SAR获取双时相的原始SAR图像，然后对原始SAR图像进行预处理，得到增强图像，再从增强图像中提取出边缘信息，计算边缘的方向特征，结合灰度信息、边缘信息、方向信息采用粒子滤波方法和SNAKE方法提取道路中心线，最后，以道路中心线为基础提取出待检测图像的道路损毁信息。本发明是利用道路在SAR图像中的辐射特征、几何特征等特征信息，有效地实现了道路中心线的提取和道路损毁信息提取。

Description

一种基于SAR图像的道路损毁信息提取方法

技术领域

本发明属于图像处理技术领域，更为具体地讲，涉及一种基于SAR图像的道路损毁信息提取方法。

背景技术

中国是自然灾害多发国家，洪涝地震等自然灾害作为不可抗因素，通常会对人们的生命财产安全造成巨大威胁，因此抗灾减灾的需求仍然是当今人类发展中需要面对的重要议题。

道路作为一种典型的人造地物，是构成现代交通体系的主要部分，具有重要的地理、政治、经济和军事意义。道路信息在应急响应和灾害救援领域中有着重要的价值，其中道路损毁信息对灾害评估和灾害救援具有重要的指导意义。目前在道路损毁信息提取领域，主要采用光学遥感图像，光学图像具有直观性，很容易识别道路信息，但光学图像的获取受天气和时间因素影响较大，不能保证随时获得清晰和有效的图像，而灾害发生时，一般都伴随恶劣天气，一定程度上影响了数据获取的实时性。SAR图像具有全天时、全天候和大面积成像的优点，在应用中体现为快速响应和高实时性特点，且近年SAR图像分辨率不断提高，如TerraSar X聚束模式距离向分辨率可以达到0.5m，道路在SAR图像中体现出暗黑条带状特征，因此基于高分辨率SAR图像的道路损毁信息提取具有很高的研究和实用价值。

现阶段基于SAR图像的道路损毁信息提取技术仍鲜见报道，如何利用SAR图像的优势和特点，进行道路特征信息提取和道路损毁信息提取，是目前研究的热点。

发明内容

本发明的目的在于克服现有技术的不足，提供一种基于SAR图像的道路损毁信息提取方法，以实现具有快速响应和高实时性特点的道路损毁信息提取。

为实现上述发明目的，本发明一种基于SAR图像的道路损毁信息提取方法，其特征在于，包括以下步骤：

(1)、获取双时相原始SAR图像

在同一监测点，分别获取道路损毁发生前、后两幅双时相原始SAR图像，标记为基准原始SAR图像和待检测原始SAR图像；

(2)、对双时相原始SAR图像进行预处理

对双时相原始SAR图像进行图像配准，得到配准后的基准SAR图像和待检测SAR图像，再对基准SAR图像进行增强处理，得到增强图像；

(3)、边缘信息提取

通过双边缘检测算子Duda对增强图像进行边缘信息检测，从该增强图像中提取出边缘信息图像；

(4)、边缘方向信息获取

对边缘信息图像进行方向离散值计算，获取到边缘信息图像中逐像元点的方向信息，从而得到边缘方向信息图像；

(5)、道路中心线提取

根据增强图像、边缘信息图像以及边缘方向信息图像中的特征，采用粒子滤波方法和SNAKE方法提出道路的中心线；

(6)、提取道路损毁信息

以步骤(5)中提取的中心线为基础，获取待检测SAR图像中道路中心线上各点的边缘损毁指数和灰度损毁指数，再将边缘损毁指数和灰度损毁指数求并集，得到道路的损毁结果。

本发明的发明目的是这样实现的：

本发明基于SAR图像的道路损毁信息提取方法，通过SAR卫星或航拍SAR获取获取双时相的原始SAR图像，然后对原始SAR图像进行预处理，得到增强图像，再从增强图像中提取出边缘信息，计算边缘的方向特征，结合灰度信息、边缘信息、方向信息采用粒子滤波方法和SNAKE方法提取道路中心线，最后，以道路中心线为基础提取出待检测图像的道路损毁信息。本发明是利用道路在SAR图像中的辐射特征、几何特征等特征信息，有效地实现了道路中心线的提取和道路损毁信息提取。

附图说明

图1是本发明基于SAR图像的道路损毁信息提取方法流程图；

图2是双时相原始SAR图像；

图3是增强图像；

图4是边缘信息图像；

图5是边缘方向信息提取结果图像；

图6是粒子滤波方法提取道路中心线的流程图；

图7是道路中心线提取结果图像；

图8是道路损毁信息提取结果图像；

具体实施方式

下面结合附图对本发明的具体实施方式进行描述，以便本领域的技术人员更好地理解本发明。需要特别提醒注意的是，在以下的描述中，当已知功能和设计的详细描述也许会淡化本发明的主要内容时，这些描述在这里将被忽略。

实施例

图1是本发明基于SAR图像的道路损毁信息提取方法流程图。

在本实施例中，如图1所示，本发明基于SAR图像的道路损毁信息提取方法，包括以下步骤：

S1、获取双时相原始SAR图像

本发明需要在同一监测点，分别获取道路损毁发生前、后两幅双时相原始SAR图像，标记为基准原始SAR图像，如图2(a)所示，以及待检测原始SAR图像，如图2(b)所示；

在本实施例中，如图2所示，通过卫星遥感和航拍SAR成像系统在同一监测点拍摄的两幅SAR图像。其中，图2(a)可以辨认出较为清晰道路区域，即呈暗黑色条的线性条带区域；图2(b)中道路区域几乎不可见。

S2、对双时相原始SAR图像进行预处理

对双时相原始SAR图像进行图像配准，得到配准后的基准SAR图像和待检测SAR图像，再对基准SAR图像进行增强处理，得到增强图像，如图3所示；

在本实施例中，对基准SAR图像先进行DMFF滤波器，使低灰度线形结构增强，再平滑基准SAR图像，降低相干斑噪声影响，同时增强道路与周围环境的对比度，如图3所示，增强图像中道路区域均匀度及对比度均高于原基准图像。

S3、边缘信息提取

通过双边缘检测算子Duda对增强图像进行边缘信息检测，从该增强图像中提取出边缘信息图像。

在本实施例中，如图4所示，边缘信息图像是由基准增强图像经Duda算子检测得到，图像中像素点的灰度值越高，则该点处的边缘特征越明显，而道路的真实边缘又隐含在具有一定灰度值的像素点集合中。

S4、边缘方向信息获取

S4.1、对边缘信息图像通过预定的阈值进行二值化处理，得到二值分割图像；

S4.2、在本实施例中，采用如图5(b)所示的八方向计算模板逐点遍历二值分割图像中每个像元点的边缘方向，即计算每个像元点的方向离散值；

其中，Direction为边缘方向信息图；I为二值分割图像；d_k为方向k上的方向模板窗口，k∈[1,8]；为方向信息计算时的当前像元点坐标，N为边缘信息图像中像元点的个数；

通过上述处理后，可以获取到每个像元点的方向信息，所有像元点的方向信息再构成边缘方向信息图像，如图5(a)所示，图中的每一个像素点的灰度值表示该点的边缘方向值，局部区域内的像素点的灰度值一般较为接近。图5(b)为八方向离散值示意图，以22.5°为间隔将180度范围内的方向以1到8进行表示；在本实施例中，d_k为大小为7*7的矩形窗口，其中经过窗口中心点、方向为k的直线上所有点的加权值为1，窗口中其余点加权值均为0。

S5、道路中心线提取

S5.1、先用利用粒子滤波方法提取道路的中心线，在本实施例中，如图6所示，包括以下步骤：

S5.1.1、先构建一个道路中心点空间其中，m为道路中心线上像元点的个数；同时定义迭代次数为t，t∈[0,m]，为一下步骤使用；

S5.1.2、在基准SAR图像中的道路区域任意选取一个像元点point₀，将point₀作为道路种子点，分别构建两个大小为n的粒子集空间，即当前粒子集PC＝{pc(i)}，预测粒子集PP＝{pp(i)}，其中，i＝1,2,…，n， {xⁱ,yⁱ}代表序号为i的粒子的位置坐标，θⁱ为序号为i的粒子的方向，wⁱ为两粒子集中各粒子的粒子权重，且初始粒子权重设置为1；

在后续过迭代程中，第t次迭代中的当前粒子集PC等价于第t-1次迭代完成时的预测粒子集PP；在本实施例中，n＝1500。

S5.1.3、对种子点point₀进行采样，得到n个粒子元素，并将这n个粒子元素作为当前粒子集PC中的元素，同时将初始化迭代次数t设置为0，再进入步骤S5.1.4；

S5.1.4、如果t＜m，则令t＝t+1，再进入步骤S5.1.5；如果t＞m-1，则停止迭代，跳转至步骤S5.1.8；

S5.1.5、更新预测粒子集PP

通过状态转移方程计算第t次迭代时当前粒子集PC中每个粒子pc(i)的预测状态，并将该结果赋值给对应的预测粒子pp(i)；

$\left\{\begin{array}{c}{x}_t=x_{t-1}+d\cos \left({\mathrm{\θ}}_{t-1}\right)+u_1\\ y_t=y_{t-1}+d\sin \left({\mathrm{\θ}}_{t-1}\right)+u_2\\ {\mathrm{\θ}}_t={\mathrm{\θ}}_{t-1}+u_3\end{array}\right.$

其中，表示第t次迭代时预测粒子pp(i)的坐标，表示第t次迭代时预测粒子pp(i)的方向，预设的步长间隔d＝2，u₁,u₂,u₃分别为用于修正状态转移方程的随机噪声；

S5.1.6、对预测粒子集PP中逐粒子进行粒子权重计算，并进行粒子权重归一化，再利用得到的归一化权重估计预测道路点的状态；

其中，粒子权重计算的方法为：

T1)、获取预测粒子在第t次迭代时的灰度值边缘强度值r_t ⁱ和边缘方向信息离散值

设第i个预测粒子的坐标为(X_i,Y_i)，通过该坐标值在基准增强图像中查找出该预测粒子在第t次迭代时的灰度值在边缘信息图像中查找出该预测粒子的边缘强度值r_t ⁱ，在边缘方向信息图像中查找出该预测粒子的边缘方向信息离散值

T2)、计算预测粒子的最终权重

先根据r_t ⁱ，求出第i个预测粒子的灰度权重边缘权重方向权重

$w_g^i=\left\{\begin{array}{cc}1/|g_{t-1}-g_t^i|& \mathrm{ThGrayL}<g_t^i<\mathrm{ThGrayH}\\ 0& \mathrm{Others}\end{array}\right.$

$w_e^i=\frac{1}{\sqrt{2\mathrm{\&pi;}}\mathrm{\&sigma;}}{e}^{-\frac{{(1-r_t^i)}^2}{2{\mathrm{\&sigma;}}^2}}$

$w_d^i=1/|{\mathrm{\&theta;}}_{t-1}-{\mathrm{\&theta;}}_t^i|$

再计算出预测粒子的最终权重

$w_t^i=w_e^i+\mathrm{\&alpha;}{w}_g^i+\mathrm{\&beta;}{w}_d^i$

其中，ThGrayL、ThGrayH分别为预设的低灰度阈值和高灰度阈值，σ为常数，它描述了粒子集中边缘权重值的散度，α、β为权重比例调整参数；

粒子权重归一化的方法为：

$\mathrm{SumWeight}=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{w}_t^i$

$w_t^i=w_t^i/\mathrm{SumWeight}$

$\left\{\begin{array}{c}{x}_t=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_t^i{w}_t^i/\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{w}_t^i\\ y_t=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{y}_t^i{w}_t^i/\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{w}_t^i\\ {\mathrm{\&theta;}}_t=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&theta;}}_t^i{w}_t^i/\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{w}_t^i\end{array}\right.$

其中，表示预测粒子集中第i个预测粒子的状态，其中， 分别表示第t次迭代时预测粒子i的坐标、方向、以及粒子权重；(x_t,y_t,θ_t)则组成了估计后的预测道路点Cpoint(t)的状态。在本实施例中，ThGrayL＝0.125，ThGrayH＝0.325，σ＝0.3，α＝1，β＝1。

S5.1.7、采用经典粒子滤波重采样方法，用当前粒子集PC对预测粒子集PP进行重采样，完成对当前粒子集PC的更新，并返回步骤S5.1.4；

在本实施例中，当前粒子集PC对预测粒子集PP进行重采样的步骤为：

1)、对预测粒子集PP根据权重值的大小进行粒子排序，排序方向由大到小，并清零当前粒子集PC；

2)、在预测粒子集PP中，对每一个预测粒子PP(i)计算再取整得到其中i＝1,....,n；

3)、由预测粒子集PP到当前粒子集PC的复制循环过程

定义计数变量k＝1，k'＝k，同时令i＝1；

(a)、令如果且k′＜n+1，则将PP(i)复制次至当前粒子集PC中，对应被赋值粒子为PC(k)～PC(k′-1)，进入步骤b)；如果且k′＞n+1，则将PP(i)复制(n-k+1)次至当前粒子集PC中，对应被赋值粒子为PC(k)～PC(n)，进入步骤4)；

(b)、令k＝k′,i＝i+1，对i进行判断，如果i＜n，则返回步骤(a)；如果i＞n，将未被赋值的PC(k)～PC(n)均赋值为PP(1)，进入步骤4)；

4)、将当前粒子集PC中的所有粒子的粒子权重设置为1/n，完成当前粒子集PC的更新；

S5.1.8、输出道路中心点空间Cpoint作为粒子滤波结果，如图7(a)所示，粒子滤波结果描述了道路的大概位置和方向，但在细节上与真实的道路位置存在一定的差异，平滑度欠佳；

S5.2、利用SNAKE方法对道路中心线的粒子滤波结果进行校正，如图7(b)所示，校正后道路中心线的平滑性更好，与道路真实位置的符合程度也更高。

在本实施例中，SNAKE方法对粒子滤波结果进行校正如下：

1)、构建道路中心线的轮廓控制点序列空间V＝{v(s),s∈[1,m]}，其中，v(s)＝{(x_vs,y_vs),(x_vs,y_vs)∈cpoint(s)}；构建能量矩阵序列E＝{E(sj),sj∈[1,8]}，其初始值均为0。定义迭代次数t′＝1；进入步骤2)；

2)、对空间V中的点从v(1)开始到v(m)依次进行能量遍历计算，并根据逐点的能量遍历结果对空间V的各点根据序号从v(1)到v(m)依次进行更新，其遍历计算过程如下：

以当前遍历点v(s)为中心，计算其八邻域内每个像元点的三个能量参数E_elastic,E_bending,E_ext，并得到该点的能量值E(sj)，设邻域像元点集合为Point8＝{point8(sj),sj∈[1,8]}，其中point8(sj)＝(x_psj,y_psj)，对Point8(sj)的能量参数和能量值的计算方法如下：

E_elastic＝α′|Point8(sj)-v_s-1|²

Ebending＝β′|v_s+1-2Point8(sj)+v_s-1|²

Eext＝γ(Point8(sj))

E(sj)＝E_elastic+Ebending-Eext

在本实施例中，应力参数α′＝0.5，刚度β′＝0.75，γ(Point8(sj))为边缘信息图像中点Point8(sj)对应的边缘强度值。特别地，当s＝1时，设v_s-1＝v_s；

待八邻域所有像元点计算完成后，对E中的元素进行排序，找出能量最小值E(sj′)＝min(E(sj))，并选择Point8(sj′)对空间V中的点v(s)进行更新，如下所示：

v(s)＝Point8(sj′)

更新完成后，进入步骤3)；

3)、令t′＝t′+1，若t′≤T_end，则返回步骤2)；若t′＞T_end，则输出道路中心线的轮廓控制点空间V＝{v₁,v₂,....,v_m}作为校正后的道路中心线；在本实施例中，T_end＝20。

S6、提取道路损毁信息

S6.1、获取道路中心线上各点的边缘损毁指数

以道路中心线为基准，对道路中心线上的像元点依次进行遍历，得到道路中心线上各点的边缘损毁指数；下面对具体的获取方法进行详细说明，如下：

(T6.1)、确定待检测SAR图像中道路中心线上的像元点

设待检测SAR图像中道路中心线上的像元点集合为其中， $\mathrm{ROAD}\left(\hat{i}\right)=\{x_{\hat{i}},y_{\hat{i}}\},$ 表示道路中心线上第个像元点在待检测SAR图像中的坐标，为道路中心线上像元点的总数；以道路中心线为基准，对道路中心线上的像元点依次进行遍历：

如果当前待检测点不为道路中心线的最后一个像元点，则进入步骤(T6.2)；如果当前待检测点为待检测SAR图像中道路中心线的最后一个像元点，则进入步骤(T6.7)；

(T6.2)、以当前待检测点Q为中心，选取长度为L的道路中心点序列S；在本实施例中，L＝9。

(T6.3)、计算当前待检测点Q的水平区域均值比γ_Q

(T6.3.1)、将序列S在水平方向上分别向左和向右平移(3*W-1)，在本实施例中，W＝3；

从而构成了长度为L＝9、宽度为W＝3的三块窗口区域，即左背景区域、中心区域和右背景区域；

(T6.3.2)、计算当前待检测点Q的水平区域均值比γ_Q

${\mathrm{\&gamma;}}_Q=\min (\frac{u_{\mathrm{lQ}}}{u_{\mathrm{cQ}}},\frac{u_{\mathrm{cQ}}}{u_{\mathrm{lQ}}},\frac{u_{\mathrm{rQ}}}{u_{\mathrm{cQ}}},\frac{u_{\mathrm{cQ}}}{u_{\mathrm{rQ}}})$

其中，u_lQ,u_cQ,u_rQ分别代表左背景区域、中心区域和右背景区域内的区域灰度均值；

(T6.4)、计算当前待检测点Q的垂直区域均值比

(T6.4.1)、将序列S在垂直方向上分别向上和向下平移(3*W-1)，构成长度为L＝9、宽度为W＝3的三块窗口区域，即上背景区域、中心区域和下背景区域；

(T6.4.2)、计算当前待检测点Q的垂直区域均值比

${\mathrm{\&gamma;}}_Q^{*}=\min (\frac{u_{\mathrm{uQ}}}{u_{\mathrm{cQ}}},\frac{u_{\mathrm{cQ}}}{u_{\mathrm{uQ}}},\frac{u_{\mathrm{dQ}}}{u_{\mathrm{cQ}}},\frac{u_{\mathrm{cQ}}}{u_{\mathrm{dQ}}})$

其中，u_uQ,u_cQ,u_dQ分别代表上背景区域、中心区域和下背景区域内的区域灰度均值；

(T6.5)、计算当前待检测点Q的边缘损毁参数I_Q

$I_Q=\sqrt{({\mathrm{\&gamma;}}_Q^2+{\left({\mathrm{\&gamma;}}_Q^{*}\right)}^2)}$

(T6.6)、根据当前待检测点Q的边缘损毁参数I_Q得到该点处的边缘损毁指数I_eQ；

其中，Te为分割道路损毁区域的高边缘阈值；在本实施例中，Te＝0.8。

(T6.7)、当时，则退出遍历，输出边缘损毁指数集合

S6.2、获取道路中心线上各点的灰度损毁指数

以道路中心线为基准，对道路中心线上的像元点依次进行逐点遍历：

如果当前待检测点不为道路中心线的最后一个像元点，则进行如下计算，得到当前待检测点Q的灰度损毁指数I_gQ；

其中，Tg和Te'分别是灰度损毁特征阈值和低边缘损毁特征阈值，G(Q)为当前待检测点Q在待检测SAR图像中的邻域灰度均值；在本实施例中Tg＝0.35,Te′＝0.3；

如果当前待检测点为待检测SAR图像中道路中心线的最后一个像元点，即时，则退出遍历，输出灰度损毁指数集合

S6.3、将道路中心线上各点的边缘损毁指数和灰度损毁指数求并集I_e∪I_g，可以得到道路损毁结果。

在本实施例中，如图8所示，白色点或线段标注了待检测SAR图像中道路的损毁段。

尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

Claims (6)

1.一种基于SAR图像的道路损毁信息提取方法，其特征在于，包括以下步骤：

(1)、获取双时相原始SAR图像

在同一监测点，分别获取道路损毁发生前、后两幅双时相原始SAR图像，标记为基准原始SAR图像和待检测原始SAR图像；

(2)、对双时相原始SAR图像进行预处理

对双时相原始SAR图像进行图像配准，得到配准后的基准SAR图像和待检测SAR图像，再对基准SAR图像进行增强处理，得到增强图像；

(3)、边缘信息提取

通过双边缘检测算子Duda对增强图像进行边缘信息检测，从该增强图像中提取出边缘信息图像；

(4)、边缘方向信息获取

对边缘信息图像进行方向离散值计算，获取到边缘信息图像中逐像元点的方向信息，从而得到边缘方向信息图像；

(5)、道路中心线提取

根据增强图像、边缘信息图像以及边缘方向信息图像中的特征，采用粒子滤波方法和SNAKE方法提出道路的中心线；

(6)、提取道路损毁信息

以步骤(5)中提取的中心线为基础，获取道路中心线上各点的边缘损毁指数和灰度损毁指数，再将边缘损毁指数和灰度损毁指数求并集，得到道路的损毁结果。

2.根据权利要求1所述的基于SAR图像的道路损毁信息提取方法，其特征在于，所述的步骤(4)中，边缘方向信息图像的获取方法为：

(2.1)、对边缘信息图像通过预定的阈值进行二值化处理，得到二值分割图像；

(2.2)、采用八方向计算模板逐点遍历二值分割图像中每个像元点的边缘方向，即计算每个像元点的方向离散值；

其中，Direction为边缘方向信息图；I为二值分割图像；d_k为方向k上的方向模板窗口，k∈[1,8]；为方向信息计算时的当前像元点坐标，N为边缘信息图像中像元点的个数。

3.根据权利要求1所述的基于SAR图像的道路损毁信息提取方法，其特征在于，所述步骤(5)中，道路中心线提取的方法为：

(3.1)、利用粒子滤波方法提取道路的中心线

(3.1.1)、定义迭代次数为t，t∈[0,m]；构建道路中心点空间 $\mathrm{Cpoint}=\{\mathrm{Cpoint}\left(\tilde{i}\right),\tilde{i}\&Element;[1,m\left]\right\},$ 其中， $\mathrm{Cpoint}\left(\tilde{i}\right)=(x_{\tilde{i}},y_{\tilde{i}},{\mathrm{\&theta;}}_{\tilde{i}}),$ m为道路中心线上像元点的个个数；

(3.1.2)、在基准SAR图像中的道路区域任意选取一个像元点point₀，将作为道路种子点，分别构建两个大小为n的粒子集空间，即当前粒子集PC＝{pc(i)}，预测粒子集PP＝{pp(i)}，其中，i＝1,2,…,n， {xⁱ,yⁱ}代表序号为i的粒子的位置坐标，θⁱ为序号为i的粒子的方向，wⁱ为两粒子集中各粒子的粒子权重，且初始粒子权重设置为1；

(3.1.3)、对种子点point进行采样，得到n个粒子元素，并将这n个粒子元素作为当前粒子集PC中的元素，同时将初始化迭代次数t设置为0，再进入步骤(3.1.4)；

(3.1.4)、如果t＜m，则令t＝t+1，再进入步骤(3.1.5)；如果t＞m-1，则停止迭代，跳转至步骤(3.1.8)；

(3.1.5)、更新预测粒子集PP

通过状态转移方程计算第t次迭代时当前粒子集PC中每个粒子pc(i)的预测状态，并将该结果赋值给对应的预测粒子PP(i)；

$\left\{\begin{array}{c}{x}_t=x_{t-1}+d\cos \left({\mathrm{\&theta;}}_{t-1}\right)+u_1\\ y_t=y_{t-1}+d\sin \left({\mathrm{\&theta;}}_{t-1}\right)+u_2\\ {\mathrm{\&theta;}}_t={\mathrm{\&theta;}}_{t-1}+u_3\end{array}\right.$

其中，表示第t次迭代时预测粒子PP(i)的坐标，表示第t次迭代时预测粒子PP(i)的方向，d为预设的步长间隔，u₁,u₂,u₃分别为用于修正状态转移方程的随机噪声；

(3.1.6)、对预测粒子集PP中逐粒子进行粒子权重计算，并进行粒子权重归一化，再利用得到的归一化权重估计预测道路点的状态；

(3.1.7)、采用经典粒子滤波重采样方法，用当前粒子集PC对预测粒子集PP进行重采样，完成对当前粒子集PC的更新，并返回步骤(3.1.4)；

(3.1.8)、输出道路中心点空间Cpoint作为粒子滤波结果；

(3.2)、利用SNAKE方法对道路中心线的粒子滤波结果进行校正，得到准确的道路中心线。

4.根据权利要求3所述的基于SAR图像的道路损毁信息提取方法，其特征在于，所述的步骤(3.1.6)中，对预测粒子集PP中逐粒子进行粒子权重计算和权重归一化的方法为：

权重计算：

(4.1)、获取预测粒子在第t次迭代时的灰度值边缘强度值r_t ⁱ和边缘方向信息离散值

设第i个预测粒子的坐标为(X_i,Y_i)，通过该坐标值在基准增强图像中查找出该预测粒子在第t次迭代时的灰度值在边缘信息图像中查找出该预测粒子的边缘强度值r_t ⁱ，在边缘方向信息图像中查找出该预测粒子的边缘方向信息离散值

(4.2)、计算预测粒子的最终权重

先根据求出第i个预测粒子的灰度权重边缘权重方向权重

$w_g^i=\left\{\begin{array}{cc}1/|g_{t-1}-g_t^i|& \mathrm{ThGrayL}<g_t^i<\mathrm{ThGrayH}\\ 0& \mathrm{Others}\end{array}\right.$

$w_e^i=\frac{1}{\sqrt{2\mathrm{\&pi;}}\mathrm{\&sigma;}}{e}^{-\frac{{(1-r_t^i)}^2}{2{\mathrm{\&sigma;}}^2}}$

$w_d^i=1/|{\mathrm{\&theta;}}_{t-1}-{\mathrm{\&theta;}}_t^i|$

再计算出预测粒子的最终权重

$w_t^i=w_e^i+\mathrm{\&alpha;}{w}_g^i+\mathrm{\&beta;}{w}_d^i$

其中，ThGrayL、ThGrayH分别为预设的低灰度阈值和高灰度阈值，σ为数，α、β为权重比例调整参数；

权重归一化：

$\mathrm{SumWeight}=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{w}_t^i$

$w_t^i=w_t^i/\mathrm{SumWeight}$

$\left\{\begin{array}{c}{x}_t=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_t^i{w}_t^i/\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{w}_t^i\\ y_t=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{y}_t^i{w}_t^i/\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{w}_t^i\\ {\mathrm{\&theta;}}_t=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&theta;}}_t^i{w}_t^i/\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{w}_t^i\end{array}\right.$

其中，表示预测粒子集中第i个预测粒子的状态，其中， 分别表示第t次迭代时预测粒子i的坐标、方向、以及粒子权重；(x_t,y_t,θ_t)则组成了估计后的预测道路点Cpoint(t)的状态。

5.根据权利要求1所述的基于SAR图像的道路损毁信息提取方法，其特征在于，所述的步骤(6)中，提取道路损毁信息的方法为：

(5.1)、获取道路中心线上各点的边缘损毁指数

以道路中心线为基准，对待检测SAR图像中道路中心线上的像元点依次进行遍历，得到道路中心线上各点的边缘损毁指数；

(5.2)、获取道路中心线上各点的灰度损毁指数

以道路中心线为基准，对待检测SAR图像中道路中心线上的像元点依次进行遍历，得到道路中心线上各点的灰度损毁指数；

(5.3)、将道路中心线上各点的边缘损毁指数和灰度损毁指数求并集，得到道路损毁结果。

6.根据权利要求5所述的基于SAR图像的道路损毁信息提取方法，其特征在于，所述的步骤(5.1)中，道路中心线上各点的边缘损毁指数获取方法为：

(6.1)、确定待检测SAR图像中道路中心线上的像元点

设待检测SAR图像中道路中心线上的像元点集合为 $\mathrm{ROAD}=\{\mathrm{ROAD}\left(\hat{i}\right),\hat{i}\&Element;[1,\hat{n}\left]\right\},$ 其中， $\mathrm{ROAD}\left(\hat{i}\right)=\{x_{\hat{i}},y_{\hat{i}}\},$ 表示道路中心线上第个像元点在待检测SAR图像中的坐标，为道路中心线上像元点的总数；以道路中心线为基准，对道路中心线上的像元点依次进行遍历：

如果当前待检测点不为道路中心线的最后一个像元点，则进入步骤(6.2)；如果当前待检测点为待检测SAR图像中道路中心线的最后一个像元点，则进入步骤(6.7)；

(6.2)、以当前待检测点Q为中心，选取长度为L的道路中心点序列S；

(6.3)、计算当前待检测点Q的水平区域均值比γ_Q

(6.3.1)、将序列S在水平方向上分别向左和向右平移(3*W-1)，构成长度为L、宽度为W的三块窗口区域，即左背景区域、中心区域和右背景区域；

(6.3.2)、计算当前待检测点Q的水平区域均值比γ_Q

${\mathrm{\&gamma;}}_Q=\min (\frac{u_{\mathrm{lQ}}}{u_{\mathrm{cQ}}},\frac{u_{\mathrm{cQ}}}{u_{\mathrm{lQ}}},\frac{u_{\mathrm{rQ}}}{u_{\mathrm{cQ}}},\frac{u_{\mathrm{cQ}}}{u_{\mathrm{rQ}}})$

其中，u_lQ,u_cQ,u_rQ分别代表左背景区域、中心区域和右背景区域内的区域灰度均值；

(6.4)、计算当前待检测点Q的垂直区域均值比

(6.4.1)、将序列S在垂直方向上分别向上和向下平移(3*W-1)，构成长度为L、宽度为W的三块窗口区域，即上背景区域、中心区域和下背景区域；

(6.4.2)、计算当前待检测点Q的垂直区域均值比

${\mathrm{\&gamma;}}_Q^{*}=\min (\frac{u_{\mathrm{uQ}}}{u_{\mathrm{cQ}}},\frac{u_{\mathrm{cQ}}}{u_{\mathrm{uQ}}},\frac{u_{\mathrm{dQ}}}{u_{\mathrm{cQ}}},\frac{u_{\mathrm{cQ}}}{u_{\mathrm{dQ}}})$

其中，u_uQ,u_cQ,u_dQ分别代表上背景区域、中心区域和下背景区域内的区域灰度均值；

(6.5)、计算当前待检测点Q的边缘损毁参数I_Q

$I_Q=\sqrt{({\mathrm{\&gamma;}}_Q^2+{\left({\mathrm{\&gamma;}}_Q^{*}\right)}^2)}$

(6.6)、根据当前待检测点Q的边缘损毁参数I_Q得到该点处的边缘损毁指数I_eQ；

其中，Te为分割道路损毁区域的高边缘阈值；

(6.7)、当时，则退出遍历，输出边缘损毁指数集合I_e，

所述的步骤(5.2)中，道路中心线上各点的灰度损毁指数获取方法为：

对道路中心线上的像元点依次进行逐点遍历：

如果当前待检测点不为道路中心线的最后一个像元点，则进行如下计算，得到当前待检测点Q的灰度损毁指数I_gQ；

其中，Tg和Te′分别是灰度损毁特征阈值和低边缘损毁特征阈值，G(Q)为当前待检测点Q在待检测SAR图像中的邻域灰度均值；

如果当前待检测点为待检测SAR图像中道路中心线的最后一个像元点，即时，则退出遍历，输出灰度损毁指数集合I_g，