海水表面温度是一项最基本的海洋水文要素。海洋表面温度不仅影响全球气候变化，还和飓风以及海洋动植物的生长和分布都具有密切的联系，海洋表面温度异常会影响海洋动植物群落，改变海洋鱼群和植物的分布与数量，实践表明海洋渔场与特定的洋面温度分布有直接的关系。一些海洋动力变化特征和灾害性的海洋现象均与海洋表面温度的变化有着密切关系，因此，海表温度的遥感测量作为海表面温度的重要获取手段有着重要的应用价值和科学意义。

目前，温度反演主要通过遥感的方法，常用的遥感温度反演算法有劈窗算法（又称分裂窗算法）、单窗算法和多通道算法等^[1]。从1970起国外开始利用热红外遥感进行地表/海表温度反演研究，最早的方法是由Anding等^[2]提出的分裂窗方法，该方法成功应用于海洋表面温度反演具有较高的精度。McMillin根据大气在两个相邻的波段窗口具有不同的吸收特性，提出通过这两个通道辐射亮温的某种组合来消除大气影响的劈窗算法^[3]。Price^[4]利用大气辐射传输理论，对大气的影响作简化处理后改进劈窗算法，该算法通过两个通道对大气的不同吸收响应来部分地消除大气的影响。Becker和Li^[5]提出针对AVHRR数据的局地劈窗算法，Wan等^[7]进一步提出针对MODIS数据的通用劈窗算法，随后Becker等^[8]在Price 算法的基础上将劈窗算法进行改进^[8]。覃志豪等^[9]提出一种新的劈窗算法，该算法引入大气透过率和发射率两个重要参数，基于普朗克方程推导，该算法模型简单，而且可进一步提高温度反演精度。Trigo等采用通用劈窗算法实现MSG-SEVIRI数据的地表温度反演。覃志豪等^[10,11]针对Landsat 6 TM数据只有一个热红外波段的特点，根据地表热辐射传导理论提出了一个更简易的单窗算法，只需要考虑地表发射率、大气透射率和大气平均作用温度3个基本参数。Jimenez-Munoz等^[12]提出了一种普适性单通道算法，适用于不同的卫星传感器，所需的参数少，可用来从Landsat 6 TM数据中反演地表温度。Wan等^[13]针对 MODIS 数据提出了具有代表性的多通道温度反演算法。单窗算法是针对当时只存在一个热红外通道的TM数据提出的，需要的大气实时剖面数据难以准确获得，多通道算法需要白天晚上两景影像数据，所需参数多，运算复杂，反演难度比较大。劈窗算法所需参数少，运算速度快，精度相对较高，是目前研究最为成熟的算法，但劈窗算法需要地表发射率作为输入数据，对地表发射率误差非常敏感^[14]。翟劭燚等将劈窗算法应用于Landsat 8 TIRS数据进行地表温度反演，克服了传统地表温度反演依赖MODIS水汽产品外部数据造成的局限，实现利用Landsat 8热红外数据反演得到较高空间分辨率的地表温度产品^[15]。吴亮等^[16]利用劈窗算法对Landsat8遥感影像进行地表温度反演，对比了根据MODIS L1B数据和基于经验值的Landsat 8数据两种不同方式的LST反演结果，结果发现劈窗算法得到的反演温度更接近实际温度，精度较高且优于MODIS产品。毛克彪等^[17]针对MODIS的波段特点，通过对热辐射强度和温度之间的关系计算，对Planck方程进行线性简化，然后从MODIS的近红外波段反演大气水汽含量，通过建立大气水汽含量与热波段透过率的关系，从而计算透过率和比辐射率。通过MODIS的可见光波段、近红外和中红外波段数据，获得地表温度反演所需要的基本参数，从而提出了适合于MODIS数据的地表温度反演的劈窗算法。

本文利用一种形态学扩张的多阈值云检测方法和加权时空相关性去云方法来修复MODIS影像，使用简易的劈窗算法进行海表温度反演研究，依据大气辐射传输方程构建地表/海表温度反演方程，基于水体MODIS数据热红外通道特点、MODTRAN模拟的大气透过率与大气水汽含量关系以及DN值与亮度温度的转换方程分别求解出第31、32波段的比辐射率、大气透过率和亮度温度值，以实现海表温度反演；通过对同一时段的海表反演温度、ARGO浮标海表实测温度以及MODIS海表温度产品数据进行对比验证本算法反演精度；最后利用本算法反演得到研究区多时相海表温度数据，并对其时空变化特征进行深入分析。

研究区地理位置为8°~12° N，114°~120° E，如图1所示。该研究区被巴拉望岛分为东西两部分，其中巴拉望岛西部为南海海域范畴，东部为苏禄海海域范畴。南海是中国最大的外海，依据大陆架和中央海盆的分布，可以将中国南海分为东南西北四个海区，本研究区域主要是南海南部海域，靠近巴拉望岛。南海南部海域是连接西太平洋和东印度洋的重要通道，整个海域环绕众多岛屿，自北向南有菲律宾群岛、巴拉望群岛和加里曼丹岛等。由于其特殊的地理位置特征，以及海表温度的变化对周围岛屿和国家可能造成的气候影响，因此南海海表温度的变化规律一直受到各国气象学家重点关注，同时也是国内外学者重点研究内容。

本研究主要数据来源于搭载在Terra和Aqua两颗卫星上中分辨率成像光谱仪（MODIS）所获取的数据，MODIS具有36个中等分辨率水平（0.25~1 μm）的光谱波段，每1~2 d对地球表面观测一次，获取陆地和海洋温度、初级生产率、陆地表面覆盖、云、气溶胶、水汽以及火情等目标的图像，波谱范围覆盖0.4~14 μm，热红外波段有8个，MODIS第1~2波段分辨率为250 m，3~7波段分辨率为500 m，其他波段分辨率为1 000 m，其中第31和32波段与AVHRR的第4和5波段的波谱范围十分接近。因此其空间分辨率和时间分辨率都非常适用于地表/海表温度反演研究。Terra卫星过境时间为上午10:30，Aqua卫星的过境时间为下午13:30，一天最多可以获得4条过境轨道资料，虽然遥感观测为瞬时观测，闵文彬等将MODIS卫星反演的瞬时地表温度与月平均温度和年平均温度进行对比分析，发现其相关系数都在0.85以上，相关性较好，因此可用卫星反演温度来合成月平均温度^[18]。

采用的实测数据为ARGO浮标数据，ARGO浮标数据是指用于建立全球观测网的一种专用测量设备浮标数据，ARGO浮标可以在海洋中自由漂移，自动测量海面到2 000 m水深的海水温度、盐度以及深度, 并可通过移动轨迹，获取海水的移动速度和方向。ARGO浮标的温度测量范围为-3~32 ℃，测量精度在±0.005 ℃之间，温度分辨率为0.001 ℃。通过下载研究区对应影像范围内的浮标数据(http://www.argodatamgt.org)，利用Matlab软件读取海洋表面温度值，用于后续验证反演精度。

MODIS数据的预处理主要包括几何校正和影像去云。由于MODIS数据空间分辨率较低，并且研究区域主要是海域，因此无法选取控制点进行精确的几何校正，但是MODIS数据轨道参数已知，根据其轨道参数可建立图像校正的RPC模型并对图像进行几何校正。该方法的原理是首先求得每个像元的经纬度坐标，然后对图像进行重采样得到几何校正结果。通过对比几何校正后的MODIS影像与Google地图影像，两者地理位置信息相差在MODIS数据的2个像元以内，可以证明利用MODIS数据自带的经纬度信息进行几何校正较为准确。为了提高MODIS影像的质量，需要对几何校正后的影像进行去云处理。去云处理主要采用参考影像校正待去云影像的方法，对待去云影像应用形态学扩张的双阈值云识别算法检测影像中的云。首先使用一个较高阈值K1检测每个像元是否属于云，然后在所有云像元周围寻找大于低阈值K2的像元作为新的云像元。利用地物之间的时空相关性对影像进行去云处理。时空相关性影像去云方法前提假设是遥感影像上的地物在地理上存在很强的时间和空间相关性。MODIS影像含云量通常较高，一般情况下待去云影像和参考影像上都会被云覆盖，只是覆盖区域有所不同，同时MODIS影像分辨率低，不容易进行影像之间的配准，因此对MODIS影像去云会更加复杂。厚云层会使影像造成像元缺失问题，可以用空间插值的方法，如反距离加权法估算云覆盖下的像元值。反距离加权插值法是根据相近相似原则对估测点进行空间插值的，将像元缺失的区域像元视作估测点，将像元缺失区域周边的无云像元视作已知的样本点，以缺失像元和临近无缺失像元之间的距离D作为权重因子，公式如下所示：

其中:T₀是云覆盖像元缺失的像元估测值；T_i是无云像元i的像元值；n为空间插值所需的周边无云像元数；p是距离的幂函数，一般情况下可取p=2。修复后的影像薄云污染的像元明显较少，如图1所示。

基于MODIS数据的温度反演研究主要包括数据预处理、辐射定标处理、水汽含量反演、劈窗算法温度反演、精度验证和结果评价几个步骤。MODIS L1B数据是经过了辐射定标处理且获取定标后的反射率数据 、辐射亮度数据和发射率数据。利用MODIS可见光-近红外的第2和19波段进行水汽含量反演，并计算大气透过率。基于亮度温度和大气透过率数据开展劈窗算法温度反演，最后与实测数据及产品数据进行精度对比分析。

研究主要采用覃志豪提出的劈窗算法，该算法以大气辐射传输方程为基础，通过选取MODIS数据两个相邻的热红外通道31和32波段来反演地表/海表温度，算法需要估算比辐射率、大气水汽含量、大气透过率以及MODIS第31和32波段的亮度温度等，该算法地表/海表温度Ts用以下公式表示：

其中：T₃₁和T₃₂为MODIS第31和32波段的亮度温度。A₀、A₁、A₂为系数，由下式求出:

其中:a31、a32、b₃₁和b₃₂为常数，当地表或海表温度在0~50 ℃范围内时，a31=-64.60363，a32=-68.72575， b₃₁=0.440817，b₃₂=0.473453，C₃₁、C₃₂、D₃₁、D₃₂为中间参数，由下式求出：

其中:ε_i为第31波段或第32波段的地表比辐射率，τ_i为第31波段或第32波段的大气透过率。

比辐射率是物体与黑体在同温度和波长下的辐射出射度的比值，它受物体表面状态、介电常数、含水量、温度、物体辐射能的波长等多种因素的影响。水体在热红外波段范围内的比辐射率很高，接近于黑体，根据ASTER提供的常用地表比辐射率光谱库和MODIS的第31和32波段区间的10.780~11.280 μm和11.770~12.270 μm，在MODIS第31和32 波段的比辐射率可分别取ε31=0.996，ε32=0.992。

大气透过率与大气水汽含量之间存在一定函数关系，大气透过率τ及大气辐射强度与大气中的水汽(w)有关，遵循如下辐射传输方程^[19]：

其中:R_s为大气辐射强度，R_g是地面直接反射和散射的辐射强度，ρ是地面反射率。

大气水汽含量估算一般用MODIS的第2和第19波段，估算公式如下：

其中:ω表示大气水汽含量，α和β是常量，分别取值为0.02和0.651。ρ2和ρ19分别表示MODIS第2和第19波段的反射率数据。大气透过率与大气水汽含量之间的关系可以通过MODTRAN软件模拟得到。应用MODTRAN软件还需要TIGR数据，TIGR数据包含全球不同区域的2311条大气廓线，内容包括经纬度、近地面气温、大气气压廓线、气温廓线、水汽含量廓线和臭氧含量廓线等。图2为模拟得到的MODIS第31、32波段的大气透过率与大气水汽含量之间的关系。

由于MODIS传感器边缘区域的视角与星下点视角差别较大，而视角的变化会改变电磁波在大气中的传播路程，从而改变大气的散射、吸收等作用，直接影响到大气透过率的大小。因此，在计算大气透过率参数时，需要考虑传感器的视角校正问题，校正函数如下：

其中:δτ31(θ)、δτ32(θ)分别为第31、32波段的大气透过率视角校正值；θ为MODIS传感器的天顶视角。

MODIS数据是用DN值表示的，因此，要计算星上亮温，必须先将DN值转换成相应的辐射强度值，然后再用Planck函数求解星上亮温，计算公式为:

其中：T_i为MODIS第31和32波段的亮度温度，i=31，32；rad为热辐射强度；scales、offsets为辐射标定系数；K_i,l、K_i,2为常数，分别为1.191 043 56×10^-16W·m²和1.438 768 5×10⁴μm·K。

利用劈窗算法对MODIS数据进行温度反演后，再将研究区部分数据裁剪出来进行制图。研究区包含巴拉望岛岛屿，在制图过程中将该区域剔除，避免巴拉望岛岛屿温度对海洋表面温度分析构成干扰。2013年3月18日研究区温度反演结果如图3所示。

本研究的验证数据为MODIS海洋表面温度产品数据以及部分实测数据。其中MODIS海洋表面温度产品数据为MOD28，是海洋2、3级标准数据产品，内容为海面温度，空间分辨率为1 km，拥有日、旬、月的数据。为了验证劈窗算法对MODIS数据反演温度的精度，选取研究区内42个点反演温度同时期对应的实测海洋表面温度数据和MODIS海洋表面温度产品数据作为验证数据。将实测温度与MODIS产品温度、实测温度与MODIS反演温度作精度验证，结果表明实测温度的平均温度为27.60 ℃，MODIS产品温度的平均温度为26.92 ℃，MODIS反演温度的平均温度为27.69 ℃，3组数据的散点图如图4所示。

为了衡量反演结果的准确性，采用标准决定系数（R²）、和方差（SSE）、均方根误差（RMSE）进行精度评价：

其中：y_i为各散点图的横坐标参数，y_j为各散点图的纵坐标参数，y¯为各散点图横坐标参数的平均值，n为测试的样本数目。对3组数据进行精度评价，如表1所示。

R²即相关系数的平方，相关系数的大小决定了相关的密切程度，当R²越接近1，表示其相关性越好。SSE是拟合数据和原始数据对应点的误差平方和，SSE越接近0，说明拟合更好，数据反演也越成功。RMSE也叫回归系统的均方根误差，RMSE越小，表示反演结果的误差的波动越小，反演结果也越好。从图4和精度评价表可以看出，MODIS产品温度与实测温度拟合效果良好，标准相关系数R²在0.8以上，SSE和RMSE都较小，显示精度良好，在没有实测数据的情况下，MODIS产品温度可以用来作为精度验证数据。实测温度与MODIS反演温度的拟合效果最好，标准相关系数R²在0.9以上，且它们的SSE和RMSE最小，温度误差范围在1 ℃左右，平均温度误差为0.09 ℃，造成反演温度有些误差的原因主要有3个方面：MODIS影像上的云的影响，虽然影像经过去云处理，但是毕竟不能真实还原影像信息，因此，会造成估算大气水汽含量时出现一定的误差，从而导致估算MODIS数据第31、32波段的大气透过率也存在一定的误差；海水的比辐射率受海水盐度、风速、浪高和太阳高度角的影响，会在一定范围浮动，由于在实验过程中将它设定为一固定值，因此也会对温度反演造成一定误差。

2014年MODIS L1B遥感影像云覆盖较少，选取具有代表性的2014年不同时相数据，采用劈窗算法进行多时相温度反演，如图5所示， 从图中可以看出，研究区内温度范围为27.6 ~ 31.3 ℃之间，具有明显的区域和季节变化特征，冬季（12月至2月）温度低，春季（3月至5月）温度逐渐回升，5月份至全年最高，夏季（6月至8月）从最高开始逐渐下降，秋季（9月至11月）降至平稳低温期。其中1月份平均海表温度最低，为27.6 ℃，极大值出现在5月，5月份平均海表温度最高，为31.3 ℃，自5月份以后，月平均海表温度呈缓慢下降趋势，这与内陆地区每年7月或8月温度最高有所不同。研究区月平均温度变化趋势表明利用遥感方法获取的温度与南海海表温度的实际情况一致。

反演得到从2007年至2014年各月的月平均温度，研究海域全年的温度均值范围为27.5~31.5 ℃之间，温度最高的月份为5月份，最低温度一般为1月或2月。为更加清楚地看到2007~2014年的变化趋势，将2007~2014年共96个月的温度采用不记周期连续方法表示，如图6所示，从图中可以明显看出8年相应出现8个波峰，波峰对应于每年5月份，这一变化趋势与研究海域的实际情况是一致的。为了分析研究区海表温度分布特征及变化趋势，对反演后的日温度图进行角投影和数据的格点化，生成反演格点日温度图，基于一个日历月内的30 d左右日温度图进行简单合成和平均，得到月值的月平均温度图，从而制作2014年月平均海表温度等值线图，如图7所示，从图可以看出研究区海域内温度空间分布主体特征，研究区巴拉望岛左侧为南海海域，右侧为苏禄海海域，巴拉望岛附近海域温度略高，在全年都维持一个较高水平。巴拉望岛右侧苏禄海温度低于南海海域，整体上温度变化范围不大。温度从离近岸向中心海域方向递减，海盆中心温度低，这一现象与实际情况相符。

本文以巴拉望岛附近海域为研究区域，主要采用MODIS数据，基于劈窗算法构建研究区海洋温度反演模型，并对反演结果与产品数据及ARGO浮标数据进行对比分析和精度评价。MODIS数据反演标准绝对系数为0.923，均方根误差为0.158。反演精度较高，反演模型具有较好的适用性。本研究结果表明经过精确辐射定标和影像去云等处理的MODIS热红外遥感数据具有很好的海洋表面温度反演能力，简易的劈窗算法有利于反演海表温度。

通过分析研究区海水温度的反演结果，从宏观上揭示了巴拉望岛附近海域温度时空间变化特征，空间上从离近岸向中心海域方向递减，海盆中心温度较低，温度受气候的影响，海表温度的变化与厄尔尼诺/拉尼娜现象存在较强的关联，在出现厄尔尼诺现象的年份研究区相应时间段的海表温度呈偏高趋势，而出现拉尼娜现象的年份研究区相应时间段的海表温度偏低。

当然，本文仅应用MODIS数据进行了海表温度反演，没有开展多源影像数据集成反演算法研究。验证所用的ARGO数据是有限的，今后应该采用更多实测站点的数据对结果进行验证，同时还可以采用其他卫星的海表温度产品进行交叉验证。一般对海洋表面温度的研究倾向于大范围、长时间序列的，而本文应实际要求，只对巴拉望岛附近海域海温变化趋势进行分析，后续可以在前面研究基础上向大范围、长时间序列的海温变化趋势进行研究，获取更多更有价值的信息。