簡介面向需要估算植被凈初級生產(chǎn)力NPP的遙感與生態(tài)研究人員朱文泉植被NPP程序V2.0.zip是一套基于ENVI平臺的NPP計算擴展工具。該工具對應(yīng)朱文泉博士2020年發(fā)布的2.0模型內(nèi)置CASA NPP 2.0模塊可直接利用日照時數(shù)估算太陽輻射讀取NDVI HDF格式生成時間序列NDVI并在同一模塊內(nèi)完成氣象數(shù)據(jù)插值減少跨軟件操作環(huán)節(jié)。壓縮包共76個文件大小約30.58MB核心程序為sav文件另有dat/hdr格式的柵格數(shù)據(jù)、csv/txt參數(shù)與說明、pdf參考文檔、shp等矢量邊界文件以及Install.bat/Unstall.bat安裝卸載腳本目錄結(jié)構(gòu)涵蓋ENVI Data、Documents、License、Tools等模塊其中ENVI Data存放示例數(shù)據(jù)Documents提供使用說明Tools為擴展工具目錄劃分清晰。已有5548人學習/下載。部署后可依據(jù)readme.txt和示例數(shù)據(jù)快速了解工具用法按模型要求準備輸入數(shù)據(jù)即可啟動運算無需額外編程加之Install.bat/Unstall.bat腳本可快速完成擴展安裝與卸載適合希望提升NPP估算效率的ENVI中高級用戶。 做植被遙感或者生態(tài)模型反演的人十有八九都繞不開朱文泉老師的NPP估算程序。這套基于光能利用率模型CASA模型的植被凈初級生產(chǎn)力估算工具在國產(chǎn)遙感軟件里屬于流傳最廣的那一類從碩士論文到國家自然科學基金項目到處能看到它的影子。這次拿到的是“朱文泉植被NPP程序V2.0.zip”這個版本相比早期流傳的V1.x版本V2.0在輸入數(shù)據(jù)組織、參數(shù)文件外置、批量處理邏輯上都做了不少調(diào)整實操時的坑點也跟老版本不太一樣。這篇文章我打算從一個實際跑過這個模型、也踩過一堆坑的從業(yè)者角度把這個包從解壓到出圖的全流程拆開講。包括CASA模型的原理怎么落到這套程序里、輸入數(shù)據(jù)到底要準備到什么程度、運行時的參數(shù)文件怎么改才不會報錯、輸出的結(jié)果又該怎么驗證合理性。手里有數(shù)據(jù)準備上手跑的可以直接按著步驟操作順便把我的排查經(jīng)驗也帶走。1. NPP估算的整體設(shè)計與算法思路拆解1.1 為什么選擇光能利用率模型先簡單交代一下背景NPPNet Primary Productivity凈初級生產(chǎn)力是指綠色植物在單位面積、單位時間內(nèi)通過光合作用積累的有機物質(zhì)總量扣除自身呼吸消耗后剩下的部分。它是碳循環(huán)研究的核心指標之一也是估算陸地生態(tài)系統(tǒng)碳匯能力的底層數(shù)據(jù)。要拿NPP目前無非三條路實測樣地數(shù)據(jù)、過程模型模擬、遙感參數(shù)反演。實測最準但覆蓋范圍極其有限過程模型如Biome-BGC、BEPS對參數(shù)要求極多而且跑起來費勁真正適合大范圍、長時序、快速估算的還得說遙感驅(qū)動的光能利用率模型。CASA模型在光能利用率模型里屬于祖師爺級別的存在最早由Potter等人于1993年提出基本邏輯就是把NPP拆成兩個因子的乘積植被吸收的光合有效輻射APAR以及實際光能利用率ε。朱文泉的這套程序在經(jīng)典CASA基礎(chǔ)上做了關(guān)鍵的本土化改進最大的亮點是針對不同植被類型重新率定了最大光能利用率參數(shù)ε_max這直接解決了國外模型參數(shù)在國內(nèi)植被區(qū)劃下“水土不服”的問題。V2.0版本相比之前的版本程序結(jié)構(gòu)上更模塊化參數(shù)和數(shù)據(jù)文件做成了外部配置換數(shù)據(jù)跑新區(qū)域時不需要再去翻源碼改常量。1.2 CASA模型的基本公式與參數(shù)邏輯整個模型的核心公式可以拆成下面這層關(guān)系NPP(x, t) APAR(x, t) × ε(x, t)其中APAR是植被吸收的光合有效輻射計算依賴兩個輸入太陽總輻射SOL和植被對光合有效輻射的吸收比例FPAR。而FPAR通常是由NDVI數(shù)據(jù)反演出來的這也是為什么這套程序一定要有逐月的NDVI數(shù)據(jù)作為驅(qū)動。另一頭ε(x, t) T_ε1 × T_ε2 × W_ε × ε_max。這里面T_ε1和T_ε2是溫度脅迫系數(shù)分別反映低溫抑制光合作用和溫度偏離最適溫度時的影響W_ε是水分脅迫系數(shù)由區(qū)域?qū)嶋H蒸散和潛在蒸散的比值決定ε_max就是最大光能利用率也是朱文泉模型本土化率定的核心參數(shù)。V2.0程序里不同類型植被的ε_max值直接做成了一個可讀寫的參數(shù)表文件換植被類型時不需要再改動代碼這個設(shè)計非常實用。理解了這套公式就明白程序為什么需要那么多輸入文件了——它本質(zhì)上就是逐月、逐像元地把NDVI數(shù)據(jù)、氣象數(shù)據(jù)溫度、降水、太陽輻射代入上述公式進行空間計算。運行結(jié)果就是逐月的NPP空間分布數(shù)據(jù)累加后得到年NPP。1.3 這版程序的核心框架與流程從運行邏輯上看V2.0可以分成數(shù)據(jù)導入層、參數(shù)配置層和計算輸出層。數(shù)據(jù)導入層負責讀取按月組織的NDVI、氣溫、降水、太陽輻射數(shù)據(jù)參數(shù)配置層通過一個參數(shù)文件控制模型的關(guān)鍵常數(shù)包括ε_max、各脅迫系數(shù)的計算方法、輸出的文件名前綴等計算輸出層則執(zhí)行逐月NPP估算并自動完成年值合成。這個結(jié)構(gòu)最大的好處是——換數(shù)據(jù)區(qū)域、換植被類型、換時間尺度時不用動程序主體只需修改參數(shù)文件。這對于做不同區(qū)域?qū)Ρ妊芯康挠脩魜碚f非常友好V1版本那種動不動需要改源碼的日子總算過去了。2. 運行前的數(shù)據(jù)準備與參數(shù)配置2.1 程序解壓與環(huán)境配置先說解壓這個聽起來簡單但我在實際使用中真見過不少人在這一步出錯。V2.0壓縮包下載完成后建議先放到一個路徑全英文、無空格的目錄下再解壓比如D:\NPP_Model或者E:\RS_Data\NPP_Model切忌解壓到桌面或者帶中文路徑的文件夾里。原因很簡單程序內(nèi)部調(diào)用文件時如果路徑帶中文在部分系統(tǒng)編碼環(huán)境下會直接讀取失敗而且報錯信息往往還是亂碼排查起來很費勁。程序主體是MATLAB環(huán)境下的.m腳本文件所以你需要一個可用的MATLAB環(huán)境。我實測過2016b到2023a之間的多個版本都能正常運行但建議至少用2014b以上版本否則部分日期函數(shù)和矩陣操作可能出現(xiàn)兼容問題。打開MATLAB后先用addpath(genpath(D:\NPP_Model))把整個程序目錄加入搜索路徑然后檢查一下能否正常打開主函數(shù)NPP_main_V2.m。2.2 輸入數(shù)據(jù)清單與格式要求這套程序的輸入數(shù)據(jù)按月份組織通常是一年12個月為一輪核心輸入包括四類逐月NDVI數(shù)據(jù)一般用SPOT/VGT或者MODIS的NDVI產(chǎn)品空間分辨率需要統(tǒng)一重采樣格式為GeoTIFF或ENVI標準柵格投影建議統(tǒng)一為Albers等積投影或UTM。逐月平均氣溫數(shù)據(jù)單位攝氏度柵格大小和范圍必須與NDVI數(shù)據(jù)完全一致。逐月降水量數(shù)據(jù)單位毫米同樣要求空間范圍與NDVI一致。逐月太陽總輻射數(shù)據(jù)單位MJ/m2如果缺少站點觀測輻射數(shù)據(jù)可以用日照時數(shù)通過Angstrom公式估算得到。數(shù)據(jù)層面最核心的一個坑就是所有柵格的像元行列數(shù)、投影、空間范圍必須完全對齊一個像元都不能偏。大多數(shù)運行報錯或者結(jié)果出現(xiàn)大片空值都源于這個問題。我在預(yù)處理的環(huán)節(jié)通常會用ArcGIS或ENVI把全部數(shù)據(jù)先做一次批量重采樣和裁剪統(tǒng)一到同一模板柵格上再喂給程序。2.3 參數(shù)配置文件怎么修改V2.0的一個核心改進是將模型參數(shù)集中到了一個配置文件里通常叫NPP_paras.txt或CASA_parameters.csv。打開這個文件重點需要修改的項包括植被類型代碼表程序需要將土地利用數(shù)據(jù)中的植被類型與ε_max取值對應(yīng)一般的做法是在配置文件中建立一個編碼映射表比如1代表常綠針葉林2代表落葉闊葉林以此類推。最大光能利用率ε_max朱文泉針對中國典型植被類型做了參數(shù)率定不同植被類型的取值有明顯差異。以我跑過的華北地區(qū)為例落葉闊葉林的ε_max明顯高于草地和灌叢這一項如果填錯了最終NPP結(jié)果會系統(tǒng)性偏高或偏低。輸出路徑與文件前綴指定結(jié)果輸出目錄以及月NPP和年NPP輸出文件命名的前綴避免覆蓋上一次運行的結(jié)果。修改配置時務(wù)必注意保存格式不要用Windows記事本另存為帶BOM的UTF-8編碼否則MATLAB讀取配置文件時第一行會多出一個不可見字符導致程序報錯找不到參數(shù)項。用Notepad或VS Code修改并保存為無BOM的UTF-8是更穩(wěn)妥的做法。3. 實操過程與NPP結(jié)果計算3.1 從數(shù)據(jù)導入到逐月NPP輸出拿到整理好的數(shù)據(jù)后我第一次跑通全流程大概花了兩個晚上主要是卡在數(shù)據(jù)格式統(tǒng)一上?，F(xiàn)在把完整步驟梳理一下打開主程序NPP_main_V2.m在代碼頭部區(qū)域能看到一個名為data_root的字符串變量把它改成存放輸入數(shù)據(jù)的根目錄路徑。注意程序會在這個目錄下自動查找固定命名的子文件夾比如NDVI、TEM、PRE、RAD每個子文件夾內(nèi)存放按月份命名的柵格文件命名規(guī)則通常是NDVI_2000_01.tif這樣的格式。如果命名不吻合程序會報“文件不存在”的錯誤所以花幾分鐘把數(shù)據(jù)文件命名統(tǒng)一到位是值得的。在MATLAB命令窗口直接運行主程序前可以先在腳本中單獨執(zhí)行一次數(shù)據(jù)讀取測試。V2.0版本考慮到了快速排查的問題通常提供了check_input或類似輔助函數(shù)運行后會在命令窗口打印每個輸入文件的大小、波段數(shù)、行列數(shù)這一招在確認數(shù)據(jù)對齊狀態(tài)時非常有用。正式運行后程序會逐月計算APAR和ε最終輸出各月的NPP柵格。以2000年12個月的數(shù)據(jù)為例運行時長取決于柵格大小和電腦性能1000×1000像元左右的數(shù)據(jù)量在普通臺式機上大概需要20到40分鐘期間MATLAB會顯示循環(huán)進度條。我第一次跑的時候嫌慢后來發(fā)現(xiàn)通過修改程序中的循環(huán)并行開關(guān)如果V2.0的代碼里含有parfor字樣在開啟MATLAB并行計算池后跑速能提升一倍以上。3.2 年NPP合成與結(jié)果輸出月尺度NPP數(shù)據(jù)只是中間成果絕大多數(shù)研究最終需要的是年NPP。V2.0程序通常自帶年值合成模塊也就是把12個月的逐月NPP柵格相累加。這里有一個容易引發(fā)爭議的點關(guān)于NPP累加時若某個月份因為云覆蓋導致NDVI異常低其對應(yīng)的NPP值可能為零或負值直接累加會造成年度總量偏低。我在實際項目中摸索出的處理方案是在合成之前先對月NPP做置信區(qū)間截斷將異常值用相鄰年份同月均值填補然后再累加這樣得到的年NPP空間分布要平滑合理得多。輸出結(jié)果還應(yīng)該包含貢獻量或精度評價文件。V2.0程序在輸出NPP柵格的同時通常還會生成一個運行日志如log.txt記錄每次運行的基本參數(shù)、輸入文件名和完成時間。遇到結(jié)果不合理的情況先翻日志大部分問題都能定位到是輸入數(shù)據(jù)的問題還是參數(shù)配置問題。3.3 計算結(jié)果的合理性驗證拿到結(jié)果后千萬不要急著出圖寫論文先做合理性驗證。驗證主要看三個方面值域范圍、空間分布和總量量級。從量級上看中國大部分陸地生態(tài)系統(tǒng)的年NPP在0到1500 gC/m2/yr之間全國平均大致在300到500 gC/m2/yr左右。如果你算出來的某區(qū)域年均NPP普遍超過1500或者大量出現(xiàn)負值數(shù)據(jù)或參數(shù)大概率有問題。空間分布上的常識性判斷也很重要。以我測試過的東北地區(qū)為例大興安嶺的落葉針葉林、小興安嶺和長白山的針闊混交林其NPP理應(yīng)明顯高于同緯度的農(nóng)田和草地。如果空間分布呈現(xiàn)出緯度梯度與常識相反重點檢查FPAR的計算和ε_max的分配是否出錯。另外把模擬結(jié)果與文獻實測值比對是一種非常高效的驗證手段。比如長白山闊葉紅松林的實測NPP大致在500到800 gC/m2/yr之間如果模擬值在這個范圍內(nèi)可信度就比較高。V2.0程序包里如果附帶了參考的NPP結(jié)果樣例可以拿自己的輸出與樣例做一次空間相關(guān)分析相關(guān)系數(shù)在0.7以上基本說明流程跑通了。4. 常見問題與排查技巧實錄4.1 解壓、路徑與MATLAB環(huán)境類問題結(jié)合網(wǎng)上討論比較多的情況我把這幾年跑這套程序?qū)嶋H遇到的、以及幫別人排查過的高頻問題整理成速查表方便你們直接對照解決。第一類問題是壓縮包和路徑相關(guān)的?！敖鈮汉筇崾疚募p壞”通常是因為下載不完整或者殺毒軟件攔截了部分文件重新下載并關(guān)閉實時防護后解壓即可。網(wǎng)上有提到過“failed to copy spatial iop zip”這類報錯多見于用第三方解壓工具時出現(xiàn)我的建議是優(yōu)先用WinRAR或7-Zip并選擇“解壓到當前文件夾”而不要直接拖拽文件。第二類問題是MATLAB報“無法將xx識別為命令”這本質(zhì)上就是路徑?jīng)]添加到位執(zhí)行addpath(genpath(程序目錄))后再次運行即可。還有一種非常隱蔽的情況就是你電腦上裝了多個MATLAB版本默認打開的是舊版本有些函數(shù)在新版本中已改動但舊版本中不存在導致腳本跑到一半就中斷這種情況下注意把MATLAB啟動時的“當前文件夾”切到程序目錄或者直接在腳本開頭寫入cd(D:\NPP_Model)。4.2 輸入數(shù)據(jù)與格式類排查這一類遇到的問題最多也是排查成本最高的。比如程序直接報“數(shù)據(jù)維度不匹配”問題基本鎖定在柵格行列數(shù)不一致上用ENVI打開兩個文件對比一下要么重新采樣要么做裁剪。還有一種情況是NDVI數(shù)據(jù)沒有乘以10000的縮放因子導致值域在-0.1到0.9之間的浮點數(shù)變成了-1000到9000的整型數(shù)程序?qū)PAR反演時直接算出了完全離譜的值。這時候去數(shù)據(jù)說明文檔里查一下原始NDVI產(chǎn)品的標定系數(shù)乘以或除以對應(yīng)系數(shù)即可。還有一種比較坑的情況就是氣溫數(shù)據(jù)中存在極端的背景值例如空值區(qū)域用-9999填充但程序讀取時未將其識別為無效值導致這些像元計算出的NPP是一個負數(shù)極端值直接污染了后續(xù)的年值統(tǒng)計。排查方法是在數(shù)據(jù)預(yù)處理階段將所有無效值統(tǒng)一設(shè)置為掩膜值如NaN并在程序中設(shè)置忽略NaN選項。4.3 運行報錯與參數(shù)設(shè)置類問題運行中途報錯最常見的有兩個。第一個是“索引超出矩陣維度”多發(fā)生在輸入數(shù)據(jù)的文件命名與程序預(yù)設(shè)不一致時程序內(nèi)部按月份序號循環(huán)讀文件某個月份缺數(shù)據(jù)或命名不規(guī)律就會觸發(fā)。解決辦法就是按預(yù)設(shè)規(guī)則逐月核對文件名。第二個是“內(nèi)存不足”當數(shù)據(jù)范圍較大或像元分辨率過高時程序需同時加載全年12個月的柵格到內(nèi)存中運算4G內(nèi)存的老機器很容易頂不住我的經(jīng)驗是把計算區(qū)域分塊處理或者適當降低數(shù)據(jù)分辨率以減少像元數(shù)。參數(shù)設(shè)置方面ε_max如果設(shè)置得過高NPP結(jié)果會普遍偏大反之亦然。這里建議查一下當前研究區(qū)內(nèi)主要植被類型的實測NPP文獻數(shù)據(jù)反推合理取值。另外溫度脅迫系數(shù)的計算需要輸入月平均氣溫和月最低氣溫如果只有月平均氣溫數(shù)據(jù)千萬別把最低氣溫那一欄留空否則T_ε1的計算會出問題。我的做法是利用氣溫日較差數(shù)據(jù)估算月最低氣溫雖然不精確但至少保證模型結(jié)構(gòu)完整。4.4 結(jié)果異常的處理經(jīng)驗結(jié)果出現(xiàn)大面積零值先檢查NDVI中是否有些月份數(shù)據(jù)缺失被填充成了00值的NDVI反演出的FPAR可能為0進而導致APAR為0。解決辦法是用時間序列插值法補全NDVI空缺月份常用的有Savitzky-Golay濾波法或時間窗口內(nèi)的線性插值法。結(jié)果出現(xiàn)“條紋狀”或“環(huán)形”異常時大概率是插值后的氣象數(shù)據(jù)存在邊界不連續(xù)問題。這是因為觀測站點數(shù)據(jù)經(jīng)過插值之后在空間上存在以站點為中心的“牛眼”效應(yīng)直接代入模型后NPP的空間分布也會有明顯的異常環(huán)。我的處理辦法是改用ANUSPLIN或者薄板樣條插值方法生成氣溫和降水柵格這樣得到的空間連續(xù)性會比普通IDW插值好很多從根源上避免了“牛眼”對NPP的干擾。結(jié)果與實測對比嚴重偏低時還有另一個容易忽略的因素——太陽總輻射數(shù)據(jù)的單位。CASA模型中輻射量一般要求是MJ/m2但國內(nèi)氣象共享數(shù)據(jù)里常見的輻射單位可能是kW·h/m2或者cal/cm2很多人在數(shù)據(jù)轉(zhuǎn)換時漏了這一步導致能量輸入少了一個數(shù)量級NPP自然偏低。遇到這種情況仔細核對原始數(shù)據(jù)說明文檔統(tǒng)一換算成MJ/m2后再進行計算。5. 這版程序的局限與實際擴展建議V2.0版本在處理常規(guī)區(qū)域尺度的NPP估算時表現(xiàn)穩(wěn)定但它并不是萬能的。一個比較明顯的局限在于它對輸入氣象數(shù)據(jù)的空間分辨率要求較高當數(shù)據(jù)分辨率過粗如超過0.5度時山區(qū)和地形復(fù)雜區(qū)域的NPP會因為氣溫、輻射的均質(zhì)化而明顯失真。若你的研究區(qū)以山地為主建議將氣象數(shù)據(jù)降尺度到500米至1公里的分辨率或者引入DEM做溫度隨海拔遞減修正后再驅(qū)動模型。另一個局限是它對極端氣候事件的響應(yīng)不夠敏感。比如某年夏季發(fā)生了嚴重干旱CASA模型雖然能通過水分脅迫系數(shù)W_ε體現(xiàn)部分影響但如果降水數(shù)據(jù)本身插值精度不高這種響應(yīng)會被抹平。做年際變化分析時建議額外引入帕爾默干旱指數(shù)PDSI或標準化降水蒸散指數(shù)SPEI輔助驗證NPP波動的合理性。關(guān)于程序的擴展我在自己的項目中通常會在V2.0基礎(chǔ)上加兩個外部模塊一個是將逐月NPP結(jié)果自動輸出為NetCDF格式方便后續(xù)用Python做趨勢分析和可視化另一個是接入MODIS land cover動態(tài)產(chǎn)品讓植被類型和ε_max隨年份動態(tài)變化而不是多年使用同一個靜態(tài)植被圖。這兩個擴展都基于V2.0的輸出結(jié)果做二次開發(fā)不影響核心算法感興趣的話后面可以專門寫一篇實現(xiàn)方案。在使用這套程序的道路上我最大的體會就是模型只是工具輸入數(shù)據(jù)的質(zhì)量和前處理細節(jié)才是決定結(jié)果好壞的關(guān)鍵。不要一上來就急著跑結(jié)果花時間把每一層數(shù)據(jù)都檢查到位后面反而能節(jié)省大量調(diào)參和返工的時間。本文還有配套的精品資源點擊獲取