回歸與克里金插值:空間統(tǒng)計(jì)下的空氣質(zhì)量預(yù)測實(shí)戰(zhàn))
簡介一份完整的大學(xué)課程論文內(nèi)容圍繞地理加權(quán)回歸GWR模型與克里金插值法在空氣質(zhì)量指數(shù)預(yù)測中的應(yīng)用適合地理信息科學(xué)、環(huán)境統(tǒng)計(jì)及相關(guān)專業(yè)的學(xué)生與研究者參考。包內(nèi)為單個(gè)docx文檔大小約1.08MB包含論文正文與附錄中的三份Python源代碼分別對應(yīng)GWR模型預(yù)測、克里格插值預(yù)測及地圖可視化便于對照學(xué)習(xí)與復(fù)現(xiàn)實(shí)驗(yàn)。資源重點(diǎn)闡述了GWR如何通過局部參數(shù)估計(jì)刻畫空間異質(zhì)性以及克里金插值如何基于變異函數(shù)生成連續(xù)污染分布圖同時(shí)給出了模型校準(zhǔn)、帶寬優(yōu)化、預(yù)測評估等關(guān)鍵步驟的完整論述。目前已有1258人學(xué)習(xí)下載。對于希望掌握空間統(tǒng)計(jì)分析方法或完成相關(guān)課程設(shè)計(jì)、論文寫作的讀者這份資源提供了從理論到編程實(shí)現(xiàn)的系統(tǒng)參考具有較強(qiáng)的實(shí)用價(jià)值。1. 地理加權(quán)回歸為什么同一個(gè)氣象變量在不同點(diǎn)位的影響方向可能完全相反同樣一份空氣監(jiān)測數(shù)據(jù)直接丟進(jìn)普通線性回歸氣溫每升高 1℃對空氣質(zhì)量指數(shù)的影響被壓縮成一個(gè)固定的全局系數(shù)換成地理加權(quán)回歸之后同一個(gè)氣溫變量的回歸系數(shù)在北部山區(qū)可能低到 -224在南部城區(qū)卻接近 0 甚至為正方向都能反過來。這種「系數(shù)跟著坐標(biāo)走」的現(xiàn)象就是空間非平穩(wěn)性。地理加權(quán)回歸模型的思路很直白不要用一條全局?jǐn)M合線概括整個(gè)城市而是讓每個(gè)監(jiān)測點(diǎn)根據(jù)自己的地理位置從鄰近樣本中估計(jì)出一套專屬的局部回歸系數(shù)。配合克里格插值方法補(bǔ)全未布設(shè)站點(diǎn)的預(yù)測值就構(gòu)成了一套完整的空氣污染空間分析閉環(huán)。把這兩套底層機(jī)制拆開看你才會明白為什么空間統(tǒng)計(jì)不能只靠換一個(gè)復(fù)雜模型來解決。2. GWR核心機(jī)制空間權(quán)重矩陣、核函數(shù)與帶寬搜索的數(shù)學(xué)邏輯2.1 從全局參數(shù)到局部參數(shù)OLS為什么在空間數(shù)據(jù)上失靈普通最小二乘回歸的隱含假設(shè)是所有樣本共享同一個(gè)總體回歸關(guān)系。對空氣質(zhì)量數(shù)據(jù)來說這個(gè)假設(shè)幾乎必然不成立城市北部山區(qū)植被覆蓋率高、污染物擴(kuò)散條件好氣溫每升高 1℃對 AQI 的邊際影響可能很小而南部城區(qū)建筑密度大、尾氣排放集中同樣的氣溫升幅會顯著推高污染物濃度。如果用 OLS 估計(jì)得到的只是這些局部效應(yīng)的加權(quán)平均既不能代表北部也不能代表南部。GWR 的一般形式可以寫為y_i β0(ui, vi) Σ βk(ui, vi) · x_{ik} ε_i其中 (ui, vi) 是第 i 個(gè)觀測點(diǎn)的坐標(biāo)β0 和 βk 都隨地理位置變化。系數(shù)估計(jì)采用的是加權(quán)最小二乘形式β?(ui, vi) (X? W(ui, vi) X)?1 X? W(ui, vi) y這里的 W 是一個(gè)對角權(quán)重矩陣對角線上的每個(gè)元素代表每個(gè)觀測點(diǎn)對回歸點(diǎn) i 的空間權(quán)重。權(quán)重矩陣怎么構(gòu)造決定了 GWR 是能夠抓住局部細(xì)節(jié)還是退化成一種帶裝飾的 OLS。有一點(diǎn)很多人初次接觸時(shí)會忽略GWR 估計(jì)的不是一個(gè)系數(shù)向量而是 n 個(gè)系數(shù)向量——每個(gè)觀測點(diǎn)各有一套。你的輸出不是一張系數(shù)表而是一組按坐標(biāo)組織的系數(shù)曲面。每個(gè)變量的系數(shù)在每個(gè)觀測位置都有一個(gè)估計(jì)值伴隨的還有局部標(biāo)準(zhǔn)差、局部 t 值這些都可以直接落到地圖上做可視化。2.2 空間核函數(shù)bi-square 權(quán)重如何實(shí)現(xiàn)距離衰減要構(gòu)造權(quán)重矩陣 W首要任務(wù)是定義一個(gè)空間核函數(shù)。核函數(shù)回答的問題是站在當(dāng)前回歸點(diǎn)上周圍哪些觀測點(diǎn)有資格參與計(jì)算以及它們各自的話語權(quán)有多大。論文選用的核函數(shù)是 bi-square 函數(shù)其表達(dá)形式是當(dāng) d_ij ≤ b 時(shí)w_ij [1 - (d_ij / b)2]2 當(dāng) d_ij b 時(shí)w_ij 0其中 d_ij 是第 i 個(gè)回歸點(diǎn)和第 j 個(gè)觀測點(diǎn)之間的歐氏距離b 是帶寬。bi-square 核的特點(diǎn)是「軟硬結(jié)合」在帶寬范圍內(nèi)權(quán)重隨距離以二次方下降但不會突然截?cái)嘣趲掃吔缣幥『媒禐榱阈纬梢环N平滑截?cái)?。相比同屬常見選項(xiàng)的高斯核bi-square 核的計(jì)算效率更高因?yàn)槌鰩挼挠^測點(diǎn)權(quán)重直接置零矩陣運(yùn)算只需關(guān)注帶寬內(nèi)的鄰近樣本。這里有一個(gè)經(jīng)常被忽略的操作細(xì)節(jié)d_ij 的計(jì)算必須基于投影坐標(biāo)不能直接用經(jīng)緯度的度數(shù)。經(jīng)度 1 度和緯度 1 度代表的地面距離不同直接用經(jīng)緯度歐氏距離會讓權(quán)重矩陣產(chǎn)生方向性偏差進(jìn)而污染全部系數(shù)估計(jì)。數(shù)據(jù)預(yù)處理階段必須統(tǒng)一投影坐標(biāo)系這一步?jīng)]做對后面所有結(jié)果都不可信。2.3 帶寬搜索與模型復(fù)雜度CV、AICc、帽子矩陣的三角關(guān)系帶寬 b 的選擇是 GWR 模型校準(zhǔn)的核心關(guān)節(jié)。帶寬過小每個(gè)回歸點(diǎn)只用很少的鄰近樣本擬合局部系數(shù)方差巨大模型只是在背誦訓(xùn)練數(shù)據(jù)帶寬過大權(quán)重矩陣趨近均勻分布模型退化為 OLS空間細(xì)節(jié)全部丟失。論文采用自適應(yīng)帶寬策略不固定一個(gè)距離值而是讓每個(gè)回歸點(diǎn)根據(jù)自身周圍的站點(diǎn)密度調(diào)整鄰域范圍再統(tǒng)一搜索最優(yōu)參數(shù)。適應(yīng)度評估使用交叉驗(yàn)證分?jǐn)?shù)CV(b_s) Σ [y_i - ?_{-i}(b_s)]2其中 ?_{-i} 是在校準(zhǔn)過程中剔除第 i 個(gè)觀測點(diǎn)后用剩余樣本對該點(diǎn)的預(yù)測值。CV 值越小說明留一預(yù)測能力越強(qiáng)。但 CV 有一個(gè)眾所周知的缺陷它只評價(jià)訓(xùn)練數(shù)據(jù)內(nèi)部的留一表現(xiàn)并不顯式懲罰模型復(fù)雜度。數(shù)據(jù)量足夠大時(shí)帶寬略微偏小帶來的 CV 下降會被誤認(rèn)為擬合更好實(shí)際上模型已經(jīng)開始過擬合。這時(shí)需要 AICc 出場。論文使用的 AICc 公式是AICc 2n ln(σ?) n ln(2π) n(n tr(S)) / (n - 2 - tr(S))其中 n 是觀測點(diǎn)總數(shù)σ? 是估計(jì)的標(biāo)準(zhǔn)誤差tr(S) 是帽子矩陣 S 的跡。帽子矩陣 S 是一個(gè) n×n 的投影矩陣第 i 行可以表示為S_i x_i (X? W_i X)?1 X? W_itr(S) 在 GWR 中被稱為有效參數(shù)個(gè)數(shù)ENP。普通回歸的有效參數(shù)個(gè)數(shù)就是自變量個(gè)數(shù)加 1但 GWR 中每個(gè)樣本點(diǎn)都有自己的局部權(quán)重模型實(shí)際使用的自由度遠(yuǎn)高于 p1。帽子矩陣的跡恰恰度量了這種實(shí)際自由度帶寬越小權(quán)重矩陣越尖銳每個(gè)回歸點(diǎn)獨(dú)立貢獻(xiàn)的約束越強(qiáng)ENP 就越大。AICc 通過 tr(S) 對復(fù)雜度進(jìn)行懲罰避免只優(yōu)化訓(xùn)練誤差導(dǎo)致的過擬合。帶寬搜索的實(shí)現(xiàn)策略在論文中寫得很清楚先指定空間核函數(shù)再利用 CV 或 AICc 值做最優(yōu)帶寬搜索搜索算法采用 Golden Section黃金分割方法用搜索到的最優(yōu)帶寬構(gòu)造權(quán)重矩陣最后求解局部回歸系數(shù)矩陣 β?_k(ui, vi)。每次評估一個(gè)候選帶寬都要重新計(jì)算一遍完整的帽子矩陣因此這個(gè)環(huán)節(jié)的計(jì)算量在樣本量大的時(shí)候會顯著上升。CV 和 AICc 的選擇也值得斟酌論文提到有時(shí)候?yàn)榱吮苊膺^擬合并提高運(yùn)算效率可以直接計(jì)算 AICc 值來優(yōu)化帶寬——AICc 尋優(yōu)得到的帶寬通常比 CV 尋優(yōu)的帶寬偏大泛化穩(wěn)定性更好。3. GWR校準(zhǔn)全流程從變量篩選、七步校準(zhǔn)到系數(shù)t值空間可視化3.1 數(shù)據(jù)預(yù)處理全零變量為什么會拆掉回歸矩陣論文的原始自變量有 7 個(gè)temp氣溫、hpa氣壓、wet濕度、speed風(fēng)速、dir風(fēng)向距離值、mm、height高度。數(shù)據(jù)預(yù)處理階段第一個(gè)動作就是把 mm 剔除因?yàn)?mm 變量在所有觀測點(diǎn)上的值都是 0。一個(gè)全零變量放在回歸設(shè)計(jì)矩陣?yán)锏葍r(jià)于給模型添加了一個(gè)常數(shù)倍的截距列。設(shè)計(jì)矩陣 X 的列之間產(chǎn)生完全共線性矩陣 X?WX 變成奇異矩陣無法求逆。即使某些數(shù)值計(jì)算庫用偽逆強(qiáng)行求解得到的系數(shù)也會因?yàn)楣簿€性而劇烈震蕩。因此在進(jìn)入 GWR 之前做方差檢查幾乎是一條鐵律任何方差為零或接近零的變量都應(yīng)當(dāng)剔除。在真實(shí)項(xiàng)目中某個(gè)監(jiān)測因子整列為 0 的原因可能是傳感器斷電、單位換算錯誤或者該因子在特定時(shí)期確實(shí)沒有被檢測到。無論是哪種情況報(bào)告里都應(yīng)該明確寫出剔除的變量和原因否則復(fù)現(xiàn)者拿著原始 7 變量直接跑代碼會在矩陣求逆時(shí)收到一堆 NaN 警告。3.2 七步校準(zhǔn)流程從核函數(shù)到加權(quán)最小二乘求解GWR 模型的校準(zhǔn)流程可以歸納為六個(gè)基本步驟論文中從 1 到 7 編號數(shù)據(jù)預(yù)處理把數(shù)據(jù)整理成模型可接受的形式必要時(shí)對某些屬性值做插值補(bǔ)全。指定空間核函數(shù)選擇 bi-square 核或高斯核。利用 CV 或 AICc 值通過 Golden Section 搜索尋找最優(yōu)空間帶寬 b_s*。用最優(yōu)帶寬構(gòu)造權(quán)重矩陣 W。利用加權(quán)最小二乘公式求解局部回歸系數(shù)矩陣 β?_k(ui, vi)。輸出結(jié)果并在模型結(jié)果的基礎(chǔ)上進(jìn)行預(yù)測分析或顯著性測試。第 5 步中系數(shù)矩陣的求解并不一定一次到位。論文提到使用迭代加權(quán)最小二乘法IRLS來求解原因在于 GWR 的權(quán)重矩陣既包含空間核權(quán)重還可能包含穩(wěn)健估計(jì)中的殘差權(quán)重。IRLS 的常規(guī)做法是先用普通最小二乘得到初始系數(shù)計(jì)算殘差用殘差調(diào)整下一輪權(quán)重矩陣中每個(gè)觀測點(diǎn)的權(quán)重對殘差大的點(diǎn)降權(quán)再重新擬合如此往復(fù)直至收斂。這相當(dāng)于給空間權(quán)重疊了一層抗差保護(hù)避免個(gè)別異常監(jiān)測點(diǎn)主導(dǎo)局部擬合。3.3 局部系數(shù)顯著性檢驗(yàn)把 t 值畫回地圖上模型校準(zhǔn)完以后每個(gè)觀測點(diǎn)都有一組局部系數(shù)估計(jì)。論文對每個(gè)變量逐一進(jìn)行 t 檢驗(yàn)顯著水平取 α0.05并把各監(jiān)測點(diǎn)的顯著性結(jié)果落到地圖上這一步用到的關(guān)鍵計(jì)算是局部方差。對第 i 個(gè)回歸點(diǎn)的模型系數(shù)估計(jì)方差由下面兩個(gè)式子給出C_i (X? W_i X)?1 X? W_i var(β?_i) diag(C_i C_i?) · σ?2其中 σ?2 由估計(jì)的標(biāo)準(zhǔn)誤差平方而來diag 表示取矩陣對角元素。于是每個(gè)位置的局部 t 值為t_i β?_i / sqrt(var(β?_i))這個(gè) t 值可以和標(biāo)準(zhǔn)正態(tài)分布臨界值 1.96 比較α0.05 時(shí)所有 |t| 大于 1.96 的點(diǎn)位標(biāo)記為顯著。論文正是這樣把全部監(jiān)測點(diǎn)的 t 值分布畫到地圖上形成一張帶地理指向性的顯著性圖。以溫度為例在 α0.05 水平上溫度對 AQI 的影響在研究區(qū)域幾乎全域?yàn)樨?fù)但從南部到北部存在明顯的強(qiáng)度梯度南部的負(fù)效應(yīng)更強(qiáng)北部的負(fù)效應(yīng)減弱個(gè)別點(diǎn)位甚至接近 0。再看濕度正面影響集中在西北部和南部的幾個(gè)城區(qū)點(diǎn)位負(fù)面影響則以北部區(qū)域?yàn)樽畲竽骋槐辈奎c(diǎn)位擬合出的系數(shù)低至 -343是整個(gè)區(qū)域負(fù)向影響最極端的值。風(fēng)速的影響方向同樣存在分區(qū)東北部區(qū)域?yàn)檎绊懼鞒菂^(qū)及部分外圍區(qū)域?yàn)樨?fù)影響。高度變量的正面影響集中在中南部城區(qū)負(fù)面影響較大區(qū)域集中在北部及其周邊。這些空間異質(zhì)性細(xì)節(jié)是 OLS 的全局系數(shù)完全給不出的。每個(gè)系數(shù)都變成了一個(gè)可追問的空間變量——哪個(gè)區(qū)影響最強(qiáng)、哪個(gè)點(diǎn)位方向與全局相反。附錄里提供的地圖可視化程序代碼就是把這些系數(shù)和 t 值按地理坐標(biāo)著色成圖的關(guān)鍵環(huán)節(jié)。預(yù)測階段論文利用第 23 期數(shù)據(jù)進(jìn)行預(yù)測輸出 16 個(gè)監(jiān)測點(diǎn)的 AQI 預(yù)測值。預(yù)測結(jié)果與真實(shí)值對比均方誤差為 449.4411776。這里的預(yù)測邏輯值得注意GWR 預(yù)測依賴協(xié)變量的未來值——你需要先給出預(yù)測時(shí)刻的氣溫、氣壓、濕度等才能利用局部系數(shù)算 AQI。這在真實(shí)業(yè)務(wù)里意味著你需要一套數(shù)值天氣預(yù)報(bào)來驅(qū)動空氣質(zhì)量預(yù)報(bào)而這也正是 GWR 和克里金需要搭配使用的根本原因。4. 克里金插值變異函數(shù)、泛克里金與模型診斷指標(biāo)4.1 區(qū)域化變量與變異函數(shù)空間相關(guān)性的量化克里格插值方法是地統(tǒng)計(jì)學(xué)的核心內(nèi)容之一。它的出發(fā)點(diǎn)是一個(gè)假設(shè)未采樣點(diǎn)的值不是孤立隨機(jī)的而是與空間中已采樣點(diǎn)的值存在相關(guān)性——距離越近相關(guān)性越強(qiáng)。空氣質(zhì)量數(shù)據(jù)在某種程度上滿足這個(gè)假設(shè)污染物濃度在幾百米范圍內(nèi)往往是平滑變化的只有在強(qiáng)排放源附近才會出現(xiàn)突變。這種空間相關(guān)性通過半變異函數(shù)來量化。對于距離檔 h半變異函數(shù)定義為γ(h) (1 / 2N(h)) · Σ [Z(x_i) - Z(x_j)]2其中 N(h) 是距離落在該檔位內(nèi)的樣本點(diǎn)對數(shù)量。實(shí)際計(jì)算時(shí)先把樣本點(diǎn)兩兩配對按距離分檔在每個(gè)檔位內(nèi)求屬性差平方的平均值得到經(jīng)驗(yàn)變異函數(shù)再用理論模型指數(shù)、高斯、球狀等去擬合這條經(jīng)驗(yàn)曲線??死锝鸬暮诵墓绞前汛傈c(diǎn)的值寫成周邊觀測值的線性組合Z*(x0) Σ λ_i · Z(x_i)權(quán)重 λ_i 的求解有兩個(gè)約束。第一是無偏性要求所有權(quán)重之和等于 1這樣估計(jì)的期望才能等于真實(shí)期望。第二是最小方差要求在無偏約束下使插值方差最小。這是一個(gè)帶等式約束的優(yōu)化問題用拉格朗日乘子法可以解出一組線性方程方程的輸入是樣本點(diǎn)之間的變異函數(shù)值和樣本點(diǎn)到待估點(diǎn)的變異函數(shù)值輸出是每個(gè) λ_i??死锝鸩皇呛唵蔚姆淳嚯x加權(quán)而是把空間相關(guān)性用變異函數(shù)顯式建模后再通過最優(yōu)化計(jì)算權(quán)重這保證了估計(jì)是無偏且方差最小的。4.2 泛克里金當(dāng)區(qū)域化變量攜帶漂移趨勢普通克里金要求區(qū)域化變量滿足二階平穩(wěn)假設(shè)也就是均值恒定協(xié)方差只與距離有關(guān)。但空氣質(zhì)量數(shù)據(jù)往往不滿足這個(gè)條件城市熱島效應(yīng)、地形高度變化、人口密度分布都會讓 AQI 的均值在空間上出現(xiàn)明顯趨勢。如果直接用普通克里金預(yù)測趨勢會被吸收進(jìn)殘差導(dǎo)致變異函數(shù)在近距離處表現(xiàn)異常。泛克里金把問題拆成兩部分來處理Z(x) Σ f_l(x) · β_l R(x)其中 Σ f_l(x)β_l 是確定性漂移成分通常取低階多項(xiàng)式R(x) 是剩余項(xiàng)滿足二階平穩(wěn)假設(shè)。無偏條件因此發(fā)生變化權(quán)系數(shù) λ 在求和等于 1 之外還需要滿足一組與漂移基函數(shù)相關(guān)的約束方程。泛克里金用一句話概括把空間趨勢和隨機(jī)波動分開建模然后只對隨機(jī)波動的空間自相關(guān)求解權(quán)重。論文在空氣質(zhì)量案例中選擇的正是泛克里金。預(yù)處理階段同樣剔除了 mm使用 temp、hpa、wet、speed、dir、height 共 6 個(gè)自變量抽取 time22 期的樣本擬合模型選用高斯模型對半變異函數(shù)進(jìn)行擬合然后對 time23 期的坐標(biāo)點(diǎn)做預(yù)測。這里需要強(qiáng)調(diào)克里金與 GWR 的一個(gè)關(guān)鍵差異克里金的預(yù)測不需要提供未采樣點(diǎn)位置的協(xié)變量值它完全憑借空間相關(guān)性和已采樣數(shù)據(jù)插值。這在業(yè)務(wù)上的價(jià)值是可以直接把稀疏的監(jiān)測站點(diǎn)數(shù)據(jù)擴(kuò)展成連續(xù)污染濃度分布圖為整個(gè)行政區(qū)提供覆蓋式預(yù)測。4.3 高斯模型擬合與診斷指標(biāo)Q1、Q2、cR 怎么看論文選用的高斯模型其理論變異函數(shù)形式為γ(h) c0 c · [1 - exp(-h2 / a2)]擬合得到的三個(gè)參數(shù)分別是塊金常數(shù)為 3.88034880e01基臺值為 2.95014344e-02變程為 7.59106034e00。這三個(gè)參數(shù)的物理含義是塊金常數(shù)代表距離為零時(shí)仍然存在的測量誤差和微觀尺度變異基臺值代表變異函數(shù)趨于平穩(wěn)時(shí)的總方差變程表示空間相關(guān)性存在的尺度范圍超過這個(gè)距離變量之間基本不具備空間關(guān)聯(lián)。在這組數(shù)據(jù)里塊金常數(shù)明顯大于基臺值這是一個(gè)值得警覺的信號。它說明在現(xiàn)有監(jiān)測點(diǎn)間距尺度上樣本之間的空間自相關(guān)性很弱絕大多數(shù)方差來自測量誤差或更微觀尺度的過程。這也解釋了為什么論文給出的三項(xiàng)診斷指標(biāo)只能算「效果一般」。三項(xiàng)指標(biāo)中Q1 越接近 0 越好Q2 越接近 1 越好cR 越小越好。Q1 為 0.63 意味著預(yù)測存在系統(tǒng)性高估傾向Q2 為 1.126 說明模型對預(yù)測不確定性的估計(jì)略偏保守cR 為 6.68 說明理論模型與經(jīng)驗(yàn)變異函數(shù)的匹配度一般。三項(xiàng)指標(biāo)組合起來看高斯模型對這個(gè)數(shù)據(jù)集只能算「能用」談不上「優(yōu)秀」。預(yù)測階段利用 time23 期的坐標(biāo)點(diǎn)做預(yù)測16 個(gè)點(diǎn)的輸出值與真實(shí)值對比預(yù)測均方誤差為 423.2909487。值得對比的是這個(gè)數(shù)值比 GWR 的 449.44 還要略低一點(diǎn)。GWR 的預(yù)測誤差取決于協(xié)變量未來值是否準(zhǔn)確克里金則不依賴協(xié)變量只靠空間自相關(guān)。在樣本量小、協(xié)變量不確定性高的場景下克里金的預(yù)測反而更穩(wěn)。5. 避坑指南GWR與克里金實(shí)戰(zhàn)中的四個(gè)典型踩坑點(diǎn)5.1 坐標(biāo)系不統(tǒng)一經(jīng)緯度直接算距離導(dǎo)致系數(shù)方向反轉(zhuǎn)現(xiàn)象模型跑通了但回歸系數(shù)的符號和常識完全相反。比如溫度系數(shù)在某區(qū)域顯示為強(qiáng)正影響而氣象常識告訴你溫度升高應(yīng)該促進(jìn)擴(kuò)散、降低 AQI。原因在計(jì)算權(quán)重矩陣 W 時(shí)用了 WGS84 經(jīng)緯度坐標(biāo)的度數(shù)距離。經(jīng)緯度不是等距坐標(biāo)系——經(jīng)度 1 度在不同緯度對應(yīng)的地面距離差別很大直接用歐氏距離等于給不同方向分配了錯誤的權(quán)重最終系數(shù)方向和大小全部失真。解決所有距離計(jì)算必須在投影坐標(biāo)系下完成。城市尺度研究建議用 UTM 或高斯-克呂格投影單位是米。在代碼里加一步坐標(biāo)轉(zhuǎn)換進(jìn)入核函數(shù)之前統(tǒng)一檢查一遍各監(jiān)測點(diǎn)兩兩距離是否符合直覺。這個(gè)小習(xí)慣能擋掉后面一大半的臟數(shù)據(jù)問題。5.2 帶寬搜索陷入局部極小別把 Golden Section 當(dāng)萬能藥現(xiàn)象CV 曲線在帶寬軸上出現(xiàn)多個(gè)波谷Golden Section 搜索最終停在其中一個(gè)波谷上導(dǎo)致回歸系數(shù)估計(jì)出現(xiàn)斷崖式突變。原因Golden Section 是單峰函數(shù)搜索算法假設(shè)目標(biāo)函數(shù)只有一個(gè)極小值。當(dāng)監(jiān)測站點(diǎn)分布極不均勻時(shí)——中心城區(qū)密集、外圍稀疏——CV 函數(shù)可能出現(xiàn)多峰形態(tài)搜索結(jié)果的正確性完全取決于初始區(qū)間的位置。解決搜索之前先做一次網(wǎng)格粗掃。從帶寬范圍的最小值到最大值均勻取 20 個(gè)點(diǎn)逐個(gè)計(jì)算 AICc 或 CV把曲線畫出來。單峰再交給 Golden Section 細(xì)化多峰就換全局優(yōu)化策略。這件事花不了十分鐘卻能避免系數(shù)震蕩帶來的無窮麻煩。把優(yōu)化器當(dāng)黑匣子直接用是很多空間統(tǒng)計(jì)項(xiàng)目翻車的起點(diǎn)。5.3 變異函數(shù)塊金常數(shù)過大先分解趨勢再談插值現(xiàn)象克里金的變異函數(shù)擬合效果不好塊金常數(shù)達(dá)到 3.88e01遠(yuǎn)超基臺值 2.95e-02Q1 和 cR 明顯偏離理想值。原因原始數(shù)據(jù)里存在空間趨勢比如 AQI 從北部山區(qū)到南部城區(qū)整體遞增。趨勢沒有被分離就進(jìn)入變異函數(shù)擬合會被當(dāng)成「空間隨機(jī)波動」的一部分導(dǎo)致變異函數(shù)在距離接近零時(shí)就表現(xiàn)出很大的方差。解決先用最小二乘擬合一階或二階趨勢面把趨勢項(xiàng)和殘差項(xiàng)分離再對殘差項(xiàng)擬合變異函數(shù)。如果分離趨勢后塊金常數(shù)仍然偏大再考慮更換變異函數(shù)模型比如 Matérn 或指數(shù)型而不是一上來就斷定數(shù)據(jù)沒有空間相關(guān)性。5.4 只看訓(xùn)練 R 平方高興太早預(yù)測 MSE 才是硬指標(biāo)現(xiàn)象GWR 的訓(xùn)練決定系數(shù)從 OLS 的 0.417 提升到 0.916你覺得模型已經(jīng)接近完美結(jié)果用 time23 期數(shù)據(jù)一預(yù)測MSE 高達(dá) 449.44誤差并不小。原因GWR 的有效參數(shù)數(shù)量遠(yuǎn)高于 OLS。帶寬選擇偏小時(shí)局部系數(shù)會產(chǎn)生嚴(yán)重的逐點(diǎn)擬合傾向訓(xùn)練 R 平方高不代表泛化能力好換到新時(shí)刻新位置空間關(guān)系是否仍然成立完全是另一回事。解決嚴(yán)格做時(shí)間切片驗(yàn)證。用 time22 的數(shù)據(jù)擬合用 time23 的數(shù)據(jù)評估以預(yù)測 MSE 為最終對比指標(biāo)。論文里的 GWR 和克里金都做了這一步這是最值得借鑒的做法。任何人給一個(gè)孤零零的 R 平方你都應(yīng)該追問是訓(xùn)練集還是測試集。提示GWR 和克里金的結(jié)果輸出都帶空間坐標(biāo)預(yù)測驗(yàn)證時(shí)務(wù)必保證訓(xùn)練集和驗(yàn)證集的坐標(biāo)系完全一致否則誤差分析本身也會出偏差。6. 設(shè)計(jì)取舍解釋性模型與預(yù)測性模型如何互補(bǔ)6.1 GWR 的強(qiáng)項(xiàng)是解釋克里金的強(qiáng)項(xiàng)是預(yù)測論文同時(shí)采用 GWR 和克里金并不是重復(fù)建設(shè)而是兩種方法回答的問題不同。GWR 的空間非平穩(wěn)系數(shù)告訴你氣溫、氣壓、濕度這些因素對 AQI 的影響在哪個(gè)區(qū)域強(qiáng)、在哪個(gè)區(qū)域弱方向是否隨位置反轉(zhuǎn)??死锝饎t在采樣點(diǎn)之間填充未知位置把離散監(jiān)測點(diǎn)擴(kuò)展成連續(xù)的污染濃度面。從結(jié)果看克里金的預(yù)測 MSE423.29略低于 GWR449.44這不能說明克里金永遠(yuǎn)更好。背后的邏輯是GWR 做預(yù)測需要外部協(xié)變量的未來值協(xié)變量本身的預(yù)報(bào)誤差會傳導(dǎo)到 AQI 預(yù)測中克里金省略了協(xié)變量只借助空間自相關(guān)在協(xié)變量不可靠或缺失的場景里反而更穩(wěn)。但克里金回答不了「為什么北部濕度系數(shù)是 -343」這類機(jī)理問題這個(gè)解釋只能由 GWR 給出。兩種方法對這個(gè)案例而言是互補(bǔ)關(guān)系先用 GWR 建模出各影響因子的空間異質(zhì)性看清驅(qū)動機(jī)制再用克里金生成未采樣點(diǎn)的連續(xù)預(yù)測面服務(wù)實(shí)際決策。如果你只跑 GWR你會在沒有監(jiān)測站覆蓋的區(qū)域無所適從如果只跑克里金你對污染成因一無所知。從這份課程設(shè)計(jì)論文里我學(xué)到的最大教訓(xùn)是空間統(tǒng)計(jì)必須把「解釋」和「預(yù)測」分開思考。從那以后我每次拿到帶經(jīng)緯度坐標(biāo)的回歸任務(wù)都會強(qiáng)制先做三件事——第一畫一張因變量的空間分布圖確認(rèn)是否存在明顯的空間分層或梯度第二跑一遍 OLS 記錄 AICc作為 GWR 的對照基線第三先剔除方差接近零的冗余變量再統(tǒng)一投影坐標(biāo)系。三道工序做完GWR 和克里金的路才算真正鋪好。希望幫到你。本文還有配套的精品資源點(diǎn)擊獲取