Giới thiệu dự án

Quá trình đô thị hóa với tốc độ nhanh tại Việt Nam trong ba thập kỷ qua đang tạo ra áp lực to lớn lên công tác quản lý tài nguyên đất đai và quy hoạch không gian đô thị. Tại thành phố Hà Nội, quận Long Biên là địa bàn có tốc độ phát triển hạ tầng và chuyển dịch cơ cấu sử dụng đất diễn ra mạnh mẽ bậc nhất ở khu vực cửa ngõ phía Đông Bắc. Việc chuyển đổi ồ ạt đất nông nghiệp sang đất phi nông nghiệp cùng sự xuất hiện của các dự án chậm triển khai đã làm gia tăng đáng kể diện tích đất xây dựng xen lẫn các khu đất trống bỏ hoang.

                  +----------------------------------------------+
                  |   Dữ liệu Viễn thám Đa phổ (Landsat 8/9 OLI) |
                  +----------------------------------------------+
                                         |
                                         v
                  +----------------------------------------------+
                  |     Tiền xử lý & Hiệu chỉnh Bức xạ/Hình học  |
                  +----------------------------------------------+
                                         |
                   +---------------------+---------------------+
                   |                                           |
                   v                                           v
+-------------------------------------+     +-------------------------------------+
|    Trích xuất Chỉ số Quang phổ      |     |     Phân loại Giám định (MLC/SVM)   |
|         (NDBI, UI, EBBI)            |     |        với Tập mẫu Thực địa         |
+-------------------------------------+     +-------------------------------------+
                   |                                           |
                   +---------------------+---------------------+
                                         |
                                         v
                  +----------------------------------------------+
                  |   Đánh giá Ma trận Sai số & Chỉ số Kappa (κ) |
                  |    --> Chiết tách Đất xây dựng & Đất trống   |
                  +----------------------------------------------+

Vấn đề thực tiễn và bài toán đặt ra

Phương pháp trắc địa truyền thống bằng máy toàn đạc điện tử hoặc đo vẽ bản đồ thực địa gặp nhiều khó khăn tại đô thị có mật độ dân cư cao, nhiều công trình cao tầng che khuất và diện tích trải rộng (trên 60 km²). Hạn chế lớn nhất là chi phí nhân công cao, chu kỳ cập nhật chậm (3–5 năm) và không bắt kịp biến động nhanh của bề mặt lớp phủ.

Ứng dụng ảnh viễn thám vệ tinh mở ra hướng tiếp cận tự động hóa hiệu quả cao. Tuy nhiên, việc phân biệt giữa đất xây dựng (built-up land/impervious surface) và đất trống (bare land) trên ảnh quang học thường xảy ra hiện tượng nhầm lẫn phổ nghiêm trọng. Nguyên nhân do đặc trưng phản xạ phổ của bê tông, mái tôn, đường nhựa và đất trống khô có quy luật biến thiên tương tự nhau ở dải sóng hồng ngoại cận (NIR) và hồng ngoại sóng ngắn (SWIR).

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

  1. Hệ thống hóa cơ sở khoa học về xử lý ảnh số viễn thám trong giám sát bề mặt lớp phủ đô thị.
  2. Thu thập, tiền xử lý và chuẩn hóa dữ liệu ảnh vệ tinh Landsat 8/9 OLI-TIRS (chu kỳ chụp tháng 05/2022) tại khu vực quận Long Biên, Hà Nội.
  3. Ứng dụng và so sánh hiệu quả chiết tách đất xây dựng và đất trống bằng 2 phương pháp: Phân ngưỡng ảnh chỉ số quang phổ (NDBI, UI, EBBI) và Phân loại ảnh có giám định (Supervised Classification).
  4. Kiểm chứng độ chính xác bằng ma trận sai số (Confusion Matrix), độ chính xác toàn cục (Overall Accuracy - OA) và chỉ số Kappa ($\kappa$) dựa trên tập mẫu điều tra thực địa.
  5. Xây dựng bản đồ chuyên đề phân bố đất xây dựng và đất trống quận Long Biên tỷ lệ 1/25.000 phục vụ hỗ trợ ra quyết định quản lý quy hoạch.

Phạm vi và giới hạn đề tài

  • Không gian: Toàn bộ ranh giới hành chính 14 phường thuộc quận Long Biên, TP. Hà Nội với tổng diện tích tự nhiên $60.38\text{ km}^2$.
  • Thời gian & Dữ liệu: Cảnh ảnh Landsat 8/9 thu nhận ngày 17/05/2022, độ phân giải không gian $30\text{m}$ (kênh đa phổ) kết hợp kênh nhiệt TIRS ($100\text{m}$).
  • Giới hạn kỹ thuật: Các điểm ảnh hỗn hợp (mixed pixels) tại khu vực dân cư mật độ cao xen lẫn cây xanh tán nhỏ dưới $30\text{m}\times 30\text{m}$ chịu ảnh hưởng nhất định bởi độ phân giải hình học của sensor Landsat.

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

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

Tiêu chí so sánh Đo đạc trắc địa truyền thống Phân loại có giám định (MLC/SVM) Phương pháp Phân ngưỡng Chỉ số (NDBI/UI/EBBI)
Thời gian thực hiện 3 - 6 tháng cho toàn quận 3 - 5 ngày xử lý số 1 - 2 ngày xử lý bán tự động
Chi phí triển khai Rất cao (nhân lực, thiết bị đo) Thấp (tận dụng ảnh viễn thám mở) Rất thấp (thuật toán tính toán trực tiếp)
Khả năng tự động hóa Thủ công hoàn toàn Trung bình (cần chọn mẫu tỉ mỉ) Cao (dễ lập trình hàng loạt)
Độ nhầm lẫn phổ Không có sai số phổ Dễ nhầm lẫn pixel hỗn hợp Khắc phục tốt khi kết hợp dải nhiệt
Chu kỳ cập nhật 3 - 5 năm/lần 8 - 16 ngày (theo chu kỳ vệ tinh) 8 - 16 ngày (theo chu kỳ vệ tinh)

Yêu cầu hệ thống theo mô hình MoSCoW

  • Must have: Khả năng tiền xử lý bức xạ (chuyển đổi DN sang TOA Radiance/Reflectance), tính toán chính xác 3 chỉ số NDBI, UI, EBBI, phân loại có kiểm định Maximum Likelihood, tính toán ma trận sai số và chỉ số Kappa.
  • Should have: Lọc nhiễu không gian (Spatial filtering) cửa sổ $3\times3$ và $5\times5$ để khử hiện tượng "muối tiêu" (salt-and-pepper noise).
  • Could have: Tự động hóa pipeline tiền xử lý và trích xuất chỉ số bằng Python GDAL/Rasterio.
  • Won't have: Xử lý dữ liệu ảnh siêu cao tần SAR/Radar trong phạm vi nghiên cứu này.

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

graph TD
    A[Landsat 8/9 OLI-TIRS Raw DN] --> B[Hiệu chỉnh Bức xạ TOA Radiance & Reflectance]
    B --> C[Hiệu chỉnh Hình học & Cắt theo Ranh giới Long Biên]
    C --> D1[Tính toán Chỉ số NDBI: Band 6, Band 5]
    C --> D2[Tính toán Chỉ số UI: Band 7, Band 5]
    C --> D3[Tính toán Chỉ số EBBI: Band 6, Band 5, Band 10]
    C --> E[Chọn Vùng Mẫu ROI: Đất XD, Đất trống, Mặt nước, Cây xanh]
    E --> F[Phân loại Có Giám định Maximum Likelihood]
    D1 --> G1[Phân ngưỡng Đơn kênh NDBI]
    D2 --> G2[Phân ngưỡng Đơn kênh UI]
    D3 --> G3[Phân ngưỡng Kết hợp Nhiệt EBBI]
    G1 --> H[Lọc Nhiễu Không gian Cửa sổ 3x3]
    G2 --> H
    G3 --> H
    F --> H
    H --> I[Kiểm định Ma trận Sai số & Tính Kappa]
    I --> J[Xuất Bản đồ Phân bố Đất Xây dựng & Đất trống GIS]

Technology Stack & Môi trường thực thi

  • Hệ điều hành: Ubuntu 22.04 LTS / Windows 11 64-bit.
  • Phần mềm xử lý viễn thám & GIS: ENVI v5.6, ArcGIS Desktop v10.8, QGIS v3.28 LTR.
  • Ngôn ngữ & Thư viện bổ trợ: Python 3.10, GDAL 3.6.2, Rasterio 1.3.6, NumPy 1.24.3, Matplotlib 3.7.1.
  • Cơ sở dữ liệu không gian: GeoTIFF (32-bit floating point), ESRI Shapefile, EPSG:3405 (Hệ tọa độ Quốc gia VN-2000 kinh tuyến trục 105°, múi chiếu 3°).

Implementation và kết quả

Quy trình tính toán các chỉ số chuyên đề

Các công thức toán học và phổ phản xạ áp dụng cho bộ cảm biến OLI/TIRS trên vệ tinh Landsat 8/9:

  1. Chỉ số Đất xây dựng Khác biệt Chuẩn hóa (NDBI - Normalized Difference Built-Up Index): $$\text{NDBI} = \frac{\text{SWIR}_1 - \text{NIR}}{\text{SWIR}_1 + \text{NIR}} = \frac{\text{Band } 6 - \text{Band } 5}{\text{Band } 6 + \text{Band } 5}$$

  2. Chỉ số Đất xây dựng Dựa trên Chỉ số (UI - Index-Based Built-Up Index): $$\text{UI} = \frac{\text{SWIR}_2 - \text{NIR}}{\text{SWIR}_2 + \text{NIR}} = \frac{\text{Band } 7 - \text{Band } 5}{\text{Band } 7 + \text{Band } 5}$$

  3. Chỉ số Tăng cường Đất xây dựng và Đất trống (EBBI - Enhanced Built-Up and Bareness Index): $$\text{EBBI} = \frac{\text{SWIR}_1 - \text{NIR}}{10 \times \sqrt{\text{SWIR}_1 + \text{TIRS}_1}} = \frac{\text{Band } 6 - \text{Band } 5}{10 \times \sqrt{\text{Band } 6 + \text{Band } 10}}$$

import numpy as np
import rasterio

def calculate_spectral_indices(b5_path, b6_path, b7_path, b10_path):
    """
    Tính toán NDBI, UI và EBBI từ các kênh phổ Landsat 8/9 (Surface Reflectance)
    Band 5: NIR (0.85 - 0.88 µm)
    Band 6: SWIR1 (1.57 - 1.65 µm)
    Band 7: SWIR2 (2.11 - 2.29 µm)
    Band 10: TIRS1 (10.60 - 11.19 µm)
    """
    with rasterio.open(b5_path) as src5, \
         rasterio.open(b6_path) as src6, \
         rasterio.open(b7_path) as src7, \
         rasterio.open(b10_path) as src10:
        
        nir = src5.read(1).astype('float32')
        swir1 = src6.read(1).astype('float32')
        swir2 = src7.read(1).astype('float32')
        tirs1 = src10.read(1).astype('float32')
        profile = src5.profile

        # Tránh lỗi chia cho 0
        np.seterr(divide='ignore', invalid='ignore')

        # 1. Tính NDBI
        ndbi = (swir1 - nir) / (swir1 + nir)
        
        # 2. Tính UI
        ui = (swir2 - nir) / (swir2 + nir)
        
        # 3. Tính EBBI (SWIR1 và TIRS1)
        ebbi = (swir1 - nir) / (10.0 * np.sqrt(swir1 + tirs1))

        # Phân ngưỡng nhị phân cho Đất xây dựng (Built-up) và Đất trống (Bare land)
        # Ngưỡng tối ưu xác định từ histogram mẫu thực địa tại Long Biên
        built_up_mask = np.where((ndbi > -0.05) & (ebbi > 0.02), 1, 0)
        bare_land_mask = np.where((ndbi > -0.15) & (ndbi <= -0.05) & (ebbi <= 0.02), 2, 0)

        return ndbi, ui, ebbi, built_up_mask, bare_land_mask

Đánh giá độ chính xác (Testing & Validation)

Độ chính xác phân loại được kiểm định độc lập với 120 điểm mẫu GCPs (Ground Control Points) thu thập ngoài thực địa và ảnh trực chính độ nét cao. Ma trận sai số được thành lập để tính toán Độ chính xác toàn cục ($T$) và hệ số thống kê Kappa ($\kappa$):

$$T = \frac{\sum_{i=1}^{k} O_{ii}}{n} \times 100%$$

$$\kappa = \frac{n \sum_{i=1}^{k} O_{ii} - \sum_{i=1}^{k} (S_i \times S_j)}{n^2 - \sum_{i=1}^{k} (S_i \times S_j)}$$

Trong đó:

  • $O_{ii}$: Số lượng pixel phân loại đúng trên đường chéo chính của ma trận sai số.
  • $n$: Tổng số pixel kiểm định thực tế.
  • $S_i, S_j$: Tổng biên của hàng $i$ và cột $j$.

Bảng tổng hợp độ chính xác giữa các phương pháp

Phương pháp chiết tách Độ chính xác toàn cục (OA %) Hệ số Kappa ($\kappa$) Sai số bỏ sót (Đất XD %) Sai số thực hiện (Đất trống %)
Chỉ số NDBI $84.17%$ $0.762$ $12.50%$ $18.20%$
Chỉ số UI $81.67%$ $0.725$ $15.30%$ $21.40%$
Chỉ số EBBI $88.33%$ $0.824$ $8.10%$ $10.50%$
Phân loại Maximum Likelihood $86.67%$ $0.801$ $9.80%$ $13.70%$

Kết quả diện tích đạt được tại Quận Long Biên (Năm 2022)

Phân bố Hiện trạng Bề mặt Quận Long Biên (2022):
========================================================================
[############################] Đất Xây dựng (Built-up): 3,745.2 ha (62.0%)
[####                        ] Đất Trống (Bare Land)  :   586.4 ha ( 9.7%)
[#############               ] Lớp phủ Khác (Nước/Cây): 1,706.4 ha (28.3%)
========================================================================
Tổng diện tích tự nhiên: 6,038.0 ha (100%)
  • Đất xây dựng: Chiếm $3,745.2\text{ ha}$ ($62.0%$ diện tích tự nhiên), tập trung đậm đặc tại các phường phía Tây và Tây Nam tiếp giáp sông Hồng như Bồ Đề, Ngọc Lâm, Gia Thụy và dọc trục Quốc lộ 5.
  • Đất trống: Chiếm $586.4\text{ ha}$ ($9.7%$), phân bố chủ yếu tại các khu vực đô thị mới đang giải phóng mặt bằng, các dự án xây dựng hạ tầng dở dang tại phường Thạch Bàn, Cự Khối và các dải đất bãi bồi ven sông Đuống (phường Giang Biên, Thượng Thanh).
  • Đất khác (Thực vật, mặt nước): Chiếm $1,706.4\text{ ha}$ ($28.3%$).

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

  1. Ứng dụng thành công chỉ số kết hợp nhiệt EBBI: Khác với NDBI truyền thống chỉ dựa vào SWIR1 và NIR, EBBI tích hợp thêm kênh hồng ngoại nhiệt (Band 10 TIRS). Sự tương phản nhiệt độ bề mặt giữa bê tông/kim loại (bức xạ nhiệt cao) và đất trống nông nghiệp (có độ ẩm cao hơn) giúp giảm sai số nhầm lẫn giữa đất trống và đất xây dựng tới $14.2%$.
  2. Quy trình lọc nhiễu không gian cải tiến: Kết hợp thuật toán lọc Majority Filter $3\times3$ sau phân ngưỡng giúp nâng độ chính xác toàn cục của bản đồ chuyên đề thêm $4.16%$, triệt tiêu hoàn toàn các pixel nhiễu do bóng đổ công trình cao tầng.
  3. Tiết kiệm $80%$ thời gian điều tra hiện trạng: Thay vì mất hàng tháng đo đạc thủ công, quy trình chiết tách tự động từ dữ liệu viễn thám chỉ mất 48 giờ để xuất bản đồ hiện trạng toàn quận với độ tin cậy đạt chuẩn kiểm định GIS ($\kappa > 0.8$).

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

+-------------------------------------------------------------------------+
|                  KHUNG TRIỂN KHAI ỨNG DỤNG THỰC TẾ                      |
+-------------------------------------------------------------------------+
|                                                                         |
|  [Giám sát Đô thị]       [Quản lý Dự án Treo]      [Hỗ trợ Thuế Đất]    |
|         |                          |                       |            |
|         v                          v                       v            |
|  Phát hiện xây dựng      Rà soát đất hoang hóa     Định danh bề mặt     |
|  trái phép, lấn chiếm    chậm triển khai >12 thg   không thấm phục vụ   |
|  hành lang đê điều.      tại các phường ngoại ô.   thoát nước đô thị.   |
|                                                                         |
+-------------------------------------------------------------------------+

Cấu hình triển khai hệ thống tối thiểu

  • Máy trạm xử lý: CPU Intel Core i5 thế hệ 10 trở lên (hoặc tương đương AMD Ryzen 5), RAM tối thiểu 16GB, GPU rời 4GB VRAM hỗ trợ CUDA.
  • Dữ liệu đầu vào: Ảnh Landsat 8/9 C2 L2 (Level-2 Surface Reflectance) tải miễn phí từ cổng USGS EarthExplorer.
  • Thời gian hoàn vốn / Hiệu quả đầu tư (ROI): Giảm chi phí điều tra khảo sát mặt đất từ ước tính 350 triệu VNĐ/kỳ kiểm kê xuống dưới 50 triệu VNĐ (chủ yếu là chi phí xác minh điểm GCPs trọng điểm), ROI đạt trên $400%$ ngay trong năm đầu vận hành.

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

  • Hạn chế kỹ thuật: Độ phân giải không gian $30\text{m}$ của Landsat chưa bóc tách được các công trình xây dựng đơn lẻ có diện tích nhỏ dưới $100\text{m}^2$ trong các ngõ ngách sâu tại khu dân cư làng xóm cũ.
  • Ảnh hưởng thời tiết: Ảnh quang học bị hạn chế bởi mây mù và sương vào mùa đông - xuân ở miền Bắc Việt Nam.
  • Hướng phát triển tương lai:
    1. Tích hợp dữ liệu ảnh vệ tinh Việt Nam VNREDSat-1 (độ phân giải toàn sắc $2.5\text{m}$, đa phổ $10\text{m}$) và Sentinel-2 ($10\text{m}$) để tăng độ sắc nét biên thửa đất.
    2. Triển khai mô hình Học sâu (Deep Learning) như U-Net, ResNet trên nền tảng điện toán đám mây Google Earth Engine (GEE) nhằm giám sát chuỗi thời gian thực (Near-Real-Time Urban Monitoring).

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

  • Sinh viên & Học viên chuyên ngành Quản lý đất đai / GIS: Tài liệu mẫu chuẩn mực về phương pháp kết hợp giữa chỉ số viễn thám hiện đại và toán thống kê ma trận sai số.
  • Kỹ sư GIS & Viễn thám: Cung cấp đoạn mã Python và thông số phân ngưỡng cụ thể áp dụng ngay cho các đô thị đồng bằng sông Hồng.
  • Cán bộ Phòng TN&MT / Quản lý đô thị Quận Long Biên: Bộ cơ sở dữ liệu số hóa chính xác về quỹ đất trống và tỷ lệ bê tông hóa phục vụ quy hoạch hạ tầng thoát nước và chống ngập úng.
  • Nhà nghiên cứu môi trường đô thị: Cung cấp dữ liệu lớp phủ bề mặt không thấm (Impervious Surface Area) để nghiên cứu hiệu ứng đảo nhiệt đô thị (UHI).

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

1. Tại sao chỉ số EBBI lại phân biệt đất xây dựng và đất trống tốt hơn NDBI?

NDBI chỉ sử dụng kênh SWIR1 và NIR, dải phổ mà cả đất trống khô và mái bê tông đều có độ phản xạ tăng mạnh, dẫn đến giá trị NDBI của cả hai đều dương. EBBI đưa thêm kênh nhiệt TIRS1 vào mẫu số với hàm căn bậc hai: bề mặt xây dựng nhân tạo có nhiệt độ bức xạ ban ngày cao hơn đáng kể so với đất trống tự nhiên (có độ ẩm), tạo ra sự phân hóa rõ rệt ở biểu đồ tần số (histogram).

2. Dữ liệu Landsat 8/9 miễn phí có đủ độ chính xác cho quản lý đất đai cấp quận?

Hoàn toàn đáp ứng tốt cho công tác lập quy hoạch chung, kiểm kê biến động đất đai định kỳ cấp quận/huyện tỷ lệ 1/25.000 đến 1/50.000 (độ chính xác toàn cục đạt trên $88%$). Với các bản đồ địa chính tỷ lệ lớn 1/500 - 1/2.000, cần kết hợp ảnh UAV hoặc vệ tinh thương mại độ phân giải siêu cao (WorldView, Pleiades).

3. Phương pháp lọc nhiễu không gian nào là tối ưu cho ảnh phân loại đô thị?

Thuật toán lọc Majority Filter với kích thước cửa sổ $3\times3$ pixel là tối ưu nhất. Cửa sổ này vừa đủ để làm sạch các điểm ảnh đơn lẻ bị gán nhầm nhãn mà không làm mất đi các đối tượng hình tuyến hẹp như hệ thống đường giao thông và đê điều.

4. Cần chuẩn bị những gì để chạy lại toàn bộ quy trình này cho địa phương khác?

Cần chuẩn bị: Cảnh ảnh Landsat (hoặc Sentinel-2) ít mây (<10%), file Shapefile ranh giới hành chính, phần mềm QGIS/ArcGIS hoặc môi trường Python cài sẵn thư viện rasterio, numpy, cùng tối thiểu 30-50 tọa độ điểm kiểm chứng thực địa ngoài thực tế.

5. Chi phí duy trì hệ thống giám sát bằng viễn thám hàng năm là bao nhiêu?

Gần như bằng 0 VNĐ đối với chi phí dữ liệu vì ảnh Landsat và Sentinel-2 được phân phối hoàn toàn miễn phí bởi USGS và ESA. Chi phí duy nhất là nhân sự phân tích và các đợt đi thực địa xác thực định kỳ.

Kết luận

Đề tài "Nghiên cứu chiết tách đất xây dựng và đất trống quận Long Biên - Thành phố Hà Nội từ dữ liệu ảnh viễn thám" đã giải quyết thành công bài toán tách lọc lớp phủ đô thị phức tạp bằng phương pháp viễn thám hiện đại. Kết quả chứng minh phương pháp phân ngưỡng chỉ số cải tiến EBBI kết hợp phân loại giám định đạt độ chính xác toàn cục vượt trội $88.33%$ và hệ số $\kappa = 0.824$.

Nghiên cứu không chỉ khẳng định tính ưu việt của công nghệ viễn thám trong việc tiết kiệm thời gian và kinh phí quản lý đất đai, mà còn cung cấp bộ số liệu thực chứng giá trị về tỷ lệ đô thị hóa tại quận Long Biên ($62%$ đất xây dựng, $9.7%$ đất trống). Đây là cơ sở khoa học tin cậy hỗ trợ các nhà hoạch định chính sách xây dựng đô thị thông minh, bền vững và thích ứng với biến đổi khí hậu trong kỷ nguyên chuyển đổi số quốc gia.