13. 迁移U-Net以实现遥感图像变化检测

一、项目背景

1. 应用领域概述：

掌握国土资源利用和土地覆盖类型，是地理国情普查与监测的重要内容。高效获取准确、客观的土地利用情况，监测国土变化情况，可以为国家和地方提供地理国情信息决策支撑。随着遥感、传感器技术的发展，特别是多时相高分辨率遥感图像数据的普及，使我们可以足不出户，就能掌握全球任一地表的细微变化。
目前，我国遥感领域已步入了高分辨率影像的快车道，对遥感数据的分析应用服务的需求也与日俱增。传统方式对高分辨率卫星遥感图像的对特征刻画能力差且依赖人工经验工作量巨大。随着人工智能技术的兴起，特别是基于深度学习的图像识别方法获得了极大的发展，相关技术也推动了遥感领域的变革。相对于传统基于人海战术的目视解译方法，基于深度学习的遥感图像识别技术可以自动分析图像中的地物类型，在准确率和效率方面展现出极大的潜力。

2. 解决的问题：

变化检测——聚焦于建筑的变化检测任务，即给定某地区的双时相遥感图像，要求获得该区域的建筑变化。数据主要来源于北航 LEVIR 团队的公开论文，利用提供的训练数据，实现对多时相图像中的建筑变化检测。具体而言，多时相遥感图像建筑物变化检测任务是给定两张不同时间拍摄的相同位置（地理配准）的遥感图像，要求定位出其中建筑变化的区域。

使用的数据中包括有标注训练样本637对，无标注的测试样本363对（提供测试图像，但不提供对应标签）。每个样本包括前时向遥感图像，后时相遥感图像以及对应的建筑变化标签图。

示例数据——本数据集中的遥感图像大小为1024×1024，像元分辨率为0.5米。标签图与图像等大，像素值为0或255，其中0代表没有建筑变化，255代表存在建筑变化。提供的训练数据集样例如下图所示，每个样本包括前时向图（Image 1），后时相图（Image 2）和标签图（Ground Truth）。

结果输出规范：
输入：给定一对前后时相遥感图像，要求模型能够定位出其中建筑变化的区域。输入图像为RGB三通道的png格式，尺寸为1024×1024像素。
输出：输出图像为单通道png格式的灰度图（8bit），其中255代表变化的建筑区域，0代表未变化建筑区域。输出图像尺寸与输入图相同。

3. 效果：

最终训练数据显示，在已知标签的数据集划分出的预测集上，best_model得到的校验参数分别如下：

miou=0.834496, category_iou=[0.98285981 0.68613132],
oacc=0.983479, category_acc=[0.98986214 0.84118711],
kappa=0.805218, category_F1-score=[0.99135582 0.81385277]；

在无标签的测试数据上，只能得到提交预测结果后系统给出的F1-score信息，F1-score=0.81108。

$$ F1=\frac{TP}{TP+\frac{1}{2(FP+FN)}} $$

某次训练过程曲线如下所示：

多图旋转打分操作效果如下所示：

最终模型的预测效果如下所示：

二、网络架构与主程序流程

1. 网络架构：

本项目采用迁移学习的方法，迁移了卷积神经网络(CNN)中的预训练方式U-Net模型：unet_coco_v3。该语义分割模型在 MSCOCO 验证集上测试得到，大小为 107.3MB。

U-Net 整个流程为U型，左边为下采样过程，右边为上采样过程，中间的灰色箭头是将特征图进行跳层联结，其原理和 DenseNet 相同，即 concatenate，torch.cat([x1,x2])可以将浅层的定位信息和高层的像素分类判定信息进行融合，从而得到更佳的结果。

除此之外，U-Net 有几点值得注意的地方：
1）卷积层 Padding = 0，所以每次做卷积，特征图的大小都会 -2；
2）特征提取的卷积层都为 3×3 大小；
3）下采样使用 max-pooling；
4）上采样使用步长为 2 的反卷积；
5）最后的分类使用 1×1 的卷积层。

2. 主程序流程：

所用的训练环境为国内百度 AI Studio 平台，配置为4核CPU，RAM:32GB，Tesla V100 32GB GPU，磁盘大小：100GB。调参训练20次后得到的最优结果F1-score=0.81108训练时长可控制在6小时以内，滑动窗口参数调小后可达4小时以内，但是分数会有相应减少。

主程序为main.ipynb，需要调用prepare_my_data.py、train.py和predict.py。可以在此处下载，也可以直接运行公开的项目，具体代码如下

三、PaddleX实现代码

1、需要用到的.py文件

1-1. prepare_my_data.py

import os
import os.path as osp
import numpy as np
import cv2
import shutil
import random

import glob
import re

# 为保证每次运行该脚本时划分的样本一致，故固定随机种子
random.seed(0)

import paddlex as pdx

# 定义训练集切分时的滑动窗口大小和步长，格式为(W, H)
train_tile_size = (1024, 1024)
train_stride = (512, 512)
# 定义验证集切分时的滑动窗口大小和步长，格式(W, H)
val_tile_size = (512, 512)
val_stride = (512, 512)
# 训练集和验证集比例
train_ratio = 0.95
val_ratio = 0.05
# 切分后的数据集保存路径
tiled_dataset = './dataset/train'

# 切分后的图像文件保存路径
tiled_image_dir = osp.join(tiled_dataset, 'JPEGImages')
# 切分后的标注文件保存路径
tiled_anno_dir = osp.join(tiled_dataset, 'Annotations')

# 解压数据集
pdx.utils.decompress('/home/aistudio/data/data134796/train_data.zip')
pdx.utils.decompress('/home/aistudio/data/data134796/test_data.zip')
change_det_dataset = '/home/aistudio/data/data134796/train'
image1_dir = osp.join(change_det_dataset, 'A')
image2_dir = osp.join(change_det_dataset, 'B')
label_dir = osp.join(change_det_dataset, 'label')

if not osp.exists(tiled_image_dir):
    os.makedirs(tiled_image_dir)
    print('mkdir' + tiled_image_dir + ' finished!~')
if not osp.exists(tiled_anno_dir):
    os.makedirs(tiled_anno_dir)
    print('mkdir' + tiled_anno_dir + ' finished!~')


# 划分数据集
im1_file_list = os.listdir(image1_dir)
im2_file_list = os.listdir(image2_dir)
label_file_list = os.listdir(label_dir)
im1_file_list = sorted(
    im1_file_list, key=lambda k: int(k.split('_')[-1].split('.')[0]))
im2_file_list = sorted(
    im2_file_list, key=lambda k: int(k.split('_')[-1].split('.')[0]))
label_file_list = sorted(
    label_file_list, key=lambda k: int(k.split('_')[-1].split('.')[0]))

file_list = list()
for im1_file, im2_file, label_file in zip(im1_file_list, im2_file_list,
                                          label_file_list):
    im1_file = osp.join(image1_dir, im1_file)
    im2_file = osp.join(image2_dir, im2_file)
    label_file = osp.join(label_dir, label_file)
    file_list.append((im1_file, im2_file, label_file))
random.shuffle(file_list)
train_num = int(len(file_list) * train_ratio)

# 将大图切分成小图
for i, item in enumerate(file_list):
    if i < train_num:
        stride = train_stride
        tile_size = train_tile_size
    else:
        stride = val_stride
        tile_size = val_tile_size
    set_name = 'train' if i < train_num else 'val'

    # 生成原图的file_list（保存在 data134796/train 里）
    im1_file, im2_file, label_file = item[:]
    mode = 'w' if i in [0, train_num] else 'a'
    with open(
            osp.join(change_det_dataset, '{}_list.txt'.format(set_name)),
            mode) as f:
        f.write("A/{} B/{} label/{}\n".format(
            osp.split(im1_file)[-1],
            osp.split(im2_file)[-1], osp.split(label_file)[-1]))

    # 读取要切分的文件：
    im1 = cv2.imread(im1_file)
    im2 = cv2.imread(im2_file)
    # 将三通道的label图像转换成单通道的png格式图片
    # 且将标注0和255转换成0和1
    label = cv2.imread(label_file, 0)
    label = label != 0
    label = label.astype(np.uint8)

    # 准备切分与保存
    H, W, C = im1.shape
    tile_id = 1
    # 读取原图的文件名列表
    im1_name = osp.split(im1_file)[-1].split('.')[0]
    im2_name = osp.split(im2_file)[-1].split('.')[0]
    label_name = osp.split(label_file)[-1].split('.')[0]
    # 开始切分与保存文件
    for h in range(0, H, stride[1]):
        for w in range(0, W, stride[0]):
            left = w
            upper = h
            right = min(w + tile_size[0], W)
            lower = min(h + tile_size[1], H)
            tile_im1 = im1[upper:lower, left:right, :]
            tile_im2 = im2[upper:lower, left:right, :]
            cv2.imwrite(
                osp.join(tiled_image_dir,
                         "{}_{}_{}.bmp".format('A' ,im1_name, tile_id)), tile_im1)
            cv2.imwrite(
                osp.join(tiled_image_dir,
                         "{}_{}_{}.bmp".format('B', im2_name, tile_id)), tile_im2)
            cut_label = label[upper:lower, left:right]
            cv2.imwrite(
                osp.join(tiled_anno_dir,
                         "{}_{}.png".format(label_name, tile_id)), cut_label)
            mode = 'w' if i in [0, train_num] and tile_id == 1 else 'a'
            with open(
                    osp.join(tiled_dataset, '{}_list.txt'.format(set_name)),
                    mode) as f:
                f.write(
                    "JPEGImages/{}_{}_{}.bmp JPEGImages/{}_{}_{}.bmp Annotations/{}_{}.png\n".
                    format('A', im1_name, tile_id, 'B', im2_name, tile_id, label_name,
                           tile_id))
            tile_id += 1

# 生成labels.txt
label_list = ['unchanged', 'changed']
for i, label in enumerate(label_list):
    mode = 'w' if i == 0 else 'a'
    with open(osp.join(tiled_dataset, 'labels.txt'), 'a') as f:
        name = "{}\n".format(label) if i < len(
            label_list) - 1 else "{}".format(label)
        f.write(name)

1-2. train.py

# 环境变量配置，用于控制是否使用GPU
# 说明文档：https://paddlex.readthedocs.io/zh_CN/develop/appendix/parameters.html#gpu
import os
os.environ['CUDA_VISIBLE_DEVICES'] = '0'

import paddlex as pdx
from paddlex.seg import transforms

# 定义训练和验证时的transforms
# API说明 https://paddlex.readthedocs.io/zh_CN/develop/apis/transforms/seg_transforms.html
train_transforms = transforms.Compose([
    transforms.ResizeStepScaling(
        min_scale_factor=0.5, max_scale_factor=2.,
        scale_step_size=0.25), transforms.RandomRotate(
            rotate_range=45,
            im_padding_value=[127.5] * 6), transforms.RandomPaddingCrop(
                crop_size=512, im_padding_value=[127.5] * 6),
    transforms.RandomDistort(), transforms.RandomHorizontalFlip(),
    transforms.RandomVerticalFlip(), transforms.Normalize(
        mean=[0.5] * 6, std=[0.5] * 6, min_val=[0] * 6, max_val=[255] * 6)
])

eval_transforms = transforms.Compose([
    transforms.Padding(
        target_size=512, im_padding_value=[127.5] * 6), transforms.Normalize(
            mean=[0.5] * 6, std=[0.5] * 6, min_val=[0] * 6, max_val=[255] * 6)
])

# 定义训练和验证所用的数据集
# API说明：https://paddlex.readthedocs.io/zh_CN/develop/apis/datasets.html#paddlex-datasets-changedetdataset
train_dataset = pdx.datasets.ChangeDetDataset(
    data_dir='dataset/train',
    file_list='dataset/train/train_list.txt',
    label_list='dataset/train/labels.txt',
    transforms=train_transforms,
    num_workers=8,
    shuffle=True)

eval_dataset = pdx.datasets.ChangeDetDataset(
    data_dir='dataset/train',
    file_list='dataset/train/val_list.txt',
    label_list='dataset/train/labels.txt',
    num_workers=8,
    transforms=eval_transforms)

# 初始化模型，并进行训练
# 可使用VisualDL查看训练指标，参考https://paddlex.readthedocs.io/zh_CN/develop/train/visualdl.html
num_classes = len(train_dataset.labels)

# API说明：https://paddlex.readthedocs.io/zh_CN/develop/apis/models/semantic_segmentation.html#paddlex-seg-unet
model = pdx.seg.UNet(num_classes=num_classes, input_channel=6)

# 再次训练时加载模型
# str = './output/EzXxY_11_unet/best_model'

# API说明：https://paddlex.readthedocs.io/zh_CN/develop/apis/models/semantic_segmentation.html#train
# 各参数介绍与调整说明：https://paddlex.readthedocs.io/zh_CN/develop/appendix/parameters.html
model.train(
    num_epochs=50,
    train_dataset=train_dataset,
    train_batch_size=2,
    eval_dataset=eval_dataset,
    learning_rate=0.02,
    save_interval_epochs=1,
    pretrain_weights='CITYSCAPES',
    save_dir='output/EzXxY_9_unet',
    use_vdl=True)

1-3. predict.py

# 环境变量配置，用于控制是否使用GPU
# 说明文档：https://paddlex.readthedocs.io/zh_CN/develop/appendix/parameters.html#gpu
import os
os.environ['CUDA_VISIBLE_DEVICES'] = '0'
import cv2
import numpy as np

import paddlex as pdx

model_dir = 'output/EzXxY_9_unet/best_model'
data_dir = 'dataset/test'
file_list = 'dataset/test/test.txt'
save_dir = 'output/EzXxY_9_unet/pred_the_test'

if not os.path.exists(save_dir):
    os.makedirs(save_dir)
# color = [0, 0, 0, 255, 255, 255]

model = pdx.load_model(model_dir)

with open(file_list, 'r') as f:
    for line in f:
        items = line.strip().split()
        img_file_1 = os.path.join(data_dir, items[0])
        img_file_2 = os.path.join(data_dir, items[1])

        # 原图是tiff格式的图片，PaddleX统一使用gdal库读取
        # 因训练数据已经转换成bmp格式，故此处使用opencv读取三通道的tiff图片
        #image1 = transforms.Compose.read_img(img_file_1)
        #image2 = transforms.Compose.read_img(img_file_2)
        image1 = cv2.imread(img_file_1)
        image2 = cv2.imread(img_file_2)
        image = np.concatenate((image1, image2), axis=-1)

        # API说明：https://paddlex.readthedocs.io/zh_CN/develop/apis/models/semantic_segmentation.html#overlap-tile-predict    
        pred = model.overlap_tile_predict(
            img_file=image,
            tile_size=(512, 512),
            pad_size=[256, 256],
            batch_size=2)
        
        # 查看模型预测结果字典
        # 二维图像（0-1之间）
        # print(pred['label_map'])
        # 概率值（0-1之间）
        # print(pred['score_map'].shape)

        # img = pred['label_map']
        # print(img)
        img = np.array(255 * pred['label_map'])
        # img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
        cv2.imwrite(save_dir + '/' + name, img)
        print('save: ' + save_dir + '/' + name + ' finished!~')

        # pdx.seg.visualize(
        #     img_file_1, pred, weight=0., save_dir=save_dir, color=color)

**接下来进入主程序（main.ipynb）**

2、训练模型与初预测

2-1. 安装相关包并进入工作目录

# 约 30 秒
!pip install "paddlex<2.0.0" -i https://mirror.baidu.com/pypi/simple

!pip install --upgrade tqdm

# 进入工作目录
%cd /home/aistudio/work/Change

2-2. 准备数据集并生成训练列表

# 准备数据集（需要大约 2 分钟时间）
!python prepare_my_data.py

2-3. 加载预训练模型开始训练

# 开始训练(已不是 EzXxY_9_unet 的训练过程，此项为小 batch 的模型训练过程)
!python train.py

2-4. 解压测试集并生成测试列表

# 解压测试集
!unzip -o -d ./dataset /home/aistudio/data/data134796/test_data.zip > /dev/null

# 生成测试集列表

# 导入需要的包
import os

# 生成数据列表
def create_list(data_path):
    with open(os.path.join(data_path, 'test.txt'), 'w') as vf:
        for i in range(1, 364):
            img1 = 'A/test_' + '%d'%i + '.png'
            img2 = 'B/test_' + '%d'%i + '.png'
            vf.write(img1 + ' ' + img2 + '\n')

    print('数据列表生成完成')

data_path = './dataset/test'
create_list(data_path)  # 生成数据列表

2-5. 最优kappa系数模型初预测

# 模型预测（预测已有数据，大约需要 18 分钟）
!python predict.py

# 压缩初次预测的测试数据
!zip -r test.zip output/EzXxY_9_unet/pred_the_test

# 压缩训练得到的 best_model 
!zip -r best_model.zip output/EzXxY_9_unet/best_model

3、多图操作打分预测

3-1. 等大旋转生成多图

import os
import glob
import cv2
import numpy as np
import re

# 定义旋转测试集 test 的函数 test_rotate() 
def test_rotate(file_path, file_savepath):
    print(' <rotate the test data> '+file_path+'-->'+file_savepath+':')
    n = 0
    for files in glob.glob(file_path+"/*.png"): #遍历 ./dataset/test/A & B 每个文件
        k = 0
        n += 1
        a = []
        a = re.findall("\d+\.?\d*", files)  #正则表达式提取文件中的数字

        for i in range(7):
            k += 1
            x = cv2.imread(files)

            h, w = x.shape[:2]
            # 第一个参数旋转中心，第二个参数旋转角度，第三个参数：缩放比例, 生成一个 2×3 的矩阵
            center = (w/2, h/2)  # 以图像中心为旋转中心
            angle = 45 * i             # 每次顺时针多转45°
            scale = 1                  # 等比例旋转，即旋转后尺度不变    
            M = cv2.getRotationMatrix2D(center, angle, scale)

            cos = np.abs(M[0, 0])
            sin = np.abs(M[0, 1])
        
            # compute the new bounding dimensions of the image
            nW = int((h * sin) + (w * cos))
            nH = int((h * cos) + (w * sin))
            # adjust the rotation matrix to take into account translation
            M[0, 2] += (nW / 2) - w / 2
            M[1, 2] += (nH / 2) - h / 2

            q = cv2.warpAffine(x, M, (nW, nH))
            cv2.imwrite(file_savepath+'/' + 'test_' + a[0][0:-1] + '_' + '%s' % k + '.png', q)

    print(' new_data_number = ' + '%d'%(n*7))

# 新建文件夹，得到旋转后的测试图像数据（大约 8 分钟）
!mkdir ./dataset/test/optimize
!mkdir ./dataset/test/optimize/A
!mkdir ./dataset/test/optimize/B

# 旋转 A 
file_path = './dataset/test/A'
file_savepath = './dataset/test/optimize/A'
test_rotate(file_path, file_savepath)

# 旋转 B 
file_path = './dataset/test/B'
file_savepath = './dataset/test/optimize/B'
test_rotate(file_path, file_savepath)

3-2. 生成多图测试集列表

# 生成多图测试集列表

# 导入需要的包
import os

# 生成数据列表
def create_list(data_path):
    with open(os.path.join(data_path, 'optimize.txt'), 'w') as vf:
        for i in range(1, 364):
            for j in range(1, 8):
                img1 = 'A/test_' + '%d'%i + '_' + '%d'%j + '.png'
                img2 = 'B/test_' + '%d'%i + '_' + '%d'%j + '.png'
                vf.write(img1 + ' ' + img2 + '\n')

    print('数据列表生成完成')

data_path = './dataset/test'
create_list(data_path)  # 生成数据列表

3-3. 导入模型开始预测

# 环境变量配置，用于控制是否使用GPU
# 说明文档：https://paddlex.readthedocs.io/zh_CN/develop/appendix/parameters.html#gpu
import os
os.environ['CUDA_VISIBLE_DEVICES'] = '0'
import cv2
import numpy as np

import paddlex as pdx

model_dir = 'output/EzXxY_9_unet/best_model'

model = pdx.load_model(model_dir)

data_dir = 'dataset/test/optimize'
file_list = 'dataset/test/optimize.txt'
save_dir = 'output/EzXxY_9_unet/pred_optimize'

if not os.path.exists(save_dir):
    os.makedirs(save_dir)
# color = [0, 0, 0, 255, 255, 255]

with open(file_list, 'r') as f:
    for line in f:
        items = line.strip().split()
        img_file_1 = os.path.join(data_dir, items[0])
        img_file_2 = os.path.join(data_dir, items[1])
        name = img_file_1.split('/')[-1]

        # 原图是tiff格式的图片，PaddleX统一使用gdal库读取
        # 因训练数据已经转换成bmp格式，故此处使用opencv读取三通道的tiff图片
        #image1 = transforms.Compose.read_img(img_file_1)
        #image2 = transforms.Compose.read_img(img_file_2)
        image1 = cv2.imread(img_file_1)
        image2 = cv2.imread(img_file_2)
        image = np.concatenate((image1, image2), axis=-1)

        # API说明：https://paddlex.readthedocs.io/zh_CN/develop/apis/models/semantic_segmentation.html#overlap-tile-predict
        pred = model.overlap_tile_predict(
            img_file=image,
            tile_size=(512, 512),
            pad_size=[256, 256],
            batch_size=2)
        
        # 查看模型预测结果字典
        # 二维图像（0-1之间）
        # print(pred['label_map'])
        # 概率值（0-1之间）
        # print(pred['score_map'].shape)

        # img = pred['label_map']
        # print(img)
        img = np.array(255 * pred['label_map'])
        # img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
        cv2.imwrite(save_dir + '/' + name, img)
        print('save: ' + save_dir + '/' + name + ' finished!~')

        # pdx.seg.visualize(
        #     img_file_1, pred, weight=0., save_dir=save_dir, color=color)

3-4. 等大反旋转预测结果

import os
import glob
import cv2
import numpy as np
import re

# 定义旋转预测标签 pred_optimize 的函数 optimize_rotate() 
def optimize_rotate(file_path, file_savepath):
    print(' <rotate the optimize data> '+file_path+'-->'+file_savepath+':')
    n = 0
    for files in glob.glob(file_path+"/*.png"): #遍历 ./output/.../pred_optimize 每个文件
        k = 0
        n += 1
        a = []
        a = re.findall("\d+\.?\d*", files)  #正则表达式提取文件中的数字
        rora = a[2][0:-1]
        
        k += 1
        x = cv2.imread(files, 0)

        h, w = x.shape[:2]
        # 第一个参数旋转中心，第二个参数旋转角度，第三个参数：缩放比例, 生成一个 2×3 的矩阵
        center = (w/2, h/2)  # 以图像中心为旋转中心
        if rora == '1':
            angle = 0
        elif rora == '2':
            angle = -45
        elif rora == '3':
            angle = -45 * 2
        elif rora == '4':
            angle = -45 * 3
        elif rora == '5':
            angle = -45 * 4
        elif rora == '6':
            angle = -45 * 5
        elif rora == '7':
            angle = -45 * 6

        scale = 1            # 等比例旋转，即旋转后尺度不变    
        M = cv2.getRotationMatrix2D(center, angle, scale)
    
        # compute the new bounding dimensions of the image
        nW, nH = 1024, 1024
        # adjust the rotation matrix to take into account translation
        M[0, 2] += (nW / 2) - w / 2
        M[1, 2] += (nH / 2) - h / 2

        q = cv2.warpAffine(x, M, (nW, nH))
        cv2.imwrite(file_savepath+'/' + 'test_' + a[1] + '_' + rora + '.png', q)

    print(' new_data_number = ' + '%d'%(n))

# 新建文件夹，得到旋转后的预测图像数据（大约 8 分钟）
!mkdir ./output/optimize

# 旋转回去
file_path = './output/EzXxY_9_unet/pred_optimize'
file_savepath = './output/optimize'
optimize_rotate(file_path, file_savepath)

!zip -r optimize.zip ./output/optimize

3-5. 多图打分评价预测结果

# 新建文件夹保存结果
!mkdir ./output/submisson

# 打分评价预测结果
for i in range(1, 364):
    # 获取图像路径
    str1 = './output/optimize/'
    file = str1 + 'test_' + '%d'%i + '_'
    name = './output/submisson/test_' + '%d'%i + '.png'
    path1 = file + '1.png'
    path2 = file + '2.png'
    path3 = file + '3.png'
    path4 = file + '4.png'
    path5 = file + '5.png'
    path6 = file + '6.png'
    path7 = file + '7.png'
    
    img1 = cv2.imread(path1, 0)
    img2 = cv2.imread(path2, 0)
    img3 = cv2.imread(path3, 0)
    img4 = cv2.imread(path4, 0)
    img5 = cv2.imread(path5, 0)
    img6 = cv2.imread(path6, 0)
    img7 = cv2.imread(path7, 0)
    
    img = np.zeros((1024,1024))
    
    for j in range(1024):
        for k in range(1024):
            m = 0;
            if img1[j][k] == 255:
                m += 1
            if img2[j][k] == 255:
                m += 1
            if img3[j][k] == 255:
                m += 1
            if img4[j][k] == 255:
                m += 1
            if img5[j][k] == 255:
                m += 1
            if img6[j][k] == 255:
                m += 1
            if img7[j][k] == 255:
                m += 1
            
            if m >= 4:
                img[j][k] = 255
    
    img = np.array(img, dtype='uint8')
    cv2.imwrite(name, img)
    print('save: ' + name + ' finished!~')

# 保存结果
!zip -r submisson.zip ./output/submisson