とらりもんHOME  Index  Search  Changes  Login

とらりもん - 土地被覆機械学習の手法間比較 Diff

  • Added parts are displayed like this.
  • Deleted parts are displayed like this.

(2016.05 水落裕樹; 2017.02追記)
{{toc}}
{{br}}古典的な機械学習手法をいくつか試してみたので、結果と解析方法を書きます。機械学習にはR、GISにはGRASSを使ってみました(それぞれの使い方は解説しません)。
(2019.01)
環境設定等も含めて加筆修正したコンテンツを[[こちらに|https://sites.google.com/site/mizuochipublic/%E5%AE%9F%E8%B7%B5%E3%82%B3%E3%83%B3%E3%83%86%E3%83%B3%E3%83%84-remote-sensing-tutorials/%E6%A9%9F%E6%A2%B0%E5%AD%A6%E7%BF%92%E3%81%AB%E3%82%88%E3%82%8B%E5%9C%9F%E5%9C%B0%E8%A2%AB%E8%A6%86%E5%88%86%E9%A1%9E]]作成中です。

また、あくまで解析方法の紹介と、手法間の比較が目的のため、最終プロダクトの精度にはこだわってません。
{{br}}なお、リモセン解析や機械学習の理論的・手法的なマチガイにお気づきの方は、ぜひご教示ください。

!衛星データと前処理
つくば界隈(140.0E-140.5E,36.0N-36.5N)の土地被覆分類を公開衛星データでやってみましょう。
USGSの[[Landsat8の地表面反射率プロダクト|http://landsat.usgs.gov/about_LU_Special_Issue_3.php]]、JAXAの公開している[[PRISM DEM|http://www.eorc.jaxa.jp/ALOS/aw3d/index.htm]]を入力データとして使ってみます。
[[ここ|http://pen.agbi.tsukuba.ac.jp/~ryuiki_plus/sample_data/Landsat/LULC]]|http://pen.envr.tsukuba.ac.jp/~ryuiki_plus/sample_data/Landsat/LULC]]から以下のものをDLし、tarボールは展開して下さい。
landsat8.tar.gz #Landsat8
N036E140_AVE_DSM.tif #DEMデータ
N035E140_AVE_DSM.tif
N036E139_AVE_DSM.tif
N035E139_AVE_DSM.tif
tsukuba_saclaj_extract.csv #教師&検証データ
landsat8データは、earthexplorerからDLした晴天時の地表面反射率プロダクト
LC81070352015090-SC20160520024926.tar.gz #2015年DOY090のLandsat8
LC81070352015298-SC20160520025130.tar.gz #2015年DOY298のLandsat8
を、正距円筒図法に投影法変換し、更に、国土基盤情報と重ねて幾何補正したものです(何故か元データが大きく位置ずれしてました。landsat8を使う際は最初に位置ずれしてないか確認したほうがいいかも)。
最後のデータは、つくば周辺のGIS検証写真合計4356点のデータを次の10カテゴリに整理したものです:
DF(落葉林), EF(常緑林), bamboo, barren, croplands, grasslands, orchard, rice, urban, water
カテゴリ間でデータサイズが相当ばらついてますが、そのあたりの検討はスルーです。

もし上記フォルダにアクセスできない場合は、USGSおよびJAXAの本家をあたってください。教師データは[[github|https://github.com/hmizuochi/LULC/blob/master/tsukuba_saclaj_extract.csv]]にも置いてあります。

{{br}}GRASS(EPSG:4326)の環境を作成し、データを読み込んでください。
r.in.gdal input=N036E140_AVE_DSM.tif output=N036E140_AVE_DSM
r.in.gdal input=N035E140_AVE_DSM.tif output=N035E140_AVE_DSM
r.in.gdal input=N036E139_AVE_DSM.tif output=N036E139_AVE_DSM
r.in.gdal input=N035E139_AVE_DSM.tif output=N035E139_AVE_DSM
for NAME in LC81070352015090LGN00 LC81070352015298LGN00;do
     NAME2=`echo $NAME | cut -c 1-16`
     for i in `seq 1 7`;do
         r.in.gdal input=${NAME}_sr_band${i}_ll_modified.tif output=${NAME2}_sr_band${i}_ll_modified --o
     done
done
regionの設定(140.0E-140.5E, 36.0N-36.5N, 分解能はLandsatデータにあわせる)をしてください。できたら、取り敢えずtrue color画像を見てみましょう。
for i in 2 3 4;do r.colors map=LC81070352015090_sr_band${i}_ll_modified color=grey.eq; done
d.rgb red=LC81070352015090_sr_band4_ll_modified green=LC81070352015090_sr_band3_ll_modified blue=LC81070352015090_sr_band2_ll_modified
*true color画像(みえない場合は[[ここからDL|https://sites.google.com/site/mizuochipublic/_/rsrc/1486745186539/yan-jiu-nei-rong-research-themes/landsat.png?height=200&width=200]])
[[true color画像|http://pen.agbi.tsukuba.ac.jp/~ryuiki_plus/sample_data/Landsat/LULC/landsat.png]]|http://pen.envr.tsukuba.ac.jp/~ryuiki_plus/sample_data/Landsat/LULC/landsat.png]]
{{br}}筑波山と霞ヶ浦あたりを範囲にいれておきました。
{{br}}では各種特徴量を計算しましょう。まずLandsat8の2015年DOY090データから
band1〜band7の反射率、NDVI、NDWI、NDVIのテクスチャ値(asm,contrast,entropy,variance)
を計算します。
for i in `seq 1 7`;do
     r.mapcalc "LC81070352015090_sr_band${i}_C=LC81070352015090_sr_band${i}_ll_modified/10000.0"
     r.mapcalc "LC81070352015090_sr_band${i}_C=if(LC81070352015090_sr_band${i}_C>=1.0,1.0,LC81070352015090_sr_band${i}_C)"
     r.mapcalc "LC81070352015090_sr_band${i}_C=if(LC81070352015090_sr_band${i}_C<=0,0,LC81070352015090_sr_band${i}_C)"
done
r.mapcalc "LC81070352015090_NDVI=(LC81070352015090_sr_band5_C-LC81070352015090_sr_band4_C)/(LC81070352015090_sr_band5_C+LC81070352015090_sr_band4_C)"
r.mapcalc "LC81070352015090_NDWI=(LC81070352015090_sr_band3_C-LC81070352015090_sr_band7_C)/(LC81070352015090_sr_band3_C+LC81070352015090_sr_band7_C)"
r.mapcalc "LC81070352015090_NDVI_rescale=int((255-0)*(LC81070352015090_NDVI-(-1.0))/(1.0-(-1.0))+0)" #テクスチャ計算前に、0〜255の離散データにスケール変更
r.texture input=LC81070352015090_NDVI_rescale prefix=LC81070352015090_NDVI_asm -a
r.texture input=LC81070352015090_NDVI_rescale prefix=LC81070352015090_NDVI_contrast -c
r.texture input=LC81070352015090_NDVI_rescale prefix=LC81070352015090_NDVI_entropy -e
r.texture input=LC81070352015090_NDVI_rescale prefix=LC81070352015090_NDVI_variance -v
つづいて、DEMデータから傾斜と方位角を計算します。regionの縁がDEMデータの縁にちょうどかぶっているため、まずちょっと広めのregionにしてモザイクしてから処理しましょう。
g.region save=tsukuba.rgn
d.zoom #ミドルクリックでアンズーム
r.mapcalc "N036E140_AVE_DSM_dummy=if(isnull(N036E140_AVE_DSM),-1,N036E140_AVE_DSM)"
r.mapcalc "N035E140_AVE_DSM_dummy=if(isnull(N035E140_AVE_DSM),-1,N035E140_AVE_DSM)"
r.mapcalc "N035E139_AVE_DSM_dummy=if(isnull(N035E139_AVE_DSM),-1,N035E139_AVE_DSM)"
r.mapcalc "N036E139_AVE_DSM_dummy=if(isnull(N036E139_AVE_DSM),-1,N036E139_AVE_DSM)"
r.mapcalc "DSM_mosaic=max(N036E140_AVE_DSM_dummy,N035E140_AVE_DSM_dummy,N035E139_AVE_DSM_dummy,N036E139_AVE_DSM_dummy)"
r.colors map=DSM_mosaic color=srtm
r.slope.aspect elevation=DSM_mosaic slope=DSM_mosaic_slope aspect=DSM_mosaic_aspect
g.region tsukuba.rgn
最後に、Landsat8の2015年DOY298とDOY090の差分画像も計算します。全特徴量についてやってもいいのですが、今回はNDVI,NDWIだけにしました(DOY298のNDVI,NDWIも上記にならって予め計算しておいてください)。
r.mapcalc "NDVI_diff=LC81070352015298_NDVI-LC81070352015090_NDVI"
r.mapcalc "NDWI_diff=LC81070352015298_NDWI-LC81070352015090_NDWI"
では、教師データと特徴量の対応表を作りましょう。機械学習的にはこの段階で第一の勝負が決まってる感があります。
echo "B1,B2,B3,B4,B5,B6,B7,DEM,NDVI,NDWI,SLOPE,ASM,CONTRAST,ENTROPY,VARIANCE,NDVIDF,NDWIDF,CAT" > mldata.csv
for l in `seq 1 4356`;do
     N=`cat tsukuba_saclaj_extract.csv | head -n ${l} | tail -n 1 | cut -d"," -f 4`
     E=`cat tsukuba_saclaj_extract.csv | head -n ${l} | tail -n 1 | cut -d"," -f 3`
     CAT=`cat tsukuba_saclaj_extract.csv | head -n ${l} | tail -n 1 | cut -d"," -f 1`
     B1=`r.what input=LC81070352015090_sr_band1_C east_north=${N},${E} | cut -d"|" -f 4`
     B2=`r.what input=LC81070352015090_sr_band2_C east_north=${N},${E} | cut -d"|" -f 4`
     B3=`r.what input=LC81070352015090_sr_band3_C east_north=${N},${E} | cut -d"|" -f 4`
     B4=`r.what input=LC81070352015090_sr_band4_C east_north=${N},${E} | cut -d"|" -f 4`
     B5=`r.what input=LC81070352015090_sr_band5_C east_north=${N},${E} | cut -d"|" -f 4`
     B6=`r.what input=LC81070352015090_sr_band6_C east_north=${N},${E} | cut -d"|" -f 4`
     B7=`r.what input=LC81070352015090_sr_band7_C east_north=${N},${E} | cut -d"|" -f 4`
     DEM=`r.what input=DSM_mosaic east_north=${N},${E} | cut -d"|" -f 4`
     NDVI=`r.what input=LC81070352015090_NDVI east_north=${N},${E} | cut -d"|" -f 4`
     NDWI=`r.what input=LC81070352015090_NDWI east_north=${N},${E} | cut -d"|" -f 4`
     SLOPE=`r.what input=DSM_mosaic_slope east_north=${N},${E} | cut -d"|" -f 4`
     ASM=`r.what input=LC81070352015090_NDVI_asm_ASM_0 east_north=${N},${E} | cut -d"|" -f 4`
     CONTRAST=`r.what input=LC81070352015090_NDVI_contrast_Contr_0 east_north=${N},${E} | cut -d"|" -f 4`
     ENTROPY=`r.what input=LC81070352015090_NDVI_entropy_Entr_0 east_north=${N},${E} | cut -d"|" -f 4`
     VARIANCE=`r.what input=LC81070352015090_NDVI_variance_Var_0 east_north=${N},${E} | cut -d"|" -f 4`
     NDVIDF=`r.what input=NDVI_diff east_north=${N},${E} | cut -d"|" -f 4`
     NDWIDF=`r.what input=NDWI_diff east_north=${N},${E} | cut -d"|" -f 4`
     echo "$B1,$B2,$B3,$B4,$B5,$B6,$B7,$DEM,$NDVI,$NDWI,$SLOPE,$ASM,$CONTRAST,$ENTROPY,$VARIANCE,$NDVIDF,$NDWIDF,$CAT" >> mldata.csv
done
指定した座標の画像値を取得するr.whatコマンドをfor文で回しているのですが(5分かかった)、スキルが足りない感が否めません。
{{br}}なお、計算しておいてなんですが、方位角 (ASPECT) はあまり関係ないかな?と思ったので、特徴ベクトルからは外しました。

!Rの準備
{{br}}機械学習を始める前に、GRASSで吐き出した対応表を整形しましょう。
cat mldata.csv | sed "s/DF/1/g" | sed "s/EF/2/g" | sed "s/bamboo/3/g" | sed "s/barren/4/g" | \
sed "s/croplands/5/g" | sed "s/grasslands/6/g" | sed "s/orchard/7/g" | sed "s/rice/8/g" | \
sed "s/urban/9/g" | sed "s/water/10/g" > mldata_d.csv #カテゴリ名が文字列だと面倒なので番号に変更
また、幾つかのライブラリを使うので、インストールしておいてください。以下は筆者の環境での例です(ubuntu 14.04)。
sudo apt-get install r-cran-car
R #GRASS上で実行すること。
> install.packages("caret")
> install.packages("VGAM")
> install.packages("e1071")
> install.packages("randomForest")
続いて、Rにデータを読み込みます。
> mldata<-read.csv("mldata_d.csv")
> mldata.scale<-data.frame(scale(mldata[,-18]),CAT=as.factor(mldata$CAT)) #一応スケーリングするとともに、CAT列をfactor型に変更
ちゃんと読み込めたかどうか、head(mldata)やdim(mldata)などで確認してください。
{{br}}機械学習においては、教師データを学習・検証・チューニングにどう割り振るか考える必要がありますが、今回は[[交差検証|http://musashi.osdn.jp/tutorial/mining/xtclassify/accuracy.html]]を利用しようと思います。交差検証ではデータセットを何分割するか決める必要があります(データサイズぶん分割するLOOというのもありますが、さすがに時間かかりそうなのでやめます)。それは、[[とあるページによれば|http://d.hatena.ne.jp/hoxo_m/20110324/p1]]、データサイズをn、分割数をkとすると
k=1+log(n)/log(2)
とするのがいいそうです(スタージェスの公式)。交差検証はRではcaretパッケージで実現できます。
fold<-1+log(nrow(mldata.scale),2) #確認してみると13.09、つまり13分割がよさそう
group<-createFolds(mldata.scale$CAT,k=13)
これでようやく機械学習の準備が整いました。余談ですが、昨今は第三次AIブームと言われているそうですが、人工「知能」というからには、このへんの(ある意味)面倒くさいデータ集め&前処理こそ機械自らやってほしい、というかやってくれない限り、(巷で懸念されているように)人間の仕事が無くなるということはないのでは、と思います。

!多項ロジスティック回帰
2値ロジスティック回帰を多クラス分類に一般化したものです。詳細は触れませんが、回帰曲線をロジスティック関数ではなくソフトマックス関数に変えることで、特徴ベクトルから各クラスに属する事後確率を計算し、事後確率最大のクラスに分類します。基本的に線形な識別器なので、衛星データにどこまで適応できるかは心許ないです。
{{br}}まあ、とりあえずやってみましょう。VGAMライブラリを使います。
library("VGAM")
OA.ave<-0 #全体精度を格納
for(i in 1:fold){
    training<-mldata.scale[-group[[i]],]
    validation<-mldata.scale[group[[i]],]
    mldata.scale.vglm<-vglm(CAT~.,training,family=multinomial) #VGAMパッケージのvglmで多項ロジット回帰
    out.vglm<-predict(mldata.scale.vglm,newdata=validation[,-18],type="response") #validationデータの特徴ベクトルで事後確率予測
    classify<-apply(out.vglm,1,which.max) #クラス分類
    classify<-factor(classify,levels=c("1","2","3","4","5","6","7","8","9","10")) #とちゅうで勝手にnumericになるので、factor型に読み直し
    errormat<-table(classify,validation$CAT)
    OA<-sum(diag(errormat))/sum(errormat)
    OA.ave<-OA.ave+OA
}
OA.ave<-OA.ave/fold
print("multilogit OA")
print(OA.ave)
8分くらいかかりました。交差検証による全体精度は0.75。イマイチです。ちなみに最後のループのエラーマトリクスを貼り付けておきます。
classify  1  2  3  4  5  6  7  8  9 10
       1  23  2  0  0  0  2  0  0  0  0
       2   0 37  2  0  1  1  0  0  0  0
       3   1  0  7  0  0  0  2  0  0  0
       4   0  0  0  4  0  0  0  0  1  0
       5   2  3  0  1 58 10  4 10  3  0
       6   0  0  0  0  4 17  0  3  1  0
       7   0  0  0  0  0  0  0  0  0  0
       8   0  0  0  0 20  3  0 80  3  0
       9   1  0  0  0  0  0  0  2 19  0
       10  0  0  0  0  1  0  0  0  0  8
左端の列が作成者の分類結果を表しています。これをみると、cropland(5)をrice(8)に分類、rice(8)をcropland(5)に分類、grassland(6)をcropland(5)に分類、常緑林(2)をcropland(5)に分類、といった誤分類が目立ちます。また果樹園(7)に至っては、その多くがcroplandに誤分類された結果、一つも正しく分類できていません。croplandが鬼門のようですね。ちなみに番号とカテゴリの対応は
DF/1 /EF/2 bamboo/3 barren/4 croplands/5 grasslands/6 orchard/7 rice/8 urban/9 water/10
となっています。

!ニューラルネットワーク
隠れ層を追加した[[パーセプトロン|http://hokuts.com/2015/11/25/ml2_perceptron/]]です。出力層でソフトマックス関数を用いることで、先ほどと同様、事後確率が明示的に計算できます。Rだとcaretライブラリの中に入ってます。
{{br}}今回は、ハイパーパラメータのチューニングも行うことにします。caretパッケージのニューラルネットワークでは、隠れ層のユニット数sizeと、重みの正則化パラメータdecayが指定でき、これをいじるだけでも結果がだいぶ変わったりするためです(汎化を優先するか、フィッティングを優先するか)。そこで、これ以降は、13分割したデータセットのうち、1つのグループをチューニングに用い、残りの12個のグループで交差検証をする、という考え方で進めてみます。
tuning<-mldata.scale[group[[1]],]
nnet.tune<-train(CAT~.,data=tuning,method="nnet",maxit=100,tuneLength=4,trace=F) #チューニングに数分かかった
print(nnet.tune)
s<-nnet.tune$finalModel$tuneValue[[1]] #隠れ層のユニットサイズ
d<-nnet.tune$finalModel$tuneValue[[2]] #重みの正則化パラメータ
確認してみると、size=7,decay=0.00316がいいそうですので、そうしてみます。
OA.ave<-0
for(i in 2:fold){
    training<-mldata.scale[-c(group[[1]],group[[i]]),]
    validation<-mldata.scale[group[[i]],]
    mldata.scale.nnet<-nnet(CAT~.,training,size=s,decay=d) #ちなみにmaxitでiterationの回数を変えることもできる
    out.nnet<-predict(mldata.scale.nnet,validation[,-18])
    classify<-apply(out.nnet,1,which.max)
    classify<-factor(classify,levels=c("1","2","3","4","5","6","7","8","9","10")) #しょっちゅうnumericになるので、factor型に読み直し
    errormat<-table(classify,validation$CAT)
    OA<-sum(diag(errormat))/sum(errormat)
    OA.ave<-OA.ave+OA
}
OA.ave<-OA.ave/(fold-1.0)
print("neuralnet OA")
print(OA.ave)
精度は…0.75でした。精度だけ見れば多項ロジットとあまり変わらないですね。エラーマトリクスを貼り付けておきます。
classify  1  2  3  4  5  6  7  8  9 10
       1  23  3  0  0  0  2  0  0  0  0
       2   0 34  0  0  1  0  0  0  0  0
       3   1  2  8  0  1  0  0  0  0  0
       4   0  0  0  0  0  0  0  0  0  0
       5   1  2  0  2 62 10  3 11  2  1
       6   0  1  1  0  2 19  2  1  1  0
       7   1  0  0  0  0  0  1  0  0  0
       8   0  0  0  0 17  2  0 79  5  0
       9   0  0  0  3  0  0  0  4 19  0
       10  1  0  0  0  1  0  0  0  0  7

!ランダムフォレスト
弱識別器として決定木を用いる[[バギング|https://www1.doshisha.ac.jp/~mjin/R/32/32.html]]です。軽くて汎化性能のすぐれた手法として、リモセン界でも人気があります。決定木の各分岐において、各特徴量が平均的に不純度をどのくらい下げたかを指標として、変数重要度なども出してくれます。RだとずばりrandomForestライブラリというのがありますのでそれを使います。
{{br}}なお、ランダムフォレストが既にして交差検証に似た(ブートストラップ)枠組みで動いているので、更に交差検証をする必要もないのかもしれませんが、他の識別器と比較手法を合わせるという意味で、先ほどまでの検証の枠組みを踏襲します。
{{br}}ランダムフォレストのチューニングパラメータは、決定木に与える特徴量の数です(分岐のたびに、全特徴量のなかからランダムに幾つか候補を選ぶ。ゆえにランダムフォレスト)。[[ものの本|http://www.amazon.co.jp/%E3%81%AF%E3%81%98%E3%82%81%E3%81%A6%E3%81%AE%E3%83%91%E3%82%BF%E3%83%BC%E3%83%B3%E8%AA%8D%E8%AD%98-%E5%B9%B3%E4%BA%95-%E6%9C%89%E4%B8%89/dp/4627849710]]によれば、特徴ベクトルの次元(全特徴量の数)をdとするとき、√dとするのがよい、とのことです。しかしせっかくなのでチューニングしてみましょう。
> library("randomForest")
> tuneRF(tuning[,-18],tuning[,18],doBest=T)
Call:
randomForest(x = x, y = y, mtry = res[which.min(res[, 2]), 1])
               Type of random forest: classification
                     Number of trees: 500
No. of variables tried at each split: 4
        OOB estimate of  error rate: 26.29%
というわけで与える特徴量の数(mtry)は4とするのがよさそうです。特徴ベクトルの次元が17ですから、√17≒4で経験的にもそのぐらいが良さそうですね。
OA.ave<-0
for(i in 2:fold){
    training<-mldata.scale[-c(group[[1]],group[[i]]),]
    validation<-mldata.scale[group[[i]],]
    mldata.scale.rf<-randomForest(CAT~.,training,mtry=4)
    out.rf<-predict(mldata.scale.rf,validation[,-18])
    errormat<-table(out.rf,validation$CAT)
    OA<-sum(diag(errormat))/sum(errormat)
    OA.ave<-OA.ave+OA
}
OA.ave<-OA.ave/(fold-1.0)
print("randomforest OA")
print(OA.ave)
1分くらいで終わりました。早い。精度は0.82。強い。確かにこれはランダムフォレスト使いたくなります。
{{br}}エラーマトリクスは以下。
out.rf  1  2  3  4  5  6  7  8  9 10
    1  24  3  0  0  1  2  0  0  0  0
     2   0 37  0  0  1  0  0  0  0  0
     3   1  0  9  0  0  0  0  0  0  0
     4   0  0  0  4  0  0  0  0  0  0
     5   1  2  0  0 68  6  2  9  3  0
     6   0  0  0  1  3 21  1  0  1  0
     7   0  0  0  0  0  0  3  0  0  0
     8   1  0  0  0 11  2  0 85  1  0
     9   0  0  0  0  0  1  0  1 22  0
     10  0  0  0  0  0  1  0  0  0  8
やはりcropland(5)とrice(8)の見分けがしんどそうです。しかし、果樹園(7)が割と見分けられるようになりました。
{{br}}ちなみに、
> importance(mldata.scale.rf)
         MeanDecreaseGini
B1              129.03139
B2              150.97811
B3              193.81596
B4              164.04862
B5              215.51119
B6              210.49034
B7              204.33100
DEM             280.63575
NDVI            303.59199
NDWI            276.06097
SLOPE           200.76071
ASM              38.54532
CONTRAST         90.16030
ENTROPY          35.21386
VARIANCE        106.49409
NDVI1           207.83371
NDWI1           159.61049
で変数重要度が出せます(ジニ係数の平均低下量)。予想通りと言いますか、DEMやSLOPE、NDVI、NDWI、二時期の差分あたりが効いています。テクスチャの中でも重要度に開きがあるのも面白いですが、テクスチャに関しては算出時のウィンドウサイズやGLCM計算方向などによっても変わってくるので、話半分に聞いておくことにします。

!サポートベクタマシン
ランダムフォレストと並んで人気のある分類器です。マージン最大化とカーネルトリックにより高速かつ非線形な分類を実現します。分類だけなので、事後確率や変数重要度などの副産物は出ません。チューニングの依存度がだいぶ高いという印象が個人的にはありますので、[[このあたりを参考に|http://d.hatena.ne.jp/hoxo_m/20110325/p1]]グリッドサーチでしっかりチューニングしてみます。
library(e1071)
t<-tune.svm(CAT~.,data=tuning,gamma=10^(-5:5),cost=10^(-2:2),tunecontrol=tune.control(sampling="cross",cross=9)) #交差検定数はスタージェスの公式より決定
これがすごく時間かかりました(1〜数時間)。summary(t)とすると、gamma=0.1,cost=10がベストと出ました。gammaはガウスカーネルの指数部の係数、costはソフトマージンのコスト係数です。
OA.ave<-0
for(i in 2:fold){
    training<-mldata.scale[-c(group[[1]],group[[i]]),]
    validation<-mldata.scale[group[[i]],]
    mldata.scale.svm<-svm(CAT~.,data=training,gamma=0.1,cost=10) #4秒くらい
    out.svm <- predict(mldata.scale.svm,validation)
    errormat<-table(out.svm,validation$CAT)
    OA<-sum(diag(errormat))/sum(errormat)
    OA.ave<-OA.ave+OA
}
OA.ave<-OA.ave/(fold-1.0)
print("supportvector OA")
print(OA.ave)
精度は0.80でした。ランダムフォレストには劣りますがそれなりかと思います。エラーマトリクスを貼っておきます。
out.svm  1  2  3  4  5  6  7  8  9 10
      1  24  3  0  0  1  1  0  0  0  0
      2   1 36  1  0  0  0  0  0  0  0
      3   1  1  8  0  0  0  2  0  0  0
      4   0  0  0  4  0  0  0  0  0  0
      5   1  1  0  0 62  3  2  8  1  1
      6   0  1  0  0  4 25  1  0  2  0
      7   0  0  0  0  0  0  1  0  0  0
      8   0  0  0  1 16  3  0 85  5  0
      9   0  0  0  0  0  1  0  2 19  0
      10  0  0  0  0  1  0  0  0  0  7
やはりcropland(5)とrice(8),grassland(6)とcropland(5)が問題です。ところでSVMの場合はカーネルを変更する(多項式カーネルとかシグモイドカーネルとか)ことでも結果がかなり変わるはずですが、マニアック過ぎるのでスルーします。

!マッピング
精度比較で結果の良かったランダムフォレストとSVMで実際にマッピングしてみましょう。
{{br}}ちなみに、本来は交差検証のループごとにマップを作成して多数決をとる…みたいにするべきなのかもしれませんが、今回は最後のループで学習したマシンで分類を行いたいと思います。まずは特徴ベクトルとなるラスタマップを読み込みます。spgrass6ライブラリを使いますので、各自でインストールしておいてください(環境によって結構違うようなので、インストール部分は割愛)。
library(spgrass6)
name<-c("LC81070352015090_sr_band1_C","LC81070352015090_sr_band2_C","LC81070352015090_sr_band3_C","LC81070352015090_sr_band4_C","LC81070352015090_sr_band5_C","LC81070352015090_sr_band6_C","LC81070352015090_sr_band7_C","DSM_mosaic","LC81070352015090_NDVI","LC81070352015090_NDWI","DSM_mosaic_slope","LC81070352015090_NDVI_asm_ASM_0","LC81070352015090_NDVI_contrast_Contr_0","LC81070352015090_NDVI_entropy_Entr_0","LC81070352015090_NDVI_variance_Var_0","NDVI_diff","NDWI_diff") #名前ベクトルを用意
system("g.region -p") #総ピクセル数を確認。2776608ピクセルだそうです
feature<-data.frame(matrix(0,2715904,17)) #格納用のデータフレームを用意。この位の大きさなら楽勝
names(feature)<-c("B1","B2","B3","B4","B5","B6","B7","DEM","NDVI","NDWI","SLOPE","ASM","CONTRAST","ENTROPY","VARIANCE","NDVI1","NDWI1")
for(i in 1:length(name)){
     feature[,i]<-readRAST6(name[i])[[1]]
}
feature.scale<-scale(feature)
では、ランダムフォレストで分類マップを作ってみましょう。
lulc.rf<-predict(mldata.scale.rf,feature.scale) #3分くらいで終わりました
system("r.mapcalc \"box=null()\"")
box<-readRAST6("box") #出力用の配列を用意
box[[1]]<-as.numeric(lulc.rf) #numeric型に変更して格納
writeRAST6(box,"lulc_rf",overwrite=T) #書き出し
続いてSVMもやってみます。
lulc.svm <- predict(mldata.scale.svm,feature.scale) #8分くらいかかりました
box[[1]]<-as.numeric(lulc.svm)
writeRAST6(box,"lulc_svm",overwrite=T) #書き出し
save.image("ml.RData")
quit() #Rから抜ける
ではいよいよ可視化してみます。
r.colors map=lucl_rf color=rules rule=lulc_color.txt
d.rast lulc_rf
d.legend map=lulc_rf use=1,2,3,4,5,6,7,8,9,10
r.colors map=lulc_svm color=rules rule=lulc_color.txt
d.rast lulc_svm
d.legend map=lulc_svm use=1,2,3,4,5,6,7,8,9,10
*ランダムフォレストのマップ(みえない場合は[[こちら|https://sites.google.com/site/mizuochipublic/_/rsrc/1486745211947/yan-jiu-nei-rong-research-themes/landsat_rf.png?height=199&width=200]],[[カラーテーブルはこちら|https://sites.google.com/site/mizuochipublic/_/rsrc/1486745311698/yan-jiu-nei-rong-research-themes/saclass_legend.png?height=200&width=37]])
[[ランダムフォレストのマップ|http://pen.agbi.tsukuba.ac.jp/~ryuiki_plus/sample_data/Landsat/LULC/landsat_rf.png]]|http://pen.envr.tsukuba.ac.jp/~ryuiki_plus/sample_data/Landsat/LULC/landsat_rf.png]]
*SVMのマップ(みえない場合は[[こちら|https://sites.google.com/site/mizuochipublic/_/rsrc/1486745198186/yan-jiu-nei-rong-research-themes/landsat_svm.png?height=199&width=200]])
[[SVMのマップ|http://pen.agbi.tsukuba.ac.jp/~ryuiki_plus/sample_data/Landsat/LULC/landsat_svm.png]]|http://pen.envr.tsukuba.ac.jp/~ryuiki_plus/sample_data/Landsat/LULC/landsat_svm.png]]
{{br}}上記は筆者の作ったマップです。人によって結果は多少異なるはずですが、大雑把な雰囲気は一致していることと思います。
{{br}}番号とクラスの対応をもう一度記しておきます。
DF/1 /EF/2 bamboo/3 barren/4 croplands/5 grasslands/6 orchard/7 rice/8 urban/9 water/10
ちなみに、使ったカラーテーブル(lulc_color.txt)の中身は次のようなものです:
1 128:128:0 #DF,olive
2 0:128:0 #EF,green
3 0:0:0 #bamboo,black
4 255:255:255 #barren,white
5 255:128:0 #croplands,orange
6 128:255:0 #grasslands,lawngreen
7 255:0:0 #orchaard,red
8 255:255:0 #rice,yellow
9 128:128:128 #urban,grey
10 0:0:255 #water,blue
マップをみてみると、ランダムフォレストとSVMでだいぶ趣きが違います。ざっくりした見た目はランダムフォレストの方が綺麗ですね。SVMでは霞ヶ浦の中にcroplandがあったりして、全体的にノイジーです。また、ランダムフォレストでは平地の多くが水田に分類されていますが、SVMではcroplandに分類されています。どちらが正しいのかは、google-earth等と対比しながら細かく見ていくなり、現地におもむくなり、つくば界隈の土地利用区分とかを調べればわかりそうですが、そこまではやりません(どなたか加筆願います)。
{{br}}少し細かく見てみましょう。霞ヶ浦の北西に伸びていく水域(高浜入)周辺にズームしました。
*ランダムフォレストでズーム(みえない場合は[[こちら|https://sites.google.com/site/mizuochipublic/yan-jiu-nei-rong-research-themes/materials/rf_zoom.png?attredirects=0]])
[[ランダムフォレストでズームしてみた|http://pen.agbi.tsukuba.ac.jp/~ryuiki_plus/sample_data/Landsat/LULC/rf_zoom.png]]|http://pen.envr.tsukuba.ac.jp/~ryuiki_plus/sample_data/Landsat/LULC/rf_zoom.png]]
*SVMでズーム(みえない場合は[[こちら|https://sites.google.com/site/mizuochipublic/yan-jiu-nei-rong-research-themes/materials/svm_zoom.png?attredirects=0]])
[[SVMでズームしてみた|http://pen.agbi.tsukuba.ac.jp/~ryuiki_plus/sample_data/Landsat/LULC/svm_zoom.png]]|http://pen.envr.tsukuba.ac.jp/~ryuiki_plus/sample_data/Landsat/LULC/svm_zoom.png]]
*対応するtrue color画像(みえない場合は[[こちら|https://sites.google.com/site/mizuochipublic/yan-jiu-nei-rong-research-themes/materials/landsat_zoom.png?attredirects=0]])
[[対応するtrue color画像|http://pen.agbi.tsukuba.ac.jp/~ryuiki_plus/sample_data/Landsat/LULC/landsat_zoom.png]]|http://pen.envr.tsukuba.ac.jp/~ryuiki_plus/sample_data/Landsat/LULC/landsat_zoom.png]]
{{br}}恋瀬川が石岡方面に伸びていくのですが、30m分解能ではさすがに途中で消えてしまい、grasslandと分類されています。また、いずれのマップでも、もう少しurbanが多くてもいいのでは?という感じです(ランダムフォレストではgrasslandと水田が、SVMではcroplandが多すぎる)。森林は割とよく当たっているようですが、常緑林か落葉林かまではもう少し分解能が高くないとはっきりしません。
{{br}}
{{br}}もっと精度を上げるためにはどうすべきか?
*検証情報をカテゴリ間でできるだけ均等にとる、空間的にランダムにとる、精度よくとる、たくさんとる
*使用する衛星データや特徴量をもっと増やす(LSTや夜間光やSAR、もっと多時期、もっと高分解能)
*ディープラーニングなど、アルゴリズムに凝ってみる
*職人芸(エキスパートシステムとも言う)と併用する
などなどの点を改良していけば、もう少しよくなっていくと思います。
*検証情報については、人海戦術と気合で増やすか、クラウドなどから自動で探してくる仕組みを考えてみるか、さもなくばカテゴリ間が均等になるように人工的に調整する方法もあるそうです(多いカテゴリはランダムに減らし、足りないカテゴリは人工的に作るらしい)。
*使用データは、もちろん独立なデータセットをたくさん使うほうがいいのですが、データ取得にかかるコストやマシンスペックとの兼ね合いです。次元を増やすと、必要な教師データサイズや計算負荷が指数的に膨張していきます(次元の呪い)。また当然、たくさんのデータを使えば、前処理にかかる手間も増えます。
*ディープラーニングについては、そのうち試してみたいと思います。
*ここでエキスパートシステムと言っているのは、「標高◯◯m以上にこんな植生があるはずがない」といった専門知から、「湖の真ん中に耕作地があるわけないよね?」といった常識にいたるまで、あらかじめ分類ルールを人間が定めておくというやり方です。そうやって職人芸的にルールを決めるのは汎用性にかける(そもそも全てのルールを論うことなど不可能)ということで第二次AIブームは下火になったのかもしれませんが、別に併用してはいけない決まりはないので、機械学習と職人芸両方使えばいいのではないかと思います。例えば機械学習を走らせた後、「隣接するピクセルが全て水域なら、当該ピクセルも水域とする」というルールを追加しておけば、SVMの霞ヶ浦croplandは低減できるはずです。
*また、今回は取り上げませんでしたが、そういう先見知を取りこみやすいという意味で、ベイズ的学習器も機会があれば試してみたいです。あとオブジェクトベース分類とかも近年流行りっぽいので興味のあるかたは試していただければ。
*Rのベクトル・行列演算がわかりやすくてなかなか手放せないのですが、pythonやC++などに書き換える企画も機会があれば。