潘建平,趙瑞淇,蔡卓言,袁雨馨,李鵬霞
(重慶交通大學智慧城市學院,重慶 400074)
Sentinel-1A數(shù)據(jù)從2014年10月對全球用戶開放以來[1],已經(jīng)被廣泛應用于地震與地質構造監(jiān)測[2-3]、火山與冰川監(jiān)測[4-5],以及其他形變監(jiān)測領域[6-8],現(xiàn)已成為合成孔徑雷達干涉測量(interferometric synthetic aperture radar,InSAR)領域內重要的C波段數(shù)據(jù)源(波長為5.6 cm)。但是其短波長的特性給其應用帶來了諸多限制,特別是在植被區(qū),C波段的雷達信號通常難以穿透植被的表層或頂層(如樹枝和樹葉),衛(wèi)星重訪周期內植被的理化特性變化,甚至是枝葉的微小運動,都會導致干涉失相干現(xiàn)象的發(fā)生,并將誤差引入InSAR地表形變監(jiān)測。
以往的研究表明,植被區(qū)的失相干程度可以在兩景單視復數(shù)影像(single look complex,SLC)發(fā)生干涉后由相干系數(shù)(coherence)進行評估[9],相干系數(shù)反映兩幅SLC影像的相似程度,可以衡量干涉條紋的質量與SLC影像中像素間的相位和振幅信息的變化[10]。但是干涉后評估必然會導致SAR數(shù)據(jù)存在選擇的盲區(qū),使得研究者不能選擇最有效的干涉對,從而降低InSAR對地表形變監(jiān)測的效率,增加了時間與經(jīng)濟成本。對于短波長的Sentinel-1A數(shù)據(jù),如果失相干能在干涉前得到精確評估,將十分有助于克服上述干涉后評估的缺點,從而提高監(jiān)測效率。
作為應用最廣泛的植被指數(shù),歸一化差分植被指數(shù)(normalized difference vegetation index,NDVI)與植被覆蓋度密切相關,可用于監(jiān)測植被的動態(tài)變化,對于地表覆被的可靠表征為建立其與InSAR失相干之間的定量模型提供了可能。
以往針對NDVI與Sentinel-1A相干系數(shù)的研究主要集中于定性說明兩者的關系[11-12],其定量關系尚不明確?;诖?本文選擇位于四川的一個植被覆蓋區(qū)作為研究區(qū),在同極化(VV)與交叉極化(VH)模式下利用免費光學影像Landsat 8計算得到的NDVI,分別建立其與Sentinle-1A相干系數(shù)之間的定量模型,并在鄰近區(qū)域對模型進行驗證。
本文的研究區(qū)位于四川省,其西部為川西高原,東部緊鄰成都平原,面積約為340 km2,如圖1紅色實線框所示。研究區(qū)內植被覆蓋茂盛,主要分布有林地和耕地(主要作物為水稻),并伴隨有少量村落分布,屬于InSAR地表形變監(jiān)測中的低相干地區(qū),在該區(qū)域建立模型有助于提高模型的普適性與應用價值;同時選擇位于研究區(qū)西南的鄰近區(qū)域(圖1黑色實線框所示)作為驗證區(qū),以驗證模型的可靠性,以及區(qū)內地形和植被覆蓋狀況與研究區(qū)類似。
圖1 研究區(qū)位置與影像覆蓋范圍
本文使用的光學影像為Landsat 8 OLI影像,成像時間為2019年8月11日,云量為0.89%,空間分辨率為30 m,覆蓋范圍如圖1中黑色虛線框所示。對光學影像進行大氣校正和輻射校正兩步處理,并按照研究區(qū)和驗證區(qū)矢量范圍進行裁剪,最后計算NDVI,公式為
(1)
式中,PNIR、PRED分別代表近紅外和紅光波段的反射率。NDVI的值介于[-1,1],隨著地表植被覆蓋度的增加而增大。NDVI為負值時代表地表覆蓋有水或雪,NDVI為0表示地表為巖石或裸土,NDVI為正值表示地表有植被覆蓋。研究區(qū)和驗證區(qū)的NDVI分布如圖2所示。
圖2 研究區(qū)和驗證區(qū)NDVI分布
Sentinel-1A影像分別為成像時間為2019年7月3日和8月20日的干涉寬幅模式下的SLC影像,時間基線為48 d,分別作為干涉對的主影像和輔影像,其覆蓋范圍如圖1中紅色虛線框所示。SLC影像的預處理包括3部分: ①使用雙過差分對主、輔影像進行差分干涉,并分別使用精密軌道數(shù)據(jù)和SRTM DEM(30 m)消除軌道誤差和模擬地形相位; ②對干涉對進行多視處理以消除斑點噪聲并使SAR影像與Landsat 8影像保持相同的分辨率;③使用文獻[13]提出的公式消除空間失相干的影響,并計算研究區(qū)和驗證區(qū)的相干系數(shù),公式為
(2)
式中,γ代表相干系數(shù);E[·]代表數(shù)學期望;μ*代表SLC影像的共軛復數(shù)。相干系數(shù)的值介于[0,1],當主輔影像完全相干時,相干系數(shù)的值為1;當兩者完全失相干時(如干涉圖完全被噪聲污染),其值為0。
最后使用SRTM DEM(30 m)對NDVI分布圖和相干系數(shù)圖進行配準,配準后完成數(shù)據(jù)預處理,數(shù)據(jù)預處理流程如圖3所示。
圖3 數(shù)據(jù)預處理流程
文獻[14]提出了對相干系數(shù)產生影響的各個因素,主要有空間失相干、時間失相干及熱噪聲失相干。其中,熱噪聲失相干可以忽略不計;空間失相干由于研究區(qū)較為平坦的地形(平均坡度小于15°)并不明顯,且已通過文獻[13]提出的公式消除;通過建立時間基線與相干系數(shù)的模型削弱時間失相干對相干系數(shù)的影響。
基于多主影像策略,利用覆蓋研究區(qū)的Sentinel-1A影像生成了若干干涉對,干涉對的時間基線為12~696 d,時間基線與相干系數(shù)的關系如圖4所示。
圖4 時間基線與相干系數(shù)的關系
圖4(a)顯示,VV極化下相干系數(shù)在時間基線≤216 d時隨著時間基線的增加呈現(xiàn)指數(shù)衰減,其擬合曲線為y=A1·exp(-x/t1),其中A1=0.743,t1=206 d。在216 d之后沒有明顯的衰減趨勢,相干系數(shù)在0.370左右浮動。圖4(b)顯示,VH極化下相干系數(shù)在時間基線小于168 d時呈現(xiàn)指數(shù)衰減,其擬合曲線為y=A2·exp(-x/t2),其中A2=0.560,t2=222 d。在168 d之后沒有明顯的衰減趨勢,在0.325左右浮動。本文所使用的干涉對時間基線為48 d,因此有必要考慮時間失相干對相干系數(shù)指數(shù)衰減的作用。
試驗表明,將配準后的NDVI分布圖與相干系數(shù)圖的所有像素進行擬合分析會引入非常大的全局誤差,從而無法建立起NDVI和相干系數(shù)之間的可靠模型,因此本文考慮使用文獻[15]中提出的窗口采樣方法。該方法包括3個步驟。
(1)對NDVI分布圖設置一個采樣窗口,同時對相干系數(shù)圖設置一個相同大小的采樣窗口,使用采樣窗口對兩張圖像同時進行采樣,直至窗口遍歷整張圖像。
(2)對兩個窗口中的NDVI像素值和相干系數(shù)像素值計算皮爾遜相關系數(shù)[16],公式為
(3)
(3)設置一個閾值,若兩個采樣窗口中的像素值皮爾遜相關系數(shù)大于該閾值,則將窗口中的像素值保留,參與后續(xù)分析;否則,將窗口中的像素值舍棄。
該方法被用于建立NDVI與L波段ALOS數(shù)據(jù)相干系數(shù)之間的關系,但是ALOS數(shù)據(jù)由于其長波段特性失相干現(xiàn)象較弱,NDVI與相干系數(shù)之間的負相關關系并不明顯。初步試驗表明,Sentinel-1A數(shù)據(jù)的相干系數(shù)與NDVI之間呈現(xiàn)較強的負線性相關關系,直接使用該方法將會忽略掉大量像素值,導致模型精度大幅降低。因此,本文對該方法的相似性測度進行了改進,在步驟(3)中計算皮爾遜相關系數(shù)的絕對值,以顧及采樣窗口中NDVI與Sentinel-1A相干系數(shù)大量呈負線性相關的像素值。最佳采樣窗口大小和預設閾值大小的設定需要通過試驗確定,改進后的窗口采樣方法如圖5所示。
圖5 改進相似性測度的窗口采樣方法
圖5上下兩張圖像分別為研究區(qū)NDVI分布圖和相干系數(shù)圖,T為皮爾遜相關系數(shù)的絕對值,T0為預設的閾值。大量試驗表明,對于VV極化,最佳的采樣窗口大小為5×5,閾值大小為0.7,VH極化下的最佳采樣窗口和閾值大小為9×9和0.6。將滿足閾值條件的保留像素值進行分析,得到了VV和VH極化下NDVI和相干系數(shù)的關系,如圖6所示。
圖6 NDVI和相干系數(shù)的關系
圖6中藍色星星為真實像素值,兩條直線分別為VV和VH極化下的擬合曲線。由于時間失相干對相干系數(shù)的衰減作用,因此需要在模型中加入時間失相干指數(shù)衰減因子,在95%的置信區(qū)間下采用最小二乘擬合得到了VV和VH極化模式下NDVI和Sentinel-1A相干系數(shù)的定量模型,分別為
(4)
(5)
基于式(4)和式(5),利用NDVI分布圖模擬得到了VV和VH極化下的研究區(qū)相干系數(shù)圖,與研究區(qū)真實相干系數(shù)的對照如圖7所示。
圖7 研究區(qū)相干系數(shù)
圖7(a)和圖7(c)分別為研究區(qū)VV和VH極化下的真實相干系數(shù),圖7(b)和圖7(d)分別為VV和VH極化下的模擬相干系數(shù)。利用真實相干系數(shù)減去模擬相干系數(shù)得到了圖7(e)和圖7(f)所示的誤差分布圖。由誤差分布圖可以看出,VV極化下誤差分布除研究區(qū)內河流分布區(qū)域外其他區(qū)域不明顯,VH極化下誤差分布除了河流區(qū)域外主要集中在圖7(f)中的村落分布區(qū)域。對誤差分布進行統(tǒng)計,如圖8所示。研究區(qū)VV極化下誤差均值為-0.037,標準差為0.205,大部分誤差集中在-0.3~0.3之間,誤差整體基本服從正態(tài)分布。研究區(qū)VH極化下誤差均值為0.045,標準差為0.200,大部分誤差同樣集中分布在-0.3~0.3之間,雖然誤差整體高于VV極化,但是仍然基本服從正態(tài)分布。
圖8 研究區(qū)誤差統(tǒng)計
為了驗證模型的可靠性,在驗證區(qū)進行了相同的試驗,數(shù)據(jù)處理方式與研究區(qū)保持一致,驗證區(qū)的相干系數(shù)如圖9所示。圖9(a)和圖9(c)分別為驗證區(qū)VV和VH極化下的真實相干系數(shù),圖9(b)和圖9(d)分別為VV和VH極化下的模擬相干系數(shù),利用真實相干系數(shù)減去模擬相干系數(shù)得到了圖9(e)和圖9(f)所示的誤差分布。由于驗證區(qū)分布有兩條明顯的河流,因此由誤差分布圖可以看出,驗證區(qū)的誤差主要分布在河流分布區(qū)域及區(qū)域南部的村落分布區(qū)。對驗證區(qū)誤差分布進行統(tǒng)計,如圖10所示。驗證區(qū)VV極化下的誤差均值為-0.067,標準差為0.256,VH極化下的誤差均值為-0.065,標準差為0.230。兩種極化下的誤差均集中分布在-0.3~0.3之間,且服從正態(tài)分布。綜合圖7(e)、圖7(f)、圖9(e)和圖9(f)的差值及圖8和圖10的誤差分布情況來看,利用式(4)和式(5)模擬得到的相干系數(shù)沒有隨機誤差分布,模型可靠。
圖9 驗證區(qū)相干系數(shù)
圖10 驗證區(qū)誤差統(tǒng)計
從研究區(qū)和驗證區(qū)的角度來看,兩種極化模式下驗證區(qū)模擬相干系數(shù)的誤差分布更離散且整體大于研究區(qū),主要原因有兩點:①盡管驗證區(qū)的地形與研究區(qū)相似,但仍存在微小差異,雖已去除空間失相干的影響,但是微小地形差異引起的空間失相干會將誤差引入模型中;②相較于研究區(qū),驗證區(qū)分布有兩條明顯的河流,河流區(qū)域的誤差成為驗證區(qū)的主要誤差來源。
從極化方式來看,無論研究區(qū)還是驗證區(qū),VV極化下的相干系數(shù)整體均高于VH極化下的相干系數(shù),這是由于研究區(qū)和驗證區(qū)的地表覆被主要為耕地(主要作物為水稻),其次是林地,已有研究表明水稻的理化參數(shù)與反射特性與VH極化的關系比VV極化更強[16],導致其對VH極化的雷達波更敏感,相同時間基線下的失相干相較于VV極化更嚴重。同時研究區(qū)和驗證區(qū)VV極化下的整體誤差均低于VH極化下的整體誤差,這是由于盡管經(jīng)過大量試驗得到了兩種極化方式下的最佳采樣窗口與閾值,但是VH極化下的像素值相較于VV極化分布更離散(如圖6所示),整體擬合精度低于VV極化,從而導致VH極化下模型精度低于VV極化下模型精度。
綜合上述模型的誤差源與缺點,在未來的研究中還需考慮如何完全去除空間失相干對模型的影響,另外水系不可避免地會出現(xiàn)在植被覆蓋區(qū),因此應考慮如何消除水系對模型的影響,如建立光學遙感水體指數(shù)與相干系數(shù)的定量模型并融入已建立的模型),從而提高模型的普適性與可靠性。其他諸如研究區(qū)的天氣狀況、地表植被類型及影像不同的成像視角(升軌與降軌)等因素也應當在模型中加以考慮。除了本文建立的NDVI與C波段Sentinel-1A相干系數(shù)之間的定量模型之外,NDVI與同極化L波段的ALOS-1/PALSAR-1相干系數(shù)之間的定量模型已經(jīng)建立[15],因此在未來的工作中需要建立NDVI與交叉極化的L波段數(shù)據(jù)的相干系數(shù)及與X波段數(shù)據(jù)的相干系數(shù)之間的定量模型。
本文建立了兩個基于Landsat 8影像計算得到的NDVI與Sentinel-1A不同極化方式下干涉相干系數(shù)之間的二階線性模型,從而得到了InSAR失相干與地表植被覆蓋之間的定量關系。相干性隨著植被覆蓋度的增加而呈線性減小,且線性減小的趨勢因極化方式的不同而有所區(qū)別,同時這兩個模型在驗證區(qū)得到了驗證,表明了模型的可靠性?;诒疚慕⒌哪P?可以利用免費光學影像計算得到的NDVI在進行InSAR處理之前(即SLC影像干涉之前)對研究區(qū)的InSAR失相干現(xiàn)象做出定量評估;同時在干涉前對SAR數(shù)據(jù)選擇提供指導,從而提升研究者的效率并減少時間與經(jīng)濟成本。