李彬, 孫小龍*, 蔣姍姍, 盧士慶, 李云鵬, 李靜

(1.內(nèi)蒙古自治區(qū)生態(tài)與農(nóng)業(yè)氣象中心, 呼和浩特 010051; 2.高分辨率對地觀測系統(tǒng)內(nèi)蒙古數(shù)據(jù)與應(yīng)用中心, 呼和浩特 010051;3.通遼市氣象局, 通遼 028000; 4.滿都拉氣象站, 包頭 014000)

沙塵作為吸收性氣溶膠對太陽輻射有較強的吸收作用,加熱大氣改變大氣穩(wěn)定度及其與云微物理的相互作用從而減少云量抑制降水。同時,作為環(huán)境污染物,對人類健康有一定危害。而且對社會經(jīng)濟(jì)發(fā)展,生態(tài)文明以及水循環(huán)等一系列活動有著直接或間接的影響。沙塵天氣一般在春季爆發(fā),對中國北方大部地區(qū)的生產(chǎn)和生活有較大影響[1-2]。

地基觀測資料的限制使得沙塵氣溶膠的研究還有很大的不確定性,而且沙塵爆發(fā)一般具有較大的影響范圍,一般發(fā)生在自然條件惡劣的地區(qū),而這些地區(qū)地基觀測站點少且分散,給監(jiān)測、預(yù)報和研究沙塵暴帶來很大困難。而衛(wèi)星遙感監(jiān)測范圍廣、波段多的優(yōu)勢可以彌補這些不足,因此,針對沙塵遙感監(jiān)測人們開展了許多研究,發(fā)展了許多算法。

目前,科學(xué)家在利用衛(wèi)星資料識別沙塵爆發(fā)方面主要通過多通道組合,主要是可見光、近紅外通道及熱紅外通道。Ackerman[3]較早提出分裂窗插值識別沙塵,利用紅外亮溫(brightness temperature,BT)監(jiān)測沙塵,為以后的沙塵定量識別提供了思路。Xu等[4]利用中分辨率成像光譜儀(moderate-resolution imaging spectroradiometer,MODIS)第3和第7波段構(gòu)建歸一化差值沙塵指數(shù) (normalized difference dust index,NDDI)來判識沙塵。詹奕哲等[5]結(jié)合成熟的紅外分裂窗差值監(jiān)測沙塵及衛(wèi)星導(dǎo)風(fēng)技術(shù)。Romano等[6]利用基于可見光和紅外波段的多通道方法對二代氣象衛(wèi)星MSG(meteosat second generation)所搭載的旋轉(zhuǎn)增強型可見光及紅外線成像輻射計(spinning enhanced visible and infra-red imager,SEVIRI)的數(shù)據(jù)進(jìn)行了沙塵識別。羅敬寧等[7]對多通道沙塵光譜特征進(jìn)行分析,利用風(fēng)云三號衛(wèi)星的近紅外1.6 μm、中紅外3.7 μm以及熱紅外分裂窗進(jìn)行了一次全球沙塵的實例監(jiān)測。Di等[8]提出一種包含可見光、中紅外、熱紅外波段的增強粉塵指數(shù)(enhanced dust index,EDI)來進(jìn)行INSAT-3A/D衛(wèi)星的沙塵監(jiān)測。

研究表明紅外分裂窗在沙塵識別上具有較大優(yōu)勢,特別是在高反射地表和夜間,因此紅外亮溫差(brightness temperature difference,BTD)也經(jīng)常被用于沙塵的監(jiān)測。延昊等[9]利用熱紅外亮溫差(T11 μm～T12 μm)針對NOAA/AVHRR資料進(jìn)行沙塵監(jiān)測。Zhang等[10]利用BTD進(jìn)行了MODIS數(shù)據(jù)的沙塵監(jiān)測,并利用輻射傳輸模式建立了不同光學(xué)厚度、BTD、BT11 μm與有效粒子半徑的查找表。Chaboureau等[11]利用BTD對沙塵和卷云的區(qū)域氣象預(yù)測模型進(jìn)行了評估。Schepanski等[12]采用8.7、10.8 、12 μm三通道亮溫差對撒哈拉沙漠地區(qū)的沙塵源區(qū)進(jìn)行了分析。Huang等[13]也利用雷達(dá)資料提出了微波極化亮溫差 (microwave polarization temperature difference,MPTD) 微波極化指數(shù) (microwave polarization index,MPI)來識別云下沙塵,取得了一定效果。Kluser等[14]結(jié)合BTD與時間序列方法進(jìn)行沙塵監(jiān)測,Sang等[15]在BTD的基礎(chǔ)上引入了亮溫比(temperature difference ratio,BTR)來消除地表溫度變化造成的影響。也有研究利用地物比輻射率特征(比輻射率是指物體的輻射出射度與同溫度下黑體輻射出射度的比值),結(jié)合熱輻射波段特征,討論了沙塵暴發(fā)時比輻射率在熱紅外波段處的變化來識別沙塵[16]。比輻射率在云和沙塵混合的區(qū)域識別效果較差,因為水汽云在11 μm的散射能力強于12 μm,而吸收能力弱于12 μm,表現(xiàn)出與沙塵相似的比輻射率特征,特別是弱沙塵區(qū)域[5]。Yang等[17]利用NPP/VIIRS的3.7、11、12 μm的BTD組合進(jìn)行沙塵判識。

可見光波段難與沙漠戈壁等亮地表區(qū)分;熱紅外波段對強沙塵的判識效果較好,但對于弱沙塵則效果較差;微波傳感器目前多搭載于極軌衛(wèi)星,數(shù)據(jù)存在時間和空間分辨率都較低的局限,無法對沙塵發(fā)生發(fā)展進(jìn)行實時有效的監(jiān)測預(yù)警。因此如何利用高時空分辨率數(shù)據(jù)進(jìn)行更加有效的沙塵識別是當(dāng)前主要的研究方向。

Himawari-8是世界上第一顆可拍攝真彩色圖像的靜止氣象衛(wèi)星,觀測頻率提高到每10 min一次,對沙塵天氣過程及動向的持續(xù)觀測性能也顯著提升,從真彩色圖像上還能解譯其他圖像中與云層混雜的沙塵、弱浮塵。因此,研究者針對Himawari-8沙塵監(jiān)測也開展了許多研究,而利用亮溫差BTD是目前應(yīng)用最廣泛的監(jiān)測方法。She等[18]利用Himawari-8的3.9、8.6、11、12 μm的BTD組合進(jìn)行沙塵判識,并針對亮、暗地表進(jìn)行了閾值區(qū)分。邱昀等[19]利用Himawari-8熱紅外BTD判識沙塵并結(jié)合地面氣象數(shù)據(jù)及激光雷達(dá)進(jìn)行了三維動態(tài)觀測。李彬等[20]結(jié)合熱紅外BTD及可見光灰度熵方法對Himawari-8純沙塵及混合性沙塵進(jìn)行了判識。劉方偉等[21]通過對典型地物及沙塵光譜特性和多波段亮溫差分析,提出組合閾值進(jìn)行沙塵監(jiān)測。以上方法閾值的確定多以強度較大的沙塵暴進(jìn)行經(jīng)驗統(tǒng)計得到,且在不同場景(如晝夜)下的不確定性較大。Himawari-8的11 μm與12 μm的BTD隨時空變化有明顯差異,單一閾值容易產(chǎn)生漏判。張海香等[22]基于Himawari-8的11 μm和12 μm的BTD背景場,提出基于動態(tài)閾值的沙塵判識方法。Legrand等[23]也曾利用紅外亮溫粉塵差異指數(shù)(infrared dust difference index,IDDI)來識別沙塵。但是此方法本身也存在著一些不足,對強度大持續(xù)時間長的沙塵效果較差,動態(tài)參考閾值的確定也是制約因素。

一次持續(xù)時間長,影響范圍大的沙塵過程一般都伴隨著不同強度的沙塵,而局地沙塵過程通常強度較弱,而沙塵起始階段的弱沙塵對沙源識別具有重要意義;特別是在干旱半干旱地表判識較為困難,且云沙混合區(qū)及沙漠戈壁易發(fā)生誤判[24-25]。因此,現(xiàn)使用Himawari-8衛(wèi)星數(shù)據(jù)提出一種新的基于動態(tài)亮溫差背景場的沙塵覆蓋度指數(shù)對不同地表條件及強度下的沙塵進(jìn)行判識,并據(jù)此進(jìn)行了沙塵強度估算。

1 數(shù)據(jù)介紹

Himawari-8是日本氣象廳的第二代靜止氣象衛(wèi)星,有16個通道,分布于可見光-近紅外-熱紅外,搭載高級成像儀(advanced himawari imager,AHI)能掃描區(qū)域包括:Full Disk(全球范圍)、Japan Area(兩個日本地區(qū)的范圍)、特定區(qū)域(兩個范圍)、地標(biāo)性區(qū)域(兩個范圍)。用于拍攝全球范圍和日本地區(qū)的掃描范圍是固定的,其他的兩個特定區(qū)域和地標(biāo)性區(qū)域是可以靈活調(diào)整的,不論是哪種拍攝模式,時間分辨率最少都能達(dá)到10 min,分辨率0.5～2 km。

本文中使用了2017—2019年春夏季發(fā)生于中國內(nèi)蒙古及其周邊地區(qū)的幾次典型的沙塵天氣過程,所用數(shù)據(jù)的日期和時次如表1所示,結(jié)合目視解譯與地面觀測對方法進(jìn)行了分析驗證。研究中所使用影像范圍為北緯41°～54°、東經(jīng)97°～127°和北緯41°～54°、東經(jīng)110°～127°。為了分析方法對強度較弱沙塵過程的監(jiān)測效果,選取了三處典型干旱半干旱區(qū)域弱沙塵天氣事件,即2019年4月15日、2019年3月27日、2019年4月20日,分布標(biāo)記為S1、S2和S3。2017年5月3日的沙塵過程波及的范圍較大,有不同強度分布,因此利用這一天的數(shù)據(jù)對本文算法在大范圍進(jìn)行了分析和驗證。地面觀測數(shù)據(jù)來自全國綜合氣象信息共享平臺和地方氣象臺站觀測的天氣現(xiàn)象。分析過程中,還使用了2020年全球30 m地表覆蓋精細(xì)分類產(chǎn)品V1.0,MODIS的植被指數(shù)NVDI產(chǎn)品(MYD13A3)、MODIS的氣溶膠光學(xué)厚度產(chǎn)品(MCD19A2)及SRTM的DEM數(shù)據(jù),分辨率均統(tǒng)一為1 km。

表1 實驗所用Himawari-8衛(wèi)星影像列表Table 1 List of dates of Himawari-8/AHI date used for training

2 判識方法

2.1 云、水體、冰雪檢測

由于云、水體、冰雪會對判識所需底圖及待判識圖像歸一化等處理產(chǎn)生影響,因此首先需要對數(shù)據(jù)中的上述目標(biāo)物進(jìn)行掩膜處理。具體方法如參考了葛邦宇等[26]針對Himawari-8的算法。

2.2 熱紅外亮溫差底圖

對于長波紅外輻射而言,它們在空中組成類似云的物質(zhì)表面,足以改變長波發(fā)射輻射的分布。在波長11 μm和12 μm的熱紅外波段,水汽的吸收存在較大差異。沙塵氣溶膠折射率虛部的差異會影響沙塵的大氣輻射,導(dǎo)致亮溫差ΔT<0[9]。近紅外和中紅外波段,對于強沙塵的識別效果較好,對于弱浮塵,受到地表反射輻射影響較大,識別效果較差。熱紅外分裂窗利用兩個通道的差異來識別,這種對強弱沙塵都存在。BTD方法還可以有效剔除大范圍連續(xù)云區(qū)。因此,本文中使用BTD值構(gòu)建亮溫差背景場底圖,其表達(dá)式為

BTD=T11 μm-T12 μm

(1)

Himawari-8作為一顆靜止氣象衛(wèi)星具有全天時觀測能力,而地表溫度具有日變化特征,因此不同于極軌衛(wèi)星數(shù)據(jù),靜止衛(wèi)星觀測亮溫和亮溫差也具有一定日變化特征。本文中利用2019年4月2日08:00—17:00的晴空數(shù)據(jù)生成了各時次的亮溫差產(chǎn)品,并依據(jù)中國科學(xué)院空天信息創(chuàng)新研究院發(fā)布的2020年全球30 m地表覆蓋精細(xì)分類產(chǎn)品提取了城市、農(nóng)田、林地和沙地(包含沙漠和戈壁)四種典型地物在各時次的亮溫差值(10個代表性像元的平均值)得到其日變化規(guī)律,如圖1所示。

圖1 不同地物BTD值日變化Fig.1 Diurnal variation of BTD for different land covers

可以看出,四種地物類型均表現(xiàn)出先升后降的特征,一般都在中午前后到達(dá)最大值,在晨昏時間為最小值。植被覆蓋較高的林地的日波段最小,植被稀疏的沙地日波段最為劇烈,日較差近3 K。 為了研究BTD在不同季節(jié)之間的變化,我們選取了能表征地表狀態(tài)變化的植被覆蓋度指數(shù)(fractional vegetation cover,FVC)來進(jìn)行分析比較。從圖2可以看到BTD的分布與植被覆蓋度之間沒有特別明顯的規(guī)律性,其中還可能受到海拔和緯度的影響,但是在植被覆蓋度低的區(qū)域,BTD的高值也有較為集中的分布,而在植被覆蓋度高的區(qū)域也有BTD低值的集中分布。

圖2 植被覆蓋度和海拔高度與BTD值分布散點圖Fig.2 Scatter plot of BTD with fractional vegetation cover and altitude

She等[18]在研究中發(fā)現(xiàn)當(dāng)海拔低于3 000 m時,BTD值通常大于2.5 K。本文中選取了研究區(qū)內(nèi)海拔2 500 m以下的樣本進(jìn)行統(tǒng)計,結(jié)果基本一致;在1～2.5 K也有分布,總體隨著海拔升高BTD值呈現(xiàn)下降趨勢,受植被狀態(tài)、緯度等影響在高海拔 區(qū)也有BTD高值的分布。基于以上分析,BTD值分布在不同區(qū)域、不同時段具有較大的不確定性,因此背景場的構(gòu)建需基于動態(tài)。

構(gòu)建背景場?；谏鲜鰧τ贐TD在不同時間,不同植被覆蓋度及海拔等分布規(guī)律的討論,需要針對判識區(qū)構(gòu)建動態(tài)背景場。構(gòu)建晴空亮溫差背景場是整個沙塵監(jiān)測的重要一步,背景場的純凈程度將決定沙塵判識的準(zhǔn)確性,特別是對一些強度較弱的沙塵。

首先對云、水體、冰雪等進(jìn)行掩膜處理。掩膜處理后選擇同一時刻的多天數(shù)據(jù),保證判識區(qū)域內(nèi)的每一個像元為晴空背景,并把多天的數(shù)據(jù)合成為一個最小值矩陣,設(shè)為背景場矩陣。表達(dá)式為

(2)

式(2)中:n為構(gòu)建動態(tài)背景場的參考數(shù)據(jù)的個數(shù),本文算法選取待判識日期前7～10 d的數(shù)據(jù),保證期間植被狀態(tài)改變不大,其中有沙塵天氣的要剔除;i、j為像元的行列號;ΔT為亮溫差值。如圖3所示為2019年4月下午16:00的亮溫差背景場合成圖?？梢钥闯鲈诓煌貐^(qū)表現(xiàn)出較大差異,在高海拔的林區(qū)BTD值較低,在低海拔的城市區(qū)和沙漠戈壁等地區(qū)BTD則相對較高。

圖3 2019年4月16:00亮溫差背景圖Fig.3 Background of BTD at 16:00 on April,2019

2.3 沙塵覆蓋度指數(shù)(FDCI)

目前大部分沙塵監(jiān)測都是基于亮溫、亮溫差等閾值進(jìn)行的,而由于衛(wèi)星觀測的這些參量在不同時間、地表狀態(tài)、海拔等具有較大差異,因此使用固定閾值存在局限,使用動態(tài)閾值也具有較大波動性。因此提出一種基于動態(tài)背景場的沙塵覆蓋度指數(shù)來進(jìn)行沙塵判識,該方法既能通過構(gòu)造動態(tài)背景彌補亮溫等變化產(chǎn)生的固定閾值局限,也能通過指數(shù)化的方法將判斷閾值進(jìn)行固定化,避免動態(tài)閾值選取的不確定性。

基于11 μm和12 μm的熱紅外波段亮溫差建立動態(tài)背景場,并對該背景場進(jìn)行歸一化。由于BTD的值域差別較大,且存在負(fù)值,歸一化之后可以把數(shù)值統(tǒng)一為0～1,這樣便于進(jìn)行指數(shù)化的計算。歸一化公式為

(3)

φnorm0=1-φnorm

(4)

式中:φ和φnorm別為歸一化前后的變量值;φmax、φmin和φmean分別為變量BTD的最大值、最小值和平均值。由于BTD負(fù)值越小的沙塵強度越大,為了讓歸一化值與之對應(yīng),用1減去φnorm得到最后的φnorm0值[27]。圖4所示為一次沙塵過程中BTD<0的像元分布情況和其歸一化φnorm0的分布情況,歸一化后二者呈現(xiàn)出鏡像分布。

圖4 一次沙塵過程熱紅外亮溫差及其歸一化值分布Fig.4 Distribution for BTD of heavy dust and its normalized value

在遙感生態(tài)評估業(yè)務(wù)經(jīng)常使用植被覆蓋度FVC來表征植被的生長狀態(tài),Wang等[28]也提出積雪覆蓋指數(shù)(fractional snow cover,FSC)來進(jìn)行積雪的判識。沙塵覆蓋度指數(shù)也是基于以上思想,利用該指數(shù)來反映像元沙塵覆蓋的程度,并據(jù)此判斷是否為沙塵。 將待判識數(shù)據(jù)的BTD與背景場的BTD分別進(jìn)行歸一化后,再計算沙塵覆蓋度(fractional dust cover index,FDCI)指數(shù)值,公式為

(5)

式(5)中:φnorm0,det為待判識像元的歸一化值;φnorm0,back為像元背景場歸一化值;φnorm0,max為沙塵歸一化值統(tǒng)計分布中的最大值,這里取0.75。

經(jīng)過實驗統(tǒng)計發(fā)現(xiàn),將該閾值設(shè)定為大于0.2能夠覆蓋絕大部分沙塵區(qū)。圖5、圖6所示為從2017年5月3日、2019年3月27日、4月15日、4月20日、3月26日和2018年4月29日的六次沙塵過程中結(jié)合地面觀測選取不同強度的沙塵像元FDCI值的分布情況。其中有6.8%的值小于0.2,59.3%在0.2～0.6之間,33.9%大于0.6。不同強度沙塵FDCI和BTD的分布情況表明,強度較弱沙塵的亮溫差值有約90%大于0,而FDCI值則絕大部分大于0.2。

圖5 7 166個弱及中等強度沙塵像元FDCI和BTD分布Fig.5 Distribution of FDCI and BTD of 7 166 thin and moderate dust pixels

圖6 9 980個強沙塵像元的FDCI值和BTD值分布Fig.6 Distribution of FDCI and BTD of 9 980 heavy dust pixels

為了分析FDCI對不同強度沙塵的區(qū)分效果,分別選取了100個弱沙塵和300個強沙塵(主要是目視解譯的純像元)像元的FDCI進(jìn)行統(tǒng)計,得到如圖7所示的箱線圖?？梢钥闯鰪娚硥m的FDCI區(qū)間為0.4～0.8,集中分布于0.55～0.65,算法對此類沙塵的區(qū)分作用明顯。弱沙塵的FDCI區(qū)間為0.1～0.4,集中分布于0.2～0.3,與強沙塵的FDCI分布具有一定連續(xù)性。因此FDCI指數(shù)對于不同強度沙塵區(qū)分效果明顯,特別是對弱沙塵。

圖7 弱沙塵與強沙塵的FDCI分布箱線圖Fig.7 Boxplot of FDCI for thin dust and heavy dust

2.4 孤立區(qū)的剔除與填充

由于沙塵天氣過程通常都會呈現(xiàn)出連續(xù)分布特征,極少出現(xiàn)孤立分布,因此這些零散孤立區(qū)域很大可能是誤判,應(yīng)予濾除。如圖8所示,要進(jìn)一步濾除這些誤判的孤立像元,只需將不連續(xù)的細(xì)碎區(qū)域濾除,保留具有連續(xù)特征的沙塵區(qū)即可。

圖8 灰度熵處理結(jié)果Fig.8 The results before and after using entropy of brightness

本文中采用基于圖像灰度熵的團(tuán)塊目標(biāo)檢測方法。該方法假設(shè)n×n的窗口內(nèi)像素灰度分別為a1,a2,…,an,則該窗口內(nèi)的灰度熵H定義為

(6)

(7)

灰度熵(0≤H≤1)反映了窗口內(nèi)像素灰度的差異程度,值越大表明窗口內(nèi)能量的空間分布越均一,當(dāng)窗口內(nèi)所有像素灰度值相同時,圖像灰度熵為1;反之,表明灰度差異大。

首先將經(jīng)過RDI(reflectivity difference index)判識之后的結(jié)果轉(zhuǎn)化為二值圖像,并設(shè)定3×3像元的窗口進(jìn)行灰度熵檢測[20]。這樣窗口內(nèi)的灰度熵值就只有0和 1兩個值,將值為0的窗口判斷為非沙塵區(qū),反之為沙塵。

另外,在一片連續(xù)的沙塵區(qū)中也會有受到云干擾而被漏判的像元,對于這類像元應(yīng)該進(jìn)行補充判識。本文中對待填充區(qū)域(云覆蓋區(qū))首先進(jìn)行邊緣提取,然后逐一判斷區(qū)域邊緣像元是否為已判識沙塵像元,如果全部為沙塵像元則說明該區(qū)域為連續(xù)沙塵區(qū)中的漏判部分并予以填充。

2.5 基于FDCI的沙塵強度指數(shù)(EDII)

沙塵強度作為沙塵氣象監(jiān)測的重要指標(biāo)也是遙感定量反演的主要對象,一般是在沙塵判識結(jié)果的基礎(chǔ)上進(jìn)行參數(shù)估算。不同研究者從需求角度出發(fā),有不同的強度等級劃分標(biāo)準(zhǔn)。2006年國家標(biāo)準(zhǔn)《沙塵暴天氣等級》(GB/T 20480—2006),規(guī)定了沙塵天氣過程等級,將沙塵天氣強度等級劃分為5個等級:浮塵、揚沙、沙塵暴、強沙塵暴、特強沙塵暴。 遙感反演沙塵強度主要依據(jù)反射率、亮溫和光學(xué)厚度幾個指標(biāo),羅敬寧等[29]利用1.6 μm波段構(gòu)建了可比沙塵強度指數(shù),隨后又利用11 μm和12 μm亮溫差和1.6 μm波段反射率建立了FY-3的沙塵強度指數(shù)。Di等[8]利用INSAT衛(wèi)星數(shù)據(jù)提出一種增強型沙塵指數(shù)(enhanced dust intensity index,EDII),該指數(shù)綜合了歸一化的反射率、亮溫及光學(xué)厚度。公式為

(8)

式(8)中:A=1,B=10,C=0.1;R1.6表示1.6 μm波段反射率,R0.47表示0.47 μm波段反射率。AOD由延昊等[9]利用SBDART輻射傳輸模式模擬的沙塵BTD(11～12 μm)與0.55 μm光學(xué)厚度之間的關(guān)系計算得到。公式為

(9)

但是上述方法計算沙塵強度主要適用于強度較大的沙塵天氣,而弱沙塵通常對于地表有一定透明度,地表信息對沙塵的光譜信息存在干擾,導(dǎo)致算法對于強度弱的沙塵適用性較差,計算強度一般偏大。

因此對于強度較大的沙塵使用式(8)進(jìn)行計算,對于強度較弱的沙塵(0.2

EDII=1.67e-0.25VIS+0.47

(10)

式(10)中:VIS為能見度,km。能見度計算根據(jù)余娟等[30]高斯曲線優(yōu)化的能見度與氣溶膠光學(xué)厚度轉(zhuǎn)換模型。公式為

(11)

式(11)中:τα為氣溶膠光學(xué)厚度;M為月份參數(shù),當(dāng)觀測月份大于6時,M取值為月份減6,反之加6。a、b、c為模型經(jīng)驗系數(shù)。

3 結(jié)果與分析

本節(jié)主要對前文所述實驗數(shù)據(jù)和區(qū)域,首先針對弱沙塵利用FDCI進(jìn)行識別,同時也在相同區(qū)域使用多波段分裂窗亮溫差[18]方法進(jìn)行判識,并將結(jié)果進(jìn)行比較分析。圖9～圖11分別為區(qū)域S1～S3對應(yīng)的真彩色影像、FDCI算法判識結(jié)果、多波段BTD所得結(jié)果。最后在大范圍區(qū)域內(nèi)進(jìn)行判識及沙塵強度估算,并與觀測數(shù)據(jù)進(jìn)行對比驗證。

圖9 區(qū)域S1真彩色圖像與兩種方法的判識結(jié)果Fig.9 True color image of area S1 and detective results by two methods

區(qū)域S1為一個典型的弱沙塵與強沙塵都存在的過程,從圖9(a)的真彩色圖像可以明顯看出上方的沙塵強度較弱,能一定程度得辨識出地表,這樣強度的沙塵結(jié)合沙漠地表使得判識難度大大增加,導(dǎo)致BTD方法不能對此類情況進(jìn)行有效識別。如圖9(c)所示,BTD所得結(jié)果沒有判識出該區(qū)域的沙塵。如圖9(b)所示,FDCI方法則準(zhǔn)確地識別出了區(qū)域內(nèi)這段強度較弱的沙塵,且兩者對下方強度較大沙塵的判識結(jié)果基本一致。表2中包含了氣象站的經(jīng)緯度和當(dāng)時次所觀測的天氣現(xiàn)象,可以看出觀測與判識結(jié)果一致。

表2 不同時間地面氣象觀測與算法識別結(jié)果Table 2 Observation on the ground and detective results by FDCI of various time points

由圖10(a)的真彩色圖像可以看到區(qū)域S2為一個典型的沙塵與積云混合狀態(tài),碎積云分布細(xì)碎而分散,一般由對流生成。沙塵區(qū)的下部及邊緣處的沙塵較薄,因此BTD綜合判識方法未能對此類情況進(jìn)行有效識別,如圖10(c)所示。FDCI方法則有效地識別出了區(qū)域內(nèi)這些強度較弱的沙塵,圖10(b)所示。

圖10 區(qū)域S2真彩色圖像與兩種方法的判識結(jié)果Fig.10 True color image of area S2 and detective results of two methods

從表2中觀測的天氣現(xiàn)象與判識結(jié)果一致。為了在面上進(jìn)一步驗證算法的準(zhǔn)確性,還使用了MODIS的氣溶膠產(chǎn)品進(jìn)行比對。如圖10(d)所示,通過添加的地理邊界線可以看出判識結(jié)果與氣溶膠光學(xué)厚度AOD大于0.25的區(qū)域(橙色)基本一致,MODIS的氣溶膠產(chǎn)品將沙塵區(qū)內(nèi)的積云進(jìn)行了最大程度的濾除,因此出現(xiàn)了一些孤立區(qū),本文算的予以了填充。

根據(jù)阿拉善盟氣象臺2019年4月20日發(fā)布的氣象信息,4月20日8時49分至13時57分阿拉善盟額濟(jì)納旗拐子湖地區(qū)出現(xiàn)揚沙天氣,區(qū)域S3即為這次局地性沙塵過程。

由圖11(a)可以看到沙塵區(qū)西部的厚度較大,其他部分則多呈現(xiàn)絮狀分布,加之沙漠地表的干擾導(dǎo)致多波段BTD幾乎不能識別,如圖11(b)所示。而本文方法識別出的區(qū)域與MODIS產(chǎn)品AOD值大于0.25[圖11(d)中的橙色區(qū)域AOD>0.25,黃色區(qū)域AOD>0.12],綠色區(qū)域為無沙塵AOD<0.12)有較高一致性,紅圈即為拐子湖站所在位置。

圖11 區(qū)域S3真彩色圖像與兩種方法的判識結(jié)果Fig.11 True color image of area S3 and detective results of two methods

通過對三個典型的沙塵過程的判識結(jié)果的分析可以發(fā)現(xiàn),相比傳統(tǒng)方法在弱沙塵識別上的局限性,本文算法在各種不同強度條件下的沙塵判識中取得了較理想的效果,提升了沙塵遙感判識的準(zhǔn)確性和適用性。

為了進(jìn)一步論證判識及強度算法的適用性,選取了2017年5月3日16時30分的一幅影像進(jìn)行綜合判識,并用進(jìn)行沙塵強度估算。2017年5月3日發(fā)生了一次持續(xù)時間長、波及范圍廣的強沙塵過程,主要是由一個較強的冷渦和蒙古氣旋共同影響所導(dǎo)致的,影響了中國內(nèi)蒙古及北方大部分地區(qū)。如圖12所示為真彩色圖像和本文算法及沙塵強度的估算結(jié)果。表3列出了2017年5月3日17時的27個地面氣象站觀測的天氣現(xiàn)象與本文算法識別的結(jié)果,這27個站均識別有沙塵,有浮塵、揚沙、沙塵暴三種不同強度類型的觀測。本文根據(jù)能見度標(biāo)準(zhǔn)劃分的沙塵強度等級與EDII值進(jìn)行對應(yīng),給出EDII值對不同類型沙塵的強度等級劃分:浮塵:0.4≤EDII≤0.8;揚沙:0.8

圖12 2017年5月3日16時30分影像與本文方法沙塵識別結(jié)果及沙塵強度估算結(jié)果 Fig.12 The images of 16:30 on May 3, 2017 and detective results and intensity by the method in this paper

表3 2019年5月3日17時站點觀測與算法識別結(jié)果Table 3 The results by observation and method of 17:00 on May 3, 2019

可以看出22個站的地面觀測與判識及估算結(jié)果一致,有5個站沙塵強度的估算結(jié)果與地面觀測有偏差。其中53 149觀測為沙塵暴,估算為揚沙;53 336、53 419、53 502、53 513觀測為揚沙,估算為沙塵暴。從這5個站的結(jié)果可以看出,誤差中高估是主要的,原因可能是沙塵過程往往存在含沙量垂直分布不均的情況,在高空輸送較強時上層含沙量較高,地面含沙量相對較低,造成地面觀測較衛(wèi)星估算結(jié)果偏低。通過與地面觀測結(jié)果比較可以發(fā)現(xiàn),本文算法明顯優(yōu)于傳統(tǒng)BTD算法,特別是對弱沙塵,伴隨著積云的沙塵,云系邊緣沙塵等情況。對于這些復(fù)雜條件下的沙塵,傳統(tǒng)方法都表現(xiàn)出了一定的局限性,而本文方法則能夠進(jìn)行有效識別和較高精度的強度估算。

4 結(jié)論與討論

由于衛(wèi)星觀測參量在不同時間、地表狀態(tài)、海拔等具有較大差異,因此使用固定閾值存在局限,使用動態(tài)閾值也具有較大波動性。采用Himawari-8衛(wèi)星數(shù)據(jù),提出一種針對性的識別方法,基于動態(tài)亮溫差底圖引入了沙塵覆蓋度指數(shù)FDCI。該方法既能通過構(gòu)造動態(tài)背景彌補亮溫等變化產(chǎn)生的固定閾值局限,也能通過指數(shù)化的方法將判斷閾值進(jìn)行固定化,避免動態(tài)閾值選取的不確定性。

結(jié)果表明該指數(shù)在一個較確定范圍內(nèi)對不同強度沙塵都可以進(jìn)行有效判識,并能提升對弱沙塵的識別效果。在此基礎(chǔ)上本文還構(gòu)建了沙塵強度指數(shù)EDII模型,其中對2017年5月3日沙塵強度的地面驗證表明,位于觀測區(qū)的27個站中有22個站的地面觀測沙塵強度與估算結(jié)果一致,具備業(yè)務(wù)應(yīng)用價值。適用于各種地表及強度條件下的沙塵監(jiān)測,該方法是對分裂窗亮溫差法的有效改進(jìn)。將為中國新一代靜止氣象衛(wèi)星FY-4B(已發(fā)射,在軌測試中)的沙塵遙感業(yè)務(wù)提供更多有益參考。由于國產(chǎn)衛(wèi)星在通道設(shè)置、輻射定標(biāo)等方面存在一定差異,因此國產(chǎn)衛(wèi)星的適用性研究仍需做針對性研究。

但是,對于一些極弱的地面性浮塵天氣,衛(wèi)星監(jiān)測條件下的FDCI對其敏感性低,仍有較大的判識難度;另一方面,在大范圍復(fù)雜云系影響下判識的準(zhǔn)確性也會受到影響。大范圍沙塵過程有時會有一些云系與沙塵混合的情況,使整體亮溫值降低、亮溫差值變大,對沙塵判識結(jié)果造成一定影響,特別是一些薄云下的沙塵,容易被當(dāng)云掩膜處理掉。此外,由于強沙塵過程常涉及夜間判識,而本文基于FDCI構(gòu)建了EDII的沙塵強度,但夜間沒有能見度,因此并未對算法夜間的適用性進(jìn)行檢驗。研究過程中發(fā)現(xiàn)方法對于沙塵強度存在一定的高估,在以后的研究中或可結(jié)合激光雷達(dá)CALIPOSO垂直結(jié)構(gòu)做進(jìn)一步修正的探究。