CN106709465A - 基于条件随机场的极化sar图像道路提取方法 - Google Patents

Abstract

本发明提供一种基于条件随机场的极化SAR图像道路提取方法，构建一个条件随机场框架实现道路提取，包括以利用多尺度线状目标检测算子MLFD对输入SAR图像构建金字塔，获得多尺度的图像金字塔；通过CRF进行建模，引入上下文信息，其中单元势函数通过logistic函数建模，根据成对势函数设定采用Beamlet分解方式寻找最优划分；采用相关约束来对道路的基元进行编组和标记，描述成对势函数，合并从全局数据中得到的信息，对特征矩阵进行统一的归一化处理。应用本发明技术方案进行道路提取的效率高，精度好，适于推广使用。

Description

基于条件随机场的极化SAR图像道路提取方法

技术领域

本发明属于图像处理技术领域，特别涉及一种基于条件随机场的极化SAR图像道路提取方法。

背景技术

合成孔径雷达(Synthetic Aperture Radar,SAR)图像具有极强的现实应用，如测绘，遥感，城市规划，农业和灾害预防。与原来的单一极化的SAR相比，极化SAR显示出在数据的应用程序上的能力，并且能够获得更丰富的目标信息，以及对完全偏振光散射机制的鉴定。它采用主动成像原理，实现全天候对目标的连续跟踪，不受天气和照明的影响。其中SAR图像解译应用，道路提取具有十分重要的研究意义，因为线性指标(包括道路，桥梁，脊线，沿海线)由于奇散射将在SAR图像中呈现为较暗的区域。

从SAR图像进行道路提取一般分为两步：局部道路候选段检测和全局道路网络优化。对于局部路段检测，可用不同的检测算子从局部像素获得候选段，如传统的边缘检测、形态学算子。对于全局道路网优化，可由先验信息约束提供全局选择技术。如Tupin等人提出通过D1和D2操作的两步技术提取局部候选段，然后在局部候选段建立马尔可夫随机场(MRF)实现全局道路的网络优化。

传统的基于线性检测算子的方法，是通过SAR图像的局部特征来获得道路的片段。由于SAR图像的斑点噪声，需通过全局优化以改进局部结果。候选道路段的全局选择本质上是利用上下文信息类别的标识过程，即每条道路的候选段标记为道路或噪声。马尔科夫随机场(MRF)是利用上下文信息，最大限度地降低全局损失函数，这是一种对分类结果产生积极影响的方法。然而，MRF需要估计标签和数据的联合分布，并涉及数据的分布，因此总是难以实现。

发明内容

本发明的目的在于针对极化SAR图像道路提取问题，提出一种基于条件随机场模型的方法，实现了多尺度的道路提取，并在统一的框架内进行推断学习。

本发明的技术方案提供一种基于条件随机场的极化SAR图像道路提取方法，构建一个条件随机场框架实现道路提取，包括以下步骤，

步骤1，利用多尺度线状目标检测算子MLFD对输入SAR图像构建金字塔，获得多尺度的图像金字塔；

步骤2，通过CRF基于步骤1所得结果进行建模，引入上下文信息，实现如下，

设给定图像I，x表示观测数据，y表示标号集合，记I＝{x₁,x₂,...,x_M}；设i表示图像所有块集合S＝{1,2,...,M}中的一个子块，观测数据对应为x_i；设j表示块i的邻域集合N_i中的一块，给定数据x时标号y的后验概率由指数形式给出如下，

式中，φ(x)表示对观测数据x进行特征映射，A_i(.)和I_ij(.)分别是单元势函数和成对势函数，A_i(.)表示给定观测数据x_i时相应块i被标记为y_i的期望，I_ij(.)描述块i和块j之间的相互影响；φ_i(x)将块i的观测数据x_i映射到特征空间向量；μ(φ_i(x),φ_j(x))是成对块(i,j)的新特征向量；Z(φ(x))是归一化常数，表示如下，

步骤3，单元势函数A_i(y_i,φ_i(x))通过logistic函数来建模，表示如下，

式中，h_i(x)表示第i个二进方块的特征矩阵，由MLFD的响应值构成，T_n(m_i)表示该方块在尺度n下的MLFD最大响应值；w＝[w₁,w₂,…,w_K,α₁]^T是模型参数，包含特征h_i(x)的各项权重；

步骤4，根据成对势函数I_ij(y_i,y_j,μ(φ_i(x),φ_j(x)))，设定采用Beamlet分解方式寻找最优划分；

步骤5，采用相关约束来对道路的基元进行编组和标记，描述成对势函数，合并从全局数据x中得到的信息，表达如下式，

I_ij(y_i,y_j,μ(φ_i(x),φ_j(x)))＝y_iy_jv^Tμ_ij(x) (1)

式中，v＝[v₁,v₂,…,v_m]^T表示特征矩阵μ_ij(x)中每维特征的权重参数，m表示联合特征μ_ij(x)的维度，μ_ij(x)是块i和块j的联合特征；

步骤6，对特征矩阵进行统一的归一化处理。

而且，步骤5中，所述相关约束包括夹角距离、横向距离、端点距离、临近度和曲率连接。

而且，步骤6中，归一化处理采用阈值局部化方式。

而且，采用L-BFGS优化算法训练单元势特征的参数w和成对势特征的参数v。

而且，所述多尺度线状目标检测算子MLFD，将图像迭代地分割成2×2块，构建图像金字塔，在每一层上提取图像特征，合并不同层次的特征用于模型优化。

与现有技术相比，本发明的有益效果是：本发明构建了一个条件随机场框架实现道路提取，可实现先验与似然之间的交互，同时可实现相互学习互动和事先可能性信息相互学习。在此框架中：1)多尺度线状目标检测算子(MLFD)提供道路段中有效的似然信息来获取道路结构的稀疏表示；2)局部相邻段连接用若干个几何先验关系来描述；3)图像金字塔的构建引入多尺度信息，使单元势和成对势可在统一框架中进行推断和优化。应用本发明技术方案进行道路提取的效率高，精度好，适于推广使用。

附图说明

图1是本发明实施例的原理图；

图2是本发明实施例中步骤1的示意图，其中图2(a)为构建金字塔示意图、图2(b)为MLFD算法掩膜示意图，图2(c)为多尺度应用示意图；

图3是本发明实施例步骤4的5种相关约束示意图，其中图3(a)为夹角距离约束示意图，图3(b)为横向距离约束示意图，图3(c)为端点距离约束示意图，图3(d)为临近度约束示意图，图3(e)为曲率连接约束示意图。

具体实施方式

本发明针对极化SAR图像，结合极化特征、多分辨率信息以及统计特性，提出一种基于条件随机场的极化SAR图像道路提取方法。

以下根据附图和实施例详细说明本发明的技术方案。

条件随机场(CRF)模型是判别图模型，该模型训练和推理过程使用整个图像现有标签可能性和后验概率分布，具有实现精确和稳健的标记结果的优势，本发明基于条件随机场模型实现道路检测。本发明实施例提供的技术方案原理如图1所示，包括输入SAR图像，经多尺度检测算子得到多尺度的图像金字塔，进行单元势提取得到MLFD框架及单元势中的似然；然后进行成对势提取，基于子束分析、先验约束和成对势，构建模型，获得条件随机场模型，基于模型提取道路，获取道路提取结果。

以下分步说明具体实现方式：

步骤1利用多尺度线状目标检测算子MLFD对输入SAR图像构建金字塔，获得多尺度的图像金字塔。

优选地，基于Ratio相关检测器(Ratio表示比例)和Wedgelet多尺度分析(Wedgelet表示楔波变换)，使用多个掩膜滤波器将图像分成不同的匀质区域来检测局部线条片段，利用多分辨率方法扩展掩膜宽度的选择范围，并自适应改变掩膜大小和中心区域的宽度与方向，综合考虑局域上下文的统计和几何属性，如图2(c)所示。

MLFD将图像迭代地分割成2×2块，构建图像金字塔，如图2(a)所示。在每一层上提取图像特征，合并不同层次的特征用于模型优化。其中特征提取细节是将局部图像块分成三块区域，如图2(b)所示，包括区域1、2、3，掩膜表示为(v₁,v₂,w,μ₁,μ₂,μ₃)，其中v₁,v₂是中间区域的中心线端点，w是中间区域的宽度，μ₁,μ₂和μ₃分别表示3个不同区域的均值，即各区域超像素块数据的平均值region1、region2、region3，由下式计算得到

式(2)中，p为区域i的像素点，i＝1,2,3；A_p是像素点p的幅值，n_i为区域i内的像素点个数。用l表示中间区域的长度，η表示统一性系数，则MLFD的响应值为

式(3)中，m表示待检测区域的模板中心线，r和ρ分别是D1算子和D2算子的响应值，计算公式由下式给出

式(4)中，参数c_ij＝μ_i/μ_j，γ_i和γ_j分别表示区域i和j的方差系数。参数η用于评估中间区域的连续性。此时，i＝1，j＝2,3；即中间参数r_ij包括r₁₂,r₁₃，中间参数ρ_ij包括ρ₁₂,ρ₁₃。

MLFD线状目标检测算法是通过将图像分解成不同尺度且不重叠的二进区域完成的，分解后的二进区域集合C满足下式

式(5)中，m:s表示检测模板m在二进区域s上，#C表示集合C中的对象数目，λ为复杂性惩罚系数，表征二进区域的分解精度。采用树裁剪方法求解上式，将图像分解为一系列边长相等的二进方块，尺度n的方块被分解为2×2的尺度为n+1方块，直到最大尺度n_max，最终形成一颗完全四叉树。例如用一个2×2的窗，以步长为1进行滑动，那么n_max＝n-2+1。所有子块使用式(5)计算响应值

式(6)中，T_n-1(m_p)表示尺度为n-1的某一方块的最大响应值，T_n(m_i)表示该方块的4个子块中第i个子块的最大响应值，其尺度为n。保存最大值以及对应的检测模板，得到道路基元若上式成立，则保留四个子块的响应值和对应的道路基元，否则删除子块，保留父二进块，并更新其响应值为四个子块的最大响应值之和。重复裁剪直到根，则保留下来的子块满足上式。

步骤2与所有的图模型一样，CRF能根据提供的样本和样本特征等数据给分类结果提供精确决策，并为标号域计算出分配概率。CRF同时结合描述模型和判别模型，构建描述样本之间关系的判别函数，只对在给定样本数据时的标号场进行后验概率建模。

经步骤1对原图像进行分割后，步骤2对步骤1中的结果进行建模，即通过CRF引入上下文信息。

给定图像I，x表示观测数据，y表示标号集合，则I＝{x₁,x₂,...,x_M}，即图像由M个超像素块的观测数据x_i(一小部分像素点的观测数据集合)构成，x＝1,2,...,M；为便于标记起见，设i表示图像所有块集合S＝{1,2,...,M}中的一个子块(超像素块)，其观测数据对应为x_i；j表示块i的邻域集合N_i中的一个块。给定数据x时标号y的后验概率可由指数形式给出

式(7)中，φ(x)表示对观测数据x进行特征映射，即对x进行特征提取；A_i(.)和I_ij(.)分别是单元势函数和成对势函数，A_i(.)表示给定观测数据x_i时块i被标记为y_i的期望，I_ij(.)描述块i和块j之间的相互影响。φ_i(x)将块i的观测数据x_i映射到特征空间向量；μ(φ_i(x),φ_j(x))是成对块(i,j)的新特征向量。Z(φ(x))是归一化常数，可表示如下，

CRF用于建模分类的条件概率即P(y|φ(x))，该条件概率中的参数通过训练样本学习得到；得到该条件概率模型后，输入需要分类的样本，分类标签为使得该条件概率最大的标签。

步骤3本发明中单元势函数A_i(y_i,φ_i(x))通过logistic函数来建模，表示如下

式(9)中_，h_i(x)表示第i个二进方块的特征矩阵，由多尺度MLFD的响应值构成，T_n(m_i)表示该方块在尺度n下的MLFD最大响应值。w＝[w₁,w₂,…,w_K,a₁]^T是模型参数，包含特征h_i(x)的各项权重。

步骤4因为需要首先将输入图像进行划分成方块，然后将方块间的相关信息通过CRF建模，本发明针对成对势函数的建立，提出方块的划分是通过Beamlet实现的：

成对势函数I_ij(y_i,y_j,μ(φ_i(x),φ_j(x)))决策块i和块j之间的空间相互作用，影响的强度取决于观测数据x。在已有的图模型研究中，统计物理学中的Ising或Potts模型在区域建模中得到广泛应用，通过参数调整可实现对局部属性的较好描述。

在CRF道路提取框架中，成对势函数所需的先验信息使用子束波分析框架获取，并引入该框架下局部邻域内道路基元之间的多种约束信息，充分利用CRF灵活设计成对势函数特征矩阵的特点，块与块之间的局部上下文关系利用道路基元的长度、曲率和交叉距离等先验信息来设计成对势函数，以用于本发明对道路和噪声的二分类任务。

子束波(Beamlet)是一种多分辨率图像分析框架，在此框架下，线段扮演的角色类似于点在小波分析中的扮演的角色。Beamlet链提供对平面上曲线的最优稀疏逼近表示，Beamlet图结构以图像的像素为顶点，连接每一对像素作为边，其基于邻域Beamlet系数寻找最优的全局优化算法。

Beamlet基元基于二进方块，连接方块边界上的任意两点形成，边界的长度则决定基元的尺度。产生过程如下：首先对样本图像块进行二进剖分，在尺度0≤n≤N的二进方形集合中，对于K ≥0，固定一个分辨率量δ＝2^-N-K，即将每个正方形一直分割到最小尺度；然后在每个二进方块中，从右上角开始顺时针沿边界旋转，相隔距离δ标记出等间隔的定点，则在边长为2^-n的方块中，存在M_n＝4·2^K·2^N-n个顶点；最后在正方形中连接任意一对顶点，即形成一条Beamlet基元。

离散Beamlet变换的响应由下式给出

式(10)中，g(p)是像素p的灰度值，b为Beamlet基元，l(b)为该基元的长度，即基元上像素点的个数。给定一幅图像，将图像划分为一系列二进方块，每个方块用一个Beamlet基表示，Beamlet分解则是寻找最优划分

式(11)中，P表示所有可能的划分，λ为复杂性惩罚系数，#p表示该划分中二进方块的数目，b～s表示beamlet基b在方块s上，λ决定分割的精细程度，值越小能得到越多的细节，较大λ则可降低噪声，获取提取目标的轮廓。

步骤5SAR图像路网提取中，道路基元之间的上下文联系，如交叉角度距离、曲率、端点距离、临近度等参数，为全局优化提供有效的先验信息。本发明实施例使用了如下的几种相关约束来对道路的基元进行编组和标记。

(1)夹角距离：该关系检测两条基元的交叉夹角(可能是线段的延长线)，如图3(a)所示，由下式定义

D^a＝min(θ,π-θ) (12)

式(12)中，θ₁和θ₂是两条基元的切线角，可由线条的两端点计算得到，取值范围在(-π/2,π/2]。

(2)横向距离：道路的局部区域内的基元线条应该有较近的横向距离，该距离测量短线条的中点到长线条的垂直距离。横向距离如图3(b)所示，可表示为

D^t＝D/L

(13)

式(13)中，L是长线条的长度，D是短线条的中点到长线条的垂直距离。

(3)端点距离：与横向距离类似，道路基元应与其邻域内基元之间有较小的端点距离。设d是两条线条的端点之间最小的交叉距离，L是长线条的长度，如图3(c)所示，则端点距离可写成

D^e＝d/L

(14)

(4)临近度：两条基元线条的端点之间的临近度反映基元之间的感知灵敏度，决策毗邻基元应该邻接的程度。计算示意图如图3(d)，可用式子表示为

P＝L²/2πR²

(15)

式(15)中，L是两条道路片段中长度最小者，R是两条线条的端点之间最小的交叉距离，类似端点距离计算公式中的距离d。

(5)曲率连接：曲率连接描述片段之间的结构信息，该关系决定道路片段被连接成光滑曲线的权重。如图3(e)，可用式子表示为

式(16)中，α和β分别是连接片段在连接点处与两片段的夹角，ε₁和ε₂分别是连接片段的平滑度的权重和对连接点之间的距离R的权重。

本发明使用泛化的Ising模型来描述成对势函数，合并从全局数据x中得到的信息，可写成下式

I_ij(y_i,y_j,μ(φ_i(x),φ_j(x)))＝y_iy_jv^Tμ_ij(x) (17)

式(17)中，v＝[v₁,v₂,…,v_m]^T表示特征矩阵μ_ij(x)中每维特征的权重，m表示联合特征μ_ij(x)的维度，μ_ij(x)是块i和块j的联合特征。

步骤6已提取道路的不同特征和约束信息，这些信息获取道路的不同属性，但对道路的特性具有统一的描述。为了避免CRF在训练过程中出现数值不平稳，需要对特征矩阵进行统一的归一化处理。一般而言，归一化过程需要选择阈值，该值可能被设置为全局常数，但是更为高效的算法是阈值局部化，在测试点的邻域内寻找最佳阈值参数，局部化阈值使得特征与邻域相似的点被标记和保留。

本发明实施例中对夹角距离D^a、横向距离D^t和端点距离D^e进行二值化处理，每个二值化过程均与对应的局部邻域均值(4连接邻域)比较。块i与其邻域块内道路基元之间的夹角距离平均值由下式计算

式(18)中，N_i表示块i的邻域块集合，块k是N_i中的块，是块i和块k内道路基元的夹角距离，num(N_i)是集合N_i内块的数目。

考虑到道路的延展平滑性，局部块间道路基元的夹角距离D^a、横向距离D^t和端点距离D^e应该具有较小值，因而具体实施时可设置由实验得到的经验阈值T_a，如果块i与块j的夹角距离与局部平均值的绝对差小于T_a，则保留该夹角信息，否则删除，该夹角距离可表示为

根据经验阈值T_t和T_e，可求得二值化的横向距离和端点距离

对临近度P_ij和曲率连接C_ij，由于图像中存在如立交桥等曲线道路，而曲线和直线的曲率变化存在较大差异，因而本发明对P_ij和C_ij的后处理，使用邻域平均值对两者进行归一化，而不进行二值化处理，如下式

这里给出成对势函数式中的联合特征μ_ij(x)，其包含块i和块j的特征矩阵h_i(x)和h_j(x)，以及步骤5提到的5种先验约束信息，则μ_ij(x)可表示为

式(22)中，0<n≤N，N的大小指方块边界的长度，该长度决定了基元的尺度；x_nij＝f(h_ni(x),h_nj(x))，其中h_ni(x)和h_nj(x)分别是块i和块j的第n维特征，x_nij是h_ni(x)和h_nj(x)的函数关系f(.,.)，由两者的比率决定，同时设置一个固定的由实验得到的b，将x_nij归一化到[0,1]，可表示如下

式(23)中，参数r₁,r₂可表示为

本发明使用现有L-BFGS优化算法来训练CRF模型中单元势特征参数w和成对势特征参数v，推理过程则由现有迭代信度传递LBP算法来测试样本的标号。具体实施时可采用计算机软件技术实现自动道路提取。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (6)

1.一种基于条件随机场的极化SAR图像道路提取方法，其特征在于：构建一个条件随机场框架实现道路提取，包括以下步骤，

步骤1，利用多尺度线状目标检测算子MLFD对输入SAR图像构建金字塔，获得多尺度的图像金字塔；

步骤2，通过CRF基于步骤1所得结果进行建模，引入上下文信息，实现如下，

设给定图像I，x表示观测数据，y表示标号集合，记I＝{x₁,x₂,...,x_M}；设i表示图像所有块集合S＝{1,2,...,M}中的一个子块，观测数据对应为x_i；设j表示块i的邻域集合N_i中的一块，给定数据x时标号y的后验概率由指数形式给出如下，

$P\left(y\right|\mathrm{\&phi;}\left(x\right))=\frac{1}{Z\left(\mathrm{\&phi;}\right(x\left)\right)}\exp \left\{\underset{i\&Element;S}{\&Sigma;}{A}_i\left(y_i,{\mathrm{\&phi;}}_i\left(x\right)\right)+\underset{i\&Element;S}{\&Sigma;}\underset{j\&Element;N_i}{\&Sigma;}{I}_{ij}\left(y_i,y_j,\mathrm{\&mu;}\left({\mathrm{\&phi;}}_i\left(x\right),{\mathrm{\&phi;}}_j\left(x\right)\right)\right)\right\}$

式中，φ(x)表示对观测数据x进行特征映射，A_i(.)和I_ij(.)分别是单元势函数和成对势函数，A_i(.)表示给定观测数据x_i时相应块i被标记为y_i的期望，I_ij(.)描述块i和块j之间的相互影响；φ_i(x)将块i的观测数据x_i映射到特征空间向量；μ(φ_i(x),φ_j(x))是成对块(i,j)的新特征向量；Z(φ(x))是归一化常数，表示如下，

$Z\left(\mathrm{\&phi;}\right(x\left)\right)=\underset{y}{\&Sigma;}\exp \left\{\underset{i\&Element;S}{\&Sigma;}{A}_i\left(y_i,{\mathrm{\&phi;}}_i\left(x\right)\right)+\underset{i\&Element;S}{\&Sigma;}\underset{j\&Element;N_i}{\&Sigma;}{I}_{ij}\left(y_i,y_j,\mathrm{\&mu;}\left({\mathrm{\&phi;}}_i\left(x\right),{\mathrm{\&phi;}}_j\left(x\right)\right)\right)\right\}$

步骤3，单元势函数A_i(y_i,φ_i(x))通过logistic函数来建模，表示如下，

$\left\{\begin{array}{c}{A}_i(y_i,{\mathrm{\&phi;}}_i(x\left)\right)=log(1+e^{-y_i{w}^T{h}_i\left(x\right)})\\ h_i\left(x\right)={\&lsqb;T_1\left(m_i\right),T_2\left(m_i\right),...,T_n\left(m_i\right),1\&rsqb;}^T\end{array}\right.$

式中，h_i(x)表示第i个二进方块的特征矩阵，由MLFD的响应值构成，T_n(m_i)表示该方块在尺度n下的MLFD最大响应值；w＝[w₁,w₂,...,w_K,α₁]^T是模型参数，包含特征h_i(x)的各项权重；

步骤4，根据成对势函数I_ij(y_i,y_j,μ(φ_i(x),φ_j(x)))，设定采用Beamlet分解方式寻找最优划分；

步骤5，采用相关约束来对道路的基元进行编组和标记，描述成对势函数，合并从全局数据x中得到的信息，表达如下式，

I_ij(y_i,y_j,μ(φ_i(x),φ_j(x)))＝y_iy_jv^Tμ_ij(x) (1)

式中，v＝[v₁,v₂,…,v_m]^T表示特征矩阵μ_ij(x)中每维特征的权重参数，m表示联合特征μ_ij(x)的维度，μ_ij(x)是块i和块j的联合特征；

步骤6，对特征矩阵进行统一的归一化处理。

2.根据权利要求1所述基于条件随机场的极化SAR图像道路提取方法，其特征在于：步骤5中，所述相关约束包括夹角距离、横向距离、端点距离、临近度和曲率连接。

3.根据权利要求2所述基于条件随机场的极化SAR图像道路提取方法，其特征在于：步骤6中，归一化处理采用阈值局部化方式。

4.根据权利要求1或2或3所述基于条件随机场的极化SAR图像道路提取方法，其特征在于：采用L-BFGS优化算法训练单元势特征的参数w和成对势特征的参数v。

5.根据权利要求1或2或3所述基于条件随机场的极化SAR图像道路提取方法，其特征在于：所述多尺度线状目标检测算子MLFD，将图像迭代地分割成2×2块，构建图像金字塔，在每一层上提取图像特征，合并不同层次的特征用于模型优化。

6.根据权利要求4所述基于条件随机场的极化SAR图像道路提取方法，其特征在于：所述多尺度线状目标检测算子MLFD，将图像迭代地分割成2×2块，构建图像金字塔，在每一层上提取图像特征，合并不同层次的特征用于模型优化。