A TRM-based attribution analysis of differences in land surface temperatures in rubber plantations in the Southeast Asian countries

橡胶树（Hevea brasiliensis）又称巴西橡胶树、三叶橡胶树起源于南美洲的亚马逊盆地，现主要产地位于东南亚，其余部分产自南美洲和中美洲，最近西非和中非也开始生产^[1]，是落叶乔木，在每年旱季落叶^[2]。橡胶树一般需要生长5到7年开始生产胶乳，当树龄到30年左右产胶能力开始下降直至无法再生产乳胶^[3]。橡胶林的传统种植区位于赤道附近南北纬10°之间全年潮湿多雨的地区，如泰国、马来西亚和印度尼西亚等；而热带其他地区通常被称为橡胶非传统种植区，如老挝境地区^[4]。由于国际市场对天然橡胶的强劲需求，使得热带地区的非传统种植国家纷纷发展橡胶种植业。特别是上世纪50年代中国在云南西双版纳地区对橡胶林的成功种植和推广，打破了橡胶种植的传统界线，从技术上保障了橡胶林的跨区域种植^[5]。橡胶种植园的快速扩张正持续侵占并破碎化原始森林生境，改变景观连通性与物种栖息地格局，进而触发一系列生态与气候效应^[6]。在生态过程层面，森林转换为单一树种且高强度经营的人工林往往伴随物种多样性衰退、群落结构简化与生态系统功能弱化，包括授粉与天敌网络瓦解、土壤团聚体稳定性下降及微生物群落重组等^[7,8]；在地表能量与水分收支层面，地表反照率、粗糙度与蒸散能力的系统性改变会重塑感热与潜热通量分配，影响局地地表温度（LST）^[9,10]；在碳循环层面，原始森林被替代不仅直接释放生物量与土壤碳，同时降低作为碳汇能力，并可能通过光合能力、呼吸代谢与凋落物输入的改变，影响土壤有机碳形成与周转^[11,12]；在化学通量层面，森林通过交换二氧化碳、挥发性有机物（VOCs）与反应性氮，调节对流层化学与气溶胶生成，从而影响云滴凝结核数目，导致地表粗糙度与冠层湍流条件的改变，又会进一步反馈至臭氧干沉降效率与近地层空气质量^[13-15]。相关研究表明，LST在全球范围内呈上升态势，且橡胶林和其他林有明显温度差异^[16]，这一变化既受到大尺度气候变暖和大气CO₂浓度升高的背景影响，也与区域尺度上的土地利用/土地覆盖变化和植被管理方式密切相关（如造林、退化、农林转换等）。当一种植被类型被另一种类型取代时（例如天然林向人工林），其生物物理属性会发生系统性调整，从而改变地表能量收支并引起温度差异^[17]。有研究表明橡胶种植园相较周边常绿阔叶林常呈系统性升温^[18,19]，气温升高会伴随土壤含水率下降，从而降低土壤团聚体稳定性与胞外聚合物（EPS）生成，抑制微生物碳周转与氮矿化效率，进而增加地表疏松化、溅蚀与片蚀风险，扬尘与粒径选择性流失加剧导致土壤有效养分与有机质下移或输出^[20-22]；而夜间增温与热岛效应会提高可燃物温度并降低平衡含水率，抬升细尺度燃料可燃性与点火概率，若与区域热浪同期发生，极端火险日的频次与持续时间可能非线性上升^[23,24]；林下微气候趋热趋干削弱阴湿生态位，导致地被植物、蕨类、苔藓及节肢动物等小型生物群落的丰度与β多样性下降，传粉与分解等关键生态功能受损^[16,25,26]。经营层面上，白天高温与空气干燥压缩适宜割胶时窗，影响凝固动力学与乳胶干胶率，诱发年际与季内产量与品质波动^[27,28]；为维持树体水分需增加灌溉、喷雾降温与防虫投入，单位产品的能耗上升^[29]；劳动者在高温工况下面临热应激、脱水与肾功能损伤等健康风险上扬，对作业时间安排、个人防护与职业卫生提出更高要求^[30,31]。然而，现有工作多聚焦反照率与蒸散差异，对空气动力学阻力（ra）的定量贡献缺乏月尺度解析，致使橡胶林全年持续增温的物理图像不完整^[9,32]。

尽管多项研究已指向ra在调控地表能量分配与LST对植被变化响应中的核心地位^[17,33,34]，通过改变下垫面粗糙度、位移高度与冠层-边界层耦合强度来影响感热（H）与潜热（LE）的分配^[35,36]，但针对橡胶林和周边其他林这一特定对照体系的机制性量化仍明显不足。此外，LST差异并不能用遥感观测到的表面温度变化曲线做简单替代。评估橡胶林相对其他林的LST差异需要采用更具因果识别力的框架。为此，我们采用一种基于物理的归因方法，温度响应模型和TRM框架相结合^[35,37]，根据地表能量平衡方程中地表长波辐射项的线性化，用于估算土地利用和覆盖变化对地表温度的直接响应。该温度响应模型广泛应用于地表温度研究中，以区分与土地覆盖转化相关的各个生物物理因子对模型预测的温度变化。

东南亚位于亚洲大陆的东南部（92°E–140°E，10°S–28°26′N），涵盖泰国、中国南部、越南、老挝、马来西亚和印度尼西亚等国家与地区。该区域处在亚洲与大洋洲以及印度洋与太平洋的交汇处，既是生物地理学上的过渡带，也是欧亚大陆与大洋洲之间的重要桥梁。其南部与东部为马来群岛，横贯太平洋与印度洋的广阔海域，由两万余座岛屿组成。该地区深受季风影响，表现为大尺度的季节性风向反转，通常在冬季降水较多。东南亚亦是全球重要的生物多样性热点区。其卓越的物种丰富度以及复杂而动态的地形、岛屿体系与海洋环境，反映出多样的地质构造，这些因素共同塑造了区域的生物多样性演化与生物地理分布格局^[38,39]。

本研究所使用的橡胶林样本点（图1）综合整合了多源信息，主要包括三类数据来源：（1）海南岛实地调查获取的高可信度样本点；（2）基于区域已有的橡胶林分布数据中筛选出的可靠橡胶林像元^[1,40,41]，以及（3）通过高分辨率遥感影像进行目视解译所获得的橡胶林样本。其他林（非橡胶林）样本的构建则由两部分组成：一部分来自海南岛实地调查的非橡胶林样本点；另一部分来源于LUCU数据集中森林类别所对应的非橡胶林像元。此外，为确保与橡胶林样本的空间环境保持一致性，我们进一步从LUCU数据中识别出与橡胶林像素相邻的森林地块，将其作为周边森林的代表性样本，以增强非橡胶林样本在结构、物候和环境条件上的可比性。主要使用了MODIS数据包括：8天平均的1 km昼/夜地表温度数据（MOD11A2）^[42] 和下行短波（DSR）数据MCD18A1^[43]、每日500 m的反照率数据（VNP43DNBA3）^[44]、8天的500 m潜热数据（MOD16A2）^[45]；除此之外还有月度的感热数据（GLEAM4.2a）^[46]，每日1 km的下行长波（SLDR）数据（ELITE）^[47]。由于部分地区云影或者其他因素影响出现的空值，这部分数据均被剔除。其中，反照率数据同时包含短波（0.3～5.0 µm）宽谱段的定向半球反射率（黑天反照率）与双向半球反射率（白天反照率），为了更接近实际条件的情况下，我们取两者的平均进行分析^[44,48]。其中，所有遥感数据均被重采样至1 km空间分辨率，用于统计分析和成对比较。此外，本研究还采用了来自全球气候变化研究数据集（ISIMIP）的数据2m温度数据^[49]，并通过温度–植被指数方法（TVX）将分辨率降为1km^[50]。样本点则是东南亚范围内共计6510个橡胶，及其周围的15441个其他林，栅格大小为1 km。

Figure 1.  Distribution of sampling points

为保证后续1km分辨率下橡胶林气候效应归因分析的尺度一致性，本研究首先对橡胶林分布数据进行空间聚合处理。具体而言，我们以1km×1km网格作为目标尺度，在每个网格中分别统计森林与橡胶林的像元比例。为确保目标网格在生态属性上的主导性和代表性，首先筛选森林覆盖比例超过50%的网格作为森林主导地块；随后在这些地块中进一步筛选橡胶林面积占比超过50%的网格，识别为1km尺度的橡胶林主导区域。且为避免空间分辨率升尺度至1km过程中产生统计偏差，本研究进一步考虑了由非线性函数形式引起的尺度效应，采用泰勒二阶展开来矫正尺度效应（公式如下），使1km网格在保持主导植被类型一致性的同时，最大程度减少非线性聚合所导致的估计偏差。通过上述逐级面积阈值筛选策略，本研究有效减弱了由空间分辨率差异引起的尺度效应，确保升尺度后的橡胶林数据在空间一致性与生态意义上与原始高分辨率数据保持匹配，从而为后续气候效应归因分析提供可靠的基础数据。（代码开源：升尺度与二阶校正算法作为补充材料上传。）

其中，μ是NDVI的均值，σ是NDVI的标准差。

温度响应模型基于地表能量平衡方程中对地表长波辐射项的线性化处理^[51]，用于估算土地利用/覆盖变化（LUCC）引起的LST响应，已被广泛用于LST相关研究，可以分离相互关联的各个生物物理因子对模型预测的ΔT的贡献。为使计算得到的ΔTs更接近观测值，我们采用了经修改后的温度响应模型^[35]。相较于传统能量平衡方法，该模型不仅能够保持明确的物理意义，还具有可分解性强、参数需求低、适用于大尺度遥感研究等优势。为进一步揭示非辐射过程在地表温度变化中的作用机制，本研究在温度响应模型的基础上融入了部分扩展的TRM框架^[52]，用以构建同时量化辐射与非辐射贡献的综合归因体系。以地表能量平衡方程为基础，将净辐射项分解为感热通量、潜热通量以及地热通量（式3）。其中，感热通量进一步结合空气动力学阻力（ra）（式4），从而使非辐射过程直接纳入温度变化的数学表达。又对地表温度的控制方程进行一阶泰勒展开，得到地表温度对各生物物理因子的解析敏感度，包括反照率、短波辐射、向下长波辐、空气阻力、空气温度与地热通量变化的偏导项（式5），反映了LST对不同能量分量及大气背景条件变化的响应强度。

具体公式如下：

其中，$ {\lambda }_{0} $：月平均Planck响应，计算方式为1/4$ {\varepsilon }_{s}\delta T_{s}^{3} $，单位为K·m²·W⁻¹；$ {\varepsilon }_{s} $是地表的月平均发射率，近似为常数0.97；$ \Delta R_{n}^{*} $：净辐射变化，单位为W·m⁻²；$ \Delta G $：土壤热通量变化，单位为W·m⁻²；$ f $：能量再分配因子；$ {T}_{S} $：地表温度，单位为°C；$ {T}_{a} $：空气温度，单位为°C；$ \alpha $：反照率；$ S $：短波辐射，单位为W·m⁻²；$ {L}_{\downarrow } $：向下长波辐射，单位为W·m⁻²；$ \delta $ ：Stefen-Boltzmann常数，为5.670374419×10⁻⁸W/m²K⁻⁴；ρ表示空气密度，约为1.2 kg/m³；$ {C}_{P} $为空气在定压条件下的比热容，1004J/（kg∙K⁻¹）；$ {r}_{a} $：空气动力阻力，单位为s·m⁻¹；H：显热通量，单位为W·m⁻²；LE：潜热通量，单位为W·m⁻²。

其中，Δ表示变化量（即橡胶林和周边其他林差异）。偏导数表示LST对各地表生物物理因子变化的敏感性，这些敏感性的解析形式可由$ {T}_{S} $的表达式推导得到^[52]。

在选取满足“森林占地块≥50%且橡胶林占森林≥50%”的1 km网格后，对样本进行了尺度效应校正，最终共获得7879个符合条件的1km橡胶林样本点（图1）。采用这些样本点对高分辨率NDVI与1kmNDVI进行一致性评估，结果显示：皮尔逊相关系数达到0.9023，均方误差（MSE）为0.0032，平均绝对误差（MAE）为0.0484，结构相似性（SSIM）为1；T检验的结果显示：t=−0.9518、p=0.341，说明了两组数据的差异相对于样本波动很小，即升尺度过程没有显著偏高或偏低；进行纬度与坡向分组统计的结果也显示，不同纬度区间内的误差均值保持在±0.001以内，且未出现随纬度升高而增大的趋势，在坡向分类（N/E/S/W）下的误差均值介于−1.98×10⁻⁵至4.63×10⁻⁵之间，说明升尺度并未导致空间结构性偏移，可用于后续分析。整体而言，两种尺度NDVI之间表现出极高的一致性。基于上述统计结果，可判断所选取的7879个样本点能够有效代表橡胶林，其植被特征受橡胶林覆盖主导。因此，本研究后续的相关性分析均基于这7879个橡胶林样本点开展。

橡胶林与其他林地之间的地表温度差（ΔT）呈现一定的季节性变化特征（图2）。全年ΔT在2～6月维持较高水平，并在5月和6月达到最大值（接近1.2°C）。其中，2～4月ΔT持续上升，3月和4月分别达到0.87～1.03°C。进入7～9月后，ΔT显著下降并出现全年最低值，最低值出现在7月（约0.07°C）。随后ΔT自10月起再次升高，在11～12月恢复至0.67～0.73°C的水平。ΔT在不同纬度呈现一定的空间异质性：低纬区（0°～5°N）ΔT波动较小，主要集中在−1℃至1℃之间；中低纬区（5°～10°N）ΔT的变幅明显增大，局地可达2℃以上；中纬区（10°～15°N）则出现负向偏移，有ΔT落在−2℃以下；而在更高纬区（15°～20°N），ΔT再次出现正向上升并伴随高频波动，峰值可达3～4℃。总体而言，ΔT的全年分布表现为：春末至初夏偏高、雨季最低、秋季后再度升高，并呈现出明显的月际梯度变化且在空间上呈现出明显的起伏。

Figure 2.  Monthly temperature differences between rubber plantations and other forests across regions in 2021 (observations).

为深化对橡胶林与周边其他林地表温度（LST）差异的物理机理认知，ΔLST进行了逐月分解，量化了Δrₐ、ΔL_↓、ΔS与Δα等路径对温度差异的边际贡献。结果显示：

（1）Δrₐ主导升温，其正向信号，且抵消了ΔL_↓的负向信号；

（2）ΔL_↓对ΔLST呈降温作用，但其逐月变化与ΔLST的相关性并不显著；

（3）两类贡献在11～2月差异尤为突出，体现为rₐ路径增温与L_↓路径降温共同塑造，但净效应仍由rₐ所决定。

以1月为例，路径分解给出的边际贡献为：rₐ、α 与S分别带来约1.20℃、0.03℃ 与0.07℃的升温；L_↓贡献约0.90℃，体现显著的降温效应。若将各项求和，可得到与模型ΔTs近似但非完全一致的值，我们对贡献极小项（如ΔG、$ \Delta {T}_{a} $均小于0.01℃）做了阈值截断以及线性化近似与栅格尺度下时空异质性共同引入的小幅残差。因此，模型分解的ΔTs与MODIS反演ΔTs在数值上存在一定偏差，但趋势与量级保持一致（表1），支持上述机理判断的稳健性。

全年ΔT显示出一定的季节差异，反照率项贡献始终较小（多为0.02～0.04°C），对总体变化影响有限；短波辐射项在6月达到峰值（0.47°C），在7～9月转为负值，呈现季节反转。长波辐射项全年均表现为冷却效应，典型值为–0.30°C，并在1月与11月达到显著低值（–0.90与–0.92°C）。相比之下，空气动力阻力项在所有月份均为正值，贡献幅度最大（多数月份≥1.0°C），为全年ΔTs的主导因子。整体来看，增温主要由r_a驱动，S在光照强的月份提供附加强化，而L_↓和α在多数月份表现为次要或抵消作用。

为系统阐明橡胶林与周边其他林之间LST差异的物理成因，在温度响应模型及TRM框架下，对全域样本的年度尺度数据（图3）和月尺度数据（图4）分别进行了归因分析，将ΔLST细分为若干可观测、可解释的物理路径，包括α、r_a、L_↓与S等分量的贡献。在夜间条件下，r_a的贡献为56.84%，在白天增加至66.82%，全天平均仍维持在64.39%；L_↓贡献范围介于18.75%～24.02%，全天为22.73%；S夜间24.02%、白天18.75%、全天10.61%；α贡献最弱（<3%），在各时间尺度均处于边缘化水平，对ΔT的影响相对有限。

Figure 3.  Proportion of the main factors influencing LST

Figure 4.  Monthly relative contributions of individual factors to ΔTs

结果表明（图4），r_a在全年均保持最高的贡献水平，多数月份占比超过50%，在3月、4月、10月与12月达到60～95%的主导地位，相比之下，L_↓的贡献呈中等水平，在1月、7月、8月、11月与12月等湿度相对偏低或辐射条件发生转折的月份贡献占比可达30～40%，体现其在部分月份具有次要但不可忽视的影响。短波辐射S的贡献存在季节差异，在5～6月达到全年峰值（30～38%），与雨季来临前后地表净辐射变化一致；其余月份贡献较低（5～15%）。反照率α的贡献始终保持在较低水平（约1～6%），在全年温度差异解释中作用相对有限。总体而言，月尺度ΔTs的贡献结构呈现出高度一致的主导因子格局：r_a为全年第一主导因子，L_↓为次要贡献因子，而S与α的贡献相对较弱但具有明显的时序波动性。

在像元尺度的归因分析中（图5），各主导生物物理因子对ΔLST的贡献呈现出极为不均衡的分布格局。统计结果显示，绝大多数像元的ΔLST差异主要由r_a所主导；相比之下，S仅在3个像元中表现为主导因子，α仅在7个像元中占据主导地位，L_↓则在27个像元中成为主要贡献者。整体而言，r_a不仅是最普遍的主导因子，也是解释区域尺度ΔLST的关键生物物理机制。

Figure 5.  distribution of pixels dominated by different attribution factors

在将土壤湿度划分为干旱（<20%）、正常（20～35%）和潮湿（>35%）后可以发现（图6），在土壤干旱条件下，日间ΔLST为0.25°C，而夜间转为−0.33°C；在正常水分条件下，日间ΔLST升高至1.44°C，夜间为0.80°C；在潮湿条件下，日间ΔLST为0.73°C，夜间为−0.14°C。总体来看，在干旱和正常水分区域，夜间ΔLST相较日间普遍降低约0.5～0.6°C，而在潮湿区域，这一昼夜差幅扩大至0.8～0.9°C，橡胶林增温效应在日间更为明显，在夜间则出现一定的差异化，其原因可能在于冠层结构和空间异质性差异在夜间在湍流减弱和辐射冷却的条件下，使得橡胶林相对于周边林地更易出现温度差异的减弱甚至反向变化。同时，在干旱情境下，r_a的贡献占比为59.8%，为主导性因子，表明此阶段地表能量过程主要受动力限制控制；当土壤湿度处于正常范围时，r_a的贡献进一步上升至73.9%，动力过程的主控地位得到强化；然而，在进入潮湿阶段后，主导因子由r_a转向长波辐射L_↓，其贡献升至78.4%，在这一阶段，r_a的贡献降虽至12.7%，表现为次级调节项，但并未完全消失，而是由核心驱动退化为背景调控机制。

Figure 6.  Dominant contributors to ΔLST across soil moisture regimes (from inner to outer: dry, normal, wet)

本研究通过对尺度一致化后的ΔLST逐月分解与物理路径归因，揭示了橡胶林与周边其他林之间地表温度差异主导机制和不同季节的差异性。整体上，橡胶林全年LST均高于其他林，且差异在旱季末至雨季初最为显著，而在雨季中期显著减弱。该季节性差异不仅与区域水热条件调控的能量分配相关，也反映出橡胶林独特的物候节律对地表能量平衡的影响。首先，从能量分配角度来看，旱季末与雨季初往往伴随较高的辐射输入，但土壤水分与冠层导度却未能同步恢复至雨季水平，使得LE受限，能量优先分配给H，从而推动LST升高^[53]。这一机制可能与橡胶林干季落叶—返叶的物候节律导致的叶面积指数降低、冠层未完全封闭、气孔调节能力下降等因素相关，使得橡胶林在此阶段更易向感热主导的能量分配模式偏移^[54]。另一方面，TRM与温度响应模型的路径分解结果表明，r_a是全年ΔLST的首要贡献因子，在多数月份提供主导性的升温信号。橡胶林具有较大的粗糙度长度和位移高度，这种冠层形态在结构上赋予其更强的湍流交换潜力，使其在相同能量输入条件下更易发生动量与热量的垂直输送^[55]。然而，该潜力并未完全转化为有效的蒸散调节：由于LAI配置、冠层通风效率及气孔导度受水力条件限制，导致了在rₐ偏高/湍流传输受限的情况下，能量更倾向于以H向外传导，进而推升地表温度^[56]。相比之下，辐射相关项对ΔLST的贡献虽然方向明确，表明橡胶林与其他林之间的辐射差异更多受云量、降水与区域性大气背景调控。因此，辐射项在解释局地ΔLST差异时呈现一定的被动响应特征，而非主驱动因子^[57]。ΔLST的形成并非由单一物理过程决定，当水分受限时，蒸散减弱和气动阻力升高驱动感热通量显著增强，使地表温度对ra敏感性变高^[58]；而在湿润环境下，水汽的累积提升大气光学厚度，下行长波辐射反馈加强，使温度变化由辐射主导^[59]。表明ΔLST调控机制具有水分依赖性，在干旱区可优先关注降低r_a的管理策略，而湿润区的调控策略则更依赖辐射收支平衡。

同时，雨季出现温差最小的现象也进一步印证了这一点：在云量增加、净辐射被抑制的雨季背景下，植被物候差异的影响被显著削弱，而辐射项的调节作用加强，从而压低了ΔLST的季节波动^[60]。值得关注的是，模型分解得到的逐月边际贡献项与实际ΔLST之间存在轻微偏差，主要来源于阈值截断、线性化与栅格异质性等因素。然而，从趋势与量级上看，两者保持一致，说明TRM路径分解结果具有较高的物理合理性与统计稳健性。尤其是ra路径的主导效应，与热力边界层理论、森林种植园湍流结构差异的经典认识高度一致，进一步强化了“空气动力学结构差异”而非“辐射差异”是橡胶林更热的核心原因。在泰国-老挝边境已划定的橡胶限制区内，推广“宽行密株+保留次生带”的林分结构优化策略，通过改善林冠层通风效率和增加地表粗糙度来降低r_a。根据本研究的情景推算，ra可下降约15%，并在实现约0.3℃的降温效应。同时，将r_a纳入东盟橡胶碳汇项目（AS-REDD+）的协同效益评价体系，并与区域碳价机制挂钩，可实现碳汇增益与微气候调节的双重收益。这些结果不仅深化了对橡胶林与其他林能量交换过程的认识，也为解释热带种植园的气候效应提供了新的证据，对区域气候评估、生态管理与可持续种植政策具有重要参考价值。

本研究以温度响应模型结合TRM框架，将橡胶林和其他林的ΔLST分解为辐射项与非辐射项两类可观测路径。结果显示，2021年东南亚橡胶林全年较周边其他林更热（年均+0.63℃），6月差值最大（1.16℃），7月最小（0.07℃），且体现出季节差异：旱季（2～4月）与早雨季（5～6月）“更热”信号更强，而雨季盛期（7～9月）有所收敛。归因路径上表明ra是解释ΔLST的首要因子，其贡献占比达64.4%，显著高于α、S与L_↓等辐射相关项。并逐月分解进一步揭示了季节差异背后的物理基础：在11～2月，Δrₐ的正向（增温）更高，能够抵消ΔL_↓的系统性降温效应，解释了干冷季节ΔLST仍以正值为主的现象。但在干旱与正常水分条件下，温度调控以ra为主导，而在湿润条件下，L_↓的影响会更显著。值得注意的是，TRM中以ra表征的非辐射效应，本质上是ra通过影响感热与潜热通量、鲍文比值以及总传输阻力等多路径的耦合影响导致了地表增温。