结合先验概率估计的GF-3影像水体概率估计方法

孟令奎，毛旭东，魏祖帅，张 文

1. 武汉大学遥感信息工程学院，湖北 武汉 430079； 2. 地球空间信息技术协同创新中心，湖北 武汉 430079

用遥感影像提取水体信息是水利行业解决洪涝灾害监测、水资源管理等需求的重要方式。SAR数据以其全天候、全天时的成像特点，适用于洪涝灾害信息的提取和监测[1]，在水利遥感领域发挥着重要作用。作为我国自主研发的首颗分辨率高达1 m的C波段SAR卫星[2]，GF-3的投入运行为水利遥感带来了更为广阔的发展前景。

经典的遥感影像水体提取方法包括阈值分割[3-4]、结合多种特征的分类器分割[5]、基于水平集理论的分割[1,6]、基于能量函数的图像分割[7-9]等。当研究区域环境复杂时，水体提取过程常常出现虚警、漏警现象，导致水体提取结果精度不高。这种环境复杂性带来的精度问题，一定程度上源自像元分类的不确定性[10]。根据某些属性对像元分类时，可能存在将某些与水体像元散射特征或是光谱特征相似，但实际上不是水体的像元错分为水体的现象，使得像元分类结果存在不确定性；当数据的分辨率不高时，影像中某些像元反映的实际区域可能同时包含了多种地物，在分类上存在着不确定性。此外，在灰度值分割中阈值的选取同样会造成分类结果的不确定性。显然，二值水体分布图式的水体提取结果无法反映上述不确定性。

不同于传统的图像分割思路，水体分布概率估计[11]基于贝叶斯推断，根据SAR影像后向散射系数值计算水体的概率分布图。其最大的特点在于结果图上每一点取值范围不再是离散的{0,1}集合，而是连续的[0,1]区间，该值表示该点属于水体的概率。与常规水体提取方法得到二值水体分布图形式的结果不同，水体分布概率估计方法得到的是信息量更丰富的概率图，能够定量反映SAR影像中水体像元分类存在的不确定性。水体分布概率图充分描述了影像中水体像元的分布信息，可以用索引图的形式将水体概率值映射为不同颜色、进而直观表示水体像元的分布，也可以通过阈值分割退化为二值水体提取结果，还可以应用于水位估计[12-13]、水动力学模型校准[14-16]等场景。

水体分布先验概率是水体概率估计方法中的一个重要参数。在先验概率估计方面，现有的概率估计方法[11,17]将先验概率默认置为0.5，这会使推导出的地物后向散射系数统计模型不够准确，降低概率估计结果的精度。针对先验概率估计不充分的问题，本文提出了改进的水体概率估计方法，借助k-means聚类分析分割影像的结果估计先验概率，并利用这一先验概率估计值优化了后向散射系数分布参数估计方法。本文试验以GF-3影像为数据源，结果表明改进的水体概率估计方法能够有效提升水体概率估计精度。

1 水体概率估计理论

(1)

(2)

对于影像中的某个像元，希望根据其后向散射系数σ0判断其类别，即计算p(W|σ0)。根据条件概率公式，有

pW|σ0p(σ0)=p(W)pσ0|W

(3)

则p(W|σ0)可以写作

pW|σ0=Ppriorpσ0|W/p(σ0)

(4)

(5)

(6)

相应的，结合式(2)，边际分布p(σ0)表现为含有两个簇的混合Gauss分布形式。理论上其概率密度函数恰好具有双峰曲线的形式，这与实际分析得到的直方图形态特征一致。

2 方法设计

2.1 先验概率估计

关于水体分布先验概率估计，在文献[11]的方法中，先验概率值默认设置为0.5，称为默认先验方法。其假设的情况是研究区域中水体、背景像元数相同的情况，然而该条件在大多数情况下远不成立。注意到含有水体的影像后向散射系数边际分布p(σ0)的混合Gauss分布模型形式，考虑采用k-means算法对各像元进行聚类分析，并用聚类的结果估计水体分布先验概率。

k-means算法以平方误差为标准，通过迭代对给定样本集寻找簇内样本相似度最高的簇划分，具体算法步骤不再赘述。根据p(σ0)的混合Gauss分布模型形式，结合统计分布直方图h(σ0)呈现出的双峰形态，研究区域所有像元的σ0组成的样本集合按取值大小聚集为两簇。并且根据水体像元后向散射系数普遍低于背景像元这一特点，可以直接判定聚类得到的低σ0簇为水体像元簇，高σ0簇为背景像元簇。得到聚类结果后，计算水体像元比例作为水体分布先验概率Pprior的估计值。此外，统计两簇样本的均值和标准差用于初始化后续参数估计中的迭代步骤。

2.2 σ0分布参数估计

(1) 计算研究区域影像后向散射系数统计分布直方图。统计分布直方图h(σ0)计算需要的参数为其每个条带的宽度，称为带宽。本次试验中带宽的计算公式为

(7)

式中，n为统计样本个数；IQR为样本的四分位差。研究表明[21]，由上式给出的带宽适用于Gauss分布样本数据。

(2) 叠加理论概率密度函数p(σ0)。直方图h(σ0)的取值是像元个数，概率密度函数p(σ0)取值是概率值，二者在曲线拟合之前需经过叠加[21]操作。具体方法为先对p(σ0)乘以系数A，再用Ap(σ0)对h(σ0)作曲线拟合，其中A表示直方图的面积(带宽×高度)。

(3) 非线性拟合。使用Levenberg-Marquardt算法进行非线性最小二乘拟合[11,22]，迭代计算分布参数。迭代初始化时，用聚类步骤中得到的水体、背景样本均值和标准差作为算法初值。

2.3 方法流程

综上所述，本文方法的主要流程如图1所示。由于引入了先验概率估计的步骤，理论上本文方法有两点优势：

(1) 精度更高。先用聚类估计先验概率、再用非线性拟合估计分布参数的流程会估计得到的使p(σ0)更贴合h(σ0)，即更能贴合σ0分布的实际统计规律，理论上会获得更高的概率估计精度。

(2) 效率更高。在估计模型参数初始化值设置时，默认先验的方法采用人工设置经验初值，当初值选取不合适时可能会出现求解过程中迭代次数过多、甚至找不到局部最优解的情况。本文方法用聚类结果数据初始化模型参数，效率更高，同时也避免了因初值设置不当而可能出现的找不到最优解的情况。

图1 本文方法流程Fig.1 Workflow of the algorithm

3 试验结果

3.1 试验数据

试验使用了GF-3 L2级产品，包含FSI、FSII、QPSI等3种成像模式，分别对应3、10与8 m分辨率。验证数据由研究区域的GF-1影像人工标注得来，影像的拍摄时间与试验用GF-3影像相距在15 d以内。在开始数据分析前对影像应用了窗口大小为5×5的Gamma滤波，研究表明[23-24]这一选择可以在尽可能地保持影像空间细节的情况下提高信噪比。

3.2 结果可靠性评价

选取研究区域的GF-1影像，采用人工标注的方法，得到二值水体区域图作为真值，采用可靠性图(reliability diagram)[11,25]的方式验证结果。将概率取值的[0,1]区间均匀分割为N=10个小区间Ik，即

(8)

(9)

3.3 结合先验概率估计的水体概率估计方法试验结果

以河北省鹿泉市黄壁庄水库周边区域和湖北省黄冈市龙感湖周边区域为例开展试验。这些区域内水体集中分布于水库库区或湖区及周边河道内，其他部分则均为典型的背景地物，成分简单且分布互不混淆，适用于初步测试本文方法的效果，试验结果如图2、图3和表1。图2(a)为河北鹿泉地区经辐射定标得到的GF-3后向散射系数影像。可以看出与非水地物相比，水体因其低σ0值而呈现为明显的暗色区域。图2(b)为用作参考的GF-1影像，分辨率为8 m。图2(c)为人工标注的水体区域真值图。图2(d)为水体概率估计的结果，由颜色的深浅变化对应条件概率p(W|σ0)不同大小。结合图2(c)和图2(d)来看，概率估计有效提取出了水体区域。水库库区及其周边河道的水体条件概率取值很高(pW|σ0>0.9)，与周围地物区分明显。湖北黄冈地区试验结果与鹿泉地区情况相似，如图3所示。

表1 参数估计结果

图2 鹿泉地区试验结果Fig.2 The experimental result in Luquan area

图3 黄冈地区试验结果Fig.3 The experimental result in Huanggang area

图4 水体概率分析图Fig.4 Analysis of probabilistic water body map

3.4 与默认先验方法的结果精度对比

对比试验研究区域地处湖北咸宁地区长江沿岸。该区域包含水田、沟渠、河道等零散细小水体以及大范围的滩涂、湿地等难以与水体区分的地物，可用于检验本文方法在复杂地物情况下的效果，数据及结果见图5。同时，使用默认先验概率估计方法用作对照，对比分析的结果见图6。

图5 咸宁地区试验结果Fig.5 The experimental result in Xianning area

图6 水体概率图分析Fig.6 Analysis of probabilistic water body map

3.5 先验概率对概率估计精度的影响

前文的试验表明结合先验概率估计的水体概率估计较默认先验的估计方法在精度上的提升，并且推测这一结果差异与先验概率估计准确与否有关。对于先验概率对结果带来的影响，可以作进一步的解释。

从之前的研究区域中截取7个大小为200×200像元的子区域。这些子区域中水体面积所占区域总面积的比例各不相同，即各区域的水体分布先验概率各不相同。而在这些区域实施水体概率估计时，默认先验的方法统一将先验概率估计值置为0.5，这就造成了估计值与实际值的偏差，且各子区域偏差值各不相同。根据7个子区域的试验结果，尝试分析先验概率对结果带来的影响，试验结果见表2及图7。

表2 各子区域的试验结果

图7 可靠性度量比较分析Fig.7 Comparison of reliability

默认先验方法的Re值按1到7的顺序整体呈现先减小再增大的“U”形走势。这与默认先验概率值与各区域先验概率实际值偏差大小的变化规律一致。说明Re与先验概率偏差呈正相关，意味着准确的先验概率可以提升水体概率估计精度。

根据可靠性度量的减小比例，除去第5个区域的结果，在其他6个区域上本文方法精度均优于默认先验的方法，且精度提升比例按1到7的顺序也呈两端大中间小的“U”形走势。这一现象是合理的，根据表第2、4行，本文方法对先验概率的估计较为准确，相应的先验概率偏差值保持在较小的范围内，而默认先验方法的先验概率偏差值则会出现“U”形波动。这一结果也再次验证了水体概率估计精度与先验概率偏差呈正相关关系，结合先验概率估计的水体概率估计方法精度更高。

3.6 从概率图到二值图

即97.51%的像元被正确分类。剩余2.49%分类错误像元在各概率区间的分布情况如表3所示。结合图9，分类错误像元在[0,0.1]与[0.9,1]区间占比较高，但这两个区间内像元所占像元总数比例也高，因此分类错误像元占区间像元总数的比例反而较低。与之相反的是，8个中间概率区间只含有影像10.3%的像元，分类错误部分却占了错误分类像元总数的57.5%，占各区间像元数目的比例也远远高于[0,0.1]和[0.9,1]两个区间。这说明具有中间概率值的像元更容易被错误分类。特别取出这部分像元(水体概率取值于0.1到0.9)，如图8(b)，发现普遍位于水体边界地带。该类区域多为湿地、沼泽，常覆盖有浅层水面，呈现出与水体类似的低后向散射特征，易于出现虚警、漏警现象，仅依据后向散射系数无法准确分类此种地物。这是阈值划分方法不可避免的缺陷，取0.5作为水体概率的分割阈值恰好是一种折中的思路。

图8 南昌地区试验结果Fig.8 The experimental result in Nanchang area

对水体概率图做阈值划分本质上是对后向散射系数影像做阈值划分。本节所述的由水体概率图退化得到水体二值图的方法属于阈值法提取水体的范畴。选择水体概率取0.5时的后向散射系数值作为分割阈值，对应理论分布概率密度曲线上水体与背景两类像元两组密度函数曲线的交点处，也即直方图两峰值之间的“波谷”处，与一些阈值划分方法中通过寻找直方图“波谷”确定阈值的思路一致。不过与之不同的是，本节所述水体提取方法根据两类像元理论分布推导而来，在像元分布估计准确的前提下，阈值选取的效果必然是好的。与直接在直方图上寻找“波谷”的方法相比，该方法具有更强的理论基础，同时会减少直方图中噪声带来的误差。

表3 各概率区间像元分布

注：n为区间像元数；N为像元总数；nerror为区间分类错误像元数；Nerror为分类错误像元总数。

图9 各概率区间像元分布Fig.9 The distribution of pixels in each interval

4 结 语

本文针对水体概率估计中水体分布先验概率估计不充分的问题，使用k-means聚类估计先验概率并优化了后续的参数估计步骤，改进了水体概率估计方法。在水体概率估计方法中，先验概率估计的偏差会带来模型对地物后向散射系数分布规律描述的偏差，进而影响水体概率估计精度。本文方法由于对先验概率估计更充分，进而结果精度更高。此外，本文方法使用聚类结果数据初始化模型参数让参数迭代求解过程能够更快地收敛至最优值，同时也简化了人工迭代初值设置的步骤。此外，由于该方法使用的是单极化数据，因此同样适用于多时相或是多种极化方式的影像，具有很好的推广性。

需要指出的是，本文方法的推导建立在研究区域后向散射系数统计分布直方图呈双峰形这一前提条件上，根据水域与背景间后向散射系数的差异提出模型假设。但是这一条件未必总是成立，例如当影像中水体与背景像元中的一类数目远远大于另一类时，其后向散射系数分布会被一类像元主导；或是风浪导致的水体后向散射系数变高并与背景像元混淆，直方图也无法呈现双峰形走势，影响概率估计效果。对于此类问题及本文方法适用条件的定量描述，是今后研究工作的重点。在实际应用中，可以通过设置河流、湖泊等水体目标缓冲区等方式避免研究区域水体比例出现极端值的情况，保证概率估计的效果。