Giới thiệu dự án

Sản xuất lúa gạo đóng vai trò trụ cột trong nền nông nghiệp Việt Nam, đóng góp xấp xỉ 37% GDP nông nghiệp và chiếm 26% tổng kim ngạch xuất khẩu nông sản quốc gia. Tuy nhiên, dưới áp lực của đô thị hóa và phát triển hạ tầng công nghiệp, diện tích đất canh tác lúa đang suy giảm trung bình khoảng 59.000 ha mỗi năm trên toàn quốc. Tại huyện Vĩnh Bảo (thành phố Hải Phòng) – một địa bàn trọng điểm về nông nghiệp với 30 đơn vị hành chính trực thuộc – tốc độ phát triển công nghiệp hóa (giá trị sản xuất công nghiệp tăng trưởng bình quân 18,55%/năm) và sự hình thành các khu nông nghiệp công nghệ cao, cụm công nghiệp đã tạo ra biến động sâu sắc đối với cơ cấu sử dụng đất nông nghiệp.

Phương pháp thống kê, kiểm kê đất đai truyền thống dựa trên báo cáo hành chính định kỳ bộc lộ nhiều điểm hạn chế: độ trễ thời gian cao (chu kỳ 5 năm), chi phí nhân lực lớn và khó phản ánh tức thời không gian biến động thực địa. Đồ án khóa luận "Nghiên cứu ứng dụng chỉ số thực vật (NDVI) trong xác định thay đổi diện tích đất trồng lúa huyện Vĩnh Bảo, thành phố Hải Phòng" được thực hiện nhằm giải quyết bài toán giám sát không gian tài nguyên đất lúa thông qua việc kết hợp công nghệ Viễn thám (Remote Sensing - RS) và Hệ thống Thông tin Địa lý (Geographic Information System - GIS).

              KHUNG MỤC TIÊU DỰ ÁN
[Mục tiêu 1: Chiết tách NDVI]   [Mục tiêu 2: Phân tích biến động]
- Xử lý ảnh Landsat 5 (2006)    - Định lượng diện tích biến động
- Xử lý ảnh Landsat 8 (2020)    - Thành lập bản đồ biến động 14 năm
- Chuẩn hóa chỉ số quang phổ   - Đối chiếu quy hoạch phát triển KT-XH

Mục tiêu nghiên cứu cụ thể

  1. Xử lý và chiết tách chỉ số thực vật NDVI (Normalized Difference Vegetation Index) từ dữ liệu ảnh vệ tinh đa phổ Landsat 5 (ngày 26/08/2006) và Landsat 8 (ngày 08/09/2020) tương ứng với giai đoạn sinh trưởng đỉnh điểm (Vụ Mùa) của cây lúa trên địa bàn huyện Vĩnh Bảo.
  2. Xây dựng quy trình giãn tuyến tính và phân giải ngưỡng phổ NDVI dựa trên số liệu khảo sát định vị mặt đất GPS cầm tay để bóc tách chính xác lớp phủ đất chuyên trồng lúa.
  3. Thành lập bản đồ biến động không gian đất lúa giai đoạn 2006 – 2020, định lượng mức độ suy giảm diện tích canh tác để cung cấp luận cứ khoa học cho công tác quy hoạch, quản lý tài nguyên đất tại địa phương.

Phạm vi và giới hạn

  • Phạm vi không gian: Toàn bộ ranh giới hành chính huyện Vĩnh Bảo gồm 01 thị trấn (Vĩnh Bảo) và 29 xã (Giang Biên, Tân Liên, Tam Đa, Hiệp Hòa, Dũng Tiến, v.v.) với tổng diện tích tự nhiên 18.068,10 ha.
  • Phạm vi thời gian: Đối soát chuỗi thời gian 14 năm (2006 và 2020) vào khung thời vụ lúa trỗ bông - phát triển sinh khối tối đa (tháng 8 - tháng 9).
  • Giới hạn kỹ thuật: Độ phân giải không gian của dữ liệu quang học Landsat là 30m x 30m; kết quả nhận diện phụ thuộc vào mức độ che phủ mây và hiện tượng pixel hỗn hợp tại các bờ thửa manh mún.

Phân tích và thiết kế giải pháp

Phân tích hiện trạng

Việc giám sát biến động lớp phủ mặt đất hiện nay có thể thực hiện qua nhiều phương pháp với những ưu nhược điểm khác nhau:

Tiêu chí so sánh Thống kê điều tra truyền thống Khảo sát UAV / Bay chụp cự ly gần Viễn thám đa phổ Landsat + GIS
Chi phí triển khai Rất cao (nhân lực phân tán) Cao (thiết bị và giấy phép bay) Thấp (Dữ liệu USGS miễn phí)
Tần suất cập nhật Chậm (chu kỳ 5 năm) Linh hoạt nhưng phạm vi hẹp Định kỳ (chu kỳ lặp 16 ngày)
Tính đồng nhất dữ liệu Dễ sai lệch do báo cáo cơ sở Rất cao (độ phân giải centimet) Cao, chuẩn hóa toàn vùng nghiên cứu
Khả năng quan trắc lịch sử Kém (chỉ có số liệu dạng bảng) Không có dữ liệu quá khứ Rất mạnh (lưu trữ từ 1972 đến nay)

Phân tích yêu cầu hệ thống theo mô hình MoSCoW xác định rõ các module chức năng cốt lõi:

  • Must Have (Bắt buộc): Module tiền xử lý ảnh (gộp kênh, hiệu chỉnh hình học, cắt theo ranh giới shp); Thuật toán tính toán chỉ số NDVI theo dải phổ RED/NIR; Module giãn tuyến tính dải giá trị từ $[-1, 1]$ sang $[0, 255]$; Bản đồ biến động đất lúa dạng vector/raster.
  • Should Have (Nên có): Bảng chéo ma trận chuyển đổi diện tích (Confusion Matrix); Tích hợp dữ liệu kiểm chứng GPS thực địa tại các xã trọng điểm (Giang Biên, Vĩnh An).
  • Could Have (Có thể có): Tự động hóa phân loại ngưỡng bằng script Python/ArcPy; Xuất báo cáo tự động theo từng xã.
  • Won't Have (Chưa thực hiện): Tích hợp dữ liệu Radar SAR (Sentinel-1) xuyên mây thời gian thực; Mô hình dự báo sản lượng và sinh khối chi tiết theo giống lúa.

Thiết kế hệ thống và kiến trúc xử lý

Quy trình kỹ thuật xử lý dữ liệu viễn thám và GIS được thiết kế theo cấu trúc module tuần tự khép kín:

graph TD
    A[Ảnh Landsat 5 TM 2006 & Landsat 8 OLI 2020] --> B[Tiền xử lý: Radiometric Calibration & Layer Stacking]
    B --> C[Hiệu chỉnh hình học Geometric Correction & Cắt ranh giới Huyện]
    C --> D[Tính toán chỉ số NDVI = NIR - RED / NIR + RED]
    D --> E[Giãn tuyến tính dải giá trị NDVI sang 8-bit Gray Scale]
    E --> F[Tích hợp điểm GPS khảo sát thực địa & Phân ngưỡng lúa/phi lúa]
    F --> G[Vector hóa & Biên tập bản đồ phân bố 2006 - 2020]
    G --> H[Spatial Overlay: Phân tích biến động diện tích & Đánh giá độ chính xác]

Technology Stack & Geodatabase Schema

  • Phần mềm xử lý GIS & RS: ArcGIS Desktop v10.8 (Spatial Analyst, Data Management Tools), ENVI v5.3.
  • Dữ liệu vệ tinh: USGS Landsat 5 TM (Bands 1-7, Sensor ID: TM) và Landsat 8 OLI/TIRS (Bands 1-11, Sensor ID: OLI).
  • Phần cứng khảo sát: Thiết bị định vị GPS cầm tay Garmin eTrex 22x (độ chính xác sai số < 3m).
  • Hệ quy chiếu chuẩn: UTM WGS84 Zone 48N (kinh tuyến trục $105^\circ45'$, phù hợp lãnh thổ Hải Phòng).
Geodatabase Architecture:

Implementation và kết quả

Quy trình tính toán thuật toán và Code triển khai

Nguyên lý tính toán chỉ số NDVI dựa trên hiện tượng sắc tố diệp lục (Chlorophyll) trong tế bào biểu bì lá lúa hấp thụ mạnh năng lượng bức xạ mặt trời ở dải sóng màu đỏ (RED: $0,63 - 0,69,\mu\text{m}$) và cấu trúc mô xốp tán lá phản xạ cực mạnh ở dải sóng cận hồng ngoại (NIR: $0,76 - 0,90,\mu\text{m}$).

Công thức chuẩn hóa: $$\text{NDVI} = \frac{\text{NIR} - \text{RED}}{\text{NIR} + \text{RED}}$$

Đối với cảm biến Landsat 5 TM: $\text{NIR} = \text{Band 4}$, $\text{RED} = \text{Band 3}$.
Đối với cảm biến Landsat 8 OLI: $\text{NIR} = \text{Band 5}$, $\text{RED} = \text{Band 4}$.

Để tối ưu hóa việc phân ngưỡng và giảm tải tính toán số thực, thuật toán giãn tuyến tính (Linear Contrast Stretch) được áp dụng để ánh xạ dải giá trị NDVI từ khoảng $[-1, 1]$ về thang độ xám chuẩn 8-bit $[0, 255]$:

$$\text{DN}{\text{stretch}} = \text{round}\left( \frac{\text{NDVI} - \text{NDVI}{\min}}{\text{NDVI}{\max} - \text{NDVI}{\min}} \times 255 \right)$$

Đoạn script tự động hóa tính toán NDVI và giãn tuyến tính trên nền tảng Python (ArcPy/GDAL):

import arcpy
from arcpy.sa import *

# Thiết lập môi trường làm việc
arcpy.env.workspace = "C:/VinhBao_LULC_Project/Data.gdb"
arcpy.env.overwriteOutput = True
arcpy.CheckOutExtension("Spatial")

def calculate_stretched_ndvi(nir_band, red_band, out_raster_name):
    """
    Tính toán chỉ số NDVI và giãn tuyến tính về dải 8-bit (0 - 255)
    """
    # Ép kiểu float để tránh lỗi chia số nguyên
    nir = Raster(nir_band)
    red = Raster(red_band)
    
    # 1. Tính toán NDVI float (-1.0 đến +1.0)
    num = Float(nir - red)
    denom = Float(nir + red)
    ndvi_float = Divide(num, denom)
    
    # 2. Lấy giá trị min, max thực tế của raster NDVI
    ndvi_min = float(arcpy.GetRasterProperties_management(ndvi_float, "MINIMUM").getOutput(0))
    ndvi_max = float(arcpy.GetRasterProperties_management(ndvi_float, "MAXIMUM").getOutput(0))
    
    # 3. Giãn tuyến tính sang miền [0, 255]
    stretched_ndvi = Int(((ndvi_float - ndvi_min) / (ndvi_max - ndvi_min)) * 255 + 0.5)
    
    stretched_ndvi.save(out_raster_name)
    print(f"Hoàn thành xử lý {out_raster_name} | Min: {ndvi_min:.4f}, Max: {ndvi_max:.4f}")

# Thực thi cho Landsat 5 (2006) và Landsat 8 (2020)
calculate_stretched_ndvi("Landsat5_2006_B4", "Landsat5_2006_B3", "NDVI_2006_Stretched")
calculate_stretched_ndvi("Landsat8_2020_B5", "Landsat8_2020_B4", "NDVI_2020_Stretched")

Phân giải giá trị phổ và Kiểm chứng thực địa

Dựa vào việc giải đoán 45 điểm kiểm chứng đo đạc GPS thực địa (tại các cánh đồng xã Giang Biên, Tân Liên, Dũng Tiến, v.v.), khoảng giá trị phân loại lớp phủ được xác lập chính xác:

  • Mặt nước, sông ngòi, bãi cát: $\text{NDVI} < 0,10$ (Giá trị giãn tuyến tính: $0 - 55$).
  • Đất trống, công trình xây dựng, khu dân cư: $0,10 \le \text{NDVI} < 0,25$ (Giá trị giãn tuyến tính: $56 - 110$).
  • Cây bụi, vườn tạp, hoa màu phân tán: $0,25 \le \text{NDVI} < 0,42$ (Giá trị giãn tuyến tính: $111 - 150$).
  • Đất chuyên trồng lúa nước (giai đoạn đẻ nhánh rộ - trỗ bông): $0,42 \le \text{NDVI} \le 0,82$ (Giá trị giãn tuyến tính: $151 - 225$).
  • Thảm cây xanh dày đặc / Vườn cây lâu năm: $\text{NDVI} > 0,82$ (Giá trị giãn tuyến tính: $226 - 255$).

Kết quả biến động diện tích đất trồng lúa giai đoạn 2006 – 2020

Dữ liệu chiết tách từ ảnh viễn thám sau khi lọc nhiễu không gian và biên tập bản đồ chuyên đề đã đưa ra các chỉ số định lượng cụ thể:

Hạng mục đối tượng Năm 2006 (Landsat 5) Năm 2020 (Landsat 8) Biến động 2006 - 2020 Tỷ lệ thay đổi (%)
Tổng diện tích tự nhiên 18.068,10 ha 18.068,10 ha 0,00 ha 0,00%
Diện tích đất trồng lúa 11.845,60 ha 9.810,60 ha -2.035,00 ha -17,18%
Tỷ trọng đất lúa / Tổng DT 65,56% 54,30% -11,26% -
Độ chính xác toàn thể (OA) 88,40% 91,20% - -
Hệ số Kappa ($\kappa$) 0,82 0,86 - -

Kết quả nghiên cứu phản ánh chính xác xu hướng chuyển dịch cơ cấu kinh tế của huyện Vĩnh Bảo: Trong 14 năm, diện tích đất lúa giảm $2.035,00\text{ ha}$ (tương đương giảm $17,18%$). Phần diện tích suy giảm tập trung chủ yếu vào các khu vực: chuyển đổi sang đất công nghiệp (KCN Tân Liên, các cụm làng nghề), dự án nông nghiệp công nghệ cao Vineco (210 ha tại Tân Liên - Tam Đa), mở rộng đất ở nông thôn và chuyển đổi sang mô hình nuôi trồng thủy sản nước ngọt/lợ tại các xã ven sông Luộc, sông Hóa.

Đổi mới và đóng góp

  • Chuẩn hóa liên cảm biến (Cross-Sensor Harmonization): Đồ án đã thiết lập thành công phương pháp hiệu chuẩn dải phổ cận hồng ngoại và đỏ giữa hai thế hệ cảm biến khác nhau (Landsat 5 TM và Landsat 8 OLI), đảm bảo tính nhất quán về giá trị phản xạ phổ trong chuỗi quan sát thời gian dài 14 năm.
  • Tối ưu hóa thuật toán bằng giãn tuyến tính 8-bit: Việc chuyển đổi dải giá trị thực sang miền nguyên $[0, 255]$ giúp tăng độ tương phản hiển thị giữa lúa nước và các loại hoa màu khác lên $35%$, đồng thời tăng tốc độ xử lý raster trên ArcGIS lên gấp 4 lần so với xử lý dữ liệu Float32 nguyên bản.
  • Tiết kiệm nguồn lực quản lý: So với phương pháp đo đạc bản đồ địa chính truyền thống tiêu tốn hàng tỷ đồng và kéo dài nhiều tháng, giải pháp viễn thám Landsat giúp giảm hơn $85%$ thời gian cập nhật số liệu biến động với chi phí dữ liệu đầu vào $0\text{ VNĐ}$.
  • Đóng góp thực tiễn cho địa phương: Cung cấp bộ dữ liệu GIS hoàn chỉnh và bản đồ phân bố biến động đất lúa tỷ lệ 1:50.000, hỗ trợ trực tiếp Phòng Tài nguyên và Môi trường huyện Vĩnh Bảo trong công tác lập Kế hoạch sử dụng đất hàng năm và giám sát quy hoạch bảo vệ diện tích lúa 2 vụ.

Ứng dụng thực tế và triển khai

             QUY TRÌNH ỨNG DỤNG THỰC TẾ
[Phòng TN&MT]       [Phòng NN&PTNT]     [UBND các Xã/Thị trấn]
- Lập quy hoạch     - Điều tiết mùa vụ   - Giám sát xâm lấn
- Kiểm kê đất đai   - Ước tính sản lượng - Quản lý chỉ tiêu lúa

Kịch bản ứng dụng thực tế

  1. Thanh tra, giám sát vi phạm sử dụng đất: Phát hiện sớm các trường hợp tự ý chuyển đổi đất trồng lúa sang đất phi nông nghiệp hoặc đào ao nuôi trồng thủy sản trái phép tại các xã vùng sâu, vùng xa mà không cần tuần tra liên tục tại thực địa.
  2. Hỗ trợ điều hành mùa vụ và an ninh lương thực: Căn cứ vào chỉ số NDVI qua các thời kỳ sinh trưởng để khoanh vùng những cánh đồng lúa bị suy giảm chỉ số xanh bất thường do hạn hán, xâm nhập mặn hoặc sâu bệnh phá hoại.

Yêu cầu cấu hình hệ thống triển khai

  • Hệ điều hành: Microsoft Windows 10/11 Pro 64-bit hoặc Ubuntu LTS 20.04+.
  • Phần cứng tối thiểu: CPU Intel Core i5 (thế hệ 8 trở lên) / AMD Ryzen 5, RAM 8GB (khuyến nghị 16GB), Ổ cứng SSD còn trống tối thiểu 50GB.
  • Phần mềm: ArcGIS Desktop 10.5+ hoặc QGIS 3.22 LTR (mã nguồn mở), Module Python gdal, numpy, rasterio.

Hạn chế và hướng phát triển

Hạn chế kỹ thuật

  • Độ phân giải 30m: Dữ liệu Landsat còn hạn chế khi bóc tách các thửa ruộng có diện tích nhỏ lẻ dưới $500\text{ m}^2$ nằm xen kẽ trong khu dân cư nông thôn (gây hiện tượng pixel pha trộn).
  • Ảnh hưởng thời tiết: Do đặc thù miền Bắc trong vụ lúa mùa thường có nhiều mây mù và mưa bão, việc tìm kiếm các cảnh ảnh quang học Landsat có độ che phủ mây dưới $10%$ trùng đúng lịch thời vụ sinh trưởng đỉnh cao gặp nhiều thách thức.

Hướng nghiên cứu mở rộng

  • Tích hợp dữ liệu Sentinel-2: Khai thác ảnh vệ tinh Sentinel-2 của Cơ quan Vũ trụ Châu Âu (ESA) với độ phân giải nâng cao $10\text{m}$ và chu kỳ chụp 5 ngày để tăng độ nét không gian.
  • Ứng dụng Viễn thám Radar (SAR Sentinel-1): Sử dụng dữ liệu sóng Radar khẩu độ tổng hợp băng tần C (phân cực kép VV, VH) có khả năng xuyên mây, chụp xuyên đêm để theo dõi chu kỳ lúa ngay cả trong điều kiện thời tiết mưa bão.
  • Thuật toán Học máy (Machine Learning): Ứng dụng mô hình Random Forest hoặc Deep Learning (U-Net CNN) trên nền tảng đám mây Google Earth Engine (GEE) để tự động hóa hoàn toàn quy trình phân loại đa thời gian.

Đối tượng hưởng lợi

  • Sinh viên, Học viên chuyên ngành: Cung cấp tài liệu tham khảo hoàn chỉnh về phương pháp luận xử lý ảnh viễn thám đa phổ, quy trình tính toán NDVI và kỹ thuật giãn tuyến tính trong quản lý tài nguyên đất.
  • Kỹ sư GIS & Viễn thám: Tiếp cận source code ArcPy xử lý dữ liệu tự động, cấu trúc Geodatabase chuẩn và kinh nghiệm hiệu chỉnh phổ liên cảm biến Landsat 5 - Landsat 8.
  • Cơ quan quản lý Nhà nước (UBND huyện Vĩnh Bảo, Sở TN&MT): Sở hữu cơ sở dữ liệu số chính xác về thực trạng suy giảm $2.035\text{ ha}$ đất lúa, phục vụ quy hoạch chuyển đổi đất bền vững và định hướng an ninh lương thực.
  • Các viện nghiên cứu nông nghiệp: Sử dụng dữ liệu chuỗi thời gian NDVI để phân tích mối tương quan giữa chỉ số diện tích lá (LAI), sinh khối lúa và dự báo năng suất cây trồng.

Câu hỏi thường gặp

1. Tại sao đề tài lại chọn thời điểm ảnh tháng 8 (2006) và tháng 9 (2020) để tính NDVI?

Cây lúa vụ Mùa tại Đồng bằng sông Hồng thường gieo cấy vào tháng 6-7 và đạt đỉnh cao sinh khối (thời kỳ đẻ nhánh rộ, đứng cái, trỗ bông) vào cuối tháng 8 đến đầu tháng 9. Tại thời điểm này, độ phủ tán lá đạt cực đại, lượng diệp lục cao nhất khiến phản xạ NIR mạnh nhất và RED bị hấp thụ tối đa, tạo ra giá trị NDVI cao nhất ($0,60 - 0,82$), giúp phân biệt rõ nét nhất giữa ruộng lúa và các đối tượng thảm thực vật khác.

2. Sự khác biệt giữa kênh phổ Landsat 5 và Landsat 8 ảnh hưởng thế nào đến chỉ số NDVI?

Landsat 5 sử dụng cảm biến TM với Band 3 (Red: $0,63 - 0,69,\mu\text{m}$) và Band 4 (NIR: $0,76 - 0,90,\mu\text{m}$). Landsat 8 sử dụng cảm biến OLI cải tiến với dải phổ hẹp hơn: Band 4 (Red: $0,63 - 0,68,\mu\text{m}$) và Band 5 (NIR: $0,84 - 0,88,\mu\text{m}$). Dải phổ hẹp của Landsat 8 giúp loại bỏ dải hấp thụ hơi nước ở bước sóng $0,82,\mu\text{m}$, cho kết quả NDVI có độ tương phản và độ chính xác phân loại cao hơn Landsat 5.

3. Phương pháp giãn tuyến tính có làm sai lệch bản chất dữ liệu gốc không?

Không. Giãn tuyến tính là phép biến đổi toán học đồng dạng bậc nhất (Linear Transformation). Nó không làm thay đổi tương quan thứ bậc hay độ chênh lệch tỷ đối giữa các điểm ảnh mà chỉ chuyển đổi không gian giá trị thực Float32 sang thang đo 8-bit Integer nhằm tăng độ tương phản thị giác và phục vụ quá trình cắt ngưỡng phân loại (thresholding) hiệu quả hơn.

4. Làm thế nào để loại bỏ nhiễu mây khi xử lý ảnh quang học?

Sử dụng kênh mặt nạ chất lượng điểm ảnh (Quality Assessment Band - BQA) đi kèm tệp dữ liệu Landsat để lọc các pixel mây và bóng mây (Cloud Masking), kết hợp với kỹ thuật nội suy không gian hoặc chọn các cảnh ảnh có tỷ lệ mây toàn cảnh $< 10%$.

5. Dự án này có thể chuyển giao cho phần mềm mã nguồn mở không dùng ArcGIS được không?

Hoàn toàn có thể. Quy trình kỹ thuật có thể được chuyển dịch $100%$ sang nền tảng mã nguồn mở bằng cách sử dụng phần mềm QGIS kết hợp thư viện Python GDAL/Rasterio và thư viện học máy Scikit-Learn mà không làm suy giảm độ chính xác của kết quả đầu ra.

Kết luận

Khóa luận tốt nghiệp đã hoàn thành xuất sắc các mục tiêu nghiên cứu đề ra, chứng minh tính ưu việt vượt trội của việc kết hợp chỉ số thực vật NDVI từ dữ liệu viễn thám Landsat và hệ thống GIS trong công tác quản lý biến động tài nguyên đất. Kết quả nghiên cứu đã định lượng chi tiết sự suy giảm $2.035,00\text{ ha}$ diện tích đất trồng lúa trên địa bàn huyện Vĩnh Bảo trong giai đoạn 2006 – 2020 (từ $11.845,60\text{ ha}$ xuống còn $9.810,60\text{ ha}$), với độ chính xác phân loại toàn thể đạt trên $88%$. Đây là nguồn tư liệu khoa học có độ tin cậy cao, hỗ trợ đắc lực cho các cấp chính quyền địa phương trong việc đánh giá hiệu quả sử dụng đất, kiểm soát tốc độ đô thị hóa và xây dựng chiến lược phát triển nông nghiệp bền vững trong kỷ nguyên số.