本发明涉及一种基于多光谱遥感影像的海图养殖区识别方法，用于多光谱遥感影像分析与目标识别，属于图像处理和计算机视觉领域。
背景技术：
随着我国水产养殖业的发展，沿海地区养殖区范围不断扩大，在产生了巨大经济效益的同时，也对生态环境、海上交通造成了不少负面影响。每年各地都有多起运输商船误入养殖区的事件。商船进入养殖区，一方面会对水产养殖产生影响或局部破坏，产生经济纠纷；另一方面也可能发生渔网绞缠螺旋桨事故，影响船舶安全航行。目前，海图上并未对近海养殖区信息进行标注，船舶驾驶员仅能根据以往的经验来避开养殖区。针对这种航行安全现状，监视近海养殖区的分布与位置范围变化，对保障船舶安全航行具有重要的现实意义。
传统的海水养殖监测方式，包括，无人机巡视，严格地被天气制约，而且只有白天才能有效的监视部分海水养殖区；再有利用船舶进行海面巡视，一方面误闯率大，另一方面航行视线受海面气象环境的影响，降低监测识别的效率和准确性；还有进行航海雷达探测，其范围受限，只能进行小面积的养殖区监测。针对这一问题，本发明提出了一种基于多光谱遥感影像的海图养殖区识别方法及工作流程，能够客观、准确、及时地对近海养殖区域进行动态监测，在海图上及时显示养殖区域的变化更新状态，同时可实现海洋资源的合理开发利用，减少对环境破坏以及海上交通事故的发生。
目前已有的基于卫星遥感的养殖区检测方法均是在限定水域范围内针对特定类型的养殖区纹理和光谱特性提出的。养殖区的提取分为两类。一类是采用光谱特征，针对养殖用海信息的提取。比如，马艳娟，赵冬玲等在基于aster数据的近海养殖区的提取方法对比研究中，利用养殖区在影像上不同波段上的不同光谱反射特征，提出了一种区分水体和非水体的指数表达式[1]；徐京萍等以spot5卫星遥感数据为数据源，通过多尺度图像分割，结合对象的光谱、形状和语义特征，建立分类规则集，实现池塘养殖用海信息提取[2]。另一类是采用纹理特征，针对养殖区纹理的信息提取。比如，朱长明等利用光谱求倒以增强遥感图像纹理信息，借助岸线掩膜，通过纹理分析和阈值分割实现养殖区的识别[3]；
这两类方法可基本识别养殖区信息，但针对近海区域不同类型养殖区的一整套自动提取方法还少有提及。
技术实现要素：
本发明的目的在于针对目前无法动态获取大范围近海养殖区的分布与位置范围信息，并且目前海图上并未对近海养殖区信息进行标注的问题，提出一种基于多光谱遥感影像的海图养殖区识别方法。通过遥感卫星影像对近海养殖区进行提取，并形成海图养殖区数据产品。该方法可以提高对大范围近海养殖区的变化监测的工作效率，并可以减少由于缺少此类的航行障碍物信息标识所导致的航行事故。
为实现这样的目的，本发明提供一种基于多光谱遥感影像的海图养殖区识别方法，其特征在于，该识别方法包括以下步骤：
步骤一：对海岸带浅海水域的低分辨率多光谱遥感影像与高分辨率全色影像进行图像融合。为了保证融合效果，首先对全色影像和多光谱影像进行大气校正、空间配准等预处理，将多光谱波段的大小达到与全色影像一致。
(一)全色影像是传感器获取的单波段的灰度影像，多光谱影像是传感器对地物辐射中包含多个波段光谱信息的影像，图像融合时选用全色影像的单波段以及多光谱影像的蓝、绿、红光三个波段。
(二)采用施密特正交化变换的图像融合方法，流程如下：
(1)利用原始低分辨率多光谱影像的波段进行光谱重采样，模拟出一幅全色影像作为施密特正交化变换的第一个分量，对模拟的全色影像与低分辨率的多光谱影像进行施密特正交化变换，具体变换公式如下所示：
其中，gst表示经施密特正交化变换后的第t个正交分量，bt表示原始低分辨率多光谱影像的第t个波段，μt为原始低分辨率多光谱影像的第t波段像元灰度值的均值；
根据公式(2)、(3)计算得到的gs1的均值和标准差对原始高分辨率全色影像的均值和标准差进行调整，得到优化后高分辨率全色影像。
(2)用优化后高分辨率全色影像的单波段替换施密特正交化变换的第一个分量，即gs1分量。
(3)最后通过公式(4)对替换后的分量进行施密特逆变换，完成低分辨率多光谱影像与高分辨率全色影像的图像融合，得到高分辨率多光谱影像。
其中，公式(4)所有符号含义同公式(1)、(2)、(3)。
步骤二：采用局部空间相关统计指数对步骤一所得的融合后高分辨多光谱影像的养殖区纹理进行描述。
(一)根据遥感影像，对近海海域两种不同类型养殖区影像纹理特征定义如下：网箱养殖区呈规则的长方形网格状分布，在蓝、绿、红光三波段合成影像中显示为棕色和暗棕色；浮筏养殖区呈规则的带状分布，密度较大，在蓝、绿、红光三波段合成影像中显示为深黑色，边缘有白色浮球。
(二)在一阶四邻域规则下，对步骤一所得的融合后影像的绿光波段采用局部空间统计热点分析指数进行纹理特征的分析和计算。一阶四邻域(rook’scase)是指某个(x,y)点的上、下、左、右四个紧邻点，该邻域规则通过上下左右四个方向的像元进行比较；采用局部空间统计热点分析(getis-ordgi)指数对地物类别进行自相关分析，减少同物异谱和异物同谱对分类精度带来的影响。该指数可以计算每一个空间单元与邻近单元的集聚程度。集聚程度代表单元间的性质相似与否，高值集聚代表单元间性质相似，低值集聚代表单元间性质相异；局部空间统计热点分析指数相关公式为：
其中，若gi(d)值为负，表明位置i周围的值相对较小，即低于均值，表现为低的观测值的集聚；若gi(d)值为正，表明位置i周围的值相对较大，即高于均值，表现为高的观测值的集聚，养殖区属于低值集聚。公式(5)中，d为影像中空间单元的间隔距离，wij为空间权重矩阵的元素，xj为影像中空间单元j的属性值。空间权重矩阵是对空间关系的建模和理解，将影像中n个空间单元的邻接关系用1个二元的空间权重矩阵w表示，如公式(6)所示：
式中，wij为影像中空间单元i，j的邻接关系，邻接则取值为1，不邻接则取值为0。
步骤三：建立支持向量机分类模型，对近海养殖区进行分类、识别。支持向量机(supportvectormachines，svm)是一种定义在特征空间上的使两类样本之间距离最大化的二类线性分类器。算法通过核函数(kernel)技术将原始输入空间映射到更高维度的特征空间，使原本线性不可分的数据样本在新的核空间内变为线性可分。然后基于这种变换，根据有限的样本信息找到预测可能的两种分类之间的最佳边界，即最优决策超平面，使得该平面两侧距离该平面最近的两类样本之间的距离最大化。
设给定的训练集为(x1,y1)，(x2,y2)，…，(xn,yn)，其中xi∈rn为输入向量，输出向量为yi∈{-1,1}，如果该训练集可被一个超平面线性划分，则该超平面为
wx+b＝0(7)
其中，w是特征空间中分类超平面的系数向量；b是分类面的阈值；w和b决定了超平面的位置，wx为两个向量的内积。为了得到最优化的划分，则该问题就转化为求最优化超平面的问题：
s.t.yi[(w·x1)+b]≥1-ξi,ξi≥0(9)
其中，ξi是考虑分类误差而引入的松弛因子；c是对于错分样本的惩罚因子，构造出来的最优化分类器为：
f(x)＝wx+b(10)
由此，建立养殖区支持向量机分类模型的三个主要步骤为：
(一)选取训练样本
训练样本的质量和数量对分类效果和分类精度有较大影响，选取七种地物类别包括：裸地、盐湖、林地、房屋、浮筏养殖区、网箱养殖区以及水体；选取过程中反复对训练样本进行筛选，剔除相似度高的样本，确保养殖区与其他地物类别的样本可分离度达到峰值。其中，裸地呈不规则的条纹状分布，在蓝、绿、红光三波段合成影像中显示为橙棕色；盐湖呈规则的网格状分布，在蓝、绿、红光三波段合成影像中显示为浅蓝或深蓝色；林地呈不规则形状分布，边界较为模糊，在蓝、绿、红光三波段合成影像中显示为深绿或浅绿色；房屋呈不规则的多边形分布，边界较为清晰，在蓝、绿、红光三波段合成影像中显示为深灰或深紫色；水体形状自然弯曲或局部平直，色彩均匀，在蓝、绿、红光三波段合成影像中显示为孔雀蓝色；浮筏、网箱养殖区特征定义同步骤二。
(二)确定分类模型的核函数及相关参数
选择高斯核函数作为支持向量机分类模型的核函数，相关公式为：
k(xi,xj)＝exp(-γ||xi-xj||)2(11)
其中，γ>0为核参数，可以通过调节γ来匹配维度的大小，γ越大，维度越低。
使用高斯核函数需要确定公式(11)中的核参数γ和公式(8)中的惩罚因子c的取值，通过交叉对比验证的方法确定γ、c的取值，即将选取的训练样本分为n个子集，每次将其中n-1个子集作为分类模型的训练样本，剩下的一个子集作为确定模型参数γ、c的检验样本，利用检验样本验证n-1个子集的分类精度，根据公式(12)不断改变γ、c的值以获取最高的样本分类精度，最终确定γ、c的取值；
c＝2-5,2-3,…,215,γ＝2-5,2-3,…,215(12)
其中，惩罚因子c和核参数γ是独立增长的；
(三)将步骤二所得的运算结果作为纹理波段和融合后高分辨率多光谱波段一起输入到支持向量机分类模型中进行分类，实现近海浮筏和网箱养殖区的分类与识别。
设融合后高分辨率多光谱影像中蓝光波段为band1、绿光波段为band2、红光波段为band3、步骤二中所得的运算结果作为纹理波段为band4。根据公式(13)计算裸地、盐湖、林地、房屋、浮筏养殖区、网箱养殖区以及水体等七种样本类别的最优化超平面分别为：
其中，i＝1,2,3,4,5,6,7，裸地i取值为1，盐湖i取值为2，林地i取值为3，房屋i取值为4，浮筏养殖区i取值为5，网箱养殖区i取值为6，水体i取值为7；xi1、xi2、xi3分别代表第i个地物在band1、band2、band3的光谱强度特征，xi4代表第i个地物在band4的纹理特征。
步骤四：将步骤三所得的分类后的遥感影像与海图进行配准，消除两幅图像之间的几何形变，确保不存在重叠区。图像配准就是将不同时间、不同传感器或不同条件下获取的两幅或多幅图像进行匹配、叠加的过程。根据配准后的分类结果图，将养殖区提取结果转换为矢量文件，对照海图进行图斑镶嵌，实现海图上养殖区的标识。
附图说明：
图1为本发明的流程图。
图2为基于施密特正交化变换的融合方法流程图。
具体实施方式：
以下结合附图，具体说明实施方法，对本发明做进一步阐述。
步骤一：对海岸带浅海水域的低分辨率多光谱遥感影像与高分辨率全色影像进行图像融合。
为了保证融合效果，首先对海岸带浅海水域的多光谱遥感影像与高分辨率全色影像进行大气校正、空间配准等预处理工作，将多光谱波段的大小达到与全色影像一致；然后，采用施密特正交化变换方法将全色影像与多光谱影像进行融合，增强影像的空间分辨率，提升影像精确度，增加信息量，突出养殖区内部纹理结构。
(一)全色影像是传感器获取的单波段的灰度影像，多光谱影像是传感器对地物辐射中包含多个波段光谱信息的影像，图像融合时选用全色影像的单波段以及多光谱影像的蓝、绿、红光三个波段。
(二)采用施密特正交化变换的图像融合方法，流程图如图2所示。
(1)利用原始低分辨率多光谱影像的波段进行光谱重采样，模拟出一幅全色影像作为施密特正交化变换的第一个分量，对模拟的全色影像与低分辨率的多光谱影像进行施密特正交化变换，具体变换公式如下所示：
其中，gst表示经施密特正交化变换后的第t个正交分量，bt表示原始低分辨率多光谱影像的第t个波段，μt为原始低分辨率多光谱影像的第t波段像元灰度值的均值；
根据公式(2)、(3)计算得到的gs1的均值和标准差对原始高分辨率全色影像的均值和标准差进行调整，得到优化后高分辨率全色影像。
(2)用优化后高分辨率全色影像的单波段替换施密特正交化变换的第一个分量，即gs1分量。
(3)最后通过公式(4)对替换后的分量进行施密特逆变换，完成低分辨率多光谱影像与高分辨率全色影像的图像融合，得到高分辨率多光谱影像。
其中，公式(4)所有符号含义同公式(1)、(2)、(3)。
步骤二：采用局部空间相关统计指数对步骤一所得的融合后高分辨多光谱影像的养殖区纹理进行描述。
近海养殖区是在浅海海域通过浮筏、网箱等进行养殖生产的区域，在遥感影像上的养殖信息与浅海岩石、盐湖、水体等地物斑块相混淆。
(一)根据遥感影像，对近海海域两种不同类型养殖区影像纹理特征定义如下：网箱养殖区呈规则的长方形网格状分布，在蓝、绿、红光三波段合成影像中显示为棕色和暗棕色；浮筏养殖区呈规则的带状分布，密度较大，在蓝、绿、红光三波段合成影像中显示为深黑色，边缘有白色浮球。
(二)在一阶四邻域规则下，对步骤一所得的融合后影像的绿光波段采用局部空间统计热点分析指数进行纹理特征的分析和计算，局部空间统计热点分析(getis-ordgi)指数相关公式为：
其中，若gi(d)值为负，表明位置i周围的值相对较小，即低于均值，表现为低的观测值的集聚；若gi(d)值为正，表明位置i周围的值相对较大，即高于均值，表现为高的观测值的集聚，养殖区属于低值集聚。在公式(5)中，d为影像中空间单元的间隔距离，wij为空间权重矩阵的元素，xj为影像中空间单元j的属性值。空间权重矩阵是对空间关系的建模和理解，将影像中n个空间单元的邻接关系用1个二元的空间权重矩阵w表示，如公式(6)所示：
式中，wij为影像中空间单元i，j的邻接关系，邻接则取值为1，不邻接则取值为0。
步骤三：建立支持向量机分类模型，对近海养殖区进行分类、识别。
设给定的训练集为(x1,y1)，(x2,y2)，…，(xn,yn)，其中xi∈rn为输入向量，输出向量为yi∈{-1,1}，如果该训练集可被一个超平面线性划分，则该超平面为
wx+b＝0(7)
其中，w是特征空间中分类超平面的系数向量；b是分类面的阈值；w和b决定了超平面的位置，wx为两个向量的内积。为了得到最优化的划分，则该问题就转化为求最优化超平面的问题：
s.t.yi[(w·x1)+b]≥1-ξi,ξi≥0(9)
其中，ξi是考虑分类误差而引入的松弛因子；c是对于错分样本的惩罚因子，构造出来的最优化分类器为：
f(x)＝wx+b(10)
由此，建立养殖区支持向量机分类模型的三个主要步骤为：
(一)选取训练样本
训练样本的质量和数量对分类效果和分类精度有较大影响，选取七种地物类别包括：裸地、盐湖、林地、房屋、浮筏养殖区、网箱养殖区以及水体；选取过程中反复对训练样本进行筛选，剔除相似度高的样本，确保养殖区与其他地物类别的样本可分离度达到峰值。
其中，裸地呈不规则的条纹状分布，在蓝、绿、红光三波段合成影像中显示为橙棕色；盐湖呈规则的网格状分布，在蓝、绿、红光三波段合成影像中显示为浅蓝或深蓝色；林地呈不规则形状分布，边界较为模糊，在蓝、绿、红光三波段合成影像中显示为深绿或浅绿色；房屋呈不规则的多边形分布，边界较为清晰，在蓝、绿、红光三波段合成影像中显示为深灰或深紫色；水体形状自然弯曲或局部平直，色彩均匀，在蓝、绿、红光三波段合成影像中显示为孔雀蓝色；浮筏、网箱养殖区特征定义同步骤二。
选取训练样本时需要遵循以下规则：
(1)对于光谱特征较为单一的区域，根据目视解译预判只有一类地物，则选择均匀分布的多个地块；
(2)对于光谱特征变化较大的地物，为了使训练样本满足建立判别函数的要求，每类地物选取的训练样本的光谱特征分布均匀，得到较为可靠的统计数据；
(3)在选择样本的时候，选择纯净的像元，提高样本可分离性。
(4)对于检验样本的选择，其位置要与训练样本的位置不同，类别要与训练样本的类别一致。
(二)确定分类模型的核函数及相关参数
选择高斯核函数作为支持向量机分类模型的核函数，相关公式为：
k(xi,xj)＝exp(-γ||xi-xj||)2(11)
其中，γ>0为核参数，可以通过调节γ来匹配维度的大小，γ越大，维度越低。
使用高斯核函数需要确定公式(11)中的核参数γ和公式(8)中的惩罚因子c的取值，通过交叉对比验证的方法确定γ、c的取值，即将选取的训练样本分为n个子集，每次将其中n-1个子集作为分类模型的训练样本，剩下的一个子集作为确定模型参数γ、c的检验样本，利用检验样本验证n-1个子集的分类精度，根据公式(12)不断改变γ、c的值以获取最高的样本分类精度，最终确定γ、c的取值；
c＝2-5,2-3,…,215,γ＝2-5,2-3,…,215(12)
其中，惩罚因子c和核参数γ是独立增长的。
(三)将步骤二所得的运算结果作为纹理波段和融合后高分辨率多光谱波段一起输入到支持向量机分类模型中进行分类，实现近海浮筏和网箱养殖区的分类与识别。
设融合后高分辨率多光谱影像中蓝光波段为band1、绿光波段为band2、红光波段为band3、步骤二中所得的运算结果作为纹理波段为band4。根据公式(13)计算裸地、盐湖、林地、房屋、浮筏养殖区、网箱养殖区以及水体等七种样本类别的最优化超平面分别为：
其中，i＝1,2,3,4,5,6,7，裸地i取值为1，盐湖i取值为2，林地i取值为3，房屋i取值为4，浮筏养殖区i取值为5，网箱养殖区i取值为6，水体i取值为7；xi1、xi2、xi3分别代表第i个地物在band1、band2、band3的光谱强度特征，xi4代表第i个地物在band4的纹理特征。
步骤四：将步骤三所得的分类后的遥感影像与海图进行配准，消除两幅图像之间的几何形变，确保不存在重叠区。图像配准就是将不同时间、不同传感器或不同条件下获取的两幅或多幅图像进行匹配、叠加的过程。根据配准后的分类结果图，将养殖区提取结果转换为矢量文件，对照海图进行图斑镶嵌，实现海图上养殖区的标识。
综上所述，本发明公开了一种基于多光谱遥感影像的海图养殖区识别方法，抑制遥感影像中的其他地物，能够有效识别并区分近海地区两种不同类型养殖区，为进一步开展整个航行监测和沿海水产养殖区分布调查工作提供方法和技术支撑。
尽管本发明的内容已经通过上述实施例作了详细介绍，但应当认识到上述的描述不应被认为是对本发明的限制。在本领域技术人员阅读了上述内容后，对于本发明的多种修改和替代都将是显而易见的。因此，本发明的保护范围应由所附的权利要求来限定。
参考文献：
[1]马艳娟，赵冬玲。基于aster数据的近海水产养殖区提取方法对比研究[j]。测绘通报，2011，20(1)：59-63。
[2]徐京萍，赵建华，张丰收等.面向对象的池塘养殖用海信息提取[j]。国土资源遥感，2013，25(1)：82-85。
[3]朱长明，骆剑承，沈占锋等.高分辨率遥感影像近海养殖区自动识别方法[j]。大连海事大学学报，2011，37(3)：66-69。