CN116185616A - 一种fy-3d mersi l1b数据自动化再处理方法 - Google Patents
一种fy-3d mersi l1b数据自动化再处理方法 Download PDFInfo
- Publication number
- CN116185616A CN116185616A CN202310095470.XA CN202310095470A CN116185616A CN 116185616 A CN116185616 A CN 116185616A CN 202310095470 A CN202310095470 A CN 202310095470A CN 116185616 A CN116185616 A CN 116185616A
- Authority
- CN
- China
- Prior art keywords
- data
- pixel
- mersi
- radiation
- 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.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 41
- 238000012958 reprocessing Methods 0.000 title claims abstract description 26
- 238000012545 processing Methods 0.000 claims abstract description 56
- 238000002310 reflectometry Methods 0.000 claims abstract description 20
- 238000004364 calculation method Methods 0.000 claims abstract description 18
- 238000001514 detection method Methods 0.000 claims abstract description 15
- 238000005516 engineering process Methods 0.000 claims abstract description 15
- 230000003595 spectral effect Effects 0.000 claims abstract description 12
- 238000006243 chemical reaction Methods 0.000 claims abstract description 8
- 238000012544 monitoring process Methods 0.000 claims abstract description 8
- 230000005855 radiation Effects 0.000 claims description 65
- 238000012937 correction Methods 0.000 claims description 26
- 238000001228 spectrum Methods 0.000 claims description 25
- 238000013139 quantization Methods 0.000 claims description 7
- 230000002159 abnormal effect Effects 0.000 claims description 6
- PXFBZOLANLWPMH-UHFFFAOYSA-N 16-Epiaffinine Natural products C1C(C2=CC=CC=C2N2)=C2C(=O)CC2C(=CC)CN(C)C1C2CO PXFBZOLANLWPMH-UHFFFAOYSA-N 0.000 claims description 5
- 101100096319 Drosophila melanogaster Spc25 gene Proteins 0.000 claims description 5
- 238000007781 pre-processing Methods 0.000 claims description 4
- 101100382321 Caenorhabditis elegans cal-1 gene Proteins 0.000 claims description 3
- 238000004458 analytical method Methods 0.000 claims description 3
- 101150014174 calm gene Proteins 0.000 claims description 3
- 238000013507 mapping Methods 0.000 claims description 3
- 210000001747 pupil Anatomy 0.000 claims description 3
- 238000005520 cutting process Methods 0.000 claims description 2
- 238000007689 inspection Methods 0.000 claims description 2
- 230000005540 biological transmission Effects 0.000 abstract description 7
- 238000011161 development Methods 0.000 abstract description 4
- 230000008569 process Effects 0.000 abstract description 4
- 239000007789 gas Substances 0.000 description 6
- 230000008859 change Effects 0.000 description 5
- 238000009826 distribution Methods 0.000 description 4
- 238000011282 treatment Methods 0.000 description 4
- 230000018109 developmental process Effects 0.000 description 3
- 238000003860 storage Methods 0.000 description 3
- 101000573444 Gallus gallus Multiple inositol polyphosphate phosphatase 1 Proteins 0.000 description 2
- 230000032683 aging Effects 0.000 description 2
- 230000008878 coupling Effects 0.000 description 2
- 238000010168 coupling process Methods 0.000 description 2
- 238000005859 coupling reaction Methods 0.000 description 2
- 238000007405 data analysis Methods 0.000 description 2
- 230000006870 function Effects 0.000 description 2
- 238000012417 linear regression Methods 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000012546 transfer Methods 0.000 description 2
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 1
- 206010063385 Intellectualisation Diseases 0.000 description 1
- 102100029469 WD repeat and HMG-box DNA-binding protein 1 Human genes 0.000 description 1
- 101710097421 WD repeat and HMG-box DNA-binding protein 1 Proteins 0.000 description 1
- 238000010521 absorption reaction Methods 0.000 description 1
- 239000005427 atmospheric aerosol Substances 0.000 description 1
- 239000005441 aurora Substances 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000013500 data storage Methods 0.000 description 1
- 238000013135 deep learning Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 238000007667 floating Methods 0.000 description 1
- 239000005431 greenhouse gas Substances 0.000 description 1
- 238000011221 initial treatment Methods 0.000 description 1
- 239000005433 ionosphere Substances 0.000 description 1
- 238000007726 management method Methods 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 238000003032 molecular docking Methods 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 238000013441 quality evaluation Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000012876 topography Methods 0.000 description 1
- 238000011144 upstream manufacturing Methods 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Chemical compound O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F9/00—Arrangements for program control, e.g. control units
- G06F9/06—Arrangements for program control, e.g. control units using stored programs, i.e. using an internal store of processing equipment to receive or retain programs
- G06F9/46—Multiprogramming arrangements
- G06F9/50—Allocation of resources, e.g. of the central processing unit [CPU]
- G06F9/5005—Allocation of resources, e.g. of the central processing unit [CPU] to service a request
- G06F9/5027—Allocation of resources, e.g. of the central processing unit [CPU] to service a request the resource being a machine, e.g. CPUs, Servers, Terminals
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
- G06F16/10—File systems; File servers
- G06F16/16—File or folder operations, e.g. details of user interfaces specifically adapted to file systems
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
- G06F16/10—File systems; File servers
- G06F16/17—Details of further file system functions
- G06F16/172—Caching, prefetching or hoarding of files
-
- 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)
- Theoretical Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Software Systems (AREA)
- Data Mining & Analysis (AREA)
- Databases & Information Systems (AREA)
- Human Computer Interaction (AREA)
- Image Processing (AREA)
- Radar Systems Or Details Thereof (AREA)
Abstract
本发明针对风云三号D(FY‑3D)气象卫星接收和预处理的MERSIL1B轨道观测数据，提出了一种面向FY‑3DMERSIL1B数据的自动化再处理方法，其采用Python开发平台和GDAL栅格处理技术，集约管理数据检测传输、辐射定标、光谱反射率/亮温计算、等经纬投影、区域裁剪、数据格式转换及归档等模块，实现各处理模块紧密衔接和自动化运行。处理结果表明，本方法具有数据读取效率高和处理耗时少的优点，能自动化加工和归档1km和250m空间分辨率等经纬网格多通道图像，图像质量与原处理技术的图像基本一致，能够应用于卫星遥感监测系统中，从而有效支撑本区域气象遥感服务业务。
Description
技术领域
本发明涉及FY-3D MERSI L1B数据再处理技术方法领域，具体涉及一种FY-3DMERSI L1B数据自动化再处理方法。
背景技术
风云三号D星(FY-3D)于2017年11月在太原卫星发射中心发射升空，是我国第二代极轨气象卫星，与风云三号C星(FY-3C)、风云三号E星(FY-3E)组成我国极轨气象卫星黎明、上午、下午观测业务网络，在每6小时内为数值预报同化业务提供一次全球完整覆盖的遥感观测资料。FY-3D搭载了第二代中分辨率光谱成像仪(MERSI-II)、第二代微波温度计(MWTS-II)、微波成像仪(MWRI)、第二代微波湿度计(MWHS-II)、空间环境检测仪(SEM)、全球导航卫星掩星探测仪(GNOS)、红外高光谱大气探测仪(HIRAS)、近红外高光谱温室气体监测仪(GAS)、广角极光成像仪(WAI-I)、电离层光度计(IPM)共10套遥感观测仪器，观测数据广泛应用于空间天气、大气辐射、陆表生态、数值预报等领域。MERSI-II广泛用于区域大气气溶胶、云、火情监测和陆表生态的监测评估产品制作，目前国内已经开展了基于FY-3D MERSI数据的遥感应用研究，实现了基于FY-3D MERSI-II远红外数据火情监测研究(郑伟、陈洁等，2020)和AOD遥感反演(陈辉等，2022)，基于深度学习研究了FY-3D MERSI的云检测模型(瞿建华等，2019)，针对FY-3D MERSI NDVI产品进行了质量评估(王圆圆等，2022)。省级气象部门主要采用FY3D MERSI-II数据产品开展本地气象遥感实时监测服务业务。
因FY-3D MERSI L1B数据没有等经纬投影信息，光谱波段数据主要以整型辐射亮度灰度值(Digital Number,DN)存储，非遥感技术专业人员解析及应用该数据仍需对数据做辐射定标、像元反射率/亮温计算、重投影、几何校正、区域裁剪和图像格式转等再处理加工，输出易于理解和应用的多通道栅格遥感图像数据。目前本地业务上，FY-3D MERSI L1B数据的再处理加工方式需大量人工干预和设置，没有形成数据处理的自动化处理工程模式，再处理效率不高，不能满足气象遥感监测业务智能化发展要求。
发明内容
针对上述存在的问题，本发明提供一种FY-3D MERSI L1B数据自动化再处理方法，基于网络文件系统传输协议(NetworkFile System，NFS)，能够与风云三省级利用站预处理平台数据传输进行对接，采用Python和GDAL栅格技术对轨道数据进行再处理，实现FY-3DMERSI L1B数据文件检测、元数据和地理数据解析、辐射定标、光谱反射率/亮温计算、重投影、几何校正和区域裁剪等主要模块耦合连接和触发式异步自动运行，提高FY-3D MERSI数据处理效率和自动化水平，自动生成和归档250m和1km两种空间分辨率多通道栅格格式遥感图像。
实现本发明目的的技术解决方案为：
一种FY-3D MERSI L1B数据自动再处理方法，其特征在于，包括以下步骤：
步骤1：通过风云三省级直收站实时接收并预处理为FY-3D MERSI L1B过境区域轨道观测数据；
步骤2：自动检测新接收的FY-3D MERSI L1B数据，并进行自动完整性检查和读取；
步骤3：解析检查后的观测数据属性，得到数据文件光谱波段像元的元数据信息；对元数据信息进行特征分析得到观测角、订正系数以及地理坐标信息；
步骤4：基于得到的观测角、订正系数计算光谱波段通道的大气表观反射率或热辐射波段的亮温；
步骤5：基于得到的像元地理坐标，通过设置仿射参数，对图像数据进行等经纬网格投影和几何校正；
步骤6：根据业务观测区域经纬度范围对监测图像数据进行区域裁剪，利用GDAL栅格数据处理技术生成TIF格式区域图像产品。
进一步地，步骤2的具体操作步骤包括：
步骤21：定期扫描FY-3D MERSI L1B轨道时序观测数据目录下所接收到的FY-3DMERSI L1B数据文件，进行标识并生成标识文件；
步骤22：判断是否有新的标识文件，若无则返回步骤21继续扫描；若有则进入步骤23进行完整检查；
步骤23：根据文件大小、数据文件数目以及标识文件时间属性，基于计划任务服务自动检查数据标识文件完整性，若文件完整则可进行数据读取，若文件不完整则返回步骤21。
进一步地，步骤4的具体步骤包括：
步骤41：计算大气表观反射率
步骤411：对每种波段像元数字量化值进行异常值和缺测值检测后，通过式(1)对合理的像元DN值进行辐射定标，得到辐射定标后的定标后辐射量化dn，即：
dn＝DN*slope+intercept (1)
其中，slope和intercept为辐射定标所需增益和偏移系数，可从数据集属性中获取；DN为通道像元辐射数字量化值，dn为辐射定标后的辐射值；
步骤412：根据订正系数计算该通道波段订正后的像元辐射R：
R＝Cal2*dn2+Cal1*dn+Cal0 (2)
其中，Cal2、Cal1、Cal0为订正系数；
步骤413：通过(3)式对大气顶层平均光谱辐射开展太阳高度角订正并计算光谱波段通道大气表观反射率vtoa：
其中，Em为大气顶层平均光谱辐射，d为日地天单位，θ为太阳天顶角；Lsensor为卫星传感器上入瞳辐亮度；
步骤414：通过(4)式计算Lsensor，并将式(4)代入式(3)得到式(5)，通过式(5)快速计算大气表观反射率：
Lsensor＝(R*Em)/π (4)
步骤42：计算亮温
步骤421：对热辐射波段异常值和缺测值进行检测后，对红外通道辐射L进行辐射定标，计算辐射定标后波段辐射亮度Lλ：
Lλ＝L*slope+intercept (6)
其中，L为通道未定标的辐射亮度量化值，slope为定标的增益，intercept为定标的偏移，Lλ为定标后的辐射亮度；
步骤422：根据普朗克函数转换公式(7)，计算热红外波段等效黑体亮温Tf：
其中，C1和C2为整合后的计算系数，h＝6.62606876e-34J.s为普朗克常量，c＝2.99792458e+8m/s为光速，k＝1.3806503e-23J/K玻尔兹曼常量，υ为该通道波数，Lλ为辐射定标后辐射亮度；
步骤423：通过(8)式对Tf进行订正，得到该波段观测到的亮温Tb：
Tb＝TBB_a*Tf+TBB_b (8)
其中，TBB_a和TBB_b为亮温订正系数，Tf热红外通道等效黑体亮温，Tb波段亮温，卫星遥感观测仪器观测的波段亮温。
进一步地，步骤5的具体步骤包括：
步骤51：根据读取的不同空间分辨率的像元地理坐标，获取不同分辨率图像的最大纬度Latmax、最小经度Lonmin和图像像元行高Res_Line、像元列宽Res_pixel；
步骤52：根据Latmax、Lonmin、Res_Line、Res_pixel设置卫星图像仿射参数，构建一个等经纬投影网格，通过式(9)-(10)构建原始图像像元地理坐标到投影网格行列号之间的关系：
Latvar＝Latmax+lines*Res_line (9)
Lonvar＝Lonmin+cols*Res_pixel (10)
其中，Latmax为卫星图像最大纬度，Latvar为卫星图像像元纬度，Res_line为投影网格行高，lines为投影网格行数；中Lonmin为卫星图像最小经度，Lonvar为卫星图像像元经度，Res_pixel为投影网格列宽，cols为投影网格列数；
步骤53：基于所述映射关系通过gdalWarp进行卫星图像几何校正和重投影。
进一步地，步骤6的具体步骤包括：
步骤61：获取FY-3D-MERSI待监测图像区域的经纬度坐标范围信息，对FY-3DMERSI图像数据进行裁剪；
步骤62：在裁剪后的区域内，基于GDAL栅格数据处理技术构建三维栅格数据集；
步骤64：将获取的像元太阳和卫星观测角度、计算出的光谱波段通道大气表观反射率和热红外波段像元亮温，根据图像数据行列维度大小，依次写入到三维栅格数据集中；
步骤65：输出1km、250m两种空间分辨率多通道栅格图像数据。
本方法与原有技术相比，具有以下有益效果：
本发明对轨道数据进行再处理，实现FY-3D MERSI L1B数据文件检测、元数据和地理数据解析、辐射定标、光谱反射率/亮温计算、重投影、几何校正和区域裁剪等主要模块耦合连接和触发式异步自动运行，提高FY-3D MERSI数据处理效率和自动化水平，自动生成和归档250m和1km两种空间分辨率多通道波段栅格遥感图像；通过实验表明经本发明方法处理加工后的多通道TIF图像的数据量整体小于处理前数据量，便于用户读取和调用；并且由于本发明采用网络文件传输协议NFS与FY-3D储存设备建立实时传输链接，实现了数据自动化检测和快速读取能力，耗时明显减少，从而在一定程度上提高数据再处理效率。
附图说明
图1为数据自动检测流程图；
图2为本发明所提自动再处理方法；
图3为自动再处理耗时时效变化趋势图；
图4为band13(0.709μm)反射率图像；其中图4(a)为原处理技术所得图像，图4(b)为本发明得到的图像；
图5为热红外波段band23(8.550μm)图像；其中图5(a)为原处理技术所得图像，图5(b)为本发明得到的图像；
图6为band1-band19可见光红外RSB波段的MAE和RMSE；其中图6(a)为MAE，图6(b)为RMSE；
图7为Band20-Band25中长波热红外辐射TEB波段亮温MAE和RMSE；其中图7(a)为MAE，图7(b)为RMSE；
图8为将原处理方法和本发明处理方法进行拟合的结果；其中图8(a)为波段band17线性回归拟合结果，8(b)为波段band25线性回归拟合结果；
图9为250m空间分辨率遥感图像；其中图9(a)为band3(0.650μm)，图9(b)为band24(10.8μm)。
具体实施方式
为了使本领域的普通技术人员能更好的理解本发明的技术方案，下面结合附图和实施案例对本发明的技术方案做进一步描述。
1、FY-3D MERSI L1B数据介绍
本发明要处理的数据为风云三省级直收站系统接收和预处理的FY-3D MERSI L1B数据，FY-3D的MERSI是第二代中分辨率光谱成像仪，光谱通道共有25个(见表1)，与上一代MERSI-I相比，增加了卷云、水汽和大气窗口陆表温度等共5个光谱通道。
风云三省级直收站接收和预处理的FY-3D MERSI L1B轨道过境观测数据主要有四个HDF5格式数据文件，分别是空间分辨率1km的定标观测数据集、空间分辨率1km的地理坐标数据、空间分辨率250m的定标观测数据集、空间分辨率250m的地理坐标数据(见表2)。从国家卫星气象中心数据服务网上可查询到数据说明，空间分辨率1km定标观测数据集主要包含4个子集，空间分辨率250m数据集含有6个子集以及不同分辨率地理坐标文件。表1为MERSI-I和MERSI-II光谱波段信息，表2为FY-3D MERSI L1B数据文件信息。
表1 MERSI-I和MERSI-II光谱信息
表2 FY-3D MERSI L1B文件信息
下面再从反射太阳波段(Reflective Solar Bands,RSB)反射率计算、热辐射波段(Thermal Emissive Bands,TEB)亮温计算、等经纬投影和栅格数据生成描述本发明所提的再处理方法，其基于Python开发平台和gdal处理模块从而实现数据在处理模块的集约管理和自动化运行。
2、反射率计算
FY-3D MERSI L1B级数据不同分辨率(1km，250m)文件数据存储类型均是UInt16像元光谱亮度DN值，不便直接应用开展遥感监测服务，需将各通道数据DN辐射定标后，计算相对应的反射率或辐射亮温，再应用于植被指数和陆表温度等基础观测因子反演模型中。
反射率是指卫星传感器接收到被遥感物体反射辐亮度，与被遥感物体吸收太阳辐亮度的比值，表征被遥感物体对太阳辐射的吸收和反射能力，本文提到的反射率为大气表观反射率，是卫星传感器接收到的包含大气和被遥感地面物体的光谱反射率总和。FY-3D的MERSI仪器的反射太阳波段(RSBs，Reflective Solar Bands)可见光近红外通道共有19个，其中250m空间分辨率通道4个，1km空间分辨率通道15个。计算光谱波段大气表观反射率，需消耗大量计算时间。每种波段像元数字量化(Digital Number，DN)值，经异常值和缺测值检测处理后，通过(1)式对合理DN值进行辐射定标，得到定标后辐射量化dn，如下：
dn＝DN*slope+intercept (1)
其中，slope和intercept为辐射定标所需增益和偏移系数，可从数据集属性中获取；DN为通道像元辐射数字量化值，dn为辐射定标后的辐射量化；
通过辐射定标后，再读取文件订正系数，计算该通道波段订正后的像元辐射R，如下：
R＝Cal2*dn2+Cal1*dn+Cal0 (2)
其中，Cal2、Cal1、Cal0为订正系数，可从L1数据文件中获取；
计算获得波段通道订正辐射R后，通过(3)式对大气顶层平均光谱辐射Em开展太阳高度角订正并计算光谱波段通道大气表观反射率ρtoa：
其中，Em为大气顶层平均光谱辐射，d为日地天单位，θ为太阳天顶角，这些变量可从L1B文件数据集属性中获取；Lsensor为卫星传感器上入瞳辐亮度，
Lsensor和订正辐射R、大气顶层平均光谱辐射Em存在如(4)式所示的关系，通过这种关系可计算Lsensor：
Lsensor＝(R*Em)/π (4)
将(4)式带入(3)式简化得到通道大气表观反射率计算公式(5)，将大气表观反射率ρtoa和Lsensor的关系转换成和通道定标订正辐射R之间的关系。通过式(5)可快速计算大气表观反射率ρtoa：
3、亮温计算
当一物体的光谱辐射亮度与某一黑体的光谱辐射亮度相等时，可利用该黑体的物理温度表征该物体的亮度温度(即亮温)，亮温虽有温度量纲和单位，但不具备温度的物理含义，是物体辐射亮度的代名词。FY-3D MERSI-II的热红外辐射光谱TEB波段共6个，其中具备2个250m空间分辨率热红外波段。该热红外波段数据存储为未定标辐射亮度整型量化值L。在完成热辐射波段异常值和缺测值检测处理后，对红外通道辐射L进行辐射定标，计算辐射定标后波段辐射亮度Lλ，求解过程如下：
Lλ＝L*slope+intercept (6)
其中，L为通道未定标的辐射亮度量化值，slope为定标的增益，intercept为定标的偏移，Lλ为定标后的辐射亮度；
此处的热红外波段辐射亮度数据的单位是mW/(m2*cm-1*sr),但大气校正或陆表温度反演算法模型中，辐射亮度单位有时为W/(m2*μm*sr)或μW/(cm2*sr*nm),根据辐射亮度具体应用场景的单位，必要时需要进行波数υ和波长λ的转换和单位变换，此单位变换采用现有技术中已公开方式，不再赘述。
根据普朗克函数转换公式(7)，计算热红外波段等效黑体亮温Tf：
其中，C1和C2为整合后的计算系数，h＝6.62606876e-34J.s为普朗克常量，c＝2.99792458e+8m/s为光速，k＝1.3806503e-23J/K玻尔兹曼常量，为该通道波数，Lλ为辐射定标后辐射亮度。
通过式(7)完成热红外通道等效黑体亮温Tf计算后，通过(8)式对Tf进行订正，从而得到该波段观测到的亮温Tb：
Tb＝TBB_a*Tf+TBB_b (8)
其中，TBB_a和TBB_b为通道亮温订正系数，Tf热红外通道等效黑体亮温，Tb波段亮温，卫星遥感观测仪器观测的波段亮温。
反射太阳波段反射率和热辐射波段亮温的计算主要代码为：
4、等经纬度重投影
由于FY-3D MERSI L1B数据图像没有等经纬投影信息，在光谱波段数据的反射率/辐射亮温计算处理后，需对图像数据进行重投影，基于原轨道数据经纬度信息，对图像进行等经纬投影。分别读取FY-3D L1B 1km和250m分辨率图像数据地理坐标，获取不同空间分辨率图像的最大纬度Latmax和最小经度Lonmin和图像像元行高Res_Line和像元列宽Res_pixel，根据这些参数设置卫星图像仿射参数，构建一个等经纬投影网格，如下。
Latvar＝Latmax+lines*Res_line (9)
Lonvar＝Lonmin+cols*Res_pixel (10)
其中，Latmax为卫星图像最大纬度，Latvar为卫星图像像元纬度，Res_line为投影网格行高，lines为投影网格行数；中Lonmin为卫星图像最小经度，Lonvar为卫星图像像元经度，Res_pixel为投影网格列宽，cols为投影网格列数；
通过(9)式和(10)式建立了原卫星图像地理坐标到新投影网格图像行列号之间映射关系，再利用gdalWarp即可完成卫星图像几何校正和重投影。
5、多通道栅格图像生成
本发明设置的FY-3D MERSI图像数据裁剪区范围为104.0E～112.0E和26.0N～34.0N，覆盖重庆长江上游地区及邻近省份部分区域，通过GDAL栅格数据处理技术，将太阳和卫星观测角度等辅助数据集、可见光近红外波段反射率、热红外波段辐射亮温数据集重新组合到一起，构建1km分辨率大小为29*800*800和250m分辨率大小为10*3200*3200的三维数据集，处理输出1km、250m两种空间分辨率多通道栅格图像数据。表3显示了处理输出图像文件信息，其中PRJ表示经过重投影处理，L2为数据处理等级，YYYYMMDDHHmm为卫星轨道数据接收日期，GLL为等经纬度投影，还有卫星、载荷和空间分辨率等信息也可以在文件中发现。
表3TIF图像产品信息
6、数据自动处理技术
关于FY-3D MERSI L1数据自动再处理技术，其首先要实现L1数据自动检测和数据完整性检查，该步骤主要基于Linux系统网络文件传输协议(NFS)和计划任务(Crond)服务，从而实现FY-3D MERSI L1B数据多频次自动化检测，所述自动化数据检测流程见附图1。当直收站有L1B数据更新时，部署直收站平台端的程序会定期检测，生成数据标识文件，自动检测模块扫描到新数据标识文件时，通过文件大小、文件数等逻辑条件检查数据文件完整性，自动读取数据文件数据信息。
FY3D MERSI L1数据自动再处理技术能够与省级直收站接收预处理平台紧密衔接，采用Python开发平台，结合HDF5文件和GDAL栅格数据处理模块，将数据检测和读取、辐射定标、反射率/亮温计算、等经纬投影、几何校正、区域裁剪和栅格图像产品生成等处理模块耦合连接和集约化管理，实现各模块间异步条件触发方式运行模式，提高FY-3D MERSIL1B级数据再处理的智能化水平。数据自动再处理技术流程见图2，输入端为FY-3D MERSIL1B过境轨道数据，输出为再处理加工后的FY-3D MERSI多通道栅格图像。
实施例
为了进一步验证本发明所提自动处理方法的有效性，选取了接收完整的2020年8月14日13时33分的FY-3D MERSI L1B数据，该时间点重庆中西部以晴空为主，云量较少，可见光波段和热红外波段图像数据较清晰，便于精确分析图像质量。对比分析本发明方法的图像质量。
1、处理前后数据量
经本发明方法处理前和处理后的平均数据大小情况见表4，处理加工后的多通道TIF图像的数据量整体要小于处理前数据量，1km空间分辨率数据处理后数据量减少了约117M，250m空间分辨率处理后数据量减少约206M。本发明为了便于技术人员快速读取，未能像已有的原处理方式那样，将数据转换为整型进行二进制存储，光谱波段数据以真实浮点类型存储，总体上处理后数据量没有原处理方式减少明显，但便于用户读取和调用。
表4处理前后数据大小(单位：MB)
Table 1 Data size before and after processing(unit:MB)
2、数据处理耗时
FY-3D MERSI L1数据自动再处理耗时时效变化趋势见图3，可以看出本发明处理耗时范围在79～128秒之间，平均运行耗时103.14秒，因服务器实际性能和网络波动因素影响，数据处理耗时呈现一定的变化波动。
为了深入分析再处理过程的耗时，对比了5个数据处理模块的耗时(表5)，原处理方式中卫星数据文件需人工手动检测和手动传输，耗时较大。本发明采用网络文件传输协议NFS与FY-3D储存设备建立实时传输链接，实现了数据自动化检测和快速读取能力，耗时明显减少。并且本发明其它处理节点耗时也基本优于原处理方式，总体上本发明耗时比原处理方式减少了约111秒，一定程度上提高数据再处理效率。
表5原处理方式和新处理技术主要处理模块平均耗时(单位：秒)
Table 7 the main models average runtime between old method and newtechnology(unit：seconds)
3、图像质量对比
以2020年8月14日13时33分的FY-3D MERSI L1B数据为例开展图像质量对比分析，因该时间点重庆中西部以晴天为主，上空云量较少，可见光通道和热红外通道的图像数据比较清晰，便于对比图像质量。
band13(0.709μm)是可见光波段，可用于植被指数NDVI计算，band23(8.550μm)是热红外辐射波段，可用于大气水汽反演，本发明选取这两个波段开展遥感图像质量直观比较。从原处理方式和本发明技术的band13(0.709μm)两幅图像的空间分布上看(如图4所示)，观测云的分布均集中在渝东部区域，重庆主城区的位置均清晰可见，主城区域地理位置指定准确，新处理方式在太阳高度角订正后，图像更亮度更高，整体图像更加清晰。原处理方式和本发明技术的热红外波段band23(8.550μm)图像对比来看(如图5所示)，受西南地区山地地形影响，发现图像通道辐射亮温变化范围明显。遥感图像在重庆中西部区域的亮温变化基本一致，在渝中西部地区红外通道辐射亮温基本地物亮温，且亮温值偏高。
通过计算光谱反射太阳光波段RSB和热辐射波段TEB的像元平均绝对误差MAE和均方根误差RMSE，进一步分析两种处理方式的图像数据误差。从band1-band19可见光红外RSB波段的MAE和RMSE(如图6所示)发现band5的MAE和RMSE最小，band8-band15光谱波段误差较大，总体上RSB波段反射率MAE不超过0.07，RMSE在0.08以内。从Band20-Band25中长波热红外辐射TEB波段亮温MAE和RMSE(如图7所示)，发现辐射亮温误差整体不大，其MAE不超过0.53，RMSE不超过0.8，其中band22的误差最小，在0.01以下。由于计算处理过程中数据精度的损失和插值算法差异，新处理技术与原处理方式的图像数据存在一定误差，但各光谱波段误差均控制在可接受范围内。这里选取了Band17和Band25两个波段，拟合原处理方式和新处理技术两个方法的结果，发现线性相关性，相关系数都在0.98左右(如图8所示)。
本发明也能够输出250m空间分辨率波段栅格遥感图像。以band3(0.650μm)红光和band24(10.8μm)热红外辐射波段图像数据为例，分析250m空间分辨率图像数据质量。图9为250m空间分辨率可见光的反射率图像和热红外波段辐射亮温图像。整体上看，可见光图像纹理更加精细，云量和云空间分布的精度更高，重庆主城区山脉和河流地理位置指向准确。通过两个波段的方差和标准差(表6)，band3标准差为0.256，band24标准差14.134，因未有效去除大气和云的影响，导致数据最大最小值差值较大，但红外波段辐射亮温标准差分布能够基本反映重庆及周边山地地形变化大的特征。
表6band3和band24数据信息
本说明书中未作详细描述的内容属于本领域专业技术人员公知的现有技术。尽管参照前述实施例对本发明专利进行了详细的说明，对于本领域的技术人员来说，其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
Claims (5)
1.一种FY-3D MERSI L1B数据自动再处理方法，其特征在于，包括以下步骤：
步骤1：通过风云三省级直收站实时接收并预处理为FY-3D MERSI L1B过境区域轨道观测数据；
步骤2：自动检测新接收的FY-3D MERSI L1B数据，并进行自动完整性检查和读取；
步骤3：解析检查后的观测数据属性，得到数据文件光谱波段像元的元数据信息；对元数据信息进行特征分析得到观测角、订正系数以及地理坐标信息；
步骤4：基于得到的观测角、订正系数计算光谱波段通道的大气表观反射率或热辐射波段的亮温；
步骤5：基于得到的像元地理坐标，通过设置仿射参数，对图像数据进行等经纬网格投影和几何校正；
步骤6：根据业务观测区域经纬度范围对监测图像数据进行区域裁剪，利用GDAL栅格数据处理技术生成TIF格式区域图像产品。
2.如权利要求1所述的一种FY-3D MERSI L1B数据自动再处理方法，其特征在于，步骤2的具体操作步骤包括：
步骤21：定期扫描FY-3D MERSI L1B轨道时序观测数据目录下所接收到的FY-3D MERSIL1B数据文件，进行标识并生成标识文件；
步骤22：判断是否有新的标识文件，若无则返回步骤21继续扫描；若有则进入步骤23进行完整检查；
步骤23：根据文件大小、数据文件数目以及标识文件时间属性，基于计划任务服务自动检查数据标识文件完整性，若文件完整则可进行数据读取，若文件不完整则返回步骤21。
3.如权利要求2所述的一种FY-3DMERSIL1B数据自动再处理方法，其特征在于，步骤4的具体步骤包括：
步骤41：计算大气表观反射率
步骤411：对每种波段像元数字量化值进行异常值和缺测值检测后，通过式(1)对合理的像元DN值进行辐射定标，得到辐射定标后的定标后辐射量化dn，即：
dn＝DN*slope+intercept (1)
其中，slope和intercept为辐射定标所需增益和偏移系数，可从数据集属性中获取；DN为通道像元辐射数字量化值，dn为辐射定标后的辐射值；
步骤412：根据订正系数计算该通道波段订正后的像元辐射R：
R＝Cal2*dn2+Cal1*dn+Cal0 (2)
其中，Cal2、Cal1、Cal0为订正系数；
步骤413：通过(3)式对大气顶层平均光谱辐射开展太阳高度角订正并计算光谱波段通道大气表观反射率ρtoa：
其中，Em为大气顶层平均光谱辐射，d为日地天单位，θ为太阳天顶角；Lsensor为卫星传感器上入瞳辐亮度；
步骤414：通过(4)式计算Lsensor，并将式(4)代入式(3)得到式(5)，通过式(5)快速计算大气表观反射率：
Lsensor＝(R*Em)/π (4)
步骤42：计算亮温
步骤421：对热辐射波段异常值和缺测值进行检测后，对红外通道辐射L进行辐射定标，计算辐射定标后波段辐射亮度Lλ：
Lλ＝L*slope+intercept (6)
其中，L为通道未定标的辐射亮度量化值，slope为定标的增益，intercept为定标的偏移，Lλ为定标后的辐射亮度；
步骤422：根据普朗克函数转换公式(7)，计算热红外波段等效黑体亮温Tf：
其中，C1和C2为整合后的计算系数，h＝6.62606876e-34J.s为普朗克常量，c＝2.99792458e+8m/s为光速，k＝1.3806503e-23J/K玻尔兹曼常量，υ为该通道波数，Lλ为辐射定标后辐射亮度；
步骤423：通过(8)式对Tf进行订正，得到该波段观测到的亮温Tb：
Tb＝TBB_a*Tf+TBB_b (8)
其中，TBB_a和TBB_b为亮温订正系数，Tf热红外通道等效黑体亮温，Tb波段亮温，卫星遥感观测仪器观测的波段亮温。
4.如权利要求3所述的一种FY-3D MERSI L1B数据自动再处理方法，其特征在于，步骤5的具体步骤包括：
步骤51：根据读取的不同空间分辨率的像元地理坐标，获取不同分辨率图像的最大纬度Latmax、最小经度Lonmin和图像像元行高Res_Line、像元列宽Res_pixel；
步骤52：根据Latmax、Lonmin、Res_Line、Res_pixel设置卫星图像仿射参数，构建一个等经纬投影网格，通过式(9)-(10)构建原始图像像元地理坐标到投影网格行列号之间的关系：
Latvar＝Latmax+lines*Res_line (9)
Lonvar＝Lonmin+cols*Res_pixel (10)
其中，Latmax为卫星图像最大纬度，Latvar为卫星图像像元纬度，Res_line为投影网格行高，lines为投影网格行数；中Lonmin为卫星图像最小经度，Lonvar为卫星图像像元经度，Res_pixel为投影网格列宽，cols为投影网格列数；
步骤53：基于所述映射关系通过gdalWarp进行卫星图像几何校正和重投影。
5.如权利要求4所述的一种FY-3D MERSI L1B数据自动再处理方法，其特征在于，步骤6的具体步骤包括：
步骤61：获取FY-3D-MERSI待监测图像区域的经纬度坐标范围信息，对FY-3DMERSI图像数据进行裁剪；
步骤62：在裁剪后的区域内，基于GDAL栅格数据处理技术构建三维栅格数据集；
步骤64：将获取的像元太阳和卫星观测角度、计算出的光谱波段通道大气表观反射率和热红外波段像元亮温，根据图像数据行列维度大小，依次写入到三维栅格数据集中；
步骤65：输出1km、250m两种空间分辨率多通道栅格图像数据。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310095470.XA CN116185616A (zh) | 2023-02-10 | 2023-02-10 | 一种fy-3d mersi l1b数据自动化再处理方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310095470.XA CN116185616A (zh) | 2023-02-10 | 2023-02-10 | 一种fy-3d mersi l1b数据自动化再处理方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN116185616A true CN116185616A (zh) | 2023-05-30 |
Family
ID=86443904
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310095470.XA Pending CN116185616A (zh) | 2023-02-10 | 2023-02-10 | 一种fy-3d mersi l1b数据自动化再处理方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116185616A (zh) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117115621A (zh) * | 2023-10-24 | 2023-11-24 | 中国海洋大学 | 基于改进U-Net网络的卫星云图预测方法 |
-
2023
- 2023-02-10 CN CN202310095470.XA patent/CN116185616A/zh active Pending
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117115621A (zh) * | 2023-10-24 | 2023-11-24 | 中国海洋大学 | 基于改进U-Net网络的卫星云图预测方法 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110750904B (zh) | 一种基于遥感数据的区域碳储量空间格局监测系统和方法 | |
| CN110348107B (zh) | 一种云下像元真实地表温度重建方法 | |
| CN114019579B (zh) | 高时空分辨率近地表空气温度重构方法、系统、设备 | |
| CN105678085A (zh) | 一种pm2.5浓度的估算方法及系统 | |
| CN113447137B (zh) | 一种面向无人机宽波段热像仪的地表温度反演方法 | |
| CN116519557B (zh) | 一种气溶胶光学厚度反演方法 | |
| CN104636608A (zh) | 一种modis卫星数据的直接同化方法 | |
| CN114564767A (zh) | 一种基于太阳-云-卫星观测几何的云下地表温度估算方法 | |
| CN115204691B (zh) | 基于机器学习和遥感技术的城市人为热排放量估算方法 | |
| CN117077437B (zh) | 基于多源卫星的极区海表净辐射模型构建和确定方法 | |
| CN116519913B (zh) | 基于星载和地基平台融合的gnss-r数据土壤水分监测方法 | |
| CN111191380B (zh) | 一种基于地基光谱仪测量数据的大气气溶胶光学厚度估算方法和装置 | |
| CN114970214A (zh) | 一种气溶胶光学厚度反演方法及装置 | |
| CN113408111A (zh) | 大气可降水量反演方法及系统、电子设备和存储介质 | |
| CN116185616A (zh) | 一种fy-3d mersi l1b数据自动化再处理方法 | |
| CN116822141A (zh) | 利用卫星微光遥感反演夜间大气气溶胶光学厚度的方法 | |
| Zhao et al. | A Nonlinear Split-Window Algorithm for Retrieving Land Surface Temperatures from Fengyun-4B Thermal Infrared Data | |
| Li et al. | Retrieving nighttime aerosol optical depth using combined measurements of satellite low light channels and ground-based integrating spheres | |
| Peng et al. | GOES-R Advanced Baseline Imager (ABI) algorithm theoretical basis document for surface albedo | |
| CN113740263A (zh) | 一种气溶胶光学厚度反演方法及大气颗粒物遥感反演方法 | |
| CN116558652A (zh) | 一种气象卫星热红外通道在轨辐射定标方法 | |
| CN115452167A (zh) | 基于不变像元的卫星遥感器交叉定标方法和装置 | |
| CN114296061B (zh) | 基于多元变量检测和不同辐射传输模型的交叉定标方法 | |
| CN115933007A (zh) | 基于多时空尺度卫星气象的天气监测方法及装置 | |
| KR20100028337A (ko) | 구름을 통한 통신해양기상위성의 복사 보정 방법 |
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 |