樹實(shí)現(xiàn)遙感巖性識別:特征構(gòu)建與調(diào)參實(shí)戰(zhàn))
簡介針對遙感圖像巖性識別中人工特征依賴與參數(shù)調(diào)優(yōu)繁瑣的痛點(diǎn)采用極端隨機(jī)樹模型并集成布谷鳥、粒子群等優(yōu)化算法實(shí)現(xiàn)端到端自動調(diào)參所有工程代碼基于Python編寫是遙感地學(xué)與機(jī)器學(xué)習(xí)交叉方向的高分畢設(shè)工程。項(xiàng)目面向計(jì)算機(jī)、人工智能、遙感地信等相關(guān)專業(yè)在校生及企業(yè)算法人員尤其適合作為課程設(shè)計(jì)、畢業(yè)設(shè)計(jì)或算法對比實(shí)驗(yàn)的參照方案。包體共10個(gè)文件含8個(gè)Python源碼、1個(gè)訓(xùn)練好的RF模型pickle文件及1份Markdown說明文檔代碼模塊覆蓋數(shù)據(jù)預(yù)處理、TXT轉(zhuǎn)CSV、模型構(gòu)建訓(xùn)練、布谷鳥與PSO尋優(yōu)等完整流程壓縮包僅63KB方便快速部署與二次開發(fā)。目前已有159人學(xué)習(xí)下載資源描述顯示測試運(yùn)行穩(wěn)定、答辯評價(jià)高既可直接復(fù)現(xiàn)巖性識別流程也便于替換數(shù)據(jù)以拓展至其他遙感分類任務(wù)。整體輕量且模塊化適合逐步跟蹤調(diào)試。1. 用極端隨機(jī)樹做遙感巖性識別為什么這個(gè)方案值得你親手跑一遍地質(zhì)填圖、礦產(chǎn)勘查甚至工程選址里最耗時(shí)的環(huán)節(jié)之一就是把遙感影像上的色調(diào)、紋理差異翻譯成巖性界線。傳統(tǒng)目視解譯靠肉眼看波段組合一個(gè)圖幅解譯下來要兩三天而且不同人圈出來的界線經(jīng)常對不上。極端隨機(jī)樹模型把這件事變成“選特征、點(diǎn)樣本、跑模型、出圖”的自動化流水線每個(gè)像元被當(dāng)成一個(gè)樣本光譜特征加地形特征送進(jìn)樹模型輸出巖性類別。這篇文章給你完整落地路徑巖性識別為什么能用遙感數(shù)據(jù)做、極端隨機(jī)樹比隨機(jī)森林強(qiáng)在哪、布谷鳥和粒子群兩種優(yōu)化算法怎么給模型調(diào)參以及全程 Python 代碼和寫文檔說明時(shí)該突出的重點(diǎn)。適合手里有遙感影像、想用一份能復(fù)用代碼直接跑通的人群新手能跟著做熟手能避開幾個(gè)常見大坑。2. 巖性識別的光譜基礎(chǔ)與極端隨機(jī)樹模型這對組合為什么靠譜2.1 遙感巖性識別的物理基礎(chǔ)我們從影像里能看到什么巖石種類不同造巖礦物組合不同反射光譜就有差異。碳酸鹽巖在 2.3 μm 附近有明顯的吸收特征黏土礦物在 2.2 μm 附近有 Al-OH 吸收帶鐵氧化物在可見光紅波段到藍(lán)波段之間吸收強(qiáng)烈。Landsat 8 OLI 的第二個(gè)到第七個(gè)波段恰好覆蓋這些區(qū)間所以可以用短波紅外兩個(gè)波段的比值識別黏土化用紅/藍(lán)比值識別鐵染再疊上 DEM 派生的坡度和坡向一個(gè)像元就能組成十幾維的特征向量。這就是整個(gè)巖性識別方案能成立的物理前提。要澄清一個(gè)定位遙感巖性識別和遙感圖像目標(biāo)檢測不是一回事。目標(biāo)檢測在影像里找“哪里有礦坑、哪里有水體”輸出外接框巖性識別做的是像元級分類每個(gè)像素被分到某個(gè)巖性類別輸出是一張柵格圖。既然是分類任務(wù)隨機(jī)森林、支持向量機(jī)也都能做真正的問題是它們在二三十維光譜地形特征、十幾萬甚至幾十萬像元樣本、類別又高度不平衡時(shí)精度和效率不好兼顧。這也是極端隨機(jī)樹在這類任務(wù)里被頻繁選用的原因。2.2 極端隨機(jī)樹 vs 隨機(jī)森林Extra-Trees 到底“極端”在哪里極端隨機(jī)樹Extremely Randomized Trees簡稱 Extra-Trees與隨機(jī)森林的差別從名字就能看出比隨機(jī)森林更隨機(jī)。隨機(jī)森林做兩件隨機(jī)的事——用 bootstrap 抽樣構(gòu)造每棵樹的訓(xùn)練子集在每個(gè)節(jié)點(diǎn)從隨機(jī)特征子集中尋找最優(yōu)分裂閾值。Extra-Trees 把隨機(jī)推到極致第一不抽樣每棵樹直接用全量訓(xùn)練集第二分裂閾值也是隨機(jī)生成的算法在幾個(gè)隨機(jī)閾值里挑一個(gè)讓不純度下降最大的而不是遍歷所有特征值找全局最優(yōu)。這兩處改動換來的是更低的方差和更快訓(xùn)練速度代價(jià)是單棵樹的偏差略微變大靠多棵樹集成把偏差補(bǔ)回來。在遙感巖性識別場景里這個(gè)特點(diǎn)非常合適。遙感像元樣本在空間上高度自相關(guān)同一巖性單元里相鄰像素光譜幾乎一樣bootstrap 抽樣對降低方差的幫助有限特征又多——原始波段、比值、地形、紋理疊起來能到二三十維其中很多特征互相相關(guān)隨機(jī)閾值分裂反而讓樹在特征選擇上更均勻不容易被一兩個(gè)強(qiáng)特征帶偏。我拿同一份數(shù)據(jù)對比過隨機(jī)森林和 Extra-Trees總體精度接近時(shí)極端隨機(jī)樹的訓(xùn)練時(shí)間大約少三分之一對特征噪聲也更不敏感。2.3 為什么必須參數(shù)調(diào)優(yōu)Extra-Trees 的門檻在超參數(shù)sklearn 里 ExtraTreesClassifier 的默認(rèn)參數(shù)是通用經(jīng)驗(yàn)值不是為遙感數(shù)據(jù)設(shè)計(jì)的。n_estimators 默認(rèn) 100對像元級分類動輒幾十萬樣本來說偏少max_features 默認(rèn) sqrt高維特征下每個(gè)節(jié)點(diǎn)只隨機(jī)考察四五個(gè)特征樹與樹之間相關(guān)性變高集成效果打折扣min_samples_leaf 默認(rèn) 1遙感數(shù)據(jù)里巖性邊界處有大量混合像元標(biāo)簽本身帶噪聲葉子上只有一個(gè)樣本很容易過擬合。這幾個(gè)參數(shù)的合理取值和數(shù)據(jù)集規(guī)模、類別數(shù)、特征維度強(qiáng)相關(guān)手調(diào)很容易陷入“調(diào)一個(gè)、壞一個(gè)”的循環(huán)。用網(wǎng)格搜索又太貴一個(gè) 5 折交叉驗(yàn)證的 Extra-Trees 在中等規(guī)模數(shù)據(jù)上要跑幾十秒三參數(shù)網(wǎng)格隨便一搜就是上千次實(shí)驗(yàn)還都是獨(dú)立實(shí)驗(yàn)完全浪費(fèi)計(jì)算量。布谷鳥搜索和粒子群優(yōu)化這類元啟發(fā)式算法反而更現(xiàn)實(shí)——它們不要求目標(biāo)函數(shù)可導(dǎo)把 Extra-Trees 的交叉驗(yàn)證精度當(dāng)黑匣子用種群迭代的方式在參數(shù)空間里找好點(diǎn)。接下來兩章先解決數(shù)據(jù)和樣本問題第四章給出完整調(diào)參實(shí)現(xiàn)。3. 遙感圖像標(biāo)注與特征棧構(gòu)建精度天花板其實(shí)在這一步先把結(jié)論放前面模型再強(qiáng)也救不了爛樣本和弱特征。巖性識別項(xiàng)目的精度上限早在你加載影像、圈樣本的那一刻就定死了。這一章把數(shù)據(jù)側(cè)的問題一次講透。3.1 數(shù)據(jù)源與特征層組合Landsat 8 DEM 常見搭配怎么搭我常用的特征層組合如下注意特征順序必須固定第五章會講為什么順序錯一位預(yù)測結(jié)果就面目全非。特征層來源/計(jì)算方式主要響應(yīng)的巖性信息OLI 波段 B2-B7Landsat 8 表面反射率鐵氧化物、黏土礦物的基本光譜輪廓鐵氧化物比值B4 / B2火成巖區(qū)氧化鐵含量差異黏土礦物比值B6 / B7黏土化、熱液蝕變帶識別坡度、坡向DEM 派生巖性抗風(fēng)化能力導(dǎo)致的地形差異GLCM 紋理對第一主成分做 3×3 窗口區(qū)分厚層塊狀砂巖與薄層泥巖NDVI 掩膜植被指數(shù)不直接入特征用于剔除植被干擾Landsat 8 OLI 的 B6、B7 兩個(gè)短波紅外波段是巖性識別的核心因?yàn)榇蠖鄶?shù)造巖礦物的診斷性吸收特征都在 1.6-2.4 μm 區(qū)間。B2-B5 提供鐵氧化物和植被的上下文信息。DEM 的作用是輔助尤其在沉積巖地區(qū)陡坎、緩坡與巖性差異有明顯的對應(yīng)關(guān)系坡度特征能把“光譜相似但地形表現(xiàn)不同”的兩類巖石分開。遙感圖像標(biāo)注的常見做法是打開影像套已有地質(zhì)圖把每個(gè)巖性單元手動圈幾個(gè)多邊形。這樣能得到大量樣本點(diǎn)但里面藏著一個(gè)大問題空間自相關(guān)。同一巖性單元內(nèi)部相鄰像素光譜高度相似它們并不是獨(dú)立樣本。如果把訓(xùn)練區(qū)和驗(yàn)證區(qū)從同一批多邊形里隨機(jī)切開驗(yàn)證集里全是訓(xùn)練樣本的“近親”交叉驗(yàn)證精度能報(bào)到 96% 以上把模型換到另一圖幅預(yù)測直接掉到 70% 出頭。解決思路是訓(xùn)練多邊形和驗(yàn)證多邊形在地理上物理隔開。比如某巖性單元在研究區(qū)有 5 個(gè)出露區(qū)塊用其中 3 個(gè)區(qū)塊訓(xùn)練剩下 2 個(gè)區(qū)塊完全不出現(xiàn)在訓(xùn)練集里只做驗(yàn)證。這個(gè)邏輯一定要寫進(jìn)文檔說明里答辯時(shí)老師最常問的“驗(yàn)證精度怎么來的”就在這里。3.3 樣本均衡與掩膜先把水體、植被和陰影剔掉巖性類別天然不平衡。大面積沉積巖類——砂巖、泥巖——樣本量輕松上萬小面積侵入巖脈可能只有幾百個(gè)像元。Extra-Trees 對不平衡的敏感度比 SVM 低但不處理照樣會出問題。我一般做兩件事先做掩膜把水體、濃密植被、云陰影對應(yīng)的像元全部置為無效因?yàn)檫@些地物的光譜和巖石無關(guān)留著只會讓模型學(xué)習(xí)“區(qū)分水和石頭”這類無效規(guī)則再做類別加權(quán)設(shè)置 class_weight 參數(shù)讓少數(shù)類在分裂時(shí)獲得更高的權(quán)重。數(shù)據(jù)加載和特征矩陣構(gòu)建的代碼長這樣import numpy as np import pandas as pd import rasterio from sklearn.preprocessing import LabelEncoder # 1. 讀入訓(xùn)練區(qū)樣本點(diǎn)坐標(biāo)與標(biāo)簽 train_pts pd.read_csv(train_samples.csv) # 列: x, y, lithology # lithology 是巖性字符串例如 limestone / granite / sandstone # 2. 按坐標(biāo)從多波段影像中提取光譜向量 def extract_pixels(src_path, pts): with rasterio.open(src_path) as src: rows, cols rasterio.transform.rowcol(src.transform, pts[x].values, pts[y].values) data src.read() # 形狀 (bands, H, W) samples data[:, rows, cols].T # 轉(zhuǎn)置為 (n_points, n_bands) return samples X_spectral extract_pixels(feature_stack.tif, train_pts) # 3. 拼接 DEM 派生的地形特征 X_dem train_pts[[slope, aspect]].values X np.hstack([X_spectral, X_dem]) # 4. 巖性字符串編碼為整數(shù)類別 le LabelEncoder() y le.fit_transform(train_pts[lithology]) print(特征矩陣形狀:, X.shape) # (樣本數(shù), 波段數(shù) 地形特征數(shù)) print(巖性類別:, le.classes_)邏輯說明第一步讀樣本點(diǎn) CSV坐標(biāo)必須和影像使用同一投影坐標(biāo)系否則提取出的光譜全是錯的。第二步用 rasterio.transform.rowcol 把投影坐標(biāo)換算成影像的行列號一次性取出所有波段對應(yīng)位置的像元值。這里用 numpy 數(shù)組索引代替 for 循環(huán)樣本量上萬時(shí)速度差異非常明顯。第三步把坡度、坡向橫向拼到光譜后面特征順序從這一步就固定下來后面訓(xùn)練和預(yù)測必須保持一致。第四步用 LabelEncoder 把巖性字符串變成 0、1、2 等整數(shù)。參數(shù)說明src.read() 不帶 window 會把整幅影像載入內(nèi)存訓(xùn)練區(qū)一般圖幅不大可以這么做如果影像有幾個(gè) GB必須用 windowed reading第六章給出分塊方案。rowcol 返回的是整數(shù)行列號numpy 的整數(shù)數(shù)組索引可以直接用于 data[:, rows, cols]。注意數(shù)據(jù)清洗要在特征提取之前完成NDVI 掩膜最好在生成特征棧時(shí)就把無效值替換為 NaN 或統(tǒng)一標(biāo)記不要留到建模階段再處理。4. Python 實(shí)現(xiàn)極端隨機(jī)樹訓(xùn)練 布谷鳥/粒子群調(diào)參全流程到這里數(shù)據(jù)已經(jīng)變成特征矩陣 X 和標(biāo)簽 y。這一章給出三塊可直接抄的代碼Extra-Trees 基線與交叉驗(yàn)證設(shè)置、布谷鳥搜索調(diào)參、粒子群調(diào)參。調(diào)參時(shí)注意只用訓(xùn)練區(qū)數(shù)據(jù)驗(yàn)證區(qū)碰都不能碰。4.1 極端隨機(jī)樹基線模型參數(shù)怎么設(shè)才不是默認(rèn)值先搭建基線模型目的是拿到一個(gè)“調(diào)參前的分?jǐn)?shù)”后面所有優(yōu)化結(jié)果都要跟它對比from sklearn.ensemble import ExtraTreesClassifier from sklearn.model_selection import StratifiedKFold, cross_val_score # 基線 Extra-Trees et_base ExtraTreesClassifier( n_estimators300, max_features0.5, # 每個(gè)節(jié)點(diǎn)隨機(jī)考察 50% 的特征 min_samples_leaf5, # 葉子節(jié)點(diǎn)至少 5 個(gè)樣本 min_samples_split10, bootstrapFalse, # Extra-Trees 標(biāo)準(zhǔn)做法不抽樣 n_jobs-1, random_state42 ) cv StratifiedKFold(n_splits5, shuffleTrue, random_state42) base_score cross_val_score(et_base, X_train, y_train, cvcv, scoringbalanced_accuracy).mean() print(f基線 balanced accuracy {base_score:.4f})邏輯說明這里用 balanced_accuracy 而不是普通 accuracy因?yàn)閹r性類別不平衡普通精度會被大面積類別綁架少數(shù)類分錯多少都看不出來。balanced accuracy 先逐類算召回率再取平均每一類權(quán)重相同。cross_val_score 內(nèi)部會自動訓(xùn)練和評估返回 5 折的分?jǐn)?shù)列表。參數(shù)說明n_estimators300 是考慮到像元級樣本量通常很大默認(rèn) 100 棵樹不夠穩(wěn)max_features0.5 是關(guān)鍵改動遙感特征里波段和比值高度相關(guān)sqrt 模式下每棵樹只看四五個(gè)特征容易漏掉波段比值之間的組合信息改成 0.5 讓每棵樹考察一半特征min_samples_leaf5 防止巖性邊界上的混合像元把樹逼到過擬合bootstrapFalse 是 Extra-Trees 區(qū)別于隨機(jī)森林的典型設(shè)置改成 True 反而改變算法性質(zhì)效果通常更差。4.2 布谷鳥搜索調(diào)參Lévy 飛行和 pa 丟棄率怎么實(shí)現(xiàn)布谷鳥搜索的核心是 Lévy 飛行新解不是朝固定方向走而是頻繁出現(xiàn)短步長加偶爾長跳躍的行走模式長跳躍負(fù)責(zé)跳出局部最優(yōu)。宿主鳥發(fā)現(xiàn)外來蛋后以概率 pa 丟棄部分解相當(dāng)于每輪淘汰一批差解保持種群多樣性。import numpy as np from math import gamma # 參數(shù)空間統(tǒng)一映射到 [0,1]解碼時(shí)還原成真實(shí)參數(shù) # 位置 pos[0] - n_estimators, pos[1] - max_features, pos[2] - min_samples_leaf def decode(pos): n_est int(pos[0] * (800 - 100) 100) # 100 ~ 800 max_feat pos[1] * 0.8 0.1 # 0.1 ~ 0.9 min_leaf int(pos[2] * 29 1) # 1 ~ 30 return n_est, max_feat, min_leaf def fitness(pos): n_est, max_feat, min_leaf decode(pos) model ExtraTreesClassifier( n_estimatorsn_est, max_featuresmax_feat, min_samples_leafmin_leaf, bootstrapFalse, n_jobs-1, random_state42 ) return cross_val_score(model, X_train, y_train, cv5, scoringbalanced_accuracy).mean() def levy_flight(beta1.5): sigma (gamma(1 beta) * np.sin(np.pi * beta / 2) / (gamma((1 beta) / 2) * beta * 2 ** ((beta - 1) / 2))) ** (1 / beta) u np.random.normal(0, sigma) v np.random.normal(0, 1) return u / (abs(v) ** (1 / beta)) def cuckoo_search(n_nests10, n_iter30, pa0.25): nests np.random.uniform(0, 1, (n_nests, 3)) best_score -np.inf best_pos None for it in range(n_iter): # 對所有巢做 Lévy 飛行更新 for i in range(n_nests): step levy_flight() other nests[np.random.choice(n_nests)] new_pos nests[i] 0.01 * step * (nests[i] - other) new_pos np.clip(new_pos, 0, 1) if fitness(new_pos) fitness(nests[i]): nests[i] new_pos # 宿主鳥發(fā)現(xiàn)布谷鳥蛋隨機(jī)丟棄一部分解 for i in range(n_nests): if np.random.rand() pa: nests[i] np.random.uniform(0, 1, 3) # 更新全局最優(yōu) for i in range(n_nests): s fitness(nests[i]) if s best_score: best_score s best_pos nests[i].copy() print(f迭代 {it1}: 最優(yōu) balanced accuracy {best_score:.4f}) return decode(best_pos), best_score best_params, best_score cuckoo_search() print(布谷鳥最優(yōu)參數(shù):, best_params, 得分:, best_score)邏輯說明每個(gè)巢就是一個(gè)參數(shù)組合decode 把 0-1 之間的連續(xù)值映射到參數(shù)真實(shí)區(qū)間這樣做是為了讓粒子/巢穴在同一個(gè)量綱的空間里移動避免 n_estimators 范圍幾百、min_samples_leaf 范圍只有幾十導(dǎo)致搜索效率失衡。fitness 函數(shù)每次調(diào)用都做一次 5 折交叉驗(yàn)證返回 balanced accuracy這是整個(gè)優(yōu)化過程的“黑匣子評估器”。參數(shù)說明步長縮放因子 0.01 控制 Lévy 跳躍的幅度太大容易在參數(shù)空間亂飛太小收斂慢pa0.25 是布谷鳥搜索的經(jīng)典取值表示每輪約四分之一的解被隨機(jī)重置這個(gè)值一般不用調(diào)。n_nests10、n_iter30 是兼顧計(jì)算量的常見設(shè)置注意一次 fitness 調(diào)用要訓(xùn)練 5 棵樹模型30 輪乘 10 個(gè)巢就是 300 次評估跑完全量數(shù)據(jù)可能要幾小時(shí)。實(shí)操建議先用下采樣到 2 萬樣本的數(shù)據(jù)把調(diào)參流程跑通拿到最優(yōu)參數(shù)區(qū)間后再放全量精修。4.3 粒子群優(yōu)化調(diào)參慣性權(quán)重衰減和速度更新粒子群走的是另一條路每個(gè)粒子記住自己的歷史最優(yōu)位置 pbest種群共享全局最優(yōu) gbest速度和位置按公式更新。我把慣性權(quán)重 w 從 0.9 線性衰減到 0.4前期大權(quán)重讓粒子廣泛探索后期小權(quán)重讓粒子在最優(yōu)附近精細(xì)搜索避免一上來就收斂到局部最優(yōu)。def pso_tune(n_particles12, n_iter30, w_start0.9, w_end0.4, c11.5, c21.5): dim 3 pos np.random.uniform(0, 1, (n_particles, dim)) vel np.random.uniform(-0.1, 0.1, (n_particles, dim)) pbest_pos pos.copy() pbest_score np.array([fitness(p) for p in pos]) gbest_idx np.argmax(pbest_score) gbest_score pbest_score[gbest_idx] gbest_pos pbest_pos[gbest_idx].copy() for it in range(n_iter): w w_start - (w_start - w_end) * it / n_iter # 慣性權(quán)重線性衰減 for i in range(n_particles): r1, r2 np.random.rand(2) # 速度更新慣性 個(gè)體認(rèn)知 群體認(rèn)知 vel[i] (w * vel[i] c1 * r1 * (pbest_pos[i] - pos[i]) c2 * r2 * (gbest_pos - pos[i])) # 位置更新并限制在 [0, 1] 邊界內(nèi) pos[i] np.clip(pos[i] vel[i], 0, 1) s fitness(pos[i]) if s pbest_score[i]: pbest_score[i] s pbest_pos[i] pos[i].copy() if s gbest_score: gbest_score s gbest_pos pos[i].copy() print(f迭代 {it1}: gbest {gbest_score:.4f}) return decode(gbest_pos), gbest_score pso_params, pso_score pso_tune() print(粒子群最優(yōu)參數(shù):, pso_params, 得分:, pso_score)邏輯說明速度更新三項(xiàng)分別代表三個(gè)方向的力——慣性項(xiàng)保留上一輪運(yùn)動趨勢認(rèn)知項(xiàng)把粒子拉向自己歷史最優(yōu)社會項(xiàng)把粒子拉向種群全局最優(yōu)。通過調(diào)節(jié) c1 和 c2 控制兩股力的強(qiáng)弱c1 太大會讓每個(gè)粒子只顧自己亂跑c2 太大會讓種群過早抱團(tuán)。位置更新后 np.clip 把粒子限制在 [0,1] 空間里解碼函數(shù)的區(qū)間邊界在那里等著。參數(shù)說明w 線性衰減是粒子群調(diào)參的經(jīng)典技巧從 0.9 到 0.4 是經(jīng)驗(yàn)區(qū)間c1c21.5 是比較保守的配置不會讓粒子飛太猛。粒子數(shù) n_particles12 三輪迭代的評估次數(shù)是 360 次和布谷鳥差不多。兩種算法跑完后對比分?jǐn)?shù)如果差異在 0.005 以內(nèi)取更簡單的那個(gè)參數(shù)組合不要迷信“看起來更高的那一位”。表布谷鳥搜索與粒子群調(diào)參設(shè)置對比項(xiàng)目布谷鳥搜索 CS粒子群 PSO種群大小10 巢12 粒子迭代輪數(shù)3030核心機(jī)制Lévy 飛行 pa 丟棄慣性權(quán)重 pbest/gbest 吸引關(guān)鍵參數(shù)步長 0.01pa0.25w 0.9→0.4c1c21.5典型適應(yīng)性跳出局部最優(yōu)能力強(qiáng)收斂速度快后期精細(xì)搜索4.4 調(diào)參結(jié)果回讀與模型保存優(yōu)化算法輸出的只是一組參數(shù)坐標(biāo)拿到后要重新訓(xùn)練一個(gè)完整模型并用空間隔離的驗(yàn)證區(qū)做最終評估from sklearn.metrics import balanced_accuracy_score import joblib best_model ExtraTreesClassifier( n_estimatorsbest_params[0], max_featuresbest_params[1], min_samples_leafbest_params[2], bootstrapFalse, n_jobs-1, random_state42 ) best_model.fit(X_train, y_train) train_score balanced_accuracy_score(y_train, best_model.predict(X_train)) val_score balanced_accuracy_score(y_val, best_model.predict(X_val)) print(f訓(xùn)練集 balanced acc {train_score:.4f}) print(f驗(yàn)證集 balanced acc {val_score:.4f}) joblib.dump(best_model, lithology_et_best.joblib) joblib.dump(le, label_encoder.joblib)邏輯說明train_score 和 val_score 的差距是判斷過擬合的直接證據(jù)。兩者相差超過 0.1說明參數(shù)搜出來的組合過分依賴訓(xùn)練區(qū)細(xì)節(jié)優(yōu)先檢查調(diào)參時(shí)是不是不小心把驗(yàn)證區(qū)樣本卷進(jìn) fitness 評估了。joblib 保存的模型文件里包含完整樹結(jié)構(gòu)后面逐像元預(yù)測直接加載不用重新訓(xùn)練。參數(shù)說明這里用的是隨機(jī)種子 42 固定復(fù)現(xiàn)如果你在其他論文里看到“以 0.87 為最優(yōu)參數(shù)”但無法復(fù)現(xiàn)大概率是沒固定 random_state。文檔說明里寫清隨機(jī)種子和交叉驗(yàn)證折數(shù)是評閱老師最看重的可復(fù)現(xiàn)性細(xì)節(jié)。5. 巖性識別避坑5 個(gè)讓分類結(jié)果翻車的常見問題這一章全部是踩過的坑按“現(xiàn)象 → 原因 → 解決”寫每條都能對應(yīng)到實(shí)際運(yùn)行時(shí)的報(bào)錯或異常結(jié)果。5.1 特征順序錯位訓(xùn)練跑得很好整幅圖預(yù)測時(shí)全是噪聲現(xiàn)象模型在訓(xùn)練集和驗(yàn)證集上的精度都很正常但用整幅影像預(yù)測時(shí)輸出圖像完全錯亂花崗巖分布區(qū)變成了砂巖紋理像是被打亂的馬賽克。原因預(yù)測階段重建特征向量時(shí)特征順序和訓(xùn)練時(shí)不一致。訓(xùn)練時(shí)你可能用“B2-B7 坡度 坡向”的順序拼接預(yù)測時(shí)用了“坡度 坡向 B2-B7”或者漏了某個(gè)波段比值。Extra-Trees 對特征順序不敏感但特征矩陣的列順序必須和訓(xùn)練時(shí)完全一致樹節(jié)點(diǎn)上記錄的分裂特征索引是按列位置存的。解決把特征棧的構(gòu)建寫成一個(gè)函數(shù)訓(xùn)練和預(yù)測都調(diào)用同一個(gè)函數(shù)生成特征矩陣不要在兩處各寫一遍。我在代碼里習(xí)慣把 band_list、terrain_features 定義成全局列表預(yù)測前打印一下 X.shape[1]和訓(xùn)練時(shí)對不上就直接停下排查。5.2 隨機(jī)切分驗(yàn)證集精度虛高到 97%一換圖幅就崩現(xiàn)象交叉驗(yàn)證精度 96% 以上模型在訓(xùn)練區(qū)附近看起來完美一旦預(yù)測到相鄰圖幅精度掉到 70% 左右。原因訓(xùn)練和驗(yàn)證樣本來自同一批巖性出露區(qū)空間上強(qiáng)相關(guān)的像素被隨機(jī)分到兩側(cè)驗(yàn)證集沒有獨(dú)立代表性。遙感像元不是獨(dú)立樣本同一巖性單元的相鄰像素光譜幾乎一樣模型記住的是局部光譜模式而不是巖性規(guī)律。解決按地理區(qū)塊劃分訓(xùn)練和驗(yàn)證同一巖性單元的不同出露區(qū)塊必須嚴(yán)格隔離。我見過的最簡單做法是在 GIS 里按多邊形編號切分區(qū)塊 id 為奇數(shù)的進(jìn)訓(xùn)練偶數(shù)的進(jìn)驗(yàn)證。這樣驗(yàn)證精度才有實(shí)際意義也禁得住答辯追問。5.3 調(diào)參時(shí)把驗(yàn)證區(qū)卷進(jìn)來布谷鳥“優(yōu)化”出的參數(shù)是假高分現(xiàn)象布谷鳥搜索報(bào)出的最優(yōu)分?jǐn)?shù) 0.93按最優(yōu)參數(shù)重新訓(xùn)練后驗(yàn)證區(qū)分?jǐn)?shù)只有 0.81差距大得離譜。原因fitness 函數(shù)里做交叉驗(yàn)證時(shí)用的是全量 X 而不是 X_train驗(yàn)證區(qū)樣本混進(jìn)調(diào)參過程優(yōu)化算法等于提前“看過答案”。這比隨機(jī)切分的坑更隱蔽因?yàn)楸砻嫔峡唇徊骝?yàn)證流程完整實(shí)際上數(shù)據(jù)泄漏已經(jīng)發(fā)生。解決調(diào)參前嚴(yán)格區(qū)分 X_train 和 X_valfitness 內(nèi)部只接收 X_train。一個(gè)檢查技巧把調(diào)參后的最優(yōu)參數(shù)用固定的隨機(jī)種子重跑一遍直接看驗(yàn)證區(qū)分?jǐn)?shù)如果明顯低于調(diào)參時(shí)報(bào)出的分?jǐn)?shù)立刻檢查數(shù)據(jù)切割順序。5.4 少數(shù)巖性類全滅花崗巖脈幾百個(gè)像元全被劃成砂巖現(xiàn)象混淆矩陣?yán)锎竺娣e砂巖類精度 90%花崗巖類召回率只有 8%模型幾乎沒把它識別出來。原因類別不平衡加上 Extra-Trees 的分裂準(zhǔn)則偏向樣本量大的類。雖然基線模型用了 balanced_accuracy 評估但算法本身的分裂過程仍然按基尼系數(shù)最小化少數(shù)類提供的分裂收益太小經(jīng)常被忽略。解決兩層處理。第一用 class_weightbalanced_subsample 給少數(shù)類加權(quán)sklearn 的 Extra-Trees 直接支持第二對少數(shù)類做樣本擴(kuò)增——不是合成少數(shù)類樣本而是對少數(shù)類多邊形內(nèi)的像元做輕度平移、旋轉(zhuǎn)生成更多訓(xùn)練樣本。擴(kuò)增幅度控制在 1.5-2 倍以內(nèi)過度擴(kuò)增會引入虛假空間模式。5.5 整幅影像預(yù)測內(nèi)存爆炸一次 read() 吃掉 16 GB 內(nèi)存現(xiàn)象訓(xùn)練和調(diào)參都順利到了整幅影像預(yù)測那一步程序直接卡死或報(bào) MemoryError。原因rasterio 的 src.read() 把整幅影像一次性載入內(nèi)存。以 7000×7000 像元、7 個(gè)波段的 GeoTIFF 為例float32 數(shù)據(jù)就是 1.3 GB 起步加上預(yù)測輸出數(shù)組、特征工程中間結(jié)果16 GB 內(nèi)存很快見底。解決分塊預(yù)測。用 rasterio 的 Window 按 512×512 像元讀取和寫出第六章給出完整分塊代碼。這個(gè)坑最容易在最后一步翻車提前做好分塊規(guī)劃輸出時(shí)再設(shè)定合理的壓縮參數(shù)可以完全避免。6. 逐像元預(yù)測出圖與 Kappa 驗(yàn)證把成果做成能上答辯臺的交付物6.1 分塊預(yù)測整幅影像rasterio 窗口讀取的正確姿勢import numpy as np import rasterio from rasterio.windows import Window def predict_map(model, src_path, out_path, block_size512): 按 512x512 窗口分塊預(yù)測整幅影像避免內(nèi)存溢出 with rasterio.open(src_path) as src: profile src.profile profile.update(dtypeuint8, count1, compressdeflate) height, width src.height, src.width with rasterio.open(out_path, w, **profile) as dst: for row in range(0, height, block_size): for col in range(0, width, block_size): win Window(col, row, min(block_size, width - col), min(block_size, height - row)) data src.read(windowwin) # (bands, h, w) n_bands, h, w data.shape # 每個(gè)像元展開成一個(gè)樣本特征順序與訓(xùn)練一致 flat data.reshape(n_bands, -1).T pred model.predict(flat).astype(uint8) dst.write(pred.reshape(h, w), 1, windowwin)邏輯說明外層循環(huán)按行、列步長 512 切窗口邊界處用 min 限制窗口尺寸防止索引越界。每個(gè)窗口讀取后先 reshape 成 (h×w, bands) 的樣本矩陣預(yù)測后再 reshape 回二維柵格寫入輸出文件。模型是之前用 joblib 加載的特征矩陣的波段順序必須與訓(xùn)練一致這也是第三章多次強(qiáng)調(diào)順序的原因。參數(shù)說明block_size512 是內(nèi)存與 IO 的折中窗口越大 IO 次數(shù)越少但內(nèi)存峰值越高compressdeflate 能顯著減小輸出 GeoTIFF 體積巖性分類圖類別少壓縮比通常能達(dá)到 5:1 以上。如果影像有無效值預(yù)測前先做掩膜把無效像元直接賦 0不要送進(jìn)模型。6.2 用空間隔離的驗(yàn)證區(qū)計(jì)算混淆矩陣和 Kappa出圖只是交付物的骨架精度驗(yàn)證才是答辯時(shí)能拿出手的硬數(shù)據(jù)。對驗(yàn)證區(qū)樣本計(jì)算混淆矩陣和 Kappa 系數(shù)from sklearn.metrics import confusion_matrix, cohen_kappa_score, classification_report y_val_pred best_model.predict(X_val) cm confusion_matrix(y_val, y_val_pred) kappa cohen_kappa_score(y_val, y_val_pred) print(混淆矩陣:\n, cm) print(fKappa 系數(shù) {kappa:.4f}) print(classification_report(y_val, y_val_pred, target_namesle.classes_))邏輯說明混淆矩陣要按驗(yàn)證區(qū)的真實(shí)標(biāo)簽逐類看重點(diǎn)看對角線外的錯誤去向——花崗巖被誤判成什么、砂巖和泥巖之間是否互相混淆。Kappa 系數(shù)在類別不平衡時(shí)比總體精度可靠得多一般認(rèn)為 0.8 以上屬于高度一致。如果 Kappa 比 balanced accuracy 低很多說明某些類的錯誤分布很不均勻模型可能在“用巖性大類掩蓋小類錯誤”。參數(shù)說明classification_report 輸出的 per-class precision、recall、F1 是文檔說明里最該放的三列數(shù)據(jù)。答辯時(shí)老師問“哪兩類最容易混淆”指著混淆矩陣說“砂巖和泥巖光譜相似錯誤主要發(fā)生在二者之間”比任何泛泛的話都有說服力。6.3 出圖的三個(gè)細(xì)節(jié)固定類別顏色、圖例和空間參照最后一個(gè)技巧是出圖。巖性分類結(jié)果圖最容易犯的錯是類別顏色不固定——同一套數(shù)據(jù)畫兩張圖繪圖庫自動配色花崗巖第一次是紅色第二次變成綠色圖例直接錯位。正確做法是給每個(gè)類別固定一個(gè) RGB 顏色字典全程復(fù)用。圖例上標(biāo)注類別名和對應(yīng)地質(zhì)年代空間參照信息從源影像的 profile 里繼承保證輸出 GeoTIFF 能直接疊加到 GIS 里和地質(zhì)圖對位置。我交這類作業(yè)前有個(gè)習(xí)慣出完圖一定把結(jié)果和原始影像在 GIS 里疊著看一遍從目視解譯的角度抽查幾個(gè)巖性邊界是否符合地形和色調(diào)的直覺。這個(gè)習(xí)慣救過我很多次——有一次模型把線性構(gòu)造帶上的花崗巖全錯分成閃長巖就是因?yàn)橛?xùn)練樣本里閃長巖多邊形不干凈混入了接觸變質(zhì)帶上的混合像元光看精度指標(biāo)完全看不出問題。邊界不符合地質(zhì)常識時(shí)先別懷疑模型回去檢查訓(xùn)練樣本問題多半在那里。希望幫到你。本文還有配套的精品資源點(diǎn)擊獲取