陳建國 ，王晨輝 ，徐緒堪 ，嵇慶才

（1. 寧夏水利電力工程學校，寧夏 銀川 750006；2. 河海大學商學院，江蘇 常州 213022；3. 常州工業(yè)大數(shù)據(jù)挖掘與知識管理重點實驗室，江蘇 常州 213022；4. 鎮(zhèn)江新區(qū)城鄉(xiāng)建設(shè)局，江蘇 鎮(zhèn)江 212132）

0 引言

隨著城鄉(xiāng)建設(shè)的快速發(fā)展，城市現(xiàn)代化進程不斷加快，大量河湖違章建筑影響城鄉(xiāng)水域公共空間，河湖流域違法的圈圩和建設(shè)嚴重影響城市防洪排澇，給城市社會經(jīng)濟穩(wěn)定發(fā)展帶來較大潛在危害。城市河湖違建整治是城市全面推進河（湖）長制、建設(shè)人水和諧生態(tài)的重要舉措之一，根據(jù)江蘇省委辦公廳要求，全省各地對轄區(qū)內(nèi)河湖流域逐一排查，確保各類違法圈圩和建設(shè)不遺漏，2020 年底將全面整治到位。

在過去的河湖“清四亂”專項活動中，水域違章建筑檢測主要依賴人工對 0.3 m 分辨率遙感影像進行比對排查，為解決人力資源有限、人工比對耗時長的問題，提出一種基于遙感影像變化技術(shù)的水域違章建筑識別模型，并通用實例驗證識別的效果，為水域違章建筑檢測提供科學快捷的方法。

1 相關(guān)研究

眾多學者在違章建筑檢測方面提出了很多方案，詹金瑞[1]提出了以歷年衛(wèi)星影像為基礎(chǔ)，結(jié)合實地測繪和規(guī)劃信息的違章建筑檢測總體架構(gòu)；朱建偉等[2]提出利用無人機低空環(huán)拍影像構(gòu)建三維模型進而識別違章建筑的系統(tǒng)；黃蓉等[3]以青島市為例，設(shè)計了基于高分辨率遙感影像變化檢測技術(shù)的違章建筑檢測技術(shù)路線，解決了實地測繪的低效率和無人機環(huán)拍的高技術(shù)門檻問題，成為了目前主流的違建檢測方法。

遙感影像變化檢測技術(shù)是指通過圖像處理等技術(shù)對不同時期同一區(qū)域的遙感影像數(shù)據(jù)進行分析以確定 2 個時期影像間的變化[4]。依據(jù)處理單元粒度大小不同，可以將現(xiàn)有的變化檢測算法分為兩大類：1）基于像素點的變化檢測。以像素為最小單位進行處理，試圖將所有像素歸為發(fā)生和未發(fā)生變化2 類的方法，例如：Malila 最早提出變化矢量分析法（CVA）[5]，利用多波段的遙感影像數(shù)據(jù)將每個像素描述為 1 個一維列向量，進而計算前后時期同位置像素之間的向量差值，代表變化的強度；陳晉等[6]將 CVA 算法用于國內(nèi)土地利用/覆蓋變化檢測，并提出一種雙窗口變步長閾值搜尋的方法，在變化檢測的自動程度上更進一步；黃維等[7]結(jié)合主成分分析（PCA）和 CVA 算法，對多波段影像數(shù)據(jù)提取第一主成分后進行差值運算和閾值劃分，降低了圖像噪聲的影響；趙秋菊[8]進一步改進 PCA 算法，以差值圖像的一倍標準差作為變化閾值，自動確定出變化區(qū)域。這類算法一般對輸入影像向量表示并求差值后，人工設(shè)定或自動搜索閾值，確定每個像素是否發(fā)生變化。2）面向?qū)ο蟮淖兓O(jiān)測。常常采用先分類后檢測變化的方法，首先對某一類型像素集合看作對象，研究改對象范圍內(nèi)是否發(fā)生變化。例如：李亮等[9]通過影響分割獲取像斑，結(jié)合 EM 算法和貝葉斯最小錯誤理論確定每個像斑的變化/非變化類別歸屬；張亞亞等[10]提出一種基于對象的語義網(wǎng)遙感圖像知識分類框架通過制定規(guī)則解決對象分割問題。近年來，隨著計算機運算性能的提高，基于深度學習的監(jiān)督分類方法往往在實際中取得更好的效果，張曉東等[11]將目標檢測領(lǐng)域的主流網(wǎng)絡(luò)Faster R-CNN 應用到高分辨率遙感圖像變化檢測中取得了理想的效果。

綜上所述，國內(nèi)外學者為違章建筑檢測貢獻了有效的系統(tǒng)架構(gòu)和技術(shù)方法，然而傳統(tǒng)的基于像素的變化檢測算法容易出現(xiàn)“椒鹽”現(xiàn)象，難以滿足違章建筑粗篩選的精度需求，而基于深度學習的運算時間高度依賴計算機的性能，難以適用實時的檢測系統(tǒng)。因此，本研究提出基于改進主成分分析和k-Means 聚類的遙感影像快速變化檢測算法以實現(xiàn)水域違章建筑粗篩選。

2 遙感影像快速變化檢測算法設(shè)計

2.1 問題定義

給定當前時期某區(qū)域遙感影像，欲篩選出違建區(qū)域，首先獲取過去某時期同區(qū)域遙感影像；接著將 2 個時期影像數(shù)據(jù)代入算法，獲得變化區(qū)域圖斑，并在當前時期遙感影像中標記出變化位置；最后將標注有變化區(qū)域的遙感影像推送至核查人員處，審查變化區(qū)域是否發(fā)生違章建筑。其中，獲取變化區(qū)域是核心算法，X1，X2是同一地區(qū)不同時期的 2 張遙感影像圖片，滿足以下條件：

式中：H和W分別表示圖像的高度和寬度；x1(i，j)和x2(i，j) 表示像素點 (i，j) 的灰度值。

目標是生成一張變化區(qū)域圖片Xc= {xc(i，j) | 1 ≤i≤H，1 ≤j≤W}，滿足：

2.2 算法流程

算法流程主要包括計算差分圖像、構(gòu)造圖像變化特征矩陣和像素點二分類 3 步。1）差分圖像通過2 個時期遙感影像灰度值差值計算；2）圖像變化特征通過 2 次劃分像素塊，將差分圖像的每個像素點變化信息表示成向量，并利用 PCA 算法降維；3）像素點二分類是利用改進的k-Means 聚類將圖像像素點分為發(fā)生和未生變化 2 類。算法流程如圖 1 所示。

圖 1 算法流程圖

2.2.1 計算差分圖像

在變化檢測模型中有 2 種主流方法：合并 2 個不同時間的圖像；用后一時期影像減去前一時期影像得到差分圖像。

本研究提出的快速檢測的算法是基于差分圖像的，差分圖像Xd表示圖像X1和X2每個像素點灰度值的絕對差值。

2.2.2 構(gòu)造圖像變化特征矩陣

將差分圖像Xd劃分成大小為h×h(h≥2) 的像素塊，則

式中：xd(m，n) 表示位于m行n列的像素塊，xd(m，n) 簡化寫作將h×h大小的像素塊依據(jù)行列順序展開為h2×1 的向量表示第k個像素塊。

式中：K為劃分后像素塊的總個數(shù)每個向量距離平均向量Ψ的距離

PCA 算法試圖尋找一組N(N≤h2) 個線性無關(guān)的綜合指標代替原來的h2個指標，對集合Δk使用PCA 算法，首先計算協(xié)方差矩陣其中i和j的協(xié)方差cij= Cov(Δi，Δj)，i，j=1，2，…，h2。然后計算協(xié)方差矩陣C的特征向量與特征值，由于矩陣C= (cij)h2×h2，所以共有h2個特征向量es與特征值λs，按照特征值遞減順序排列，選取前N(N≤h2) 個特征向量組成特征向量空間EVS= [e1，e2，…，eN]T，N≤h2。

最后重新將差分圖像Xd劃分成大小為h×h(h≥2) 的像素塊，這次劃分采用重疊的方式，即同一行或同一列相鄰的 2 個像素塊之間的間隔為定長，用stride表示，設(shè)置每個stride為 1 個像素，劃分后Xd={xd(m，n) |m+h=W，n+h=H}，每個像素塊用來表示其中心位置像素點的變化特征，同樣按行列展開為h2×1 的向量，記作(H-h)，按式 (6) 和 (7) 將映射到EVS中去：

最終通過 PCA 算法，將差分圖像大部分區(qū)域像素的變化特征轉(zhuǎn)化為維度間線性無關(guān)的N維向量，稱為圖像變化特征矩陣。

2.2.3 像素點二分類

使用聚類算法依據(jù)上一步得到的圖像變化特征矩陣對差分圖像像素點進行二分類，相較于傳統(tǒng)的k-Means 聚類（k= 2）[12]，做出如下改進：

1）比較使用多種k值，由于不同時期的衛(wèi)星圖片在整體上區(qū)別不大，而改變總在細微處，所以可以認為包含像素點最少的一類，是發(fā)生改變的一類。

2）為選擇初始聚類中心，優(yōu)化k-Means 聚類方法，獲得更快的迭代速度和更好的聚類效果，采用如下算法：a. 輸入數(shù)據(jù)集X，聚類中心個數(shù)k，在數(shù)據(jù)集中隨機選取 1 個樣本作為初始聚類中心C1。b. 計算每個樣本到與之最近的 1 個聚類中心的距離D(x)，x∈X；計算每個樣本被選為下一聚類中心的概率選擇概率最大的樣本作為下一聚類中心。c. 重復步驟 b 直到選擇出k個聚類中心。

設(shè)定“1”表示某像素點發(fā)生改變，“0”表示該像素點沒有發(fā)生改變，因此分類結(jié)果圖中，黑色區(qū)域表示沒有變化，白色區(qū)域表示發(fā)生變化。

3 實驗與分析

3.1 數(shù)據(jù)來源與預處理

用一組拍攝時間分別為 2015 年 6 月和 2018 年6 月的江蘇省常州市某區(qū)域遙感影像數(shù)據(jù)進行實驗，該區(qū)域遙感影像圖片如圖 2 所示。

圖 2 常州市某區(qū)域遙感影像圖片

圖片大小為 8 000 像素×7 000 像素，圖片格式為 TIFF；影像信息格式為 TFW。圖 2 所包含影像信息內(nèi)容及含義如表 1 所示，其中：像素分辨率表示單位像素寬度代表的實際距離，單位為 m/像素，在 0.3 m/像素的分辨率尺度下，每張圖片的實際地理長寬分別是H= 0.3 ×8 000 = 2 400 m，W= 0.3 ×7 000 = 2 100 m。

數(shù)據(jù)預處理主要包括：

表 1 某遙感影像圖片信息

1）成對圖片位置校正。由于獲取的遙感影像可能存在圖片大小誤差、地理位置偏差等異常情況，誤差示意圖如圖 3 所示。

圖 3 誤差示意圖

若A，B點的像素坐標分別為 (x1，x2)，(y1，y2)，地理坐標分別為 (x1'，x2')，(y1'，y2')，2 個圖像素分辨率都是α，那么A，B像素坐標及地理信息滿足下式：

分別對 2 個時期遙感影像圖片剪裁出圖 3 中黃色區(qū)域，代替原始數(shù)據(jù)，完成圖片位置校正。

2）切割圖像。由于原始圖像過大，容易造成內(nèi)存溢出，所以將原始圖像切割成長寬都為 1 000 像素的矩形。

3）轉(zhuǎn)化成灰度圖。將紅、綠、藍三通道圖像轉(zhuǎn)換為單通道灰度圖。

預處理完成后，每次輸入算法的都是 2 張1 000 像素×1 000 像素的單通道同區(qū)域不同時期遙感影像圖片。

3.2 對比實驗

在本研究提出的方法中，主要超參數(shù)包括像素塊大小h和聚類中心個數(shù)k，選用如表 2 所示的 6 組超參數(shù)進行實驗。

在實驗中，前 3 個主成分的累積貢獻率達到了99.7%，故取主成分個數(shù)N= 3，對比實驗結(jié)果如圖4 所示。

對比方案 1，2，3 或 4，5，6，發(fā)現(xiàn)像素塊大小h對實驗結(jié)果的影響不顯著；對比方案 1，2 或3，4 或 5，6，發(fā)現(xiàn)聚類中心個數(shù)k決定了變化區(qū)域的閾值，k值越小，結(jié)果圖像的椒鹽化現(xiàn)象就越顯著，誤分類率也隨之提高。圖 5 所示是圖像上部一塊區(qū)域的實際圖片與實驗結(jié)果，從實際圖像中可以看出，該區(qū)域變化主要是下方農(nóng)田改建為道路的施工用地，在k= 3 個時，算法將上方一塊池塘誤歸類為變化區(qū)域，這是由于前后 2 個時期遙感影像光線不同導致的圖像顏色深淺不一致；而k= 5 個的方案很好地解決了這一誤分類。

表 2 實驗超參數(shù)選擇

圖 4 各參數(shù)方案提取的變化圖像

圖 5 某區(qū)域原始圖片與實驗結(jié)果

在定性分析的基礎(chǔ)上，保持像素塊大小h= 5 像素恒定，選用 2～6 這 5 種不同k值，使用樣點檢驗的方式進行定量分析。依據(jù)特定原則選取 500 個樣本點比較各算法的精度[4]。人工判斷這 300 個樣本點的集合歸屬，其中變化樣本 166 個，未變化樣本334 個，計算各個方法的召回率R、準確率P和總體精度結(jié)合傳統(tǒng) CVA 方法比較結(jié)果如表 3所示。

表 3 算法精度比較

從表 3 可以看出：隨著k值增加，準確率提高至收斂，召回率卻不斷降低；同時，本研究算法普遍優(yōu)于傳統(tǒng)的 CVA 算法，在k= 4 個取得最優(yōu)的總體精度，通過上述分析，實驗的最優(yōu)參數(shù)選擇為h=5 像素，k= 4 個。

3.3 結(jié)果分析

從圖 5 中可以看出：常州市某區(qū) 2015—2018 年的主要變化體現(xiàn)在區(qū)域上部新修一條主干道，其他的一些典型變化列舉如圖 6 所示。

上述用地變化是否涉及違章建筑需要人工進一步比對審查，但在總體上，遙感影像快速變化檢測算法提供了可能發(fā)生違建的大致位置，節(jié)省了人力的比對。

4 結(jié)論

本研究為解決違章建筑的粗篩選問題，以常州市某區(qū)域遙感影像為例，改進了基于 PCA 和k-Means 聚類的變化檢測算法，通過選擇不同的分割像素塊邊長和聚類中心個數(shù)的超參數(shù)組合進行比對實驗。結(jié)果表明：在精度水平上該算法優(yōu)于傳統(tǒng)CVA 算法，并且在像素塊邊長為 5 像素、聚類中心個數(shù)為 4 個時，既能有效抑制噪聲，減小偽變化的出現(xiàn)，又能保證一定的召回率，使得整體精度達到理想效果，并在測試用例中篩選出池塘新修小道、池塘變?yōu)檗r(nóng)田等變化區(qū)域。

圖 6 典型變化示例