CN114114358B - 一种基于多源卫星数据融合的北极海冰厚度空间分辨率改进方法 - Google Patents
一种基于多源卫星数据融合的北极海冰厚度空间分辨率改进方法 Download PDFInfo
- Publication number
- CN114114358B CN114114358B CN202111400996.1A CN202111400996A CN114114358B CN 114114358 B CN114114358 B CN 114114358B CN 202111400996 A CN202111400996 A CN 202111400996A CN 114114358 B CN114114358 B CN 114114358B
- Authority
- CN
- China
- Prior art keywords
- data
- product
- sea ice
- products
- observation
- 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
Links
- 230000004927 fusion Effects 0.000 title claims abstract description 85
- 238000000034 method Methods 0.000 title claims abstract description 45
- 230000006872 improvement Effects 0.000 title claims description 5
- 238000005259 measurement Methods 0.000 claims abstract description 41
- 238000001228 spectrum Methods 0.000 claims abstract description 33
- 230000000694 effects Effects 0.000 claims abstract description 19
- 235000021581 juice product Nutrition 0.000 claims abstract description 4
- 238000007781 pre-processing Methods 0.000 claims abstract description 4
- 238000004458 analytical method Methods 0.000 claims description 35
- 238000005070 sampling Methods 0.000 claims description 24
- 238000004364 calculation method Methods 0.000 claims description 8
- 238000012935 Averaging Methods 0.000 claims description 6
- 230000003595 spectral effect Effects 0.000 claims description 6
- 238000001514 detection method Methods 0.000 claims description 4
- 238000011156 evaluation Methods 0.000 claims description 3
- 230000003287 optical effect Effects 0.000 claims description 3
- 230000008569 process Effects 0.000 claims description 3
- 238000012216 screening Methods 0.000 claims description 3
- 239000002689 soil Substances 0.000 claims description 3
- 238000012545 processing Methods 0.000 claims 1
- 230000002123 temporal effect Effects 0.000 claims 1
- 230000008901 benefit Effects 0.000 abstract description 3
- 238000011161 development Methods 0.000 description 4
- 238000010586 diagram Methods 0.000 description 3
- 230000006978 adaptation Effects 0.000 description 2
- 238000012937 correction Methods 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 239000000700 radioactive tracer Substances 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- ZYJDYPZRRJTBMY-UHFFFAOYSA-N Troupin Natural products COc1c(O)cc2C(=COC(=O)c2c1OC)C ZYJDYPZRRJTBMY-UHFFFAOYSA-N 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000010168 coupling process Methods 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000005553 drilling Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 238000013213 extrapolation Methods 0.000 description 1
- 230000004907 flux Effects 0.000 description 1
- 239000013505 freshwater Substances 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 230000007774 longterm Effects 0.000 description 1
- 238000007500 overflow downdraw method Methods 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 230000004044 response Effects 0.000 description 1
- 239000013535 sea water Substances 0.000 description 1
- 238000010183 spectrum analysis Methods 0.000 description 1
- 230000007480 spreading Effects 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
- G01S19/38—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system
- G01S19/39—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system the satellite radio beacon positioning system transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
- G01S19/42—Determining position
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
- G01B7/00—Measuring arrangements characterised by the use of electric or magnetic techniques
- G01B7/02—Measuring arrangements characterised by the use of electric or magnetic techniques for measuring length, width or thickness
- G01B7/06—Measuring arrangements characterised by the use of electric or magnetic techniques for measuring length, width or thickness for measuring thickness
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/25—Fusion techniques
Landscapes
- Engineering & Computer Science (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Data Mining & Analysis (AREA)
- Artificial Intelligence (AREA)
- Life Sciences & Earth Sciences (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Evolutionary Biology (AREA)
- Evolutionary Computation (AREA)
- General Engineering & Computer Science (AREA)
- Computer Networks & Wireless Communication (AREA)
- Radar Systems Or Details Thereof (AREA)
Abstract
本发明公开了一种基于多源卫星数据融合的北极海冰厚度空间分辨率改进方法，首先获取海冰厚度产品多源数据和海冰厚度实测资料，包括PIOMAS、CS2SMOS、APP‑x和CPOM四种海冰厚度产品和ULS、IMB、IceBridge三种海冰厚度实测资料；然后对获取的海冰资料进行预处理，将月资料展开为日数据，对冗余的实测资料进行降采样后与海冰产品进行时空匹配；基于ULS的长时间序列实测资料，选取精度最高的日产品与空间分辨率高的月产品作为待融合资料，代入Diva模型进行数据融合；设计多组参数方案，选取效果最好的参数确定融合模型的参数；最后通过波数谱对比各产品的空间分辨率，检验融合效果；本发明提供的方法通过Diva插值融合了不同产品的优点，可以显著提高产品的时空分辨率。
Description
技术领域
本发明涉及遥感地学技术领域，主要涉及一种基于多源卫星数据融合的北极海冰厚度空间分辨率改进方法。
背景技术
北极作为北半球的冷空气源地和热量的“汇”，对北半球大气、海洋环流和能量体系的平衡都有重要意义。海冰厚度在海冰气耦合过程中起到了很重要的作用。它的变化更能影响海气通量，淡水收支，能量平衡等。除了气候方面，海冰厚度也是实际作业的一个重要因素，如船舶冰区航行和冰上钻井等。准确的厚度观测不仅影响海冰本身的建模预测，还对气候长期的全球响应以及航道利用、资源开发等社会经济领域影响显著。
至今为止，虽然海冰厚度观测手段很多，但结果受限于时间跨度、空间分辨率和观测精度，观测记录是不充分的，且可靠性有限。当前，实测数据获取代价高、观测时空连续性差，海冰厚度资料主要包括卫星遥感和数值模式资料，但是这些产品在精度和时空分辨率上都或多或少存在缺陷。单纯的偏差修正都过度依赖于实测数据，在实测资料不足的情况下，训练出的模型的适用性受实测数据的取样时间和取样地点所限制，只能证明其在取样多的北极西半球和冬半年修正效果较佳，由于实测值在北极东半球和夏半年的缺失，无法验证模型在取样点缺少的区域和时间的修正效果。海冰厚度产品种类很多，多源数据融合分析可以优势互补，减少不确定性、提高时空分辨率，其可以将多种精度的厚度资料融合提高整体的产品精度，实测充分时可以融合精度最高的实测产品，在实测数据不足时，也可以将不同精度的产品融合。采用第二种方式时，面对实测数据不足的情况效果也不会太片面，因为实测数据仅用于评估融合方案。
因为海冰厚度日产品数据的空间分辨率有限，前人的工作多集中于提高微波辐射计产品精度，或者提高卫星高度计产品的空间覆盖率，目前已经广泛使用的CS2SMOS产品便是微波辐射计SMOS和卫星高度计CS2融合后的产品，该产品有效地提高了CS2的空间覆盖率和薄冰区的精度，但是其空间分辨率仍然是25km左右，无法满足航行保障或者中小尺度分析的需求。目前并没有相关工作是提高日产品的空间分辨率。
发明内容
发明目的：针对上述背景技术中存在的问题，本发明提供了一种基于多源卫星数据融合的北极海冰厚度空间分辨率改进方法，通过基于变分分析的数据插值方法融合高精度的低空间分辨率海冰厚度日(周)产品和低精度高空间分辨率的海冰厚度月产品能有效提升海冰厚度产品的精度和空间分辨率，得到高精度、高时空分辨率的海冰厚度产品，填补了相关产品的空缺，为进行气候分析和北极航道开发、北极资源开发等经济活动领域提供保障。
技术方案：为实现上述目的，本发明采用的技术方案为：
一种基于多源卫星数据融合的北极海冰厚度空间分辨率改进方法，包括以下步骤：
步骤S1、数据获取；准备海冰厚度数据并对数据进行提取和筛选；获取的数据包括海冰厚度产品多源数据和海冰厚度实测资料；
所述海冰厚度产品多源数据包括海冰厚度日/周产品和海冰厚度月产品；所述海冰厚度周产品包括CryoSat-2和土壤水分海洋盐度卫星SMOS的融合产品CS2SMOS、日产品包括依据能量预算模型反演光学卫星AVHRR得到的产品APP-x和同化海表温度和海冰密集度的冰-海模式产品PIOMAS；所述海冰厚度月产品包括CPOM中心得到的CryoSat-2的月平均资料；所述海冰厚度产品均为格点数据；
所述海冰厚度实测资料包括飞机航拍观测产品IceBridge、物质质量平衡浮标观测产品IMB和水下定点浮标向上声呐探测产品ULS；
采用PIOMAS、CS2SMOS、APP-x和CPOM四种海冰厚度产品进行后续融合，寻找最佳的融合方案，并以ULS、IMB、IceBridge作为实测值进行融合效果评估；
步骤S2、对步骤S1获取的海冰厚度产品数据进行预处理；采用时间平均或扩展的方式进行时间匹配；采用小区域平均法和降采样的方式将格点数据与实测值进行空间匹配；
步骤S3、基于ULS的长时间序列实测资料，选取精度最高的日产品与空间分辨率高的月产品作为待融合资料，代入Diva模型进行数据融合；设计多组参数方案，选取融合结果的偏差和均方根偏差最小、相关系数最大时的参数确定为最佳融合模型的参数；
步骤S4、融合效果检验；将得到的最佳数据融合模型扩展至北极全域，通过波数谱对比各产品的空间分辨率，检验融合效果。
进一步地，所述步骤S2中分别采用小区域平均法和降采样的方式将格点数据与实测值进行时间匹配和空间匹配，具体包括：
所述步骤S2中对格点数据与实测值进行时间匹配和空间匹配，具体包括：
步骤S2.1、对各海冰厚度产品进行时间匹配；
针对CS2SMOS，选取观测当天、观测前3天及观测后3天的所有观测数值，取均值作为当天的观测值；将CPOM由月平均化的格点数据展开为日数据；选取ULS的日平均值作为该点的海冰厚度；针对IceBridge和IMB，选取每个采样点采样时间所在的日期为该点的采样日期；
步骤S2.2、基于小区域平均和降采样方式进行空间匹配；针对ULS观测值，以取样点为圆心，针对CS2SMOS、APP-x和PIOMAS产品，以25km为半径，针对CPOM产品，以5km为半径，圆内的格点观测值的平均值作为对应ULS观测处的产品观测值；针对CS2SMOS、PIOMAS和APP-x产品，以取样点为圆心，将半径25km内的格点观测值的平均值作为对应观测处的产品观测值；针对CPOM产品，以取样点为圆心，将半径5km内的格点观测值的平均值作为对应观测处的产品观测值；针对IMB和IceBridge，实测值空间分布密集观测冗余，为减少计算量，先确定IMB或IceBridge实测值的空间覆盖区域，取覆盖区域外扩1°为待计算区域，以待计算区域的所有产品观测格点为中心，对区域内的实测值以反比例加权的方式与海冰产品进行空间匹配，通过对冗余实测值加权平均，实现对实测数据的降采样。
进一步地，所述步骤S2.2中空间匹配具体数据处理方法包括：
步骤S2.2.1、针对IceBridge或IMB资料，以四类海冰厚度产品在北极的空间格点(Xpi,Ypi)k为基准，计算与当日每个实测点(Xoj,Yoj)k的球面距离Dijk如下：
其中，i＝1,…,m，j＝1,…,n，k＝1,…,t；m为海冰厚度产品在北极空间覆盖的点数，n为IceBridge或IMB在第k天探测到的空间点数，t为有实测数据的日期数目；
当球面距离Dijk满足以下条件：
Dijk≤L
则认为(Xoj,Yoj)k处于第k日以(Xpi,Ypi)k为圆心，以L为半径的圆内；其中当产品选择CS2SMOS、APP-x或PIOMAS时，L＝25km，当产品选择CPOM时，L＝5km；根据不同海冰厚度产品，选取不同L，将满足上述条件的Dijk对应的海冰厚度实测值依据Dijk的大小进行反距离加权平均，将均值结果作为(Xpi,Ypi)k产品格点处的实测值，实现降采样过程；
步骤S2.2.2、针对ULS，以实测点(Xoj,Yoj)k，j＝1,2,3为圆心，依据步骤S2.2.1所述球面距离计算公式，寻找半径为25km内包含的s个空间点(Xpi,Ypi)k，i＝1,…,s；将对应的s个观测值求均值，作为第k天在(Xoj,Yoj)k处的四类产品的观测值。
进一步地，所述步骤S2.2.1中各类海冰厚度产品在北极覆盖的空间点数m选取分别如下：
对于PIOMAS产品，m＝360×120；对于APP-x产品，m＝361×361，对于CS2SMOS产品，m＝432×432，对于CPOM产品，m＝432×432；
根据不同算力情况，基于(Xoj,Yoj)k的覆盖范围，可以减少经纬度范围，m取值随k的取值变化。
进一步地，所述步骤S3中，将四类产品分布与连续定点观测的ULS实测值计算差值均值和差值方差并对比，选取最高精度的海冰厚度日产品作为观测场，选取高空间分辨率的CPOM产品作为背景场，进行Diva插值，得到分析场，进而获取高精度插值场，获取数据融合模型；具体地，
基于Diva插值得到的分析场的成本函数J包括分析场与插入值之间的距离和分析场的规律，具体表示如下：
其中代表分析场与插入值之间的距离，即观测约束，N代表插值点(xi，yi)的个数，di为观测值的异常值，/>为分析场中对应点的取值，μi为待确定的各项权重参数；
代表分析场的规律，及光滑约束，用于计算分析场在Ω域上的空间变异性；其中▽为水平梯度，“:”代表双点乘运算，α0和α1为待确定的参数；
将所述三个待确定的参数简化为两个有物理意义的参数，包括相关长度len和信噪比SNR；len表示插值点对周围分析场产生影响的距离，单位为千米；SNR为背景误差协方差和观测误差协方差之比，SNR愈大，观测场相对背景场更加准确，需要对比待融合的资料，由外界输入。
进一步地，步骤S4中通过波数谱对比各产品的空间分辨率，具体包括：
使用经向谱，沿区域内每条经度线上的序列数据求取波数谱，并进行平均；首先将序列数据去除趋势性，然后运用Black-Harris窗，防止谱能量泄露，通过Burg法估计融合模型的波数谱。
有益效果：
本发明提供的基于多源卫星数据融合的北极海冰厚度空间分辨率改进方法，主要目的在于借助变分分析的方法融合多源海冰厚度产品的优点，得到高分辨率的厚度日产品。具体做法是面对现有的、数据质量层次不齐的多源海冰厚度数据，通过产品质量评估挑选出合适的观测场和背景场，将传统的Diva插值算法迁移至数据融合，整合多源数据，调整信噪比，得到参数不定的数据融合模型，再对比实测数据确定最佳参数和最佳融合模型，将该精度高模型推广至整个北极，选取浮冰区的海冰厚度应用波数谱进行验证，证明最终得到的产品的有效空间分辨率可以达到5km，验证北极海冰厚度日产品的空间分辨率也将得到显著的提升。
本发明使用了高空间分辨率的CPOM产品，避开了融合缺乏代表性的海冰厚度实测数据，而是融合多源数据提高产品精度。相对于现有多源数据融合的方法来说，本方法创新性的融合了高空间分辨率低时间分辨率的CPOM产品和低空间分辨率高时间分辨率的CS2SMOS产品。二者优势互补且使用的该采用的Diva算法较传统的最优插值算法计算效率高且效果类似，结果显著提高了产品的时空分辨率。本发明相对于其他方法简单高效地得到高空间分辨率的海冰厚度日产品，为提高产品时空分辨率提供了简单快速的思路，可以满足未来航行保障、资源开发或者中小尺度分析的需求。
附图说明
图1是本发明提供的基于多源卫星数据融合的北极海冰厚度空间分辨率改进方法流程图；
图2是本发明实施例中2010.11.18-2018.3.31海冰厚度产品观测值平均场空间分布示意图；
图3是四种海冰厚度产品及ULS A在210.42°E，75°N位置的海冰厚度观测值时间序列；
图4是2021年3月21日海冰厚度实测资料的空间分布图；
图5是2013年3月21日采用Diva算法，在不同SNR取值下CS2SMOS与CPOM融合产品的海冰厚度示意图；
图6是不同SNR取值时CS2SMOS与CPOM融合产品的误差曲线图；
图7是2013年3月21日CPOM产品、CS2SMOS产品以及当SNR为25时数据融合模型的海冰厚度空间分布示意图；
图8是CPOM产品、CS2SMOS产品以及CS2SMOS与CPOM融合产品在[84°N-87°N，353.5°E-357.5°E]海区内的经向波数谱。
具体实施方式
下面结合附图对本发明作更进一步的说明。显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
本发明提供的基于多源卫星数据融合的北极海冰厚度空间分辨率改进方法，具体如图1所示，包括以下步骤：
步骤S1、准备海冰厚度数据并对数据进行提取和筛选；获取的数据包括海冰厚度产品多源数据和海冰厚度实测资料；
海冰厚度产品多源数据包括海冰厚度日产品和海冰厚度月产品；所述海冰厚度日产品包括CryoSat-2和土壤水分海洋盐度卫星SMOS的融合产品CS2SMOS、依据能量预算模型反演光学卫星AVHRR得到的产品APP-x和同化海表温度和海冰密集度的冰-海模式产品PIOMAS；所述海冰厚度月产品包括CPOM中心得到的CryoSat-2的月平均资料；所述海冰厚度产品均为格点数据；
海冰厚度实测资料包括飞机航拍观测产品IceBridge、物质质量平衡浮标观测产品IMB和水下定点浮标向上声呐探测产品ULS；
本实施例中收集到的资料来自官网公开发布，资料为2010.11.18-2018.3.31的北极圈以北的所有的海冰厚度及其时间、经纬度资料。
步骤S2、对步骤S1获取的海冰厚度产品数据进行预处理；分别采用小区域平均法和降采样的方式将格点数据与实测值进行时间匹配和空间匹配；具体地，
步骤S2.1、对各海冰厚度产品进行时间匹配；针对CS2SMOS，选取观测当天、观测前3天及观测后3天的所有观测数值，取均值作为当天的观测值；将CPOM由月平均化的格点数据展开为日数据；选取ULS的日平均值作为该点的海冰厚度；针对IceBridge和IMB，选取每个采样点采样时间所在的日期为该点的采样日期；
步骤S2.2、基于小区域平均和降采样方式进行空间匹配；针对ULS观测值，以取样点为圆心，针对CS2SMOS、APP-x和PIOMAS产品，以25km为半径，针对CPOM产品，以5km为半径，圆内的格点观测值的平均值作为对应ULS观测处的产品观测值；针对CS2SMOS、PIOMAS和APP-x产品，以取样点为圆心，将半径25km内的格点观测值的平均值作为对应观测处的产品观测值；针对CPOM产品，以取样点为圆心，将半径5km内的格点观测值的平均值作为对应观测处的产品观测值；针对IMB和IceBridge，实测值空间分布密集观测冗余，为减少计算量，先确定IMB或IceBridge实测值的空间覆盖区域，取覆盖区域外扩1°为待计算区域，以待计算区域的所有产品观测格点为中心，对区域内的实测值以反比例加权的方式与海冰产品进行空间匹配，通过对冗余实测值加权平均，实现对实测数据的降采样。
针对IceBridge或IMB资料，以四类海冰厚度产品在北极的空间点(Xpi,Ypi)k为基准，计算与当日每个实测点(Xoj,Yoj)k的球面距离Dijk如下：
其中，i＝1,…,m，j＝1,…,n，k＝1,…,t；m为海冰厚度产品在北极空间覆盖的点数，n为IceBridge或IMB在第k天探测到的空间点数，t为有实测数据的日期数目；对于PIOMAS产品，m＝360×120；对于APP-x产品，m＝361×361，对于CS2SMOS产品，m＝432×432，对于CPOM产品，m＝432×432.
根据不同算力情况，基于(Xoj,Yoj)k的覆盖范围，可以减少经纬度范围，m取值随k的取值变化。
当球面距离Dijk满足以下条件：
Dijk≤L
则认为(Xoj,Yoj)k处于第k日以(Xpi,Ypi)k为圆心，以L为半径的圆内；其中当产品选择CS2SMOS、APP-x或PIOMAS时，L＝25km，当产品选择CPOM时，L＝5km；根据不同海冰厚度产品，选取不同L，将满足上述条件的Dijk对应的海冰厚度实测值依据Dijk的大小进行反距离加权平均，将均值结果作为(Xpi,Ypi)k产品格点处的实测值，实现降采样过程；
针对ULS，以实测点(Xoj,Yoj)k，j＝1,2,3为圆心，依据步骤S2.2.1所述球面距离计算公式，寻找半径为25km内包含的s个空间点(Xpi,Ypi)k，i＝1,…,s；将对应的s个观测值求均值，作为第k天在(Xoj,Yoj)k处的四类产品的观测值。
步骤S3、基于ULS的长时间序列实测资料，选取精度最高的日产品与空间分辨率高的月产品作为待融合资料，代入Diva模型进行数据融合；设计多组参数方案，选取融合结果的偏差和均方根偏差最小、相关系数最大时的参数确定为最佳融合模型的参数；。
本实施例中选取高精度海冰厚度日产品作为观测场，选取高空间分辨率的CPOM产品作为背景场，进行Diva插值，得到分析场，进而获取高精度插值场，获取数据融合模型。Diva是Troupin等2012年公开发布的基于变分分析的数据插值方法，同时发布的是相应的MATLAB工具包，在海洋温度、盐度等数据的插值上得到了成功应用，其初步应用时将信度高的离散实测资料作为观测场，将信度低的网格场作为背景场，得到的分析场可以减少外插的误差，得到精度较高的插值场。本发明将低分辨率但是精度高的日产品作为观测场，高分辨率但是精度低的CPOM作为背景场，此处的插值可以理解为两种资料的融合，可以较为合理的提高产品的空间分辨率。
基于Diva插值得到的分析场的成本函数J包括分析场与插入值之间的距离和分析场的规律，具体表示如下：
其中代表分析场与插入值之间的距离，即观测约束，N代表插值点(xi，yi)的个数，di为观测值的异常值，一般为原始值减去气候态，在没有气候态参考值时可以减去均值代替气候态。/>为分析场中对应点的取值，μi为待确定的各项权重参数；
代表分析场的规律，及光滑约束，用于计算分析场在Ω域上的空间变异性；其中▽为水平梯度，“:”代表双点乘运算，α0和α1为待确定的参数；
如果只考虑观测约束，结果将是纯插值，并不适用于包含很多噪音的气象水文要素，因为观测值除了观测误差，观测值无法完全代表气候态这也带来了代表性误差，加上光滑约束，可以提高插值水平。
将所述三个待确定的参数简化为两个有物理意义的参数，包括相关长度len和信噪比SNR；len表示插值点对周围分析场产生影响的距离，单位为千米；SNR需要对比待融合的资料，由外界输入。关于参数简化操作可以参考Troupin C,Barth A,Sirjacobs D,etal.Generation of analysis and consistent error fields using the DataInterpolating Variational Analysis(DIVA).Ocean Modelling.2012；52-53:90-101.
进行Diva插值时需选取合适的产品作为背景场和观测场，因为融合精度较高的依赖于观测值，需先依据实测资料对产品质量进行评估，选取合适的融合产品，参数的确定需要依据实测资料手动调参，选取精度最好的时的参数确定最佳融合方案。
步骤S4、将得到的最佳数据融合模型扩展至北极全域，通过波数谱对比各产品的空间分辨率，检验融合效果。
波数谱是以波数为横轴绘制的功率谱，最简单直观的波数谱通过傅立叶变换得到的波数和谱能量(即傅立叶空间内振幅的平方)进行表示。本发明使用的经向谱是沿区域内每条经度线上的序列数据求取波数谱，然后进行平均，参考前人方式，首先将序列数据去趋势，然后运用Black-Harris窗防止谱能量泄露，采用Burg法(MATLAB自带函数)估计其波数谱。波数谱反映了从大尺度到中小尺度的能量衰减规律，前人研究显示由准地转湍流理论计算的理想的波数谱的斜率大约为-2，沿等密度线运动的被动示踪物的波数谱斜率也是约等于-2，海冰并不属于严格的被动示踪物，但是海冰厚度与海水温度、盐度等都关系密切，二者的波谱斜率处于-3和-2之间，所以此处本发明依据此认为海冰厚度合理的波谱斜率也是处于-3和-2之间。在谱线斜率合理的情况下，谱能量高则证明含有的中小尺度能量更为丰富，即数据产品含有更多的中小尺度结构且这些结构是与大尺度合理衔接的，通常说明数据的质量更高。
下面结合实验给出详细实施方式，进一步验证数据融合模型对分辨率改进的效果。
以2010.11.18-2018.3.31的北极海冰厚度为出发点，选择实测值较多的2021年3月21日的多源海冰厚度产品融合为例，结合附图1所示步骤，进行融合试验，研究区域包括整个北极，CS2SMOS等日产品空间分辨率为25km，CPOM空间分辨率为5km。
收集PIOMAS、CS2SMOS、APP-x、CPOM四种卫星和模式产品，以及ULS、IMB、IceBridge等实测值以备后期融合效果评估。其中除CPOM产品为月数据外，其余产品皆为日(周)产品。但是CPOM是空间分辨率最高的产品，产品分辨率可达5km(在极地≈0.03°)。2010.11.18-2018.3.31的多年平均场如图2所示，可见四种产品形态差异很大。
将实测数据与厚度产品进行时空匹配。首先为寻找精度最高的一类产品，计算四种产品在ULS A、ULS B、ULSD的观测值并与实测值对比，依据上述球面距离计算公式寻找以ULS为圆心、半径固定的圆内所包含的空间点平均值，为保证区域圆尽可能小且内有产品观测值，此处CPOM的区域圆半径定义为5km，其他三类产品定为25km。以ULS A为例，四者在浮标处的时间序列和浮标的实测值如图3所示。在(210.42°E，75°N)处，实测海冰厚度(SIT)平均值为1.02m，产品观测值皆偏高，相差最少的产品是1.14m的CS2SMOS，其次是1.31m的CPOM。实测方差较大，四类产品的观测值波动较小，尤其是CPOM，时间曲线呈现方波状，无法达到实测值的时间分辨率。对比产品差值可以得到差值平均值最小的是CS2SMOS和CPOM，差值方差最小的是PIOMAS和CS2SMOS，差值方差最大的是CPOM。由ULS A处的观测可知，质量最好的产品是CS2SMOS。
将CS2SMOS(半径为25km)作为观测场，CPOM(半径为5km)作为背景场，通过Diva插值可得到空间分辨率为5km的产品。为保证融合产品的精度，需确定SNR、len两个参数的取值。因为选取的观测场的空间分辨率为25km，所以相关长度len定义为25km，SNR无法确定，分别取0.1,0.5,1,5,10,25,50七个值进行试验。为对比不同参数的产品融合效果，需要尽量多的实测值进行评估，具体实验结果如图4所示，＝IceBridge在2013年3月21日进行了长距离的观测，IMB在该天也有观测资料，故而选取该天的融合产品进行对比分析。
如图5所示，实测值曲线波动明显，空间分辨率最高，CPOM次之，CS2SMOS波动最小，空间分辨率最低。两种产品对1m以下的薄冰探测误差较大，正如1-80的取样点，二者整体呈现薄冰厚估的问题，CS2SMOS偶尔有几个低值点观测较为准确，可能是融合了SMOS的缘故。
Diva插值融合产品的取值并没有限制于两类原始产品的取值之间，融合的SIT误差曲线在不同SNR取值时趋势类似，如图6所示，可见融合效果受CS2SMOS原始产品质量限制较为明显。
下表1计算了原始CS2SMOS、CPOM产品和不同SNR取值时融合产品的均值、方差，并参照实测值计算偏差、均方根误差和相关系数。以发现当日CS2SMOS质量明显强于CPOM，使用Diva融合后的产品偏差、均方根误差均减小，R增加，且SNR越大，产品质量越好，在SNR取25或50时，效果提升到相似的程度。此处我们取SNR为25，len为25km的模型作为最佳融合模型，得到融合前后的产品对比如图7所示。由表1可知，融合后的产品精度得到了很大的提升，视觉上融合产品空间细微结构更明显，但是该分辨率是否有效，还需要进一步采用波数谱的方式进行衡量。
表1CS2SMOS和CPOM多源卫星数据融合效果对比结果
本发明对选定SIT不缺测、经纬度范围为[84°N-87°N，353.5°E-357.5°E]的海区作为研究区域计算CPOM和CS2SMOS及其融合产品的经向谱，结果如图8所示。波数谱反映了从大尺度到中小尺度的能量衰减规律，SIT合理的波谱斜率也是处于-3和-2之间。融合产品和CPOM的波数谱在150km至10km完全重合，且在50km至10km时其波谱斜率处于-2和-3之间，可见该产品可以识别SIT在10km上的中尺度变化，根据nyquist折叠频率，融合产品和CPOM的空间分辨率皆可达到5km。而CS2SMOS的波谱曲线自50km后突然上折，说明其无法识别25km以下的中尺度变化。证明了融合产品5km的空间分辨率是有效的。至于CPOM和融合产品在50km的突然凸起可能是因为海冰厚度不是连续变量，其观测值在此处出现了跃值，致使波谱能量也突然变化，但是其小尺度的观测还是相对连续的。
依据定点浮标实测值对CPOM月产品(5km)、PIOMAS日产品(≈25km)、CS2SMOS日产品(≈25km)、APP-x(≈25km)的日产品在SIT日尺度上的精度进行评估，选出精度最好的CS2SMOS产品和空间分辨率最高、时间分辨率稍差的CPOM作为待融合的产品带入Diva和加权平均模型，以2021年3月21日的实测值为参考选取合适的参数值，最终得到相关长度为25km、信噪比为25的Diva模型得到的产品精度最好，偏差减少至0.08m，均方根偏差也减少至0.46，结果比任一融合前的产品质量都好。通过波数谱分析，可以得见融合后的产品有效空间分辨率可以达到5km，证明了此次多源卫星数据融合的科学性。
将该模型继续推广至所有时间段，本发明成功制造出了精度和空间分辨率皆提高的5km逐日海冰厚度产品，这是本方法的一大创新点，通过多源卫星数据融合成功弥补了卫星观测在海冰厚度产品上时间分辨率、空间分辨率和空间覆盖率不能同时兼顾的缺点，且有效提高了产品的精度。
以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。
Claims (5)
1.一种基于多源卫星数据融合的北极海冰厚度空间分辨率改进方法，其特征在于，包括以下步骤：
步骤S1、数据获取；准备海冰厚度数据并对数据进行提取和筛选；获取的数据包括海冰厚度产品多源数据和海冰厚度实测资料；
所述海冰厚度产品多源数据包括海冰厚度日/周产品和海冰厚度月产品；所述海冰厚度周产品包括CryoSat-2和土壤水分海洋盐度卫星SMOS的融合产品CS2SMOS、日产品包括依据能量预算模型反演光学卫星AVHRR得到的产品APP-x和同化海表温度和海冰密集度的冰-海模式产品PIOMAS；所述海冰厚度月产品包括CPOM中心得到的CryoSat-2的月平均资料；所述海冰厚度产品均为格点数据；
所述海冰厚度实测资料包括飞机航拍观测产品IceBridge、物质质量平衡浮标观测产品IMB和水下定点浮标向上声呐探测产品ULS；
采用PIOMAS、CS2SMOS、APP-x和CPOM四种海冰厚度产品进行后续融合，寻找最佳的融合方案，并以ULS、IMB、IceBridge作为实测值进行融合效果评估；
步骤S2、对步骤S1获取的海冰厚度产品数据进行预处理；采用时间平均或扩展的方式进行时间匹配；采用小区域平均法和降采样的方式将格点数据与实测值进行空间匹配；
步骤S3、基于ULS的长时间序列实测资料，选取精度最高的日产品与空间分辨率高的月产品作为待融合资料，代入Diva模型进行数据融合；设计多组参数方案，选取融合结果的偏差和均方根偏差最小、相关系数最大时的参数确定为最佳融合模型的参数；
步骤S4、融合效果检验；将得到的最佳数据融合模型扩展至北极全域，通过波数谱对比各产品的空间分辨率，检验融合效果；
所述步骤S3中，将四类产品分布与连续定点观测的ULS实测值计算差值均值和差值方差并对比，选取最高精度的海冰厚度日产品作为观测场，选取高空间分辨率的CPOM产品作为背景场，进行Diva插值，得到分析场，进而获取高精度插值场，获取数据融合模型；具体地，
基于Diva插值得到的分析场的成本函数J包括分析场与插入值之间的距离和分析场的规律，具体表示如下：
其中代表分析场与插入值之间的距离，即观测约束，N代表插值点(xi，yi)的个数，di为观测值的异常值，/>为分析场中对应点的取值，μi为待确定的各项权重参数；
代表分析场的规律，即光滑约束，用于计算分析场在Ω域上的空间变异性；其中/>为水平梯度，“:”代表双点乘运算，α0和α1为待确定的参数；
将三个待确定的参数简化为两个有物理意义的参数，包括相关长度len和信噪比SNR；len表示插值点对周围分析场产生影响的距离，单位为千米；SNR为背景误差协方差和观测误差协方差之比，SNR愈大，观测场相对背景场更加准确，需要对比待融合的资料，由外界输入。
2.根据权利要求1所述的一种基于多源卫星数据融合的北极海冰厚度空间分辨率改进方法，其特征在于，所述步骤S2中对格点数据与实测值进行时间匹配和空间匹配，具体包括：
步骤S2.1、对各海冰厚度产品进行时间匹配；
针对CS2SMOS，选取观测当天、观测前3天及观测后3天的所有观测数值，取均值作为当天的观测值；将CPOM由月平均化的格点数据展开为日数据；选取ULS的日平均值作为该点的海冰厚度；针对IceBridge和IMB，选取每个采样点采样时间所在的日期为该点的采样日期；
步骤S2.2、基于小区域平均和降采样方式进行空间匹配；针对ULS观测值，以取样点为圆心，针对CS2SMOS、APP-x和PIOMAS产品，以25km为半径，针对CPOM产品，以5km为半径，圆内的格点观测值的平均值作为对应ULS观测处的产品观测值；针对CS2SMOS、PIOMAS和APP-x产品，以取样点为圆心，将半径25km内的格点观测值的平均值作为对应观测处的产品观测值；针对CPOM产品，以取样点为圆心，将半径5km内的格点观测值的平均值作为对应观测处的产品观测值；针对IMB和IceBridge，实测值空间分布密集观测冗余，为减少计算量，先确定IMB或IceBridge实测值的空间覆盖区域，取覆盖区域外扩1°为待计算区域，以待计算区域的所有产品观测格点为中心，对区域内的实测值以反比例加权的方式与海冰产品进行空间匹配，通过对冗余实测值加权平均，实现对实测数据的降采样。
3.根据权利要求2所述的一种基于多源卫星数据融合的北极海冰厚度空间分辨率改进方法，其特征在于，所述步骤S2.2中空间匹配具体数据处理方法包括：
步骤S2.2.1、针对IceBridge或IMB资料，以四类海冰厚度产品在北极的空间格点(xpi,ypi)k为基准，计算与当日每个实测点(xoj,yoj)k的球面距离Dijk如下：
其中，i＝1,…,m，j＝1,…,n，k＝1,…,t；m为海冰厚度产品在北极空间覆盖的点数，n为IceBridge或IMB在第k天探测到的空间点数，t为有实测数据的日期数目；
当球面距离Dijk满足以下条件：
Dijk≤L
则认为(xoj,yoj)k处于第k日以(xpi,ypi)k为圆心，以L为半径的圆内；其中当产品选择CS2SMOS、APP-x或PIOMAS时，L＝25km，当产品选择CPOM时，L＝5km；
根据不同海冰厚度产品，选取不同L，将满足上述条件的Dijk对应的海冰厚度实测值依据Dijk的大小进行反距离加权平均，将均值结果作为(xpi,ypi)k产品格点处的实测值，实现降采样过程；
步骤S2.2.2、针对ULS，以实测点(xoj,yoj)k，j＝1,2,3为圆心，依据步骤S2.2.1所述球面距离的计算公式，寻找半径为25km内包含的s个空间点(xpi,ypi)k，i＝1,…,s；将对应的s个观测值求均值，作为第k天在(xoj,yoj)k处的四类产品的观测值。
4.根据权利要求3所述的一种基于多源卫星数据融合的北极海冰厚度空间分辨率改进方法，其特征在于，所述步骤S2.2.1中各类海冰厚度产品在北极覆盖的空间点数m选取分别如下：
对于PIOMAS产品，m＝360×120；对于APP-x产品，m＝361×361，对于CS2SMOS产品，m＝432×432，对于CPOM产品，m＝432×432；
根据不同算力情况，基于(xoj,yoj)k的覆盖范围，可以减少经纬度范围，m取值随k的取值变化。
5.根据权利要求1所述的一种基于多源卫星数据融合的北极海冰厚度空间分辨率改进方法，其特征在于，步骤S4中通过波数谱对比各产品的空间分辨率，具体包括：
使用经向谱，沿区域内每条经度线上的序列数据求取波数谱，并进行平均；首先将序列数据去除趋势性，然后运用Black-Harris窗，防止谱能量泄露，通过Burg法估计融合模型的波数谱。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111400996.1A CN114114358B (zh) | 2021-11-24 | 2021-11-24 | 一种基于多源卫星数据融合的北极海冰厚度空间分辨率改进方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111400996.1A CN114114358B (zh) | 2021-11-24 | 2021-11-24 | 一种基于多源卫星数据融合的北极海冰厚度空间分辨率改进方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114114358A CN114114358A (zh) | 2022-03-01 |
| CN114114358B true CN114114358B (zh) | 2024-05-28 |
Family
ID=80440793
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202111400996.1A Active CN114114358B (zh) | 2021-11-24 | 2021-11-24 | 一种基于多源卫星数据融合的北极海冰厚度空间分辨率改进方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114114358B (zh) |
Families Citing this family (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114239704A (zh) * | 2021-12-02 | 2022-03-25 | 集美大学 | 基于海冰密集度数据的观测站点融化和冻结时间提取方法 |
| CN115326208B (zh) * | 2022-07-29 | 2024-07-02 | 武汉大学 | 基于被动微波多系点反演海冰厚度的方法及系统 |
| CN117892055B (zh) * | 2024-01-18 | 2024-07-09 | 黄河水利委员会水文局 | 一种基于曲线拟合和奈奎斯特定理的水温数据计算方法 |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2005265733A (ja) * | 2004-03-19 | 2005-09-29 | National Institute Of Information & Communication Technology | Ｓａｒによる氷厚推定方法 |
| RU2460968C1 (ru) * | 2011-03-22 | 2012-09-10 | Федеральное государственное бюджетное учреждение "Арктический и Антарктический научно-исследовательский институт" (ФГБУ "ААНИИ") | Способ определения высоты снежного покрова на льду акваторий |
| CN109765550A (zh) * | 2019-01-17 | 2019-05-17 | 中国人民解放军61741部队 | 海冰厚度反演方法、系统及电子设备 |
| CN112434814A (zh) * | 2020-12-07 | 2021-03-02 | 中国人民解放军国防科技大学 | 一种基于多源异构信息融合算法对航运经济潜力的分析方法 |
| CN112504144A (zh) * | 2020-12-04 | 2021-03-16 | 南京大学 | 一种海冰表面积雪厚度的遥感估算方法 |
| CN113063360A (zh) * | 2021-03-15 | 2021-07-02 | 上海工程技术大学 | 一种基于单光子激光测高数据的海冰厚度估算方法 |
-
2021
- 2021-11-24 CN CN202111400996.1A patent/CN114114358B/zh active Active
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2005265733A (ja) * | 2004-03-19 | 2005-09-29 | National Institute Of Information & Communication Technology | Ｓａｒによる氷厚推定方法 |
| RU2460968C1 (ru) * | 2011-03-22 | 2012-09-10 | Федеральное государственное бюджетное учреждение "Арктический и Антарктический научно-исследовательский институт" (ФГБУ "ААНИИ") | Способ определения высоты снежного покрова на льду акваторий |
| CN109765550A (zh) * | 2019-01-17 | 2019-05-17 | 中国人民解放军61741部队 | 海冰厚度反演方法、系统及电子设备 |
| CN112504144A (zh) * | 2020-12-04 | 2021-03-16 | 南京大学 | 一种海冰表面积雪厚度的遥感估算方法 |
| CN112434814A (zh) * | 2020-12-07 | 2021-03-02 | 中国人民解放军国防科技大学 | 一种基于多源异构信息融合算法对航运经济潜力的分析方法 |
| CN113063360A (zh) * | 2021-03-15 | 2021-07-02 | 上海工程技术大学 | 一种基于单光子激光测高数据的海冰厚度估算方法 |
Non-Patent Citations (2)
| Title |
|---|
| Feasibility of the Northeast Passage: The role of vessel speed, route planning, and icebreaking assistance determined by sea-ice conditions for the container shipping market during 2020–2030;Yangjun Wang et al.;《Transportation Research Part E》;20210311;第1-24页 * |
| 遥感海表盐度分析产品的有效分辨率研究;陈建;《海洋学研究》;20181231;第36卷(第4期);第17-26页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114114358A (zh) | 2022-03-01 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN114114358B (zh) | 一种基于多源卫星数据融合的北极海冰厚度空间分辨率改进方法 | |
| Ishida et al. | Use of cokriging to estimate surface air temperature from elevation | |
| Chatterjee et al. | A new atlas of temperature and salinity for the North Indian Ocean | |
| CN111401602B (zh) | 基于神经网络的卫星以及地面降水测量值同化方法 | |
| Mao et al. | An empirical orthogonal function model of total electron content over China | |
| Li et al. | Mapping near-surface air temperature, pressure, relative humidity and wind speed over Mainland China with high spatiotemporal resolution | |
| CN112505068A (zh) | 一种基于gnss-ir的地表土壤湿度多星组合反演方法 | |
| CN112861072B (zh) | 一种星地多源降水自适应动态融合方法 | |
| White et al. | Space/time statistics of short‐term climatic variability in the western North Pacific | |
| CN116644379A (zh) | 多源海面物理要素的机器学习融合方法、设备及介质 | |
| CN113553766A (zh) | 一种使用机器学习反演北极积雪深度的方法 | |
| CN113281754A (zh) | 一种雨量站融合雷达定量估测降雨的WRF-Hydro关键参数率定方法 | |
| Abushandi et al. | Rainfall estimation over the Wadi Dhuliel arid catchment, Jordan from GSMaP_MVK+ | |
| Beaucage et al. | Synthetic aperture radar satellite data for offshore wind assessment: A strategic sampling approach | |
| CN111597692B (zh) | 一种地表净辐射估算方法、系统、电子设备和存储介质 | |
| González-Gambau et al. | First SMOS Sea surface salinity dedicated products over the Baltic Sea | |
| Verron et al. | Assimilation of geosat data into a quasi-geostrophic model of the north-atlantic between 20-degrees and 50-degrees-N-preliminary-results | |
| Fioravanti et al. | Interpolating climate variables by using INLA and the SPDE approach | |
| CN109886497B (zh) | 基于纬度改进的反距离加权法的地面气温插值方法 | |
| CN112444825A (zh) | 一种基于北斗geo的电离层tec同化建模方法 | |
| CN112561140A (zh) | 基于东亚副热带急流和极锋急流协同变化的中国四季降水预测方法 | |
| Peng et al. | Validation of wind speeds from brown-peaky retracker in the Gulf of Mexico and east coast of north America | |
| Handoko et al. | The ENSO’s Influence on the Indonesian Sea Level Observed Using Satellite Altimetry, 1993–2016 | |
| CN114004426B (zh) | 一种短时暴雨预报释用模型的动态调整方法 | |
| Vujec et al. | Kalman filter sensitivity tests for the NWP and analog-based forecasts post-processing |
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 |