CN110046415A - 一种时空精细化的土壤有机质含量遥感动态反演方法 - Google Patents

Abstract

本发明涉及一种时空精细化的土壤有机质含量遥感动态反演方法。选择Sentinel‑2和Sentinel‑3作为遥感影像源，利用Gram‑Schmidt光谱锐化图像融合方法进行影像融合，保证了该方法使用的影像数据兼具高空间、高时间、高光谱分辨率；构建包含融合后影像多波段反射率、多植被指数的反演模型框架，通过关系拟合、遴选比对、自相关校验等方法在指标遴选池中选取最优指标组合方式，构建形成一种基于影像融合、最能捕捉土壤有机质时空变化特征的影像反演方法。本发明的土壤有机质含量遥感反演方法兼顾高光谱分辨率和高时空精度，可获得区域尺度内时空精细化的土壤有机质含量变化图，为探究农田产能的变化趋势、从而精准指导农业生产提供重要的基础数据与决策支撑。

Description

一种时空精细化的土壤有机质含量遥感动态反演方法

技术领域

本发明属于遥感技术领域，具体涉及一种兼顾高光谱分辨率和高时空精度的土壤有机质含量遥感反演方法。

背景技术

土壤有机质是指土壤中的动植物残体、微生物及其分解与合成的物质，它是反映土壤肥力的重要指标同时也是响应生态环境变化的重要参数。一方面，利用多光谱遥感数据对有机质进行定量反演具有低成本，省时高效且时效性强的优势；另一方面，由于遥感影像在时间、空间、光谱分辨率上难以同时达到最优，空间分辨率高的影像往往时间分辨率较低，不能将研究时间细化，时间分辨率较高的影像通常空间分辨率较低，小范围区域反演效果较差。因此目前小区域有机质的研究时间相对静态，以反演一期为主。而有机质易受施肥、湿度、温度等因素影响，变化敏感，所以时间尺度上的精细化动态监测是必要的，这为推进精准农业的发展提供重要的数据资料成果。

遥感影像融合是将在空间、时间、波谱上冗余或互补的多源遥感数据按照一定的算法进行运算处理，获得比任何单一数据更精确、更丰富的信息，生成具有新的空间、波谱、时间特征的合成影像数据。影像通过融合既可以提高多光谱影像空间分辨率，又保留其多光谱特性。因此，它不仅是数据间的简单复合，而且强调信息的优化，以突出有用的专题信息，消除或抑制无关的信息，改善目标识别的影像环境，从而增强解译的可靠性，减少不确定性及误差，提高分类精度，扩大应用范围和效果。融合方法包括主分量变换、乘积变换、BROVEY变换融合法、Gram-schmidt等。

土壤的发育是一个由成土母质、气候条件、微生物环境、地形条件等成土因素以及人类活动影响综合作用的过程。在这一过程中，成土时间的长短和外部环境也会不同程度地影响土壤的发育。不同的成土因素会使土壤的物理和化学性质，如土壤有机质含量产生各种差异，这种差异性决定了每一类土壤都有其独特的光谱特征，它是土壤属性定量反演的物理基础^[1,2]。研究发现，有机质对土壤光谱的影响主要在红光波段与近红外波段。刘焕军利用Landsat8和Sentinel-2A卫星红波段建立基于BP神经网络的黑龙江黑土带SOM遥感反演模型^[3]。孙建英选取中国北方潮土作为研究对象，分析结果表明土壤有机质含量可以利用田间土样的近红外光谱特性进行分析和检测，建立的线性模型是有效的^[4]。而植被指数可削弱环境背景对光谱的干扰，比单波段具有更好灵敏性^[5～7]，它与有机质的相关性也一直被探讨，李燕丽利用AVHRR和MODIS NDVI遥感影像数据研究发现，广西南宁、崇左地区NDVI变化与土壤有机质变化呈显著正相关，NDVI时间序列数据能够在一定程度上反映土壤有机质的变化^[8]。余博应用地统计学方法，研究发现伊犁绢蒿荒漠草地土壤有机质和NDVI的异质性具有相关性^[9]。Ayoubi S等利用人工神经网络发现伊朗中西部的牧场和林区Landsat ETM波段1、2以及NDVI与土壤有机质相关性较高^[10]。Pavan Kumar利用Sentinel-2A卫星的可见光和近红外反射率数据计算得到NDVI、SAVI、RDVI，再通过回归分析得到通过卫星数据测量的SOC与叶面积指数(LAI)之间的关系，相关系数R²为0.95(p值<0.05)^[11]。

目前土壤有机质的遥感反演方法与研究较多，但在中小区域尺度而言，都面临一个共性问题：现有方法多以Landsat影像或Modis传感器为数据源开展研究，其中Modis影像空间分辨率较低，在中小区域尺度难以有效把握空间异质性信息，而Landsat虽然空间分辨率较高，但受云覆盖、过境时间等影响，在一定的时间范围内难以同时获取多期可用的数据，且波段数较少，这就造成了当前的有机质反演方法在时间上相对静态，对于时空异质性的把握不足。而实际上，土壤有机质是体现土壤肥力最重要、最直接的指示性指标，其相比于氮磷等营养元素而言，对施肥、湿度、温度等因素影响更为敏感，因此在不同尺度上的变化也更为剧烈。由此可见，开展土壤有机质时空精细化信息的提取和空间制图十分必要，而现有的遥感反演方法集中于单一的数据源(时间、空间、光谱分辨率难以兼顾)与单一的波段信息(无法削弱环境背景对光谱的干扰，灵敏性不高)。本发明提供的一种基于影像融合的区域尺度土壤有机质反演方法，可以有效解决上述问题。

参考文献

[1]高永光,胡振琪.高光谱遥感技术在土壤调查中的应用闭.矿业研究与开发,2006,(01):44-46.

[2]吕国楷,洪启旺,郝允充,等.遥感概论[M].北京:高等教育出版社,1995.

[3]刘焕军,潘越,窦欣,等.黑土区田块尺度土壤有机质含量遥感反演模型[J].农业工程学报,2018,34(01):127-133.

[4]孙建英,李民赞,郑立华,等.基于近红外光谱的北方潮土土壤参数实时分析[J].光谱学与光谱分析,2006(03):426-429.)

[5]Jan Pisek,Jing M Chen.Comparison and validation of MODIS andVEGETATION global LAI products over four BigFoot sites in North America[J]Remote Sensing of Environment,2007,109(1):81-94

[6]Sims,D.and Gamon,J.A.Relationships between leaf pigment contentand spectral reflectance across a wide range of species,leaf structures anddevelopmental stages[J].Remote Sensing of Environment,2002 81(2):337-354.

[7]郑曼迪,熊黑钢,乔娟峰,等.基于宽波段与窄波段综合光谱指数的土壤有机质遥感反演[J].激光与光电子学进展,2018,55(07):457-465.

[8]李燕丽,潘贤章,周睿,等.长期土壤肥力因子变化及其与植被指数耦合关系[J].生态学杂志,2013,32(03):536-541.

[9]余博,朱进忠,范燕敏,等.伊犁绢蒿荒漠草地土壤养分和植被指数的空间异质性研究[J].新疆农业科学,2009,46(06):1294-1300.

[10]Ayoubi S,Shahri AP,Karchegani PM,et al.2011.Application ofartificial neural network(ANN)to predict soil organic matter using remotesensing data in two ecosystems.Biomass and Remote Sensing of Biomass.Croatia:InTech:181－196.

[11]Kumar,P.,Sajjad,H.,Tripathy,B.R.,Ahmed,R.,&Mandal,V.P..2017.Prediction of spatial soil organic carbon distribution usingsentinel-2a and field inventory data in sariska tiger reserve.NaturalHazards.

[12]Bruno Aiazzi,Stefano Baronti,Massimo Selva.Enhanced Gram-SchmidtSpectral Sharpening Based[J].IEEE Geoscience and Remote Sensing Symposium,2006,6:3806-3809.

[13]C.A.Laben,B.V.Brower.Process for enhancing the spatial resolutionof multispectral imagery using Pan-sharpening[M].U.S.Patent Tech.Rep.,EastmanKodak Company,2000,4.

[14]李存军,刘良云,王纪华,等.两种高保真遥感影像融合方法比较[J].中国图像图形学报,2004，9(11):1376-1385.

发明内容

本发明提供一种兼顾高光谱分辨率和高时空精度的土壤有机质含量遥感反演方法。可获得区域尺度内时空精细化的土壤有机质含量变化图，为探究农田产能的变化趋势、从而精准指导农业生产提供重要的基础数据与决策支撑。

为达成上述目的，本发明所采用的技术方案如下：

一种土壤有机质含量遥感反演方法，包括如下步骤：

(1)采集遥感影像并进行预处理，所述遥感影像为无云或少云遮盖的Sentinel-2和Sentinel-3影像；

(2)采用Gram-Schmidt光谱锐化方法对Sentinel-2和Sentinel-3影像进行融合，获取空间分辨率增强的多光谱影像；

(3)采集土壤样品，测定各样品溶解性有机碳TOC含量，作为训练数据；

(4)采取由多植被指数与多波段反射率构成的最优化指标进行反演模型构建，如下所示：

Y＝AX₁+BX₂+C (1)

式中X₁为波段反射率，X₂为植被指数；

在X₁与X₂构成的指标池中进行指标遴选，获取与TOC含量相关系数最高的指标，作为最优指标；

(5)从训练数据中选取实测数据作为计算集，选用步骤(4)确定的最优指标代入式(1)，利用SPSS建立多元逐步线性回归模型，获取反演模型。

作为本发明的进一步改进，所述步骤(1)中，所述Sentinel-2影像产品类型选择S2MSI1C；Sentinel-3影像产品类型选择OL_1_EFR_，探测器选择OLCI。

作为本发明的进一步改进，所述步骤(1)中，Sentinel-2影像经过ENVI 5.5辐射定标和FLAASH大气校正处理；Sentinel-3影像经过ENVI 5.5几何定位、辐射定标和FLAASH大气校正处理。

作为本发明的进一步改进，所述步骤(2)中，高分辨率全色波段图像选择Sentinel-2影像第8波段，低分辨率多光谱图像选择Sentinel-3影像21个波段，利用多光谱波段的平均值来模拟低分辨率的全色波段。

作为本发明的进一步改进，所述步骤(3)中，在过境影像当天，采用梅花形布点法垂直采集0-20cm的表层土壤样品，土壤样品经室内自然风干至恒重，去除杂质，敲碎研磨过2mm尼龙筛，保存备用；土壤TOC的测定采用重铬酸钾容量法。

作为本发明的进一步改进，所述步骤(4)中，X₁从中心波长涵盖400nm～1020nm的21个波段中进行选择；X₂从由近红外与红外波段组合而成的归一化植被指数NDVI、比值植被指数RVI、差值植被指数DVI、优化土壤调节植被指数OSAVI、重归一化植被指数RDVI和增强型植被指数EVI中进行选择。

通过遍历算法对需要选取的植被指数、植被指数中近红外与红外具体波段号进行遴选，与对应样点土壤样品TOC含量进行决定系数R²与绝对误差RMSE比较选择最优指标。

指标池最优指标构建原则为：

a)X₁与X₂中选取的波段无重复；

b)X₁与X₂中选取的波段间不具有共线性；

c)波段尽可能多的反映波谱信息。

作为本发明的进一步改进，所述步骤(5)中，选择训练数据的部分数据作为计算集建立反演模型，其余部分用于对构建的反演模型进行验证，通过R²与RMSE进行精度评价。进一步的，所述步骤(5)还包括：将融合后影像基于上述模型计算得到TOC空间分布图，并通过已去除反演地类的土地利用数据对TOC空间分布图进行掩膜去除。

4)时间推演

在2016年11月至2018年11月间形成了可用于东海县土壤TOC含量反演的52期影像。

本发明的土壤有机质含量遥感反演方法在现有土壤有机质遥感动态估算的基础上，兼顾了时空精细度与光谱分辨率的提升，可获得区域尺度内时空精细化的土壤有机质含量变化图，为探究农田产能的变化趋势、从而精准指导农业生产提供重要的基础数据与决策支撑。

应当理解，前述构思以及在下面更加详细地描述的额外构思的所有组合只要在这样的构思不相互矛盾的情况下都可以被视为本公开的发明主题的一部分。另外，所要求保护的主题的所有组合都被视为本公开的发明主题的一部分。

从下面的描述中可以更加全面地理解本发明教导的前述和其他方面、实施例和特征。本发明的其他附加方面例如示例性实施方式的特征和/或有益效果将在下面的描述中显见，或通过根据本发明教导的具体实施方式的实践中得知。

附图说明

图1是东海县某研究区土壤采集点示意图；

图2是Gram-Schmidt图像融合流程图；

图3是影像融合前后效果图；

图4是模型选择示意图；

图5是建模集与验证集拟合结果；

图6是掩膜前后反演效果图；

图7是东海县土壤TOC含量反演的20期典型分布图；

图8是东海县土壤TOC含量反演的52期波动曲线。

具体实施方式

为详细了解本发明的技术内容，特举具体实施例并配合所附图式说明如下。

在本公开中参照附图来描述本发明的各方面，附图中出示了许多说明的实施例。本公开的实施例不必定义在包括本发明的所有方面。应当理解，上面介绍的多种构思和实施例，以及下面更加详细地描述的那些构思和实施方式可以以很多方式中任意一种来实施，这是因为本发明所公开的构思和实施例并不限于任何实施方式。另外，本发明公开的一些方面可以单独使用，或者与本发明公开的其他方面的任何适当组合来使用。

实施例1

本实施例以图1所示的东海县某研究区土壤TOC含量反演方法为例，具体说明本发明遥感监测方法的具体实现方式。

本实施例基于融合遥感影像数据对土壤有机质含量进行定量反演，实现方式如下：在基于Gram-Schmidt(GS)融合同时期Sentinel-2的高空间分辨率信息和Sentinel-3丰富的光谱信息基础上，结合东海县两期实测土壤TOC含量数据，其中92个点建模，利用融合后影像中提取的采样点的反射率数据计算采样点像元的六个植被指数(NDVI、RVI、DVI、OSAVI、RDVI、EVI)，选取与TOC相关性最高的RVI以及与TOC相关性较好的第20波段利用多元逐步线性回归分析方法建立模型，再将(AX₁+BX₂+C)的值与土壤TOC分别进行线性、指数、对数、幂函数拟合选择最优模型。以15个TOC实测数据对模型进行验证。将该模型推广到其他时间的土壤TOC含量反演，得到多期土壤TOC空间分布图，来反演东海县土壤TOC含量，最后得到从2016.11年至2018.11共52期土壤TOC空间分布图，实现精细化时空动态反演，可获取区域土壤年际、月际、季度等不同时间尺度下的有机质变化规律及其精细化的空间分布。为精准指导农业生产提供重要的基础数据。

作为示例性的描述，下面结合图表所示，对前述方法的实施进行具体说明。

一、兼顾高时间频率与高光谱分辨率的Sentinel同步影像融合

1)影像采集与预处理

从欧洲航天局(ESA)官网上(https://scihub.copernicus.eu/)获取影像数据，选择覆盖研究区的无云或尽可能少云遮盖的Sentinel-2和Sentinel-3影像。Sentinel-2产品类型选择S2MSI1C；Sentinel-3产品类型选择OL_1_EFR_，探测器选择OLCI。下载了从2016年11月至2018年11月共104景包含东海县行政区的影像；

Sentinel-2影像经过ENVI 5.5辐射定标和FLAASH大气校正处理；Sentinel-3影像经过ENVI 5.5几何定位、辐射定标和FLAASH大气校正处理。气溶胶模型选择Urban，气溶胶反演方法选择2-band(K-T)，研究区平均高程为0.25km。

2)影像融合

本发明融合方法为Gram-Schmidt光谱锐化方法(GS)，它是由基于代数运算的方法——Gram-Schmidt正交化发展而来。，正交化过程如公式1所示：

其中，{u₁，u₂，…，u_n}是一组相互独立的向量，{v₁，v₂，…，v_n}为通过Gram-Schmidt构造的正交向量。

Gram-Schmidt光谱锐化融合在Sentinel-2和Sentinel-3影像中应用的具体步骤如图2(Laben等，2000)：

Gram-Schmidt光谱锐化图像融合的关键步骤可表述如下：

使用多光谱低空间分辨率影像对高分辨率波段影像进行模拟(图1(1))。

为达到较好的空间信息增强效果，又保证较高的光谱保真度，改进融合图像的质量，利用高光谱图像的相关波段进行光谱重采样，以模拟全色的分辨率图像。

模拟的高分辨率波段影像信息量特性与高分辨率全色波段影像的信息量特征比较接近。模拟的高分辨率波段影像在后面的处理中被作为Gram-Schmidt第一分量进行GS变换。由于在GS变换中第一分量GS₁没有变换，故模拟的高分辨率波段影像将被用来与高分辨率全色波段影像进行交换，这样可使信息失真少。

利用模拟的高分辨率波段影像作为GS变换的第一个分量来对模拟的高分辨率波段影像和低分辨率波段影像进行GS变换(图1(2))。

该算法在运用到图像处理时，对GS变换进行了修改，将第T个GS分量由前T-1个GS分量构造，具体修改如公式(2)(李存军等，2004)：

其中，GS_T是GS变换后产生的第T个分量，B_T是原始多光谱影像的第T个波段影像，u_T是第T个原始多光谱波段影像灰度值的均值。

通过调整高分辨率波段影像的统计值来匹配GS变换后的第一个分量GS₁，以产生经过修改的高分辨率波段影像。修改方式见(4)(5)。该修改有助于保持原多光谱波段影像的光谱特征。

将经过修改的高分辨率波段影像替换GS变换后的第一个分量，产生一个新的数据集。将新的数据集进行反GS变换(图1(3))，即可产生空间分辨率增强的多光谱影像。

GS反变换的公式如下(李存军等，2004)：

GS方法有效去除信息冗余，能保持较好的空间纹理信息，不受波段限制，信噪比高，信息失真小，可在较大程度上保持低空间分辨率影像的光谱特征。融合日期相同或相近的Sentinel-2影像与Sentinel-3影像，融合操作在ENVI Classic 5.3Gram-SchmidtSpectral Sharpening模块中进行，在Select Low Spatial Resolution Multi BandInput File对话框中选择Sentinel-3影像，在Select High Spatial Resolution PanInput Band对话框中选择Sentinel-2影像第8波段，利用多光谱波段的平均值来模拟低分辨率的全色波段。融合后影像重访周期为5天，空间分辨率为10米，涵盖21个波段，形成了面向土壤TOC反演光谱信息丰富、重访时间较短的高空间分辨率影像。影像融合前后效果如图3所示。

二、适宜土壤有机质反演的关键指标遴选

实地参数获取(实地采样、试验指标)

土壤样点采集工作分两期进行，一期于2017年12月底开展，土壤样品采集位于整治项目实施区温泉镇，采集样品15件；二期于2018年3月上旬开展，土壤样品采集分布于全县，共采集样品92件。采样期间天气晴朗，采样点土地利用类型主要是耕地、林地、园地等，主要植被类型为小麦、玉米、桃树。

采用梅花形布点法垂直采集0-20cm的表层土壤。在每个采样点采样时，先把采样区域划分为九宫格，然后采集该区域内的四个角以及中心点的土壤样品，接着把采集的5个土壤样品等量混合成一份，这个过程去除体积较大的地表砾石以及动植物残体，其余土壤用四分法弃去，留下1kg左右样品。最后，将土壤样品装入干净的布袋，记好编号。土壤样品经过室内自然风干至恒重，去除杂质，敲碎研磨过2mm尼龙筛，研磨好的土壤样品装袋保存以备用。土壤TOC的测定采用重铬酸钾容量法，在过量的硫酸存在下，用氧化剂重铬酸钾氧化有机碳，剩余的氧化剂用标准硫酸亚铁溶液回滴，从消耗的氧剂量来计算有机碳量。

1)构建适用于土壤TOC反演的模型框架

利用ENVI 5.3Input Points from ASCII模块从融合后影像中获取采样点的反射率。模型中若能包含植被指数以及与TOC相关性较好的特征波段，则可以容纳更多的表征TOC的波段信息，提高模型的反演精度。对融合后影像采取由多植被指数与多波段反射率构成的最优化指标进行反演模型构建，其表达形式为：

Y＝AX₁+BX₂+C (7)

模型中X₁为波段反射率，X₂为植被指数；模型构建原则：

a)模型内波段不能重复；

b)波段间不能具有共线性；

c)波段应尽可能多的反映波谱信息。

2)遥感指标集的选择

模型中波段反射率X₁从中心波长涵盖400nm～1020nm的21个波段中进行选择，融合后影像的第3波段为蓝光波段，第8～12波段为红光波段，第13～21波段为近红外波段。植被指数X₂从归一化植被指数(NDVI)、比值植被指数(RVI)、差值植被指数(DVI)、优化土壤调节植被指数(OSAVI)、重归一化植被指数(RDVI)、增强型植被指数(EVI)中进行选择，植被指数计算采用遍历算法。利用SPSS 20Pearson相关性分析统计采样点像元的21个波段以及七个植被指数与TOC的相关性，选择相关性最高的植被指数与波段建模。

表1遥感指标集选择

通过波段反射率和植被指数两个指标进行各种模型组合方法的量测，以东海县土壤92个样点为示例。结果发现，以B₁₈作为近红外代表，B₉作为红外代表组成的RVI(即B₁₈/B₉)与土壤TOC相关系数最高，为0.641，RMSE最小，为1.3926，将RVI_18/9作为X₂代表。同时8,9,10,20四个波段反映土壤TOC的波谱信息最为丰富，相关系数分别为-0.520、-0.535、-0.531、-0.43，但是考虑到波段间不能重复且应尽量避免共线性的原则，而8，9，10具有极强的互相关，存在共线性，相关系数达到1左右，故选取B₂₀作为X₁代表。选用RVI_18/9与B₂₀波段利用SPSS建立多元逐步线性回归模型得到AX₁+BX₂+C。

表2 21个波段与TOC的相关性

表3 RVI与TOC的相关性

表4波段8、9、10相关性

3)关系比对与遴选(关系对比)

将(AX₁+BX₂+C)的值与土壤TOC分别进行线性、指数、对数、幂函数拟合，结果如图4所示。线性拟合R²＝0.6396，RMSE＝0.3761；指数拟合R²＝0.6149，RMSE＝0.4001；幂函数拟合R²＝0.5446，RMSE＝0.3810；对数拟合R²＝0.5270，RMSE＝0.3888，线性模型拟合R²最大，RMSE最小，拟合效果最好。

表5关系比对与遴选

三、构建基于Sentinel-2与Sentinel-3同步卫星影像融合数据的土壤有机质时空精细化反演方法

1)模型构建

选取92个点作为计算集，选用RVI_18/9与B₂₀波段利用SPSS建立多元逐步线性回归模型，再将(AX₁+BX₂+C)的值与土壤TOC分别进行线性、指数、对数、幂函数拟合，结果显示线性模型拟合R²最大，RMSE最小，拟合效果最好，故确定有机质反演模型。模型如以下公式所示：

2)精度验证

通过决定系数(R²)、均方根误差(RMSE)来评价模型反演TOC的优劣。R²越大，说明模型越稳定；RMSE越小，则模型精度越高。选取2017年12月的15个样点进行验证，本发明模型：建模集R²＝0.6396，RMSE＝0.4018；验证集R²＝0.6595，RMSE＝0.1495，模型拟合效果较好，结果如图5所示。

3)掩膜去除

将52期融合后影像均利用ENVI 5.3Band Math模块使用上述模型计算得到TOC空间分布图，并通过已去除水体、建设用地等无法反演地类的2015年土地利用数据对TOC空间分布图进行掩膜，过程下图所示。最后利用ArcGIS制图工具绘制52期研究区TOC空间分布图，从而实现土壤有机质含量遥感精细化定量反演，图6右图为52期均值效果图。

4)时间推演

在2016年11月至2018年11月间形成了可用于东海县土壤TOC含量反演的52期影像。从52期中选取典型的20期分布图如图7所示，52期波动曲线如图8所示。

虽然本发明已以较佳实施例揭露如上，然其并非用以限定本发明。本发明所属技术领域中具有通常知识者，在不脱离本发明的精神和范围内，当可作各种的更动与润饰。因此，本发明的保护范围当视权利要求书所界定者为准。

Claims (10)

1.一种时空精细化的土壤有机质含量遥感动态反演方法，其特征在于，包括如下步骤：

(1)采集遥感影像并进行预处理，所述遥感影像为无云或少云遮盖的Sentinel-2和Sentinel-3影像；

(2)采用Gram-Schmidt光谱锐化方法对Sentinel-2和Sentinel-3影像进行融合，获取空间分辨率增强的多光谱影像；

(3)采集土壤样品，测定各样品溶解性有机碳TOC含量，作为训练数据；

(4)采取由多植被指数与多波段反射率构成的最优化指标进行反演模型构建，如下所示：

Y＝AX₁+BX₂+C (1)

式中X₁为波段反射率，X₂为植被指数；

在X₁与X₂构成的指标池中进行指标遴选，获取与TOC含量相关系数最高的指标，作为最优指标；

(5)从训练数据中选取实测数据作为计算集，选用步骤(4)确定的最优指标代入式(1)，利用SPSS建立多元逐步线性回归模型，获取反演模型。

2.根据权利要求1所述的方法，其特征在于，所述步骤(1)中，所述Sentinel-2影像产品类型选择S2MSI1C；Sentinel-3影像产品类型选择OL_1_EFR_，探测器选择OLCI。

3.根据权利要求1所述的方法，其特征在于，所述步骤(1)中，Sentinel-2影像经过ENVI5.5辐射定标和FLAASH大气校正处理；Sentinel-3影像经过ENVI5.5几何定位、辐射定标和FLAASH大气校正处理。

4.根据权利要求1所述的方法，其特征在于，所述步骤(2)中，高分辨率全色波段图像选择Sentinel-2影像第8波段，低分辨率多光谱图像选择Sentinel-3影像21个波段，利用多光谱波段的平均值来模拟低分辨率的全色波段。

5.根据权利要求1所述的方法，其特征在于，所述步骤(3)中，在过境影像当天，采用梅花形布点法垂直采集0-20cm的表层土壤样品，土壤样品经室内自然风干至恒重，去除杂质，敲碎研磨过2mm尼龙筛，保存备用；土壤TOC的测定采用重铬酸钾容量法。

6.根据权利要求1所述的方法，其特征在于，所述步骤(4)中，X₁从中心波长涵盖400nm～1020nm的21个波段中进行选择；X₂从由近红外与红外波段组合而成的归一化植被指数NDVI、比值植被指数RVI、差值植被指数DVI、优化土壤调节植被指数OSAVI、重归一化植被指数RDVI和增强型植被指数EVI中进行选择。

7.根据权利要求1或6所述的方法，其特征在于，所述步骤(4)中，通过遍历算法对需要选取的植被指数、植被指数中近红外与红外具体波段号进行遴选，与对应样点土壤样品TOC含量进行决定系数R²与绝对误差RMSE比较选择最优指标。

8.根据权利要求1所述的方法，其特征在于，所述步骤(4)中，指标池最优指标构建原则为：

a)X₁与X₂中选取的波段无重复；

b)X₁与X₂中选取的波段间不具有共线性；

c)波段尽可能多的反映波谱信息。

9.根据权利要求1所述的方法，其特征在于，所述步骤(5)中，选择训练数据的部分数据作为计算集建立反演模型，其余部分用于对构建的反演模型进行验证，通过R²与RMSE进行精度评价。

10.根据权利要求1所述的方法，其特征在于，所述步骤(5)还包括：将融合后影像基于上述模型计算得到TOC空间分布图，并通过已去除反演地类的土地利用数据对TOC空间分布图进行掩膜去除。