吳鳳敏 余 靜 王素偉 王 琛 鄭稚棚 李卓錕 張靈犀
(1. 重慶市地理信息和遙感應用中心, 重慶 401147;2. 重慶市地質礦產研究院, 重慶 400042)
礦產資源開發在國家社會經濟發展中具有不得替代的作用,但隨著礦山的開采,會造成生態環境破壞,導致各種各樣的環境問題發生,因此需要開展綜合治理[1]。我國礦山地質環境恢復和綜合治理任務要求到2025年,全面建立動態監測體系[3],構建礦山恢復治理的遙感監測體系是礦山生態環境保護與治理恢復任務的重要內容,是加強礦山治理過程動態監管和治理效果評價的重要依據[4]。
近年來,針對礦山修復治理遙感監測已有大量研究,部分學者根據遙感影像,建立不同傳感器、不同分辨率礦山及內部結構(尾礦庫、煤矸石堆等)解譯標志庫[5-7],便于礦山的識別。還有學者利用中分辨率遙感影像,如Landsat TM(thematic mapper)等對礦區礦產開發對土地、植被等生態環境的影響進行研究[8-9]。后續越來越多學者利用高分辨率遙感影像對礦山目標和應用進行相關研究[10-11],如杜培軍分析了高分辨率遙感圖像在礦山的應用情況,并研究了開展應用的關鍵問題[12],薛衛寧對礦山治理修復前后的地質環境變化和治理修復效果進行了解譯[13],鄧錦山等人采用2016—2018年高分辨率遙感影像建立解譯標志,構建礦山地質環境評價指標體系并進行綜合評價[14],周英杰等人利用遙感衛星圖像對山東省尾礦庫的遙感監測與分析[15]。此外,還有開展生態環境恢復治理規劃設計實景三維平臺建設研究,利用三維地理信息系統(geographic information system,GIS)技術建立礦山真三維場景,露天礦山環境恢復治理規劃方案、設計評審、政府決策提供智能化支撐[16]。
已有對礦山恢復治理的監測研究中,有大尺度中低分辨率影像礦山識別方法,有礦山高分辨率遙感影像人機交互式解譯方法,還有實景三維傾斜攝影應用等多個方面,但對于礦山恢復治理遙感監測技術沒有進行綜合運用。基于以上背景,根據重慶實際情況,結合中分辨率遙感、高分辨率等遙感數據資源,提出重慶市礦山修復治理遙感動態監測的技術方法,以便為后續礦山修復工程動態監管提供基礎。
對于礦山恢復治理過程中遙感動態監測技術研究,首先是礦山范圍的識別研究,通過利用高分辨率遙感影像,建立礦山遙感解譯樣本庫,對礦山范圍采用人機交互式遙感解譯提取礦山損毀土地范圍;其次對礦山植被狀況進行動態監測研究,采用中分辨率遙感影像,基于植被在近紅外波段的光譜響應特征構建植被指數,對礦山植被變化進行動態監測;研究采用哨兵數據,通過干涉分析對礦山恢復治理過程中地表形變狀況進行監測分析;同時,通過無人機采集的高分辨率遙感影像,進行三維模型構建,為礦山恢復治理三維實景效果提供數據支撐,主要技術流程如下(圖1)。

圖1 露天開采礦山恢復治理遙感動態監測技術體系
研究采用國外衛星影像如WorldView、Pleiades,國產衛星影像如高分2號、高景1號等高分辨率遙感影像進行礦山損毀土地范圍的精確識別。不同遙感影像在數據處理方面具有一定差異,衛星遙感影像一般處理流程包括:影像配準、影像融合、正射校正、影像增強(影像調色)、去云處理、影像鑲嵌等(圖2)。

圖2 衛星遙感影像數據處理流程
影像配準以選取同名點對多光譜影像進行糾正,多光譜糾正模型及數字高程模型(digital elevation model,DEM)數據選擇應與全色波段一致,兩景影像之間配準精度不得大于1個像素,典型地物和地形特征(如山谷、山脊)不能有重影。影像融合目的是保持影像紋理清晰、色彩自然,無虛像和重影現象,研究采用Pansharpen或高通濾波融合方法。正射糾正采用雙線性內插或立方卷積插值重采樣方式,糾正后正射影像應盡量避免影像拉伸和扭曲現象。影像坐標投影方式采用標準規定的坐標系統和投影方式以及分帶,中央經線的選取按照面積最大原則。
影像增強包括對比度/色彩飽和度調整、勻光勻色處理以及銳化處理等,使影像能真實反映地表特征和自然色調,影像對比度和色彩飽和度調整采用濾波和直方圖拉伸方法。影像鑲嵌時,應保持接邊處色彩過渡自然,地物合理接邊,無重影和發虛現象。如鑲嵌區內有人工地物時,拼接線應繞開人工地物,使鑲嵌結果保持人工地物的完整性和合理性。
礦山損毀土地范圍識別可結合高分辨率遙感影像,通過建立礦山的遙感樣本庫,獲取其光譜和紋理特征,通過人機交互式人工目視解譯,形成礦山損毀土地范圍成果。重慶露天開采類礦山主要為露天石灰巖、建筑石料用灰巖等,其中石灰巖開采面在遙感影像上為亮白色,一般呈現不規則的團狀或長條狀,與周邊的色調、紋理有明顯差異[5]。恢復治理過程中關閉石材廠,開采面主要成暗灰色,紋理也較為光滑,周邊一般植被較好,且部分區域由于植被恢復,呈現淺綠色低植被覆蓋或者淺棕色土壤裸露等(圖3)。此外,礦山建筑一般位于采砂采石場周邊區域,建筑數量不多或呈零星分布,遙感影像上多為灰色房屋或藍色鋼棚。

圖3 礦山影像特征示例——2020年北碚區天府礦務局劉家灣煤礦
礦山損毀土地范圍識別采用人機交互式目視解釋方式,通過建立遙感影像特征樣本,認為經驗判斷礦山,并勾畫礦山范圍,同時結合礦山相關數據,如礦業權數據,實地調查數據以及其他專題資料數據(地理國情監測、第三次國土調查)等,進一步對礦山判斷和識別礦山損毀土地范圍。
礦山恢復治理植被遙感動態監測,由于需要多年數據,考慮到監測成本、分辨率和監測時效性,研究采用landsat數據用于植被監測。為更好研究礦山恢復治理植被狀況,選取植被主要生長季節數據(5~10月)進行分析,重慶由于部分數據夏季云量較大,選擇鄰近夏季影像。數據處理主要包括輻射定標、大氣校正、正射校正等流程。輻射定標處理,主要指將影像的灰度(D)值轉化成輻射亮度值或表觀反射率。研究將影像灰度值轉化為輻射亮度值,具體公式為
(1)
式中,L為波段輻射能量;D為像元亮度值;Gain為增量校正系數;Offest為校正偏差量。
大氣校正主要是將輻射定標后的表觀反射率或輻射亮度值轉化為地表反射率的過程,研究采用FLAASH進行大氣校正,計算公式為
(2)
式中,L為總輻射亮度;ρ為像素表面反射率;ρe為像素周圍平均表觀反射率;S為大氣球面反照率;La為大氣球面反射;A、B取決于大氣條件和幾何條件。
正射校正需基于DEM來對影像中所有的像元執行地形變形的校正,讓圖像滿足正射投影的要求,研究可根據影像自帶的遠程過程調用(remote procedure call,RPC)參數和遙感圖像處理平臺(the environment for visualizing images,ENVI)中的DEM信息,對影像數據進行正射校正。
研究采用目前應用較為廣泛的哨兵數據(Sentinel-1)對礦山地表形變進行監測,該數據干涉模式幅寬為250 km×250 km,分辨率為30 m,單個Sentinel-1衛星每12 d覆蓋全球一次,兩個衛星星座的重復周期可縮短至6 d,且數據能夠免費獲取,有利于開展時間序列的礦山形變監測。Sentinel-1數據一般利用GAMMA軟件編寫數據預處理腳本,在數據預處理基礎上,采用小基線集技術手段,進行時間系列干涉分析,精確反演線性形變速率(圖4)。

圖4 哨兵數據主要處理技術流程
為最大化所有干涉對象間的相干性,需要選擇一幅主影像。已有研究表明,干涉像對間的相干值依賴于時間基線(T)、空間垂直基線(B⊥)、多普勒中心基線(FDC)和熱噪聲四部分乘積[19],干涉對總體相關性計算公式為
(3)

在選擇主影像的基礎上,對雷達數據進行高精度配準研究。由于Sentinel-1數據獨特的循序掃描地形觀測(terrain observation by progressive scans,TOPS)成像模式,獲取到的每組Burst的多普勒頻率沿方位向上不斷變化,為避免方位向上相位不連續帶來的影響,采用基于增強光譜屬性的Sentinel-1數據配準方法,提高配準精度。為獲得相干相位中的形變相位,必須消除參考相位和地形相位成分,其中參考相位可利用幾何姿態和相關成像參數,通過多項式擬合去除;地形相位根據合成孔徑雷達(synthetic aperture radar,SAR)觀測數據和已有數字高程模型(DEM),通過二軌差分法去除。二軌差分法主要利用形變前后的SAR數據組成干涉對,干涉得到包含形變值的干涉相位,然后依據已有高程信息和SAR成像參數模擬參考相位和地形相位,將干涉相位和模擬的相位進行差分得到地表形變相位,進而將形變相位轉換為形變值(圖5)。

圖5 二軌差分干涉測量處理流程圖
由于定軌數據精度有限,不可能完全去除軌道誤差影響,因此在實際應用中必須去除軌道誤差影響,軌道誤差影響在整個干涉區域內呈現系統性的趨勢,研究采用二次曲面擬合方法進行去除[20-21],設二次曲面模型為
(4)
式中,φ為軌道殘差相位;i,j為干涉圖像元位置;a0,a1,a2,a3,a4,a5為二次曲面模型系數。已知模型系數就可求出干涉圖中每點的軌道殘差相位,其中,Φ=MA。
解纏后的差分干涉相位計算公式為

相干點雖然在時間序列中保持高的相干性,但是不同時期相干值不一樣,對應解纏的干涉相位也存在差異,且可能存在解纏錯誤情況,研究采用L1范求解,即:
(12)
假定A滿秩,采用下式進行迭代求解:
(13)
其中P為權矩陣,采用下式確定
(14)
假設地表二維形變分量為d=(Dv,Dh)T,分別代表局部參考系的垂直和平面的形變分量,u=(Dv,Dh)T為局部參考系的雷達視線方向或者方位方向的單位形變量。根據SAR衛星的成像幾何關系,當u為雷達視線方向單位形變量時可表示為u=(cosθ,sinθ)T,合成孔徑雷達干涉(interferometric synthetic aperture radar, InSAR) 監測結果可表示為
(15)
為獲取監測區域地面形變速率和時序累計形變量,研究將計算結果導入ArcGIS軟件進行解譯分析。地面變化分為兩類:沉陷與上升,并結合礦山范圍總體計算地表形變時間序列情況。
研究采用Smart三維(three-dimensional,3D)建模軟件對無人機采集傾斜攝影數據進行處理,主要處理流程包括:數據預處理、空三加密、模型構建、三維模型修改及整理(圖6)。數據預處理針對采集數據情況,檢查數據質量,對照片進行勻色云光處理等。傾斜攝影獲取正視和斜視多個角度影像,空三加密較傳統攝影測量更復雜,多視角影像的空中三角測量以初始全球定位系統(global positioning system,GPS)數據為基礎,采取由粗到精的金字塔影像匹配策略,在每級影像上進行同名點匹配及光束發平差,得到同名點匹配結果,為確保平差精度,可建立連接點、控制點及慣性測量單元(inertial measurement unit ,IMU)輔助數據的多視角影像聯合解算。基于軟件進行空三加密為全自動解算過程,在導入照片和定位定姿系統(position and orientation system, POS)信息后,設置照片定位、地理參考及空三解算參數,開展傾斜影像空三解算。

圖6 實景三維模型技術流程圖
根據任務需求選擇分塊計算,自動選擇不同視角上的最佳像對模型,生成三維尺度的密集點云;基于三維重建算法,三維密集點云自動構建為不規則三角網(triangulated irregular network,TIN)模型,實現TIN模型進行平滑和優化;根據三角網TIN的空間位置信息,獲取最佳視角影像紋理,自動賦予模型紋理,輸出三維模型成果(圖7)。

圖7 中梁山片區傾斜攝影數據成果示例(osgb格式)
重慶市露天礦山具有面積小、分不廣等特征,目前礦山恢復治理過程中遙感動態監測還未有成熟的技術體系,研究結合重慶市礦山實際情況,在數據選擇、指標選取以及遙感監測技術方法實施等方面,提出了露天礦山恢復治理遙感動態監測技術方法體系,希望對未來全市礦山修復治理遙感監測提供基礎支撐。研究總結出以下幾種監測方式:
(1)礦山范圍的識別可通過高分辨率遙感影像,建立礦山遙感解譯樣本庫,對礦山范圍采用人機交互式遙感解譯提取礦山損毀土地范圍。
(2)礦山植被恢復狀況的遙感動態監測,可結合應用多時間序列的中分辨率或高分辨率遙感影像,利用植被在近紅外波段的光譜響應特征,構建植被指數NDVI,對礦山恢復治理植被動態變化進行分析。
(3)礦山地表形變監測,采用哨兵數據,通過干涉分析對礦山恢復治理過程中地表形變狀況進行動態分析。
本研究提出技術方法具有一定局限性,例如在植被遙感監測數據的選擇方面,雖然結合重慶礦山情況,沒有選擇中分辨率成像光譜儀(moderate-resolution imaging spectroradiometer, MODIS)等低分辨率遙感數據,但采用landsat數據,對于小型礦山植被監測也具有較大誤差。此外,未來隨著高分衛星數據的廣泛應用,長時間序列的礦山植被狀況監測可采用高分1號、高分6等數據,地表形變監測可采用高分3號數據。