1. 從“綠了”到“綠多久”物候提取的工程化視角每年春天當(dāng)?shù)谝荒ㄐ戮G悄然爬上枝頭或者秋天第一片黃葉飄落我們感知到的就是“物候”。但在生態(tài)學(xué)、農(nóng)業(yè)、氣候研究乃至商業(yè)遙感領(lǐng)域我們需要的遠(yuǎn)不止這種模糊的感性認(rèn)知。我們需要精確的、可量化的、大范圍的“植被物候”數(shù)據(jù)這片森林什么時(shí)候開始生長(zhǎng)生長(zhǎng)季持續(xù)了多久什么時(shí)候達(dá)到生長(zhǎng)頂峰什么時(shí)候進(jìn)入休眠回答這些問(wèn)題就是“植被物候提取”的核心任務(wù)。這聽起來(lái)像是個(gè)純粹的科研問(wèn)題但它的應(yīng)用場(chǎng)景早已滲透到我們生活的方方面面。農(nóng)業(yè)保險(xiǎn)公司需要根據(jù)作物生長(zhǎng)季的長(zhǎng)短來(lái)評(píng)估災(zāi)害風(fēng)險(xiǎn)林業(yè)部門需要監(jiān)測(cè)森林健康狀況預(yù)警病蟲害氣候變化研究者需要全球植被的生長(zhǎng)周期數(shù)據(jù)來(lái)驗(yàn)證模型甚至城市綠化管理部門也需要知道公園里草坪的返青和枯黃時(shí)間以優(yōu)化養(yǎng)護(hù)計(jì)劃。所有這些需求都指向一個(gè)共同的技術(shù)動(dòng)作從海量的、看似雜亂的時(shí)間序列遙感數(shù)據(jù)主要是植被指數(shù)中提取出那幾個(gè)關(guān)鍵的物候期節(jié)點(diǎn)。作為一名長(zhǎng)期與遙感數(shù)據(jù)打交道的從業(yè)者我處理過(guò)從區(qū)域農(nóng)田到全球森林的各種物候提取任務(wù)。市面上論文和教程很多但真正把各種方法的原理、坑點(diǎn)、適用場(chǎng)景講透并能讓你直接“抄作業(yè)”的內(nèi)容卻很少。今天我就拋開那些復(fù)雜的數(shù)學(xué)公式從工程實(shí)踐的角度系統(tǒng)梳理一下幾種最常用、最核心的植被物候提取方法。我會(huì)重點(diǎn)講清楚每種方法“為什么”要那么做參數(shù)該怎么設(shè)以及我最常踩的那些“坑”。無(wú)論你是剛?cè)腴T的學(xué)生還是需要快速應(yīng)用的研究員或工程師這篇文章都能給你一套清晰的“導(dǎo)航圖”。2. 基石與原料理解植被指數(shù)與時(shí)間序列數(shù)據(jù)在討論任何提取方法之前我們必須先徹底理解我們工作的“原料”——時(shí)間序列植被指數(shù)數(shù)據(jù)。這是所有物候提取算法的輸入原料的質(zhì)量直接決定了成品物候參數(shù)的可靠性。2.1 植被指數(shù)的選擇NDVI、EVI 與它們的“表親”植被指數(shù)是通過(guò)衛(wèi)星傳感器捕獲的不同波段主要是紅光和近紅外反射率計(jì)算而來(lái)的用于量化植被的“綠度”或生物物理參數(shù)。最常用的兩位“主角”是NDVI歸一化差異植被指數(shù)(NIR - Red) / (NIR Red)。這是物候研究中的“萬(wàn)金油”歷史數(shù)據(jù)豐富對(duì)中等至高密度植被敏感。但它有個(gè)著名的缺點(diǎn)在植被覆蓋度高的地區(qū)容易飽和且對(duì)大氣影響和土壤背景比較敏感。EVI增強(qiáng)型植被指數(shù)在NDVI的基礎(chǔ)上引入了藍(lán)光波段進(jìn)行大氣校正并加入了土壤調(diào)節(jié)因子。它的優(yōu)勢(shì)在于減少了大氣和土壤背景的噪音在高生物量區(qū)不易飽和對(duì)植被結(jié)構(gòu)變化更敏感。但EVI對(duì)傳感器性能和數(shù)據(jù)預(yù)處理的要求更高。實(shí)操心得對(duì)于大多數(shù)溫帶和熱帶森林、農(nóng)田的物候研究我首選EVI。尤其是在生長(zhǎng)季植被茂密NDVI容易達(dá)到飽和值接近1無(wú)法有效區(qū)分生長(zhǎng)峰值的細(xì)微變化。而EVI的動(dòng)態(tài)范圍更寬能更好地刻畫整個(gè)生長(zhǎng)過(guò)程。但對(duì)于稀疏植被如草原、荒漠或歷史數(shù)據(jù)對(duì)比Landsat系列NDVI仍然是可靠的選擇。一個(gè)簡(jiǎn)單的原則植被越密、研究越關(guān)注生長(zhǎng)峰值細(xì)節(jié)越傾向用EVI數(shù)據(jù)連續(xù)性要求高、植被較稀疏可用NDVI。除了這兩位還有一系列針對(duì)特定場(chǎng)景的指數(shù)NDWI歸一化差異水指數(shù)用于監(jiān)測(cè)植被水分含量在干旱區(qū)物候或與水分脅迫相關(guān)的研究中很有用。GCC綠度色譜坐標(biāo)來(lái)自近地面攝影如物候相機(jī)與NDVI高度相關(guān)但數(shù)據(jù)來(lái)源和尺度完全不同。LAI葉面積指數(shù)、FPAR光合有效輻射吸收比例這些是更直接的生物物理參數(shù)但通常由植被指數(shù)反演而來(lái)本身噪聲可能更大。選擇哪個(gè)指數(shù)首先取決于你的科學(xué)問(wèn)題其次是數(shù)據(jù)的可用性和質(zhì)量。沒(méi)有“最好”只有“最適合”。2.2 時(shí)間序列數(shù)據(jù)的預(yù)處理濾波與平滑是必修課原始的衛(wèi)星時(shí)間序列數(shù)據(jù)是無(wú)法直接用于物候提取的。它充滿了噪聲云、氣溶膠、傳感器誤差、太陽(yáng)高度角變化等等。這些噪聲點(diǎn)會(huì)偽裝成植被的突然生長(zhǎng)或枯萎導(dǎo)致提取算法“誤判”。因此時(shí)間序列重構(gòu)濾波與平滑是物候提取前絕對(duì)不可或缺的一步。這一步的目標(biāo)是去除噪聲還原植被真實(shí)的生長(zhǎng)曲線同時(shí)盡可能保留真實(shí)的物候信號(hào)。常用方法有Savitzky-Golay濾波這是我個(gè)人最常用也最推薦給新手的方。它本質(zhì)上是一個(gè)移動(dòng)窗口多項(xiàng)式擬合。對(duì)于窗口內(nèi)的數(shù)據(jù)點(diǎn)用一個(gè)多項(xiàng)式去擬合然后用擬合出的多項(xiàng)式中心點(diǎn)的值替代原始值。它的優(yōu)點(diǎn)是能有效平滑噪聲同時(shí)較好地保留曲線的峰值、寬度等形態(tài)特征——這些特征對(duì)物候提取至關(guān)重要。關(guān)鍵參數(shù)窗口大小和多項(xiàng)式階數(shù)。窗口越大平滑力度越強(qiáng)但可能過(guò)度平滑抹掉真實(shí)的快速變化如農(nóng)作物收割。我的一般起手式是對(duì)于16天合成的MODIS數(shù)據(jù)窗口大小取4即前后各4個(gè)點(diǎn)共9個(gè)點(diǎn)多項(xiàng)式階數(shù)取2或3。你需要用眼睛觀察平滑后的曲線是否去掉了明顯的“毛刺”同時(shí)又沒(méi)有把生長(zhǎng)季的“駝峰”壓成“平原”。雙邏輯斯蒂函數(shù)擬合這種方法用一個(gè)預(yù)設(shè)的函數(shù)模型兩個(gè)邏輯斯蒂函數(shù)的組合去強(qiáng)行擬合整個(gè)年度的生長(zhǎng)曲線。它基于一個(gè)強(qiáng)假設(shè)植被的生長(zhǎng)和衰亡過(guò)程符合S型曲線。擬合出的函數(shù)非常光滑能直接給出物候參數(shù)如拐點(diǎn)即對(duì)應(yīng)物候期。但缺點(diǎn)是如果某年的物候曲線不符合雙邏輯模型例如受到干旱、火災(zāi)等干擾擬合會(huì)失敗或產(chǎn)生巨大偏差。HANTS時(shí)間序列諧波分析這種方法將時(shí)間序列分解為不同頻率的諧波正弦余弦波保留代表植被生長(zhǎng)周期的低頻信號(hào)剔除代表噪聲的高頻信號(hào)。它對(duì)處理不規(guī)則采樣和數(shù)據(jù)缺失有較好效果但參數(shù)設(shè)置相對(duì)復(fù)雜。踩坑實(shí)錄我曾在一個(gè)項(xiàng)目中為了追求曲線的“光滑”使用了過(guò)大的Savitzky-Golay窗口。結(jié)果在提取北方森林的物候時(shí)把春季返青后一個(gè)短暫的低溫造成的生長(zhǎng)停滯期給平滑掉了導(dǎo)致提取出的生長(zhǎng)季起始日比地面觀測(cè)早了近10天。教訓(xùn)是平滑不是越強(qiáng)越好。一定要將平滑后的曲線與原始數(shù)據(jù)點(diǎn)疊加查看確保關(guān)鍵的轉(zhuǎn)折點(diǎn)春季快速上升、秋季快速下降被保留而不是被平滑掉。對(duì)于噪聲特別大的區(qū)域如熱帶常綠林其NDVI/EVI年際變化本身很小噪聲相對(duì)信號(hào)很強(qiáng)可能需要嘗試結(jié)合多種方法或者接受一定的不確定性。下表對(duì)比了這幾種常用預(yù)處理方法的核心特點(diǎn)方法核心原理優(yōu)點(diǎn)缺點(diǎn)適用場(chǎng)景Savitzky-Golay濾波移動(dòng)窗口局部多項(xiàng)式擬合能保留曲線形態(tài)特征參數(shù)直觀靈活性強(qiáng)對(duì)連續(xù)缺失數(shù)據(jù)敏感窗口參數(shù)需經(jīng)驗(yàn)調(diào)整通用首選尤其適用于溫帶季節(jié)性植被雙邏輯斯蒂擬合用預(yù)設(shè)的S型函數(shù)模型全局?jǐn)M合結(jié)果非常光滑可直接輸出物候參數(shù)物理意義明確假設(shè)性強(qiáng)對(duì)異常年份干旱、災(zāi)害擬合差生長(zhǎng)季規(guī)律明顯的區(qū)域如溫帶農(nóng)田、落葉林HANTS時(shí)間序列的諧波分析與重構(gòu)能處理不規(guī)則數(shù)據(jù)和缺失值分離周期信號(hào)與噪聲參數(shù)頻率數(shù)、擬合誤差容限設(shè)置復(fù)雜不易理解數(shù)據(jù)缺失嚴(yán)重或采樣不規(guī)則的時(shí)間序列預(yù)處理后的干凈時(shí)間序列才是一幅可供我們“作畫”提取物候的潔凈畫布。這一步花費(fèi)的時(shí)間往往能省去后續(xù)大量糾錯(cuò)和解釋的麻煩。3. 閾值法最直觀的“尺子”也是最易誤用的工具閾值法大概是概念上最簡(jiǎn)單、最直觀的物候提取方法了。它的邏輯直白設(shè)定一個(gè)植被指數(shù)的絕對(duì)值或相對(duì)值閾值當(dāng)時(shí)間序列曲線穿過(guò)這個(gè)閾值時(shí)對(duì)應(yīng)的日期就是物候期。例如設(shè)定NDVI0.5為生長(zhǎng)季開始閾值全年時(shí)間序列曲線從低值向上穿過(guò)0.5的那一天即為返青期。3.1 絕對(duì)閾值與相對(duì)閾值的選擇困境絕對(duì)閾值如NDVI0.5。這種方法簡(jiǎn)單粗暴但問(wèn)題巨大。不同生態(tài)系統(tǒng)、不同植被類型、甚至同一地區(qū)不同年份的植被指數(shù)基線都不同。森林的NDVI基線可能就在0.6以上草原可能只有0.3。用一個(gè)固定閾值去套全球顯然會(huì)得出荒謬的結(jié)果。相對(duì)閾值這是更常用的做法。通常定義為生長(zhǎng)季內(nèi)振幅最大值與最小值之差的一個(gè)比例。例如生長(zhǎng)季起始日SOS常被定義為曲線從最低點(diǎn)上升到“振幅的20%”所對(duì)應(yīng)的日期。生長(zhǎng)季結(jié)束日EOS則為曲線從最高點(diǎn)下降到“振幅的20%”的日期。生長(zhǎng)季峰值POS則對(duì)應(yīng)最大值點(diǎn)。相對(duì)閾值法部分解決了生態(tài)系統(tǒng)差異的問(wèn)題因?yàn)樗诒镜鼗恼穹M(jìn)行標(biāo)準(zhǔn)化。但它引入了新的參數(shù)這個(gè)比例該設(shè)多少10%20%50%3.2 閾值法的參數(shù)陷阱與實(shí)戰(zhàn)調(diào)整為什么是20%這其實(shí)沒(méi)有嚴(yán)格的生理學(xué)依據(jù)更多是經(jīng)驗(yàn)值源于早期研究發(fā)現(xiàn)在這個(gè)比例附近曲線變化較陡對(duì)閾值不敏感結(jié)果相對(duì)穩(wěn)定。但實(shí)際應(yīng)用中你需要驗(yàn)證。操作步驟示例以提取SOS為例獲取平滑后時(shí)間序列對(duì)單個(gè)像元一年的EVI數(shù)據(jù)進(jìn)行Savitzky-Golay濾波。計(jì)算本地振幅找出該年度時(shí)間序列的EVI最大值EVI_max和最小值EVI_min。振幅 EVI_max - EVI_min。定義動(dòng)態(tài)閾值閾值 EVI_min 振幅 * 比例系數(shù)。假設(shè)比例系數(shù)設(shè)為0.2。尋找交叉點(diǎn)從年初向年中搜索找到第一個(gè)EVI值大于等于該閾值的日期即為SOS。核心避坑點(diǎn)閾值法最大的敵人是“曲線形態(tài)的多樣性”。我遇到過(guò)幾個(gè)典型問(wèn)題雙峰曲線某些地區(qū)如地中海氣候或多次收割的農(nóng)作物一年內(nèi)可能出現(xiàn)兩個(gè)生長(zhǎng)峰值。閾值法會(huì)找到第一個(gè)上升沿的交叉點(diǎn)但可能會(huì)把第二個(gè)峰誤判為新的生長(zhǎng)季開始或者根本無(wú)法處理。這時(shí)需要先進(jìn)行生長(zhǎng)季分割或者使用更復(fù)雜的方法。平緩的上升沿在高緯度常綠林或熱帶雨林植被指數(shù)年際變化很小上升沿非常平緩。這意味著在閾值附近曲線可能“徘徊”很長(zhǎng)時(shí)間導(dǎo)致提取的SOS日期對(duì)閾值極其敏感。今天用20%閾值是第100天明天用25%閾值可能就是第120天結(jié)果不確定性很大。最小值定位錯(cuò)誤如果年初有積雪或云污染導(dǎo)致EVI_min異常低計(jì)算出的振幅會(huì)異常大從而使閾值虛高嚴(yán)重推遲SOS的提取。因此在計(jì)算振幅前必須合理確定“生長(zhǎng)季背景值”。我通常的做法不是直接用全年最小值而是取生長(zhǎng)季開始前一個(gè)固定窗口如冬眠期的均值作為EVI_min取生長(zhǎng)季峰值附近窗口的均值作為EVI_max這樣可以避免異常值的干擾。閾值法總結(jié)它快速、簡(jiǎn)單、易于實(shí)現(xiàn)和理解是很好的入門方法和快速普查工具。但其結(jié)果嚴(yán)重依賴于預(yù)設(shè)的閾值比例和最大值/最小值的準(zhǔn)確估計(jì)。它適用于生長(zhǎng)季信號(hào)強(qiáng)烈、曲線單峰且陡峭的地區(qū)如溫帶落葉林、一年一熟農(nóng)田。對(duì)于復(fù)雜情況需要謹(jǐn)慎調(diào)整參數(shù)并結(jié)合目視檢查。4. 導(dǎo)數(shù)法尋找變化速度的“轉(zhuǎn)折點(diǎn)”導(dǎo)數(shù)法或稱變化率法從另一個(gè)物理角度切入植被的生長(zhǎng)和衰老不是勻速的在物候期轉(zhuǎn)換時(shí)其生長(zhǎng)速度即植被指數(shù)隨時(shí)間的變化率會(huì)達(dá)到極值。簡(jiǎn)單說(shuō)返青期是“加速”最快的點(diǎn)衰亡期是“減速”最快的點(diǎn)。4.1 一階導(dǎo)數(shù)與二階導(dǎo)數(shù)的物候意義一階導(dǎo)數(shù)表示植被指數(shù)隨時(shí)間的變化速度斜率。曲線上升最快的那一點(diǎn)一階導(dǎo)數(shù)的最大值通常對(duì)應(yīng)生長(zhǎng)加速期可近似視為返青期。曲線下降最快的那一點(diǎn)一階導(dǎo)數(shù)的最小值通常對(duì)應(yīng)衰老加速期可近似視為枯黃期。二階導(dǎo)數(shù)表示變化速度本身的變化率加速度。一階導(dǎo)數(shù)的極值點(diǎn)在二階導(dǎo)數(shù)上對(duì)應(yīng)過(guò)零點(diǎn)。返青期對(duì)應(yīng)二階導(dǎo)數(shù)由正變負(fù)的過(guò)零點(diǎn)從加速增長(zhǎng)變?yōu)闇p速增長(zhǎng)這里需要糾正一個(gè)常見誤解。實(shí)際上對(duì)于一條S型生長(zhǎng)曲線生長(zhǎng)初期速度增加加速度為正。拐點(diǎn)增長(zhǎng)最快點(diǎn)速度達(dá)到最大加速度為零二階導(dǎo)數(shù)過(guò)零點(diǎn)。生長(zhǎng)后期速度減小加速度為負(fù)。 因此生長(zhǎng)季中期峰值增長(zhǎng)期才對(duì)應(yīng)二階導(dǎo)數(shù)的過(guò)零點(diǎn)。而物候期開始和結(jié)束更多與一階導(dǎo)數(shù)的極值點(diǎn)相關(guān)。在實(shí)際算法中我們通常對(duì)平滑后的時(shí)間序列計(jì)算數(shù)值微分如中心差分法然后尋找一階導(dǎo)數(shù)的局部極大值和極小值。4.2 導(dǎo)數(shù)法的實(shí)現(xiàn)細(xì)節(jié)與噪聲放大效應(yīng)數(shù)值微分會(huì)放大噪聲。即使原始數(shù)據(jù)經(jīng)過(guò)平滑微分后仍可能產(chǎn)生許多小的波動(dòng)導(dǎo)致檢測(cè)出多個(gè)虛假的極值點(diǎn)。因此導(dǎo)數(shù)法通常需要與閾值法或規(guī)則結(jié)合使用。一個(gè)常見的復(fù)合策略是用閾值法確定一個(gè)大致的生長(zhǎng)季窗口例如EVI超過(guò)振幅10%到低于振幅10%之間的時(shí)期。在這個(gè)窗口內(nèi)計(jì)算一階導(dǎo)數(shù)序列。在窗口前期尋找一階導(dǎo)數(shù)的最大值其對(duì)應(yīng)日期作為SOS。在窗口后期尋找一階導(dǎo)數(shù)的最小值負(fù)得最多其對(duì)應(yīng)日期作為EOS。實(shí)操技巧直接對(duì)離散數(shù)據(jù)做差分噪聲很大。我通常的做法是先對(duì)平滑后的時(shí)間序列進(jìn)行重采樣例如通過(guò)樣條插值生成一個(gè)更高時(shí)間分辨率如每天的連續(xù)曲線然后再對(duì)這條光滑曲線求導(dǎo)。這樣得到的導(dǎo)數(shù)曲線更干凈極值點(diǎn)更明確。Python中可以用scipy.interpolate進(jìn)行樣條插值再用scipy.misc.derivative求導(dǎo)或者直接對(duì)插值函數(shù)求導(dǎo)。導(dǎo)數(shù)法的優(yōu)勢(shì)在于它基于變化的物理意義不依賴于絕對(duì)的閾值。但它對(duì)數(shù)據(jù)平滑的質(zhì)量要求極高且容易受到生長(zhǎng)季內(nèi)短期波動(dòng)如干旱導(dǎo)致的生長(zhǎng)暫停的干擾這些波動(dòng)也會(huì)產(chǎn)生局部的導(dǎo)數(shù)極值。因此導(dǎo)數(shù)法很少單獨(dú)使用通常作為其他方法如閾值法、曲線擬合法的補(bǔ)充或驗(yàn)證工具用于在生長(zhǎng)季窗口內(nèi)精確定位變化最快的時(shí)刻。5. 曲線擬合法用數(shù)學(xué)模型“概括”生長(zhǎng)季曲線擬合法是物候提取中更為強(qiáng)大和穩(wěn)健的一類方法。其核心思想是用一個(gè)預(yù)設(shè)的、光滑的數(shù)學(xué)模型去描述整個(gè)生長(zhǎng)季的植被指數(shù)變化過(guò)程。擬合成功后物候參數(shù)可以直接從模型的數(shù)學(xué)屬性中推導(dǎo)出來(lái)例如函數(shù)的拐點(diǎn)、極值點(diǎn)、達(dá)到特定比例的日期等。5.1 雙邏輯斯蒂Double Logistic函數(shù)經(jīng)典之選這是最著名、應(yīng)用最廣泛的物候擬合模型。它用兩個(gè)邏輯斯蒂S型函數(shù)分別模擬生長(zhǎng)季的上升返青和下降衰老過(guò)程。其函數(shù)形式大致如下y(t) m1 (m2 - m1) * (1/(1exp(-k1*(t-t1))) 1/(1exp(k2*(t-t2))) - 1)其中m1,m2生長(zhǎng)季開始前和結(jié)束后背景植被指數(shù)水平。t1,k1控制上升過(guò)程返青的拐點(diǎn)日期和速率。t2,k2控制下降過(guò)程衰老的拐點(diǎn)日期和速率。t時(shí)間年積日。擬合這個(gè)模型就是找到一組最優(yōu)參數(shù)(m1, m2, t1, k1, t2, k2)使得函數(shù)曲線y(t)與觀測(cè)到的時(shí)間序列數(shù)據(jù)點(diǎn)最吻合。擬合通常使用非線性最小二乘法如Levenberg-Marquardt算法。物候提取擬合成功后生長(zhǎng)季起始日SOS通常取上升拐點(diǎn)t1或者計(jì)算曲線達(dá)到m1 (m2-m1)*比例的日期。生長(zhǎng)季結(jié)束日EOS通常取下降拐點(diǎn)t2或者計(jì)算曲線達(dá)到m2 - (m2-m1)*比例的日期。生長(zhǎng)季峰值日POS通常取曲線最大值對(duì)應(yīng)的日期約在t1和t2之間。5.2 非對(duì)稱高斯Asymmetric Gaussian函數(shù)更靈活的形態(tài)雙邏輯斯蒂函數(shù)假設(shè)上升和下降是對(duì)稱的S型但實(shí)際植被生長(zhǎng)曲線往往不對(duì)稱例如春季返青可能比秋季衰老更快。非對(duì)稱高斯函數(shù)提供了更大的靈活性它用左右寬度不同的高斯函數(shù)來(lái)擬合能更好地刻畫這種不對(duì)稱性。其函數(shù)形式基于修改的高斯函數(shù)參數(shù)包括峰值位置、峰值高度、左半寬度和右半寬度。擬合和物候提取思路與雙邏輯斯蒂類似。5.3 擬合法的優(yōu)勢(shì)、挑戰(zhàn)與關(guān)鍵步驟優(yōu)勢(shì)抗噪能力強(qiáng)模型擬合過(guò)程本身是一種全局優(yōu)化對(duì)個(gè)別噪聲數(shù)據(jù)點(diǎn)不敏感。結(jié)果物理意義明確參數(shù)t1t2直接對(duì)應(yīng)物候拐點(diǎn)。提供完整曲線擬合出的曲線是完整、光滑的便于后續(xù)分析如積分計(jì)算生長(zhǎng)季總生產(chǎn)力。挑戰(zhàn)與實(shí)操要點(diǎn)初始值猜測(cè)非線性擬合嚴(yán)重依賴于參數(shù)初始值的設(shè)置。給得不好算法可能不收斂或收斂到錯(cuò)誤的局部最優(yōu)解。我通常的初始化策略是m1,m2用生長(zhǎng)季前、后一段時(shí)間如各30天的植被指數(shù)中位數(shù)或均值。t1用閾值法如20%振幅粗略估計(jì)的SOS。t2用閾值法粗略估計(jì)的EOS。k1,k2設(shè)為經(jīng)驗(yàn)值如0.1-0.5表示變化速率。擬合失敗處理不是所有時(shí)間序列都能被完美擬合。對(duì)于常綠林曲線平坦、遭受干擾火災(zāi)、砍伐或云污染嚴(yán)重的序列擬合可能失敗。必須設(shè)置嚴(yán)格的擬合優(yōu)度檢驗(yàn)如R平方低于0.6或殘差過(guò)大則標(biāo)記該像元擬合失敗采用備用方法如閾值法或直接標(biāo)記為無(wú)效數(shù)據(jù)。參數(shù)邊界約束必須給參數(shù)設(shè)置合理的物理邊界。例如t1必須早于t2且在一年內(nèi)k1k2必須為正數(shù)m2生長(zhǎng)季峰值水平應(yīng)大于m1背景水平。這些約束能極大提高擬合的穩(wěn)定性和合理性。深度踩坑我曾用雙邏輯斯蒂函數(shù)批量擬合全球數(shù)據(jù)。在赤道熱帶雨林地區(qū)失敗率異常高。原因是熱帶雨林的植被指數(shù)年循環(huán)幅度很小曲線近乎一條直線加噪聲。雙邏輯斯蒂模型試圖去擬合一條“S型”曲線但數(shù)據(jù)中根本不存在這樣的強(qiáng)信號(hào)導(dǎo)致擬合結(jié)果完全隨機(jī)。解決方案是先計(jì)算時(shí)間序列的振幅最大值-最小值如果振幅小于一個(gè)經(jīng)驗(yàn)閾值例如對(duì)于MODIS EVI小于0.1則直接認(rèn)為該地區(qū)無(wú)顯著季節(jié)性物候跳過(guò)擬合或賦予其特殊標(biāo)識(shí)。這叫“信號(hào)強(qiáng)度檢測(cè)”是批量處理前必不可少的一步。曲線擬合法提供了更穩(wěn)健、理論上更優(yōu)美的物候提取方案尤其適合處理中等噪聲水平、具有明顯單峰季節(jié)性的數(shù)據(jù)。它是目前許多全球物候產(chǎn)品如MODIS MCD12Q2的核心算法。6. 物候提取的完整工作流與質(zhì)量評(píng)估掌握了核心方法我們需要將其串聯(lián)成一個(gè)自動(dòng)化、可批量處理、且包含質(zhì)量控制的完整工作流。這對(duì)于處理海量遙感數(shù)據(jù)至關(guān)重要。6.1 一個(gè)穩(wěn)健的物候提取流水線設(shè)計(jì)以下是一個(gè)我常用的、結(jié)合了多種方法優(yōu)勢(shì)的混合流水線以單個(gè)像元多年時(shí)間序列為例數(shù)據(jù)準(zhǔn)備與預(yù)處理輸入原始時(shí)間序列植被指數(shù)如MODIS 16天合成EVI、對(duì)應(yīng)的數(shù)據(jù)質(zhì)量標(biāo)識(shí)QA波段。利用QA波段進(jìn)行初步去云和去低質(zhì)量數(shù)據(jù)將低置信度數(shù)據(jù)標(biāo)記為缺失。使用Savitzky-Golay濾波對(duì)時(shí)間序列進(jìn)行平滑重構(gòu)填補(bǔ)缺失值生成連續(xù)光滑的曲線。生長(zhǎng)季背景值估算針對(duì)每一年定義“非生長(zhǎng)季”窗口如北半球溫帶前一年第300天至當(dāng)年第60天當(dāng)年第300天至第365天。取這些窗口內(nèi)有效數(shù)據(jù)的中位數(shù)作為該年的背景值EVI_bg。使用中位數(shù)是為了抵抗異常值。年度信號(hào)分割與振幅計(jì)算針對(duì)每一年在平滑曲線上尋找全局最大值EVI_max。計(jì)算年度振幅Amp EVI_max - EVI_bg。信號(hào)強(qiáng)度檢查如果Amp 閾值如0.1標(biāo)記該像元該年為“無(wú)顯著季節(jié)”物候參數(shù)賦空值流程結(jié)束。粗略生長(zhǎng)季窗口劃定閾值法使用相對(duì)閾值法如10%振幅確定生長(zhǎng)季大致的開始和結(jié)束范圍[SOS_rough, EOS_rough]。這個(gè)窗口用于約束后續(xù)精細(xì)提取。精細(xì)物候提取曲線擬合法為主在[SOS_rough-30天 EOS_rough30天]的擴(kuò)展窗口內(nèi)使用雙邏輯斯蒂函數(shù)或非對(duì)稱高斯函數(shù)進(jìn)行擬合。提供精心設(shè)置的參數(shù)初始值和邊界約束。計(jì)算擬合優(yōu)度R2。如果R2 0.7可調(diào)認(rèn)為擬合成功。從擬合函數(shù)中提取物候參數(shù)例如將上升拐點(diǎn)t1作為SOS下降拐點(diǎn)t2作為EOS函數(shù)最大值點(diǎn)作為POS。擬合失敗的后備方案如果擬合失敗R2過(guò)低或不收斂則回退到導(dǎo)數(shù)法或閾值法。在粗略生長(zhǎng)季窗口內(nèi)計(jì)算一階導(dǎo)數(shù)分別尋找上升段的最大值和下降段的最小值作為SOS和EOS的備選。結(jié)果后處理與過(guò)濾時(shí)間連續(xù)性檢查同一像元相鄰年份的SOS或EOS日期不應(yīng)發(fā)生劇烈跳躍如超過(guò)30天。對(duì)于異常跳躍值需要結(jié)合上下文判斷是否為真實(shí)變化如火災(zāi)或提取錯(cuò)誤?？臻g一致性檢查相鄰像元的物候日期應(yīng)具有空間連續(xù)性。孤立的、與周邊差異極大的值可能是錯(cuò)誤提取可以考慮用中值濾波等方法進(jìn)行平滑或剔除。6.2 如何評(píng)估你提取的物候數(shù)據(jù)質(zhì)量沒(méi)有評(píng)估結(jié)果就不可信。評(píng)估分為間接驗(yàn)證和直接驗(yàn)證。間接驗(yàn)證內(nèi)部一致性時(shí)間序列可視化隨機(jī)抽取一批像元將原始數(shù)據(jù)點(diǎn)、平滑曲線、擬合曲線以及提取的物候期標(biāo)記豎線畫在同一張圖上。人工目視檢查提取的點(diǎn)是否落在曲線的“合理”位置如上升沿中部、下降沿中部?？臻g分布圖將提取的SOS或EOS制成空間分布圖。檢查是否符合地理規(guī)律如緯度梯度、海拔梯度和生態(tài)系統(tǒng)分布規(guī)律落葉林早于針葉林農(nóng)田有獨(dú)特模式。出現(xiàn)大面積反常識(shí)的斑塊很可能算法在該區(qū)域失效。統(tǒng)計(jì)分布查看物候日期的直方圖。如果出現(xiàn)不合理的雙峰或多峰可能意味著算法對(duì)某些地類如農(nóng)田、混合像元處理不佳。直接驗(yàn)證與地面真值對(duì)比 這是最可靠但最困難的方式。需要獲取地面物候觀測(cè)數(shù)據(jù)如物候相機(jī)網(wǎng)絡(luò)、人工觀測(cè)記錄。尺度匹配問(wèn)題地面觀測(cè)是一個(gè)點(diǎn)衛(wèi)星像元是一個(gè)面如MODIS是500x500米。如果觀測(cè)點(diǎn)位于均質(zhì)植被內(nèi)如大片農(nóng)田中心匹配較好如果位于森林邊緣或城市匹配誤差會(huì)很大。物候定義對(duì)齊問(wèn)題地面觀測(cè)的“展葉始期”和衛(wèi)星提取的“生長(zhǎng)季起始日”在生理意義上并不完全等同。衛(wèi)星看到的是冠層整體的綠度變化。需要理解這種差異通常衛(wèi)星物候會(huì)稍晚于地面展葉期。常用指標(biāo)使用均方根誤差RMSE、偏差Bias、相關(guān)系數(shù)R來(lái)定量評(píng)估。通常在均質(zhì)植被區(qū)RMSE能控制在7-15天以內(nèi)就可以認(rèn)為算法性能不錯(cuò)。經(jīng)驗(yàn)之談在實(shí)際項(xiàng)目中目視檢查永遠(yuǎn)是最重要、最不能省略的一環(huán)。無(wú)論你的算法多么自動(dòng)化在批量運(yùn)行前一定要在不同生態(tài)系統(tǒng)、不同氣候區(qū)隨機(jī)抽取上百個(gè)像元進(jìn)行人工檢查。我經(jīng)常通過(guò)這種檢查發(fā)現(xiàn)一些意想不到的算法邊界情況比如在灌溉農(nóng)田區(qū)由于多次澆水EVI曲線出現(xiàn)多次小波動(dòng)導(dǎo)致擬合失敗。這些發(fā)現(xiàn)是優(yōu)化算法、增加規(guī)則的最直接來(lái)源。物候提取不是純數(shù)學(xué)問(wèn)題更是對(duì)生態(tài)系統(tǒng)過(guò)程的理解問(wèn)題。