CN107247935B - 深圳海域初级生产力遥感反演方法及系统 - Google Patents

深圳海域初级生产力遥感反演方法及系统 Download PDF

Info

Publication number
CN107247935B
CN107247935B CN201710392183.XA CN201710392183A CN107247935B CN 107247935 B CN107247935 B CN 107247935B CN 201710392183 A CN201710392183 A CN 201710392183A CN 107247935 B CN107247935 B CN 107247935B
Authority
CN
China
Prior art keywords
ndvi
data
sea area
shenzhen
module
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
CN201710392183.XA
Other languages
English (en)
Other versions
CN107247935A (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.)
Shenzhen Institute of Advanced Technology of CAS
Original Assignee
Shenzhen Institute of Advanced Technology of CAS
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 Shenzhen Institute of Advanced Technology of CAS filed Critical Shenzhen Institute of Advanced Technology of CAS
Priority to CN201710392183.XA priority Critical patent/CN107247935B/zh
Publication of CN107247935A publication Critical patent/CN107247935A/zh
Application granted granted Critical
Publication of CN107247935B publication Critical patent/CN107247935B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V20/00—Scenes; Scene-specific elements
    • G06V20/10—Terrestrial scenes
    • G06V20/13—Satellite images
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q10/00—Administration; Management
    • G06Q10/06—Resources, workflows, human or project management; Enterprise or organisation planning; Enterprise or organisation modelling
    • G06Q10/063—Operations research, analysis or management
    • G06Q10/0639—Performance analysis of employees; Performance analysis of enterprise or organisation operations
    • 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
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Business, Economics & Management (AREA)
  • Human Resources & Organizations (AREA)
  • General Physics & Mathematics (AREA)
  • Physics & Mathematics (AREA)
  • Economics (AREA)
  • Strategic Management (AREA)
  • Entrepreneurship & Innovation (AREA)
  • Theoretical Computer Science (AREA)
  • Development Economics (AREA)
  • Educational Administration (AREA)
  • Multimedia (AREA)
  • Remote Sensing (AREA)
  • Game Theory and Decision Science (AREA)
  • Astronomy & Astrophysics (AREA)
  • Marketing (AREA)
  • Operations Research (AREA)
  • Quality & Reliability (AREA)
  • Tourism & Hospitality (AREA)
  • General Business, Economics & Management (AREA)
  • Image Processing (AREA)

Abstract

本发明涉及一种深圳海域初级生产力遥感反演方法，包括：输入Landsat TM遥感影像数据；计算深圳海域的归一化差分植被指数数据NDVI；将NDVI数据进行线性压缩，得到NDVI′；计算上述得到的NDVI′数据的梯度值
Figure DDA0001307778690000011
在深圳四大主体海域随机布设插值控制点作为一般控制点；从上述四大主体海域上选择
Figure DDA0001307778690000012
值较高的像元作为处理对象Δp，在Δp上随机布设插值控制点作为纹理控制点；对MODIS/AQUA中的数据分别进行反距离加权插值运算；对深圳海域反演时间当日的光照周期数据，进行反距离加权插值运算；利用OVGPM模型计算得到深圳海域初级生产力遥感反演结果。本发明还涉及一种深圳海域初级生产力遥感反演系统。本发明能够更为精准、简易、快速地实现深圳海域初级生产力的反演。

Description

深圳海域初级生产力遥感反演方法及系统
技术领域
本发明涉及一种深圳海域初级生产力遥感反演方法及系统。
背景技术
深圳市位于广东省南部沿海，陆域范围为22°26＇59＂——22°51 ＇49＂N，113°45＇44＂——114°37＇21＂，东临大亚湾和大鹏湾，西濒珠江口和伶仃洋，南隔深圳河与香港相望，北接东莞和惠州；海域连接南海与太平洋，主要划分为珠江口、深圳湾、大鹏湾和大亚湾四个海区，海岸线全长257公里，海域面积1145平方公里。属南亚热带海洋性季风气候，沿岸海区年均表层水温18～28℃，近岸海域主要污染物为石油类、无机氮和活性磷酸盐，近12年共发生赤潮62次，累积发生面积621.85km2，主要发生在2～6月，已成为该海域日益突出的生态问题。通过现有技术方法对深圳市大亚湾和珠江口等海区进行现场调查后发现，叶绿素和初级生产力的周年变化明显呈现双周期型，且珠江口的初级生产力远高于大亚湾。大亚湾叶绿素a与初级生产力在春夏季变化趋势较一致，而大鹏澳海区却没有明显的季节规律，虽然其叶绿素a含量在近岸养殖区较高，但季节上的变化却十分不明显。通过以往相关技术方法的研究可知，深圳市海域的净初级生产力在自然区域上呈现出较大的差别，因而采用同一标准进行深圳市海域的整体研究往往不能得到更为可靠的结果，所以适当划分研究海区，将深圳市海域化整为零不失为一种有效的处理手段。
现有技术对海洋净初级生产力的估算主要集中在其估算模式、时空分布和影响因素等几个方面，相应的技术方法主要有较为传统的黑白瓶法、同位素固碳法、基于海洋卫星遥感的经验模型、半经验模型、 BPM(Bedford Productivity Model)、LPCM(Laboratoirede Physique et Chimie Marines)、VGPM(Vertically Generalized Production Model)等估算模型。
现有的方法大多采用传统的现场测定方法来研究海洋净初级生产力，并且其研究尺度多锁定在大亚湾等局部海湾，采样耗时长且难以做到大面积采样；另一方面，为了解决传统现场测定方法的效率问题，也有一些新的技术方法涌现，譬如利用遥感数据来反演海洋净初级生产力，虽然这些遥感反演技术的时间连续性较好且便于分析大尺度空间内海洋净初级生产力的分布特点，但是以VGPM模型为首的遥感反演手段的适用尺度往往差强人意，在面临小尺度的近海海域问题上往往不能得到很好的结果，再加上不但近海海域上情况复杂，连近海海域的不同海区也可能面临遥感反演上的差异性，除此之外，长时间序列的海洋净初级生产力估算难度较大，参与运算的参数繁多，模型十分复杂，因而对于以深圳海域为例的研究区域而言，现有的技术都存在着一定的局限性。基于以上原因，以往的相似技术方法也多选择东海、太湖、珠江口等区域为主要研究区，而关于深圳近岸海域的相关NPP遥感反演技术方法则较少。
发明内容
有鉴于此，有必要提供一种深圳海域初级生产力遥感反演方法及系统。
本发明提供一种深圳海域初级生产力遥感反演方法，该方法包括如下步骤：a.输入Landsat TM遥感影像数据；b.根据输入的Landsat TM遥感影像数据，计算深圳海域的归一化差分植被指数数据NDVI；c.将计算得到的NDVI数据进行线性压缩，得到NDVI′；d.计算上述得到的NDVI′数据的梯度值
Figure BDA0001307778670000031
e.在深圳的四大主体海域随机布设插值控制点作为一般控制点，其中：所述四大主体海域分别是珠江口海区、深圳湾海区、大鹏湾海区和大亚湾海区；f.按照
Figure BDA0001307778670000032
从高到底的顺序，从上述四大主体海域上选择
Figure BDA0001307778670000033
值较高的像元作为处理对象Δp，并在Δp上随机布设插值控制点作为纹理控制点；g.对输入的MODIS/AQUA数据中的月平均海表温度数据、月平均海水漫衰减系数数据、月平均海水光合有效辐射数据、月平均海水叶绿素a浓度数据，分别进行反距离加权插值运算；h.输入深圳海域反演时间当日的光照周期数据；i.根据步骤g中进行反距离加权插值运算的结果及光照周期数据 ，利用OVGPM模型计算得到深圳海域初级生产力遥感反演结果PPeu。
其中，所述Landsat TM遥感影像数据包括7个谱段：B1蓝绿谱段、 B2绿谱段、B3红谱段、B4近红外谱段、B5近短波红外谱段、B6热红外谱段、B7近短波红外谱段。
所述的步骤b具体包括：
采用如下公式计算深圳海域的归一化差分植被指数数据NDVI：
NDVI＝(B4-B3)/(B4+B3)
其中，NDVI表示归一化差分植被指数；B3表示Landsat TM数据红谱段的亮度值；B4表示Landsat TM数据近红外谱段的亮度值。
所述的步骤c具体包括：
将NDVI数据的数据范围由NDVI∈[-1，1]调整为NDVI∈[0，1]，线性压缩公式为：
NDVIi′＝(χi-χmin)/(χmax-χmin)，i＝1、2、3、…、n
其中，NDVI′表示线性压缩后的NDVI值；χi表示每个像元的NDVI 值；χmin表示所有像元中的NDVI的最小值；χmsx表示所有像元中的 NDVI的最大值；n表示该NDVI数据的像元数量。
所述的步骤d具体包括：
利用公式：
Figure BDA0001307778670000041
dx(i,j)＝NDVI′(i+1,j)-NDVI′(i,j)
dy(i,j)＝NDVI′(i,j+1)-NDVI′(i,j)
其中，
Figure BDA0001307778670000042
表示NDVI′的梯度值；NDVI′表示该像元的NDVI′值；i、j 表示该像元的位置。
本发明提供一种深圳海域初级生产力遥感反演系统，该系统包括该系统包括输入模块、归一化模块、线性压缩模块、梯度值计算模块、控制点模块、插值运算模块以及反演结果模块，其中：所述输入模块用于输入Landsat TM遥感影像数据；所述归一化模块用于根据输入的Landsat TM遥感影像数据，计算深圳海域的归一化差分植被指数数据NDVI；所述线性压缩模块用于将计算得到的NDVI数据进行线性压缩，得到 NDVI′；所述梯度值计算模块用于计算上述得到的NDVI′数据的梯度值
Figure BDA0001307778670000043
所述控制点模块用于在深圳的四大主体海域随机布设插值控制点作为一般控制点，其中：所述四大主体海域分别是珠江口海区、深圳湾海区、大鹏湾海区和大亚湾海区；所述控制点模块还用于按照
Figure BDA0001307778670000044
从高到底的顺序，从上述四大主体海域上选择
Figure BDA0001307778670000045
值较高的像元作为处理对象Δp，并在Δp上随机布设插值控制点作为纹理控制点；所述插值运算模块用于对输入的MODIS/AQUA数据中的月平均海表温度数据、月平均海水漫衰减系数数据、月平均海水光合有效辐射数据、月平均海水叶绿素a浓度数据，分别进行反距离加权插值运算；所述插值运算模块还用于输入深圳海域反演时间当日的光照周期数据；所述反演结果模块用于根据步骤g中进行反距离加权插值运算的结果及光照周期数据 ，利用OVGPM模型计算得到深圳海域初级生产力遥感反演结果PPeu。
其中，所述的Landsat TM遥感影像数据包括7个谱段：B1蓝绿谱段、 B2绿谱段、B3红谱段、B4近红外谱段、B5近短波红外谱段、B6热红外谱段、B7近短波红外谱段。
所述的归一化模块具体用于：
采用如下公式计算深圳海域的归一化差分植被指数数据NDVI：
NDVI＝(B4-B3)/(B4+B3)
其中，NDVI表示归一化差分植被指数；B3表示Landsat TM数据红谱段的亮度值；B4表示Landsat TM数据近红外谱段的亮度值。
所述的线性压缩模块具体用于：
将NDVI数据的数据范围由NDVI∈[-1，1]调整为NDVI∈[0，1]，线性压缩公式为：
NDVIi′＝(χi-χmin)/(χmax-χmin)，i＝1、2、3、…、n
其中，NDVI′表示线性压缩后的NDVI值；χi表示每个像元的NDVI 值；χmin表示所有像元中的NDVI的最小值；χmax表示所有像元中的 NDVI的最大值；n表示该NDVI数据的像元数量。
所述的梯度值计算模块具体用于：
利用公式：
Figure BDA0001307778670000051
dx(i,j)＝NDVI′(i+1,j)-NDVI′(i,j)
dy(i,j)＝NDVI′(i,j+1)-NDVI′(i,j)
其中，
Figure BDA0001307778670000052
表示NDVI′的梯度值；NDVI′表示该像元的NDVI′值；i、j 表示该像元的位置。
本发明针对深圳海域净初级生产力的时空分布特点，以Landsat TM 数据作为插值控制点的选择背景，生成含有小尺度纹理信息的控制点集，在此基础上对MODIS/AQUA卫星数据及其产品进行反距离加权插值操作，实现对粗空间分辨率遥感数据进行降尺度的处理，之后在VGPM模型的基础上对其必要参数进行了优化，构造出了一个在深圳海域尺度上更为精准的OVGPM模型，将只适用于一类水体遥感反演的 MODIS/AQUA数据产品扩展至以深圳近海海域为例的二级水体中去，首次尝试将MODIS/AQUA数据产品和小尺度的海区联系在一起，不仅扩展了MODIS/AQUA数据产品的适用范围，更有效地利用了现有技术手段和资源。除此之外，本发明还将深圳海域区域进行分区反演处理，最大程度的适应深圳海域净初级生产力“西高东低，北强南弱”的分布特点，最终实现更为精准、简易、快速的深圳市海域净初级生产力的反演。本发明反演结果数据精度较高，反演流程简单易操作，虽然所用参数较多，但都是可以通过MODIS/AQUA数据产品直接获取的遥感数据。
附图说明
图1为本发明深圳海域初级生产力遥感反演方法的流程图；
图2为本发明深圳海域初级生产力遥感反演系统的硬件架构图。
具体实施方式
下面结合附图及具体实施例对本发明作进一步详细的说明。
参阅图1所示，是本发明深圳海域初级生产力遥感反演方法较佳实施例的作业流程图。
步骤S1，输入Landsat TM遥感影像数据。其中：
所述Landsat TM遥感影像数据包括7个谱段，分别为：3个可见光谱段、1个近红外谱段、2个近短波红外谱段、1个热红外谱段，所述3 个可见光谱段包括：蓝绿谱段、绿谱段、红谱段；并将所述7个谱段依次命名为：B1(蓝绿谱段)、B2(绿谱段)、B3(红谱段)、B4(近红外谱段)、B5(近短波红外谱段)、B6(热红外谱段)、B7(近短波红外谱段)。
步骤S2，根据输入的Landsat TM遥感影像数据，计算研究海域的归一化差分植被指数数据NDVI。具体而言：
通过Landsat TM遥感影像数据，计算研究海域的归一化差分植被指数数据(NDVI)，计算公式如下：
NDVI＝(B4-B3)/(B4+B3)
其中，NDVI表示归一化差分植被指数；B3表示Landsat TM数据红谱段的亮度值；B4表示Landsat TM数据近红外谱段的亮度值。
步骤S3，将计算得到的NDVI数据进行线性压缩，得到NDVI′。具体而言：
将NDVI数据的数据范围由NDVI∈[-1，1]调整为NDVI∈[0，1]，线性压缩公式为：
NDVI′′＝(χi-χmin)/(χmax-χmin)，i＝1、2、3、…、n
其中，NDVI′表示线性压缩后的NDVI值；χi表示每个像元的NDVI 值；χmin表示所有像元中的NDVI的最小值；χmax表示所有像元中的 NDVI的最大值；n表示该NDVI数据的像元数量。
步骤S4，计算上述得到的NDVI′数据的梯度值
Figure BDA0001307778670000071
具体而言：
计算NDVI′数据的梯度grad(NDVI)，计算公式如下：
Figure BDA0001307778670000072
dx(i,j)＝NDVI′(i+1,j)-NDVI′(i,j)
dy(i,j)＝NDVI′(i,j+1)-NDVI′(i,j)
其中，
Figure BDA0001307778670000073
表示NDVI′的梯度值；NDVI′表示该像元的NDVI′值；i、j 表示该像元的位置。
步骤S5，在深圳的四大主体海域随机布设插值控制点作为一般控制点。其中：所述四大主体海域分别是珠江口海区、深圳湾海区、大鹏湾海区和大亚湾海区。具体而言：
本实施例中，依据深圳海域净初级生产力“西高东低，北强南弱”的分布特点，在梯度数据grad(NDVI)上划定处理区域，处理区域包括了深圳市海域的四大主体海域，分别是珠江口海区、深圳湾海区、大鹏湾海区和大亚湾海区。在划定好的四大主体海域，按照珠江口海区60、深圳湾海区45、大鹏湾海区100、大亚湾海区95的分配方式在各自区域上随机布设共计300个插值控制点来作为一般控制点。
步骤S6，按照从高到底的顺序，从上述四大主体海域上选择
Figure BDA0001307778670000083
值较高的像元作 为处理对象Δp，并在Δp上随机布设插值控制点作为纹理控制点。具体而言：
本实施例中，为了保留小尺度数据上的纹理信息，在划定好的四大主体海域上选择
Figure BDA0001307778670000082
值最大的前10％的像元作为处理对象Δp，并在Δp上随机布设100个插值控制点作为纹理控制点。本实施例中的一般控制点和纹理控制点累计共有400个，所述控制点均为所采用的IDW算法的插值点位。
步骤S7，对输入的MODIS/AQUA数据中的月平均海表温度数据、月平均海水漫衰减系数数据、月平均海水光合有效辐射数据、月平均海水叶绿素a浓度数据，分别进行反距离加权插值运算。具体包括：
步骤S71，对输入的MODIS/AQUA数据中的月平均海表温度数据，进行反距离加权插值运算，得到水体最大光合速率：
首先，输入MODIS/AQUA三级数据产品中的月平均海表温度数据 (SST)，空间分辨率为4km×4km，根据数据头文件将数据亮度值转化为摄氏度(℃)；
然后，将输入的MODIS/AQUA月平均海表温度数据(SST)进行 IDW插值运算，将该数据的空间分辨率由4km×4km提高至30m×30m；
最后，依据公式计算水体最大光合速率
Figure BDA0001307778670000091
Figure BDA0001307778670000092
f(SST)＝-3.27×10-8SST7+3.41×10-6SST6-1.35×10-4SST5
+2.46×10-3SST4-0.0205SST3+0.0617SST2
+0.2749SST+C3
其中，
Figure BDA0001307778670000093
表示水体最大光合速率，以每小时每毫克叶绿素所产生的碳的毫克数计，单位是mg/(mg·h)；SST表示月平均海表温度，单位是摄氏度(℃)；C1、C2、C3为常系数，通常情况下，C1可以取0，C2可以取0，C3可以取1.2958，其均为模型常量与实测数据的拟合值，可以根据用户需求进行调整。
步骤S72，对输入的MODIS/AQUA数据中的月平均海水漫衰减系数数据，进行反距离加权插值运算，得到真光层深度。
首先，输入MODIS/AQUA三级数据产品中的月平均海水漫衰减系数数据(Kd490)，空间分辨率为4km×4km，单位是m-1；
然后，将输入的MODIS/AQUA月平均海水漫衰减系数数据(Kd490) 进行IDW插值运算，将该数据的空间分辨率由4km×4km提高至30m×30m；
最后，依据公式计算真光层深度Zeu：
Zeu＝C1/Kd490
Figure BDA0001307778670000101
ρ＝Rrs(488)/Rrs(547)
其中，Zeu表示真光层深度，单位是m；Kd490表示海水漫衰减系数，单位是m-1；Rrs(488)表示MODIS/AQUA数据488nm波段亮度值； Rrs(547)表示MODIS/AQUA数据547nm波段亮度值；C1表示常系数，通常可以取值2.7834，可以根据用户实测数据进行调整；C2表示常系数，通常可以取值0.0166，可以根据用户实测数据进行调整。
步骤S73，对输入的MODIS/AQUA数据中的月平均海水光合有效辐射数据，进行反距离加权插值运算，得到表面日光合有效辐射。
首先，输入MODIS/AQUA三级数据产品中的月平均海水光合有效辐射数据(PAR)来代替PPeu反演公式中的E0，其空间分辨率为4km×4km，单位是Ein/(m2·d)；
然后，将输入的MODIS/AQUA月平均海水光合有效辐射数据(PAR) 进行IDW插值运算，将该数据的空间分辨率由4km×4km提高至 30m×30m。
步骤S74，对输入的MODIS/AQUA数据中的月平均海水叶绿素a 浓度数据，进行反距离加权插值运算，得到垂直叶绿素a浓度。
首先，输入MODIS/AQUA三级数据产品中的月平均海水叶绿素a 浓度数据(Chl-a)来代替PPeu反演公式中的Copt，其空间分辨率为 4km×4km，单位是mg/m3；
然后，将输入的MODIS/AQUA月平均海水叶绿素a浓度数据(Chl-a) 进行IDW插值运算，将该数据的空间分辨率由4km×4km提高至 30m×30m。
步骤S8，输入研究区域反演时间当日的光照周期数据。具体而言：
输入研究区域反演时间当日的光照周期数据(Dirr)，该数据为离散常数数值，单位是小时(h)，该数据可由当地天文台查询获得。
步骤S9，根据步骤S7中进行反距离加权插值运算的结果及光照周期数据 ，利用OVGPM模型计算得到深圳海域初级生产力遥感反演结果PPeu。具体而言：
利用OVGPM模型计算PPeu，公式如下：
Figure BDA0001307778670000111
其中，PPeu是海水表层到真光层的初级生产力，以每平方米产生的碳的毫克数计，单位是mg/m2；
Figure BDA0001307778670000112
是水体最大光合速率，以每小时每毫克叶绿素所产生的碳的毫克数计，单位是mg/(mg·h)；E0是海水表面日光合有效辐射度，单位是mol/(m2·d)；Zee是真光层深度，单位是米 (m)；Copt是
Figure BDA0001307778670000113
所在深度处的叶绿素a的浓度值，可用表层叶绿素a 的浓度值代替，单位是mg/m3；Dirr是光照周期，单位为小时(h)；α和β为常量参数，α通常取值0.6613，β通常取值4.1102，均为模型常量与实测数据的拟合值。
参阅图2所示，是本发明深圳海域初级生产力遥感反演系统10的硬件架构图。该系统包括：输入模块101、归一化模块102、线性压缩模块103、梯度值计算模块104、控制点模块105、插值运算模块106以及反演结果模块107。
所述输入模块101用于输入Landsat TM遥感影像数据。其中：
所述Landsat TM遥感影像数据包括7个谱段，分别为：3个可见光谱段、1个近红外谱段、2个近短波红外谱段、1个热红外谱段，所述3 个可见光谱段包括：蓝绿谱段、绿谱段、红谱段；并将所述7个谱段依次命名为：B1(蓝绿谱段)、B2(绿谱段)、B3(红谱段)、B4(近红外谱段)、B5(近短波红外谱段)、B6(热红外谱段)、B7(近短波红外谱段)。
所述归一化模块102用于根据输入的Landsat TM遥感影像数据，计算研究海域的归一化差分植被指数数据NDVI。具体而言：
所述归一化模块102通过Landsat TM遥感影像数据，计算研究海域的归一化差分植被指数数据(NDVI)，计算公式如下：
NDVI＝(B4-B3)/(B4+B3)
其中，NDVI表示归一化差分植被指数；B3表示Landsat TM数据红谱段的亮度值；B4表示Landsat TM数据近红外谱段的亮度值。
所述线性压缩模块103用于将计算得到的NDVI数据进行线性压缩，得到NDVI′。具体而言：
所述线性压缩模块103将NDVI数据的数据范围由NDVI∈[-1，1]调整为NDVI∈[0，1]，线性压缩公式为：
NDVIi′＝(χi-χmin)/(χmax-χmin)，i＝1、2、3、…、n
其中，NDVI′表示线性压缩后的NDVI值；χi表示每个像元的NDVI 值；χmin表示所有像元中的NDVI的最小值；χmax表示所有像元中的 NDVI的最大值；n表示该NDVI数据的像元数量。
所述梯度值计算模块104用于计算上述得到的NDVI′数据的梯度值
Figure BDA0001307778670000121
具体而言：
所述梯度值计算模块104计算得到NDVI′数据的梯度grad(NDVI)，计算公式如下：
Figure BDA0001307778670000131
dx(i,j)＝NDVI′(i+1,j)-NDVI′(i,j)
dy(i,j)＝NDVI′(i,j+1)-NDVI′(i,j)
其中，
Figure BDA0001307778670000132
表示NDVI′的梯度值；NDVI′表示该像元的NDVI′值；i、j 表示该像元的位置。
所述控制点模块105用于在深圳的四大主体海域随机布设插值控制点作为一般控制点。其中：所述四大主体海域分别是珠江口海区、深圳湾海区、大鹏湾海区和大亚湾海区。具体而言：
本实施例中，依据深圳海域净初级生产力“西高东低，北强南弱”的分布特点，在梯度数据grad(NDVI)上划定处理区域，处理区域包括了深圳市海域的四大主体海域，分别是珠江口海区、深圳湾海区、大鹏湾海区和大亚湾海区。在划定好的四大主体海域，按照珠江口海区60、深圳湾海区45、大鹏湾海区100、大亚湾海区95的分配方式在各自区域上随机布设共计300个插值控制点来作为一般控制点。
所述控制点模块105还用于按照
Figure BDA0001307778670000133
从高到底的顺序，从上述四大主体海域上选择
Figure BDA0001307778670000134
值较高的像元作为处理对象Δp，并在Δp上随机布设插值控制点作为纹理控制点。具体而言：
本实施例中，为了保留小尺度数据上的纹理信息，在划定好的四大主体海域上选择
Figure BDA0001307778670000135
值最大的前10％的像元作为处理对象Δp，并在Δp上随机布设100个插值控制点作为纹理控制点。本实施例中的一般控制点和纹理控制点累计共有400个，所述控制点均为所采用的IDW算法的插值点位。
所述插值运算模块106用于对输入的MODIS/AQUA数据中的月平均海表温度数据、月平均海水漫衰减系数数据、月平均海水光合有效辐射数据、月平均海水叶绿素a浓度数据，分别进行反距离加权插值运算。具体包括：
对输入的MODIS/AQUA数据中的月平均海表温度数据，进行反距离加权插值运算，得到水体最大光合速率；
对输入的MODIS/AQUA数据中的月平均海水漫衰减系数数据，进行反距离加权插值运算，得到真光层深度；
对输入的MODIS/AQUA数据中的月平均海水光合有效辐射数据，进行反距离加权插值运算，得到表面日光合有效辐射；
对输入的MODIS/AQUA数据中的月平均海水叶绿素a浓度数据，进行反距离加权插值运算，得到垂直叶绿素a浓度。
所述插值运算模块106还用于输入研究区域反演时间当日的光照周期数据。具体而言：
输入研究区域反演时间当日的光照周期数据(Dirr)，该数据为离散常数数值，单位是小时(h)，该数据可由当地天文台查询获得。
所述反演结果模块107用于根据所述插值运算模块106进行反距离加权插值运算的结果及光照周期数据 ，利用OVGPM模型计算得到深圳海域初级生产力遥感反演结果PPeu。具体而言：
利用OVGPM模型计算PPeu，公式如下：
Figure BDA0001307778670000141
其中，PPeu是海水表层到真光层的初级生产力，以每平方米产生的碳的毫克数计，单位是mg/m2；
Figure BDA0001307778670000142
是水体最大光合速率，以每小时每毫克叶绿素所产生的碳的毫克数计，单位是mg/(mg·h)；E0是海水表面日光合有效辐射度，单位是mol/(m2·d)；Zeu是真光层深度，单位是米 (m)；Copt是
Figure BDA0001307778670000143
所在深度处的叶绿素a的浓度值，可用表层叶绿素a 的浓度值代替，单位是mg/m3；Dirr是光照周期，单位为小时(h)；α和β为常量参数，α通常取值0.6613，β通常取值4.1102，均为模型常量与实测数据的拟合值。
本发明对VGPM模型的主要输入参数进行优化，将只适用于一类水体遥感反演的MODIS/AQUA数据产品扩展至以深圳近海海域为例的二级水体中去，首次尝试将MODIS/AQUA数据产品和小尺度的海区联系在一起，不仅扩展了MODIS/AQUA数据产品的适用范围，更有效地利用了现有技术手段和资源，为深圳市近海海域环境监测提供了廉价、高效、简易的技术服务支持。
虽然本发明参照当前的较佳实施方式进行了描述，但本领域的技术人员应能理解，上述较佳实施方式仅用来说明本发明，并非用来限定本发明的保护范围，任何在本发明的精神和原则范围之内，所做的任何修饰、等效替换、改进等，均应包含在本发明的权利保护范围之内。

Claims (10)

1.一种深圳海域初级生产力遥感反演方法，其特征在于，该方法包括如下步骤：
a.输入Landsat TM遥感影像数据；
b.根据输入的Landsat TM遥感影像数据，计算深圳海域的归一化差分植被指数数据NDVI；
c.将计算得到的NDVI数据进行线性压缩，得到NDVI′；
d.计算上述得到的NDVI′数据的梯度值
Figure FDA0002307547810000012
e.在深圳的四大主体海域随机布设插值控制点作为一般控制点，其中：所述四大主体海域分别是珠江口海区、深圳湾海区、大鹏湾海区和大亚湾海区；
f.按照
Figure FDA0002307547810000013
从高到底的顺序，从上述四大主体海域上选择
Figure FDA0002307547810000014
值最大的前10％的像元作为处理对象Δp，并在Δp上随机布设插值控制点作为纹理控制点；
g.对输入的MODIS/AQUA数据中的月平均海表温度数据、月平均海水漫衰减系数数据、月平均海水光合有效辐射数据、月平均海水叶绿素a浓度数据，分别进行反距离加权插值运算；
h.输入深圳海域反演时间当日的光照周期数据；
i.根据步骤g中进行反距离加权插值运算的结果及所述光照周期数据，利用OVGPM模型计算得到深圳海域初级生产力遥感反演结果PPeu；
其中，利用OVGPM模型计算PPeu，公式如下：
Figure FDA0002307547810000011
其中，PPeu是海水表层到真光层的初级生产力，以每平方米产生的碳的毫克数计，单位是mg/m2；
Figure FDA0002307547810000021
是水体最大光合速率，以每小时每毫克叶绿素所产生的碳的毫克数计，单位是mg；E0是海水表面日光合有效辐射度，单位是mol；Zeu是真光层深度，单位是米；Copt是
Figure FDA0002307547810000022
所在深度处的叶绿素a的浓度值，单位是mg/m3；Dirr是光照周期，单位为小时；α和β为常量参数，均为模型常量与实测数据的拟合值。
2.如权利要求1所述的方法，其特征在于，所述Landsat TM遥感影像数据包括7个谱段：B1蓝绿谱段、B2绿谱段、B3红谱段、B4近红外谱段、B5近短波红外谱段、B6热红外谱段、B7近短波红外谱段。
3.如权利要求2所述的方法，其特征在于，所述的步骤b具体包括：
采用如下公式计算深圳海域的归一化差分植被指数数据NDVI：
NDVI＝(B4-B3)/(B4+B3)
其中，NDVI表示归一化差分植被指数；B3表示Landsat TM数据红谱段的亮度值；B4表示Landsat TM数据近红外谱段的亮度值。
4.如权利要求3所述的方法，其特征在于，所述的步骤c具体包括：
将NDVI数据的数据范围由NDVI∈[-1，1]调整为NDVI∈[0，1]，线性压缩公式为：
NDVI′i＝(χi-χmin)/(χmax-χmin)，i＝1、2、3、…、n
其中，NDVI′表示线性压缩后的NDVI值；χi表示每个像元的NDVI值；χmin表示所有像元中的NDVI的最小值；χmax表示所有像元中的NDVI的最大值；n表示该NDVI数据的像元数量。
5.如权利要求4所述的方法，其特征在于，所述的步骤d具体包括：
利用公式：
Figure FDA0002307547810000023
dx(i,j)＝NDVI′(i+1,j)-NDVI′(i,j)
dy(i,j)＝NDVI′(i,j+1)-NDVI′(i,j)
其中，
Figure FDA0002307547810000031
表示NDVI′的梯度值；NDVI′表示该像元的NDVI′值；i、j表示该像元的位置。
6.一种深圳海域初级生产力遥感反演系统，其特征在于，该系统包括输入模块、归一化模块、线性压缩模块、梯度值计算模块、控制点模块、插值运算模块以及反演结果模块，其中：
所述输入模块用于输入Landsat TM遥感影像数据；
所述归一化模块用于根据输入的Landsat TM遥感影像数据，计算深圳海域的归一化差分植被指数数据NDVI；
所述线性压缩模块用于将计算得到的NDVI数据进行线性压缩，得到NDVI′；
所述梯度值计算模块用于计算上述得到的NDVI′数据的梯度值
Figure FDA0002307547810000032
所述控制点模块用于在深圳的四大主体海域随机布设插值控制点作为一般控制点，其中：所述四大主体海域分别是珠江口海区、深圳湾海区、大鹏湾海区和大亚湾海区；
所述控制点模块还用于按照
Figure FDA0002307547810000033
从高到底的顺序，从上述四大主体海域上选择
Figure FDA0002307547810000034
值最大的前10％的像元作为处理对象Δp，并在Δp上随机布设插值控制点作为纹理控制点；
所述插值运算模块用于对输入的MODIS/AQUA数据中的月平均海表温度数据、月平均海水漫衰减系数数据、月平均海水光合有效辐射数据、月平均海水叶绿素a浓度数据，分别进行反距离加权插值运算；
所述插值运算模块还用于输入深圳海域反演时间当日的光照周期数据；
所述反演结果模块用于根据步骤g中进行反距离加权插值运算的结果及所述光照周期数据，利用OVGPM模型计算得到深圳海域初级生产力遥感反演结果PPeu；
其中，利用OVGPM模型计算PPeu，公式如下：
Figure FDA0002307547810000041
其中，PPeu是海水表层到真光层的初级生产力，以每平方米产生的碳的毫克数计，单位是mg/m2；
Figure FDA0002307547810000042
是水体最大光合速率，以每小时每毫克叶绿素所产生的碳的毫克数计，单位是mg；E0是海水表面日光合有效辐射度，单位是mol；Zeu是真光层深度，单位是米；Copt是
Figure FDA0002307547810000043
所在深度处的叶绿素a的浓度值，单位是mg/m3；Dirr是光照周期，单位为小时；α和β为常量参数，均为模型常量与实测数据的拟合值。
7.如权利要求6所述的系统，其特征在于，所述的Landsat TM遥感影像数据包括7个谱段：B1蓝绿谱段、B2绿谱段、B3红谱段、B4近红外谱段、B5近短波红外谱段、B6热红外谱段、B7近短波红外谱段。
8.如权利要求7所述的系统，其特征在于，所述的归一化模块具体用于：
采用如下公式计算深圳海域的归一化差分植被指数数据NDVI：
NDVI＝(B4-B3)/(B4+B3)
其中，NDVI表示归一化差分植被指数；B3表示Landsat TM数据红谱段的亮度值；B4表示Landsat TM数据近红外谱段的亮度值。
9.如权利要求8所述的系统，其特征在于，所述的线性压缩模块具体用于：
将NDVI数据的数据范围由NDVI∈[-1，1]调整为NDVI∈[0，1]，线性压缩公式为：
NDVI′i＝(χi-χmin)/(χmax-χmin)，i＝1、2、3、…、n
其中，NDVI′表示线性压缩后的NDVI值；χi表示每个像元的NDVI值；χmin表示所有像元中的NDVI的最小值；χmax表示所有像元中的NDVI的最大值；n表示该NDVI数据的像元数量。
10.如权利要求9所述的系统，其特征在于，所述的梯度值计算模块具体用于：
利用公式：
Figure FDA0002307547810000051
dx(i，j)＝NDVI′(i+1，j)-NDVI′(i，j)
dy(i，j)＝NDVI′(i，j+1)-NDVI′(i，j)
其中，
Figure FDA0002307547810000052
表示NDVI′的梯度值；NDVI′表示该像元的NDVI′值；i、j表示该像元的位置。
CN201710392183.XA 2017-05-27 2017-05-27 深圳海域初级生产力遥感反演方法及系统 Active CN107247935B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710392183.XA CN107247935B (zh) 2017-05-27 2017-05-27 深圳海域初级生产力遥感反演方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710392183.XA CN107247935B (zh) 2017-05-27 2017-05-27 深圳海域初级生产力遥感反演方法及系统

Publications (2)

Publication Number Publication Date
CN107247935A CN107247935A (zh) 2017-10-13
CN107247935B true CN107247935B (zh) 2020-07-24

Family

ID=60017763

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710392183.XA Active CN107247935B (zh) 2017-05-27 2017-05-27 深圳海域初级生产力遥感反演方法及系统

Country Status (1)

Country Link
CN (1) CN107247935B (zh)

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109359264B (zh) * 2018-05-30 2023-06-16 深圳先进技术研究院 一种基于modis的叶绿素产品降尺度方法及装置
CN110276304A (zh) * 2019-06-25 2019-09-24 北京师范大学 基于降尺度的高分辨率植被生产力遥感估算方法
CN111125277B (zh) * 2019-11-14 2023-04-11 国家气候中心 一种基于立方体技术的Landsat遥感植被指数修复方法
CN111830211B (zh) * 2020-07-30 2021-05-04 中国水产科学研究院南海水产研究所 一种基于rs的海洋初级生产力分布可视化方法
CN114199827B (zh) * 2022-02-21 2022-05-10 中国石油大学(华东) 一种基于遥感数据反演par漫衰减系数垂直变化的方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103513290A (zh) * 2013-10-24 2014-01-15 环境保护部卫星环境应用中心 基于遥感的区域陆地生态系统呼吸监测方法
CN103886213A (zh) * 2014-04-03 2014-06-25 北京大学 一种农作物总初级生产力的遥感估算方法及系统
CN103913421A (zh) * 2012-12-28 2014-07-09 中国科学院沈阳应用生态研究所 一种用于富营养化湖泊水溶态Fe含量测定的方法
CN105115941A (zh) * 2015-09-30 2015-12-02 国家海洋局南海预报中心 一种提取复杂水体叶绿素浓度分布信息的遥感反演方法
CN106446564A (zh) * 2016-09-29 2017-02-22 华南农业大学 一种植被净初级生产力遥感估算方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103913421A (zh) * 2012-12-28 2014-07-09 中国科学院沈阳应用生态研究所 一种用于富营养化湖泊水溶态Fe含量测定的方法
CN103513290A (zh) * 2013-10-24 2014-01-15 环境保护部卫星环境应用中心 基于遥感的区域陆地生态系统呼吸监测方法
CN103886213A (zh) * 2014-04-03 2014-06-25 北京大学 一种农作物总初级生产力的遥感估算方法及系统
CN105115941A (zh) * 2015-09-30 2015-12-02 国家海洋局南海预报中心 一种提取复杂水体叶绿素浓度分布信息的遥感反演方法
CN106446564A (zh) * 2016-09-29 2017-02-22 华南农业大学 一种植被净初级生产力遥感估算方法

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
Estimation of the Primary Productivity in Pearl River Estuary Using MODIS Data;Haibin Ye等;《Estuaries and Coasts》;20150520;第38卷(第2期);506-518 *
Using primar y productivity as an index of coastal eutrophication: a case study in Daya Bay;Huaxue Liu等;《Water and Environment Journal》;20110921;235-240 *
地形调节植被指数构建及在植被覆盖度遥感监测中的应用;江洪等;《福州大学学报 (自然科学版 )》;20100828;第38卷(第4期);527-532 *
基于 VGPM 模型和 MODIS 数据估算梅梁湾浮游植物初级生产力;殷燕等;《生态学报》;20120608;第32卷(第11期);3528-3537 *
基于MODIS数据的福建海域近10年净初级生产力时空变化研究;施益强等;《中国海洋大学学报》;20150930;第45卷(第9期);61-68 *

Also Published As

Publication number Publication date
CN107247935A (zh) 2017-10-13

Similar Documents

Publication Publication Date Title
CN107247935B (zh) 深圳海域初级生产力遥感反演方法及系统
Jia et al. Validation of remotely sensed evapotranspiration over the Hai River Basin, China
CN110750904B (zh) 一种基于遥感数据的区域碳储量空间格局监测系统和方法
CN110276304A (zh) 基于降尺度的高分辨率植被生产力遥感估算方法
CN109344865B (zh) 一种多数据源的数据融合方法
Li et al. Spatiotemporal dynamics of chlorophyll-a in a large reservoir as derived from Landsat 8 OLI data: understanding its driving and restrictive factors
CN107945146B (zh) 一种基于深度卷积神经网络的时空卫星图像融合方法
Jia et al. Validation of Global LAnd Surface Satellite (GLASS) fractional vegetation cover product from MODIS data in an agricultural region
CN109903234B (zh) 一种城市热力景观定量描述及多尺度特征分析方法
Ming-An et al. Validation of satellite-derived sea surface temperatures for waters around Taiwan
Wang et al. Fractional vegetation cover estimation method through dynamic Bayesian network combining radiative transfer model and crop growth model
Chen et al. An all-sky 1 km daily surface air temperature product over mainland China for 2003–2019 from MODIS and ancillary data
CN113139901A (zh) 一种流域尺度植被净初级生产力遥感精细反演方法
Shen et al. Spatial–temporal land-use/land-cover dynamics and their impacts on surface temperature in Chongming Island of Shanghai, China
Keramitsoglou et al. An online system for nowcasting satellite derived temperatures for urban areas
Zhang et al. Phytoplankton size class in the East China Sea derived from MODIS satellite data
CN112580982A (zh) 一种基于多时相遥感和casa模型的生态保护红线实施评估
CN114019579A (zh) 高时空分辨率近地表空气温度重构方法、系统、设备
Li et al. An operational split-window algorithm for retrieving land surface temperature from geostationary satellite data: A case study on Himawari-8 AHI data
Zhu et al. Estimating sunshine duration using hourly total cloud amount data from a geostationary meteorological satellite
Feng et al. Turbidity estimation from GOCI satellite data in the turbid estuaries of China’s coast
Xu et al. A random forest-based data fusion method for obtaining all-weather land surface temperature with high spatial resolution
Kowalewska-Kalkowska et al. Combining satellite imagery and numerical modelling to study the occurrence of warm upwellings in the Southern Baltic Sea in winter
Cao et al. A new 2010 permafrost distribution map over the Qinghai-Tibet Plateau based on subregion survey maps: a benchmark for regional permafrost modeling
Huld et al. Spatial downscaling of 2-meter air temperature using operational forecast data

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