论文精读：基于深度学习和年度结果增强的东北水稻长时序制图（1985-2023）

基于深度学习和年度结果增强的东北水稻长时序制图（1985-2023）

zhang, zihui; Xia, Lang; Zhao, Fen; Gu, Yue; Yang, Jing; Zha, Yan; et al. (2024). Long history paddy rice mapping across Northeast China with deep learning and annual result enhancement method. figshare. Dataset. https://doi.org/10.6084/m9.figshare.27604839.v1
张子惠; 夏朗; 赵芬; 古越; 杨静; 查岩; 等. (2024). 基于深度学习和年度结果增强方法的长时段中国东北地区水稻制图. figshare. 数据集. https://doi.org/10.6084/m9.figshare.27604839.v1

一、研究背景与核心挑战

水稻是全球半数以上人口的主粮，准确掌握其长时序空间分布对粮食安全评估、农业政策制定和温室气体排放核算至关重要。东北三省及内蒙古东部是中国重要的水稻生产基地，近几十年来种植面积急剧扩张，但高精度、年际一致的长时序水稻分布图仍然缺乏。

现有方法的痛点：

• 传统机器学习（如XGBoost/SVM）：依赖人工特征提取，在复杂种植系统下区分水稻与其他作物能力有限
• 物候法：基于光谱指数（如NDVI）阈值判断，受云雨影响易产生时间序列中断
• 深度学习中的RNN/LSTM：依赖连续无云时序数据，而Landsat在东北地区云覆盖严重，很难获得完整物候曲线
• 语义分割模型：输入灵活，但年度制图时通常采用“多时相叠加法”（任意时相检测为水稻则最终为水稻），忽略了不同物候期水稻光谱/纹理特征的巨大差异以及模型置信度的差异，导致误差累积

本文核心创新：

1. 跨传感器训练数据集：构建包含Landsat 5 TM和Landsat 8/9 OLI共155景的高质量水稻标签，首次实现跨传感器长时序（1985–2023）水稻制图
2. FRNet（全分辨率网络）：一种轻量级语义分割网络，通过多分辨率特征融合保留空间细节，适合处理Landsat图像
3. ARE（年度结果增强方法）：利用模型输出的类别概率，自动选择年内各时相中置信度最高（即|P-0.5|最大）的时相结果作为年度最终图，有效抑制误差传播

创新点

创新点1：跨传感器训练数据集构建策略
传统问题：Landsat 5 TM和Landsat 8 OLI光谱响应函数不同，直接用TM训练的模型预测OLI影像，F1会从0.85暴跌至0.48。以往研究要么回避长时序跨传感器（只用单一传感器），要么用直方图匹配等简单校正，效果有限。

本文创新：

主动构建包含两种传感器影像的混合训练集：从TM和OLI中各选取代表性场景，通过XGBoost初步分类+目视校正，得到155景高质量标签。

系统评估9种训练/测试组合（表4），证明：混合训练（TM+OLI）测试混合数据时F1=0.84，接近同传感器精度；而仅用单传感器训练+微调的方法效果差（F1≈0.67–0.70）。

核心结论：长时序跨传感器制图，必须将不同传感器的样本在训练阶段就混合，而非靠后处理校正。这为所有长时序遥感产品提供了方法论参考。

创新点2：FRNet中的MRFU（多分辨率特征融合单元）
传统语义分割网络（如U-Net, DeepLab）通过编码器-解码器结构逐步降采样再上采样，导致空间细节损失（边缘模糊、小地块丢失）。虽然U-Net有跳跃连接，但深层特征仍丢失大量细节。

本文创新：

MRFU设计：每个单元同时维持两条流：

水平流：保持原始分辨率，用1×1卷积调整通道，保留空间细节。

垂直流：步长2卷积降采样（捕捉更大感受野的语义），通道数加倍，再通过转置卷积恢复到原分辨率。

两条流相加融合，使得输出特征图始终保留高分辨率，避免了下采样的不可逆损失。

轻量级：相比DeepLabV3+的ASPP（空洞空间金字塔池化）和U-Net的多层上采样，MRFU参数量小，计算快，适合Landsat这种中等数据量场景。

创新点3：ARE（年度结果增强）方法 —— 基于置信度最大化的时相选择
传统做法：一年内有多张影像时，采用“叠加法”（OR）——只要任一时相被分类为水稻，最终就标为水稻。问题在于：某些时相（如幼苗期）模型容易误判，错误被永久保留。

本文创新：

模型输出每个像元的类别概率P（sigmoid输出，0~1）。

定义置信度 = |P - 0.5|，该值越大表示模型越确定（接近0或1）。

对一年内的m个时相，选出置信度最大的时相t = argmax(|P_i - 0.5|)。

最终结果 = (P_t ≥ 0.5 ? 水稻 : 非水稻)。

为什么有效：

水稻在不同物候期的光谱可分性不同（抽穗期最典型，幼苗期易混）。模型在抽穗期输出的概率往往接近0.9–0.99，置信度≈0.4–0.49；在幼苗期可能只有0.6，置信度≈0.1。ARE自动选择抽穗期结果。

相比直接取平均或投票，ARE保留了单时相的最佳判别能力，同时抑制了低质量预测。

定量提升：相比叠加法，OA提升6%，F1提升5%，MCC提升13%（MCC是更严格的指标，对不平衡数据敏感）。

创新点4：长时序缺失像素填补策略
问题：Landsat历史数据中，云污染导致某些年份、某些区域在水稻生长季内没有任何一景无云影像。缺失像素不能直接留空。

本文方法（虽未在方法部分大书特书，但在结果中说明）：

对于缺失年，用前后年份的“良好观测”影像进行时间插补。具体来说，基于前后两年中同一位置的水稻类别（假设水稻种植稳定性较高），并结合邻近像元的空间信息，填补缺失值。

同时生成年际云覆盖掩模图，标注每个像素的观测来源年份，让使用者了解数据可靠性。

创新点：这是首个在东北水稻长时序制图中明确处理“多年无云影像缺失”并给出填补方法和不确定性指示的产品。

二、核心技术方法详解

2.1 数据获取与预处理

• 卫星数据：美国地质调查局（USGS）提供的Landsat Collection 2 Level-2地表反射率产品，共13,809景，覆盖1985–2023年。使用Landsat 5 TM（1985–2011）和Landsat 8/9 OLI（2013–2023），Landsat 7因条带噪声被排除。
• 波段选择：蓝、绿、红、近红外（NIR）、短波红外1（SWIR1）、短波红外2（SWIR2）共6个波段，30 m分辨率。
• 影像筛选：选取每年水稻生长季（5–9月）内云量较少的影像，每个地区每年尽量获取多时相影像。

2.2 跨传感器训练数据集构建

由于Landsat 5 TM和Landsat 8/9 OLI在光谱响应范围和辐射定标上存在差异，直接用单传感器模型预测另一传感器影像会出现精度骤降。

构建流程：

1. 从Landsat影像中人工选取水稻/非水稻的感兴趣区（ROI）
2. 用XGBoost模型初步分类，获得初始水稻分布图
3. 通过目视解译逐景人工校正，得到高精度标签
4. 共获得155景覆盖不同年份、不同季节的标签影像（图1c,d）
5. 按照3:1比例划分为训练集和验证集

跨传感器训练策略：作者对比了9种训练/测试组合（表4），发现：

• 使用单一传感器训练，测试另一传感器时F1仅0.48–0.53（组合2,3）
• 使用混合传感器（TM+OLI）共同训练，测试混合数据时F1=0.84（组合7）
• 结论：长时序跨传感器制图必须建立包含所有传感器样本的混合训练集，迁移学习效果有限

2.3 FRNet：全分辨率语义分割网络

FRNet是一种专门针对遥感图像设计的语义分割网络，核心目标是在深层网络中保留空间细节，避免传统编码-解码结构导致边缘模糊。

核心模块：MRFU（多分辨率特征融合单元）
MRFU包含两个并行分支：

• 水平流：保持原始空间分辨率，通过恒等映射保留细节
• 垂直流：空间分辨率减半（步长2卷积），通道数加倍，捕获更全局的语义信息

两个分支通过3×3卷积、批归一化（BN）和ReLU激活后融合。

网络优势：

• 结构简单，参数量少，训练速度快
• 避免梯度消失问题（无深层编码器-解码器级联）
• 输出特征图始终保持较高分辨率，适合精细农业地块边界分割

损失函数：Dice Loss，天然处理正负样本不平衡（水稻像素通常远少于非水稻）
Dice Loss = 1 − 2 ∣ A ∩ B ∣ ∣ A ∣ + ∣ B ∣ \text{Dice Loss} = 1 - \frac{2|A \cap B|}{|A|+|B|}Dice Loss=1−∣A∣+∣B∣2∣A∩B∣​

训练设置：

• 框架：Keras 2.5 + TensorFlow 2.6
• GPU：NVIDIA RTX 3090
• 优化器：Adam，学习率0.001，批大小8
• 图像切块：256×256，无重叠
• 数据增强：随机旋转5°

2.4 ARE（年度结果增强方法）：核心创新点

问题定义：
一年内同一地区可能有多个时相的Landsat影像（例如6月、7月、8月）。常规做法（式1）是：只要任一时相检测为水稻，年度结果就标为水稻（OR叠加法）。这会导致：

• 某个时相因云影、物候未完全发育等因素产生误判，错误的水稻标签被永久保留
• 不同时相模型置信度差异被忽略

ARE方法原理（式2）：
对于一年内的m张影像，模型输出每个像元属于“水稻”类别的概率P_i（范围0~1）。定义置信度C_i = |P_i - 0.5|，该值越大表示模型越确定（接近0或1）。ARE方法找出置信度最大的时相t：
t = arg ⁡ max ⁡ i ( ∣ P i − 0.5 ∣ ) , i ∈ [ 1 , m ] t = \arg\max_i \big( |P_i - 0.5| \big),\quad i \in [1, m]t=argimax​(∣Pi​−0.5∣),i∈[1,m]
然后将该时相的类别概率阈值化得到最终结果：
Result ( x , y ) = { paddy if P t ( x , y ) ≥ 0.5 non-paddy otherwise \text{Result}(x,y) = \begin{cases} \text{paddy} & \text{if } P_t(x,y) \ge 0.5 \\ \text{non-paddy} & \text{otherwise} \end{cases}Result(x,y)={paddynon-paddy​ifPt​(x,y)≥0.5otherwise​

直观理解：选择模型“最有把握”的那个物候期来做决策。例如，水稻在抽穗期（7–8月）光谱特征最典型，模型概率往往接近1；而在幼苗期（5–6月）可能混有其他湿生植被，概率在0.6左右。ARE会自动选择高置信度时相，避免低质量预测污染最终图。

效果提升：

• 与传统叠加法相比，OA提高6%，F1提高5%，MCC提高13%（表2）
• 显著消除椒盐噪声，地块边界更清晰（图4）

三、核心代码实现示例

3.1 FRNet的MRFU模块（Keras/TensorFlow）

importtensorflowastffromtensorflow.kerasimportlayers,ModeldefMRFU(input_tensor,filters):""" 多分辨率特征融合单元 input_tensor: 输入特征图 (H, W, C) filters: 输出通道数 """# 水平流：保持分辨率，1x1卷积调整通道horizontal=layers.Conv2D(filters,1,padding='same')(input_tensor)horizontal=layers.BatchNormalization()(horizontal)horizontal=layers.ReLU()(horizontal)# 垂直流：步长2下采样，通道数加倍vertical=layers.Conv2D(filters*2,3,strides=2,padding='same')(input_tensor)vertical=layers.BatchNormalization()(vertical)vertical=layers.ReLU()(vertical)# 上采样恢复分辨率vertical=layers.Conv2DTranspose(filters,3,strides=2,padding='same')(vertical)vertical=layers.BatchNormalization()(vertical)vertical=layers.ReLU()(vertical)# 融合output=layers.Add()([horizontal,vertical])returnoutputdefFRNet(input_shape=(256,256,6)):""" 简化的FRNet结构，包含3个MRFU和输出层 """inputs=layers.Input(shape=input_shape)# 初始卷积x=layers.Conv2D(32,3,padding='same')(inputs)x=layers.BatchNormalization()(x)x=layers.ReLU()(x)# 3个MRFU模块x=MRFU(x,64)x=MRFU(x,128)x=MRFU(x,256)# 输出层：1通道二分类x=layers.Conv2D(1,1,padding='same',activation='sigmoid')(x)model=Model(inputs,x)returnmodel# 编译model=FRNet()model.compile(optimizer=tf.keras.optimizers.Adam(0.001),loss=tf.keras.losses.BinaryCrossentropy(),# 配合sigmoid输出metrics=['accuracy'])

3.2 ARE方法实现

importnumpyasnpdefapply_are(prob_maps):""" 年度结果增强（ARE）方法 prob_maps: list of 2D numpy arrays，每个元素是模型输出的水稻概率图 (H, W)，值域[0,1] return: 最终二值水稻图 (H, W)，1表示水稻，0表示非水稻 """# 堆叠概率图: (m, H, W)prob_stack=np.stack(prob_maps,axis=0)# 计算每个像素每个时相的置信度 |P - 0.5|confidence=np.abs(prob_stack-0.5)# (m, H, W)# 选出置信度最高的时相索引best_idx=np.argmax(confidence,axis=0)# (H, W)# 根据最佳时相的概率决定类别# 构建与prob_maps形状相同的索引数组h,w=prob_maps[0].shape result=np.zeros((h,w),dtype=np.uint8)foriinrange(h):forjinrange(w):best_t=best_idx[i,j]ifprob_stack[best_t,i,j]>=0.5:result[i,j]=1returnresult# 示例：一年内三个时相的预测概率prob_june=np.random.rand(512,512)# 6月prob_july=np.random.rand(512,512)# 7月prob_aug=np.random.rand(512,512)# 8月final_rice=apply_are([prob_june,prob_july,prob_aug])

3.3 跨传感器训练数据准备（关键步骤）

importrasterioimportnumpyasnpfromsklearn.model_selectionimporttrain_test_splitdefprepare_cross_sensor_training(landsat5_paths,landsat8_paths,label_paths):""" 将Landsat5和Landsat8影像及其标签统一处理为训练样本 假设所有影像已重采样至30m并配准 """X_all=[]y_all=[]# 处理Landsat5影像forimg_path,lbl_pathinzip(landsat5_paths,label_paths):withrasterio.open(img_path)assrc:img=src.read([1,2,3,4,5,6])# 蓝绿红NIR SWIR1 SWIR2img=np.transpose(img,(1,2,0))# (H,W,6)withrasterio.open(lbl_path)assrc_lbl:lbl=src_lbl.read(1)# 单波段标签，1=水稻，0=非水稻# 切块256x256h,w,_=img.shapeforiinrange(0,h,256):forjinrange(0,w,256):img_patch=img[i:i+256,j:j+256]lbl_patch=lbl[i:i+256,j:j+256]ifimg_patch.shape[0]==256andimg_patch.shape[1]==256:X_all.append(img_patch)y_all.append(lbl_patch)# 同样处理Landsat8影像（可合并）# ...X_all=np.array(X_all,dtype=np.float32)y_all=np.array(y_all,dtype=np.float32)# 归一化反射率（0-1）X_all=X_all/10000.0# Landsat反射率缩放因子# 划分训练/验证X_train,X_val,y_train,y_val=train_test_split(X_all,y_all,test_size=0.25,random_state=42)returnX_train,X_val,y_train,y_val

3.4 完整训练流程

# 加载跨传感器数据train_imgs,val_imgs,train_masks,val_masks=prepare_cross_sensor_training(...)# 构建模型model=FRNet(input_shape=(256,256,6))model.compile(optimizer='adam',loss='binary_crossentropy',metrics=['accuracy'])# 训练（使用Dice Loss替代二分类交叉熵可选）fromtensorflow.kerasimportbackendasKdefdice_loss(y_true,y_pred):smooth=1e-6intersection=K.sum(y_true*y_pred)return1-(2.*intersection+smooth)/(K.sum(y_true)+K.sum(y_pred)+smooth)model.compile(optimizer='adam',loss=dice_loss)# 训练history=model.fit(train_imgs,train_masks,validation_data=(val_imgs,val_masks),batch_size=8,epochs=50,callbacks=[tf.keras.callbacks.ModelCheckpoint('frnet_best.h5',save_best_only=True)])

四、关键实验结果

4.1 ARE方法提升显著

ARE相比叠加法，MCC从0.68提升至0.81（提升13%），说明消除了大量错误预测。

4.2 最终数据集精度（表3）

• 水稻的用户精度（UA）：0.93
• 水稻的生产者精度（PA）：0.91
• 总体精度（OA）：0.91
• F1得分：0.92
• Matthews相关系数（MCC）：0.82

4.3 跨传感器训练必要性（表4）

结论：单一传感器模型无法泛化到另一传感器，必须使用跨传感器混合训练集。

4.4 长时序变化（1985–2023）

• 总种植面积从1.11万km²扩张至6.45万km²，增长4.81倍（净增5.34万km²）
• 黑龙江省贡献最大（+4.33万km²），辽宁省增长最缓慢
• 扩张主要发生在三江平原、松嫩平原等区域，呈“北扩东进”趋势

五、局限性与未来方向

当前局限：

1. 混合像元影响：30m分辨率下稻田边界存在混合像元，虽用面积比例法计算混淆矩阵，但仍有不确定性
2. 单张影像时年：如果某年某地仅有一景无云影像，ARE退化为单时相，精度提升有限
3. 早期验证数据不足：1985–2001年缺乏高分辨率影像和实地数据，仅用农业统计数据进行面积验证

未来方向：

1. 结合Sentinel-1 SAR数据填补云覆盖严重区域
2. 引入自监督/半监督学习减少对人工标注的依赖
3. 将ARE扩展为“概率加权融合”，而非硬选最大置信度
4. 构建全国乃至全球尺度的长时序水稻数据集

六、数据与代码获取

• 水稻分布图（1985–2023，30m）：Figshare https://doi.org/10.6084/m9.figshare.27604839.v1（GeoTIFF格式，1=水稻，0=非水稻，EPSG:4326）
• FRNet代码：Zenodo https://doi.org/10.5281/zenodo.17744587
• 训练标签数据集：Figshare https://doi.org/10.6084/m9.figshare.28283606
• 原始Landsat数据：USGS EarthExplorer https://earthexplorer.usgs.gov/

七、总结

本文首次实现了中国东北地区1985–2023年连续39年、30m空间分辨率的水稻年际动态制图。核心贡献在于：

1. 构建跨传感器（Landsat 5 TM + Landsat 8/9 OLI）的训练样本集，解决了不同传感器特征分布差异导致的模型失效问题
2. 提出ARE方法，利用模型输出的类别概率自动选择年内最优时相，比传统叠加法F1提升5%
3. 揭示了东北水稻种植面积近40年扩张4.8倍的时空演变规律，为农业政策制定和温室气体清单提供关键基础数据

该方法可推广至其他长时序作物制图任务（如玉米、大豆、小麦），只需替换训练标签并调整物候参数即可。