本發(fā)明涉及土地恢復(fù)分析技術(shù)領(lǐng)域，特別是一種多源時序影像的稀土礦區(qū)土地毀損與恢復(fù)分析方法。

背景技術(shù)：

離子型稀土礦于1969年在我國江西首次發(fā)現(xiàn)，具有分布廣、儲量豐富、放射性低、稀土配分齊全，富含中重稀土元素等優(yōu)勢，是我國目前特有的稀土礦產(chǎn)資源，也是世界上稀缺的礦種，廣泛分布在我國南方的江西、福建、湖南、廣東、廣西等省，其中江西省所占份額最大。其開采工藝主要有池浸、堆浸和原地浸礦三種。池浸、堆浸兩種開采工藝由于需要將礦體和表土植被剝離，導(dǎo)致礦區(qū)植被遭受了巨大破壞，并且稀土開采過程中會產(chǎn)生大量尾砂，隨著雨水流向四周，導(dǎo)致土地沙化，河道堵塞等問題；原地浸礦開采工藝相比池/堆浸開采工藝，對礦區(qū)植被破壞更小，但由于開采過程中需要向山體中注入大量的浸礦液體，且浸礦液體需長時間浸泡山體，容易導(dǎo)致土壤酸化、山體滑坡等問題。不同的稀土開采工藝導(dǎo)致的土地毀損狀況不一樣，其恢復(fù)過程也不盡相同，因此，了解不同稀土開采工藝下的土地毀損與恢復(fù)特征，為稀土礦區(qū)生態(tài)環(huán)境治理，合理有序開采提供指導(dǎo)，對稀土礦區(qū)的可持續(xù)發(fā)展來說顯得尤為重要。

目前，礦區(qū)土地毀損與恢復(fù)分析方法主要有基于實測數(shù)據(jù)的土地毀損與恢復(fù)評價模型和基于遙感數(shù)據(jù)的土地毀損與恢復(fù)過程動態(tài)監(jiān)測兩種，前一種方法不僅需要大量的實測數(shù)據(jù)，費時費力，且針對不同地點需要建立不同的模型，應(yīng)用推廣較為困難；受技術(shù)方法和經(jīng)費限制，早期的土地毀損與恢復(fù)過程動態(tài)監(jiān)測多采用少數(shù)幾幅間隔較長無云的遙感影像進行，研究方法為多時相分類后比較法，該方法能夠反映地物類型變化在空間上的分布情況，但地物類型變化的過程了解較少；遙感數(shù)據(jù)源的不斷豐富以及遙感時序分析法的發(fā)展，為改進對更廣的時空范圍生態(tài)過程的理解提供了重大機遇。時序分析法通過量化分析地物在時相上的變化趨勢和規(guī)律來獲取地物的變化情況，它對變化的判斷標(biāo)準(zhǔn)具有較高的一致性，算法效率較高。

中國發(fā)明專利申請cn105598143a公開了一種資源化礦山修復(fù)方法，將尾礦庫中的尾礦資源和/或剝巖區(qū)/排土場的礦石資源進行有價元素提??；無再提取價值的尾礦資源和/或礦石資源、有價元素提取工序產(chǎn)生的尾礦如適合深加工，則進行深加工獲得目標(biāo)產(chǎn)品；如果不適合深加工，則充填至采空區(qū)/塌陷區(qū)；對剝巖區(qū)/排土場的占用土地和/或工業(yè)場地進行無害化處理；對無害化處理后的剝巖區(qū)/排土場的占用土地和/或工業(yè)場地、充填后的采空區(qū)/塌陷區(qū)進行生態(tài)修復(fù)。

技術(shù)實現(xiàn)要素：

本發(fā)明需要解決的技術(shù)問題是一種高效的多源時序影像的稀土礦區(qū)土地毀損與恢復(fù)分析方法。

為解決上述的技術(shù)問題，本發(fā)明的一種多源時序影像的稀土礦區(qū)土地毀損與恢復(fù)分析方法，包括以下步驟，

步驟s1：研究區(qū)域數(shù)據(jù)獲取，通過數(shù)據(jù)獲取單元獲取研究區(qū)域所需的原始遙感影像；

步驟s2：數(shù)據(jù)預(yù)處理，通過對獲得的原始遙感影像進行輻射校正、大氣校正、幾何校正、圖像裁剪得到研究區(qū)域的遙感影像；

步驟s3：植被指數(shù)計算，通過植被指數(shù)計算單元計算預(yù)處理后影像的ndvi；

步驟s4：轉(zhuǎn)換方程構(gòu)建，并對轉(zhuǎn)換方程的精度進行檢驗；

步驟s5：nvdi轉(zhuǎn)換，確定時序分析法相關(guān)參數(shù)，土地毀損與恢復(fù)定量分析。

進一步的，步驟s1中所述原始遙感影像包括空間分辨率均為30m的landsat5tm、landsat8oli和hj-1bccd遙感影像。

更進一步的，所述步驟s4具體包括以下步驟，

步驟s41：轉(zhuǎn)換方程構(gòu)建，基于同時期的landsat5和hj-1bccd數(shù)據(jù)、landsat8和hj-1bccd數(shù)據(jù)，構(gòu)建landsat5和hj-1bccd數(shù)據(jù)、landsat8和hj-1bccd數(shù)據(jù)的ndvi之間的轉(zhuǎn)換方程；

步驟s42：轉(zhuǎn)換方程的精度檢驗，基于獲得的轉(zhuǎn)換方程將hj-1bccd數(shù)據(jù)的ndvi轉(zhuǎn)化為與之時相相同的landsat5或landsat8的ndvi,并與真實的landsat5或landsat8比較，采用均方根誤差的方法檢驗轉(zhuǎn)換方程的精度。

更進一步的，所述步驟s5具體包括以下步驟，

步驟51：ndvi轉(zhuǎn)換，基于獲得的轉(zhuǎn)換方程將時序影像中l(wèi)andsat8和hj-1bccd數(shù)據(jù)的ndvi轉(zhuǎn)化為landsat5標(biāo)準(zhǔn)下，并將ndvi時序影像組合成一個多光譜影像；

步驟52：時序分析法相關(guān)參數(shù)獲取，通過數(shù)據(jù)挖掘的方法獲得時序分析法相關(guān)參數(shù)；

步驟53：土地毀損與恢復(fù)定量分析，根據(jù)時序分析法相關(guān)參數(shù)對土地毀損與恢復(fù)類型進行分類，獲得土地毀損與恢復(fù)類型空間分布，用于定量分析稀土開采擾動下的土地毀損與恢復(fù)特征。

采用上述方法后，本發(fā)明能實時、快速、大面積的獲取稀土礦區(qū)的地表信息，相比傳統(tǒng)根據(jù)實測數(shù)據(jù)建模分析方法，效率更高，應(yīng)用更為簡單，推廣更容易。另外，本發(fā)明采用時序分析法分析稀土開采環(huán)境下的土地毀損與恢復(fù)過程的相關(guān)特征，使得相關(guān)特征的描述更為細(xì)化。

附圖說明

下面結(jié)合附圖和具體實施方式對本發(fā)明作進一步詳細(xì)的說明。

圖1為本發(fā)明一種多源時序影像的稀土礦區(qū)土地毀損與恢復(fù)分析方法的流程圖。

圖2為本發(fā)明一種多源時序影像的稀土礦區(qū)土地毀損與恢復(fù)分析方法的結(jié)構(gòu)示意圖。

圖3a為本發(fā)明中l(wèi)andsat5與hj-1bccd散點圖。

圖3b為本發(fā)明中l(wèi)andsat8與hj-1bccd散點圖。

圖3c為本發(fā)明中l(wèi)andsat5與模擬影像散點圖。

圖3d為本發(fā)明中l(wèi)andsat8與模擬影像散點圖。

圖4a為本發(fā)明中暫未開采和無植被覆蓋像元類型的ndvi軌跡圖。

圖4b為本發(fā)明中采后無植被和低植被覆蓋像元類型的ndvi軌跡圖。

圖4c為本發(fā)明中采后植被覆蓋達(dá)到采前水平像元類型的ndvi軌跡圖。

圖5為本發(fā)明中土地毀損與恢復(fù)類型分布圖。

圖6a-圖6f為本發(fā)明中類型1對應(yīng)的高分影像圖。

圖7a-圖7f為本發(fā)明中類型2對應(yīng)的高分影像圖。

圖8a-圖8f為本發(fā)明中類型3對應(yīng)的高分影像圖。

圖9a-圖9f為本發(fā)明中類型4對應(yīng)的高分影像圖。

圖10a-圖10f為本發(fā)明中類型5對應(yīng)的高分影像圖。

圖11為本發(fā)明中稀土開采擾動至植被恢復(fù)時長比較圖。

具體實施方式

如圖1和圖2所示，本發(fā)明一種多源時序影像的稀土礦區(qū)土地毀損與恢復(fù)分析方法，包括以下步驟，

步驟s1：研究區(qū)域數(shù)據(jù)獲取，通過數(shù)據(jù)獲取單元獲取研究區(qū)域所需的原始遙感影像；，包括空間分辨率均為30m的landsat5tm、landsat8oli和hj-1bccd遙感影像。

步驟s2：數(shù)據(jù)預(yù)處理，將步驟s1中獲得的遙感數(shù)據(jù)通過輻射定標(biāo)模塊對原始影像進行輻射定標(biāo)，輻射定標(biāo)后使用大氣校正模塊進行大氣校正，大氣校正后使用幾何校正模塊校正影像，幾何校正后使用圖像裁剪模塊獲得研究區(qū)的遙感影像。輻射定標(biāo)模塊中l(wèi)andsat5和landsat8數(shù)據(jù)使用的envi自帶radiometriccalibration工具進行，hj-1bccd數(shù)據(jù)使用插件工具進行；大氣校正模塊使用envi自帶的flaash模塊進行；幾何校正模塊中，以landsat數(shù)據(jù)為基準(zhǔn)，hj-1bccd數(shù)據(jù)為待校正，校正方法選擇二次多項式，對hj-1bccd數(shù)據(jù)進行校正；最后通過研究區(qū)域的矢量邊界對校正后的影像進行裁剪得到研究區(qū)域的遙感影像。

步驟s3：植被指數(shù)計算，通過植被指數(shù)計算單元計算預(yù)處理后影像的ndvi。獲得研究區(qū)域的ndvi影像。植被指數(shù)計算單元通過envi自帶的bandmath工具進行，計算公式如下：

式中ρnir和ρred分別表示遙感影像近紅外、紅光波段的反射率。

步驟s4：轉(zhuǎn)換方程構(gòu)建，并對轉(zhuǎn)換方程的精度進行檢驗。

步驟s41：轉(zhuǎn)換方程構(gòu)建，基于同時期的landsat5和hj-1bccd數(shù)據(jù)、landsat8和hj-1bccd數(shù)據(jù)，采用回歸分析的方法，構(gòu)建landsat5和hj-1bccd數(shù)據(jù)、landsat8和hj-1bccd數(shù)據(jù)的ndvi之間的轉(zhuǎn)換方程。

步驟s42：轉(zhuǎn)換方程的精度檢驗，基于獲得的轉(zhuǎn)換方程將hj-1bccd數(shù)據(jù)的ndvi轉(zhuǎn)化為與之時相相同的landsat5或landsat8的ndvi,并與真實的landsat5或landsat8比較，采用均方根誤差的方法檢驗轉(zhuǎn)換方程的精度。

計算公式如下:

式中xi和yi分別表示為hj-1bccd數(shù)據(jù)校正后的ndvi值、landsat5或landsat8的ndvi值，n為樣本個數(shù)。

步驟5：ndvi轉(zhuǎn)換，確定時序分析法相關(guān)參數(shù)，并根據(jù)相關(guān)參數(shù)對研究區(qū)域的土地毀損與恢復(fù)類型進行分類，得到不同類型在空間上的分布，從而對稀土開采擾動下的土地毀損與恢復(fù)特征進行分析。

步驟51：ndvi轉(zhuǎn)換。通過獲得的轉(zhuǎn)換方程將時序影像中l(wèi)andsat8、hj-1bccd影像的ndvi轉(zhuǎn)化為landsat5標(biāo)準(zhǔn)下，通過envi自帶的layerstacking工具將ndvi時序影像按照時間順序組合成一個類似多光譜影像的數(shù)據(jù)。

步驟52：確定時序分析法參數(shù)。通過目視和比較得到樣本點時序ndvi的變化軌跡，如圖4a-圖4c。研究區(qū)域范圍內(nèi)包含以下5種類型，1)暫未開采，整個時期植被覆蓋均保持在較高水平；2)池浸/堆浸開采，采前有植被，采后植被尚未恢復(fù)，始終維持在較低水平；3)始終無植被覆蓋或植被覆蓋水平較低；4)植被覆蓋水平緩慢升高并恢復(fù)至周邊同類區(qū)域采前水平；5)采前有植被，采后植被開始恢復(fù)，但仍低于原始水平。在時序ndvi影像上選擇一定數(shù)量的五種類型的樣本點，應(yīng)用cart(classificationandregressiontree)決策樹分類方法，將采集的樣本點輸入rulegenv1.02軟件中，算法選擇quest，得到有植被與基本無植被、已開采與未開采的分類閾值，分別為tr1＝0.5285,tr2＝0.3787。

步驟53：土地毀損與恢復(fù)定量分析。根據(jù)得到的時序分析法的相關(guān)參數(shù)tr1,tr2，按照表1的標(biāo)準(zhǔn)，對研究區(qū)域的土地毀損與恢復(fù)類型進行劃分，得到土地毀損與恢復(fù)類型空間分布圖及礦區(qū)毀損時間分布。

表1稀土礦區(qū)土地毀損與恢復(fù)類型劃分

表格中mini表示整個觀測期內(nèi)時序ndvi的最低值，反映i像元所在位置是否受稀土開采擾動；maxi表示整個觀測期內(nèi)像元時序ndvi的最高值，反映i像元所在位置多年植被覆蓋的峰值；maxi_pre表示受稀土開采擾動前時序ndvi的最高值，反映i像元所在位置采前植被覆蓋的峰值；maxi_post受稀土開采擾動后時序ndvi的最高值，反映i像元所在位置采后植被覆蓋的峰值；ri表示植被恢復(fù)水平，反映i像元區(qū)域采后植被覆蓋是否達(dá)到或超過采前水平，ri＝maxi_post-maxi_pre,ri≥0表示采后植被覆蓋水平達(dá)到或超過采前水平；ri＜0表示采后植被覆蓋尚未達(dá)到采前水平。

以贛州市定南縣嶺北稀土礦區(qū)為例，所采用的多源時序影像如表2所示，所采用的交互比較影像對如表3所示。

表2時序影像日期、類型和標(biāo)識符

表3交互比較影像對

如圖4a-4c所示，類型1為暫未開采，類型2為采后無植被，類型3為無植被覆蓋，類型4為采后植被覆蓋達(dá)到采前水平，類型5為采后植被得到一定程度恢復(fù)，但未達(dá)到采前水平。

通過圖5和圖11，分析稀土開采擾動下的土地毀損與恢復(fù)過程。結(jié)合googleearth高分影像發(fā)現(xiàn)，類型1多為未受稀土開采影響的林地，如圖6a-圖6f所示；類型2為稀土開采區(qū)域且未進行礦區(qū)復(fù)墾，如圖7a-圖7f所示；類型3多為農(nóng)田、建筑等，如圖8a-圖8f所示；類型4為多為退耕還林、果園種植區(qū)，如圖9a-圖9f所示；類型5為稀土開采區(qū)域但進行了礦區(qū)復(fù)墾，如圖10a-圖10f所示。

雖然以上描述了本發(fā)明的具體實施方式，但是本領(lǐng)域熟練技術(shù)人員應(yīng)當(dāng)理解，這些僅是舉例說明，可以對本實施方式作出多種變更或修改，而不背離本發(fā)明的原理和實質(zhì)，本發(fā)明的保護范圍僅由所附權(quán)利要求書限定。