協同光學與雷達遙感數據的面向對象土地覆蓋變化檢測

翟鵬飛，李世華,2※，胡月明,3,4

（1. 電子科技大學資源與環境學院，成都 611731；2. 電子科技大學長三角研究院（湖州），湖州 313001；3. 海南大學熱帶作物學院，海口 570228；4. 廣州市華南自然資源科學技術研究院，廣州 510630）

0 引 言

社會的發展和人口的急劇增加導致人類對土地資源的需求也在不斷增長，土地利用/覆蓋（Land use and Land cover）變化成為全球或區域變化領域的研究熱點[1-5]。遙感技術的快速進步使得大面積數據獲取手段更加豐富，這對土地覆蓋變化檢測提出了更高的要求和挑戰。目前，國內外很多學者利用光學數據提出了不同的方法并將這些方法應用在像素級的變化檢測任務[3-4]。基于像元的方法是對前、后時相的同一位置的像元進行比較，通常采用差值或比值等代數運算獲取差異影像，然后使用不同的閾值分割方法來確定該像元是否發生變化[6-7]，這對前、后時相遙感影像的配準、輻射校正等預處理有著較高的要求，同時會使得檢測結果出現“椒鹽現象”。隨著影像空間分辨率逐漸提高，包括米級、亞米級數據的問世，同一地物在影像中往往由多個相鄰像元組成，單個像元的光譜特性不足以反映相應地物的光譜特征，而相鄰地物的像元多呈現一定的關聯性，同時，“同物異譜”或“同譜異物”的現象極大地影響了高分辨率變化檢測的精度。隨著基于地理對象的遙感影像分析方法的發展，針對對象的變化檢測方法以影像分割對象為基礎，以具有相似光譜和空間特征的影像對象作為基本處理單元，考慮了像元及其鄰域的光譜、幾何、紋理及空間上下文特征，能夠減少圖像配準結果的影響，極大地提高了變化檢測的精度，成為高分辨率光學遙感影像變化檢測的研究重點[8-12]。

結合合成孔徑雷達（Synthetic Aperture Radar，SAR）數據和面向對象方法在土地覆蓋變化檢測研究中取得了良好的效果。尤紅建[13]采用基于方差的區域增長算法實現多時相SAR 影像聯合分割，并計算斑塊的差異指數獲得變化檢測結果。張明哲等[14]使用簡單線性迭代分割方法對單極化Radarsat-2 數據進行分割，并融合像素和對象結果，通過眾數投票的方式獲得最終結果。Zhao 等[15]提出了利用圖像回歸的方式獲取變化差異圖，在超像素分割的基礎上利用C 均值聚類獲得檢測結果。以上研究中將影像分割算法直接應用在SAR 數據上，然而由于SAR傳感器與光學傳感器不同，SAR 成像是主動相干成像，成像機理復雜，且受斑點噪聲影響嚴重，各種圖像分割的方法在SAR 影像難以獲得良好的邊界，進而影響后續的檢測結果，這使得不同時相的SAR 影像的變化檢測更具有挑戰性。同時利用SAR 數據進行變化檢測多使用單極化或雙極化數據，也在一定程度上延緩了全極化SAR數據在土地覆蓋變化檢測應用的發展。西南地區常年多云多雨，高質量的光學影像數據通常難以獲取，而全極化SAR 傳感器重訪周期較長，二者單獨應用于土地覆蓋變化檢測任務難以滿足獲取高精度變化結果的要求。針對上述不足，提出一種協同光學與雷達遙感數據的面向對象土地覆蓋變化檢測的方法，該方法在一定程度上避免了單一使用兩種數據的缺陷，充分發揮了二者的優勢，為土地覆蓋變化檢測提供了一種新的范式。2 個樣區包括多種不同地物類型并且涵蓋了該地區多數的變化類型。所有數據都進行了相關的預處理并基于光學數據進行配準。為了驗證本文方法的可行性，選擇兩個試驗區進行抽樣檢驗，樣區大小均為488×488，如圖2和圖3 所示。試驗區數據如表1 所示。

表1 試驗區數據詳情Table 1 Details of the experimental area data

1 材料與方法

1.1 研究區概況與數據

研究區位于四川省眉山市（ 30°08′25.31″～30°18′42.39″N，103°42′0.92″～103°53′50.92″E）。眉山市地處四川盆地西南邊緣，成都平原的西南部，是成都平原連通川南、川西的咽喉。眉山市轄二區四縣，境內河網密集，主要有岷江和青衣江兩支水系。眉山市地處中低緯度，總體地勢西高東低，南高北低，屬中亞熱帶濕潤季風氣候，夏季溫潤多雨，多種植水稻等糧食作物。隨著經濟的發展，城市化建設的快速推進，眉山市大量水稻種植區改種席草、果林等經濟作物；部分裸地也隨著時間的推移施工完成，形成建設用地。

分別使用2019年Sentinel-2 光學數據和2016年、2019年全極化Radarsat-2 雷達數據，影像大小為2 000×2 000像元。圖1 中紅色和藍色框圖分別表示樣區一和樣區二。

1.2 研究方法

研究路線如圖4 所示，主要包括以下步驟。首先使用分型網絡演化分割算法對2019 年Sentinel-2 影像進行分割，得到初始分割結果，然后對前、后時相Radarsat-2數據利用對數比值法獲得差異影像。利用初始分割過程中獲取的矢量邊界作為基準，再次使用分型網絡演化分割算法對雷達差異影像進行引導分割獲取最終分割結果；選取合適的變化與未變化樣本并提取樣本特征空間，將樣本間距離度量可分性特征空間優化方法（Feature Space Optimization, FSO）結果與隨機森林（Random Forest，RF）分類器有機的結合構成變化檢測框架（以下稱為FSO-RF），利用該框架獲取最終的檢測結果。

1.2.1 遙感影像分割

根據分割原理的不同，圖像分割方法分為自頂向下的區域分裂和自底向上的區域生長法，本文采用的分割方法為Baatz 等[16]在2000 年提出的分型網絡演化分割算法（Fractal Net Evolution Approach, FNEA），該算法根據區域生長的思想，基于圖像對象的光譜和幾何特征，由小對象逐步合并為較大的對象，并保證在不同對象中保持較大的異質性。對象的異質性由光譜異質性和形狀異質性確定，而緊致度因子和光滑度因子共同確定形狀異質性。異質性由以下公式定義：

對象的光滑度用來表示對象分割邊界的破碎程度，利用對象邊界的周長l與近似邊界的長度b的比值來表示。

式中hcolor_m、hcompact_m、hsmooth_m分別表示對象合并后的光譜異質性、緊致度異質性和光滑度異質性，nobj1、nobj2、lobj1、lobj2、bobj1和bobj2分別表示合并前2 個影像對象的像元個數、對象邊界周長以及對象的近似邊界長度；nmerge、lmerge、bmerge分別表示影像合并后的像元個數、合并后對象邊界的周長以及合并后對象的近似邊界的長度。

分割尺度也是分割算法的重要參數，尺度是一個抽象的概念，它控制了分割過程中形成對象的大小和最終分割獲得的對象個數，尺度參數不同，對象的大小和個數也不同。隨著尺度參數的增加，分割獲得的對象越大，最終獲得的對象個數相應越少。本文目的在于分割獲得的對象內部應僅包含變化或未變化像元，既要保證良好的分割邊界，又不要使得分割結果過于破碎。

遙感影像分割的理想結果是分割后，對象內部具有較高的同質性，而相鄰的不同對象間呈現良好的異質性[17]，并且已有研究表明，基于對象的方法能更好的抑制SAR 影像相干斑噪聲對圖像信息的破壞[18]。本文使用對數比值法獲得雷達差異影像，該方法相較差值法、比值法能夠更好的抑制噪聲[19-21]。利用分型網絡演化分割算法對Sentinel-2進行分割，分割的結果由波段權重、尺度參數、形狀參數和緊致度構成。由于所選的研究區土地覆蓋類型復雜多樣，地物十分破碎，經過多次試驗將Sentinel-2 影像的各波段權重設置為1，分割尺度設置為50，形狀參數設置為0.2，緊致度因子設置為0.3，并以初始分割結果最為基準引導雷達差異影像進行二次分割獲得最終的分割結果。

1.2.2 距離度量可分性特征空間優化

遙感影像經過分割形成對象后能夠提供更加豐富的特征，包括光譜、紋理、形狀、空間上下文等特征，可以充分利用這些特征來提高變化檢測的精度。由于特征信息的冗余等因素，當特征的數量增加到一定程度時，變化檢測的精度不僅不會繼續增加，甚至有可能降低[10]，因此有必要對提取的特征空間進行優化，以獲取最佳的特征集合進行變化檢測任務。距離度量可分性特征空間優化方法的步驟和基本思想為選取可靠的變化、未變化樣本并提取相應特征構建特征空間，其次在各個特征子空間中計算變化、未變化的樣本可分距離，理論上可分距離越大則表明二者差異越大，可分性越好。對每一個維度所有樣本的距離進行比較，獲得每個維度的最大值，代表該維度空間下的最佳可分距離，然后比較每個維度的最大值，這些值中的最大可分距離，其所代表的特征組合即為最優的特征組合。變化與未變化類的可分距離由下式定義：

式中dmax表示在最優特征空間下的最大可分距離。

使用距離度量可分性特征空間優化的方法來尋找變化檢測任務的最優特征空間。本文依據HH、HV、VV 三種極化方式，Freeman-Durden 分解獲得的雙次散射（double-bounce scattering，dbl）、體散射（volume scattering，vol）和表面散射（surface scattering，surf）參數以及Pauli 分解獲得的Pauli_r、Pauli_g 和Pauli_b 參數來計算后續所要提取的特征。在選取置信度較高的變化與未變化樣本后，提取樣本的特征，包括HH、HV、VV后向散射均值和標準差，Pauli 分解參數的對象均值和標準差以及Freeman-Durden 分解參數的對象均值和對象標準差，同時根據Haralick 提出的灰度共生矩陣（Gray-level Co-occurrence Matrix,GLCM）提取紋理特征[22]，本文利用同質性（Homogeneity）、對比度（Contrast）、熵（Entropy）、均值（Mean）、標準差（Standard Deviation）和相關性（Correlation）6 種紋理特征。利用灰度共生矩陣可以按照不同方向進行掃描求得對應參數的值，包括0°、45°、90°以及135°，本文按照所有方向掃描并求取均值的方式求得對應紋理特征的值，共提取72 維特征如表2 所示。

表2 特征參數變量Table 2 Characteristic parameter variables

根據式（8）～（10）計算各維度特征子空間的變化、未變化樣本的可分距離，并在各個維度中獲取該維度的最大值，這些值中的最大可分距離所代表的特征組合即代表最佳特征組合。

1.2.3 隨機森林面向對象變化檢測方法

隨機森林（Random Forest）是近年發展起來的一種機器學習模型，由Breiman[23]首先提出。由于其可以適用于高維特征樣本和大的數據集上，因此已經廣泛應用于遙感領域，并且在分類和土地覆蓋變化檢測的任務上獲得了良好的效果[10,11,24]。隨機森林的本質屬于集成學習，該模型的理論基礎是決策樹，它由不同的決策樹進行組合，最終的結果根據每棵決策樹生成的結果進行投票產生最優的分類結果。在隨機森林中，誤差估計是在訓練過程中進行的，訓練數據根據Bootstrap 重采樣獲得不同的訓練集，利用這些不同的訓練集分別構造決策樹；在采樣過程中使用Bagging 的原理，同一個樣本可以被多次選擇，而其他樣本可能沒有被選擇，大約由三分之一的樣本沒有被用于每個決策樹的構建，這些樣本用來測試以獲得袋外誤差（Out of Bag ,OOB）估計。

在初始樣本選擇時，規定以下樣本選擇原則：假設在分割對象內部共有Nn個像元，其中屬于變化的像元個數為Nc，屬于未變化的像元為Nu，顯然Nn=Nc+Nu。定義對象內部變化率μ如下：

當μ=0，表示對象內部均為未變化像元，則該對象選擇為未變化樣本，當μ=1，則表示對象內部均為變化樣本，相應的選擇為變化樣本。當μ>80%或μ<20%時，在該置信區間選擇置信度較高的變化與未變化樣本。圖5 為置信度較高的變化與未變化樣本示意圖。

1.2.4 檢驗方法

為了驗證所提方法的有效性，有必要對變化檢測結果進行定量的精度驗證，本文在引入混淆矩陣后采用準確率（Accuracy）、召回率（Recall）、F1、錯分誤差（Commission Errors, CE）和漏分誤差（Omission Errors,OE）進行精度評價。

準確率表示正、負類別均被正確分類的概率。

式中TP 表示將正類預測為正類數量；TN 表示將負類預測為負類數量；FP 表示將負類預測為正類數量；FN 表示將正類預測為負類的數量。

2 結果與分析

2.1 面向對象分割

以樣區一為例進行對比試驗，將分割算法直接應用于雷達差異影像，并將分割結果套合于光學影像用于可視化對比分析如圖6 所示。圖6a、6b 分別表示對光學數據分割和利用光學數據引導SAR 影像分割；圖6c 是直接將分割算法應用于雷達差異圖，將分割結果套合于光學影像如圖6d。

將圖6a 中紅色框圖內區域的分割細節作為展示如圖7，重點關注圖7 中白色框圖內的地物，其內部地物類型由草地變化為建設用地。根據分割算法獲得的分割結果來看，無論是直接對SAR 影像分割還是利用光學影像引導SAR 影像分割都可以正確提取變化的范圍，但是后者可以在雷達差異圖中獲得良好的地物邊界，獲取更加精細的變化地理實體，并且針對變化劇烈的區域分割對象內部獲取了盡可能多的變化或未變化像元。

2.2 FSO-RF 變化檢測框架結果分析

在選取置信度較高的變化和未變化樣本后，依據選的樣本計算樣本各維度特征的可分距離如表3 所示，表中顯示了各維度樣本可分距離大小和具體的特征參數，以一維、二維特征為例，樣本特征一維空間中GLCM Con_surf 特征使得72 個一維特征的變化、未變化樣本的可分距離最大，因此它的距離代表了一維特征空間的可分距離。樣本特征二維空間中GLCM Con_surf 與GLCM Mea_dbl 特征在各個二維特征組合中獲得了最大的可分距離，因此它代表了二維特征空間的可分距離，其余同理。最終在第45 維空間中獲得最大可分距離，理論上該距離所代表的特征空間即為本次任務的最佳特征組合。

表3 各維度特征可分距離Table 3 Separable distances for each dimensional feature

圖8 直觀地展示了樣本特征距離可分性大小隨著特征數量增加的變化趨勢，隨著特征的增加，樣本間的可分距離先增加，并在第45 維特征空間獲得了最大的可分距離，隨后可分距離減小，在一定程度上表明了更多的特征不一定全部有益于變化檢測，也凸顯了特征優選的必要性。

為了驗證距離度量可分性特征空間優化方法的有效性并實現FSO-RF 變化檢測框架，本文將提取的特征按照維度由小到大的方式根據樣本間可分距離的計算結果正序的將各維度的特征輸入到隨機森林分類器中，獲取各維度特征下變化檢測結果的準確率如圖9 所示。準確率的總體變化趨勢與不同維度的可分距離的計算結果基本吻合，同時在第 41 維特征獲得了最高的準確率為92.90%，此時的樣本類間可分距離為1.653，在一定程度上證明距離度量可分性特征空間優化方法的有效性。

2.3 精度評價與討論

本文采用像素的對數比值法獲得雷達影像差異圖，然后根據分型網絡演化分割方法對光學影像進行分割，并利用該分割結果引導雷達差異圖進行分割獲得最終分割結果。選取置信度較高的變化與未變化樣本并提取樣本特征，利用本文提出的FSO-RF 框架獲得最終的變化檢測結果。

將本文的方法與基于雷達差異影像分割的方法、變化向量分析（Change Vector Analysis, CVA）、主成分分析（Principal Component Analysis, PCA）、多元變化檢測（Multivariate Alteration Detection, MAD）、迭代多元變化檢測（ Iteration Multivariate Alteration Detection,IR-MAD）、以及基于像元的支持向量機（Support Vector Machines, SVM）和隨機森林（Random Forest, RF）的方法進行對比分析。其中基于雷達差異影像分割的方法將分型網絡演化分割方法直接應用于雷達差異影像，尺度參數設置為7，形狀因子和緊致度因子分別設置為0.1 和0.4，根據1.2.3 節樣本選擇原則選擇置信樣本并提取1.2.2節所示特征后，根據本文提出的FSO-RF 變化檢測流程框架獲得變化檢測結果；CVA 方法是根據計算前、后時相SAR 影像的特征的歐式距離獲得變化結果；PCA 方法是在對特征進行降維，選取貢獻度前75%的特征，然后利用CVA 獲得變化檢測結果；IR-MAD 是在MAD 的基礎上改進，根據卡方距離進行加權迭代，最大迭代次數Max Iteration 設置為100，變化收斂閾值Delta Criterion 設置為0.001；隨機森林和支持向量機方法是在選擇像素級的變化與未變化樣本進行分類獲得變化檢測結果。

同時，為了更直觀地顯示變化檢測任務的結果，對樣區一和樣區二進行定性和定量的分析。圖10 和圖11分別展示了2 個樣區的變化檢測結果。各種方法的精度評價結果如表4 所示。圖10a 和11a 是根據高分辨率遙感影像進行目視解譯結合實地考察繪制的參考變化結果。

表4 不同方法下變化檢測精度評價Table 4 Accuracy evaluation of change detection between different methods %

從整體以及2 個樣區的變化檢測結果和精度可以發現：

1）無論是傳統經典的基于像元的變化檢測方法還是有監督的SVM 和RF 方法，在研究區數據的變化檢測結果精度都相對較低。這些方法在檢測過程中只使用了前、后時相影像的后向散射特征和兩種極化分解特征；更重要的原因在于即使采用了濾波的處理手段，SAR 影像的相干噪聲也不能完全去除，這些斑點噪聲在前、后時相的差異計算中以相應特征的形式干擾算法的計算，導致針對SAR 影像基于像元的方法會造成大量的椒鹽噪聲，對于變化劇烈的區域也難以獲得良好的邊界。同時也發現，改進后的IR-MAD 相較MAD 或者其他兩種無監督的變化檢測方法在檢測結果上會減少大量的孤點；總體而言，有監督的SVM 和RF 對無監督的經典變化檢測方法在精度上有不同幅度的提升。

2）基于雷達差異影像分割的方法在各項精度評價指標下大幅度提高了變化檢測的精度，同時，利用對象分割可以極大程度的改善基于像元方法的椒鹽噪聲現象。然而，由于雷達成像原理與光學傳感器成像有著本質的區別，分型網絡演化分割方法不能分割出地物的邊界，進而難以獲得良好的變化檢測結果。根據圖10h、圖11h也可以發現，檢測結果僅能粗略地提取變化的大致范圍，而難以獲取細致、真實地理實體的變化結果。

3）2 個樣區結果在精度上有一定差異，樣區一和樣區二的準確率分別為95.08%和88.16%，原因在于初始分割僅使用了后時相光學數據。針對樣區一，前后時相的土地類型相差較大，多為建筑裸地或荒草地向建設用地轉變，因此分割形成了不同的土地覆蓋類型對象，進而分割對象內包含了變化或未變化像元；針對樣區二，精度相對較低的原因在于不同類型的農用地的轉變過程，樣區二后時相農用地多種植為果林或成片的席草等經濟作物，在使用后時相光學數據分割時，盡管已經考慮到目的是為了獲取的對象內部應僅包含變化或未變化像元，但是為保證對象內部的同質性以及盡量不將地物分割的過于破碎，將成片狀的經濟作物分割為同一對象，而對于前一時相，該對象內部存在其他地物類型，如水稻等糧食作物，因此分割對象內部未獲得盡可能多的變化或未變化像元導致在樣區二精度低于樣區一。

4）本文使用距離度量可分性特征空間優化方法對提取的72 維特征進行優化并在第45 維特征中獲得最大可分距離，在驗證的過程中表明第41 維特征中獲得了最佳的變化檢測結果。但是利用高維度特征在實際的變化檢測應用中復雜且繁瑣，而且本文方法可以計算出具體的可分距離的大小，并且在文中驗證了可分距離與準確率的變化曲線基本吻合，因此可以利用可分距離大小的變化曲線，選取可分距離大小排名靠前的對應的特征或者設定當可分距離的增長率第一次小于某閾值時，選擇當前維度特征進行變化檢測來指導實際的變化檢測應用。

5）本文將面向對象的思想與隨機森林算法相結合，利用分型網絡演化分割的方法對光學影像分割并引導分割雷達差異影像以獲取擁有良好、真實的地理邊界分割結果；同時，針對提取的多維特征根據各維度特征樣本的特征距離進行特征空間優化，驗證了距離度量可分性特征空間優化的方法進而獲得良好的變化檢測結果，各項精度指標均高于對比方法，準確率為92.90%，召回率為96.61%，F1 為96.07，錯分誤差和漏分誤差分別為4.47%和3.39%，進一步證明了該方法的有效性和可行性。

綜上所述，在多云多雨地區應用有限的光學數據引導雷達影像分割進而采用隨機森林算法進行面向對象的土地覆蓋變化檢測研究是可行的，這種方法不僅能有效地發揮光學數據豐富光譜信息的優勢，在引導雷達影像分割時能夠獲得良好的真實地物邊界，同時能將SAR 影像所提供豐富的極化信息有機的相結合來獲取良好的變化檢測結果。

3 結 論

本文提出了一種基于隨機森林算法的FSO-RF 土地覆蓋變化檢測算法框架，該框架以光學影像引導SAR 數據分割為基礎獲取良好的真實地物邊界，在選擇置信樣本提取多維特征空間后，驗證了變化檢測結果精度與不同維度的特征之間的關系，并將距離度量可分性特征空間優化方法和隨機森林算法相結合構建FSO-RF 變化檢測算法框架，該算法框架能夠有效、準確地提取變化區域，在實驗區的總體精度達到92.90%，錯分誤差和漏分誤差分別為4.47%和3.39%，相較于現有變化檢測算法的精度有一定的提升。