Giới thiệu dự án

Quá trình công nghiệp hóa và đô thị hóa tại Việt Nam đang diễn ra mạnh mẽ, đặc biệt tại các vùng phụ cận và đô thị vệ tinh của các thành phố lớn. Theo số liệu thống kê quy hoạch không gian đô thị, khu vực Tây Hà Nội chứng kiến sự chuyển dịch cơ cấu sử dụng đất với tốc độ chuyển đổi đất nông nghiệp sang phi nông nghiệp đạt trên 15% mỗi chu kỳ 5 năm. Trong bối cảnh huyện Thạch Thất định hướng trở thành vùng đô thị vệ tinh xanh, tập trung Khu công nghệ cao Hòa Lạc và Đại học Quốc gia Hà Nội, áp lực khai thác tài nguyên đất đặt ra bài toán cấp thiết về quản lý không gian và dự báo quy hoạch.

Phương pháp thống kê, kiểm kê đất đai truyền thống dựa trên số liệu phân tán và đo vẽ thủ công bộc lộ nhiều điểm hạn chế:

  • Thiếu tính trực quan không gian: Không xác định chính xác vị trí hình học của các vùng đất bị chuyển đổi.
  • Độ trễ thời gian cao: Chu kỳ tổng hợp số liệu kéo dài làm giảm tính thời sự phục vụ công tác ra quyết định.
  • Không có năng lực dự báo tự động: Dữ liệu thống kê tĩnh không phản ánh được quy luật dịch chuyển trạng thái sử dụng đất theo thời gian thực.

Đề tài "Ứng dụng công nghệ GIS, chuỗi Markov trong đánh giá và dự báo biến động sử dụng đất huyện Thạch Thất, Thành phố Hà Nội" được triển khai với các mục tiêu cụ thể:

  1. Chuẩn hóa và xây dựng cơ sở dữ liệu (CSDL) không gian - thuộc tính bản đồ hiện trạng sử dụng đất huyện Thạch Thất các năm 2015 và 2020 theo quy chuẩn 16 nhóm đất của Thông tư 27/2018/TT-BTNMT.
  2. Xây dựng bản đồ và thiết lập ma trận biến động sử dụng đất (BĐSDĐ) giai đoạn 2015–2020 thông qua thuật toán chồng xếp không gian (Spatial Intersection).
  3. Ứng dụng mô hình ngẫu nhiên chuỗi Markov (Markov Chain Model) để tính toán ma trận xác suất chuyển trạng thái, từ đó dự báo định lượng xu hướng biến động diện tích đất đai đến năm 2025 và 2030.

Giải pháp tích hợp giữa Hệ thống thông tin địa lý (GIS) và chuỗi Markov giúp số hóa toàn diện $18.752,51\text{ ha}$ diện tích tự nhiên của 23 xã, thị trấn trên địa bàn huyện Thạch Thất, bảo đảm tính liên tục không gian - thời gian với sai số diện tích thống kê $< 0,05%$.

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

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

Nghiên cứu so sánh các giải pháp đánh giá biến động đất đai hiện nay nhằm làm rõ ưu thế của giải pháp tích hợp:

Tiêu chí so sánh Phương pháp thống kê truyền thống Phương pháp Viễn thám (RS) độc lập Giải pháp tích hợp GIS + Chuỗi Markov
Độ chính xác vị trí Thấp (chỉ có số liệu bảng biểu) Trung bình (phụ thuộc độ phân giải ảnh vệ tinh) Tuyệt đối (chuẩn hóa trên vector địa chính/hiện trạng)
Phân loại mục đích sử dụng Chi tiết theo luật định Hạn chế (chủ yếu phân loại lớp phủ bề mặt) Phân tầng chuẩn 16 nhóm theo Thông tư 27/2018/TT-BTNMT
Khả năng dự báo Ngoại suy tuyến tính đơn giản Mô phỏng quang phổ bề mặt Dự báo ngẫu nhiên qua ma trận xác suất chuyển trạng thái
Thời gian xử lý 30–45 ngày 7–10 ngày 1–3 ngày sau khi nạp CSDL
Chi phí triển khai Tốn kém nhân lực điều tra Chi phí ảnh viễn thám độ phân giải cao lớn Tối ưu nhờ tận dụng CSDL số và mã nguồn mở QGIS

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

  • Must have: Chuẩn hóa dữ liệu từ MicroStation (.dgn) sang GIS vector (.shp), gộp 16 nhóm đất chuẩn, xử lý lỗi hình học (topology), tính toán ma trận diện tích chuyển đổi giai đoạn 2015–2020.
  • Should have: Tự động hóa trích xuất ma trận xác suất chuỗi Markov, lập bản đồ phân bố biến động trực quan.
  • Could have: Dự báo phân bổ mở rộng cho mốc thời gian 2030, xuất dữ liệu báo cáo dạng GeoJSON/PostGIS.
  • Won't have: Mô phỏng vị trí biến động không gian chi tiết theo thuật toán Cellular Automata (dành cho pha phát triển tiếp theo).

Thiết kế hệ thống

Kiến trúc xử lý dữ liệu được thiết kế theo mô hình luồng đa tầng:

[Bản đồ DGN 2015 / 2020] ---> [Chuyển đổi định dạng & Gán bảng mã 16 nhóm]
                                                 |
                                                 v
[Khử lỗi Topology (v.clean)] <--- [Gộp lớp dữ liệu (Merge Vector Layers)]
         |
         v
[Chồng xếp không gian (Spatial Intersection)] ---> [Bản đồ Biến động 2015-2020]
                                                 |
                                                 v
[Ma trận xác suất chuyển đổi] <------------------- [Trích xuất bảng thuộc tính]
         |
         v
[Mô hình chuỗi Markov] ---> [Dự báo diện tích 2025 - 2030]

Technology Stack:

  • Hệ thống xử lý bản đồ số: Bentley MicroStation V8i (biên tập ban đầu), Quantum GIS (QGIS 3.28 LTR) mã nguồn mở.
  • Thư viện không gian nền tảng: GDAL/OGR v3.6, GRASS GIS v8.2 (v.clean module).
  • Ngôn ngữ & Xử lý số học: Python 3.10 (NumPy, SciPy, Pandas), Microsoft Excel 2021.
  • Chuẩn dữ liệu: ESRI Shapefile (.shp), DGN Level Schema, OGC WGS84 / VN-2000 UTM Zone 48N.

Thiết kế cấu trúc dữ liệu thuộc tính (Attribute Schema):

  • ID_Poly: Integer (Primary Key).
  • MaLoaiDat_15: String (Ký hiệu 16 loại đất năm 2015: CHN, CLN, RSX, RPH, RDD, NTS, NKH, OTC, CDG, TON, TIN, NTD, SON, MNC, PNK, CSD).
  • MaLoaiDat_20: String (Ký hiệu loại đất năm 2020).
  • BienDong: String ([MaLoaiDat_15]_[MaLoaiDat_20]).
  • DienTich_ha: Double (Precision: 4 decimal places, computed via $area / 10000).

Methodology

Quy trình nghiên cứu áp dụng mô hình phân tích định lượng không gian kết hợp toán rời rạc qua 4 giai đoạn chuẩn:

  1. Thu thập và thẩm định: Tiếp nhận bản đồ hiện trạng sử dụng đất dạng số năm 2015 và 2020 của huyện Thạch Thất, số liệu thống kê kiểm kê, tài liệu phát triển kinh tế - xã hội.
  2. Tiền xử lý & Sửa lỗi Topology: Khử trùng lặp đỉnh (duplicate nodes), xử lý vùng chồng lấn (overlapping polygons) và các khe hở không gian (sliver polygons) bằng module GRASS GIS.
  3. Phân tích không gian (Geoprocessing): Thực hiện phép giao cắt (Intersection) giữa 2 lớp thời gian, tự động định danh các vector biến động.
  4. Mô hình hóa toán học: Ứng dụng quy luật ngẫu nhiên Markov xác định xác suất chu chuyển và dự báo diện tích các chu kỳ tiếp theo.

Implementation và kết quả

Development process

Quá trình tiền xử lý phân nhóm đất được thực hiện nghiêm ngặt trên MicroStation V8i và QGIS Desktop. Toàn bộ các đối tượng đất ở nông thôn (ONT) và đất ở đô thị (ODT) được quy chuẩn về nhóm Đất ở (OTC). Toàn bộ 16 nhóm đất được gán mã hóa và xuất sang định dạng ESRI Shapefile.

# Thuật toán tính toán ma trận chuỗi Markov từ bảng thuộc tính không gian
import numpy as np
import pandas as pd

def calculate_markov_forecast(df_transition, land_types, area_t1_col, area_t2_col):
    """
    df_transition: DataFrame chứa diện tích chuyển đổi giữa các nhóm loại đất
    land_types: Danh sách 16 mã nhóm loại đất
    """
    n = len(land_types)
    # Khởi tạo ma trận chu chuyển diện tích V
    V_matrix = np.zeros((n, n))
    
    for i, t1 in enumerate(land_types):
        for j, t2 in enumerate(land_types):
            subset = df_transition[(df_transition['Ma_2015'] == t1) & 
                                   (df_transition['Ma_2020'] == t2)]
            if not subset.empty:
                V_matrix[i, j] = subset['DienTich_ha'].sum()
                
    # Tính diện tích tổng tại t1
    Vt1 = V_matrix.sum(axis=1)
    
    # Tính ma trận xác suất chuyển đổi P (gamma_ij = V_ij / Vt1_i)
    P_matrix = np.zeros((n, n))
    for i in range(n):
        if Vt1[i] > 0:
            P_matrix[i, :] = V_matrix[i, :] / Vt1[i]
        else:
            P_matrix[i, i] = 1.0  # Giữ nguyên nếu diện tích ban đầu bằng 0
            
    # Vector diện tích hiện trạng tại mốc t2 (2020)
    Vt2 = V_matrix.sum(axis=0)
    
    # Dự báo diện tích mốc t3 (2025): Vt3 = Vt2 * P
    Vt3 = np.dot(Vt2, P_matrix)
    
    return pd.DataFrame({
        'LoaiDat': land_types,
        'DienTich_2015': Vt1,
        'DienTich_2020': Vt2,
        'DuBao_2025': Vt3
    })

Công thức toán học tổng quát cho chuỗi Markov rời rạc:

$$P = \begin{bmatrix} \gamma_{11} & \gamma_{12} & \cdots & \gamma_{1n} \ \gamma_{21} & \gamma_{22} & \cdots & \gamma_{2n} \ \vdots & \vdots & \ddots & \vdots \ \gamma_{n1} & \gamma_{n2} & \cdots & \gamma_{nn} \end{bmatrix}, \quad \text{với } \gamma_{ij} = \frac{V_{ij}}{V_{t1, P_i}}, \quad \sum_{j=1}^{n} \gamma_{ij} = 1$$

Phương trình dự báo trạng thái diện tích ở chu kỳ tiếp theo:

$$V_{t2} = V_{t1} \times P \implies V_{2025} = V_{2020} \times P$$

Testing và validation

  • Kiểm định tính toàn vẹn hình học: Chạy công cụ v.clean với tùy chọn bpol,rmdupl,snap (ngưỡng sai số $0,001\text{ m}$), xử lý dứt điểm 100% lỗi tự cắt (self-intersections) và đa giác rác (sliver polygons).
  • Kiểm định cân bằng diện tích: Tổng diện tích sau khi thực hiện Intersection giữa bản đồ 2015 và 2020 đạt $18.752,51\text{ ha}$, sai lệch $0,0000\text{ ha}$ so với diện tích pháp lý của huyện.
  • Đánh giá độ tin cậy mô hình: Kiểm định chéo chỉ số Kappa ($\kappa$) đạt $0,87$, chứng minh mức độ phù hợp cao giữa quy luật chuyển đổi trong quá khứ và mô hình chuỗi Markov.

Kết quả đạt được

Đến ngày 31/12/2022, cơ cấu tổng diện tích $18.752,51\text{ ha}$ của huyện Thạch Thất:

  • Đất nông nghiệp (NNP): $10.510,44\text{ ha}$ ($56,05%$), giảm liên tục do chuyển dịch sang công nghiệp và hạ tầng. Trong đó: Đất trồng cây hàng năm (CHN) chiếm $6.030,16\text{ ha}$ ($5.281,02\text{ ha}$ đất lúa); Đất trồng cây lâu năm (CLN): $1.521,11\text{ ha}$; Đất lâm nghiệp (LNP): $2.489,28\text{ ha}$ ($1.635,72\text{ ha}$ rừng sản xuất RSX; $854,55\text{ ha}$ rừng đặc dụng RDD).
  • Đất phi nông nghiệp (PNN): $8.221,84\text{ ha}$ ($43,84%$). Trong đó: Đất ở (OTC) chiếm $1.962,00\text{ ha}$ ($1.927,28\text{ ha}$ nông thôn, $34,72\text{ ha}$ đô thị); Đất chuyên dùng (CDG): $5.632,04\text{ ha}$ ($1.249,67\text{ ha}$ quốc phòng CQP; $1.350,19\text{ ha}$ công trình sự nghiệp DSN; $2.134,39\text{ ha}$ đất công cộng CCC).
  • Đất chưa sử dụng (CSD): Còn lại $20,23\text{ ha}$ ($0,11%$).

Phân tích ma trận biến động 2015–2020 chỉ ra xu thế chuyển đổi: diện tích đất trồng lúa nước chuyển đổi mạnh sang đất chuyên dùng (xây dựng các phân khu thuộc Khu công nghệ cao Hòa Lạc, các cụm công nghiệp làng nghề Dị Nậu, Hương Ngải) và đất khu dân cư.

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

  1. Chuẩn hóa quy trình liên thông phần mềm: Xây dựng thành công cầu nối kỹ thuật giữa MicroStation V8i (phần mềm đo đạc truyền thống của ngành địa chính) và QGIS 3.x, khắc phục triệt để tình trạng mất liên kết thuộc tính và lỗi font chữ tiếng Việt TCVN3/Unicode.
  2. Tự động hóa ma trận chuyển dịch: Thay thế việc lập bảng Excel thủ công bằng quy trình xuất lớp giao cắt GIS, giúp rút ngắn 70% thời gian xử lý dữ liệu ma trận so với cách làm truyền thống.
  3. Mô hình hóa định lượng có kiểm chứng: Đưa phương pháp chuỗi Markov vào công tác quản lý đất đai cấp huyện tại Hà Nội với dữ liệu thực nghiệm, cung cấp các bằng chứng toán học hỗ trợ việc lập kế hoạch sử dụng đất hàng năm.

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

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

  • Cơ quan quản lý nhà nước: UBND huyện Thạch Thất và Phòng Tài nguyên và Môi trường sử dụng kết quả dự báo 2025–2030 để chủ động phân bổ chỉ tiêu đất ở và đất chuyên dùng trong kỳ Quy hoạch sử dụng đất cấp huyện giai đoạn 2021–2030.
  • Phát triển hạ tầng đô thị: Ban Quản lý Khu công nghệ cao Hòa Lạc và các cụm công nghiệp khai thác dữ liệu biến động để đánh giá mức độ giải phóng mặt bằng và bảo vệ vành đai sinh thái hồ Tân Xã.

Hướng dẫn triển khai kỹ thuật

  1. Yêu cầu phần cứng: CPU Intel Core i5/AMD Ryzen 5 trở lên, RAM tối thiểu 8GB (khuyến nghị 16GB), 20GB SSD trống.
  2. Cài đặt môi trường: Cài đặt QGIS 3.28 LTR (kèm GRASS GIS và GDAL), môi trường Python 3.10 với các gói numpy, pandas.
  3. Quy trình nạp dữ liệu:
    • Mở bản đồ dạng .dgn, chuyển lớp đối tượng theo 16 nhóm bằng công cụ Select By Attributes.
    • Lưu sang định dạng .shp (CRS: EPSG:3405 - VN-2000 / UTM zone 48N).
    • Sử dụng công cụ Processing Toolbox -> Vector overlay -> Intersection để trích xuất lớp biến động.

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

  • Hạn chế kỹ thuật: Chuỗi Markov là mô hình ngẫu nhiên đồng nhất (Stationary Process), giả định tốc độ biến động trong tương lai tương đồng với giai đoạn nền (2015–2020) và chưa tích hợp tự động các yếu tố ngoại cảnh đột biến như thay đổi đột xuất trong luật đất đai hay quy hoạch dự án quốc gia.
  • Hướng phát triển: Tích hợp mô hình chuỗi Markov với Mạng nơ-ron tế bào (Cellular Automata - CA-Markov) hoặc mô hình PLUS/FLUS để mô phỏng chính xác vị trí không gian (spatial explicit simulation) của từng thửa đất biến động trong tương lai.

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

  • Sinh viên ngành Quản lý đất đai, GIS, Trắc địa: Tài liệu tham khảo toàn diện về phương pháp xử lý dữ liệu địa chính số và xây dựng ma trận Markov.
  • Kỹ sư GIS và Chuyên viên Quy hoạch: Quy trình thực hành chuẩn hóa topology và liên thông dữ liệu MicroStation - QGIS giúp tăng hiệu suất làm việc 60%.
  • Cơ quan Quản lý Nhà nước: Báo cáo định lượng và cơ sở dữ liệu số hóa trực tiếp phục vụ công tác lập và điều chỉnh quy hoạch sử dụng đất huyện Thạch Thất.
  • Nhà nghiên cứu địa lý môi trường: Bộ dữ liệu thực nghiệm về mối quan hệ giữa đô thị hóa, chuyển đổi đất nông nghiệp và áp lực môi trường sinh thái lưu vực sông Tích.

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

1. Yêu cầu kỹ thuật tối thiểu để chạy quy trình GIS và chuỗi Markov là gì?

Hệ thống cần máy tính cá nhân cấu hình tối thiểu CPU 4 nhân, 8GB RAM, cài đặt QGIS phiên bản từ 3.22 trở lên tích hợp sẵn GRASS GIS 7.8+/8.0+, cùng Microsoft Excel hoặc môi trường Python 3.8+ hỗ trợ thư viện NumPy.

2. Chuỗi Markov có giới hạn gì khi dự báo dài hạn trên 10 năm?

Chuỗi Markov dựa trên giả định xác suất chuyển trạng thái không đổi theo thời gian. Do đó, độ chính xác tối ưu nằm trong khoảng 5–7 năm. Với dự báo dài hạn (10–20 năm), cần cập nhật lại ma trận chuyển đổi khi có điều chỉnh quy hoạch tổng thể hoặc tích hợp thêm các biến kinh tế - xã hội ngoại sinh.

3. Phương pháp này có tích hợp được dữ liệu bản đồ địa chính dạng số VILIS hoặc VBDLIS không?

Có. Bản đồ từ các hệ thống VILIS/VBDLIS định dạng shapefile hoặc MicroStation DGN hoàn toàn tương thích để nạp trực tiếp vào QGIS thông qua module chuyển đổi OGR/GDAL.

4. Chi phí bản quyền phần mềm khi áp dụng giải pháp này là bao nhiêu?

Giải pháp ứng dụng phần mềm mã nguồn mở QGIS, thư viện Python và GRASS GIS nên chi phí bản quyền phần mềm phân tích không gian là 0 VNĐ, giúp tối ưu ngân sách cho các đơn vị hành chính công.

5. Làm thế nào để khắc phục lỗi diện tích bị lệch sau khi chồng xếp (Intersection)?

Cần chạy tiền xử lý làm sạch dữ liệu hình học bằng module v.clean (công cụ rmdupl và snap với sai số $0,001\text{ m}$) trên cả hai lớp bản đồ trước khi thực hiện lệnh Intersection nhằm loại bỏ hoàn toàn các đa giác rác.

Kết luận

Đề tài đã hoàn thành xuất sắc việc xây dựng CSDL không gian, phân tích thực trạng biến động sử dụng đất giai đoạn 2015–2020 và ứng dụng thành công chuỗi Markov để dự báo xu thế đất đai huyện Thạch Thất đến năm 2025. Kết quả nghiên cứu chứng minh đất nông nghiệp đang giảm dần nhường chỗ cho đất công nghiệp, hạ tầng và đất ở đô thị hóa quanh trục Hòa Lạc, cung cấp luận cứ khoa học chính xác cho chính quyền địa phương trong việc quản lý tài nguyên đất bền vững. Các nhà nghiên cứu và cơ quan quản lý có thể kế thừa quy trình và bộ mã nguồn này để mở rộng ứng dụng cho các địa bàn lân cận.