Segmentation of PolSAR Data Based on Mean-Shift and Spectral Graph Partitioning and Its Evaluation
-
摘要: 提出了一种基于均值漂移和谱图分割的极化SAR(PolSAR)影像分割方法。首先,通过均值漂移算法对PolSAR影像进行过分割处理,并基于Wishart统计分布和假设检验的方法构建边缘检测器,充分利用了PolSAR影像的全极化信息提取边缘信息;然后,在过分割和边缘信息的基础上构建相似性度量矩阵,并采用归一化割准则实现PolSAR影像的分割。该算法充分利用了均值漂移算法过分割的特点,降低了谱图分割算法的运算代价,并结合谱图分割算法全局优化的优点改善了PolSAR影像的分割结果;最后,利用Radar-sat-2全极化影像进行了实验,并采用改进的分割效果评价方法实现了精度评价。实验表明,该算法有效地实现了PolSAR影像的分割,显著提高了谱图分割算法的效率,分割结果优良,分割精度优于eCognition软件中的多尺度分割方法。
-
关键词:
- 极化合成孔径雷达 /
- 影像分割 /
- 均值漂移 /
- 谱图分割 /
- 归一化割中图法分类号:P237.3 /
- TP79 文献标志码:A
Abstract: A new segmentation approach of polarimetric synthetic aperture radar(PolSAR)data is pro-posed based on mean shift and spectral graph partitioning.First,Mean-shift algorithm is used to gen-erate the over-segmentation of PolSAR image.In order to extract edge information,we apply a set ofdetectors based on the Wishart distribution with a hypothesis testing method that has fully consideredthe polarization information in PolSAR images.Then,a similarity matrix is constructed based on theover-segmentation results and image edge information.The graph partitioning process is performed u-sing the normalized cut criterion.With this method,we improve the segmentation efficiency of spec-tral graph partitioning based on the over-segmentation results generated by Mean-shift.The quality ofsegmentation results is also improved as a result of the global optimization of spectral graph partitio-ning algorithm.We applied this method on Radarsat-2full polarization images and evaluated the seg-mentation results.The experiement showed that this scheme can realize PolSAR segmentation effec-tively,speed up the original algorithm,and also demonstrates a better result than eCognition’s Multi-resolution segmentation approach.-
Keywords:
- PolSAR /
- segmentation /
- Mean-Shift /
- spectral graph partition /
- normalized cut
-
数字高程模型(digital elevation model, DEM)是对地球表面地形地貌的离散化数字表达，其主要借助插值算法对不规则采样数据插值完成[1-2]。受仪器噪声、人为操作不当以及外界环境不理想等因素影响，DEM源数据不可避免地含有采样误差，即误差是空间采样数据的基本属性[3-6]。插值过程中，采样误差能够被传播和放大，进而显著影响DEM及其派生品的精度。为此，众多研究人员建议采用光滑插值算法代替准确插值用于曲面建模。例如，文献[7]使用克里金(Kriging)法值移除地球物理数据中的噪声，实现光滑地学曲面建模。文献[8]认为薄板样条算法具有较强光滑性，可以有效去除采样数据中误差。鉴于不规则三角网(triangulated irregular network, TIN)线性插值的理论缺陷，文献[9]借助二元二次多项式曲面加权平均实现了DEM光滑插值。考虑到航天飞机雷达地形测绘任务(Shuttle Radar Topography Mission, SRTM)数据受噪声影响，文献[9]利用Kriging值完成SRTM重采样。为了抑制采样数据中噪声的影响，文献[10]发展了高精度曲面建模的光滑插值算法。理论上，光滑插值算法仅降低采样数据中的高程误差分量对曲面建模的影响，而忽略了水平误差分量。尽管仅移除高程误差分量可实现DEM光滑建模，但水平误差分量会导致DEM及其派生品(如地形和水系特征线)产生平面位置误差[11]，进而降低DEM在洪水风险分析[12]、山体滑坡[13]、遥感影像配准[14]等领域数字模拟结果的精度。因此，在水平位置精度要求较高的应用中，DEM建模时需同时考虑采样数据中的水平和高程误差。
总体最小二乘算法(total least squares, TLS)是近30年发展起来的一种能同时顾及系数矩阵和观测向量误差的数据处理方法[15]。在曲线(面)拟合中，TLS通过极小化所有采样点到待定曲线(面)的正交距离之和，来降低采样数据水平和高程误差分量对模拟结果的影响。大量实例分析表明，当采样点受水平和高程误差分量影响时，TLS模拟结果精度高于最小二乘算法[16-17]。实际应用中，TLS自身并不能作为插值算法直接用于DEM建模。传统插值算法中，多面函数(digital elevation model, MQ)以其高计算精度已成为DEM建模的主要技术手段[18-19]。因此，本文以MQ为基函数，借助TLS处理误差思路，发展总体最小二乘MQ算法(MQ-T)，以降低采样点水平和高程误差分量对模拟结果的影响，进而实现DEM高精度建模。
1 MQ-T原理
1.1 传统MQ
在曲面建模中，MQ基于数值逼近理论，假设任意平滑的规则或者不规则曲面均可以由简单曲面相互叠加进而逼近到任意精度等级[20]。MQ求算未知点(xi, yi)函数值的解析表达式为[21]：
$$ f\left( {{x_i},{y_i}} \right) = \sum\limits_{j = 1}^n {{w_j}\varphi \left( {{r_{ij}}} \right)} + \sum\limits_{j = 0}^m {{b_j}{p_j}\left( {{x_i},{y_i}} \right)} $$ (1) 式中，m表示多项式阶数，本文取m=1；n表示采样点数；pj(xi, yi)和bj分别表示多项式的项和对应权重系数；rij表示待求点到第j个采样点的距离；φ(r)和w分别表示映射函数和对应权重。因此，基于式(1)求算未知点函数值f(xi, yi)的前提是已知b和w。
因采样点总数小于未知参数个数(n < n+m)，直接利用最小二乘算法无法解算未知参数值。因此，可借助岭回归推导MQ基本方程组。岭回归实质上是一种改良的最小二乘估计法，通过放弃最小二乘法的无偏性，以损失部分信息、降低精度为代价获得回归系数更为符合实际、更可靠的回归方法[22]。岭回归目标函数是在最小二乘方法目标函数基础上加一正则项。基于岭回归获得MQ目标函数为：
$$ \left\{ \begin{array}{l} \mathop {\min }\limits_{\mathit{\boldsymbol{w}},\mathit{\boldsymbol{e}}} F\left( {\mathit{\boldsymbol{w}},\mathit{\boldsymbol{e}}} \right) = \frac{1}{2}{\mathit{\boldsymbol{w}}^{\rm{T}}}\mathit{\boldsymbol{w}} + \frac{c}{2}{\mathit{\boldsymbol{e}}^{\rm{T}}}\mathit{\boldsymbol{e}}\\ {\rm{s}}.\;{\rm{t}}.\;{f_i} = {\mathit{\boldsymbol{\varphi }}_i}\mathit{\boldsymbol{w}} + {\mathit{\boldsymbol{p}}_i}\mathit{\boldsymbol{b}} + {e_i},i = 1 \cdots n \end{array} \right. $$ (2) 式中，fi=f(xi, yi)；w=[w1…wn]T，b=[b0…bm]T；e=[e1…en]T，表示采样点高程误差；φi=[φ(ri1)]…φ(rin)]; pi=[p0(xi, yi)…pm(xi, yi)]；c表示光滑参数。
为了求式(2)极小值，由拉格朗日乘数法构建辅助函数为：
$$ \begin{array}{l} L\left( {\mathit{\boldsymbol{w}},\mathit{\boldsymbol{e}},\mathit{\boldsymbol{b}},\mathit{\boldsymbol{\alpha }}} \right) = \frac{1}{2}{\mathit{\boldsymbol{w}}^{\rm{T}}}\mathit{\boldsymbol{w}} + \frac{c}{2}{\mathit{\boldsymbol{e}}^{\rm{T}}}\mathit{\boldsymbol{e}} + \\ \sum\limits_{j = 1}^n {{\alpha _j}\left( {{f_i} - {\mathit{\boldsymbol{\varphi }}_i}\mathit{\boldsymbol{w}} - {\mathit{\boldsymbol{p}}_i}\mathit{\boldsymbol{b}} - {e_i}} \right)} \end{array} $$ (3) 式中，α表示拉格朗日算子。
对L(w, e, b, α)求关于w、e、b和α的导数并等于零得：
$$ \left\{ \begin{array}{l} \frac{{\partial L\left( {\mathit{\boldsymbol{w}},\mathit{\boldsymbol{e}},\mathit{\boldsymbol{b}},\mathit{\boldsymbol{\alpha }}} \right)}}{{\partial \mathit{\boldsymbol{w}}}} = 0 \Rightarrow \mathit{\boldsymbol{w}} = \sum\limits_{i = 1}^n {{\mathit{\boldsymbol{\alpha }}_i}\mathit{\boldsymbol{\varphi }}_i^{\rm{T}}} \\ \frac{{\partial L\left( {\mathit{\boldsymbol{w}},\mathit{\boldsymbol{e}},\mathit{\boldsymbol{b}},\mathit{\boldsymbol{\alpha }}} \right)}}{{\partial {\mathit{\boldsymbol{e}}_i}}} = 0 \Rightarrow c{\mathit{\boldsymbol{e}}_i} = {\mathit{\boldsymbol{\alpha }}_i},i = 1 \cdots n\\ \frac{{\partial L\left( {\mathit{\boldsymbol{w}},\mathit{\boldsymbol{e}},\mathit{\boldsymbol{b}},\mathit{\boldsymbol{\alpha }}} \right)}}{{\partial \mathit{\boldsymbol{b}}}} = 0 \Rightarrow \sum\limits_{i = 1}^n {{\mathit{\boldsymbol{p}}_i}{\mathit{\boldsymbol{\alpha }}_i}} = 0\\ \frac{{\partial L\left( {\mathit{\boldsymbol{w}},\mathit{\boldsymbol{e}},\mathit{\boldsymbol{b}},\mathit{\boldsymbol{\alpha }}} \right)}}{{\partial {\mathit{\boldsymbol{\alpha }}_i}}} = 0 \Rightarrow {\mathit{\boldsymbol{\varphi }}_i}\mathit{\boldsymbol{w}} + \\ \;\;\;\;\;\;{\mathit{\boldsymbol{p}}_i}\mathit{\boldsymbol{b}} + {e_i} = {f_i},i = 1 \ldots n \end{array} \right. $$ (4) 式(4)中消去w和e得：
$$ \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{K}} + {\rm{diag}}\left( {\frac{1}{c}} \right)}&\mathit{\boldsymbol{P}}\\ {{\mathit{\boldsymbol{P}}^{\rm{T}}}}&\mathit{\boldsymbol{0}} \end{array}} \right]\left[ \begin{array}{l} \mathit{\boldsymbol{\alpha }}\\ \mathit{\boldsymbol{b}} \end{array} \right] = \left[ \begin{array}{l} \mathit{\boldsymbol{f}}\\ \mathit{\boldsymbol{0}} \end{array} \right] $$ (5) 式中，P=(pij)，pij=pj(xi, yi)；α=[α1…αn]T；f=[f1…fn]T；K=(kij)为MQ径向基函数，由核方法推导而来，即kij=k(rij)=φiφjT。文献[23]分析认为，三次曲面作为核函数其精度优于Hardy所采用的双曲面法，故本文采用k(r)=1+r3。
将$ \mathit{\boldsymbol{w}} = \sum\limits_{i = 1}^n {{\alpha _i}\mathit{\boldsymbol{\varphi }}_i^{\rm{T}}} $代入式(1)得：
$$ \begin{array}{l} f\left( {{x_i},{y_i}} \right) = \sum\limits_{j = 1}^m {{\alpha _j}{\mathit{\boldsymbol{\varphi }}_i}\mathit{\boldsymbol{\varphi }}_j^{\rm{T}}} + \sum\limits_{j = 0}^m {{b_j}{p_j}\left( {{x_j},{y_i}} \right)} \\ \;\;\;\;\; = \sum\limits_{j = 1}^m {{\alpha _j}{\mathit{\boldsymbol{\varphi }}_i}\mathit{\boldsymbol{\varphi }}_j^{\rm{T}}} + \sum\limits_{j = 0}^m {{b_j}{p_j}\left( {{x_j},{y_i}} \right)} \end{array} $$ (6) 其矩阵形式为：
$$ \mathit{\boldsymbol{f}} = \mathit{\boldsymbol{K\alpha + Pb}} $$ (7) 因此，基于MQ插值计算时必须首先解算方程组(5)获得α和b的值，然后基于式(6)或式(7)求算未知点函数值。由以上推导可见，MQ计算中仅考虑采样点高程误差对模拟结果的影响，而忽略了水平误差分量。当采样点同时含有水平和高程误差时，此方案势必使得MQ插值结果产生平面位置误差。
1.2 MQ-T
理论上，式(1)可看作为原始特征空间中的超平面方程，其法向量为(w －1)。因此，基于整体最小二乘算法思想，点的平面位置误差为点到超平面的垂直距离。对于MQ而言，采样点平面位置误差可表达为：
$$ {\varepsilon _i} = \frac{{{f_i} - {\mathit{\boldsymbol{\varphi }}_i}\mathit{\boldsymbol{w}} - {\mathit{\boldsymbol{p}}_i}\mathit{\boldsymbol{b}}}}{{\sqrt {1 + {\mathit{\boldsymbol{w}}^{\rm{T}}}\mathit{\boldsymbol{w}}} }} $$ (8) 因此，MQ-T的目标函数表达为：
$$ \left\{ \begin{array}{l} \mathop {\min }\limits_{\mathit{\boldsymbol{w}},\mathit{\boldsymbol{e}}} F\left( {\mathit{\boldsymbol{w}},\mathit{\boldsymbol{\varepsilon }}} \right) = \frac{1}{2}{\mathit{\boldsymbol{w}}^{\rm{T}}}\mathit{\boldsymbol{w}} + \frac{c}{2}{\mathit{\boldsymbol{\varepsilon }}^{\rm{T}}}\mathit{\boldsymbol{\varepsilon }}\\ {\rm{s}}.\;{\rm{t}}.\;{f_i} = {\mathit{\boldsymbol{\varphi }}_i}\mathit{\boldsymbol{w}} + {\mathit{\boldsymbol{p}}_i}\mathit{\boldsymbol{b}} + {\varepsilon _i}\sqrt {1 + {\mathit{\boldsymbol{w}}^{\rm{T}}}\mathit{\boldsymbol{w}}} ,i = 1 \cdots n \end{array} \right. $$ (9) 类似MQ推导过程，对式(9)由拉格朗日乘数法构建辅助函数，求其关于w、ε、b、α的导数并等于零得：
$$ \left\{ \begin{array}{l} \frac{{\partial L\left( {\mathit{\boldsymbol{w}},\mathit{\boldsymbol{e}},\mathit{\boldsymbol{b}},\mathit{\boldsymbol{\alpha }}} \right)}}{{\partial \mathit{\boldsymbol{w}}}} = 0 \Rightarrow \mathit{\boldsymbol{w}} = \sum\limits_{i = 1}^n {{\mathit{\boldsymbol{\alpha }}_i}\mathit{\boldsymbol{\varphi }}_i^{\rm{T}}} \\ \frac{{\partial L\left( {\mathit{\boldsymbol{w}},\mathit{\boldsymbol{\varepsilon }},\mathit{\boldsymbol{b}},\mathit{\boldsymbol{\alpha }}} \right)}}{{\partial {\varepsilon _i}}} = 0 \Rightarrow {\mathit{\boldsymbol{\varphi }}_i}\mathit{\boldsymbol{w}} + \sqrt {1 + {\mathit{\boldsymbol{w}}^{\rm{T}}}\mathit{\boldsymbol{w}}} {\varepsilon _i}\sqrt {1 + {\mathit{\boldsymbol{w}}^{\rm{T}}}\mathit{\boldsymbol{w}}} ,i = 1 \cdots n\\ \frac{{\partial L\left( {\mathit{\boldsymbol{w}},\mathit{\boldsymbol{e}},\mathit{\boldsymbol{b}},\mathit{\boldsymbol{\alpha }}} \right)}}{{\partial \mathit{\boldsymbol{b}}}} = 0 \Rightarrow \sum\limits_{i = 1}^n {{\mathit{\boldsymbol{p}}_i}{\mathit{\boldsymbol{\alpha }}_i}} = 0\\ \frac{{\partial L\left( {\mathit{\boldsymbol{w}},\mathit{\boldsymbol{e}},\mathit{\boldsymbol{b}},\mathit{\boldsymbol{\alpha }}} \right)}}{{\partial {\mathit{\boldsymbol{\alpha }}_i}}} = 0 \Rightarrow {\mathit{\boldsymbol{\varphi }}_i}\mathit{\boldsymbol{w}} + \sqrt {1 + {\mathit{\boldsymbol{w}}^{\rm{T}}}\mathit{\boldsymbol{w}}} = {f_i},i = 1 \cdots n \end{array} \right. $$ (10) 对式(10)消去w和ε得：
$$ \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{K}} + {\rm{diag}}\left( {\frac{1}{c}} \right)}&\mathit{\boldsymbol{P}}\\ {{\mathit{\boldsymbol{P}}^{\rm{T}}}}&\mathit{\boldsymbol{0}} \end{array}} \right]\left[ \begin{array}{l} \mathit{\boldsymbol{\alpha }}\\ \mathit{\boldsymbol{b}} \end{array} \right] = \left[ \begin{array}{l} \mathit{\boldsymbol{f}}\\ \mathit{\boldsymbol{0}} \end{array} \right] $$ (11) 因wTw=αTKα，则式(11)迭代表达式为：
$$ \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{K}} + {\rm{diag}}\left( {\frac{{1 + {{\left( {{\mathit{\boldsymbol{\alpha }}^{\rm{T}}}} \right)}^{\left( i \right)}}K{\alpha ^{\left( i \right)}}}}{c}} \right)}&\mathit{\boldsymbol{P}}\\ {{\mathit{\boldsymbol{P}}^{\rm{T}}}}&\mathit{\boldsymbol{0}} \end{array}} \right]\left[ \begin{array}{l} {\mathit{\boldsymbol{\alpha }}^{\left( {i + 1} \right)}}\\ {\mathit{\boldsymbol{b}}^{\left( {i + 1} \right)}} \end{array} \right] = \left[ \begin{array}{l} \mathit{\boldsymbol{f}}\\ \mathit{\boldsymbol{0}} \end{array} \right] $$ (12) 式中，i为迭代次数。
比较式(5)和式(12)得出，MQ为MQ-T特殊情况，即当αTKα=0时的情况。因此，式(12)计算时，应先借助式(5)求算α，然后基于式(12)迭代计算，直至α和b收敛。
2 DEM建模算法精度分析
本文分别以数值实验和实例分析验证MQ-T计算精度。数值实验中，以高斯合成曲面为研究对象，分别设计受不同误差分量污染的采样数据，基于MQ-T曲面建模，并与传统MQ方法计算结果比较，验证模型的去噪性。实例分析中，以全站仪获取的采样数据为研究对象，借助MQ-T构建测区DEM，并将其与传统插值算法(如IDW，Kriging和ANUDEM)计算结果比较，验证模型的实用性。
2.1 数值实验
高斯合成曲面表达式为：
$$ \begin{array}{l} z\left( {x,y} \right) = 3{\left( {1 - x} \right)^2}{{\rm{e}}^{ - x2 - {{\left( {y + 1} \right)}^2}}} - \\ 10\left( {\frac{x}{5} - {x^3} - {y^5}} \right){{\rm{e}}^{ - {x^2} - {y^2}}} - \frac{1}{3}{{\rm{e}}^{ - {{\left( {x + 1} \right)}^2} - {y^2}}} \end{array} $$ 计算区域为[－3, 3]×[－3, 3]，模拟区域网格数为101 × 101。从高斯合成曲面中随机采集31×31个采样点用于曲面建模，采样点分别受如下误差影响。
例1：仅受高程误差分量影响。其中，误差来源于正态分布N (0, 0.02)。
例2：仅受水平误差分量影响。其中，x和y坐标误差均来源于N (0, 0.01)。因此，例1中高程误差分量方差和例2中水平误差分量方差相等。此方案可用于比较水平和高程误差对模拟结果精度影响的程度。
例3：同时受水平和高程误差分量影响。其中，水平和高程误差分别来源于正态分布N (0, 0.01)和N (0, 0.02)。
实例中分别采用中误差(root mean squared error，RMSE)和平均误差(mean error, ME)作为精度指标比较MQ-T和MQ模拟结果精度。中误差和平均误差计算的公式为：
$$ {\rm{RMSE}} = \sqrt {\frac{{\sum\limits_{i = 1}^k {{{\left( {{z_i} - {{\hat z}_i}} \right)}^2}} }}{k}} $$ (13) $$ {\rm{ME}} = \frac{{\sum\limits_{i = 1}^k {\left( {{z_i} - {{\hat z}_i}} \right)} }}{k} $$ (14) 式中，zi和${{\hat z}_i} $分别表示第i个检核点的真实值和模拟值；k表示检核点个数。
由表 1可见，当采样点仅受高程误差分量影响时(例1)，两种方法模拟结果相当，表明MQ-T对高程误差具有较好的光滑性。当采样点仅受水平误差影响时(例2)，两种方法模拟结果精度均低于例1，且MQ-T误差小于MQ误差，表明MQ-T能较好消除采样点水平误差分量对曲面模拟的影响。当采样点水平和高程坐标分量均受误差影响时(例3)，两种方法模拟结果精度均低于例1和例2，表明采样点精度越低，模拟结果精度越低；MQ-T中误差明显小于MQ，但前者平均误差略大于后者。综上所述，当采样点含有水平误差分量时，以中误差为误差指标，MQ-T模拟结果精度高于MQ。
表 1 3种误差分布条件下，MQ和MQ-T模拟结果精度比较Table 1. Accuracy Comparison Between the Classical MQ and MQ-T Under Different Error Distributions算例 误差分布 方法 RMSE ME 水平分量 高程分量 例1 - N (0, 0.02) MQMQ-T 0.062 70.062 1 0.0070.007 例2 N (0, 0.01) - MQMQ-T 0.147 80.139 3 -0.003-0.002 例3 N (0, 0.01) N (0, 0.02) MQMQ-T 0.158 70.142 7 -0.014-0.015 2.2 实例分析
测区位于山东省济南市，测区面积为3.36 km2，最小和最大高程分别为109.6 m和364.5 m。利用全站仪获取该测区5 394个采样点，平均点间距(sample interval, SI)为25 m。基于全站仪测角和测距精度指标以及数据采集时各种观测条件，采样点平面和高程精度分别为5 cm和10 cm。为了验证采样点密度对DEM精度的影响，将原始采样数据依次抽稀，每次抽稀后数据点数和对应样点间距如表 2所示。除了MQ-T和MQ外，还采用传统插值算法如IDW、Kriging以及ANUDEM进行10 m分辨率DEM建模。为了验证DEM精度，借助实时动态定位软件RTK获取了该测区500个检核点(见图 1)。各种插值算法构建DEM精度如表 3所示。
表 2 采样点与采样点平均点间距对应关系Table 2. Relationship Between Number of Sample Points and Sample Interva采样点数 点间距/m 2 100 40 2 742 35 3 721 30 5 394 25 表 3 各种插值算法模拟结果中误差/mTable 3. RMSE of Different Interpolation Methods/m点间距 IDW Kriging ANUDEM MQ MQ-T 40 3.93 2.43 2.68 2.38 2.26 35 3.41 2.37 2.55 2.24 2.11 30 3.17 2.22 2.43 2.09 2.06 25 2.69 1.84 2.28 1.73 1.73 平均 3.30 2.22 2.49 2.11 2.04 由表 3可见，随着平均点间距增大，各种插值算法精度均有所下降。具体而言，IDW、Kriging、ANUDEM和MQ降幅分别为46%、32%、15%和38%；MQ-T精度降幅为31%。结果表明，ANUDEM模拟结果受采样点间距影响最小，而IDW影响最大。不管采用何种采样点间距，MQ-T计算精度均高于其他插值方法。MQ-T平均计算精度分别是IDW的1.61倍、Kriging法的1.09倍、ANUDEM法的1.22倍以及MQ的1.03倍。
各种插值方法基于30 m平均点间距采样点构建的10 m分辨率DEM山体阴影图如图 2所示。由图 2可见，IDW结果最粗糙，这主要是因为IDW为准确插值算法，无法光滑去噪。ANUDEM法结果最光滑，去噪效果明显，但导致很多微地形特征信息丢失(如图 2(c)中圆圈标注的区域)。Kriging法、MQ法和MQ-T法结果相似，但前两种方法模拟结果有零散坑洼和凸起(如图 2(b)和2(d)中方框标注区域)。相比较而言，MQ-T法模拟结果稍好于其他方法。但相比Kriging法，MQ-T法结果左上角的山脊有峰值削平现象(如图 2(e)中圆圈标记区域)。
![]()
图 2 各种插值算法构建DEM的山体阴影图Figure 2. Hillshades of IDW, Kriging, ANUDEM, MQ and MQ-T Algorithms3 结语
为了降低采样数据中的水平和高程误差分量对DEM建模精度的影响，本文以较高插值精度的多面函数(MQ)为基函数，借助整体最小二乘算法(TLS)思路，发展了整体最小二乘MQ算法(MQ-T)。数值实验分析表明：(1)当采样点仅受高程误差分量影响时，MQ-T法计算结果精度与MQ法相当；(2)当采样点受水平误差分量影响，以中误差为误差指标，MQ-T计算结果精度高于MQ法。实例分析中，分别基于MQ法和MQ-T法构建测区DEM，并将计算结果与传统插值算法(如IDW、Kriging、ANUDEM)比较。结果分析表明：(1)各种插值算法结果精度均随采样点密度的降低而逐步降低，其中，ANUDEM精度变化最小，IDW精度变化最大；(2)不管采用何种采样点间隔，MQ-T精度始终高于其他插值算法。
-
[1] [11]。1.2 谱图分割算法 谱图分割(spectral graph partition,SGP)是建立在谱图理论基础上的一种聚类算法,基于该理论,可以将影像分割转化为图的划分问题[6]:将待分割影像映射为加权无向图G={V,E}(V:顶点,E:边)。顶点对应影像的每个像元,带权重的边对应像元间的相似性。这样,将影像分割转化为图的最优划分问题。具体过程可以分为三步:(1)选择合适的影像特征(边缘、灰度、纹理等),建立像元间的相似性度量矩阵,完成影像到图的映射;(2)选择划分准则,将图划分为一定数量的子图;(3)将图的划分结果映射到影像空间,完成分割。因此,首先要提取PolSAR影像的边缘信息,获取分割线索。1.2.1 边缘信息提取 本文进行边缘提取的边缘检测器如图2所示,共包含4个方向的边缘检测,检测窗口可根据影像实际情况指定(3×3,5×5,7×7,…)。图2 边缘提取的检测器 Fig.2 Edge Map Calculation Detectors若以检测器中心为边缘,则两边区域应存在一定的差异性。这种差异性可基于 Wishart统计分布构建假设检验衡量: H0:Σi=Σj,H1:Σi≠Σj(2)式中,Σi、Σj分别为i、j区域的期望协方差矩阵。设Θi、Θj分别为i、j区域的样本协方差矩阵数据集。假设样本间相互独立,i、j区域间是否相似可以通过似然比检验确定。检验统计量[1]如下:LH0(^ΣΘi,Θj) ^ΣinNi^ΣjnNjQ=LH1 ^Σi,^ΣjΘi,Θj)= ^Σn(Ni+Nj)( (3)式中,L(·)为不同假设条件下的似然函数;^Σ、^Σi、^Σj分别为Σ、Σi、Σj的最大似然估计量;Ni、Nj分别为i、j区域的样本数量;n为多视视数。如果Q的值较低,则原假设H0将被拒绝。基于这个检验统计量定义两个区域的差异性度量为:D(Si,Sj)=-1lnQ=n(4)(Ni+Nj)ln^Σ-Niln^Σi-Njln^Σj 式中,Si、Sj分别代表i、j区域;i、j两个区域差异性越大,D(Si,Sj)也越大,边缘强度值也越大。采用式(4)计算每个像元上4个方向的D(Si,Sj),并保留Dmax和它的检测方向。在此基础上进行边缘优化,对于任意像元p,其边缘强度值为Dmax(p),边缘方向为θ*(p)。比较垂直于边缘方向两边像元的边缘强度值,如果Dmax(p)大于等于两边的值,则保留该值,否则,设置为零。1.2.2 相似性度量矩阵构建边缘信息作为一种有效的影像分割线索,在PolSAR影像谱图分割算法中已经有所应用[1,8]。但这些算法只适用于从像元出发的情况。本文将试图找到能够代表过分割区域空间位置的某一像元,来解决这一问题,计算方法如图3所示。图3 过分割区域中心计算Fig.3 Positioning of the Center of aOver-segmentation Object 图3中,灰色多边形为某过分割区域S,x为其内部像元,L(α)为方位角为α时像元x到该区域边界的距离,因此,像元x在过分割区域内的可扩展程度可表示为:n E(x∈S)=i∏=1L(α×i),n=2π/α (5) 根据式(5),计算区域S内的所有像元的E值,max{E(x∈S)}所对应的像元即为能够代表区域S空间位置的像元。 然后即可利用影像的边缘信息来衡量区域间的相似性程度,如图4所示。图4(a)为影像的Pauli RGB显示,S1、S2、S3为3个过分割区域,以及其相应的空间代表性像元p1、p2、p3。两个区域间的相似性度量通过代表性像元间的连线上最大的边缘强度值来体现。如图4(b)所示,p1和p3之间有明显的边缘存在,其相似性就会较低,即图4(b)中顶点间边的权重W(p1,p3)较低;同理,W(p1,p2)则较高。因此,定义差异性度量:DC(x,y)=D*(z*),z* =arg maxz∈lD*(z) (6)式中,D*(·)代表边缘强度值;l是连接像元x和y的直线;z*为该连线上具有最大边缘强度值的像元。通过高斯核函数映射即可完成权重的计算:W(x,y)=exp{-DC2(x2,y)} (7)式中,σC为核函数映射的尺度参数(2σC)。由式(7)便可完成相似性度量矩阵的构建。图4 基于边缘信息的区域间相似性权重计算Fig.4 Illustration of Extracting the DissimilarityInformation from Edge Maps1.2.3 归一化割准则 谱图分割的划分准则种类繁多,常见的划分准则及其优缺点可[6]。本文选择应用较为广泛的归一化割准则完成影像分割。假设将一个图划分为两个子图A和B,原来连接A、B后来被删去的边的权重之和在图论中称为割(cut): cut(A,B)=x∈∑A,y∈BW(x,y) (8)Shi和Malik[13]将式(8)除以表现顶点集大小的度量,完成割的归一化,即为归一化割准则:Ncut(A,B)=assoccut(A(A,B,V))+assoccut(B(B,A,V))=cut(A,B) cut(B,A)+x∈∑A,v∈VW(x,v) y∈∑B,v∈VW(y,v) (9)式中,V=A∪B;assoc(A,V)代表A中的所有顶点与该图顶点间权重的和。这时,图的划分准则即最小化归一化割Ncut(A,B)。1.3 分割评价方法 鉴于PolSAR影像特殊的统计模型,分割评价中常用的区域间对比度、区域内部一致性等评价指标不再适用。本文将以在遥感影像分割研究中应用较多的正确分割的百分数作为指标,评价PolSAR影像分割结果。陈晓秋等[14]通过累计待评价分割图中的每个区域内与参考分割图重叠的最大面积计算正确分割的百分数。该方法计算出的精度与区域个数相关,分割区域数越多,精度越高。因此,只有当分割区域数相当的情况下,上述分割精度的比较才有意义。但这种比较的公平性难以保障,而且分割精度无法同时反映过分割和欠分割的程度。本文对上述方法进行改进,计算正确分割的百分数。如图5所示,虚线区域f为参考图分割 区域,与其有重叠的待评价分割区域分别是a、b、c、d这4个区域。灰色区域e是面积最大的重叠区,其面积为Se。改进方法考虑了区域e在区域d中的面积比例(Se/Sd)。若该比例很小,则e应属于区域f的欠分割部分,相应像元属于错误分割,不应参与正确分割百分数的计算。因此,本文定义一个反映欠分割程度的参数,欠分割比例(under segmentation ratio,USR):max(Sref∩Seva(1,2,3…N))USR=1-Seva(i) (10)式中,Sref为某参考分割区域面积,Seva为待评价分割区域面积,N为重叠的Seva个数,Seva(i)为与Sref重叠最大的Seva。如图5所示,USR=1-Se/Sd。依据USR,即可通过限制欠分割程度计算正确分割的百分数。例如,设置 USR=0.2时,只有重叠区域最大且欠分割比例不大于0.2的像元才会被计入正确分割的像元数。这样计算出的分割精度对于影像过分割和欠分割的程度均有所反映,可以更客观地评价分割方法的优劣。图5 分割评价原理示意图Fig.5 Illustration Chart of SegmentationEvaluation Principle2 实验与分析2.1 实验数据 本文实验采用的数据为2013年6月16日获取的Radarsat-2全极化数据(C-波段),升轨,模式为FQ18,入射角37.56°,方位向、距离向分辨率分别为4.96m和4.73m。覆盖区域为内蒙古海拉尔垦区上库力农场依根区域,截取部分农田区域作为实验区域,影像多视处理后大小为312像元×292像元,图6为影像的Pauli RGB显示。 评价分割结果的参考图如图7所示,通过专图6 实验区Pauli RGB显示Fig.6 Experimental PolSAR Data in Pauli RGB 家知识和地面真实数据获得。主要步骤包括:(1)首先获取该地区同季节的光学影像,利用ArcGIS9.3等软件勾绘出地块边界;(2)在Radarsat-2过境当天,开展同步的土地覆盖类别调查,采用差分GPS设备根据现场实际情况修正勾绘边界,定位记录地块属性;(3)利用ArcGIS 9.3等软件将修正后的边界、地块属性等数据数字化,并进行适当后处理,例如将属于同一类但没有明显边界的地块合并,属于同一类但有明显边界的不予合并;(4)与Radarsat-2影像数据匹配,获得最终的参考分割图。该区域的土地覆盖类别主要包括油菜、小麦、灌木、裸露耕地和其他共5类。其中,裸露耕地因粗糙度、垄向、含水量等因素的影响表现为三种不同的极化反映,因此细分为三类。图7 分割评价参考图Fig.7 Reference for Segmentation Evaluation2.2 分割实验 经过多视化、Lee滤波(3×3)等数据预处理步骤,首先提取影像的边缘信息,检测器窗口设置为7×7。图8为边缘粗提取结果;图9为边缘优化结果。图10是不同参数设置下的均值漂移的分割结果,分割区域的块数分别为349、173及96。可以看出,单纯尺度的调整无法使均值漂移算法实现全局较优的分割,即便在最大尺度的情况下,图10中的椭圆黑框区域范围仍处于过分割的状态,而矩形黑框内已经出现了欠分割的情况,由此可见,分割过于破碎是均值漂移难以克服的缺点。 选择中等尺度的均值漂移(N=173)结果作为谱图分割的输入。代表性像元的空间位置如图11中黄色的点所示(α=π/4),可以看到,这些点均处于区域相对中心的位置。图8 边缘粗提取结果 图9 边缘优化结果 Fig.8 Raw Result Fig.9 Optimizedof Edge Map Result of Edge Map在上述基础上,构建过分割区域间的相似性度量矩阵,并通过归一化割完成分割。这一步骤中,可以通过设置最终的分割区域个数控制影像分割的尺度。图12是不同分割区域个数设置的分割结果,其分割细节的表现不同。当区域个数设置为42时,矩形黄色框内的两块裸土没有分割开,设置为48时则被分割开。其主要原因在于,两块裸土间边缘强度较弱,因此在区域个数设置较小时,它们将被合并为一个区域。 为了进一步验证本文方法的有效性,使用eCognition图像分割软件中的多尺度分割方法,对该实验区Pauli RGB影像进行了不同尺度的分割,参数设置及分割结果如图13所示。随着分割尺度的增大,区域个数逐渐减少,过分割现象得到有效缓解(椭圆黑框区域),但是同时,欠分割现象也在逐渐加剧(方形黑框区域)。这也说明了eCognition软件的多尺度分割与均值漂移类似(图10),在单一尺度设置上很难获得全局最优的分割结果。综合比较图12及图13的分割结果,相比之下,本文方法(图12(b))的分割结果最为优良。不同大小的农田块基本都得到了较好的分割,成为独立的对象,区域一致性好,这为下一步的影像解译奠定了良好的基础。 算法的运算代价,实验影像为312行×292列,共91 104个像元,若采用传统的基于像元的谱图分割算法,则需通过C291 104次运算建立91 104×91 104大小的相似性矩阵,而本文方法只需C2173次运算建立173×173大小的矩阵,运算空间及时间需求均降低约28万倍。由此可见,均值漂移预分割明显降低了谱图分割算法的运算代价,而谱图算法全局优化的策略也使最终的分割结果优于目前比较流行的算法。2.3 定量评价与分析 为了更客观地评价分割算法,选用定量的方法评价分割结果的优劣。首先,分别选取本文和eCognition多尺度方法较好的分割结果(图12(b)和图13(b))进行精度评价;然后,在其基础上,固定尺度外的参数不变,分析分割精度随分割区域个数的变化;最后,分析了均值漂移分割区域个数对分割精度的影响。 经计算,本文方法的分割精度为83.6%,eCognition方法为75.1%(USR=0.3)。图14是不同方法正确分割像元的分布图。其中白色区域为正确分割的像元。可以看出,eCognition方法的分割精度很大程度上受到分割过于破碎的影响。 分割精度随USR的变化趋势如图15所示。分割精度随 USR的增大而增大,但本文方法的分割精度始终大于eCognition方法的分割精度。USR=0时,即表明不允许任何欠分割现象,实际分割时几乎不可能出现,因此精度为0;USR=1时,表明对于欠分割不做任何控制,只要与参考图重叠面积最大,则计入正确分割的像元数。 图16是不同方法分割精度随分割区域个数的变化趋势(USR=0.3)。可以看出,峰值两侧的分割精度均呈下降趋势,这验证了本文改进的分割精度评价指标对过分割和欠分割现象均有反映,克服了以往文献[14]中评价指标的缺陷。分割精度的峰值高于eCognition方法,且在42～65范围内时,精度明显高于eCognition方法;而在曲线两端,eCognition方法的精度要高于本文方法,但这两侧都趋向于不理想的分割结果,精度也较低。 图17是本文方法分割精度随均值漂移过分割区域个数的变化趋势(其他参数设置与图12(b)相同)。可以看出,随着过分割区域个数的增加,分割精度整体呈增长趋势,虽然有一定的波动,但在区域个数达到一定数目时趋于稳定。这说明,在利用均值漂移进行过分割时,区域个数不宜过少,在区域个数达到一定数目时(大于320),本文方法会表现出较好的稳健性。图10 不同尺度的均值漂移分割结果Fig.10 Segmentation Results of Different Scales Based on Mean-Shift图11 过分割区域中心图12 不同尺度的谱图分割结果 Fig.11 Over-segmentation Regional CenterFig.12 Segmentation Results of Different Scales Based on SGP图13 eCognition多尺度分割结果(Sh:形状C:紧致度Sc:尺度N:区域个数)Fig.13 Segmentation Results of Multi-resolution Segmentation Approach Based on eCognition图14 不同分割方法正确分割像元数分布Fig.14 Distribution of Correct Segments Pixels图15 分割精度随USR变化趋势Based on Different MethodFig.15 Segmentation Accuracies Accordingto Different SUR图16 分割精度随分割结果区域个数变化趋势Fig.16 Segmentation Accuracies According toDifferent Region Number图17 分割精度随均值漂移过分割区域个数变化趋势Fig.17 Segmentation Accuracies of SGP Accordingto Different Over-segmentation Region NumberBased on Mean-Shift3 结论与展望 在谱图分割框架下,本文提出了一种新的基于均值漂移和谱图分割的PolSAR影像分割方法。利用Radarsat-2全极化数据,验证了提出的方法,基于改进的精度评价方法和分割结果参考图,对于本文方法和eCognition的多尺度分割方法进行了定量的精度评价。相关实验表明:(1)本文提出的算法有效实现了PolSAR影像分割,组合算法能够有效降低谱图分割算法的运算代价,得到理想的分割结果;(2)本文分割算法的分割精度要优于eCognition软件的多尺度分割方法,分割结果更趋近于参考图的分割结果;(3)本文提出的改进的分割精度评价方法,可以有效地评价分割结果的优劣,对于过分割和欠分割都能有所反映,可以独立于分割结果的区域个数评价分割精度;(4)在均值漂移过分割区域个数达到一定数目时,本文提出的算法具有较好的稳健性。本文提出的算法仍存在一定的不足,例如整个分割算法中最优参数选择,算法的适用性,以及是否可引入其他特征等方面,尚需进一步的研究。参 考 文 献[1] Liu B,Hu H,Wang H,et al.Superpixel-basedClassification with an Adaptive Number of Classesfor Polarimetric SAR Images[J].IEEE Trans.Geosci.Remote Sens.,2013,51(2):907-924[2] Zhang Jie.Segmentation of Polarimetric SyntheticAperture Radar[D].Qingdao:Shandong Universityof Science and Technology,2012(张杰.极化SAR影像的分割[D].青岛:山东科技大学,2012)[3] Yu Jie,Liu Zhenyu,Yan Qin,et al.SemiautomaticObject-oriented Classification of SAR Images onMultiscale[J].Geomatics and Information Scienceof Wuhan University,2013,38(3):253-257(余洁,刘振宇,燕琴,等.多尺度下的半自动面向对象SAR影像分类[J].武汉大学学报·信息科学版,2013,38(3):253-257)[4] Deng Shaoping,Li Pingxiang,Zhang Jixian,et al.Segmentation of Multi-polarized SAR Imagery Basedon the Theory of Evidence[J].Science of Survey-ing and Mapping,2011,36(6):32-37(邓少平,李平湘,张继贤,等.一种基于证据理论的多极化SAR图像分割方法[J].测绘科学,2011,36(6):32-37)[5] Yang Xin.Study on Segmentation and ClassificationAlgorithm of Polarimetric SAR Images[D].Cheng-du:University of Electronic Science and Technologyof China,2008(杨新.极化SAR图像的分割和分类算法研究[D].成都:电子科技大学,2008)[6] You Li.Research on Image Segmentation MethodBased on Spectral Clustering[D].Changsha:Na-tional University of Defense Technology,2011(由里.基于谱聚类的图像分割方法研究[D].长沙:国防科学技术大学,2011)[7] Yan Chengxin,Sang Nong,Zhang Tianxu.Surveyon Graph Theory Based Image Segmentation Tech-nique[J].Computer Engineering and Applica-tions,2006,5:11-15(闫成新,桑农,张天序.基于图论的图像分割研究进展[J].计算机工程与应用,2006,5:11-15)[8] Ersahin K,Cumming I G,Ward R K.Segmenta-tion and Classification of Polarimetric SAR Data U-sing Spectral Graph Partitioning[J].IEEE Trans.Geosci.Remote Sens.,2010,48(1):164-174[9] Fowlkes C,Belongie S,Chung F,et al.SpectralGrouping Using the Nystrm Method[J].IEEETrans.Pattern Anal Math Intell,2004,36(2):214-225[10] Ma Xiuli,Jiao Licheng.SAR Image SegmentationBased on Watershed and Spectral Clustering[J].Journal of Infrared and Millimeter Waves,2008,27(6):452-457(马秀丽,焦李成.基于分水岭-谱聚类的SAR图像分割[J].红外与毫米波学报,2008,27(6):452-457)[11] He W,Jaeger M,Reigber A,et al.Building Ex-traction from Polarimetricsar Data Using Mean Shiftand Conditional Random Fields[C].The 7th Euro-pean Conference on Synthetic Aperture Radar(EU-SAR),Friedrichshafen,Germany,2008[12] Zou Tongyuan,Yang Wen,Dai Dengxin,et al.AnUnsupervised Classification Method of PolSAR Im-age[J].Geomaticsand Information Science of Wu-han University,2009,34(8):90-95(邹同元,杨文,代登信,等.一种新的极化SAR图像非监督分类算法研究[J].武汉大学学报·信息科学版,2009,34(8):90-95)[13] Shi J,Malik J.Normalized Cuts and Image Seg-mentation[J].IEEE Trans.Pattern Anal.Mach.Intell.,2000,22(8):888-905[14] Chen Qiuxiao,Chen Shupeng,Zhou Chenghu.Seg-mentation Approach for Remote Sensing ImagesBased on Local Homogeneity Gradient and Its Eval-uation[J].Journal of Remote Sensing,2006,10(3):357-366(陈秋晓,陈述彭,周成虎.基于局域同质性梯度的遥感图像分割方法及其评价[J].遥感学报,2006,10(3):357-366) -
期刊类型引用(10)
1. 曾广泉，马韬，张孟希，戴妍，陈凯文，丁继辉，俞双恩，王中文. 基于无人机多光谱影像的不同施氮量水稻LAI反演方法研究. 江苏农业科学. 2024(20): 41-48 . 
百度学术
2. 高钰琪，许桂玲，冯跃华，王晓珂，任红军，由晓璇，韩志丽，李家乐. 基于冠层高光谱植被指数的水稻产量预测模型研究. 中国稻米. 2023(05): 38-44 . 百度学术
3. 彭晓伟，张爱军，王楠，赵丽，杨晓楠. 高光谱技术在土壤及适种作物的研究进展. 遥感信息. 2022(01): 32-39 . 百度学术
4. 王晓珂，刘婷婷，许桂玲，冯跃华，彭金凤，李杰，罗强鑫，韩志丽，卢苇，PHONENASAY Somsana. 基于冠层高光谱遥感的杂交水稻植被指数氮素营养诊断模型. 中国稻米. 2021(03): 21-29 . 百度学术
5. 王浩淼，宋苗语，李翔，扈朝阳，鲁任翔，王翔，马会勤. 无人机高光谱遥感监测葡萄长势与缺株定位. 园艺学报. 2021(08): 1626-1634 . 百度学术
6. 刘雅婷，龚龑，段博，方圣辉，彭漪. 多时相NDVI与丰度综合分析的油菜无人机遥感长势监测. 武汉大学学报(信息科学版). 2020(02): 265-272 . 百度学术
7. 陈晓凯，李粉玲，王玉娜，史博太，侯玉昊，常庆瑞. 无人机高光谱遥感估算冬小麦叶面积指数. 农业工程学报. 2020(22): 40-49 . 百度学术
8. 落莉莉，常庆瑞，武旭梅，杨景，李粉玲，王琦. 夏玉米叶片光合色素含量高光谱估算. 干旱地区农业研究. 2019(04): 178-183 . 百度学术
9. 张良培，刘蓉，杜博. 使用量子优化算法进行高光谱遥感影像处理综述. 武汉大学学报(信息科学版). 2018(12): 1811-1818 . 百度学术
10. 李亚妮，鲁蕾，刘勇. 基于PROSAIL模型的水稻田缨帽三角-叶面积指数模型及其应用. 应用生态学报. 2017(12): 3976-3984 . 百度学术
其他类型引用(17)
下载:
计量
- 文章访问数: 1371
- HTML全文浏览量: 77
- PDF下载量: 295
- 被引次数: 27
