利用无人机多光谱影像和环境数据，通过符号回归方法预测酸樱桃茎的水分势 Anderson L. S. Safre, Alfonso Torres-Rua, Kurt Wedegaertner, Brent Black, Brennan Bean, Burdette Barker, Matt Yost

《Remote Sensing》：Predicting Tart Cherry Stem Water Potential Using UAV Multispectral Imagery and Environmental Data via Symbolic Regression Anderson L. S. Safre, Alfonso Torres-Rua, Kurt Wedegaertner, Brent Black, Brennan Bean, Burdette Barker and Matt Yost

　　**亮点**

**主要发现是什么？**
通过结合无人机多光谱图像、气象数据和土壤湿度数据，可以使用符号回归方法准确估算酸樱桃果园的茎水势（Ψstem），其R2值可达到0.80，均方根误差（RMSE）低至0.08 MPa。从RGB波段衍生出的红色色坐标（RCC）指数是预测茎水势的最有信息量的光谱指标，尤其是在采收前期间。

**主要发现的意义是什么？**
所提出的符号回归模型提供了可解释且可转移的方程式，克服了常见机器学习方法在估算茎水势时的“黑箱”局限性。在这项研究中，符号回归能够生成对可用信息源作出响应的可解释数值表达式，同时最大化选定输入（航空、气象、土壤湿度）对描述茎水势的贡献。

**摘要**
酸樱桃是犹他州的重要水果作物，由于干旱条件，灌溉至关重要。精准灌溉需要可靠的植物水分状态指标，而茎水势（Ψstem）是最敏感的指标之一，尽管其测量过程较为繁琐且受空间限制。本研究利用高分辨率无人机（UAV）多光谱图像以及气象和土壤湿度数据，应用符号回归（SR）方法开发了茎水势估算模型。结果显示，在采收前期间，光学波段与茎水势之间存在更强的相关性。在85个植被指数中，红色色坐标（RCC）指数的表现最佳（R2 = 0.67）。针对不同的数据可用性情景，生成了六个方程式，并通过留一法（修改后的k折法）进行了验证，得到的茎水势估算R2值范围为0.67至0.80，均方根误差（RMSE）范围为0.11至0.08 MPa。值得注意的是，符号回归能够生成提高模型透明度和可转移性的可解释方程式。通过使用来自不同地点的独立数据集进一步确认了模型的稳健性。据我们所知，这是首次将符号回归应用于茎水势估算，为酸樱桃果园的灌溉管理提供了一个可扩展且可解释的工具。

**1. 引言**
酸樱桃是犹他州的重要特色作物。2022年，酸樱桃种植面积达1503公顷，其中1230公顷处于结果期，产量为10,200毫克[1]。樱桃产业管理局[2]报告称2024年的产量为19,800毫克，占全国总量的16.6%。大多数果园位于犹他县，这是一个温带半干旱地区，灌溉对于实现高产至关重要。随着人口增长和耕地有限，种植者需要采用既能保证环境可持续性又能保证经济可持续性的灌溉方式。在这种情况下，灌溉者必须就灌溉系统、策略和调度方法做出有效决策[3]。水资源日益紧张，进一步凸显了在酸樱桃生产中精确高效灌溉调度的必要性。
美国西部的灌溉管理因作物类型和灌溉方法而异[4]。大多数酸樱桃果园使用安装在固定横向管线上的微喷头，其中许多是几十年前设计的。在犹他州北部，有研究指出[5]，小规模农场经常过度灌溉，导致深层渗透损失，尤其是在生长季节初期。这些低效率问题受到灌溉水分配方式的影响[6]。有研究指出[6]，犹他县的灌溉调度通常依赖于“轮流供水”制度，即根据水资源的可用性而非作物需水量进行灌溉。虽然已经建立了喀斯喀特山脉以东甜樱桃的作物系数[7]，但酸樱桃尚未有相应的估算值。在夏季降水量低且蒸气压亏缺（VPD）高的地区，甜樱桃树在生长季节可能需要750至1000毫米的灌溉量[8]，这表明类似气候下的酸樱桃果园可能有类似的需求。

**2. 材料与方法**
**2.1. 实验地点**
实验在犹他州农业实验站Kaysville研究农场的一块1.21公顷的‘Montmorency’酸樱桃（Prunus cerasus L.）果园中进行（坐标41.0202N, ?111.9301W，海拔1314米）。该果园建于2010年，树间距为3.6米×6米。根据美国农业部土壤调查，该地点的土壤为Kidman土壤，质地为细砂壤土，排水良好，渗透性适中。该地区属于半干旱气候，夏季炎热干燥，年平均降水量为558.3毫米。5月至10月期间，参考蒸发量（ET0）通常超过降水量。灌溉采用Nelson R10喷头（Nelson Irrigation Corporation，华盛顿州Walla Walla）进行的微喷灌系统，喷头间距为8米×6米，每次灌溉持续10小时，灌溉量为约25.4毫米。气象变量（包括每小时气温（Ta）、相对湿度（RH）、太阳辐射、ET0和降水量）由位于果园200米内的犹他气候中心气象站记录（https://climate.usu.edu/），数据获取时间为2026年3月3日）。蒸气压亏缺（VPD）根据Ta和RH计算[38]。土壤含水量通过七个传感器节点监测，每个节点配备四个Teros 12传感器，连接到ZL6数据记录仪（METER Group Inc., 华盛顿州Pullman），安装深度分别为15厘米、45厘米、60厘米和76厘米。

**2.2. Ψstem采样和测量协议**
选择了14棵树进行Ψstem采样。Kaysville地点的采样树木位置、土壤湿度传感器和气象站的位置如图1所示。图1显示了实验酸樱桃果园的位置和研究区域布局。感兴趣的区域（AOI）用虚线标出，采样树木用黑色矩形表示。蓝色点表示土壤湿度传感器的位置，气象站用黄色星号标示。定位图中的红点表示犹他州的研究地点。使用配备数字压力表的Scholander压力室（Model 1505D，PMS Instruments，俄勒冈州Albany）进行Ψstem测量。在每棵采样树的树冠内选择三片靠近树干或主支架枝条的叶片。叶片在测量前至少用层压箔袋包裹一小时。在叶柄处进行干净切割后，立即将叶片放入压力室。通入氮气直至水分通过木质部排出（终点），并记录压力值。测量时间选在接近正午（12:00至14:00）以捕捉最大水分胁迫条件。共进行了24次Ψstem数据采集活动。2023年，6月18日至9月9日期间进行了六次Ψstem测量，旨在采集采收前后的数据；2024年数据采集活动增加到19次，从5月15日（叶片展开后三天）持续到9月24日，以捕捉更细致的Ψstem动态变化。

**2.3. 土壤表观电导率测绘**
测量土壤表观电导率（ECa）以表征果园内的土壤变异性。使用EM38探针（Geonics，加拿大Mississauga）进行测量，该探针工作频率为14.5 kHz，有效感应深度约为1.5米。该探针提供四相（电导率）和同相（磁化率）读数，单位为mS/m，电导率噪声水平为0.5 mS/m。EM38探针被安装在一个木制雪橇上，并使用全地形车（ATV）以11-13公里/小时的速度拖曳穿过田地。Trimble Ag-114 GNSS接收器（Trimble Inc.，美国科罗拉多州威斯敏斯特）以1秒的间隔记录位置数据[39]。场地布置如图2所示。ECa点测量数据使用ArcGIS Pro v.2.9.5（Esri，美国加利福尼亚州雷德兰兹）中的普通克里金法进行插值，拟合过程中优化了半变异图模型参数。图2. 安装在木制雪橇上的Geonics EM38探针，配备有GNSS接收器，用于土壤ECa测量。

2.4. 无人机多光谱和热成像数据采集与校准
使用配备Altum-PT多光谱和热成像相机的DJI Matrice 300 RTK（DJI Technology，中国深圳）收集高分辨率多光谱和热成像数据（Altum-PT相机由AgEagle Aerial Systems，美国堪萨斯州威奇托制造）。共进行了25次飞行。飞行遵循基于2D网格的自动任务规划，前后重叠率为85%，飞行速度为6米/秒，飞行高度为120米，地面采样距离为2.5厘米。所有飞行都在太阳正午前一小时内完成，以最小化阴影效应。Altum-PT相机在五个光谱波段捕捉图像：蓝光（475纳米，带宽32纳米）、绿光（560纳米，带宽27纳米）、红光（668纳米，带宽14纳米）、红边（717纳米，带宽12纳米）和近红外（NIR）（842纳米，带宽57纳米）。1220万像素的全色传感器（634纳米，带宽463纳米）用于图像锐化，将空间分辨率提高到每像素2.5厘米。热成像数据通过Altum-PT相机中的集成长波红外（LWIR）传感器（11微米，带宽6微米）获取，地面采样距离为30厘米。图3显示了配备Altum-PT载荷的无人机平台、Aeropoints地面控制点（GCPs）和热校准仪器。

2.5. 正射影像生成处理
使用Agisoft Metashape专业版v.2.0.3（Agisoft LLC，俄罗斯圣彼得堡）生成正射影像。根据制造商的说明进行了图像锐化和辐射度转换。每个光谱波段的反射率（ρ）使用公式（1）计算：
$$
\rho_{\lambda} = \frac{D_N\lambda(2D_N634 + 0.2D_N560 + 0.2D_N668 + 0.2D_N717 + 0.2D_N942)}{32768}
$$
其中$D_N\lambda$是特定波段的数字编号，$D_N634$是全色数字编号，32768对应于16位图像的范围。热数字编号通过以下公式转换为温度：
$$
L_WI_R(\degree C) = D_NL_WI_R \times 100 - 273.15
$$
为了校正LWIR温度估计值，将每个校准篷布中心的平均像素值与红外辐射计（IRT）测量值进行对比。得到的校准方程为：
$$
L_WI_R c_o t_t e_c d = 0.5974L_WI_R \times 100 + 9.9875
$$
使用所有25次飞行的数据，通过线性回归得到斜率和系数的拟合优度，如图4所示。图4. Altum-PT热成像与地面红外辐射计（IRT）测量值的校准。散点代表配对的温度观测值，实线红色曲线表示线性回归拟合结果。

2.6. 树冠分割、光谱信息提取和植被指数（VI）计算
开发了一个两阶段分割工作流程，从每个正射影像中分离出纯树冠像素。使用Metashape的自动地面点检测生成并分类了密集点云。参数包括单元格大小为50米、最大角度阈值为20°和最大距离阈值为1米。从未分类的点生成数字表面模型（DSM），而从已分类的点生成数字地形模型（DTM）。然后计算树冠高度模型（CHM），作为DSM和DTM之间的差值。保留CHM值≥1米的像素以生成第一个二值掩膜。
由于酸樱桃树呈开放灌木结构，中心树冠稀疏，因此应用了第二步细化步骤，使用二值NDVI掩膜去除树冠内的非叶组织和辐射极端像素。检查了树级NDVI的统计信息（最小值、平均值和最大值）以表征NDVI分布。平均NDVI值紧密聚集在主导树冠反射率范围（约0.88）附近，而较低和较高的极端值与异常值相关（附录B）。根据这些树级统计的直方图，将NDVI值≤0.78和≥0.92识别为分布尾部，并用于定义应用于每个光谱波段的二值掩膜的边界。图5显示了树冠分割过程的工作流程和结果。图5. 上半部分：从树冠像素获取平均反射率值的步骤流程图。下半部分：两阶段分割过程的结果。使用矩形网格提取了每棵样本树的平均反射率值。每个多边形内的像素值直方图显示平均值和中值之间没有显著差异（附录C），验证了使用平均值的有效性。使用Python 3.12.7脚本自动化了每次飞行的提取过程。根据“Awesome Spectral Indices”目录（https://github.com/awesome-spectral-indices/awesome-spectral-indices；2024年6月2日访问）从多光谱影像中计算了85个植被指数。指数按照决定系数（R2）和Ψstem进行排序。在模型训练之前使用完整数据集计算了这些关系。尽管这种方法可能存在信息泄露的风险，但随后使用来自不同果园的独立验证数据集评估了模型性能。
计算了前20个指数的相关矩阵，并移除了彼此高度相关的指数（R2 > 0.95），选择了前四个独立的指数作为模型预测因子。使用校正后的LWIR影像计算了树冠-空气温差（Tc ? Ta）。

2.7. Ψstem模型构建和性能评估
使用PySR库[41]开发了符号回归（SR）模型来估计Ψstem。pySR通过Python接口将遗传编程原理应用于基于Julia的优化引擎，以识别简单、人类可解释的数学表达式。算法初始化多个随机生成的候选方程种群，并通过变异和交叉操作迭代演化它们。每个进化周期后，使用损失函数评估候选方程，并根据损失和复杂性保留候选方程。PySR框架的更多实现细节可在项目仓库中找到（https://github.com/MilesCranmer/PySR；2025年9月4日访问）。使用普通最小二乘回归（OLS）标准误差评估了最终SR导出方程中系数的统计显著性。
配置了八个种群，每个种群包含50个候选方程。每次迭代进行500个周期，总共进行了40,000,000次进化迭代。操作符集包括加法、减法、乘法、除法、幂函数、平方根和指数函数。
最初根据R2和统计显著性（p值）筛选预测变量。数据集分为采前（n = 141）和采后（n = 141）两部分。在每个子集中，观测值随机分为训练集（75%）和验证集（25%）。训练集用于识别符号回归模型结构并估计系数，验证集用于内部性能评估。通过为变量数量（每个2分）、算术运算（1分）和指数项（2分）分配加权点来独立量化模型复杂性。
由于Ψstem观测值是从多棵树在多个采样日期重复收集的，因此实施了留一树法（LOTO）交叉验证程序来评估SR方程的稳健性和泛化能力。在这个过程中，每次迭代都会排除一棵树的所有观测值。每个SR导出方程的功能形式保持不变，仅使用剩余树木的数据重新拟合系数。然后为被排除的树生成预测值。这个过程重复进行，直到每棵树都作为验证单元使用过一次。使用LOTO RMSE值的配对t检验评估了方程之间的预测性能差异。
除了SR导出的方程外，还拟合了一个广义加性模型（GAM）作为灵活的非参数基准。GAM包含了所有模型中使用的预测因子，并使用R中的mgcv包实现。自动估计了每个预测因子的平滑函数，允许在没有手动指定交互项的情况下出现非线性关系。还包括一个基线比较模型，用于预测训练树的平均Ψstem。

3. 结果
3.1. 2023年和2024年生长季节的Ψstem
图6展示了2023年和2024年的Ψstem箱形图和直方图。2023年，Ψstem的范围为-0.40至-1.23 MPa，季节平均值为-0.89 MPa。2024年，值的范围为-0.41至-1.25 MPa，平均值为-0.82 MPa。两个生长季节中，Ψstem都遵循一致的季节性模式，即季节初期值较正，随着树冠发育逐渐变为更负的值。2024年的更频繁测量证实了2023年观察到的模式，包括在生长日（DOY）175附近出现的明显变化。2023年的收获发生在DOY 206-208，2024年发生在DOY 199-201。2024年的采后Ψstem测量显示出比2023年更多的短期波动。图6. 2023年和2024年生长季节Ψstem测量的分布。箱形图和直方图总结了所有采样日期的Ψstem的季节范围和变异性。近似收获日期由垂直红色虚线表示。圆圈代表异常值。图7比较了2023年和2024年相近生长日对的Ψstem值。季节初期（DOY 169-170），Ψstem较不负；而在季节中期，Ψstem变得更负，对应于果核硬化的开始和进展。后期测量没有显示出明显的单调下降；相反，Ψstem值保持在相似范围内或略低于季节中期的值。尽管两年间的总体季节模式相似，但某些生长日对的样本树之间的变异性较大，2024年的中位数略低。图7. 两个生长季节相似日期的Ψstem测量箱形图。x轴代表2023年（左侧）和2024年（右侧）的测量日期。圆圈代表异常值。

3.2. 土壤ECa和土壤含水量与Ψstem的关系
使用EM38调查数据生成的土壤表观电导率（ECa）地图如图8所示。果园内的空间变异有限，电导率（ECa）值介于14至28 mS/m之间。在监测的树木子集中，变异幅度更小（20–25 mS/m），表明该区域的土壤剖面相对均匀。图8显示了使用EM38探针测量的土壤表观电导率（ECa）的空间分布。果园的东部和北部区域观察到较高的ECa值，尽管整个研究区域的整体变异较低。图9展示了测量期间所有传感器深度的日土壤含水量（SWC）。2023年，整个季节的SWC保持相对稳定。在15厘米深度，收获前的SWC值稳定在0.20 m3/m3左右，收获后增加到0.25 m3/m3。在45厘米深度，收获前的SWC在0.21至0.22 m3/m3之间波动，收获后增加到0.25 m3/m3。更深的传感器显示出类似的趋势：60厘米深度的SWC在收获前从0.23降至0.21 m3/m3，然后再次增加到0.23 m3/m3；而在75厘米深度，SWC从0.25降至0.22 m3/m3，随后又增加到0.23 m3/m3。图9显示了整个生长季节中四个深度（15、45、60和75厘米）的日土壤含水量（SWC），数据来自七个传感器节点。灰色垂直虚线表示无人机图像采集日期。相比之下，2024年的SWC变异较大。在15厘米深度，SWC从0.23降至0.19 m3/m3，收获后稳定在0.20 m3/m3。在45厘米深度，收获前的SWC波动范围更广（0.27至0.19 m3/m3），收获后稳定在0.20 m3/m3。60和75厘米深度的SWC分别呈现季节性下降趋势，从0.26降至0.19 m3/m3和从0.27降至0.20 m3/m3。图10使用所有单个传感器测量数据展示了SWC与Ψstem之间的关系。收获前的相关性更高，15、45、60和75厘米深度的R2值分别为0.30、0.31、0.14和0.03。平均后，R2值为0.31。图10的散点图显示了收获前后四个深度（15、45、60和75厘米）的Ψstem与土壤含水量（SWC）之间的关系。每个点代表七个传感器节点中的一个单独测量值。实线回归线表示收获前的关系，虚线表示收获后的关系，点划线表示使用所有数据进行的整体拟合。相应的R2值显示在每个面板中。

2023年，累计降水量为107毫米，其中8月份降水量最高（43毫米），而累计蒸散量（ET0）达到777毫米，7月份达到峰值（196毫米）。2024年，累计降水量为177毫米，5月份的降水量最高（57毫米）。整个季节的累计蒸散量为811毫米，7月份的蒸散量再次达到最大值（193毫米）。2023年，土壤温度（Ta）的范围为4°C（第125天）至39°C（第203天），水分势差（VPD）的范围为0.03至6.04 kPa，最大值也出现在第203天。太阳辐射范围为3.8至31 MJ/m2，日蒸散量（ET0）的范围为0.8至7.9毫米。2024年，土壤温度的范围为-0.5°C（第123天）至39°C（第192天），水分势差的范围为0.02至6.0 kPa。太阳辐射范围为4.8至31.7 MJ/m2，蒸散量（ET0）的范围为1.3至8.1毫米。飞行期间的气象条件反映了温暖干燥的环境，平均土壤温度为29°C，平均水分势差为3.0 kPa。最高值出现在第18次飞行（2024年，第219天），土壤温度为36°C，水分势差为4.7 kPa；最低值出现在第24次飞行，土壤温度为22°C，水分势差为1.7 kPa。相对湿度平均为27%，通常低于32%。每次飞行的平均气象条件显示在附录D中。

使用线性回归和二次回归模型研究了Ψstem与气象变量之间的关系。评估了收获前和整个生长季节的瞬时（飞行时间）和每日气象值。由于太阳辐射和蒸散量的散点图显示非线性模式，因此包括了二次模型（附录E）。进行了回归分析，以评估四个气象参数与Ψstem之间的线性及二次关系。评估了统计显著性，并按年份、收获前和收获后时期以及整个季节对数据进行了分层。完整结果显示在附录F中，表1总结了最强的关系。表1显示了描述Ψstem与四个气象变量（空气温度（Ta）、水分势差（VPD）、太阳辐射和参考蒸散量（ET0）之间关系的最佳回归模型。这些值代表了每个因素在线性和二次模型以及瞬时和每日时间尺度上获得的最高决定系数（R2）。结果显示在观察到最强关系的年份和时期。显著性水平：T p < 0.001（***）。完整比较见表A2。收获前期间，Ψstem与气象变量之间的最强关系最为明显。两年中，使用日平均值的二次模型提供了最强的关系，解释了Ψstem变化的33%。无论使用瞬时还是每日测量数据，Ta和VPD的R2值相当，线性和二次模型之间的差异很小。

图11显示了两个生长季节中所有14棵样本树木在蓝光、绿光、红光、红边、近红外（NIR）和长波红外（LWIR）波段的反射率的时间变化。第一次飞行日期的光学波段反射率异常高，而LWIR波段的反射率低，这可能是由于当天天气阴沉所致。尽管机载光传感器可以部分补偿漫射光照，但在完全阴天的条件下收集的反射率数据不可用于定量分析[42]。因此，排除了第一次飞行的数据。两棵位于果园边缘的树木在两年中始终在光学波段表现出较高的反射率。这些树木距离碎石路5米，其较高的反射率可能受到灰尘沉积的影响。为了避免引入偏差，后续回归建模中排除了这些树木的数据。图11显示了2023年和2024年所有样本树木六个光谱波段（蓝光、绿光、红光、红边、NIR）和冠层温度（LWIR）的平均反射率时间序列。第一行显示了2023年蓝光、绿光和红光波段的平均反射率和标准差；第二行显示了2024年的相同波段。第三和第四行分别显示了2023年和2024年的红边、NIR和LWIR测量值。粗黑线代表所有14棵树木的平均反射率（或冠层温度），阴影区域代表标准差。垂直红色虚线表示每年的大致收获日期。2023年，飞行日期之间的光谱变异通常较低，红光和NIR波段除外，这些波段显示出较大的时间波动。2024年，飞行日期之间的变异在光学波段更为明显。两年中，红光和蓝光波段在整个季节呈现下降趋势，而绿光和红边波段在232天后先下降后上升。230日左右反射率明显降低，随后几个波段恢复。170日之后红光波段的反射率增加，并在收获前一天（199-201日）达到峰值。从LWIR波段得出的冠层温度（Tc）在树木间基本一致，大多数日期的差异小于2°C。2023年228日（停止灌溉超过三周）和2024年170日（灌溉与飞行同日）的变异较大。

在整个数据集（n = 330）中，Ψstem与各个光谱波段之间的相关性较弱（图12）。绿光波段（R2 = 0.20）、冠层温度（Tc = 0.20）、NIR（R2 = 0.19）和红光波段（R2 = 0.14）获得了最高的R2值。按飞行日期评估关系时，红光和绿光波段的R2值在灌溉后大约六天进行飞行时最高。图12显示了收获前后六个光谱波段（蓝光、绿光、红光、红边、NIR和LWIR）中Ψstem与反射率之间的关系。绿色点代表收获前的测量值，黄色点代表收获后的测量值。实线表示拟合收获前数据（绿色）的线性回归，虚线表示拟合收获后数据（黄色）的线性回归，点划线表示使用所有观测值进行的整体拟合。相应的R2和p值显示在每个面板中。收获前数据显示出Ψstem与反射率之间明显更强的关系。红光、绿光和NIR波段的R2值分别为0.51、0.33和0.29。收获后的Ψstem在任何光谱波段中均无显著关系，除了LWIR（R2 = 0.12）。鉴于收获后的关系较弱，回归建模仅使用了收获前的数据（n = 141）。为了评估早期测量是否由于较高的（较不负的）Ψstem值而不成比例地影响了收获前的关系，在将收获前的Ψstem范围限制在收获后的范围（-0.7至-1.2 MPa）后重新计算了R2值。结果见表2。表2显示了在限制Ψstem范围为重叠范围（-0.7至-1.2 MPa）时，多光谱波段与Ψstem之间的决定系数（R2）。如预期，限制Ψstem范围降低了所有光谱波段收获前时期的R2值。然而，对于大多数波段，收获前的关系始终高于收获后的关系，LWIR除外。这些结果表明，尽管部分较强的收获前表现可归因于较宽的早期Ψstem范围，但在相似的Ψstem值内仍存在季节性差异。

在评估的植被指数中，红色色坐标（RCC）与Ψstem之间的相关性最强（R2 = 0.68）。另外三个指数（GARI、MSR和SR）的表现也相似，R2值约为0.66。热指标Tc ? Ta也改善了模型性能，R2达到0.35。表3列出了被认为最适合估计Ψstem的植被指数和热指标。表3显示了与Ψstem关系最强的植被指数和热指标，包括指数名称、公式、R2值和主要文献来源。常用植被指数和作为Ψstem潜在估计器的指数在两年间的季节动态不同，如图13所示。2023年，NDVI、GARI和MSR在季节末逐渐增加，而RCC随时间下降。NDRE相对稳定，仅有轻微波动。相比之下，2024年表现出更强的季节内变异。NDVI和GARI在季节中期（第180-200天）急剧下降，然后部分恢复，而NDRE和RCC在季节中期出现明显峰值。图13显示了2023年和2024年所有样本树木六个植被指数（NDVI、NDRE、RCC、GARI、MSR）和冠层温度差（Tc ? Ta）的平均值时间序列。第一行显示了2023年的平均值和标准差；第二行显示了2024年的相同指数。第三和第四行分别显示了2023年和2024年的GARI、MSR和Tc ? Ta测量值。粗黑线代表所有14棵树木的平均值，阴影区域代表标准差。垂直红色虚线表示每年的大致收获日期。2023年，飞行日期之间的光谱变异通常较低，红光和NIR波段除外，这些波段显示出较大的时间波动。2024年，飞行日期之间的变异在光学波段更为明显。两年中，红光和蓝光波段在整个季节呈现下降趋势，而绿光和红边波段在232天后先下降后上升。230日左右反射率明显降低，随后几个波段恢复。170日之后红光波段的反射率增加，并在收获前一天（198日）达到峰值。从LWIR波段得出的冠层温度（Tc）在树木间基本一致，大多数日期的差异小于2°C。2023年228日（停止灌溉超过三周）和2024年170日（灌溉与飞行同日）的变异较大。

在完整的数据集中，Ψstem与各个光谱波段之间的相关性较弱（图12）。绿光波段（R2 = 0.20）、冠层温度（Tc = 0.20）、NIR（R2 = 0.19）和红光波段（R2 = 0.14）获得了最高的R2值。按飞行日期评估关系时，红光和绿光波段的最高R2值出现在灌溉后大约六天进行飞行时。图12显示了收获前后六个光谱波段（蓝光、绿光、红光、红边、NIR和LWIR）中Ψstem与反射率之间的关系。绿色点代表收获前的测量值，黄色点代表收获后的测量值。实线表示拟合收获前数据（绿色）的线性回归，虚线表示拟合收获后数据（黄色）的线性回归，点划线表示使用所有观测值进行的整体拟合。相应的R2和p值显示在每个面板中。收获前数据显示出Ψstem与反射率之间明显更强的关系。红光、绿光和NIR波段的R2值分别为0.51、0.33和0.29。Tc也有所改善（R2 = 0.23）。收获后的Ψstem在任何光谱波段中均无显著关系，除了LWIR（R2 = 0.12）。鉴于收获后的关系较弱，回归建模仅使用了收获前的数据（n = 141）。为了评估早期测量是否由于较高的（较不负的）Ψstem值而不成比例地影响了收获前的关系，在将收获前的Ψstem范围限制在收获后的范围（-0.7至-1.2 MPa）后重新计算了R2值。结果见表2。表2显示了在限制Ψstem范围为重叠范围（-0.7至-1.2 MPa）时，多光谱波段与Ψstem之间的决定系数（R2）。如预期，限制Ψstem范围降低了所有光谱波段收获前时期的R2值。然而，对于大多数波段，收获前的关系始终高于收获后的关系，LWIR除外。这些结果表明，尽管部分较强的收获前表现可归因于较宽的早期Ψstem范围，但在相似的Ψstem值内仍存在季节性差异。

在评估的植被指数中，红色色坐标（RCC）与Ψstem之间的相关性最强（R2 = 0.68）。另外三个指数（GARI、MSR和SR）的表现也类似，R2值约为0.66。热指标Tc ? Ta也改善了模型性能，R2达到0.35。表3列出了被认为是估计Ψstem的最佳候选植被指数和热指标。表3显示了与Ψstem关系最强的植被指数和热指标，包括指数名称、公式、R2值和主要文献来源。常用植被指数和作为Ψstem潜在估计器的指数的季节动态在两年间有所不同，如图13所示。2023年，NDVI、GARI和MSR在季节末逐渐增加，而RCC随时间下降。NDRE相对稳定，仅有轻微波动。相比之下，2024年表现出更强的季节内变异。NDVI和GARI在季节中期（第180-200天）急剧下降，然后部分恢复，而NDRE和RCC在季节中期出现明显峰值。图13显示了2023年和2024年所有样本树木六个植被指数（NDVI、NDRE、RCC、GARI、MSR）和冠层温度差（Tc ? Ta）的平均值时间序列。第一行显示了2023年的平均值和标准差；第二行显示了2024年的相同指数。第三和第四行分别显示了2023年和2024年的GARI、MSR和Tc ? Ta测量值。粗黑线代表所有14棵树木的平均值，阴影区域代表标准差。垂直红色虚线表示每年的大致收获日期。冠层温度差（Tc ? Ta）在两年间也有所不同。2023年，Tc ? Ta在季节中期显示出明显最小值，而2024年的值更加多变，特征是短期峰值。总体而言，2024年的指数时间变异大于2023年。使用SR方法估计Ψstem的模型性能
以下是使用SR方法找到的六个估计Ψstem的方程：
Ψ???????????=0.2364?3.4979?????????, (10)
Ψ???????????=0.7?3.2069??????????0.8164??????0, (11)
Ψ???????????=0.0297?(?????????)?2.96?????????, (12)
Ψ???????????=0.0045?(?????????)2?2.8789?????????, (13)
Ψ???????????=????????15?c?m+0.228?????????+0.0047??????1.87, (14)
Ψ???????????=11.9????+0.0272??????????????????????????????15?c?m?1.9 (15)
其中，RCC是红色色度坐标指数；ET0是参考蒸散量；Tc是冠层温度；Ta是空气温度；SWC15是15厘米深度的土壤体积含水量；MSR是修正后的简单比率指数；LWIR是长波红外波段。

为了确认符号回归保留的每个预测因子对模型性能有显著贡献，使用方程(10)–(15)的固定功能形式进行了OLS回归。所有保留的系数都高度显著（p < 0.001），证实每个变量对Ψstem的预测具有统计意义上的解释力（附录G）。

评估指标在表4中呈现。模型在不同的数据可用性场景下进行了评估，从仅使用RGB导出指数的模型到结合多光谱、热学、气象和土壤含水量变量的模型。所有六个模型的R2值范围从0.67到0.80，而RMSE仅略有变化，从0.11到0.08 MPa。MAE保持在0.09 MPa以下，MAPE值均低于12%。表4. 使用符号回归开发的六个Ψstem估计方程的复杂性和性能指标。模型适用性因数据可用性而异。当仅使用RGB图像时，仅基于RCC的方程(10)提供了强大的预测性能。当包含气象数据时，方程(9)提高了准确性。通过结合热信息，方程(12)和(13)进一步提高了准确性，而额外使用土壤湿度的方程(14)和(15)取得了最高的性能指标。这些模型涵盖了从仅RGB输入到结合光谱、热学、气象和土壤湿度测量的各种数据可用性场景。与方程(10)相比，方程(15)显著降低了RMSE（ΔRMSE = 0.03 MPa；配对t检验p < 0.001）。然而，方程(15)与较简单的模型（如方程(10)）之间的差异在统计上不显著（p > 0.05），表明模型复杂性的增加导致收益递减。

根据信息标准，方程(15)对数据的拟合最好，显示出最低的AIC（?672.81）和BIC（?661.02）值。然而，方程(13)也表现出良好的预测性能（R2 = 0.77；RMSE = 0.09 MPa）以及较低的AIC和BIC值，在模型性能和复杂性之间提供了有利的折中，用更少的输入变量实现了可比的预测准确性。

图14显示了六个提出的方程之间观测到的Ψstem与估计值的一致性。所有模型都产生了测量值和估计值之间的强线性关系，残差散布在回归线周围均匀分布。这一模式通过残差的Q-Q图得到了确认（附录H）。图14. 六个提出的估计模型（方程(8)–(13)）的观测值与预测值Ψstem。每个面板显示了拟合的回归线（蓝色）以及相应的R2和p值。所有预测因子都基于收获前的数据（n = 141）。使用留一树法交叉验证程序进一步确认了预测能力（表5）。所有模型的表现都很好，RMSE值在0.11到0.09 MPa之间，MAE值在0.08到0.06 MPa之间，MAPE值低于12%，与表3中的结果非常相似。空模型使用整体平均Ψstem作为预测，表现明显较差（RMSE = 0.19 MPa），这表明结合生理和环境预测因子的价值。在候选模型中，方程(13)在交叉验证期间的总体误差最低（RMSE = 0.09 MPa，MAPE = 8.2%），其性能与广义加性模型（GAM）相似。表5. 六个估计模型的留一树法交叉验证性能，与广义加性模型（GAM）和空模型进行比较。指标包括RMSE、MAE和MAPE。当使用2024年数据集重新校准系数并应用于2023年的观测数据时，预测性能保持在交叉验证期间观察到的相同误差范围内，证明了SR衍生结构的时间稳定性。重要的是，重新校准的系数相对于原始公式（方程(10)–(15)）保持了幅度和符号的一致性，进一步支持了所识别数学结构的时间稳健性。重新校准的系数和验证性能总结在表6中。表6. 使用2024年数据重新校准的系数及其在2023年观测数据上的验证性能。为了评估模型的泛化能力，使用来自UT州Santaquin果园的独立数据集对外部验证了六个SR衍生方程（表7）。在Ψstem测量期间，每小时平均Ta、VPD和ET0分别为27.2 °C、2.74 kPa和0.65 mm h?1。表7. 使用来自UT州Santaquin果园的外部验证数据集评估的六个Ψstem方程的性能指标。尽管站点条件和样本量较小，大多数模型仍保持了良好的预测性能。方程(12)表现出最佳的整体性能，紧随其后的是方程(13)。相比之下，方程(10)和(11)显示出更高的误差水平（RMSE > 0.12 MPa）。总体而言，结合RCC和Tc ? Ta的模型在外部果园数据集中表现更优。

图15展示了外部验证数据集中预测值和观测值Ψstem之间的一致性。方程(12)和(13)显示出强烈的线性对应关系，残差在回归线周围分散有限，预测行为在整个评估范围内一致。相比之下，方程(14)和(15)显示出更大的散布和偏离1:1关系。没有观察到模型的系统性过度或低估现象，支持了它们在独立站点条件下的稳健性。图15. 在外部验证果园获得的六个提出的估计模型（方程(10)–(15)）的观测值与预测值Ψstem。每个面板显示了拟合的回归线（蓝色）以及相应的R2和p值。灰色点代表外部验证数据集中的单个观测值。

图16展示了2023年7月2日（DOY 183）Kaysville果园中Ψstem的空间分布，由六个SR衍生方程估计得出。方程(10)和(11)产生了更异质的空间模式，特别是在果园的中央和西部区域，Ψstem值更为负。相比之下，方程(12)和(13)生成了更平滑的空间梯度，极端负值较少，整体分布较为均匀。方程(14)倾向于生成较高的Ψstem值，导致大面积的Ψstem较高。图16. 使用2023年7月2日（DOY 183）在UT州Kaysville收集的多光谱图像生成的六个提出的估计方程的Ψstem空间图。图17展示了2024年7月16日（DOY 198）生成的Ψstem图。方程(10)–(13)产生了相似的空间模式，一致地识别出果园西部区域的Ψstem值较低的区域，而东部行区域的Ψstem值相对较高。相比之下，方程(14)和(15)生成的Ψstem估计值总体上较少负，导致果园的空间分布更加均匀。在所有模型中，果园的东部部分的Ψstem估计值高于西部部分。尽管各方程之间的幅度存在差异，但总体空间模式大体一致。图17. 使用2023年7月16日（DOY 198）在UT州Kaysville收集的多光谱图像生成的六个提出的估计方程的Ψstem空间图。

4. 讨论
4.1. 影响Ψstem的环境和土壤因素
实验地点表现出最小的空间异质性，有助于分离Ψstem的生理和光谱驱动因素。已有文献记录了粘土含量与土壤ECa之间的线性关系[42,43]。研究区域内的土壤ECa值范围较窄，表明土壤质地均匀，与美国农业部土壤调查分类中的单一细砂壤土类型一致。由于ECa与粘土含量[48]密切相关，因此与持水能力相关，其空间均匀性支持了Ψstem的变化不是由土壤类型差异驱动的假设。

土壤含水量模式也与预期的根区动态一致。酸樱桃的根集中在10至50厘米深度[49]，SWC与Ψstem之间的最强关系出现在15至45厘米深度，特别是在收获前期间。土壤湿度影响调节气孔导度（gs）的根信号[50]，而gs与Ψstem高度相关[51]。在我们的研究中，SWC与Ψstem呈正线性关系。在葡萄藤中也观察到了类似的关系[13]。这种关系在收获前期间更为强烈。

气象变量与Ψstem有中等程度的相关性，其中收获前期间Ψstem与每日ET0值之间的相关性最强（R2 = 0.36）。这表明随着ET0的增加，Ψstem变得更负，直到达到一个阈值，之后Ψstem的变化变得更加渐进。参考文献[52]认为，在充分灌溉条件下，任何给定环境条件下都可以预期Ψstem的最大值。在年轻橄榄树中也报告了Ψstem与瞬时Ta之间的负线性关系[53]。这些发现共同表明，大气需求的增加触发了酸樱桃的气孔调节，与其他木本多年生作物的观察结果一致[54]。

4.2. 影响光谱和环境响应的季节性生理驱动因素
Ψstem的季节性模式和预测因子关系的时间变化可能反映了潜在的生理变化。收获是一个主要的生理事件，显著改变了源-汇关系。果实的移除是一个主要的汇，通常会减少蒸腾作用[55]，因为汇的强度影响气孔导度[56]。在酸樱桃中，机械收获可能会进一步影响收获后的水分关系。机械收获期间常见的树皮损伤[57]可能会暂时减少水分和养分的流动，直到愈伤组织形成恢复血管流动[58]。这些综合效应可能导致了收获后Ψstem与光谱和环境变量之间关系的减弱。相比之下，参考文献[59]报告了收获后甜樱桃中Ψstem与植被指数之间的相关性增加。甜樱桃是手工收获的。在其他作物中也报告了Ψstem与反射率之间季节性的变化。在柑橘类作物中，参考文献[26]发现果实细胞分裂阶段没有显著相关性，而在果实快速生长阶段则出现了更强的关系。此外，酸樱桃树产生两种类型的叶子：早季短枝叶和后期出现的新叶。叶片年龄强烈影响光谱反射率[60]，如果不考虑这一点，可能会导致遥感分析的误解[61]。

4.3. Ψstem的光谱和热指标
在各个光谱带中，红色带与Ψstem的相关性最强（R2 = 0.51）。该区域的反射率与叶绿素浓度密切相关[62]。尽管植物胁迫倾向于降低叶绿素，但由于叶对胁迫的光学响应的普遍性，使用遥感识别具体原因可能较为困难[63]。在之前的研究中也报告了Ψstem与UAV收集的光谱带反射率之间的类似关系。在缺水灌溉下的葡萄藤中，参考文献[23]发现红色和绿色带与Ψstem的相关性最强，分别为R2 = 0.64和0.62，发生在收获前期间。同样，参考文献[64]报告在柑橘类作物中，红色带与Ψstem的相关性最强，特别是在果实成熟阶段。在更粗的空间分辨率下，参考文献[23]表明Sentinel 2光谱带对解释葡萄藤中Ψstem的变异性更为重要，尤其是红色、红边和NIR带。

在评估的植被指数中，RCC在收获前期间预测Ψstem方面优于基于NIR的指数（R2 = 0.68）。RCC之前已被与冠层光合作用[65]和物候动态[66]相关联。由于RCC仅来自可见光带，并代表了RGB光谱中红色反射的相对优势，因此它对叶绿素浓度和冠层颜色动态的变化特别敏感，这些变化可能先于NIR区域可检测到的结构退化。此外，在收获前期间，红色果实的存在可能增加了冠层内的红色反射。因为果实发育和色素沉着与植物水分状态相关[67,68]，RCC可能间接捕捉到了与果实成熟相关的胁迫变化。这些因素共同解释了在本研究果园条件下观察到的RCC与Ψstem之间更强的关系。

仅热成像与Ψstem的关系较弱，与甜樱桃中的发现一致[35]。尽管Tc ? Ta相对于Tc提高了性能，但这种提升仍然较为有限。这种较低的敏感性可能反映了与成熟甜樱桃树相似的生理行为，即水分胁迫主要通过气孔关闭来调节，从而降低了Tc ? Ta的敏感性[69]。虽然作物水分胁迫指数（CWSI）仍然是标准的基于温度的指标[47]，但其应用需要定义树冠温度的上限和下限基准，而不同的研究方法在这些基准的确定上存在很大差异[70,71,72,73]，这在实际操作中具有挑战性。

4.4. 符号回归模型的性能、复杂性和可解释性
符号回归产生了准确且易于解释的方程式，用于估算Ψstem。模型性能（R2 = 0.67–0.80，RMSE = 0.11–0.09 MPa）与之前发表的关于杏仁[30]、开心果[29]、葡萄园[32]和甜樱桃[35]的机器学习模型相当。然而，正如[29]所指出的，许多机器学习模型由于可解释性有限而难以实际应用，最终用户无法使用。符号回归相对于传统机器学习方法的一个关键优势是其模型的透明度。所得方程式明确描述了变量之间的功能关系，例如当Tc ? Ta与光谱和环境预测因子结合使用时，其解释能力得到了增强。在本研究中，随着额外输入（温度、气象和土壤水分含量SWC）的加入，模型复杂性有所增加，但多光谱+温度组合之外的性能提升并不显著。这些有限的性能改进并不足以证明额外测量所带来的成本和操作负担是合理的。符号回归方程式的可解释性和可重复性增强了其在种植者和灌溉管理者中的验证、采纳和信任度。

4.5. 灌溉管理的实际意义
在果树中实施调控性亏缺灌溉的一个主要挑战是有效监测植物的水分胁迫。空间明确的Ψstem地图可以为果园内部的变异性提供有价值的见解，并支持特定地点的灌溉决策。基于无人机的图像可以在灵活的时间分辨率下获取，特别是在夏季云量较少的干旱地区，如犹他州。使用本研究开发的方程式生成的Ψstem地图可以帮助种植者通过评估水分胁迫水平来优化灌溉决策。方程式（10）在气象或土壤湿度数据不可用的情况下特别适用于实际操作。尽管方程式（12）和（13）在外部验证果园中的MAPE值（约6%）低于方程式（10）（12.5%），但RMSE的改善相对较小（0.06 MPa）。这表明使用简单得多的模型结构仍然可以实现可接受的预测性能。重要的是，方程式（10）完全依赖于RGB图像，无需获取温度数据。温度传感器会显著增加系统成本，并需要额外的校准程序，从而增加操作复杂性。相比之下，经过辐射校准的RGB图像更容易获取，成本更低，且可以通过标准的航空测量协议进行收集。例如，一个结合了RGB和辐射温度传感器的相机成本可能高达约20,000美元，而单独的RGB传感器通常价格在2,000至5,000美元之间。

4.6. 局限性
几个因素可能导致模型中存在无法解释的变异性。本研究的局限性在于，Ψstem与天气、土壤水分含量SWC、温度和多光谱图像之间的有意义关系仅在收获前能够被检测到。收获后似乎发生了生理解耦，限制了模型在收获后的适用性。未来的研究应重点关注收获后的生理过程及其对植物水分状态的影响。果园内的营养变化也会改变叶片叶绿素浓度，从而影响反射率，这种影响与水分状态无关。高氮水平通常会导致更高的叶绿素含量[74]，从而在可见光反射率中产生更大的吸收[75]。
将基于无人机的Ψstem地图用于灌溉调度时，需要为酸樱桃定义特定的水分胁迫阈值。虽然酸樱桃的阈值尚未确定，但已经为甜樱桃品种提出了参考值，例如‘Prime Giant’为-1.3 MPa[18]，‘Summit’为-1.5 MPa[76]。研究表明，适度的水分胁迫对酸樱桃[77]和甜樱桃[35]的产量没有显著影响。在解释接近水分胁迫阈值的Ψstem值时，应考虑所提模型的预测不确定性。例如，如果将操作灌溉触发点设定为-1.3 MPa，那么模型误差约为0.08 MPa，仅占阈值大小的6%。实际上，这种不确定性水平相对于Ψstem的典型季节性变化范围（本研究中约为-0.4至-1.25 MPa）来说相对较小，除非预测值非常接近选定的阈值，否则不太可能改变灌溉决策。在这种情况下，种植者可以加入安全边际或重复测量以减少决策不确定性。

5. 结论
本研究证明，可以通过符号回归结合高分辨率无人机多光谱图像和环境数据可靠地估算Ψstem。在收获前期间，观察到Ψstem与光学反射率之间存在强烈的相关性。在评估的光谱预测因子中，RCC是最具信息量的植被指数，优于单个光谱带和温度指标。符号回归生成了六个基于方程的模型，涵盖了从仅使用RGB图像到包括温度、气象（Ta和ET0）以及15厘米深度的土壤体积含水量在内的多种数据可用性场景。所有模型的预测性能都很高（R2 = 0.67至0.80；RMSE 0.11至0.08 MPa），并且通过增加模型复杂性获得了小幅的性能提升。使用独立果园数据集进行的外部验证证实了模型的稳健性和泛化能力，特别是对于结合了温度、多光谱和天气信息的方程式。这项工作的一个关键贡献是证明了符号回归可以实现与“黑箱”机器学习模型相当的性能，同时保持完全的数学可解释性。所得方程式明确揭示了光谱、温度和环境变量与Ψstem之间的功能关系，从而实现了透明性、可重复性和实际应用的可能性。这种可解释性对于将模型应用于灌溉决策支持中的Ψstem映射至关重要。能够生成果园规模的Ψstem地图为替代劳动密集型的压力室测量方法提供了有前景的方案。因此，未来的研究应探讨植物营养和Ψstem对树冠反射率的综合影响，以更好地分离营养和水分胁迫信号。最后，需要为酸樱桃确定特定的水分胁迫阈值Ψstem值，以实现实际的灌溉调度。总体而言，本研究确立了符号回归作为一种强大且易于解释的框架，用于基于植物的水分状态估算，并支持将其整合到特种作物的精准灌溉管理中。