CN113945527B - 基于卫星数据得到水质总磷参数反演最优模型的方法 - Google Patents

基于卫星数据得到水质总磷参数反演最优模型的方法 Download PDF

Info

Publication number
CN113945527B
CN113945527B CN202111349579.9A CN202111349579A CN113945527B CN 113945527 B CN113945527 B CN 113945527B CN 202111349579 A CN202111349579 A CN 202111349579A CN 113945527 B CN113945527 B CN 113945527B
Authority
CN
China
Prior art keywords
data
water body
total phosphorus
band
wave band
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202111349579.9A
Other languages
English (en)
Other versions
CN113945527A (zh
Inventor
姚千禧
陈燕婕
曹洪涛
顾行发
余涛
刘向东
黄祥志
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Jiangsu Tianhui Spatial Information Research Institute Co ltd
Original Assignee
Jiangsu Tianhui Spatial Information Research Institute Co ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Jiangsu Tianhui Spatial Information Research Institute Co ltd filed Critical Jiangsu Tianhui Spatial Information Research Institute Co ltd
Priority to CN202111349579.9A priority Critical patent/CN113945527B/zh
Publication of CN113945527A publication Critical patent/CN113945527A/zh
Application granted granted Critical
Publication of CN113945527B publication Critical patent/CN113945527B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
    • G01N21/17—Systems in which incident light is modified in accordance with the properties of the material investigated
    • G01N21/25—Colour; Spectral properties, i.e. comparison of effect of material on the light at two or more different wavelengths or wavelength bands
    • G01N21/27—Colour; Spectral properties, i.e. comparison of effect of material on the light at two or more different wavelengths or wavelength bands using photo-electric detection ; circuits for computing concentration
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
    • G01N21/17—Systems in which incident light is modified in accordance with the properties of the material investigated
    • G01N2021/1765—Method using an image detector and processing of image signal
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
    • G01N21/17—Systems in which incident light is modified in accordance with the properties of the material investigated
    • G01N2021/1793—Remote sensing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/02—Reliability analysis or reliability optimisation; Failure analysis, e.g. worst case scenario performance, failure mode and effects analysis [FMEA]
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A20/00—Water conservation; Efficient water supply; Efficient water use
    • Y02A20/152—Water filtration

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Geometry (AREA)
  • Evolutionary Computation (AREA)
  • Computer Hardware Design (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Health & Medical Sciences (AREA)
  • General Engineering & Computer Science (AREA)
  • Chemical & Material Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Biochemistry (AREA)
  • General Health & Medical Sciences (AREA)
  • Immunology (AREA)
  • Pathology (AREA)
  • Investigating Or Analysing Materials By Optical Means (AREA)

Abstract

本发明公开了基于卫星数据得到水质总磷参数反演最优模型的方法，方法包括：步骤S100：按需求提取哨兵二号L2A级数据，对数据进行数据预处理；步骤S200：依据实测点位的经纬度数据对各点位遥感波段的反射率进行提取；步骤S300：将提取到的各标定水体点位的各波段反射率与从现场各水质监测器所测得的总磷浓度数据之间进行相关性分析；从相关性计算结果中筛选出大于相关性阈值的各波段或者各组合波段；步骤S400：基于不同相关性分析分别建立不同回归模型；回归模型包括但不限于单波段回归模型、多波段组合回归模型；步骤S500：对不同回归模型分别进行误差分析，基于误差分析结果确认反演最优模型。

Description

基于卫星数据得到水质总磷参数反演最优模型的方法
技术领域
本发明涉及卫星遥感数据处理技术领域，具体为基于卫星数据得到水质总磷参数反演最优模型的方法。
背景技术
传统的水质总磷监测方法一般采取实地取样分析，然后对样本进行实验室分析；因该方法采取的数据只能代表采集瞬间采集断面的水质情况，难以获取大范围、广区域、准实时的水体水质空间分布情况以及变化趋势，获取的数据往往具有在时空尺度上不连续、范围小、水样采集和分析的数量有限的缺点，所以通过传统的水质总磷监测方法得到的监测结果只具有局部和典型的代表意义，不能满足大尺度以及实时性的水质总磷监测要求；
传统的水质总磷监测方法会消耗大量的人力、物力和财力，且该方法对站点监测仪器的依赖性较高，需要对站点监测仪器进行定期检查维修以免仪器故障造成对总磷浓度的监测影响，所以它还具有操作过程繁琐，耗费大量财力物力的缺点。
发明内容
本发明的目的在于提供基于卫星数据得到水质总磷参数反演最优模型的方法，以解决上述背景技术中提出的问题。
为了解决上述技术问题，本发明提供如下技术方案：基于卫星数据得到水质总磷参数反演最优模型的方法，方法包括：
步骤S100：将所需影像的时间、所需影像的范围作为搜索输入条件通过哨兵网站提取对应的哨兵二号L2A级数据，将提取到的哨兵二号L2A级数据作为元数据集；对元数据集进行数据预处理得到待处理数据集；
步骤S200：结合待处理数据集和元数据集中的数据并依据实测点位的经纬度数据对各点位遥感波段的反射率进行提取；
步骤S300：将提取到的各标定水体点位的各波段反射率与从现场各水质监测器所测得的总磷浓度数据之间进行相关性分析得到相关性计算结果；相关性分析包括但不限于单波段相关性分析和多波段组合相关性分析；从相关性计算结果中筛选出大于相关性阈值的各波段或者各组合波段；
步骤S400：基于不同相关性分析分别建立不同回归模型；回归模型包括但不限于单波段回归模型、多波段组合回归模型；
步骤S500：对不同回归模型分别进行误差分析，基于误差分析结果确认反演最优模型；
本发明运用遥感监测技术，具有高动态、低成本和宏观性等显著特点，在反演总磷的研究中有着传统测量不可替代的优点；遥感监测既能满足大范围的水质总磷反演的需要，也能反映总磷反演在时间上和空间上的变化情况，弥补了单一水面采样的不足，还能发现一些常规方法难以揭示的污染源和污染物的迁移特征；且多光谱遥感监测精度高、波段多、信息量大，大大提高水质总磷参数反演的监测精度，与传统监测方法相比，遥感监测效率更高、方便快速，也减少人们对仪器的检修和维护成本。
进一步的，步骤S100中的数据预处理包括：
步骤S101：选取哨兵二号L2A级数据中分辨率为10m的其中一个波段作为重采样的数据源，使用SANP软件中的Resample功能对元数据集进行重采样得到重采样数据集；对输入输出路径及文件名进行设置，将重采样数据集转换成ENVI格式，得到重采样数据结果；
将哨兵二号L2A级数据中的Band2波段作为重采样的数据源为的是将元数据集中的所有波段的空间分辨率都为Band2波段10米级的分辨率，这样影像分辨率就会变高；
步骤S102：忽略Band8A波段利用ENVI5.3中的LayerStacking功能对重采样数据结果中的其余12个波段数据进行波段融合得到融合图像；
进行波段融合的目的是为了使得到的融合图像具有兼容可见光及近红外波段信息的优势，使得融合图像色彩丰富，地表信息更加清晰；
步骤S103：对融合图像进行NDWI水体指数提取。
进一步的，步骤S200包括：
步骤S201：在Arcgis10.6中对原始水体点位数据、融合图像以及融合图像的NDWI水体指数进行加载；
步骤S202：基于融合图像的NDWI水体指数在Arcgis10.6中利用按点矢量编辑得到水体点位数据，将水体点位数据对原始水体点位数据进行修正并根据修正后得到的水体点位数据在融合图像中对水体点位进行标定得到各标定水体点位；
步骤S203：在属性表中对各标定水体点位进行经纬度计算得到各标定水体点位的经纬度数据；
步骤S204：将各标定水体点位以图像形式导出得到标定点位图像；将原始水体点位数据与各标定水体点位的经纬度数据进行汇总以表格形式进行导出得到经纬度信息汇总表；
步骤S205：在Arcgis10.6中利用多值提取至点功能，对各标定水体点位的各波段反射率进行提取并以表格形式导出得到点位各波段反射率表；将现场各水质监测器所测得的总磷浓度数据进行汇总并以表格形式导出得到实测总磷浓度表格；
由于大部分实测点位的仪器都部署在桥下，若按照原始水体点位经纬度提取反射率的话，提取的反射率大部分会变成桥的反射率，无法提取桥下水体的反射率；且由于融合图像的图像分辨率为10米，凭肉眼画出点位会不精确，因此通过对NDWI水体指数提取，就能区别出陆地、桥梁还有水体，便于得到更加精确的点位数据；
进一步的，步骤S300中的单波段相关性分析是指，利用spss将提取到的各标定水体点位的各波段反射率与从现场各水质监测器所测得的总磷浓度数据之间进行相关性计算，从相关性计算结果中筛选出大于相关性阈值的各波段。
进一步的，多波段组合相关性分析是指，先将提取到的各标定水体点位的各波段之间依次按不同的波段组合方式进行波段组合，再将波段组合的组合反射率与从现场各水质监测器所测得的总磷浓度数据之间进行相关性计算，从相关性计算结果中筛选出大于相关性阈值的各组合波段。
进一步的，波段组合方式包括两波段组合、三波段组合、四波段组合；
其中，两波段组合的组合形式包括bi/bj、bi-bj、bi+bj；三波段组合的组合形式包括：
Figure BDA0003355411540000031
bi/(bj-bk)、(bj-bk)/bi、bi/(bj+bk)、(bj+bk)/bi、
Figure BDA0003355411540000032
四波段组合的组合形式包括：
Figure BDA0003355411540000033
(bi+bj)/(bk-bh)、(bk-bh)/(bi+bj)；其中bi表示标定水体点位的bani波段；bj表示标定水体点位的banj波段；bk表示标定水体点位的bank波段；bh表示标定水体点位的banh波段；i∈(1，2，3，4，5，6，7，8，9，10，11，12)；j∈(1，2，3，4，5，6，7，8，9，10，11，12)；k∈(1，2，3，4，5，6，7，8，9，10，11，12)；h∈(1，2，3，4，5，6，7，8，9，10，11，12)；i，j，k，h出现在一种组合形式时数值不得相同；其中“/”是指将两波段相除，“+”和“-”是指将两波段相加、相减；
重采样前元数据集中有Band1波段、Band2波段、Band3波段、Band4波段、Band5波段、Band6波段、Band7波段、Band8波段、Band8A波段、Band9波段、Band10波段、Band11波段、Band12波段，其中Band8波段和Band8A波段的光谱响应函数范围重叠，中心波长接近，因此在使用时往往根据需要选取一个使用，这里我们选择Band8波段；所以是将12各波段之间依次进行两波段组合、三波段组合、四波段组合；
进一步的，单波段回归模型建立的过程包括：
步骤S401：以从相关性计算结果中筛选出的大于相关性阈值的各波段反射率作为自变量，对应实测总磷浓度作为因变量，建立单波段反演回归模型；
步骤S402：将各标定水体点位进行两种操作，一是将所有标定水体点位用作建模数据生成回归模型，二是选出其中4种点位数据用作验证数据进行回归模型准确度的检验；对大于相关性阈值的各波段建立不同模式的回归模型；回归模型的模式包括线性函数y＝kx+b；对数函数y＝logax，其中a>0且a≠1；逆函数x＝f-1(y)；二次多项式y＝ax^2+bx+c，其中a≠0,b,c,为常数；三次多项式y＝ax^3+bx^2+cx+d，其中，a≠0,b,c,d为常数；复合函数y＝f[g(x)]；幂函数y＝xa；指数函数y＝ax，其中，a为常数且a>0，a≠1；logistic函数
Figure BDA0003355411540000041
步骤S403：基于反演回归模型得到的判定系数确定最优回归模型的模式。
进一步的，多波段组合回归模型建立的过程包括：
步骤S411：以从相关性计算结果中筛选出的大于相关性阈值的各组合波段的反射率作为自变量，对应实测总磷浓度作为因变量，建立多波段反演回归模型；
步骤S412：将各标定水体点位进行两种操作，一是将所有标定水体点位用作建模数据生成回归模型，二是选出其中若干个标定水体点位用作验证数据进行回归模型准确度的检验；对大于相关性阈值的各组合波段建立不同模式的回归模型；回归模型的模式包括线性函数y＝kx+b；对数函数y＝logax，其中a>0且a≠1；逆函数x＝f-1(y)；二次多项式y＝ax^2+bx+c，其中a≠0,b,c,为常数；三次多项式y＝ax^3+bx^2+cx+d，其中，a≠0,b,c,d为常数；复合函数y＝f[g(x)]；幂函数y＝xa；指数函数y＝ax，其中，a为常数且a>0，a≠1；logistic函数
Figure BDA0003355411540000042
步骤S413：基于反演回归模型得到的判定系数确定最优回归模型的模式。
进一步的，步骤S500中误差分析包括对单波段回归模型的误差分析和对多波段组合回归模型的误差分析；
对单波段回归模型进行误差分析时，利用单波段反演回归模型，对用作验证数据用的标定水体点位数据进行计算得到对应的总磷浓度预测值，将总磷浓度预测值与总磷浓度实测值进行误差精度评定指标计算；
对多波段组合回归模型进行误差分析时，利用多波段反演回归模型，对用作验证数据用的标定水体点位数据进行计算得到对应的总磷浓度预测值，将总磷浓度预测值与总磷浓度实测值进行误差精度评定指标计算。
进一步的，误差精度评定指标包括相对误差、判定系数、平均绝对百分比误差、均方根误差；
其中，
Figure BDA0003355411540000051
判定系数表示为R^2；其中，X表示总磷浓度实测值，X′表示总磷浓度预测值，n为样本数目；基于计算得到的误差精度评定指标对不同模型的反演效果进行评定，
与现有技术相比，本发明所达到的有益效果是：本发明依据实测点位的经纬度信息对各点位遥感波段的反射率进行提取，将各波段反射率与实测点位总磷浓度进行相关性分析，通过单波段、多波段组合方式得到最优波段组合并建立回归模型，大大减少了传统水质监测所带来的耗时费力、范围小以及分析数量有限等缺点，且能够快速且有效地实现对水体总磷参数的实时监测。
附图说明
附图用来提供对本发明的进一步理解，并且构成说明书的一部分，与本发明的实施例一起用于解释本发明，并不构成对本发明的限制。在附图中：
图1是本发明基于卫星数据得到水质总磷参数反演最优模型的方法的流程示意图；
图2是哨兵2A卫星搭载推扫式多光谱成像仪(MSI)的13个覆盖光谱波段数据图；
图3是本发明基于卫星数据得到水质总磷参数反演最优模型的方法实施例中常州五牧汇水区域的标定点位图像；
图4是本发明基于卫星数据得到水质总磷参数反演最优模型的方法实施例中原始水体点位数据与各标定水体点位的经纬度数据汇总得到的经纬度信息汇总表；
图5是本发明基于卫星数据得到水质总磷参数反演最优模型的方法实施例中原始水体点位数据与各标定水体点位的经纬度数据汇总得到的经纬度信息汇总表；
图6是本发明基于卫星数据得到水质总磷参数反演最优模型的方法实施例中点位各波段反射率表；
图7是本发明基于卫星数据得到水质总磷参数反演最优模型的方法实施例中实测总磷浓度表格；
图8是本发明基于卫星数据得到水质总磷参数反演最优模型的方法实施例中与总磷浓度相关性最高的单波段的相关性分析结果图；
图9是本发明基于卫星数据得到水质总磷参数反演最优模型的方法中的多波段组合形式图；
图10是本发明基于卫星数据得到水质总磷参数反演最优模型的方法实施例中与总磷浓度相关性最高的多波段组合的相关性分析结果图；
图11是本发明基于卫星数据得到水质总磷参数反演最优模型的方法实施例中Band5波段的回归模型拟合图；
图12是本发明基于卫星数据得到水质总磷参数反演最优模型的方法实施例中Band10波段的回归模型拟合图；
图13是本发明基于卫星数据得到水质总磷参数反演最优模型的方法实施例中单波段反演最优回归模型图；
图14是本发明基于卫星数据得到水质总磷参数反演最优模型的方法实施例中(Band2-Band5)/(Band7/Band3)波段组合回归模型拟合图；
图15是本发明基于卫星数据得到水质总磷参数反演最优模型的方法实施例中Band10波段总磷浓度预测值与实测值的误差分析表；
图16是本发明基于卫星数据得到水质总磷参数反演最优模型的方法实施例中(Band2-Band5)/(Band7/Band3)波段组合总磷浓度预测值与实测值的误差分析表。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
本文中使用术语“哨兵二号L2A级数据”，哨兵2A(Sentinel-2A)卫星由欧洲航天局研制，于2015年6月23日发射升空。该数据具有空间分辨率高、多光谱成像能力强、幅宽宽以及重访周期短等优势，可用于监测地球土地覆盖变化、植被健康和水体污染情况，以及对山体滑坡、洪水等自然灾害进行快速成像，为灾害救援提供帮助；
哨兵2A卫星搭载推扫式多光谱成像仪(MSI)，覆盖13个光谱波段如图2所示；分别包括Band1波段：海岸/气溶胶波段，用来监测近岸水体和大气中的气溶胶；Band2、3、4波段：可见光波段；Band5、6、7波段：红边范围内波段，对监测植被健康信息非常有效；Band8波段：近红外波段；Band8A波段；Band9波段：水蒸气波段；Band10、11、12波段：短波红外波段；
幅宽达290km，重访周期10天，空间分辨率从可见光到近红外不等，最高可达10m，实现了短时间、大范围、高分辨率的对地观测；
哨兵二号L2A级数据空间分辨率高，能够识别河道信息，近红外波段范围(700-800nm)内光谱信息较丰富，尤其特有的红边窄波段位于水体光谱敏感波段范围，能够充分体现不同水环境状况下的水体生物光学特性；
本文中使用术语“SANP”，SNAP是由欧空局开发的一款用于哨兵卫星数据处理的开源软件；
本文使用术语“ENVI5.3”，ENVI5.3是一款十分强悍的遥感图像处理的工具，具有强大的图像信息提取、图像增强、图像计算等功能；
本文中使用术语“NDWI水体指数”，NDWI水体指数是指归一化水指数，用遥感影像的特定波段进行归一化差值处理，以凸显影像中的水体信息；
本文中使用术语“Arcgis10.6”，Arcgis10.6是一款GSI专业的电子地图信息编辑和开发软件，它主要应用于GIS访问，提供一种免费的，快速并且使用简单的方式浏览地理信息，无论是2D还是3D的信息；
本文中使用术语“spss”，spss是世界上最早采用图形菜单驱动界面的统计软件，它最突出的特点就是操作界面极为友好，输出结果美观漂亮；
请参阅图1-图16，本发明提供技术方案：基于卫星数据得到水质总磷参数反演最优模型的方法，方法包括：
步骤S100：将所需影像的时间、所需影像的范围作为搜索输入条件通过哨兵网站提取对应的哨兵二号L2A级数据；
本实施例以常州五牧汇水区域为研究对象，利用哨兵二号L2A级数据，建立与实测总磷浓度相关的单波段回归模型和多波段组合回归模型，并分析它们的的误差以及精度；在本实施例中对2020年2月18日，区域为常州五牧汇水区域的哨兵二号L2A级数据进行提取；对应编号为S2B_MSIL2A_20200218T023749_N0214_R089_T50SQA_20200218T051659；将提取到的哨兵二号L2A级数据作为元数据集；对元数据集进行数据预处理；
将提取到的哨兵二号L2A级数据作为元数据集；对元数据集进行数据预处理；
其中，数据预处理包括：
步骤S101：将哨兵二号L2A级数据中的Band2波段作为重采样的数据源，使用SANP软件中的Resample功能对元数据集进行重采样得到重采样数据集；对输入输出路径及文件名进行设置，将重采样数据集转换成ENVI格式，得到重采样数据结果；
步骤S102：忽略Band8A波段利用ENVI5.3中的LayerStacking功能对重采样数据结果中的其余12个波段数据进行波段融合得到融合图像；
步骤S103：对融合图像进行NDWI水体指数提取；
Figure BDA0003355411540000081
Green指绿波段，NIR是指近红外波段；
步骤S200：依据实测点位的经纬度数据对各点位遥感波段的反射率进行提取；
其中，步骤S200包括：
步骤S201：在Arcgis10.6中对原始水体点位数据、融合图像以及融合图像的NDWI水体指数进行加载；
步骤S202：基于融合图像的NDWI水体指数在Arcgis10.6中利用按点矢量编辑得到水体点位数据，将水体点位数据对原始水体点位数据进行修正并根据修正后得到的水体点位数据在融合图像中对水体点位进行标定得到各标定水体点位；
步骤S203：在属性表中对各标定水体点位进行经纬度计算得到各标定水体点位的经纬度数据；
步骤S204：将各标定水体点位以图像形式导出得到标定点位图像，如图3所示；将原始水体点位数据与各标定水体点位的经纬度数据进行汇总以表格形式进行导出得到经纬度信息汇总表；如图4和图5所示；
步骤S205：在Arcgis10.6中利用多值提取至点功能，对各标定水体点位的各波段反射率进行提取并以表格形式导出得到点位各波段反射率表，如图6所示；将现场各水质监测器所测得的总磷浓度数据进行汇总并以表格形式导出得到实测总磷浓度表格，如图7所示；
步骤S300：将提取到的各标定水体点位的各波段反射率与从现场各水质监测器所测得的总磷浓度数据之间进行相关性分析得到相关性计算结果；相关性分析包括但不限于单波段相关性分析和多波段组合相关性分析；从相关性计算结果中筛选出大于相关性阈值的各波段或者各组合波段；
其中，单波段相关性分析是指，利用spss将提取到的各标定水体点位的各波段反射率与从现场各水质监测器所测得的总磷浓度数据之间进行皮尔逊相关性计算，从相关性计算结果中筛选出大于相关性阈值的各波段；发现与总磷浓度相关性相对较高的波段分别为Band5波段和Band10波段，其相关性分别为-0.552、0.406，而其他波段与总磷浓度数据之间的相关性计算结果没有大于相关性阈值，即与总磷浓度没有显著的相关性；其中相关性高的单波段与总磷浓度的相关性分析结果如图8所示；
其中，多波段组合相关性分析是指，先将提取到的各标定水体点位的各波段之间依次按不同的波段组合方式进行波段组合，再将波段组合的组合反射率与从现场各水质监测器所测得的总磷浓度数据之间进行皮尔逊相关性计算，从相关性计算结果中筛选出大于相关性阈值的各组合波段；
其中，波段组合方式包括两波段组合、三波段组合、四波段组合；如图9所示，两波段组合的组合形式包括：bi/bj、bi-bj、bi+bj；三波段组合的组合形式包括：
Figure BDA0003355411540000091
bi/(bj-bk)、(bj-bk)/bi、bi/(bj+bk)、(bj+bk)/bi、
Figure BDA0003355411540000092
四波段组合的组合形式包括：
Figure BDA0003355411540000093
(bi+bj)/(bk-bh)、(bk-bh)/(bi+bj)；其中bi表示标定水体点位的bani波段；bj表示标定水体点位的banj波段；bk表示标定水体点位的bank波段；bh表示标定水体点位的banh波段；i∈(1，2，3，4，5，6，7，8，9，10，11，12)；j∈(1，2，3，4，5，6，7，8，9，10，11，12)；k∈(1，2，3，4，5，6，7，8，9，10，11，12)；h∈(1，2，3，4，5，6，7，8，9，10，11，12)；i，j，k，h出现在一种组合形式时数值不得相同；
分析各组合形式的相关性发现，与总磷浓度相关性最高的六种组合方式分别为(Band2-Band5)/(Band7/Band3)、(Band2/Band5)/Band5、(Band3/Band7)/(Band11-Band12)、(Band3+Band10)/Band5、(Band1-Band5)/Band7、(Band6+Band10)/(Band12/Band10)，其相关性均达到0.7以上，且(Band2-Band5)/(Band7/Band3)组合方式的相关性最高，为0.933；这六种组合方式的相关性如图10所示；
步骤S400：基于不同相关性分析分别建立不同回归模型；回归模型包括但不限于单波段回归模型、多波段组合回归模型；
其中，单波段回归模型建立的过程包括：
步骤S401：以从相关性计算结果中筛选出的大于相关性阈值的各波段反射率作为自变量，对应实测总磷浓度作为因变量，建立单波段反演回归模型；
步骤S402：将各标定水体点位进行两种操作，一是将所有标定水体点位用作建模数据生成回归模型，二是选出其中4种点位数据用作验证数据进行回归模型准确度的检验；对大于相关性阈值的各波段建立不同模式的回归模型；回归模型的模式包括线性函数y＝kx+b；对数函数y＝logax，其中a>0且a≠1；逆函数x＝f-1(y)；二次多项式y＝ax^2+bx+c，其中a≠0,b,c,为常数；三次多项式y＝ax^3+bx^2+cx+d，其中，a≠0,b,c,d为常数；复合函数y＝f[g(x)]；幂函数y＝xa；指数函数y＝ax，其中，a为常数且a>0，a≠1；logistic函数
Figure BDA0003355411540000101
步骤S403：基于反演回归模型得到的判定系数确定最优回归模型的模式；其中，波段反演最优回归模型及其判定系数R^2如图13所示；Band5波段与总磷浓度的最优回归模型为三次多项式，判定系数R^2为0.250；Band10波段与总磷浓度的最优回归模型为三次多项式，判定系数R^2为0.357；
其中，多波段组合回归模型建立的过程包括：
步骤S411：以从相关性计算结果中筛选出的大于相关性阈值的各组合波段的反射率作为自变量，对应实测总磷浓度作为因变量，建立多波段反演回归模型；通过对多波段组合方式与总磷浓度的相关性分析，得到了最优的多波段组合方式(Band2-Band5)/(Band7/Band3)，并且该波段组合方式较单波段与总磷浓度相关性最高；
步骤S412：将各标定水体点位进行两种操作，一是将所有标定水体点位用作建模数据生成回归模型，二是选出其中4个标定水体点位用作验证数据进行回归模型准确度的检验；对波段组合方式为(Band2-Band5)/(Band7/Band3)的组合波段建立不同模式的回归模型；回归模型的模式包括线性函数y＝kx+b；对数函数y＝logax，其中a>0且a≠1；逆函数x＝f-1(y)；二次多项式y＝ax^2+bx+c，其中a≠0,b,c,为常数；三次多项式y＝ax^3+bx^2+cx+d，其中，a≠0,b,c,d为常数；复合函数y＝f[g(x)]；幂函数y＝xa；指数函数y＝ax，其中，a为常数且a>0，a≠1；logistic函数
Figure BDA0003355411540000111
如图14所示，为(Band2-Band5)/(Band7/Band3)波段组合回归模型拟合图；
步骤S413：基于反演回归模型得到的判定系数确定最优回归模型的模式；波段组合方式(Band2-Band5)/(Band7/Band3)与总磷的最优回归模型为线性，判定系数R^2为0.926，所建立的回归方程如下：y＝1.592x+0.156R2＝0.926；
分析判定系数R^2，发现单波段Band10波段与总磷浓度相关性较高为0.357；多波段组合中(Band2-Band5)/(Band7/Band3)波段组合与实测总磷浓度相关性较高为0.926；多波段组合的相关性明显高于单波段相关性；
步骤S500：利用最佳波段、最优波段组合与实测总磷浓度建立回归模型；对不同回归模型分别进行误差分析，基于误差分析结果确认反演最优模型；
其中，误差分析包括对单波段回归模型的误差分析和对多波段组合回归模型的误差分析；
对单波段回归模型进行误差分析时，利用单波段反演回归模型，对用作验证数据用的标定水体点位数据进行计算得到对应的总磷浓度预测值，将总磷浓度预测值与总磷浓度实测值进行误差精度评定指标计算；其中，误差精度评定指标包括相对误差、判定系数、平均绝对百分比误差MAPE、均方根误差RMSE；其中，
Figure BDA0003355411540000112
Figure BDA0003355411540000113
判定系数表示为R^2；其中，X表示总磷浓度实测值，X′表示总磷浓度预测值，n为样本数目；其中，R^2越接近1，表示数据之间的相关性越强；相对误差、平均绝对百分比误差、均方根误差越小说明模型反演效果越好；
与总磷浓度单波段回归反演模型中，Band10波段三次多项式回归模型较Band5波段更优；利用单波段反演的回归模型，计算出剩余四组数据的总磷浓度预测值，并与总磷浓度实测值进行误差分析，建立单波段回归模型中Band10波段总磷浓度预测值与实测值的误差分析表，如图15所示；图中可以看出，总磷浓度预测值与实测值最大误差值为21.831％，最小误差为0.017％，R^2为0.0194，MAPE为11.366％，RMSE为2.043％；由于总磷浓度预测值与实测值相关性极低，因此，对于建立的单波段反演回归模型，可以知道该反演回归模型的反演结果不能很好的反映对总磷浓度的预测，精度低；
对多波段组合回归模型进行误差时，利用多波段反演回归模型，对用作验证数据用的标定水体点位数据进行计算得到对应的总磷浓度预测值，将总磷浓度预测值与总磷浓度实测值进行误差精度评定指标计算；其中，误差精度评定指标包括相对误差、判定系数、平均绝对百分比误差MAPE、均方根误差RMSE；其中，
Figure BDA0003355411540000121
Figure BDA0003355411540000122
判定系数表示为R^2；其中，X表示总磷浓度实测值，X′表示总磷浓度预测值，n为样本数目；其中，R^2越接近1，表示数据之间的相关性越强；相对误差、平均绝对百分比误差、均方根误差越小说明模型反演效果越好；
与总磷浓度多波段组合回归反演模型中，(Band2-Band5)/(Band7/Band3)波段组合建立的线性回归模型最优，利用多波段组合反演的回归模型，计算出剩余四组数据的总磷浓度预测值，并与总磷浓度实测值进行误差分析，建立多波段组合回归模型中(Band2-Band5)/(Band7/Band3)波段组合总磷浓度预测值与实测值的误差分析表，如图16所示；表中我们看出，总磷浓度预测值与实测值最大误差值为16.641％，最小误差为6.423％，R^2为0.9378，MAPE为12.425％，RMSE为1.688％。因此，对于建立的多波段组合反演回归模型，可以知道该反演回归模型的反演结果可以很好的反映对总磷浓度的预测，精度高；
分析两种模型的误差，发现单波段建立的三次多项式回归模型，其总磷浓度预测值与实测值不相关；多波段组合建立的线性回归模型，无论是总磷浓度预测值与实测值的相关性还是误差分析，该组合模型反演效果良好，反演结果能满足对总磷浓度估测的需求。
需要说明的是，在本文中，诸如第一和第二等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者设备所固有的要素。
最后应说明的是：以上所述仅为本发明的优选实施例而已，并不用于限制本发明，尽管参照前述实施例对本发明进行了详细的说明，对于本领域的技术人员来说，其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (6)

1.基于卫星数据得到水质总磷参数反演最优模型的方法，其特征在于，所述方法包括：
步骤S100：将所需影像的时间和所需影像的范围作为搜索输入条件通过哨兵网站提取对应的哨兵二号L2A级数据，将提取到的哨兵二号L2A级数据作为元数据集；对所述元数据集进行数据预处理得到待处理数据集；
所述步骤S100中的数据预处理包括：
步骤S101：选取哨兵二号L2A级数据中分辨率为10m的其中一个波段作为重采样的数据源，使用SANP软件中的Resample功能对所述元数据集进行重采样得到重采样数据集；对输入输出路径及文件名进行设置，将所述重采样数据集转换成ENVI格式，得到重采样数据结果；
步骤S102：忽略B8A波段利用ENVI5.3中的Layer Stacking功能对所述重采样数据结果中的其余12个波段数据进行波段融合得到融合图像；
步骤S103：对所述融合图像进行NDWI水体指数提取；
步骤S200：结合所述待处理数据集和元数据集中的数据并依据实测点位的经纬度数据对各点位遥感波段的反射率进行提取；
所述步骤S200包括：
步骤S201：在Arcgis10.6中对原始水体点位数据、融合图像以及所述融合图像的NDWI水体指数进行加载；
步骤S202：基于所述融合图像的NDWI水体指数在Arcgis10.6中利用按点矢量编辑得到水体点位数据，将所述水体点位数据对所述原始水体点位数据进行修正并根据修正后得到的水体点位数据在所述融合图像中对水体点位进行标定得到各标定水体点位；
步骤S203：在属性表中对所述各标定水体点位进行经纬度计算得到所述各标定水体点位的经纬度数据；
步骤S204：将各标定水体点位以图像形式导出得到标定点位图像；将所述原始水体点位数据与所述各标定水体点位的经纬度数据进行汇总以表格形式进行导出得到经纬度信息汇总表；
步骤S205：在Arcgis10.6中利用多值提取至点功能，对所述各标定水体点位的各波段反射率进行提取并以表格形式导出得到点位各波段反射率表；将现场各水质监测器所测得的总磷浓度数据进行汇总并以表格形式导出得到实测总磷浓度表格；
步骤S300：将提取到的各标定水体点位的各波段反射率与从现场各水质监测器所测得的总磷浓度数据之间进行相关性分析得到相关性计算结果；所述相关性分析包括单波段相关性分析和多波段组合相关性分析；从所述相关性计算结果中筛选出大于相关性阈值的各波段或者各组合波段；
所述多波段组合相关性分析是指，先将提取到的各标定水体点位的各波段之间依次按不同的波段组合方式进行波段组合，再将波段组合的组合反射率与从现场各水质监测器所测得的总磷浓度数据之间进行相关性计算，从相关性计算结果中筛选出大于相关性阈值的各组合波段；
所述波段组合方式包括两波段组合、三波段组合和四波段组合；
其中，所述两波段组合的组合形式包括bi/bj、bi-bj和bi+bj；所述三波段组合的组合形式包括：
Figure FDA0003859736190000021
bi/(bj-bk)、(bj-bk)/bi、bi/(bj+bk)、(bj+bk)/bi和
Figure FDA0003859736190000022
所述四波段组合的组合形式包括：
Figure FDA0003859736190000023
(bi+bj)/(bk-bh)和(bk-bh)/(bi+bj)；其中bi表示标定水体点位的bandi波段；bj表示标定水体点位的bandj波段；bk表示标定水体点位的bandk波段；bh表示标定水体点位的bandh波段；i∈(1，2，3，4，5，6，7，8，9，10，11，12)；j∈(1，2，3，4，5，6，7，8，9，10，11，12)；k∈(1，2，3，4，5，6，7，8，9，10，11，12)；h∈(1，2，3，4，5，6，7，8，9，10，11，12)；i，j，k，h出现在一种组合形式时数值不得相同；
步骤S400：基于不同相关性分析分别建立不同回归模型；所述回归模型包括单波段回归模型和多波段组合回归模型；通过对多波段组合方式与总磷浓度的相关性分析，得到了最优的多波段组合方式(b2-b5)/(b7/b3)，并且该波段组合方式较单波段与总磷浓度相关性最高；
步骤S500：对不同回归模型分别进行误差分析，基于误差分析结果确认反演最优模型。
2.根据权利要求1所述的基于卫星数据得到水质总磷参数反演最优模型的方法，其特征在于，所述步骤S300中的所述单波段相关性分析是指，利用SPSS将提取到的各标定水体点位的各波段反射率与从现场各水质监测器所测得的总磷浓度数据之间进行相关性计算，从相关性计算结果中筛选出大于相关性阈值的各波段。
3.根据权利要求1所述的基于卫星数据得到水质总磷参数反演最优模型的方法，其特征在于，所述单波段回归模型建立的过程包括：
步骤S401：以从相关性计算结果中筛选出的大于相关性阈值的各波段反射率作为自变量，对应实测总磷浓度作为因变量，建立单波段反演回归模型；
步骤S402：将各标定水体点位进行两种操作，一是将所有标定水体点位用作建模数据生成回归模型，二是选出其中若干个标定水体点位用作验证数据进行回归模型准确度的检验；对大于相关性阈值的各波段建立不同模式的回归模型；所述回归模型的模式包括线性函数、对数函数、逆函数、二次多项式、三次多项式、幂函数、指数函数或logistic函数；
步骤S403：基于反演回归模型得到的判定系数确定最优回归模型的模式。
4.根据权利要求1所述的基于卫星数据得到水质总磷参数反演最优模型的方法，其特征在于，所述多波段组合回归模型建立的过程包括：
步骤S411：以从相关性计算结果中筛选出的大于相关性阈值的各组合波段的反射率作为自变量，对应实测总磷浓度作为因变量，建立多波段反演回归模型；
步骤S412：将各标定水体点位进行两种操作，一是将所有标定水体点位用作建模数据生成回归模型，二是选出其中若干个标定水体点位用作验证数据进行回归模型准确度的检验；对大于相关性阈值的各组合波段建立不同模式的回归模型；所述回归模型的模式包括线性函数、逆函数、二次多项式、三次多项式、幂函数、指数函数或logistic函数；
步骤S413：基于反演回归模型得到的判定系数确定最优回归模型的模式。
5.根据权利要求1所述的基于卫星数据得到水质总磷参数反演最优模型的方法，其特征在于，所述步骤S500中误差分析包括对单波段回归模型的误差分析和对多波段组合回归模型的误差分析；
对单波段回归模型进行误差分析时，利用单波段反演回归模型，对用作验证数据用的标定水体点位数据进行计算得到对应的总磷浓度预测值，将所述总磷浓度预测值与总磷浓度实测值进行误差精度评定指标计算；
对多波段组合回归模型进行误差分析时，利用多波段反演回归模型，对用作验证数据用的标定水体点位数据进行计算得到对应的总磷浓度预测值，将所述总磷浓度预测值与总磷浓度实测值进行误差精度评定指标计算。
6.根据权利要求5所述的基于卫星数据得到水质总磷参数反演最优模型的方法，其特征在于，所述误差精度评定指标包括相对误差、判定系数、平均绝对百分比误差或均方根误差；
其中，
Figure FDA0003859736190000041
Figure FDA0003859736190000042
判定系数表示为R^2；其中，X表示总磷浓度实测值，X'表示总磷浓度预测值，n为样本数目；基于计算得到的误差精度评定指标对不同模型的反演效果进行评定。
CN202111349579.9A 2021-11-15 2021-11-15 基于卫星数据得到水质总磷参数反演最优模型的方法 Active CN113945527B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111349579.9A CN113945527B (zh) 2021-11-15 2021-11-15 基于卫星数据得到水质总磷参数反演最优模型的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111349579.9A CN113945527B (zh) 2021-11-15 2021-11-15 基于卫星数据得到水质总磷参数反演最优模型的方法

Publications (2)

Publication Number Publication Date
CN113945527A CN113945527A (zh) 2022-01-18
CN113945527B true CN113945527B (zh) 2022-11-01

Family

ID=79338196

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111349579.9A Active CN113945527B (zh) 2021-11-15 2021-11-15 基于卫星数据得到水质总磷参数反演最优模型的方法

Country Status (1)

Country Link
CN (1) CN113945527B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115424136A (zh) * 2022-08-31 2022-12-02 中国林业科学研究院资源信息研究所 一种遥感和林相图相结合的林冠健康评价方法及系统
CN116165680A (zh) * 2023-02-14 2023-05-26 中科三清科技有限公司 反演模型更新方法、装置、存储介质与芯片
CN118050329A (zh) * 2024-04-15 2024-05-17 北京四象爱数科技有限公司 一种基于国控断面和哨兵影像的流域水质反演方法及设备

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106525762A (zh) * 2016-11-07 2017-03-22 航天恒星科技有限公司 一种基于自适应模型的水质监测方法和水质监测装置
CN112051222A (zh) * 2020-08-30 2020-12-08 山东锋士信息技术有限公司 一种基于高分卫星影像的河湖水质监测方法
CN112213287B (zh) * 2020-12-07 2021-04-06 速度时空信息科技股份有限公司 基于遥感卫星影像的沿海滩涂盐度反演的方法
CN112881293A (zh) * 2021-01-08 2021-06-01 浙江工商大学 一种基于高分一号卫星的内陆湖泊清洁水体叶绿素a浓度反演方法
CN113420497B (zh) * 2021-06-01 2024-04-19 中国科学院南京地理与湖泊研究所 浑浊湖泊总磷浓度遥感估算方法

Also Published As

Publication number Publication date
CN113945527A (zh) 2022-01-18

Similar Documents

Publication Publication Date Title
CN113945527B (zh) 基于卫星数据得到水质总磷参数反演最优模型的方法
Tewabe et al. Assessing land use and land cover change detection using remote sensing in the Lake Tana Basin, Northwest Ethiopia
Levy et al. Global evaluation of the Collection 5 MODIS dark-target aerosol products over land
Matsushita et al. Integrating remotely sensed data with an ecosystem model to estimate net primary productivity in East Asia
Yang Integrated use of remote sensing and geographic information systems in riparian vegetation delineation and mapping
CN114241331B (zh) 以UAV为地面和Sentinel-2中介的湿地芦苇地上生物量遥感建模方法
CN113436153B (zh) 一种基于高光谱成像和支持向量机技术的原状土壤剖面碳组分预测方法
CN110836870B (zh) 基于gee的大区域湖泊透明度快速制图方法
Fraser et al. UAV and high resolution satellite mapping of Forage Lichen (Cladonia spp.) in a Rocky Canadian Shield Landscape
Ferrando Jorge et al. Measuring soil colour to estimate soil organic carbon using a large-scale citizen science-based approach
Dilbone et al. Spectrally based bathymetric mapping of a dynamic, sand‐bedded channel: Niobrara River, Nebraska, USA
CN113763272A (zh) 一种富营养化湖泊光合有效辐射衰减系数遥感反演方法
CN110596017B (zh) 一种基于空间权重约束和变分自编码特征提取的高光谱影像土壤重金属浓度评估方法
Costa et al. The Nutrient App: Developing a smartphone application for on-site instantaneous community-based NO3 and PO4 monitoring
He et al. Modis aerosol optical thickness product algorithm verification and analysis
Wang et al. Ammonia nitrogen monitoring of urban rivers with UAV-borne hyperspectral remote sensing imagery
Brewer et al. Evaluating generalizations of hydrography in differing terrains for The National Map of the United States
Torbick et al. Assessing invasive plant infestation and disturbance gradients in a freshwater wetland using a GIScience approach
Tong et al. Quantitative monitoring of inland water using remote sensing of the upper reaches of the Huangpu River, China
Li et al. A global historical twice-daily (daytime and nighttime) land surface temperature dataset produced by Advanced Very High Resolution Radiometer observations from 1981 to 2021
Foppa et al. Validation of operational AVHRR subpixel snow retrievals over the European Alps based on ASTER data
CN115389383B (zh) 一种高精度的近岸海域悬浮泥沙浓度的反演方法
Behn et al. Mapping forest cover, Kimberley region of Western Australia
Mittermaier et al. Using feature-based verification methods to explore the spatial and temporal characteristics of the 2019 chlorophyll-a bloom season in a model of the European Northwest Shelf
Qiao et al. A new combined atmospheric correction algorithm for GOCI-2 Data over coastal waters assessed by long-term satellite ocean color platforms

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant