技术领域
[0001] 本
发明属于湖泊
水环境模拟领域，更具体地，涉及一种基于FVCOM和遥感反演的城市湖泊富营养化模拟方法。
背景技术
[0002] 城市湖泊，作为承载着城镇供水、工业用水、土地
灌溉、渔业养殖、景点观赏、蓄洪纳污等多重功能的重要水源之一，由于粗放型经济发展模式、污
水处理系统运转不合理，
水体自净能
力差等原因，富营养化污染日益加重。湖泊的富营养化会造成许多危害，破坏生态景观、降低水体
生态系统的物种多样性和
稳定性、威胁人体健康、造成经济损失等，因此，结合
水体富营养化的因素，选择最合适的研究方法，针对城市湖泊水体富营养化进行研究是十分必要的。
[0003] 城市湖泊富营养化的形成受水域各种
生物和非生物因素的影响，当前主要的研究方法包括地面观测、遥感反演、模型模拟三种方式。地面观测取样成本高、空间代表性差；遥感反演方法普遍存在时间连续性差，且
精度依赖于遥感卫星
传感器性能的缺点；随着数据同化方法、数字流域、高性能数值模拟等现代信息技术的不断进步、渐趋实用化和多学科的交叉融合，湖泊富营养化研究总的趋势表现为从传统的单点直接观测，转向借助空间探测技术的大面积同步遥感反演、再朝着模拟与观测统一的数据同化发展方向。水质参数遥感反演是湖泊富营养化水平遥感监测的主要手段，相对于传统的反演模型，
机器学习反演模型可以更好的模拟遥感影像波段反射率与实测水质数据间复杂的非线性关系，但在内陆湖泊反演研究中，由于水体光学特征复杂，单一反演模型的应用会因为模型存在的局限性致使反演结果体现出极大的不确定特征；现有湖泊富营养化模型存在参数获取不易、复杂边界处理困难、物理机制描述不准确、数值计算方法精度不高等局限性。
[0004] 因此，目前这三种方式都无法精准表达湖泊富营养化演变变量的连续性、非线性、高度
时空异质性和不确定性特征，直接影响富营养化模拟结果的精确性。
发明内容
[0005] 针对
现有技术的以上
缺陷或改进需求，本发明提供了一种基于FVCOM和遥感反演的城市湖泊富营养化模拟方法和系统，其目的在于解决，现有湖泊富营养化研究方法不能精确表达湖泊富营养化演变变量的真实特性，造成模拟结果精确性低的技术问题。
[0006] 为实现上述目的，按照本发明的一个方面，提供了一种基于FVCOM和遥感反演的城市湖泊富营养化模拟方法，包括：
[0007] (1)利用遥感影像数据和多种机器学习
算法构建多个遥感反演模型，利用熵权法对多个遥感反演模型进行有机组合，得到最优遥感反演模型；
[0008] (2)将待模拟湖泊的遥感影像数据输入所述最优遥感反演模型，得到富营养化指标遥感反演值；
[0009] (3)将待模拟湖泊对应的边界条件、水文气象数据、水体参数和富营养化指标初始测量值输入FVCOM三维水动力水质耦合模型，得到富营养化预测值；所述FVCOM三维水动力水质耦合模型由FVCOM模型中的湖泊三维水动力模
块和湖泊三维水质模块耦合得到；
[0010] (4)根据所述富营养化预测值和富营养化指标遥感反演值之间误差，利用卡尔曼滤波算法对所述FVCOM模型参数进行更新并对输出的富营养化预测值进行校正，实现所述FVCOM模型输出与所述最优遥感反演模型输出的同化；
[0011] (5)采用经过校正和同化的FVCOM模型对待模拟湖泊富营养化过程进行动态连续模拟和富营养化预测。
[0012] 进一步地，步骤(1)具体包括：
[0013] (1.1)对待模拟湖泊的富营养化指标与遥感影像数据进行相关性分析，并对相关性最高的波段及波段组合进行标准化处理后，作为用于反演各富营养化指标的最佳波段及波段组合；
[0014] (1.2)利用多种机器学习算法构建最佳波段及波段组合数据与实测富营养化指标之间的遥感反演模型，得到不同机器学习算法对应的多个遥感反演模型；
[0015] (1.3)通过以下公式对不同机器学习算法对应的多个遥感反演模型进行组合，得到最优遥感反演模型：
[0016]
[0017] 其中， 为第i个遥感反演模型模拟结果中第k个富营养化指标的反演值，Wi为第i个遥感反演模型权重， N为遥感反演模型总个数，Di＝1-Hi，Hi为各模型熵值， Pi(k)为第i个遥感反演模型相对误差权重，
ei(k)为第i个遥感反演模型误差， y(k)为第k个富
营养化指标的实测数据。
[0018] 进一步地，所述多种机器学习算法包括多元回归、
随机森林K值近邻、
人工神经网络、
支持向量机回归和随机森林回归。
[0019] 进一步地，步骤(4)所述卡尔曼滤波算法采用的增益矩阵为：
[0020]
[0021] 其中， 为k+1时刻状态变量预测集合的方差，H为观测算子，Rk+1为k+1时刻观测误差协方差矩阵， tr[]为对矩阵主对
角线上的元素求和，N(k)＝V0(k)-H(k)V(k)Q(k-1)HT(k)-βR(k)，M(k)＝H(k)B(k)P(k-1)BT(k)HT(k)，V(k)为k时刻富营养化预测值与遥感反演最优值之间差值；H(k)为k时刻观测算子，Q(k-1)为模型误差协方差矩阵，β≥1；B(k)为系统状态转移矩阵。
[0022] 按照本发明的另一方面提供了一种基于FVCOM和遥感反演的城市湖泊富营养化模拟系统，包括：
[0023] 最优遥感反演模型构建模块，用于根据遥感影像数据和多种机器学习算法构建多个遥感反演模型，利用熵权法对多个遥感反演模型进行有机组合，得到最优遥感反演模型；
[0024] 富营养化指标遥感反演模块，用于将待模拟湖泊的遥感影像数据输入所述最优遥感反演模型，得到富营养化指标遥感反演值；
[0025] 富营养化指标FVCOM预测模块，用于将待模拟湖泊对应的边界条件、水文气象数据、水体参数和富营养化指标初始测量值输入FVCOM三维水动力水质耦合模型，得到富营养化预测值；所述FVCOM三维水动力水质耦合模型由FVCOM模型中的湖泊三维水动力模块和湖泊三维水质模块耦合得到；
[0026] 同化模拟模块，根据所述富营养化预测值和富营养化指标遥感反演值之间误差，利用卡尔曼滤波算法对所述FVCOM模型参数进行更新并对输出的富营养化预测值进行校正，实现所述FVCOM模型输出与所述最优遥感反演模型输出的同化；
[0027] 湖泊富营养化模拟模块，采用经过校正和同化的FVCOM模型对待模拟湖泊富营养化过程进行动态连续模拟和富营养化预测。
[0028] 进一步地，所述最优遥感反演模型构建模块具体包括：
[0029] 相关性分析单元，用于对待模拟湖泊的富营养化指标与遥感影像数据进行相关性分析，并对相关性最高的波段及波段组合进行标准化处理后，作为用于反演各富营养化指标的最佳波段及波段组合；
[0030] 遥感反演模型集构建单元，利用多种机器学习算法构建最佳波段及波段组合数据与实测富营养化指标之间的遥感反演模型，得到不同机器学习算法对应的多个遥感反演模型；
[0031] 遥感反演模型集组合单元，通过以下公式对不同机器学习算法对应的多个遥感反演模型进行组合，得到最优遥感反演模型：
[0032]
[0033] 其中， 为第i个遥感反演模型模拟结果中第k个富营养化指标的反演值，Wi为第i个遥感反演模型权重， N为遥感反演模型总个数，Di＝1-Hi，Hi为各模型熵值， Pi(k)为第i个遥感反演模型相对误差权重，
ei(k)为第i个遥感反演模型误差， y(k)为第k个富
营养化指标的实测数据。
[0034] 进一步地，所述多种机器学习算法包括多元回归、随机森林K值近邻、人工神经网络、支持向量机回归和随机森林回归。
[0035] 进一步地，所述卡尔曼滤波算法采用的增益矩阵为：
[0036]
[0037] 其中， 为k+1时刻状态变量预测集合的方差，H为观测算子，Rk+1为k+1时刻观测误差协方差矩阵， tr[]为对矩阵主对角线上的元素求和，N(k)＝V0(k)-H(k)V(k)Q(k-1)HT(k)-βR(k)，M(k)＝H(k)B(k)P(k-1)BT(k)HT(k)，V(k)为k时刻富营养化预测值与遥感反演最优值之间差值；H(k)为k时刻观测算子，Q(k-1)为模型误差协方差矩阵，β≥1；B(k)为系统状态转移矩阵。
[0038] 总体而言，通过本发明所构思的以上技术方案与现有技术相比，能够取得下列有益效果：
[0039] (1)本发明引入水文预报中集合建模的思想，在考虑不同模型误差信息的
基础上，集成多种水质遥感反演模型，获取水质参数反演结果的最优解，减少单模式水质遥感反演的不确定性影响，提高模型的稳定性、普适性以及水质遥感反演结果的精度，增加了水质遥感反演结果在湖泊富营养化精细化模拟中的可利用性。
[0040] (2)本发明首次将集合水质遥感反演优化值与FVCOM模型通过集合卡尔曼滤波改进算法进行耦合，用于城市湖泊的三维水动力水质模拟，实现数据在时间和空间上的扩展，为FVCOM数值计算提供具有时空一致性、连续性和更为精确的水质观测数据，修正FVCOM模型模拟结果，很好地解决了模型算子和观测算子的非线性、不连续性问题，有效提升无资料、信息不完全情况下富营养化演变动态过程的模拟与预报精度，为减少水环境数值模拟的不确定性提供了有效的途径。
附图说明
[0041] 图1是本发明的一种基于FVCOM和遥感反演的城市湖泊富营养化模拟方法
流程图；
[0044] 图4是东湖三维TP浓度场可视化模拟图。
具体实施方式
[0045] 为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及
实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。
[0046] 参考图1，本发明实施例提供了一种基于FVCOM和遥感反演的城市湖泊富营养化模拟方法，包括：
[0047] (1)利用遥感影像数据和多种机器学习算法构建多个遥感反演模型，利用熵权法对多个遥感反演模型进行有机组合，得到最优遥感反演模型；
[0048] 具体地，步骤(1)包括：
[0049] (1.1)对待模拟湖泊的富营养化指标与遥感影像数据进行相关性分析，并对相关性最高的波段及波段组合进行标准化处理后，作为用于反演各富营养化指标的最佳波段及波段组合；
[0050] 进一步具体地，分析水温、叶绿素、总氮、总磷等湖泊富营养化指标的特征
光谱，基于统计分析选择出Landsat8遥感数据中用于反演各指标浓度的最佳波段及波段组合。
[0051] (1.2)利用多种机器学习算法构建最佳波段及波段组合数据与实测富营养化指标之间的遥感反演模型，得到不同机器学习算法对应的多个遥感反演模型；
[0052] 进一步具体地，本发明实施例利用多元回归、随机森林K值近邻(KNN)、人工神经网络(ANN)、支持向量机回归(SVR)和随机森林回归(RF)构建最佳波段及波段组合数据与实测富营养化指标之间的遥感反演模型。采用上述五种机器学习算法可以不断地学习模型校正过程中的反馈误差，完善自变量与因变量之间复杂的映射关系，是解决内陆水体反演时水质参数与波段反射率间复杂非线性关系的有效方法。
[0053] (1.3)通过以下公式对不同机器学习算法对应的多个遥感反演模型进行组合，得到最优遥感反演模型：
[0054]
[0055] 其中， 为第i个遥感反演模型模拟结果中第k个富营养化指标的反演值，Wi为第i个遥感反演模型权重， N为遥感反演模型总个数，Di＝1-Hi，Hi为各模型熵值， Pi(k)为第i个遥感反演模型相对误差权重，
ei(k)为第i个遥感反演模型误差， y(k)为第k个富
营养化指标的实测数据。
[0056] (2)将待模拟湖泊的遥感影像数据输入最优遥感反演模型，得到富营养化指标遥感反演值；
[0057] (3)将待模拟湖泊对应的边界条件、水文气象数据、水体参数和富营养化指标初始测量值输入FVCOM三维水动力水质耦合模型，得到富营养化预测值；所述FVCOM三维水动力水质耦合模型由FVCOM模型中的湖泊三维水动力模块和湖泊三维水质模块耦合得到；
[0058] 具体地，FVCOM是具有非结构网格的三维有限体积海岸海洋模型，结合了有限单元法和有限差分法的优点，采用有限体积法离散方程，采用模式分裂技术和改进的半隐格式，发展高效稳定的时间推进计算方法，解决摩阻项可能引起的刚性问题，在小尺度问题、拟合复杂岸线边界、并行计算、节省计算花销等方面优势十分明显，已经在国内外河口、港湾、海域、湖泊的水动力水质模拟中得到应用，但是该模型用于内陆城市湖泊研究的案例极少。本发明首次将集合水质遥感反演优化值与FVCOM模型通过集合卡尔曼滤波改进算法进行耦合，用于城市湖泊的三维水动力水质模拟，实现数据在时间和空间上的扩展，为FVCOM数值计算提供具有时空一致性、连续性和更为精确的水质观测数据，修正FVCOM模型模拟结果，很好地解决了模型算子和观测算子的非线性、不连续性问题。
[0059] 本发明实施例收集东湖水下地形DEM数据、水文、气象资料等模型计算所需要的边界条件，运用SMS
软件生成自适应非结构网格；将水文气象数据作为驱动数据；湖泊水体参数作为模型输入参数；选取叶绿素、水温为状态变量，以k时刻对应的叶绿素和水温观测数据作为状态变量的初始值；将伪随机噪声(任意扰动)加入状态变量的初始值，产生状态变量的初始集合；基于增长模繁殖法计算最优扰动场，将最优扰动场
叠加到状态变量，计算得初始扰动场；基于初始扰动场及其它文件，运用FVCOM模型，分别计算叶绿素和水温在k+1时刻的预测值；本发明基于增长模繁殖法(BGM)计算初始扰动场，代替传统集合卡尔曼滤波中蒙特卡罗方法，大大减小了误差增长率。
[0060] (4)根据富营养化预测值和富营养化指标遥感反演值之间误差，利用卡尔曼滤波算法对FVCOM模型参数进行更新并对输出的富营养化预测值进行校正，实现FVCOM模型输出与所述最优遥感反演模型输出的同化；
[0061] 具体地，同化观测值(集合遥感反演数据)与状态变量(叶绿素、水温)，基于方差最小化原理调整模型状态变量，同时考虑模型预报误差和观测值误差，利用优化卡尔曼增益算法对模拟结果进行校正，最终得到优化的模型状态变量，进而更新状态变量；将同化更新的状态变量重新初始化模型参数，进入下一时刻的状态变量模拟，如此推进实现模拟轨迹的调整，最终得到状态变量集合的最优模拟值。
[0062] 本发明对卡尔曼滤波算法中标准增益矩阵进行改进，添加次优渐消因子，实时调整增益矩阵，对历史模拟数据误差进行渐消，以降低历史误差对当前滤波值的影响程度，提高滤波精度。k+1时刻的改进滤波增益矩阵Kk+1公式如下：
[0063]
[0064] 其中， 为k+1时刻状态变量预测集合的方差，为k+1时刻状态变量预测集合的平均
值， H为观测算子，Rk+1为k+1时刻观测误差协方差矩阵，λ(k+1)为k+1时刻的次优渐消因子， tr[]为对矩阵主对角线上的元
素求和，N(k)＝V0(k)-H(k)V(k)Q(k-1)HT(k)-βR(k)，M(k)＝H(k)B(k)P(k-1)BT(k)HT(k)，V(k)为k时刻富营养化预测值与遥感反演最优值之间差值；H(k)为k时刻观测算子，Q(k-1)为模型误差协方差矩阵，β≥1是为使状态估计值更平滑而引入的一个弱化因子；B(k)为系统状态转移矩阵，本发明实施例选择diag[]，对角矩阵。
[0065] (5)采用经过校正和同化的FVCOM模型对待模拟湖泊富营养化过程进行动态连续模拟和富营养化预测。
[0066] 为验证本发明方法有效性，本发明实施例利用2017～2018年5次东湖野外叶绿素实测数据、同步东湖区域Landsat8遥感反射率和气象水文数据构建基于集合卡尔曼滤波的东湖叶绿素同化模型，并开展了基于扰动场生成方法的集合卡尔曼滤波东湖叶绿素同化研究。同时利用已开发的基于3D GIS湖泊富营养化同化模拟系统具有的多源异构数据的一体化管理功能、湖泊富营养化同化模拟模型的可视化前处理和后处理功能，对模拟与预测结果进行三维动态展示。具体方法如下：利用MATLAB程序对基于FVCOM的水动力水质模型生成的东湖流场、叶绿素浓度等状态变量的NetCDF格式数据进行预处理，将NetCDF格式数据转化可被ArcGIS等数据分析及图表绘制软件调用的4D数组文件，以ArcGIS的Geodatabase格式存储和管理，实现多源异构数据的一体化管理、综合查询；将预处理后的状态变量数据送入湖泊三维水动力水质可视化模拟系统，利用空间信息管理与分析技术辅助完成空间数据插值、空间分析，得到湖泊富营养化情况和演变规律，并输出东湖富营养化各状态变量的全湖分布图和动态变化图。现场采样布置图如图2所示，模拟与预测动态可视化结果如图3、图4所示，将本方法的模拟结果与实测流场和富营养指标的记录比较，结果证明，总体水位符合东湖西高东低的地形走势，基本符合东湖表面自由流动的实际情况，富营养指标也比较接近；同时从模拟结果可以看出，东湖整体呈现为轻度、中度富营养化状态，且表现出明显的子湖差异性，结合东湖周边点源污染的分布现状，发现东湖点源污是造成东湖富营养化的主要原因。
[0067] 本发明实施例的另一方面提供了一种基于FVCOM和遥感反演的城市湖泊富营养化模拟系统，包括：最优遥感反演模型构建模块，用于根据遥感影像数据和多种机器学习算法构建多个遥感反演模型，利用熵权法对多个遥感反演模型进行有机组合，得到最优遥感反演模型；富营养化指标遥感反演模块，用于将待模拟湖泊的遥感影像数据输入所述最优遥感反演模型，得到富营养化指标遥感反演值；富营养化指标FVCOM预测模块，用于将待模拟湖泊对应的边界条件、水文气象数据、水体参数和富营养化指标初始测量值输入FVCOM三维水动力水质耦合模型，得到富营养化预测值；所述FVCOM三维水动力水质耦合模型由FVCOM模型中的湖泊三维水动力模块和湖泊三维水质模块耦合得到；同化模拟模块，根据所述富营养化预测值和富营养化指标遥感反演值之间误差，利用卡尔曼滤波算法对所述FVCOM模型参数进行更新并对输出的富营养化预测值进行校正，实现所述FVCOM模型输出与所述最优遥感反演模型输出的同化；湖泊富营养化模拟模块，采用经过校正和同化的FVCOM模型对待模拟湖泊富营养化过程进行动态连续模拟和富营养化预测。
[0068] 其中，最优遥感反演模型构建模块具体包括：相关性分析单元，用于对待模拟湖泊的富营养化指标与遥感影像数据进行相关性分析，并对相关性最高的波段及波段组合进行标准化处理后，作为用于反演各富营养化指标的最佳波段及波段组合；遥感反演模型集构建单元，利用多种机器学习算法构建最佳波段及波段组合数据与实测富营养化指标之间的遥感反演模型，得到不同机器学习算法对应的多个遥感反演模型；遥感反演模型集组合单元，通过以下公式对不同机器学习算法对应的多个遥感反演模型进行组合，得到最优遥感反演模型：
[0069]
[0070] 其中， 为第i个遥感反演模型模拟结果中第k个富营养化指标的反演值，Wi为第i个遥感反演模型权重， N为遥感反演模型总个数，Di＝1-Hi，Hi为各模型熵值， Pi(k)为第i个遥感反演模型相对误差权重，
ei(k)为第i个遥感反演模型误差， y(k)为第k个富
营养化指标的实测数据。
[0071] 卡尔曼滤波算法采用的增益矩阵为：
[0072]
[0073] 其中， 为k+1时刻状态变量预测集合的方差，H为观测算子，Rk+1为k+1时刻观测误差协方差矩阵， tr[]为对矩阵主对角线上的元素求和，N(k)＝V0(k)-H(k)V(k)Q(k-1)HT(k)-βR(k)，M(k)＝H(k)B(k)P(k-1)BT(k)HT(k)，V(k)为k时刻富营养化预测值与遥感反演最优值之间差值；H(k)为k时刻观测算子，Q(k-1)为模型误差协方差矩阵，β≥1；B(k)为系统状态转移矩阵。系统各模块的具体实施方式与上述方法对应，本发明不再赘述。
[0074] 本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何
修改、等同替换和改进等，均应包含在本发明的保护范围之内。
