CN107590816B - 一种基于遥感图像的水体信息拟合方法 - Google Patents

Abstract

本发明涉及水体信息提取及拟合方法，具体涉及一种基于遥感图像的水体信息拟合方法，本发明为了解决现有的常规测量很难及时掌握水体变化及水质变化情况，且常规测量有可能无法发现一些污染源和污染源的特征的缺点，而提出一种基于遥感图像的水体信息拟合方法，包括：使用水体指数法对遥感图像进行处理，得到处理后的图像；对处理后的图像使用二维Otsu阈值分割，得到阈值分割结果；将阈值分割结果与处理后的图像取交集，确定取交集后的图像的光反射率；选取函数模型，使用光反射率及函数模型对溶氧量和高锰酸盐分别进行计算，得到溶氧量的拟合程度以及高锰酸盐的拟合程度；根据预设的标准分别选取拟合程度最优的进行拟合。本发明适用于水体信息拟合。

Description

一种基于遥感图像的水体信息拟合方法

技术领域

本发明涉及水体信息提取及拟合方法，具体涉及一种基于遥感图像的水体信息拟合方法。

背景技术

常规水质分析的方法是实地采样，即通过长时间持续监测并记录水体中污染物的种类和浓度，根据所记录的数据，分析水污染变化趋势，并分级评价水质状况。为了观察不同时间不同区域需要进行多地采样多时采样，为了对观察地区水污染变化趋势有更深的了解，需要对采样后得数据点进行相应的处理，如剔除异常点等，耗费巨大的人力物力财力；随着人口增长、城区扩建、工农业污染等，水污染程度不断恶化，水污染区域持续扩散，进行采样后的部分水体样本不能代表所有区域水体污染的分布特征及水质变化情况；地表水容易受到气候变化、土地覆盖等人为因素的影响，并且水体的高度动态以及不均匀分布等特性导致通过常规测量很难及时掌握水体变化及水质变化情况；随着污染源种类的逐渐增多，通过常规测量有可能无法发现一些污染源和污染源的特征。

发明内容

本发明的目的是为了解决现有的常规测量很难及时掌握水体变化及水质变化情况，且常规测量有可能无法发现一些污染源和污染源特征的缺点，而提出一种基于遥感图像的水体信息拟合方法。

一种基于遥感图像的水体信息拟合方法，包括：

步骤一、使用水体指数法对遥感图像进行处理，得到处理后的图像；

步骤二、对处理后的图像使用二维Otsu阈值分割，得到阈值分割结果；

步骤三、将所述阈值分割结果与所述处理后的图像取交集，确定取交集后的图像的光反射率；

步骤四、从预设的函数集中依次选取函数模型，使用光反射率以及所述函数模型对溶氧量和高锰酸盐分别进行计算，得到每一种函数模型对应的溶氧量的拟合程度以及高锰酸盐的拟合程度；

步骤五、从所述溶氧量的拟合程度以及高锰酸盐的拟合程度中根据预设的标准分别选取拟合程度最优的进行拟合。

本发明的有益效果为：1、可以及时掌握水体变化及水质变化情况；2、更易提取污染源和污染源的特征；3、本发明成本低，易于实现。

附图说明

图1为本发明的基于遥感图像的水体信息拟合方法的流程图；

图2为本发明的阈值分割向量对二维直方图的划分；

图3为改进的归一化差异水体指数法、二维Otsu阈值分割法的提取结果对比图，其中图3a为RGB合成后图像，图3b为改进的归一化差异水体指数法处理后的图像，图3c为二维Otsu阈值分割处理后的图像；

图4为具体实施方式五中高锰酸盐拟合曲线图；

图5为具体实施方式五中溶氧量反演模型的曲线图；

图6为具体实施方式七中配准之后的遥感图像，其中图6a为松花江2014年5月与2015年5月配准图像；图6b为该水域2015年5月与2015年6月配准图像；图6c为该水域2015年5月与2015年7月配准图像；图6d为该水域2015年5月与2015年10月配准图像；

图7为图6中的水域的变化区域图，其中图7a为015年5月相对2014年5月变化图，图7b为2015年6月相对于2015年5月变化图；图7c为2015年7月相对2015年5月变化图，图7d为2015年10月相对于2015年5月变化图；

图8为2016年的反演结果图，其中图8a为2016年9月的反演结果图，图8b为2016年10月的反演结果图；

图9为2015年的反演结果图，其中图9a为2015年4月的反演结果图，图9b为2015年7月的反演结果图，图9c为2015年10月的反演结果图；

图10为松花江水域高锰酸盐含量分布图；

图11为松花江水域溶氧量含量分布图。

具体实施方式

具体实施方式一：本实施方式的基于遥感图像的水体信息拟合方法，包括：

步骤一、使用水体指数法对遥感图像进行处理，得到处理后的图像。

步骤二、对处理后的图像使用二维Otsu阈值分割，得到阈值分割结果。

具体而言，从不同的遥感图像中获取水体信息，通常是与水体相关的研究的第一步。传统的水体提取方法从原理上属于基于光谱特性提取水体，即利用水体与其它地物类别明显的光谱特性(水体在近红外波段吸收率高，在可见光波段反射率低)，而本发明采用水体指数法对多光谱遥感图像进行处理，增强水体并且抑制非水体，在此基础上手动选取合适的阈值进行水体提取。

本发明在改进归一化差异水体指数的基础上应用二维Otsu阈值分割进行水体提取，此方法是根据图像的二维直方图自动选取合适阈值，与手动选取阈值相比适用性更强，精度更高。

步骤三、将所述阈值分割结果与所述处理后的图像取交集，确定取交集后的图像的光反射率。

步骤四、从预设的函数集中依次选取函数模型，使用光反射率以及所述函数模型对溶氧量和高锰酸盐分别进行计算，得到每一种函数模型对应的溶氧量的拟合程度以及高锰酸盐的拟合程度。

水质监测通常是根据光谱反射率与实测水质参数浓度之间的相关性，选取一个或几个波段作为因变量，进行线性或者非线性回归，得出回归模型。进行水质监测需要的数据通常有遥感图像的光谱反射率以及实地采样点数据。其中，遥感数据最好与实地采样数据的时间保持一致，以保证时效性；其次，采样点数据不宜过少，以保证模型的准确性以及适用性。氨氮含量、水中氧含量和高锰酸盐含量和单波段反射率或波段组合反射率有一定的相关性。而在它们的敏感波段，相关系数更大。

本发明中基于此思想，采用近4年Landsat 8遥感图像以及黑龙江省肇源站点(东经124°59′20″，北纬45°28′15″)的测量数据，分析光谱反射率与实测数据的相关性，找出敏感波段，并分析在指数函数、幂函数以及多项式函数形式下，高锰酸盐和水中溶氧量的反演模型拟合情况，并依据拟合的误差平方和(Sum of the Squared Errors)、决定性系数(R-square)、均方根误差(Root Mean Squared Error)、修正后的决定系数(Adjusted R-square)选择合适的反演模型。Adjusted R-square为负数则拟合不成功，应越接近1越好；R-square区间为[-1,1]，绝对值越接近1越好；SSE和RMSE越接近0越好。

步骤五、从所述溶氧量的拟合程度以及高锰酸盐的拟合程度中根据预设的标准分别选取拟合程度最优的进行拟合。

本发明利用遥感技术进行水质监测，具有监测区域大、速度快、便于长期监测等优点，为更加准确地掌握水体及水质变化趋势，本发明在改进的归一化差异水体指数上应用二维Otsu阈值分割进行水体提取，并在水体提取的基础上根据遥感图像光谱反射率和实地参考数据之间的相关性建立水质参数的反演模型，通过建立的模型对哈尔滨地区水体进行水质分析，相比较常规水质分析的方法而言，本发明具有成本低，易于实现等优点。

具体实施方式二：本实施方式与具体实施方式一不同的是：步骤一中，所述水体指数法中归一化差异水体指数的计算方法为：

MNDWI＝(Green-SWIR1)/(Green+SWIR1)

其中MNDWI为归一化差异水体指数，SWIR1为短波红外波段。

具体而言，在各类水体提取方法中，最为广泛应用的是水体指数法。它是通过增强水体和相邻像素之间的差异来突出感兴趣的水体，同时抑制大部分噪声。

在改进的归一化差异水体指数(NDWI)的基础上，徐涵秋发现用短波红外波段代替近红外波段，可以比NDWI更加增强水体和背景之间的对比度，并提出了改进的归一化差异水体指数(MNDWI)：

MNDWI＝(Green-SWIR1)/(Green+SWIR1) (0-1)

其它步骤及参数与具体实施方式一相同。

具体实施方式三：本实施方式与具体实施方式一或二不同的是：步骤二中，Otsu阈值分割的具体过程为：

根据公式

将f(x,y)≤s*且g(x,y)＞t*；以及f(x,y)＞s*且g(x,y)≤t*的像素忽略，并设置为0或1。

二维Otsu阈值分割法是利用梯度灰度级划分二维直方图区域，并在两个投影直方图上应用传统的Otsu阈值分割两次，目的是将感兴趣区域与背景分开。该方法对噪声具有较强的鲁棒性，且运算速度较快。

假设一幅灰度图像f(x,y)(1≤x≤M,1≤y≤N)的大小为M×N，在图像中每一个像素点f(x,y)计算n×m邻域的平均灰度值，可以得到一幅平滑的g(x,y)，两者的灰度等级都为0,1，...L。其中n×m的邻域可由下面的模板表征：

L一般取255，n和m取大于1的数，它们可以不相等。在当前相关文献中邻域的选取有多种模板，通常采用8邻域模板或者4邻域模板。本发明中选取8邻域模板。

对于图像中的每个像素，我们可以获得一对(i,j)，其中i是出现在f(x,y)中的原始灰度级，j是在g(x,y)中出现的邻域平均灰度级。令c_ij表示(i,j)的概率，其联合概率为

其中

j＝0，1，...，L-1。

给定任意阈值向量(s,t)。区域A和D分别表示对象和背景。区域B和C表示边缘和噪声。让两个类C₀和C₁分别表示对象和背景；ω_i和μ_i分别表示类C_i的概率和平均值向量。

C₀和C₁的概率可以表示为：

C₀和C₁的平均值向量可以表示为：

2D直方图的总平均向量为：

在大多数情况下，远离对角线的概率可以忽略不计。很容易得出：ω₀+ω₁≈1，并且μ_T≈ω₀μ₀+ω₁μ₁。

类间离散矩阵定义为：

离散矩阵的迹表示为：

因此，最优阈值可以表示为：

满足条件：f(x,y)≤s^*且g(x,y)＞t^*；以及f(x,y)＞s^*且g(x,y)≤t^*的像素被忽略，并设置为0或1。图2显示了阈值分割向量对二维直方图的划分。

下面通过实验来对本实施方式的水体信息提取结果进行分析：

为了比较提取结果是否准确，进行彩色合成是一种比较好的方法，便于对比。本发明采用从USGS上获取的Landsat 8 1T级产品进行实验。为了增强水体，利用Green、Red以及NIR三个波段加权处理重新生成绿色波段。即R：Red，G：(Green+Red+NIR)/3，B：Green，得到新的波段，再将三个波段合成，RGB显示，合成结果如图3a所示。在Landsat 8遥感图像中，Green、Red、NIR分别对应第三、四、五波段。

改进的归一化差异水体指数法、二维Otsu阈值分割法的提取结果分别如图3b、图3c所示。

为了更加直观的比较两种的提取精度，将每种方法的提取结果与真值图中水体信息进行比较，并根据正确检出率P(TA)、虚警率P(FA)以及漏检率P(MA)的定义：

计算出正确检出率、虚警率和漏检率，分析精度，如表1-1所示。

表1-1提取精度

其它步骤及参数与具体实施方式一或二相同。

具体实施方式四：本实施方式与具体实施方式一至三之一不同的是：步骤四中，预设的函数集包括：一阶指数函数、二阶指数函数、一阶多项式、二阶多项式、三阶多项式、一阶幂函数以及二阶幂函数。

一般来说，有两种方法可以建立数学模型：第一种是机理分析方法，另一类是测试分析方法。通过机理分析方法建立的模型通常有明确的意义。测试分析是无法直接分析研究对象内部机理的情况下，将研究对象视为一个“黑匣子”，测量通过“黑匣子”的输入输出数据，并且在此基础上运用统计分析的方法，在确定好的某一类模型中选择出一个与测量数据拟合效果最好的一个模型。

本发明中反演模型的建立即属于测试分析方法。首先给出常用的反演模型函数形式，其次，分析光谱反射率与参考数据浓度之间的相关性，选取相关性最强的一个或两个波段进行回归分析，根据拟合的均方误差等选取拟合效果最好的模型。

在统计学中，回归分析是指分析两种或多种变量之间的相互依赖的关系即相关性的一种统计分析方法。按照变量的数量，可分为一元回归和多元回归分析；按照自变量和因变量之间的关系，可分为线性回归分析和非线性回归分析。

本发明中选用的水质参数反演模型函数形式如表2-1所示。

根据光谱反射率以及水质参数浓度分析出相关性，选择相关性较强的一个波段或几个波段进行回归分析，建立反演模型。其中，光谱反射率以及参考数据浓度如表2-2、表2-3所示。

表2-2光谱反射率

表2-3参考数据浓度

其它步骤及参数与具体实施方式一至三之一相同。

具体实施方式五：本实施方式与具体实施方式一至四之一不同的是：步骤四中，拟合程度具体包括：误差平方和、决定性系数、均方根误差、修正后的决定系数。

本实施方式在具体实施方式四的基础上，进一步通过误差平方和、决定性系数、均方根误差、修正后的决定系数来判断曲线的拟合程度，选取拟合程度最佳的函数模型作为最终的选择。

一个实施例的具体过程为：

将反射率以及参考数据通过相关性统计，寻找反演高猛酸盐指数与光谱反射率之间的关联程度。相关系数如表2-4所示。

表2-4反射率与参考数据间的相关系数

可以看出，前几个波段与高锰酸盐相关性较强，分别以第二波段反射率、第三波段反射率、第二波段第三波段反射率乘积、第二波段第三波段反射率和作为自变量，对应的高锰酸盐指数含量作为因变量，用幂函数、指数函数、一阶多项式、二阶多项形式等函数形式进行拟合。由于函数形式阶数过高会导致过拟合，函数变化剧烈，不符合实际水质参数变化趋势；函数形式阶数过低会导致一些点拟合不成功，精度降低。因此，选用自变量为(ρ₂+ρ₃)，即第二、第三波段反射率的和，因变量为高锰酸盐含量，函数形式为二阶指数形式，拟合效果较好。其中，高锰酸盐的拟合曲线如图4所示。

函数形式：y＝a·e^bx+c·e^dx

参数值：a＝13.17；b＝-8.742；c＝1.321；d＝4.726

拟合程度：SSE＝0.1098；R-square＝0.9596；Adjusted R-square＝0.9191；

RMSE＝0.1914

将反射率以及参考数据通过相关性统计，寻找溶氧量与光谱反射率之间的关联程度。相关系数如表2-5所示。

表2-5反射率与参考数据间的相关系数

可以看出，溶氧量与每个波段之间的相关性都不强，与第一波段相关性最强才为0.5942。可以将多个波段光谱反射率按照某种数学关系组合，再与溶氧量进行反演运算，可靠性会大幅度增加。分别以第一波段反射率、第二波段反射率、第一第二波段反射率之和、第一第二波段反射率乘积作为自变量，对应的溶氧量含量作为因变量，用幂函数、指数函数、一阶多项式、二阶多项形式等函数形式进行拟合。实验结果表明以第一、第二波段反射率的乘积作为自变量，溶氧量含量为因变量，函数类型为三阶多项式的拟合程度最好。其中，溶氧量的拟合曲线如图5所示。

函数形式：y＝p₁·x³+p₂·x²+p₃·x+p₄

参数值：p₁＝1195000；p₂＝-91810；p₃＝2183；p₄＝-8.265

拟合程度：SSE＝5.294；R-square＝0.7923；Adjusted R-square＝0.5846；RMSE＝1.328

其它步骤及参数与具体实施方式一至四之一相同。

具体实施方式六：本实施方式与具体实施方式一至五之一不同的是：

步骤五一、建立四种评价指标，具体为：

修正后的决定系数越接近于1，则拟合程度越好；

决定性系数的绝对值越接近于1，则拟合程度越好；

误差平方和误差越接近于0，则拟合程度越好；

均方根误差越接近于0，则拟合程度越好；

步骤五二、对于溶氧量和高锰酸盐含量，均计算每一种函数模型对应的的修正后的决定系数、决定性系数的绝对值、误差平方和、均方根误差；

步骤五三、统计每一种函数模型有多少个评价指标高于其他函数模型；并选出符合条件的评价指标数量最多的函数模型。

即本实施方式采用4个指标判断拟合程度是否好，选择符合指标做多的函数进行拟合。如果多个函数的符合条件指标数量同时都是最多的，就都选取，都进行拟合，然后从拟合得出的图像中选取最符合实际地形情况的作为最终选择。

本实施方式提供的是选取拟合程度最优的函数模型的一种方法，实际过程中，为了保证拟合结果更符合实际情况，可以由人工做最后的校验核准，筛选出直观上更符合水体客观条件的模型。

其它步骤及参数与具体实施方式一至五之一相同。

具体实施方式七：本实施方式与具体实施方式一至六之一不同的是：在步骤一执行之前，还包括多时相遥感图像配准步骤，具体包括：

步骤A、对多时相遥感图像进行Harris角点特征提取；

步骤B、利用归一化互相关匹配算法剔除误匹配和不匹配向量，得到配准后的图像；

步骤C、将配准后的图像作为遥感图像进行步骤一的处理。

本实施方式是针对多时相遥感图像的水质分析，本实施方式的内容是为了是为了进行摇杆图像的配准，配准后的图像可以作为步骤一种的图像输入。

具体而言，为了直观地看出水体区域的变化情况，需要对多时相遥感图像进行变化检测，在此之前需要先进行遥感图像间的配准。首先进行Harris角点特征提取，而后利用归一化互相关匹配算法(NCC)进行粗匹配，剔除误匹配和不匹配向量，基于灰度相关系数，计算配准误差从而得到配准后的图像。

(1)Harris角点检测算法

1988年Harris提出了Harris角点检测算法，如果考虑窗口内不同点贡献权重的差异，角点响应函数可以写为：

其中，

是二维高斯窗口函数。如果考虑到偏移方向的多样性，利用一阶Taylor近似，由公式(4-1)可得到Harris算法的角点响应函数：

其中，I_x和I_y是像素点(u,v)在水平和垂直方向的一阶导数：

对像素点计算三个值：

则可以得出像素点的梯度矩阵

Harris的角点响应函数(CRF)表达式由此得到：

CRF(u,v)＝det(M)-k(trace(M))²＝(AB-C²)-k(A+B)² (0-14)

其中，det(M)＝AB-C²，trace(M)＝A+B，k为常数，一般取0.04～0.06。当目标的CRF值大于等于此阈值时，该像素点即被判定为角点。

(2)NCC匹配算法

归一化互相关(Normalized Cross Correlation method,NCC)匹配算法是一种经典的统计匹配算法，通过计算模板和匹配的互相关值来确定匹配的程度。

根据向量点乘的定义：

a·b＝|a|·|b|·cosθ (0-15)

若两个向量相似，则它们的方向相同，其夹角为θ，因此可以根据cosθ的值来判断两个向量的相似性。将其推广到二维中，则

式中，R(u,v)为位置点(u,v)的归一化相关系数；N₁×N₂为匹配模板大小；x_i+u,j+v，y_i,j分别为需要匹配的两幅中(i+u,j+v),(i,j)处的灰度值。R(u,v)的值越大，证明两幅的相似性越高。

配准之后的结果如图6所示。

在已配准的每组中截取公共区域，进行水体提取，并计算变化区域。由于进行水体提取相当于将灰度进行二值化，所以，采用如下公式进行计算变化区域：

change＝abs(imt1-imt2) (0-17)

其中，imt1和imt2分别为第一幅遥感图像和第二幅遥感图像，两者做差取绝对值得出变化区域。检测结果如图7所示。

通过实地参考数据可以发现，松花江根据溶氧量分类，通常只有Ⅰ类和Ⅱ类水体，所以，为了更细致的显示分析，对于溶氧量进行更细致的分类，加以显示。

对于Ⅳ类水体，分为3-5mg/L，对于Ⅲ类水体，分为5-6mg/L，将Ⅲ类水体和Ⅳ类水体统一用绿色显示；对于Ⅱ类水体，更细致地分为6-6.5mg/L、6.5-7mg/L、7-7.5mg/L，分别用蓝色、黄色、粉色显示；对于Ⅰ类水体，更细致地分为7.5-10mg/L、10mg/L以上，分别用橙色、青色显示。其中，氨氮含量越低，水质越差，由Ⅰ类到Ⅴ类到劣Ⅴ(溶氧量浓度小于2mg/L)类，水质越来越差。以2016年为例，反演结果如图8所示。

对反演结果进行分析：由于松花江哈尔滨段水体中高锰酸盐浓度集中在第Ⅳ类水体范围内，即6-10mg/L。反演后，将高锰酸盐含量进行分类，用不同的颜色对不同高锰酸盐含量的区域加以显示。

将Ⅰ类水体、Ⅱ类水体和Ⅲ类水体统一用绿色显示；对于Ⅳ类水体，更细致地分为6-6.5mg/L、6.5-7mg/L、7-7.5mg/L、7.5-10mg/L，分别用蓝色、黄色、橙色、粉色显示；将Ⅴ类水体和劣Ⅴ类(高锰酸盐浓度大于15mg/L)水体统一用青色表示。其中，高锰酸盐含量越低，水质越好，由Ⅰ类到Ⅴ类到劣Ⅴ类，水质越来越差。以2015年为例，反演结果如图9所示。

对实验结果进行分析：

统计各类水体，计算得出不同高锰酸盐含量的各类水体所占比率以及不同氧含量的各类水体所占比率如图10、图11所示，

其它步骤及参数与具体实施方式一至六之一相同。

从图10中分析得出，松花江哈尔滨段的流域高锰酸盐含量普遍属于Ⅳ类水体，极少部分属于Ⅲ类水体，Ⅰ类水体、Ⅱ类水体、Ⅴ类水体以及劣Ⅴ类水体均不存在。Ⅳ类水体中，高锰酸盐含量大部分都属于6.5-7mg/L范围内。每年9、10月份左右，高锰酸盐大部分会集中在6-6.5mg/L范围内，水质有所提高；到夏季的时候，高猛酸盐含量会上升至7mg/L左右，水质变差。并且河道边缘的水质比水体内部的水质更差。从高锰酸盐角度分析，松花江哈尔滨流域水体污染中，高锰酸盐的污染较为严重，水质较差。

根据图11分析可知，松花江哈尔滨流域的溶氧量集中在7.5-10mg/L范围内，属于Ⅰ类水体。极少部分属于Ⅴ类水体，即浓度在3-5mg/L范围内。从含氧量角度分析，松花江哈尔滨流域的水体不存在污染，水质较好。

由图10和11可以预测得出，松花江哈尔滨流域的水体中高锰酸盐含量将会在每年的6-10月份在7mg/L左右浮动，在其余月份将会降至6mg/L左右，但总体上来看，水体仍属于Ⅳ类水体，应当采取措施改善水质；每年的10月或者11月，松花江哈尔滨流域水体中溶氧量极大部分在7.5-10mg/L范围内，水质较好，其余月份会在7mg/L左右范围内浮动，但总体而言，水体属于Ⅰ类或者Ⅱ类水体，水质较好。

本发明还可有其它多种实施例，在不背离本发明精神及其实质的情况下，本领域技术人员当可根据本发明作出各种相应的改变和变形，但这些相应的改变和变形都应属于本发明所附的权利要求的保护范围。

Claims (2)

1.一种基于遥感图像的水体信息拟合方法，其特征在于，包括：

步骤一、使用水体指数法对遥感图像进行处理，得到处理后的图像；

所述水体指数法是指归一化差异水体指数，归一化差异水体指数的计算方法为：

MNDWI＝(Green-SWIR1)/(Green+SWIR1)

其中MNDWI为归一化差异水体指数，SWIR1为短波红外波段；

步骤二、对处理后的图像使用二维Otsu阈值分割，得到阈值分割结果；

所述Otsu阈值分割的具体过程为：

根据公式

将f(x,y)≤s^*且g(x,y)>t^*；以及f(x,y)>s^*且g(x,y)≤t^*的像素忽略，并设置为0或1；

假设一幅灰度图像f(x,y)(1≤x≤M,1≤y≤N)的大小为M×N，在图像中每一个像素点f(x,y)计算n×m邻域的平均灰度值，得到一幅平滑的g(x,y)，两者的灰度等级都为0,1，...L，其中n×m的邻域由下面的模板表征：

L取255，n和m取大于1的数，邻域模板选取8；

对于图像中的每个像素，获得一对(i,j)，其中i是出现在f(x,y)中的原始灰度级，j是在g(x,y)中出现的邻域平均灰度级；令c_ij表示(i,j)的概率，其联合概率为

其中

给定任意阈值向量(s,t)，区域A和D分别表示对象和背景，区域B和C表示边缘和噪声，让两个类C₀和C₁分别表示对象和背景；ω_i和μ_i分别表示类C_i的概率和平均值向量；

C₀和C₁的概率可以表示为：

C₀和C₁的平均值向量可以表示为：

2D直方图的总平均向量为：

远离对角线的概率忽略不计，得出：ω₀+ω₁≈1，并且μ_T≈ω₀μ₀+ω₁μ₁；

类间离散矩阵定义为：

离散矩阵的迹表示为：

因此，最优阈值表示为：

满足条件：f(x,y)≤s^*且g(x,y)>t^*；以及f(x,y)>s^*且g(x,y)≤t^*的像素被忽略，并设置为0或1；

步骤三、将所述阈值分割结果与所述处理后的图像取交集，确定取交集后的图像的光反射率；

步骤四、从预设的函数集中依次选取函数模型，使用光反射率以及所述函数模型对溶氧量和高锰酸盐分别进行计算，得到每一种函数模型对应的溶氧量的拟合程度以及高锰酸盐的拟合程度；

将反射率以及参考数据通过相关性统计，寻找反演高猛酸盐指数与光谱反射率之间的关联程度，分别以第二波段反射率、第三波段反射率、第二波段第三波段反射率乘积、第二波段第三波段反射率和作为自变量，对应的高锰酸盐指数含量作为因变量，用幂函数、指数函数、一阶多项式、二阶多项形式等函数形式进行拟合，选用自变量为(ρ₂+ρ₃)，即第二、第三波段反射率的和，因变量为高锰酸盐含量，函数形式为二阶指数形式；

函数形式：y＝a·e^bx+c·e^dx

参数值：a＝13.17；b＝-8.742；c＝1.321；d＝4.726

拟合程度：SSE＝0.1098；R-square＝0.9596；Adjusted R-square＝0.9191；

RMSE＝0.1914

将反射率以及参考数据通过相关性统计，寻找溶氧量与光谱反射率之间的关联程度，分别以第一波段反射率、第二波段反射率、第一第二波段反射率之和、第一第二波段反射率乘积作为自变量，对应的溶氧量含量作为因变量，用幂函数、指数函数、一阶多项式、二阶多项形式等函数形式进行拟合，实验结果表明以第一、第二波段反射率的乘积作为自变量，溶氧量含量为因变量，函数类型为三阶多项式的拟合程度最好：

函数形式：y＝p₁·x³+p₂·x²+p₃·x+p₄

参数值：p₁＝1195000；p₂＝-91810；p₃＝2183；p₄＝-8.265

拟合程度：SSE＝5.294；R-square＝0.7923；Adjusted R-square＝0.5846；

RMSE＝1.328；

步骤五、从所述溶氧量的拟合程度以及高锰酸盐的拟合程度中根据预设的标准分别选取拟合程度最优的函数模型进行拟合；

步骤五一、建立四种评价指标，具体为：

修正后的决定系数越接近于1，则拟合程度越好；

决定性系数的绝对值越接近于1，则拟合程度越好；

误差平方和误差越接近于0，则拟合程度越好；

均方根误差越接近于0，则拟合程度越好；

步骤五二、对于溶氧量和高锰酸盐含量，均计算每一种函数模型对应的修正后的决定系数、决定性系数的绝对值、误差平方和、均方根误差；

步骤五三、统计每一种函数模型有多少个评价指标高于其他函数模型；并选出符合条件的评价指标数量最多的函数模型。

2.根据权利要求1所述的基于遥感图像的水体信息拟合方法，其特征在于，在步骤一执行之前，还包括多时相遥感图像配准步骤，具体包括：

步骤A、对多时相遥感图像进行Harris角点特征提取；

步骤B、利用归一化互相关匹配算法剔除误匹配和不匹配向量，得到配准后的图像；

步骤C、将配准后的图像作为遥感图像进行步骤一的处理。

Images