4、光学遥感反演方法：基于植被指数的土壤湿度反演（TVDI、NDWI），经验模型与半经验模型
各位同学，咱们今天聊点实在的。光学遥感反演土壤湿度，说白了就是利用卫星拍到的可见光-近红外影像，把土壤里的水分信息给“算”出来。你可能会问：光靠眼睛看照片，怎么能知道土里有多少水呢？嗯，这里面的门道，我做了这么多年项目，总结下来就两招——植被指数法和经验/半经验模型。
我个人习惯把这两种方法比作“侦探破案”。植被指数法像是看现场留下的脚印（植被状态），而经验模型则是直接套用以前破过的案子的规律。今天咱们就重点讲讲TVDI和NDWI这两个“老牌侦探”。
4.1 为什么用植被指数反演土壤湿度？
先讲个背景。纯裸土反演湿度，其实相对简单——土壤越湿，反射率越低。但一旦有植被覆盖，事情就复杂了。植被的“绿色信号”会盖住土壤的信号，就像你隔着窗帘看外面，很难看清路面的干湿。
我在项目中遇到过这种情况：用单一波段反演玉米地的土壤湿度，结果完全不准。后来才发现，植被越茂密，误差越大。所以，我们需要借助植被指数来“剥离”植被的影响，或者反过来利用植被对水分的响应来推算土壤湿度。
4.2 TVDI：温度植被干旱指数
TVDI，全称Temperature Vegetation Dryness Index。这个名字有点长，但原理很直观。它基于一个现象：在相同植被覆盖度下，土壤越干，地表温度越高。
你想想看，湿土蒸发强，带走热量，地表就凉快；干土蒸发弱，太阳一晒就烫手。TVDI就是利用这个温差来估算干旱程度。
4.2.1 原理与公式
TVDI的计算需要两个输入：地表温度（LST）和归一化植被指数（NDVI）。具体做法是：
- 在NDVI-LST散点图上，找到每个NDVI区间内的最高温度和最低温度。
- 拟合出“湿边”和“干边”两条线。
- 计算每个像元在干湿边之间的相对位置。
公式长这样：
TVDI = (Ts - Ts_min) / (Ts_max - Ts_min)其中：
Ts：像元的地表温度Ts_min：该NDVI对应的湿边温度（最低温）Ts_max：该NDVI对应的干边温度（最高温）
TVDI的取值范围是0到1。越接近1，说明越干；越接近0，说明越湿。
4.2.2 代码实现（Python示例）
下面是我常用的TVDI计算流程，用Python和GDAL实现。注意，这里假设你已经有了LST和NDVI的栅格数据。
import numpy as np
from osgeo import gdal
def calculate_tvdi(lst_path, ndvi_path, out_path):
# 读取数据
lst_ds = gdal.Open(lst_path)
ndvi_ds = gdal.Open(ndvi_path)
lst = lst_ds.GetRasterBand(1).ReadAsArray().astype(np.float32)
ndvi = ndvi_ds.GetRasterBand(1).ReadAsArray().astype(np.float32)
# 设置有效值范围（剔除异常）
valid_mask = (ndvi > 0) & (ndvi < 1) & (lst > 0)
# 按NDVI区间统计干湿边
ndvi_bins = np.arange(0, 1.01, 0.01)
ts_max = np.full_like(ndvi_bins, np.nan)
ts_min = np.full_like(ndvi_bins, np.nan)
for i in range(len(ndvi_bins)-1):
mask = (ndvi >= ndvi_bins[i]) & (ndvi < ndvi_bins[i+1]) & valid_mask
if np.sum(mask) > 100: # 至少100个像元
ts_max[i] = np.percentile(lst[mask], 98) # 用98分位数代替最大值，避免噪声
ts_min[i] = np.percentile(lst[mask], 2) # 用2分位数代替最小值
# 线性拟合干边和湿边（这里简化处理，实际可用多项式）
# 干边：Ts_max = a + b * (1 - NDVI)
# 湿边：Ts_min = c + d * (1 - NDVI)
# 计算TVDI
tvdi = np.full_like(lst, np.nan)
for i in range(len(ndvi_bins)-1):
mask = (ndvi >= ndvi_bins[i]) & (ndvi < ndvi_bins[i+1]) & valid_mask
if np.sum(mask) > 0:
tvdi[mask] = (lst[mask] - ts_min[i]) / (ts_max[i] - ts_min[i])
# 限制范围
tvdi = np.clip(tvdi, 0, 1)
# 保存结果
driver = gdal.GetDriverByName('GTiff')
out_ds = driver.Create(out_path, lst_ds.RasterXSize, lst_ds.RasterYSize, 1, gdal.GDT_Float32)
out_ds.SetGeoTransform(lst_ds.GetGeoTransform())
out_ds.SetProjection(lst_ds.GetProjection())
out_ds.GetRasterBand(1).WriteArray(tvdi)
out_ds.FlushCache()
return tvdi4.3 NDWI：归一化水体指数
NDWI大家可能更熟悉，它原本是用来提取水体的。但你可能不知道，NDWI对土壤水分也很敏感。尤其是植被冠层含水量，和土壤湿度有很强的相关性。
NDWI的公式是：
NDWI = (Green - NIR) / (Green + NIR)或者用短波红外（SWIR）版本：
NDWI = (NIR - SWIR) / (NIR + SWIR)我个人更推荐用SWIR版本。因为SWIR对水分吸收非常敏感，植被含水量稍微变化，反射率就会明显波动。我在水稻田监测项目中，就是用SWIR版本的NDWI来估算稻田土壤湿度，效果比NDVI好很多。
3.3.1 NDWI与土壤湿度的关系
NDWI和土壤湿度之间，通常呈正相关。土壤越湿，植被吸收的水分越多，叶片含水量越高，SWIR反射率越低，NDWI值就越大。
但要注意：NDWI反映的是冠层含水量，不是直接测土壤。所以它更适合植被覆盖度较高的区域。对于稀疏植被或裸土，NDWI的效果会大打折扣。
| 指数 | 适用场景 | 优点 | 缺点 |
|---|---|---|---|
| TVDI | 中低植被覆盖区 | 物理意义明确，对干旱敏感 | 需要同步LST数据，受云影响大 |
| NDWI（SWIR版） | 高植被覆盖区 | 对冠层含水量敏感，计算简单 | 不能直接测土壤，受植被类型影响 |
4.4 经验模型与半经验模型
讲完指数法，咱们再聊聊模型。经验模型和半经验模型，说白了就是“找规律”和“套公式”的区别。
4.4.1 经验模型
经验模型完全靠数据驱动。你拿一堆实测的土壤湿度数据，再配上对应的遥感反射率或指数，用回归分析拟合出一个公式。
比如：
SM = a * NDVI + b * LST + c或者：
SM = a * TVDI + b这种模型的好处是简单、快速。坏处是——换一个地方、换一个季节，参数可能就不灵了。我记得在东北黑土区做过一个项目，用华北平原拟合出来的经验模型，结果误差大得离谱。后来老老实实重新采了样，重新拟合。
4.4.2 半经验模型
半经验模型介于经验和物理模型之间。它有一定的物理基础，但部分参数还是靠经验拟合。最典型的就是植被调整后的土壤湿度指数。
举个例子，我常用的一个半经验模型是：
SM = SM_bare * (1 - f_veg) + SM_veg * f_veg其中：
SM_bare：裸土部分的湿度（用短波红外反演）SM_veg：植被部分的湿度（用NDWI或TVDI反演）f_veg：植被覆盖度（用NDVI估算）
这个模型的好处是：它考虑了植被和土壤的混合效应。你想想看，一个像元里既有土又有草，直接用一个指数肯定不准。但分开算再加权平均，精度就上去了。
4.4.3 模型对比
| 模型类型 | 输入数据 | 可移植性 | 精度 | 计算复杂度 |
|---|---|---|---|---|
| 经验模型 | 遥感指数 + 实测数据 | 低（需重新拟合） | 中等（依赖样本） | 低 |
| 半经验模型 | 遥感指数 + 先验知识 | 中等 | 较高 | 中等 |
| 物理模型 | 辐射传输模型参数 | 高 | 高 | 高 |
4.5 实战中的选择建议
说了这么多，到底该用哪个？我个人的经验是：
- 如果研究区植被稀疏（比如草原、荒漠边缘），优先用TVDI。它对土壤水分的响应更直接。
- 如果植被茂密（比如农田、森林），优先用NDWI（SWIR版）。因为TVDI在茂密植被下，地表温度信号会被冠层遮挡。
- 如果时间紧、任务急，用经验模型快速出结果。但一定要留一部分实测数据做验证。
- 如果想发文章或做业务化运行，老老实实搞半经验模型。虽然前期工作量大，但鲁棒性好。
好了，这一章的内容就到这儿。下一章咱们聊聊微波遥感反演方法，那玩意儿不怕云，很有意思。