简介这是一篇发表于《农业工程学报》2021年的学术论文PDF面向遥感、农业工程及机器学习领域的研究者聚焦微波遥感反演农田地表土壤水分时植被覆盖干扰的难题提出了基于特征选择和GA-BP神经网络的多源遥感反演方法。整包仅含1个PDF文件大小约5.54MB已有350人学习参考。论文详细阐述了Sentinel-1微波数据与Sentinel-2光学数据的预处理、21个特征参数提取、差分进化特征选择DEFS选出10个最优参数、主成分分析降维、BP神经网络建模及遗传算法优化节点权值的完整流程。通过实测数据对比该方法反演结果的决定系数为0.7893均方根误差为0.0287 cm3/cm3相比未引入DEFS与PCA的GA-BP模型决定系数提升0.2157均方根误差降低0.0295 cm3/cm3。读者可从中系统掌握多源遥感数据融合、特征优选与神经网络优化的具体实现路径为同类研究提供可直接借鉴的方法与实验依据。1. 一个反演项目最先要过的坎数据进来了但精度上不去做农田墒情监测的人应该都有这种体验卫星影像要多少有多少但真拿去反演土壤水分结果忽高忽低验证集R²迟迟上不了0.7。原因多半不在算法而在输入——单一光学影像对湿度不敏感单一雷达信号又受地表粗糙度干扰。这个标题给了一条明确路线先把光学、雷达、地形等多源遥感特征堆到一起再用特征选择筛掉冗余维度最后交给GA-BP神经网络建反演模型。遗传算法在里面做的事很具体给BP找一组好的初始权值和阈值治它容易陷进局部最优的老毛病。一套走下来常见的改善幅度是特征从三十维降到十几维R²从0.6拉到0.8附近。这篇笔记适合正在做遥感定量反演、智慧灌溉决策或者刚接手土壤水分项目的人参考。2. 多源遥感输入的组织先别调网络把数据对齐才是第一道关很多人拿到数据第一件事就是开工具箱跑模型结果反演精度上不去回头查才发现是输入数据本身对不齐。多源遥感反演土壤水分这件事里数据准备占掉整个工作量的一半以上而且这步出的问题后期基本没法补救只能重来。2.1 多源输入到底在凑哪些特征土壤水分反演用多源是因为不同传感器对土壤水分的响应机理不同互补之后模型才稳。常用组合我一般按四类收集数据源典型特征对土壤水分的响应光学多光谱蓝、绿、红、近红外、短波红外波段不同波段组合可反映地表湿度、植被覆盖和土壤类型差异热红外地表温度、温差指数干旱时地温偏高湿度大时地温偏低是强相关变量合成孔径雷达SARVV、VH后向散射强度对土壤介电常数敏感湿度越大后向散射越强但受粗糙度和植被影响地形辅助高程、坡度、坡向影响径流和水分再分布缓坡和低洼处持水性不同衍生指数NDVI、NDWI、TVDI、湿度指数大部分是波段运算结果用来表达植被水分和土壤湿度状态收集时注意一个原则宁可特征多不可缺类型。光学只能看表层SAR能穿透植被冠层一部分地形决定的是水分的空间分布规律三者在反演模型里各守一段。后续特征选择阶段会把没用的挑出去但输入阶段如果某个数据源缺失后面再怎么选也补不回来。2.2 预处理和对齐不同分辨率、不同投影怎么统一多源数据拿回来时分辨率、投影、获取时间全都不一样直接叠用就会在边界上出现“数据打架”。我一般按这套顺序处理光学影像做辐射定标和大气校正把DN值换成地表反射率这一步不做NDVI和TVDI全都不准。SAR数据做辐射地形校正和斑点滤波滤波窗口常用7×7选Refined Lee滤波既去噪又保留边缘。统一投影坐标系到UTM再按研究区裁剪掩膜。重采样到同一像元分辨率光学波段和地形数据统一到目标栅格SAR也重采样到同网格我用10米或30米看最终输出要求。时间对齐土壤水分变化快影像获取时间差控制在24小时以内最好。光学和SAR不在同一天拍的就取同一旬内云量最少的两景记录时间差作为误差来源。计算衍生指数并叠加成多波段栅格堆栈。特征矩阵的堆叠我一般用Python完成rasterio读栅格然后按波段压成一个数组示范代码如下import rasterio import numpy as np files { B02: path/to/sentinel2_B02.tif, # 蓝波段 B04: path/to/sentinel2_B04.tif, # 红波段 B08: path/to/sentinel2_B08.tif, # 近红外 VV: path/to/sentinel1_VV.tif, # 雷达同极化 VH: path/to/sentinel1_VH.tif, # 雷达交叉极化 dem: path/to/dem.tif # 高程 } stack [] for name, path in files.items(): with rasterio.open(path) as src: band src.read(1).astype(np.float32) stack.append(band) # 得到 (n_features, height, width) 的三维数组 feature_cube np.stack(stack, axis0) # 展平成样本矩阵每个像元一行每列一个特征 n_features, h, w feature_cube.shape X feature_cube.reshape(n_features, h * w).T # 样本数 × 特征数这段逻辑里关键点是所有文件必须已经过重采样否则np.stack会在数组形状不一致时报错这是最常见的第一个翻车点。另一个容易忽视的是float32转换原始GeoTIFF常常是uint16格式后续算NDVI会发生整数截断所以提前转浮点。reshape那一步把空间维度压成样本维度后面模型里每一行就是一个像元的多源特征向量和对应的土壤含水量真值配对。2.3 按窗口提取邻域统计为什么单像元经常不带信息多源影像配准误差普遍存在光学和SAR之间差半个像元是常事单像元取值容易受配准误差影响。另一个问题是单像元噪声大一个像元里可能混合着植被、裸土和阴影直接拿来做输入模型会学得很吃力。常见做法是对每个像元取邻域统计量再作为特征。窗口大小我常用5×5或7×7太小降噪效果差太大把边缘细节磨平研究区地块细碎时还会把相邻地块的异质信息混进来。每个窗口算均值、标准差作为新特征from scipy.ndimage import uniform_filter, gaussian_filter # 以7×7均值窗口为例平滑后的值作为邻域中心像元的新特征 X_mean uniform_filter(feature_cube, size(1, 7, 7), modereflect) # size(1,7,7) 表示只在空间维度平滑不跨波段 # 再算一个邻域标准差反映局部异质性 X_local_std np.sqrt(uniform_filter((feature_cube - X_mean) ** 2, size(1, 7, 7), modereflect))这段代码里modereflect很关键不设置的话边界像元会被补零研究区边缘出现一条异常的暗边模型会专门学习这条边的错误规律。窗口平滑后的特征和原始特征可以拼到一起用实践中邻域均值对SAR数据降噪特别有效标准差信息对地形复杂的区域帮助明显。做这一步时留意样本量会因此膨胀原本几百个实测点做反演扩展成邻域后每个点周围都有大量相似像元后面模型评估时要注意独立性问题这个第5章会展开讲。3. 特征选择为什么三十维特征直接进BP等于让模型瞎转数据堆到三十维以上之后很多人直接扔进BP训练结果训练集损失降得很漂亮验证集一测就崩。这不是BP不行而是高维特征里冗余和噪声信息太多。特征选择这一步是为了在进网络之前先把不相关的、重复的维度去掉。3.1 特征选择不是特征提取物理意义是遥感反演的底线做特征处理时有两条路特征提取和特征选择。特征提取典型代表是PCA和偏最小二乘PLS它们会把原始波段组合成新变量每个新变量都包含所有原始波段的贡献数学上很干净但物理意义丢了。遥感反演落地时最需要的是可解释性——一个特征对应什么物理过程换一个地区换一颗卫星还能不能用。PCA之后的特征分量换一个场景就未必可迁移。特征选择保留的是原始波段和派生指数只是去掉其中一部分。比如选出NDVI、VV、地表温度这三个特征还能说清楚NDVI反映了植被覆盖对土壤水分的遮挡VV后向散射反映介电常数地表温度反映蒸发强度。这种可解释性在写技术报告和向甲方说明模型可信度时有不可替代的价值。所以做遥感反演土壤水分我几乎不用PCA做特征提取优先走特征选择路线这也正是特征工程里说的“能选就不要变”。3.2 三个常用挑选方法的适用场景皮尔逊、方差阈值、互信息特征选择特征选择不能只靠一种方法遥感特征里同时存在线性相关和非线性相关单一筛选标准会漏掉或误删。我常用三种方法配合方差阈值是最先做的它不关心特征和土壤水分的关系只看特征本身有没有变化。有些波段在研究区整个时段内数值几乎不动对反演没有任何区分能力阈值的具体取值要看特征值的量纲标准化之后一般设0.01。方差太低的直接剔除。皮尔逊相关系数处理的是线性关系。把每个特征和土壤水分实测值算相关系数|r|小于0.15的认为线性关系太弱考虑剔除。同样也要看特征之间的相关性|r|大于0.85的两个特征比如近红外和NDVI信息高度重叠留一个就够。互信息特征选择针对的是非线性关系。遥感反演里地表温度、后向散射系数和土壤水分之间不是简单直线关系皮尔逊算出来可能很低但实际有强烈的非线性依赖。互信息不用假设任何函数形式直接衡量一个变量携带另一个变量的信息量。计算时先把连续变量离散化按经验把数值分成8到16个区间区间太少信息损失大太多则互信息估计不稳定。用Python落地时sklearn里有现成实现from sklearn.feature_selection import mutual_info_regression # X 是上面整理好的特征矩阵y 是土壤含水量实测值 # n_neighbors 是连续变量做互信息估计时用的近邻数默认 3 mi_scores mutual_info_regression(X, y, n_neighbors3, random_state42) # 每个特征得到一个互信息分数越高说明对 y 的信息贡献越大 feature_names [B02, B04, B08, VV, VH, dem, ndvi, tvdi, slope] for name, score in zip(feature_names, mi_scores): print(f{name}: {score:.4f}) # 按分数排序保留累计贡献达到 85% 的前 k 个特征 sorted_idx np.argsort(mi_scores)[::-1]参数n_neighbors是互信息估计的敏感参数样本量少时设小一点样本超过500时可以适当调到5近邻数越大估计结果越平滑但也更容易抹掉细节。这里的random_state不是摆设互信息估计里有随机性不固定种子前后两次跑结果不一样后面做重复实验就说不清了。实际项目中还要结合三套筛选结果取交集再按物理经验调整通常保留10到15个特征。3.3 特征保留数量不是越多越好用增量曲线找拐点特征数量对模型的直接影响是过拟合和高维灾难。样本量只有200个时特征加太多模型就有足够自由度去记住每个样本的噪声。业界经验值样本量和特征数的比例至少5比1最好到10比1。但这不是死规矩我实际做的时候用一个增量曲线来定按互信息分数从高到低逐个加入特征看验证集RMSE的变化。典型表现是前几个特征加进去RMSE快速下降加到中间进入平台期偶尔还会略微上升这个平台期的起点就是特征数的合理值。如果一直下降到特征全用完说明要么特征数量还不够要么样本量太大了筛选没有起到约束作用。这个操作配合5折交叉验证来做每一折单独做特征选择再评估得到的曲线才可信。曾经见过直接在全部数据上选特征再交叉验证的结果特征选择和模型评估混在一起高估了反演精度到外场验证时一塌糊涂属于典型的评估泄漏问题。4. GA-BP神经网络遗传算法到底在优化BP的什么GA-BP这个组合容易被误解成用遗传算法替代BP训练实际不是。遗传算法在这里扮演的是“预训练器”的角色它负责找BP的初始参数BP再去精调。理解这个分工后面调参才不至于跑偏。4.1 原理辨析优化的是初始权值和阈值不是网络结构BP网络的训练本质是梯度下降从一组随机的初始权值出发沿误差梯度往下走。随机初始权值的运气成分很大——同一份数据跑十次BP结果可能差出一大截差的初始点会把网络带进局部极小值出不来。遗传算法的全局搜索能力恰好补这块短板先把候选解看成一个个“个体”每个个体是一组完整权值和阈值让它们通过选择、交叉、变异不断进化找到误差较低的个体再把这个个体交给BP做精细调整。GA没有去优化网络结构——隐层几层、每层几个节点这些还是人工定好的。它只优化初始连接参数。以常见的“12个输入特征、8个隐层节点、1个输出”为例总参数数量是输入层到隐层的12×8个权值加隐层到输出层的8个权值再加上隐层8个阈值和输出层1个阈值合计113个。遗传算法编码的长度就是113每个个体是一串113维的实数向量解码后就得到BP的初始权值和阈值。4.2 关键参数表种群规模、迭代次数、交叉与变异概率GA-BP的调参和普通BP完全是两套思路常用参数范围如下参数推荐范围说明种群规模3050太少搜索不充分太多计算量爆炸迭代次数60100前期收敛快后期纠缠在局部解附近交叉概率0.70.9控制个体基因交换频率变异概率0.010.1太小易早熟太大破坏已找到的好解适应度函数用训练集MSE的倒数误差越小适应度越高BP内部训练次数100300GA每评估一个个体就要训练一次BP次数必须控制这里最大的坑是计算量。种群50、迭代100意味着要训练5000次BP每训练一次BP内部还有几百次前向和反向计算数据量大时跑一个晚上都未必能收敛。常见做法是把GA阶段看作粗搜索先用小种群40、迭代60、BP内部训练200次跑起来粗搜结束后再用精调阶段把最优个体投给完整训练的BPBP内部训练次数放到1000以上。用这个两段式思路总耗时能省掉一半以上。4.3 代码骨架GA把权值编码成个体适应度就是BP的MSEGA-BP没有统一的官方实现MATLAB里我常用这段骨架逻辑清晰且可直接改成Python版% 初始化种群每一行是一个个体一个个体就是一组完整权值阈值 % dim_N 113由网络结构计算得出12*8 8*1 8 1 NIND 40; % 种群规模 MAXGEN 60; % 最大进化代数 dim_N 113; % 编码长度 Chrom rand(NIND, dim_N) * 2 - 1; % 初始种群范围 [-1,1] for gen 1:MAXGEN % 对每个个体解码成权值训练一次BP返回MSE for i 1:NIND individual Chrom(i, :); % 解码按网络结构切分成权值矩阵和阈值向量 W1 reshape(individual(1:96), 12, 8); % 12×8 输入到隐层 B1 individual(97:104); % 隐层阈值 W2 reshape(individual(105:112), 8, 1); % 8×1 隐层到输出 B2 individual(113); % 输出层阈值 % 设置BP初始权值训练返回误差 net newff(P_train, T_train, 8, {tansig,purelin}); net.IW{1} W1; net.b{1} B1; net.LW{2} W2; net.b{2} B2; net.trainParam.epochs 200; [net, tr] train(net, P_train, T_train); mse_value tr.perf(end); % 训练结束的MSE ObjV(i) mse_value; % 作为适应度基础 end % 选择误差越小越应该保留 Fitness 1 ./ (ObjV eps); % 交叉、变异工具箱函数 SelCh select(tour, Chrom, Fitness, GGAP); SelCh recombin(xovsp, SelCh, 0.8); SelCh mut(SelCh, 0.05); Chrom reins(Chrom, SelCh, Fitness); end % 最终最优个体解码为BP初始权值重置BP做正式训练这段代码中最关键的是tr.perf(end)这一行它每一步都用真实BP训练结果衡量个体的好坏而不是用某个近似公式估算。这保证了适应度方向没错但也决定了整个方法一定慢所以前面才会强调控制GA阶段的BP训练次数。新ff创建网络时的tansig和purelin组合是隐层双曲正切、输出层线性的搭配适合土壤水分这种数值回归问题输出层不用sigmoid是因为水分含量不只在0到1之间线性输出才能给全范围预测值。4.4 网络结构和归一化隐层节点数、训练函数怎么配GA优化的是权值网络结构还是要自己定。隐含层节点数有一个经验公式h sqrt(m n) am是输入特征数n是输出节点数a是1到10之间的调节值。12个特征、1个输出时sqrt(13)约3.6加调节值后大概在5到13之间。我一般直接遍历5到13每个节点数调一次GA-BP对比验证集RMSE选最优不猜。训练函数的选择直接影响收敛质量。数据量几千行时用trainlmLM算法收敛快梯度信息利用充分但它需要计算近似海森矩阵特征多或数据量大时内存占用会翻好几倍电脑16G内存就有点吃紧。数据量大到上万行时我会换trainbr这个函数自带贝叶斯正则化相当于在误差里加了权值惩罚项能显著抑制过拟合。换训练函数后GA的适应度曲线会变原来的参数可能需要重调别指望一套参数通吃所有数据集。数据归一化这步放在模型设定之前。土壤水分的特征变量量纲差异极大——后向散射系数是浮点小数DEM高程是几百米NDVI在-1到1之间不归一化的话BP的梯度会被大数值特征吃掉。我统一归一化到[-1,1]区间因为隐层用的是tansig对应输出范围正好是[-1,1]附近的变化最灵敏区。一个容易翻车的地方是验证集也必须用训练集的归一化参数映射不能拿验证集自己算一套min-max否则训练集和验证集的数值空间不一致模型相当于换了输入分布测试结果虚高或者崩坏都从这里来。5. 避坑土壤水分反演最容易翻车的五个地方这部分是这些年做遥感反演攒下的血泪经验每一条都是真实发生过的现象写出来供参考。5.1 验证集R²高得离谱空间自相关导致的信息泄漏现象训练集R²达到0.98验证集也高达0.9但一换到邻近时段的新影像精度立刻掉到0.4。原因是数据划分时随机打乱像元而土壤水分存在强烈空间自相关——相邻像元的含水量本来就相似模型等于提前看到答案了。解决方法是按空间地块或影像分块划分训练集和验证集比如将研究区按网格划分成若干子区子区之间不相交训练集选一部分子区验证集选另一部分确保验证集像元在空间上独立。5.2 野外取样和模型输出的土壤含水量单位对不上现象模型预测值在0.3左右野外实测重量含水量却是0.2一套算下来偏差太明显。原因是田间采样常用烘干法得到的是重量含水量而遥感反演模型大多基于体积含水量建立两者相差一个土壤容重系数。解决方法是统一单位将烘干法的重量含水量乘以土壤容重转换为体积含水量再入模。如果实测点没有取容重数据那么模型训练前至少要补测或引用研究区土壤类型对应的容重经验值否则模型和验证数据根本不在同一参考系。5.3 GA参数开太大训练跑到一半想放弃现象种群50、迭代200、BP内部训练次数设成1000跑了一个小时连第一代都没进化完。原因是一个个体就要训练一次BP种群乘以迭代次数就是BP训练的调用次数计算量成指数膨胀。解决方法是分两段跑第一阶段GA阶段种群40、迭代60BP内部训练次数压到200只求找到误差量级正确的好解区域第二阶段把最优个体解码给BPBP内部训练次数放到2000以上精修。实践下来总精度几乎不降总耗时能压缩一半以上。5.4 特征选择把物理上有意义的特征删了精度反降现象按互信息特征选择排序后只保留前八个特征验证集RMSE反而比保留全部特征时更差。原因很典型互信息和相关系数衡量的只是数据里的统计相关性不代表因果机制。有些特征比如DEM高程在与土壤水分的单变量分析中相关性不高但它在模型里和坡度、雷达后向散射组合起来才发挥作用。解决方法是特征选择时做分组先按物理过程把特征分组每组至少保留一个代表再在大类内部筛细节。地形组保留高程、坡度雷达组保留VV、VH光学组保留NDVI和TVDI这样既降维又不断信息链。5.5 TVDI计算时地表温度和光学分辨率不一致现象TVDI特征在模型里贡献很大但个别像元出现异常高值模型预测值也畸变。原因是热红外波段原始分辨率比光学波段低很多比如Landsat热红外是100米可见光近红外是30米直接重采样后算TVDI像元边界处会出现假的高温值。解决方法是先分别完成光学和热红外的各自重采样统一到目标分辨率再计算TVDI不要先算TVDI再重采样那样会把错误比值扩散到邻域。算完后还要做一次异常值统计超出物理范围的值直接掩膜掉。6. 反演结果怎么验证R²、RMSE、MAE之外还缺一个时空泛化测试一套模型跑完如果只报训练集精度就是自欺欺人。土壤水分反演模型落地的底线是三个指标加一个时空泛化验证。指标公式意义判断标准R²决定系数模型解释了多少实测方差反演任务里0.8以上算可用RMSE均方根误差误差的典型幅度数值越小越好按土壤水分单位评估MAE平均绝对误差不受大误差点干扰如果MAE远小于RMSE说明存在少量大误差点NSE纳什效率系数衡量模型是否优于实测均值0.7以上可接受验证方式上除了常见的随机划分交叉验证还要加两个检验时间泛化测试和空间泛化测试。前者用前两年的观测数据训练、后一年的数据验证检验模型跨时间的稳定性后者按子地块留出交叉验证检验模型换位置的适应性。这两个测试都通过模型才算真正具备反演能力。和普通BP对比时有一个值得记住的经验GA-BP在同一份数据上的R²提升通常在0.05到0.15之间这是初始权值优化带来的合理增益。如果提升连0.05都不到先不要怀疑遗传算法的参数设置回头检查特征选择重不充分或者样本量够不够。反过来如果提升超过0.2也先别高兴很有可能是过拟合或验证集划分有问题顺手检查一下训练样本和验证样本有没有空间重叠。我现在的习惯是任何反演模型立项都先跑一版普通BP作为基线特征选择、GA优化每一步改动都和基线比较每一步增益都要能讲清楚来源。玄学式的调参在这个方向上不管用每笔投入都得有对应指标变化撑着。希望这份笔记能帮你少走几段弯路把精力花在真正决定反演精度的数据组织和验证设计上。本文还有配套的精品资源点击获取