基于多源卫星数据计算内陆水域最大叶绿素指数的方法

摘要

本发明涉及基于多源卫星数据计算内陆水域最大叶绿素指数的方法，其包括步骤：1）、数据收集；2）、获得模拟的卫星遥感反射率；3）、对共有波段的低分辨率卫星的低空间分辨率波段进行重采样；4）、基于人工神经网络算法，建立映射关系，并基于映射关系获得高精度的共有波段的低分辨率卫星的高空间分辨率波段；5）、利用所述高空间分辨率波段锐化处理非共有波段的低分辨率卫星的低空间分辨率波段，获得非共有波段的低分辨率卫星的高空间分辨率波段；6）、利用所述高空间分辨率波段计算最大叶绿素指数。其结合不同卫星的波段数据，提高其中一个卫星的波段分辨率来计算最大叶绿素指数，效果更好。

权利要求

1.基于多源卫星数据计算内陆水域最大叶绿素指数的方法，其特征在于，包括以下步骤：

1）、数据收集：收集现场实测的光谱反射率数据、遥感图像角度数据、高分辨率卫星和低分辨率卫星的波段；

2）、基于所述现场实测的光谱反射率数据和光谱响应函数获得模拟的卫星遥感反射率数据；

3）、找出高分辨率卫星和低分辨率卫星的共有波段，并对共有波段的低分辨率卫星的低空间分辨率波段进行重采样以获得共有波段的低分辨率卫星的高空间分辨率波段；

4）、基于人工神经网络算法，利用所述模拟的卫星遥感反射率数据、遥感图像角度数据、共有波段的高分辨卫星的高空间分辨率波段和重采样获得的共有波段的低分辨率卫星的高空间分辨率波段建立映射关系，并基于所述映射关系获得高精度的共有波段的低分辨率卫星的高空间分辨率波段；

5）、利用所述高精度的共有波段的低分辨率卫星的高空间分辨率波段锐化处理非共有波段的低分辨率卫星的低空间分辨率波段，获得非共有波段的低分辨率卫星的高空间分辨率波段；

6）、利用所述非共有波段的低分辨率卫星的高空间分辨率波段计算最大叶绿素指数。

2.根据权利要求1所述的基于多源卫星数据计算内陆水域最大叶绿素指数的方法，其特征在于，所述遥感图像角度数据包括太阳天顶角SOZ、太阳方位角SOA、卫星天顶角SAZ、卫星方位角SAA和相对方位角REA。

3.根据权利要求2所述的基于多源卫星数据计算内陆水域最大叶绿素指数的方法，其特征在于，收集所述高分辨率卫星和低分辨率卫星的波段包括：

1）、筛选高分辨率卫星和低分辨率卫星的遥感影像，剔除有云的遥感图像；

2）、对筛选后的遥感影像进行几何校正；

3）、对校正后的遥感图像进行掩膜处理，以获得研究区域；

4）、通过均值计算，求得高分辨率卫星和低分辨率卫星的研究区域内的相同时间的波段。

4.根据权利要求3所述的基于多源卫星数据计算内陆水域最大叶绿素指数的方法，其特征在于，所述基于所述现场实测的光谱反射率数据和光谱响应函数获得模拟的卫星遥感反射率数据具体为：利用如下公式获得模拟的卫星遥感反射率数据其中，RSR(λ)为卫星的光谱响应函数，Rrs—measured(λ)是现场实测的光谱反射率数据，Rrs(Bi)是模拟的卫星的第i个波段的遥感反射率数据，λm和λn分别是卫星的波段的下限和上限。

5.根据权利要求4所述的基于多源卫星数据计算内陆水域最大叶绿素指数的方法，其特征在于，所述锐化处理包括以下步骤：

1）、对所述高精度的共有波段的低分辨率卫星的高空间分辨率波段进行升尺度，获得升尺度后的共有波段的低分辨率卫星的低空间分辨率波段；

2）、通过线性回归对所述升尺度后的共有波段的低分辨率卫星的低空间分辨率波段和非共有波段的低分辨率卫星的低空间分辨率波段进行回归建模，以获得回归模型；

3）、估算所述回归模型的回归系数；

4）、基于所述回归模型和回归系数，利用所述高精度的共有波段的低分辨率卫星的高空间分辨率波段预测所述非共有波段的低分辨率卫星的低空间分辨率波段的高分辨率变量；

5）、计算低分辨率残差；

6）、利用所述低分辨率残差计算高分辨率残差；

7）、基于所述高分辨率变量和高分辨率残差，获得所述非共有波段的低分辨率卫星的高空间分辨率波段。

6.根据权利要求5所述的基于多源卫星数据计算内陆水域最大叶绿素指数的方法，其特征在于，高分辨率卫星的波段包括空间分辨率为10m的波长为490nm、560nm和665nm的波段和空间分辨率为20m的波长为865nm的波段，低分辨率卫星的波段包括空间分辨率为300m的波长为490nm、560nm、665nm、865nm、681.25nm、708.75nm和753.75nm的波段。

7.根据权利要求6所述的基于多源卫星数据计算内陆水域最大叶绿素指数的方法，其特征在于，所述非共有波段的低分辨率卫星的低空间分辨率波段为681.25nm、708.75nm和

753.75nm。

说明书

基于多源卫星数据计算内陆水域最大叶绿素指数的方法

技术领域

[0001] 本发明属于环保技术领域，涉及一种计算最大叶绿素指数的方法，尤其是一种基于多源卫星数据计算内陆水域最大叶绿素指数的方法。

背景技术

[0002] 蓝藻快速生长聚集形成水华，导致水体生物多样性急剧下降，破坏水体景观和生态系统平衡，给周边经济带来了巨大的障碍，严重影响到区域环境，并制约当地的经济和社会的可持续发展。快速、准确的掌握蓝藻水华分布信息，对水华预防、治理研究尤为重要。MCI（最大叶绿素指数）是指征蓝藻优势型水体叶绿素浓度的重要水色遥感指标，是内陆水体蓝藻水华遥感监测的首选方法。该方法通过叶绿素信号波段（708nm）测量值与基线的正偏距离（>0），可与水体叶绿素浓度建立良好的线性关系，避免了下限阈值不确定带来的困难。在全球各典型湖泊水华监测中应用效果较好。

[0003] 计算MCI用到了3个波段，其中心波长分别为681.25nm、708.75nm、753.75nm。具有该波长范围的水色传感器空间分辨率较低，对于较小湖泊的蓝藻水华监测有一定的局限性。

[0004] 空间分辨率越高、重访周期越短的卫星数据，越有利于小面积零星蓝藻水华区域的提取。波长为681.25nm、708.75nm、753.75nm波段的水色卫星空间分辨率较低，对于小湖泊最大叶绿素指数（MCI）计算，蓝藻水华区域的提取有很大限制。而其他高分辨卫星不含有这三个波段。因此，可以考虑将高分辨率波段与低分辨率波段相融合，以提高所求波段的分辨率。但是，目前提高波段分辨率的方法多是针对单个卫星的不同波段进行的计算，结合不同卫星波段数据，提高其中一个卫星的波段分辨率，计算MCI指数提取湖泊蓝藻水华区域的方法，目前尚无人进行研究。

发明内容

[0005] 为了克服现有技术的缺陷，本发明提出一种基于多源卫星数据计算内陆水域最大叶绿素指数的方法，其结合不同卫星的波段数据，提高其中一个卫星的波段分辨率来计算MCI指数，效果更好。

[0006] 为了实现上述目的，本发明提供如下技术方案：

[0007] 基于多源卫星数据计算内陆水域最大叶绿素指数的方法，其特征在于，包括以下步骤：

[0008] 1）、数据收集：收集现场实测的光谱反射率数据、遥感图像角度数据、高分辨率卫星和低分辨率卫星的波段；

[0009] 2）、基于所述现场实测的光谱反射率数据和光谱响应函数获得模拟的卫星遥感反射率数据；

[0010] 3）、找出高分辨率卫星和低分辨率卫星的共有波段，并对共有波段的低分辨率卫星的低空间分辨率波段进行重采样以获得共有波段的低分辨率卫星的高空间分辨率波段；

[0011] 4）、基于人工神经网络算法，利用所述模拟的卫星遥感反射率数据、遥感图像角度数据、共有波段的高分辨卫星的高空间分辨率波段和重采样获得的共有波段的低分辨率卫星的高空间分辨率波段建立映射关系，并基于所述映射关系获得高精度的共有波段的低分辨率卫星的高空间分辨率波段；

[0012] 5）、利用所述高精度的共有波段的低分辨率卫星的高空间分辨率波段锐化处理非共有波段的低分辨率卫星的低空间分辨率波段，获得非共有波段的低分辨率卫星的高空间分辨率波段；

[0013] 6）、利用所述非共有波段的低分辨率卫星的高空间分辨率波段计算最大叶绿素指数。

[0014] 优选地，所述遥感图像角度数据包括太阳天顶角SOZ、太阳方位角SOA、卫星天顶角SAZ、卫星方位角SAA和相对方位角REA。

[0015] 优选地，收集所述高分辨率卫星和低分辨率卫星的波段包括：

[0016] 1）、筛选高分辨率卫星和低分辨率卫星的遥感影像，剔除有云的遥感图像；

[0017] 2）、对筛选后的遥感影像进行几何校正；

[0018] 3）、对校正后的遥感图像进行掩膜处理，以获得研究区域；

[0019] 4）、通过均值计算，求得高分辨率卫星和低分辨率卫星的研究区域内的相同时间的波段。

[0020] 优选地，所述基于所述现场实测的光谱反射率数据和光谱响应函数获得模拟的卫星遥感反射率数据具体为：利用如下公式获得模拟的卫星遥感反射率数据

[0021]

[0022] 其中，RSR(λ)为卫星的光谱响应函数，Rrs—measured(λ)是现场实测的光谱反射率数据，Rrs(Bi)是模拟的卫星的第i个波段的遥感反射率数据，λm和λn是卫星的波段上下限。

[0023] 优选地，所述锐化处理包括以下步骤：

[0024] 1）、对所述高精度的共有波段的低分辨率卫星的高空间分辨率波段进行升尺度，获得升尺度后的共有波段的低分辨率卫星的低空间分辨率波段；

[0025] 2）、通过线性回归对所述升尺度后的共有波段的低分辨率卫星的低空间分辨率波段和非共有波段的低分辨率卫星的低空间分辨率波段进行回归建模，以获得回归模型；

[0026] 3）、估算所述回归模型的回归系数；

[0027] 4）、基于所述回归模型和回归系数，利用所述高精度的共有波段的低分辨率卫星的高空间分辨率波段预测所述非共有波段的低分辨率卫星的低空间分辨率波段的高分辨率变量；

[0028] 5）、计算低分辨率残差；

[0029] 6）、利用所述低分辨率残差计算高分辨率残差；

[0030] 7）、基于所述高分辨率变量和高分辨率残差，获得所述非共有波段的低分辨率卫星的高空间分辨率波段。

[0031] 优选地，高分辨率卫星的波段包括空间分辨率为10m的波长为490nm、560nm和665nm的波段和空间分辨率为20m的波长为865nm的波段，低分辨率卫星的波段包括空间分辨率为300m的波长为490nm、560nm、665nm、865nm、681.25nm、708.75nm和753.75nm的波段。

[0032] 优选地，所述非共有波段的低分辨率卫星的低空间分辨率波段为681.25nm、708.75nm和753.75nm。

[0033] 与现有技术相比，本发明的基于多源卫星数据计算内陆水域最大叶绿素指数的方法具有如下有益技术效果中的一者或多者：

[0034] 1、其结合不同卫星，也就是，高分辨率卫星和低分辨率卫星的波段数据，通过提高低分辨率卫星的波段分辨率，并利用其计算MCI指数，以此来提取湖泊蓝藻水华区域，使得在内陆水域水华监测中的应用效果更好。

[0035] 2、其在锐化处理时，考虑了残差对数据结果精度的影响，大大减小了误差。

[0036] 3、其在锐化处理时，采用多个高分辨波段作为类全色波段，得到更为精确的回归系数，使得预测的低分辨率卫星的高分辨率波段结果更为准确。

附图说明

[0037] 图1是本发明的基于多源卫星数据计算内陆水域最大叶绿素指数的方法的流程图；

[0038] 图2是锐化处理的流程图；

[0039] 图3是锐化处理过程中的光谱空间自适应算法的流程图。

具体实施方式

[0040] 下面结合附图和实施例对本发明进一步说明，实施例的内容不作为对本发明的保护范围的限制。

[0041] 针对现有技术中内陆水域最大叶绿素指数计算精度低的难题，本发明提供一种基于多源卫星数据计算内陆水域最大叶绿素指数的方法，其结合不同卫星的波段数据，提高计算最大叶绿素指数所需卫星波段的分辨率，效果更好。

[0042] 图1示出了本发明的基于多源卫星数据计算内陆水域最大叶绿素指数的方法的流程图。如图1所示，本发明的基于多源卫星数据计算内陆水域最大叶绿素指数的方法包括以下步骤：

[0043] 一、数据收集。

[0044] 在本发明中，收集的数据包括：现场实测的光谱反射率数据、遥感图像角度数据、高分辨率卫星和低分辨率卫星的波段。

[0045] 其中，为了获得现场实测的光谱反射率数据，需要在研究区（也就是，要进行研究或监测的内陆水域）的每个采样站点，使用便携式地物光谱仪（光谱仪波段范围为350–1050 nm，带宽增量为1.5 nm）在当地时间 10:00 到 14:00 之间采集相关数据。在采集时，采集的相关数据包括参照板、天空、水体的辐射光谱。为了有效避免船舶对水面的干扰和太阳直接辐射的影响，设置水面辐射采集仪器相对于太阳的倾斜角度为φv 135°，天顶角度为θv ~40°。测量水面辐射度后，向上旋转光谱辐射仪，以相同角度测量窗口辐射度。每个目标获~
取十个光谱，消除异常光谱后取平均值。光谱反射率数据计算公式为

[0046] （1）

[0047] 式中，Rrs(λ)为现场实测的光谱反射率数据，Lt为水体辐射度， Lsky为天空辐射度，Lp为灰色参照板的辐射度，ρp为灰色参照板的漫反射率，γ是水表面反射率，取决于风速（无‑1 ‑1风天气为2.2%，风速小于5m s 为2.5%，风速10m s 为2.6%‑2.8%）。

[0048] 遥感图像角度数据是遥感图像预处理中必不可少的参数。遥感图像角度数据包括太阳天顶角、太阳方位角、卫星天顶角、卫星方位角、相对方位角。其中，

[0049] 太阳天顶角SOZ：指像元位置处该时入射太阳光线与该像元地平面法线之间的夹角。太阳天顶角计算公式为：

[0050] cosθ = cosLatcosδcosω+sinLatsinδ （2）

[0051] 式中，θ为太阳天顶角，Lat为像元纬度，δ为太阳赤纬角，ω为太阳时角。在遥感图像角度计算中，纬度是已知的，成像时间也是已知的，太阳赤纬角δ和太阳时角ω计算公式为：

[0052]（3）

[0053] 为日角，由成像时间计算得出。

[0054] ω = 15（AST‑12） （4）

[0055] 式中，AST为像元当地时，根据成像时间的北京时，以及所求像元与北京经度之差计算，4分钟/经度变化，东加西减，去求得该项元位置的AST。

[0056] 太阳方位角SOA：由正北方向起，顺时针转至某个时刻太阳在该像元的入射光线在地表的投影所在的射线。顶点为地表像元点所成的角度。太阳方位角计算公式为

[0057] （5）

[0058] 式中，φ为太阳方位角，θ为太阳天顶角，δ为太阳赤纬角，ω为太阳时角，Lat为像元纬度。

[0059] 卫星天顶角SAZ：指卫星传感器与像元的连线和像元地平面法线所成的夹角，取值范围为0°到90°。

[0060] 卫星天顶角SAZ计算方式：

[0061] （6）

[0062] 式中，SAZ指极轨卫星图像的卫星天顶角，SatAltitude为卫星相对地面高度，单位km，EarthRadius为地球半径，取6378.135km，SatScanAngle为卫星扫描角。卫星扫描角计算公式为：

[0063] （7）

[0064] 式中，central_pixel_index为所成图像中间像元的列标，current_pixel_index为待计算当前像元的列标，FOV为传感器成像的视场角。

[0065] 卫星方位角SAA：指由正北方向起，顺时针转至某个时刻卫星和像元连线在地表的投影所在的射线，顶点为地表像元点，这个旋转的角度称为像元的卫星方位角。

[0066] 卫星方位角SAA计算公式为

[0067] SinSAD = cosEAsinLat+sinEAcosLatcosSAA （8）

[0068] 当SAH<0时，SAA=‑1*SAA

[0069] 式中，SAD为卫星赤纬角，大小为卫星星下点所在纬度，EA为地球张角，SAH为卫星时角。

[0070] 相对方位角REA：为太阳方位角与卫星方位角之差

[0071] REA=|SOZ‑SAZ| （9）

[0072] 如果计算结果REA>180,则REA=360‑REA。

[0073] 为了使得收集的数据更准确，对于高分辨率卫星和低分辨率卫星的数据，需要进行预处理。其中，预处理包括以下步骤：

[0074] 1、筛选高分辨率卫星和低分辨率卫星的遥感影像，剔除有云的遥感图像。由此，能够减少云对遥感图像的干扰。

[0075] 2、对筛选后的遥感影像进行几何校正。由此，能够获得较好的遥感影像。所述几何校正可以采用现有技术中已有的任何校正方法。

[0076] 3、对校正后的遥感图像进行掩膜处理，以获得研究区域。通过掩膜处理，能够去除掉非研究区域，减少数据处理量。

[0077] 4、通过均值计算，求得高分辨率卫星和低分辨率卫星的研究区域内的相同时间的波段。由于高分辨率卫星和低分辨率卫星的重访周期不同，通过均值计算，能够得到同一天的高分辨率卫星和低分辨率卫星的数据。

[0078] 二、基于所述现场实测的光谱反射率数据和光谱响应函数获得模拟的卫星遥感反射率数据。

[0079] 光谱响应函数SRF指传感器在每个波长处，接受的辐亮度与入射的辅亮度的比值。使用光谱响应函数可以使得地物光谱仪测得的高光谱和卫星测得的光谱对应。一般，地物光谱仪每1.5nm（350nm 1050nm）都有对应的光谱反射率值，而卫星通常只有几个波段的光~
谱反射率值，由于光谱响应的特殊性，地物光谱仪的反射率‑波段曲线趋势和卫星的反射率‑波段曲线趋势并非完全一致。将地物光谱仪的光谱反射率值通过光谱响应函数的相关运算处理才可与卫星的光谱反射率值进行比较。

[0080] 利用光谱响应函数，对地物光谱仪的光谱反射率数据进行运算，用于匹配卫星的光谱反射率数据，模拟的卫星遥感反射率数据计算公式如下：

[0081] （10）

[0082] 式中， 为卫星的光谱响应函数， 是现场实测的光谱反射率数据，Rrs(Bi)是模拟的卫星的第i个波段的反射率数据，它是由λm到λn波段范围的根据 RSR积分运算得来，λm和λn表示波段的下限和上限。

[0083] 三、找出高分辨率卫星和低分辨率卫星的共有波段，并对共有波段的低分辨率卫星的低空间分辨率波段进行重采样以获得共有波段的低分辨率卫星的高空间分辨率波段。

[0084] 以具体卫星数据实例，高分辨率卫星选Sentinel‑2（哨兵2号卫星），Sentinel‑2空间分辨率为10m，重访周期为5天，具有空间分辨率高的优点，但其没有连续的波长位于680~750nm的波段，无法利用MCI算法提取蓝藻水华区域。低分辨率卫星选择Sentinel‑3（哨兵3号卫星），Sentinel‑3空间分辨率为300m，重访周期2天，具有重访周期短，波段连续的优点，但空间分辨率较低。

[0085] 选取Sentinel‑2和Sentinel‑3相同波段红光（665nm）、绿光（560nm）、蓝光（490nm）、近红外（865nm）波段作为共同波段。将Sentinel‑3的空间分辨率为300m的490nm、665nm、560nm、490nm和865nm波段重采样为空间分辨率为10m的高分辨率波段。将Sentinel‑
2的20m空间分辨率的865nm波段重采样为10m，作为提高重采样后Sentinel‑2波段的参数。

[0086] 需要说明的是，对于大部分高分辨率卫星，其只有一个高空间分辨率，所以不需要对高分辨率卫星的波段进行重采样。而我们举例选择的Sentinel‑2（哨兵2号卫星），除了具有空间分辨率为10m的波段（665nm、560nm、490nm），还具有空间分辨率为20m的波段（865nm），所以，针对这种特殊情况，也需要对20m空间分辨率的865nm波段重采样为空间分辨率为10m的865nm波段。目的是使得高分辨率卫星的所有波段都是同一个高空间分辨率下的波段。

[0087] 四、基于人工神经网络算法，利用所述模拟的卫星遥感反射率数据、遥感图像角度数据、共有波段的高分辨卫星的高空间分辨率波段和重采样获得的共有波段的低分辨率卫星的高空间分辨率波段建立映射关系，并基于所述映射关系获得高精度的共有波段的低分辨率卫星的高空间分辨率波段。

[0088] 具体地，对于所选择的Sentinel‑2和Sentinel‑3，将300m分辨率的Sentinel‑3（490nm、560nm、665nm、865nm）波段数据直接重采样为10m分辨率，结果还不是特别的理想。因而，在本发明中，结合10m分辨率的Sentinel‑2（490nm、560nm、665nm、865nm）波段数据、模拟的卫星遥感反射率数据、遥感图像角度数据，利用神经网络算法建立的映射关系，进一步提高Sentinel‑3的10m分辨率的490nm、560nm、665nm、865nm波段的精度。

[0089] 借助人工神经网络算法建立Sentinel‑2和Sentinel‑3相同波段（490nm、560nm、665nm、865nm）相同分辨率（10m）数据和模拟的卫星遥感反射率数据的映射关系，提高Sentinel‑3的10m分辨率波段数据的精度，如公式（11）所示。

[0090] Y = ANNs(Sentinel‑2490nm,Sentinel‑2560nm,Sentinel‑2665nm,Sentinel‑2865nm,[0091] Sentinel‑3490nm, Sentinel‑3560nm, Sentinel‑3665nm, Sentinel‑3865nm, angle‑date,Rrs(Bi) （11）

[0092] 式中，ANNs表示人工神经网络，Sentinel‑2490nm表示490nm波段的Sentinel‑2数据，Sentinel‑2560nm表示560nm波段的Sentinel‑2数据，Sentinel‑2665nm表示665nm波段的Sentinel‑2数据，Sentinel‑2865nm表示重采样后的865nm波段的Sentinel‑2数据，Sentinel‑3490nm表示重采样后的490nm波段的Sentinel‑3数据，Sentinel‑3560nm表示重采样后的560nm波段的Sentinel‑3数据， Sentinel‑3665nm表示重采样后的665nm波段的Sentinel‑3数据，Sentinel‑3865nm表示重采样后的865nm波段的Sentinel‑3数据，angle‑date表示遥感图像角度数据，Rrs(Bi)表示根据现场实测的光谱反射率数据模拟的卫星的第i个波段的遥感反射率数据。

[0093] 调整人工神经网络的参数，即可得到高精度的10m分辨率的Sentinel‑3的波段数据。本发明中默认以10m分辨率的Sentinel‑2波段数据，模拟的卫星反射率数据作为特征，利用神经网络提高重采样后的10m分辨率的Sentinel‑3的490nm、560nm、665nm、865nm波段数据的精度，成为高精度的Sentinel‑3的10m分辨率数据，可替代300m分辨率的Sentinel‑3 的490nm、560nm、665nm、865nm波段数据进行后续操作。

[0094] 五、利用所述高精度的共有波段的低分辨率卫星的高空间分辨率波段锐化处理非共有波段的低分辨率卫星的低空间分辨率波段，获得非共有波段的低分辨率卫星的高空间分辨率波段。

[0095] 具体地，利用新建的Sentinel‑3的10m分辨率的490nm、560nm、665nm、865nm波段数据，锐化处理Sentinel‑3的300m分辨率的681.25nm、708.75nm、753.75nm波段数据，得到Sentinel‑3的10m分辨率的681.25nm、708.75nm、753.75nm波段数据。

[0096] 面对点回归克里金法（area‑to‑point regression kriging，ATPRK）可结合全色波段（PAN）实现低分辨波段数据的锐化处理，提高波段的分辨率。ATPRK算法完美地保留低分辨波段数据的光谱特性，克服低空间分辨率波段与PAN波段局部关系难以表征的问题，有很好的应用前景。但考虑到490nm、560nm、665nm、865nm波段和681.25nm、708.75nm、753.75nm波段之间的全局关系，从490nm、560nm、665nm、865nm波段中选择一个 10 m 分辨率PAN 相似波段，无法充分利用490nm波段、560nm波段 、665nm和865nm波段。采用光谱空间自适应面对点回归克里金法（spectral–spatial adaptive ATPRK ，SSAATPRK）可很好解决此问题，四个高空间分辨率波段（即490nm、560nm、665nm、865nm波段）都用作SSAATPRK算法输入，且在计算ATPRK算法时不用为每个低分辨率波段（681.25nm、708.75nm、753.75nm）都匹配类PAN波段。另一方面，在ATPRK算法中，回归模型是基于整幅高、低空间分辨率图像构建的，所有低分辨率像元的回归系数是固定的。单一的全局回归模型可能无法令人满意地处理局部空间变化，低分辨率波段和高分辨率波段之间的关系因像元位置而异。SSAATPRK算法针对这两点进行改进，具体锐化过程如图2所示。

[0097] 具体地，输入Sentinel‑3的300m低空间分辨率波段l(l=681.25nm、708.75nm、753.75nm)数据和Sentinel‑3的10m高空间分辨率波段k(k=490nm、560nm、665nm、865nm)数据。假设 是低空间分辨率波段l以xi为中心的低分辨率像元V的随机变量，i=1,…,N，N是波段l中低分辨率像元的数量。 是高空间分辨率波段k以xj为中心的高分辨率像元v的随机变量。SSAATPRK算法整体流程的最终目的是预测300 m低空间分辨率波段681.25nm、
708.75nm、753.75nm的随机高分辨率变量为 。

[0098] 其具体操作过程如下：

[0099] 1、对所述高精度的共有波段的低分辨率卫星的高空间分辨率波段进行升尺度，获得升尺度后的共有波段的低分辨率卫星的低空间分辨率波段。

[0100] 将高空间分辨率波段k进行升尺度操作，以匹配低空间分辨率波段l，升尺度计算后波段k(k=490nm、560nm、665nm、865nm)和波段l(l=681.25nm、708.75nm、753.75nm)的空间分辨率均为300m。为了方便求得变动的回归系数，针对局部像元块进行升尺度计算，公式如下

[0101] (12)

[0102] 式中，xp表示像元块的位置，高分辨率波段像元块 升尺度为低分辨率波段像元块 ， 是第 l个低分辨率波段的低分辨率像元块的点扩散函数，yp是在计算点扩散函数时引入的量，这里可以理解为一个常数，∗是卷积运算。

[0103] 2、通过线性回归对所述升尺度后的共有波段的低分辨率卫星的低空间分辨率波段和非共有波段的低分辨率卫星的低空间分辨率波段进行回归建模，构建回归模型。

[0104] 具体地，利用升尺度后的类PAN像元块与第l个低分辨率波段的低分辨率像元块之间的关系，构建线性回归函数。通过线性回归对升尺度得到的低分辨率像元块 与低分辨率像元块 之间的关系进行建模，回归模型如下：

[0105] (13)

[0106] 其中， ，al(xp)和bl(xp)是通过最小二乘法估算的两个像元块的线性回归系数。

[0107] 在本发明中，采用光谱空间自适应回归模型进行回归计算，用局部像元块中的共有波段和非共有波段像元进行拟合，选择与非共有低分辨率波段局部像元块相关性较大的490nm或560nm或665nm或865nm波段局部像元块作为类PAN块，如图3所示。利用公式（14）、（15）计算相关系数CClk(xp)和类 PAN 像元块标签δl(xp,k) 。

[0108] 升尺度后的低分辨率像元块 和低分辨率像元块 的相关系数CClk(xp)的计算公式为

[0109] (14)

[0110] 式中，w是xp像元块的大小， 是像元块 中的第i个像元， 是像元块中的第i个像元， 是像元块 的平均值。因此，利用相关系
数CClk(xp)，最终的类 PAN 像元块标签δl(xp,k) 可表示为

[0111] （15）

[0112] 3、估算所述回归模型的回归系数

[0113] 利用上述回归模型，依据已知的高分辨率像元块 与低分辨率像元块可以估算像元块的线性回归系数al(xp)和bl(xp)，基于所述像元块的线性回归系数al(xp)和bl(xp)可以获得整个图像的回归系数al(x)和bl(x)，类PAN像元块标签δl(xp,k) 。不同的像元块，回归系数不同，融合不同像元块的系数al(xp)、bl(xp)和δl(xp,k)，可生成图像的al(x)、bl(x)和δl(x,k)。

[0114] 4、基于所述回归模型和回归系数，利用所述高精度的共有波段的低分辨率卫星的高空间分辨率波段预测所述非共有波段的低分辨率卫星的低空间分辨率波段的高分辨率变量。

[0115] 在本发明中，可以利用公式（16）预测第l个低分辨率波段的高分辨率变量，并将其称为回归运算结果部分。

[0116] (16)

[0117] 式中，al(x)、bl(x)为回归系数，假设al(x)和bl(x)在不同的空间分辨率下通用。将原始高分辨波段 输入，可预测第l个低分辨波段的精细变量 ， 为回归模型运算结果，

[0118] 。

[0119] 5、计算低分辨率残差。

[0120] 所建立的回归模型并不严格适用于所有低分辨率像元块，因此，回归模型通常存在残差。波段l的低分辨率残差 可表示为

[0121] (17)

[0122] 融合不同像元块残差 后生成低分辨率残差 。

[0123] 6、利用所述低分辨率残差计算高分辨率残差。

[0124] ATPK算法可将低分辨率残差 的尺度降低到高分辨残差 的尺度。

[0125] （18）

[0126] 高空间分辨率残差 （在特定位置x0）是由周围低分辨残差线性组合而成。其中 ，λi是波段l 中以xi为中心的第i个低分辨残差像元的权重，N0表示以x0像元为中心的周围的像元，例如N0=5×5像元窗口。通过克里金算法可将预测的误差方差降到最低。

[0127] 7、基于所述高分辨率变量和高分辨率残差，获得所述非共有波段的低分辨率卫星的高空间分辨率波段。

[0128] 通过公式(19)可以估算第l个低分辨率波段的高分辨率随机变量 。

[0129] （19）

[0130] 高分辨率随机变量 预测由两部分组成：回归模型和ATPK模型。其中是回归模型的结果， 是ATPK模型的结果。回归模型充分利用了10m高分辨率波段的纹理信息。由于回归模型得到的高分辨率变量与300 m 低分辨率波段间存在一些误差，使用 ATPK 模型可缩小回归模型的残差。

[0131] 对每个低分辨率波段进行上述处理，即可得到681.25nm、708.75nm、753.75nm波段的高分辨率变量。

[0132] 六、利用所述非共有波段的低分辨率卫星的高空间分辨率波段计算最大叶绿素指数。

[0133] MCI是指征蓝藻优势型水体叶绿素浓度的重要水色遥感指标，是内陆水体蓝藻水华遥感监测的首选方法。针对含蓝藻水体在 650 ～ 750 nm 的反射光谱等特征，利用 2 个端点波段的辐亮度或遥感反射率，构建一条跨 2个端点波长区间的光谱基线，含叶绿素水体在信号波段具有光谱反射峰，其遥感测量值与基线波长处的内插值之差，即为 MCI，与蓝藻优势型水体的叶绿素浓度有很好的正相关。常见的微囊藻、鱼腥藻等原核细胞浮游植物水华，藻细胞在708 ～ 710 nm波长一带有很强的后向散射光谱特性，因而其MCI 遥感信号为正值。MCI 计算式如下:

[0134] (20)

[0135] 式中，MCI指最大叶绿素指数；L1、L2、L3分别指中心波长为λ1、λ2、λ3的辐亮度，其中，λ1=680.25nm,λ2=708.75nm,λ3=753.75nm。

[0136] L1b 级的 OLCI 数据为大气顶部辐亮度产品，因此，基于辐亮度计算得到的MCI具－2 －1 －1 有物理单位( w·m ·sr ·μm ) 。通过公式（21）可将研究区各像元的 OLCI 辐亮度转换为表观反射率:

[0137] （21）

[0138] 式中， rTOA(λ)为表观反射率，LroA(λ)为该波段的辐亮度，E0(λ)为太阳在该波段的光谱辐照度，θ为像元处的太阳天顶角。

[0139] 为比较MCI对研究区水体蓝藻的灵敏程度，将OLCI 影像辐亮度转换为统一的表观反射率后计算。MCI 计算公式如下:

[0140] （22）

[0141] 式中，rTOA(10)、rTOA(11)、rTOA(12)分别为 10（波长为681.25nm）、11（波长为708.75nm）、12（波长为753.75nm） 波段的表观反射率。由于公式(22) 中3个波段的中心波长均位于红光至近红外的红边区间，受到大气分子及气溶胶吸收、散射等的影响程度较为接近，加之MCI的差分计算特性，很好地消除了大气的影响。因此，大气校正处理对计算 MCI 不是必要的步骤，可以直接使用表观反射率计算。

[0142] 利用MCI算法提取研究区蓝藻水华区域，对 MCI 大于零的湖体像元进行分级显示，色阶从灰到白，代表遥感叶绿素信号强度逐渐增高，MCI为负值的湖面显示为黑色，代表表层水体叶绿素浓度极低。

[0143] 基于Sentinel‑3的OLCI 特征波段构建的 MCI 算法可以灵敏高效提取富营养、蓝藻优势型水体中叶绿素浓度状况及空间分布; MCI算法对富营养水体中叶绿素的浓度梯度、差异的检测有很好的动态范围、层次细节、线性响应。构建的高精度的水体叶绿素浓度模型，可精细的表征叶绿素浓度状况和提取蓝藻水华分布，识别叶绿素浓度相对较高的水团或羽流，为蓝藻预警防控、饮用水源安全提供更加精准的信息。

[0144] 最后应当说明的是，以上实施例仅用以说明本发明的技术方案，而非对本发明保护范围的限制。本领域的技术人员，依据本发明的思想，可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的实质和范围。