专利名称：基于回归参数变换的遥感图像专题信息提取方法
基于回归参数变换的遥感图像专题信息提取方法
技术领域：
本发明属于遥感技术领域，具体是一种基于回归参数变换的遥感图像降维压缩专题信息提取方法。
背景技术：
遥感图像专题信息提取是遥感技术应用的核心问题。所谓专题信息是指某种土地 覆盖类型，如水体、农地、道路等的专门信息。目前常用的方法可分为以下二类一是目视解 译法，指判读人员依据各种经验在打印的纸质遥感图像上提取信息或在计算机屏幕上进行 信息提取。这种方法的优点是简单易行，缺点是提取精度往往取决于判读人员的经验，工作 量大，效率低。二是基于计算机技术和统计分析技术的图像分类方法，这类方法根据人工参 与程度可分为监督分类、非监督分类以及两者结合的混合分类等。非监督分类也称为聚类 分析或点群分析，即在图像中搜寻、定义其自然相似光谱集群组的过程。它不必对影像地物 获取先验知识，仅依靠影像上不同类地物光谱(或纹理)信息进行特征提取，再统计特征的 差别来达到分类的目的，最后对已分出的各个类别的实际属性进行确认。常用的非监督分 类方法有K-均值法、迭代自组织数据分析法(IS0DATA)等。监督分类又称为训练分类法， 即用已被确认类别的样本像素去识别其他未知类别像素的过程，它需要在有先验知识的条 件下进行。监督分类一般通过两步实现，即选择训练样本和提取信息。常用的分类方法包 括平行六面体法、最大似然法、最小距离法、光谱角分类法等。在实际应用中，仅仅采用无监 督分类或监督分类常常会由于所分类型与实际类别难以匹配或者样本选择不当造成分类 精度低，而两者结合的混合分类方法能够在一定程度上改善分类结果。无论监督分类还是非监督分类，都需要利用特征空间。通常特征空间的维数较高， 给存储空间造成压力，尤其是高光谱遥感数据；由于维数高，给直观分析带来麻烦，人眼感 知在二维平面空间是最佳的，在三维空间感知能力大大下降，在高于三维的空间里根本无 法感知；同时由于维数高，一些冗余数据会对信息提取结果产生不利影响。人们为解决此类 问题提出了一些基于数据压缩方法的信息提取方法，如利用正交变换法。但这些数据压缩 方法也有缺点。如正交变换法有线性假定，当数据不满足线性假定时选择的主分量中可能 包含不需要的信息，给分类带来不利影响，而未选的主分量中也可能遗留有重要的信息。
发明内容本发明需要解决的技术问题是针对多光谱图像和高光谱图像的专题信息提取问 题，提出一个全新的数据变换和压缩方法。通过数据变换，将高维数据变换到二维平面，从 而达到压缩数据、方便直观分析、消除冗余数据不良影响、提高信息提取精度的目的。为解决上述技术问题，本发明的技术方案是按如下步骤进行(1)遥感图像预处理和主要参数设定遥感图像预处理包括常规的几何精校正和 辐射校正；主要参数设定包括确定类别数m，确定目标地类，目标地类就是要提取的主要地 物；
(2)选取有广泛代表性的学习样本采用影像目视获取和地面调查验证相结合的 方法选取学习样本；(3)选择分类波段通过光谱曲线图选择对专题信息提取最有价值的特征波段， 这些特征波段上目标地类和非目标地类的平均亮度值具有较明显的差异；(4)依据目标地类进行波段排序将选择的波段依据目标地类波段的平均亮度值 从大到小或者从小到大进行排序，排序后顺序为(1，bn)，(2，bi2)，. · ·，(n，bin)，η为选中 参与信息提取的波段数，括弧中第一个数据为序号，第二个数据为波段平均亮度值，其中 bik(l ^ k ^ η) ^aE: bn ^ bi2 bin 或 bn ^ bi2 bin ；(5)线性回归计算以排序后波段序号(1，2，…，η)为自变量χ，对应的波段平均 亮度值(bn，bi2，…，bin)为因变量y，对整幅图像逐个像素进行最小二乘法线性回归计算， 即为每个像素建立如下直线方程y = ax+b得到斜率a和截距b;(6)对a和b做线性拉伸、取整并建立二维直方图对a和b做线性拉伸，取整为 0-255的整数，以a为例，找出a的最大值amax和最小值amin，线性拉伸、取整的公式为[a] = [ “一 η χ 255]
a —a .
max min方括弧[]表示对括弧内的数据进行4舍5入取整，b的线性拉伸和取整与a相同； 分别以拉伸取整后的斜率和截距作为像素亮度值，得到斜率图和截距图；建立每个地类的 二维直方图，作为进一步分类的基础，二维直方图就是一个256女256数组h (256，256)，数 组中栅格h[i，j]中的数值就是学习样本中拉伸取整后斜率为i、截距为j的像素数，称为 密度值；(7)对二维直方图进行平滑处理对所有m个单地类的二维直方图进行平滑处 理；(8)直方图密度值转换成概率值将单地类的直方图密度值转换成单地类概率 值，先求单地类直方图所有栅格密度值之和，用这个和值去除所有单个栅格密度值，得到单 地类概率直方图P (256，256)，概率直方图的概率值之和为1，具体算法如下
255 255sum = ^^ Kh j)
i=0 j=0p(i，j) = h(i，j)/sum，对所有 i，j e
；(9)制作二维类代码查找表根据单地类概率直方图制作二维类代码查找表 1 (256，256)，查找表中栅格的数值是类代码，它是所有m个单地类概率直方图中对应栅格 的概率值中的最大值对应的地类代码，算法描述如下设max[Pl(i, j),p2(i, j), ···，Pm(i，j) ] = Pk(i，j) > 0 贝Ijl(i, j) = k, i, j e
式中1彡k彡m，若pk(i，j) = 0，则栅格l(i，j)的代码值由步骤(10)确定；(10)用数学形态学方法将查找表扩张到整个空间，使查找表中任一栅格均具有一 个唯一的地类代码；
(11)利用查找表对待分类遥感图像进行分类，从步骤(6)计算得到的斜率图和截 距图中依次读取每个像素的a和b，到查找表中查找该像素的地类代码；(12)输出结果。本方法的有益效果是将高维数据通过直线回归压缩到斜率和截距的二维空间， 在序号和波段平均亮度值线性关系成立的前提下，压缩造成的信息损失很小。压缩到二维 空间带来三方面好处一是减小特征空间的存储空间压力，尤其对高光谱遥感数据效果更 明显；二是二维平面最适合于人眼感知、空间想象和直观评价，方便通过人机交互方式在二 维平面上进行操作；三是减少冗余数据对信息提取结果产生的不利影响。
图1为本发明的流程图。图2为本发明实施例中训练样本的光谱曲线图，其中图2a、图2b、图2c分别为 118/39、118/40、119/39三景TM图像的光谱曲线图。图3为本发明实施例118/40景TM图像的3个单地类训练样本的斜率、截距二维 直方图，其中图3a、图3b、图3c分别为水体、阴影、其他的二维直方图。图4为本发明实施例中3个单地类直方图叠加的训练样本的斜率、截距二维直方 图，其中图4a、图4b、图4c分别为118/39景、118/40景、119/39景的二维直方图。图5为本发明实施例中用于二维直方图平滑的模板示意图。图6为本发明实施例中3个地类的直方图，其中图6a为完整的直方图，图6b为局 部放大图。图7为本发明实施例中信息提取结果图(局部)，其中图7a为原始图像，图7b为 提取结果。
具体实施方式下面结合实施例和附图对本发明作进一步详述。实施例采用TM图像水体提取例子。从水体信息提取考虑出发，只需将地物分为水 体和非水体2类即可。水体提取的困难之处是水体和阴影容易混淆，本实施例为了更好地 说明本发明方法对区分水体和阴影的能力，将地类分为水体、阴影(含地形阴影和云层阴 影)和其他共3类。阴影本身不是一种地物，是由于地形或云层遮挡太阳电磁辐射造成的 特殊影像。数据为接收于2007年的3个整景TM图像118/39、118/40、119/39，覆盖浙江北 部地区及上海、江苏、安徽部分地区。本方法的流程参见图1。步骤1 遥感图像预处理和主要参数设定。遥感图像预处理包括常规的几何精校 正和辐射校正；主要参数设定包括确定类别数m，本实施例m为3 ；确定目标地类。目标地类 即提取专题信息涉及的主要地物类别，本实施例目标地类为水体；步骤2 选取学习样本。采用影像目视获取和地面调查验证相结合的方法选取学 习样本。学习样本要有广泛的代表性。水体要包括各种水体上中下游不同水质的江水水 体，不同水质的湖泊水体，不同水质的海水水体，等等。阴影要包括图像上存在的不同程度 的地形阴影、云层阴影。其他地类则要包括除水体、阴影之外的所有地物；
步骤3 选择分类波段。通过光谱曲线图选择对专题信息提取最有价值的特征波 段。所谓最有价值是指在所选择的波段上目标地类和非目标地类的平均亮度值具有较明显 的差异，这些差异可以被用来区分目标地类和非目标地类。本实施例中，从图2的3景TM 图像光谱曲线图可看出，在TM2、TM3、TM4、TM5波段上水体和阴影、其他地物的光谱曲线走 势有较明显的差异，因此选择这4个波段的数据作为分类原始数据；步骤4 依据目标地类进行波段排序。将选择的波段根据目标地类波段的平均亮 度值从大到小或者从小到大进行排序。排序后顺序为(1，bn)，(2，bi2)，· · ·，(n，bin)，η为 选中参与信息提取的波段数，括弧中第一个数据为序号，第二个数据为波段平均亮度值。其 中bik(l彡k彡η)满足bn彡bi2 bin或bn彡bi2 bin。本实施例中水体按TM2、 TM3、TM4、TM5波段的顺序下降，见图2。调整后的顺序为(1，TM2)，(2，TM3)，(3，TM4)，(4， TM5)；步骤5:线性回归计算。以排序后波段序号(1，2，…，η)为自变量X，对应的波段 平均亮度值(bn，bi2，…，bin)为因变量y，对整幅图像逐个像素进行最小二乘法线性回归计 算，即为每个像素建立如下直线方程y = ax+b得到斜率a和截距b。本实施例中以排序后波段序号(1，2，3，4)为自变量X，对应 的波段平均亮度值(TM2，TM3，TM4，TM5)作为因变量y ；步骤6 对a和b做线性拉伸、取整并建立二维直方图。对a和b做线性拉伸，取整 为0-255的整数。方法以a为例说明，对b完全一样。找出a的最大值amax和最小值amin， 线性拉伸、取整的公式为⑷二[χ255] α 一 α ,
max min方括弧□表示对括弧内的数据进行4舍5入取整，如[5.3] =5，[5.5] =6。分 别以拉伸取整后的斜率和截距作为像素亮度值，得到斜率图和截距图。建立每个地类的二 维直方图，作为进一步分类的基础。二维直方图就是一个256女256数组h (256，256)，数组 中栅格h[i，j]中的数值就是学习样本中斜率为i、截距为j (均为拉伸取整后的数据)的 像素数，称为密度值。118/40景的三个地类的二维直方图如图3所示。将3个单地类的直 方图叠加在一起后的直方图如图4所示。从图4可以看出斜率和截距具有较好的分离性， 意味着它们对这3个地类有较好的区分能力；步骤7 对二维直方图进行平滑处理。对所有m个单地类的二维直方图进行平滑 处理(低通滤波)，本例选用7 X 7均勻窗口，如图5所示，图中数字为权重，黑色栅格为窗口 中心；步骤8 直方图密度值转换成概率值。将单地类直方图密度值转换成单地类概率 值，转换方法是求单地类直方图所有栅格密度值之和，用这个和值去除所有单个栅格密度 值，得到单地类概率直方图P(256，256)，概率直方图的概率值之和为1。具体算法如下
<formula>formula see original document page 7</formula>
步骤9 制作二维类代码查找表。根据单地类概率直方图制作二维类代码查找表 1 (256，256)，查找表中栅格的数值是类代码，它是所有m个单地类概率直方图中对应栅格 的概率值中的最大值对应的地类代码。取最大值对应的地类代码是因为不同类的直方图会 有一定重叠，如图6所示，其中图6a为完整的直方图，图6b为局部放大图。算法描述如下 设<formula>formula see original document page 8</formula><formula>formula see original document page 8</formula>式中1彡k彡m。若pk(i，j) =0，则栅格1(1，力的代码值由步骤10确定；步骤10 用数学形态学方法将查找表扩张到整个空间，使得查找表中任一栅格均 具有一个唯一的地类代码。本发明实施例采用3 * 3模板反复对查找表进行数学形态学扩 张，直至每个栅格都有一个唯一的类代码；步骤11 利用查找表对待分类遥感图像进行分类。从步骤6计算得到的斜率图和 截距图中依次读取每个像素的a和b，到查找表中查找该像素的地类代码；步骤12 输出结果，本发明实施例如图7所示(局部)，其中图7a为原始图像，图 7b为提取结果图像。
权利要求
一种基于回归参数变换的遥感图像专题信息提取方法，其特征在于按照以下步骤进行(1)遥感图像预处理和主要参数设定遥感图像预处理包括常规的几何精校正和辐射校正；主要参数设定包括确定类别数m，确定目标地类，目标地类就是要提取的主要地物；(2)选取有广泛代表性的学习样本采用影像目视获取和地面调查验证相结合的方法选取学习样本；(3)选择分类波段通过光谱曲线图选择对专题信息提取最有价值的特征波段，这些特征波段上目标地类和非目标地类的平均亮度值具有较明显的差异；(4)依据目标地类进行波段排序将选择的波段依据目标地类波段的平均亮度值从大到小或者从小到大进行排序，排序后顺序为(1，bi1)，(2，bi2)，...，(n，bin)，n为选中参与信息提取的波段数，括弧中第一个数据为序号，第二个数据为波段平均亮度值，其中bik(1≤k≤n)满足bi1≥bi2≥…≥bin或bi1≤bi2≤…≤bin；(5)线性回归计算以排序后波段序号(1，2，…，n)为自变量x，对应的波段平均亮度值(bi1，bi2，…，bin)为因变量y，对整幅图像逐个像素进行最小二乘法线性回归计算，即为每个像素建立如下直线方程y＝ax+b得到斜率a和截距b；(6)对a和b做线性拉伸、取整并建立二维直方图对a和b做线性拉伸，取整为0-255的整数，以a为例，找出a的最大值amax和最小值amin，线性拉伸、取整的公式为 <mrow><mo>[</mo><mi>a</mi><mo>]</mo><mo>=</mo><mo>[</mo><mfrac> <mrow><mi>a</mi><mo>-</mo><msub> <mi>a</mi> <mi>min</mi></msub> </mrow> <mrow><msub> <mi>a</mi> <mi>max</mi></msub><mo>-</mo><msub> <mi>a</mi> <mi>min</mi></msub> </mrow></mfrac><mo>×</mo><mn>255</mn><mo>]</mo> </mrow>方括弧[]表示对括弧内的数据进行4舍5入取整，b的线性拉伸和取整与a相同；分别以拉伸取整后的斜率和截距作为像素亮度值，得到斜率图和截距图；建立每个地类的二维直方图，作为进一步分类的基础，二维直方图就是一个256＊256数组h(256，256)，数组中栅格h[i，j]中的数值就是学习样本中拉伸取整后斜率为i、截距为j的像素数，称为密度值；(7)对二维直方图进行平滑处理对所有m个单地类的二维直方图进行平滑处理；(8)直方图密度值转换成概率值将单地类的直方图密度值转换成单地类概率值，先求单地类直方图所有栅格密度值之和，用这个和值去除所有单个栅格密度值，得到单地类概率直方图p(256，256)，概率直方图的概率值之和为1，具体算法如下 <mrow><mi>sum</mi><mo>=</mo><munderover> <mi>Σ</mi> <mrow><mi>i</mi><mo>=</mo><mn>0</mn> </mrow> <mn>255</mn></munderover><munderover> <mi>Σ</mi> <mrow><mi>j</mi><mo>=</mo><mn>0</mn> </mrow> <mn>255</mn></munderover><mi>h</mi><mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo></mrow> </mrow>p(i，j)＝h(i，j)/sum，对所有i，j∈
；(9)制作二维类代码查找表根据单地类概率直方图制作二维类代码查找表l(256，256)，查找表中栅格的数值是类代码，它是所有m个单地类概率直方图中对应栅格的概率值中的最大值对应的地类代码，算法描述如下设max[p1(i，j)，p2(i，j)，…，pm(i，j)]＝pk(i，j)＞0则l(i，j)＝k，i，j∈
式中1≤k≤m，若pk(i，j)＝0，则栅格l(i，j)的代码值由步骤(10)确定；(10)用数学形态学方法将查找表扩张到整个空间，使查找表中任一栅格均具有一个唯一的地类代码；(11)利用查找表对待分类遥感图像进行分类，从步骤(6)计算得到的斜率图和截距图中依次读取每个像素的a和b，到查找表中查找该像素的地类代码；(12)输出结果。
全文摘要
一种基于回归参数变换的遥感图像专题信息提取方法，主要经过选择分类波段、依据目标地类亮度值排序、分别以序号和亮度值为自变量和因变量建立直线回归方程、将斜率a和截距b线性拉伸并取整为0-255的整数、建立二维直方图并进行平滑、将单类直方图密度值转换成概率值、根据单类概率直方图制作二维类代码查找表、用数学形态学方法将查找表扩张到整个空间、计算图像每个像素的a和b并拉伸取整后到查找表中查找该像素的地类代码从而实现信息提取等步骤来实现。本方法的优点是将高维数据变换到二维平面，从而达到压缩数据、方便直观分析、消除冗余数据不良影响、提高信息提取精度的目的。
文档编号G06K9/46GK101799871SQ201010039798
公开日2010年8月11日 申请日期2010年1月18日 优先权日2010年1月18日
发明者杜华强, 王鑫, 葛宏立 申请人:浙江林学院