基于IR-MAD不變像元的風(fēng)云三號(hào)光學(xué)成像儀同平臺(tái)交叉定標(biāo)

李國(guó)榮，何玉青，胡秀清，王俊偉

1.北京理工大學(xué) 光電學(xué)院 光電成像技術(shù)與系統(tǒng)教育部重點(diǎn)實(shí)驗(yàn)室,北京 100081;2.國(guó)家衛(wèi)星氣象中心（國(guó)家空間天氣監(jiān)測(cè)預(yù)警中心）中國(guó)氣象局,北京 100081

1 引言

可見(jiàn)光紅外掃描輻射計(jì)（VIRR）是中國(guó)風(fēng)云系列極軌衛(wèi)星搭載的主要儀器之一，其積累了20年10 通道對(duì)地觀測(cè)數(shù)據(jù)，這些數(shù)據(jù)對(duì)天氣分析、氣候預(yù)測(cè)和環(huán)境災(zāi)害業(yè)務(wù)應(yīng)用發(fā)揮了重要作用。衛(wèi)星傳感器在軌運(yùn)行期間受其運(yùn)行環(huán)境和器件老化等因素的影響，輻射探測(cè)性能和穩(wěn)定性會(huì)發(fā)生改變（Xu等，2014），嚴(yán)重影響其所獲取輻射觀測(cè)信號(hào)的準(zhǔn)確性，因此需要對(duì)衛(wèi)星傳感器進(jìn)行有效的輻射定標(biāo)，為其數(shù)據(jù)產(chǎn)品的有效性和可用性提供保障。VIRR 沒(méi)有配備相應(yīng)的星載定標(biāo)設(shè)備，主要依靠每年夏季在敦煌輻射校正場(chǎng)進(jìn)行的場(chǎng)地替代定標(biāo)方法進(jìn)行定標(biāo)校正（Hu等，2010）。然而一年一次的場(chǎng)地定標(biāo)實(shí)驗(yàn)無(wú)法支持高頻次和穩(wěn)定的在軌校準(zhǔn)系數(shù)更新，缺少對(duì)VIRR 輻射響應(yīng)衰變的連續(xù)性監(jiān)測(cè)和評(píng)估。此外，敦煌輻射校正場(chǎng)地表脆弱，晴空天氣少，容易受到TSP不穩(wěn)定和沙塵天氣因素影響（王敏 等，2015），限制了VIRR 場(chǎng)地定標(biāo)的精度。

交叉定標(biāo)方法以具有高絕對(duì)定標(biāo)精度的傳感器為基準(zhǔn)，通過(guò)對(duì)同一目標(biāo)同時(shí)觀測(cè)，對(duì)待定標(biāo)的傳感器進(jìn)行相對(duì)定標(biāo)，具有不受場(chǎng)地同步觀測(cè)條件限制、有效提高定標(biāo)頻次的優(yōu)勢(shì)。Tong 等（2005）基于MODIS 敦煌輻射校正場(chǎng)的同步觀測(cè)數(shù)據(jù)，對(duì)FY-1D/VIRR 的4 個(gè)VNIR 波段進(jìn)行校準(zhǔn)，部分通道校準(zhǔn)精度達(dá)到5%；Zhong等（2018）利用MODIS反演了巴丹吉林沙漠場(chǎng)地BRDF和氣溶膠光學(xué)厚度（AOD），得到2012 年內(nèi)FY-3A/VIRR 儀器46 次交叉校準(zhǔn)結(jié)果，誤差約為5%。探測(cè)器輻射響應(yīng)存在一定目標(biāo)依賴性，基于單個(gè)場(chǎng)地展開(kāi)的場(chǎng)地定標(biāo)或交叉定標(biāo)，均存在場(chǎng)地輻射特性單一、反射率動(dòng)態(tài)范圍小的問(wèn)題。為了更好的表征儀器寬動(dòng)態(tài)范圍內(nèi)的輻射響應(yīng)情況，韋瑋（2017）利用全球定標(biāo)場(chǎng)網(wǎng)提高定標(biāo)頻次和目標(biāo)輻射范圍，得到的FY-3B/VIRR 全球定標(biāo)場(chǎng)網(wǎng)定標(biāo)結(jié)果與官方定標(biāo)結(jié)果相對(duì)偏差為5.15%；Wang 等（2018）使用多站點(diǎn)（multisite，MST）交叉校準(zhǔn)方法，將不同衛(wèi)星平臺(tái)搭載的VIRR 大氣頂反射率差異從5%—10%降低至3%?？紤]到定標(biāo)頻次和定標(biāo)成本；何靈莉等（2020）使用自動(dòng)化輻射定標(biāo)方法對(duì)FY-3C/VIRR 進(jìn)行了長(zhǎng)時(shí)間序列輻射定標(biāo)，3 年時(shí)間尺度上可定標(biāo)觀測(cè)100 天，在VNIR 具有3%的定標(biāo)精度，在SWIR 具有5%的定標(biāo)精度。探索一種自動(dòng)化、高頻次、多目標(biāo)特點(diǎn)的交叉定標(biāo)方法將有助于解決目前VIRR的現(xiàn)實(shí)定標(biāo)問(wèn)題。

隨著計(jì)算機(jī)信息技術(shù)的進(jìn)步以及統(tǒng)計(jì)分析方法在科學(xué)研究中的廣泛應(yīng)用，遙感科研工作者逐漸將一些數(shù)學(xué)方法應(yīng)用于衛(wèi)星數(shù)據(jù)的處理分析中。研究者們提出實(shí)現(xiàn)多光譜圖像相對(duì)輻射標(biāo)準(zhǔn)化的方法是基于一個(gè)理想假設(shè)下進(jìn)行的：在兩個(gè)不同時(shí)間獲取的來(lái)自恒定反射區(qū)域的傳感器輻射之間的關(guān)系在空間上是均勻的，并且可用一個(gè)線性函數(shù)來(lái)近似（Moran 等，1992；Du 等，2002）。Nielsen等（1998）利用典型相關(guān)性分析（CCA）對(duì)兩組衛(wèi)星通道數(shù)據(jù)進(jìn)行線性組合，提出了一種具有線性尺度不變性的多變量變化檢測(cè)技術(shù)（MAD），可以消除同時(shí)相傳感器通道內(nèi)部的相關(guān)性以及雙時(shí)相傳感器不同通道間的相關(guān)性。Canty 等（2004）基于MAD 方法建立了一種自動(dòng)化的變化檢測(cè)算法用于獲取衛(wèi)星觀測(cè)場(chǎng)景中的不變特征。Schmidt等（2008）應(yīng)用MAD 方法分析了AVHRR 傳感器長(zhǎng)時(shí)間序列的響應(yīng)變化，結(jié)果證明了此方法的有效性。隨后Canty 和Nielsen（2008）提出了基于迭代重加權(quán)的多變量變化檢測(cè)方法（IR-MAD），對(duì)于變化像元占主要部分的場(chǎng)景，例如季節(jié)性更替明顯的耕地牧場(chǎng)等，IR-MAD方法在檢測(cè)不變特征方面表現(xiàn)優(yōu)異。王俊偉等（2019）利用IR-MAD方法識(shí)別同一場(chǎng)景不同時(shí)相衛(wèi)星圖像中的不變像元，并基于這些不變像元的輻射變化對(duì)FY-3A/MERSI 進(jìn)行儀器響應(yīng)衰變分析，其結(jié)果與傳統(tǒng)方法獲取的儀器衰變結(jié)果具有很好的一致性。IRMAD 方法對(duì)同一傳感器不同時(shí)相場(chǎng)景中不變像元的檢測(cè)具有優(yōu)異性，對(duì)于具有相似通道的兩個(gè)傳感器，其相似通道在恒定反射區(qū)域獲取的輻射量也具有明顯的線性關(guān)系，因此此方法可以擴(kuò)展到不同傳感器同時(shí)相場(chǎng)景的不變像元檢測(cè)中，以用于交叉定標(biāo)分析。

本文針對(duì)FY-3B 上搭載的可見(jiàn)光紅外掃描輻射計(jì)（VIRR）定標(biāo)精度差、定標(biāo)頻次低的問(wèn)題，以同平臺(tái)的中分辨率光譜成像儀（MERSI）為基準(zhǔn)，基于IR-MAD 方法對(duì)兩傳感器同時(shí)相觀測(cè)場(chǎng)景中的不變像元進(jìn)行自動(dòng)化檢測(cè)，利用不變像元的輻射信息展開(kāi)交叉定標(biāo)，并對(duì)交叉定標(biāo)結(jié)果和長(zhǎng)時(shí)間序列監(jiān)測(cè)情況進(jìn)行驗(yàn)證分析。

2 研究數(shù)據(jù)與區(qū)域

2.1 FY-3B/VIRR與MERSI

FY-3B 是中國(guó)第二代極軌氣象衛(wèi)星中的第二顆衛(wèi)星，可見(jiàn)光紅外掃描輻射計(jì)（VIRR）和中分辨率成像光譜儀（MERSI）是其中的兩個(gè)主要儀器。FY-3B/VIRR 共有10 個(gè)通道波段，光譜范圍0.43—12.5 μm，其中有7 個(gè)可見(jiàn)光近紅外波段和3 個(gè)紅外發(fā)射波段，星下點(diǎn)空間分辨率為1.1 km。FY-3B/MERSI 具有19 個(gè)反射太陽(yáng)波段（0.41—2.13 μm）和1 個(gè)熱紅外波段（11.25 μm）。VIRR的1，2，6，7，8，9 通道與MERSI 的3，4，6，1，10_11，2 通道設(shè)置相似（VIRR 的8 通道對(duì)應(yīng)MERSI 的10 和11 兩個(gè)通道），具體指標(biāo)見(jiàn)表1，故將MERSI 作為參考數(shù)據(jù)，對(duì)VIRR 共計(jì)6 組對(duì)應(yīng)通道進(jìn)行交叉定標(biāo)。本研究使用的數(shù)據(jù)為FY-3B/VIRR 和MERSI 2011 年1 月21 日 至2018 年11 月14日近8年的1 km空間分辨率的L1數(shù)據(jù)。

表1 FY-3B/VIRR、MERSI 光譜波段指標(biāo)Table 1 Spectral band specifications of FY-3B/VIRR、MERSI

一般的交叉定標(biāo)研究需考慮不同衛(wèi)星傳感器由于觀測(cè)時(shí)間、觀測(cè)幾何條件、輻射光譜響應(yīng)、大氣特性差異以及目標(biāo)輻射特性帶來(lái)的不確定性（Lacherade 等，2013）。本文研究?jī)x器VIRR 及MERSI 位于同一衛(wèi)星平臺(tái)，且兩探測(cè)器視場(chǎng)掃描范圍均為±55.4°，在視場(chǎng)大小、觀測(cè)時(shí)間、觀測(cè)幾何、大氣條件等方面一致，這有利于IR-MAD 算法對(duì)場(chǎng)景內(nèi)不變像元的檢測(cè)，同時(shí)可近似認(rèn)為兩探測(cè)器的輻射標(biāo)準(zhǔn)化差異僅由儀器本身響應(yīng)特性和光譜差異引起（Obata等，2017）。

2.2 研究目標(biāo)區(qū)域

本文以中國(guó)西北地區(qū)為主要研究區(qū)域（圖1（a）），中心坐標(biāo)為（91°E、39°N）；主要包括內(nèi)蒙古自治區(qū)，青海省，甘肅省和新疆維吾爾自治區(qū)；地形以高原、盆地、山地為主，海拔1 km 以上；地貌景觀主要為黃土高原、戈壁灘、荒漠草原和戈壁沙漠。此區(qū)域?qū)儆诖箨懶詺夂颍芟柴R拉雅山脈南部的影響，冬季寒冷干燥，夏季炎熱，降水稀少，干旱是該地區(qū)的主要自然特征。研究區(qū)域地表反射率季節(jié)性變化小，符合理想的低氣溶膠量、低水汽含量和高晴空概率的交叉定標(biāo)觀測(cè)區(qū)域標(biāo)準(zhǔn)。同時(shí)為驗(yàn)證所提出方法的通用性，選用同時(shí)期北非地區(qū)（圖1（b））數(shù)據(jù)進(jìn)行比對(duì)驗(yàn)證。北非研究區(qū)域中心坐標(biāo)為（15°E、23°N），以沙漠荒漠地貌為主，氣候干旱，包含有地球觀測(cè)衛(wèi)星委員會(huì)和可見(jiàn)光遙感器定標(biāo)和驗(yàn)證工作小組選定的多個(gè)國(guó)際校準(zhǔn)目標(biāo)場(chǎng)地，其地表反射率的穩(wěn)定性和均勻性良好（Lacherade等，2013）。

3 基于IR-MAD 不變像元的交叉定標(biāo)方法

基于IR-MAD 不變像元進(jìn)行VIRR 與MERSI 交叉定標(biāo)的處理過(guò)程如圖2 所示。首先對(duì)同時(shí)相的VIRR 和MERSI 圖像數(shù)據(jù)對(duì)進(jìn)行篩選與預(yù)處理，將兩探測(cè)器的匹配通道數(shù)據(jù)進(jìn)行線性組合以構(gòu)造MAD 變量，經(jīng)IR-MAD 不變像元檢測(cè)模型分析得到場(chǎng)景內(nèi)的不變像元樣本?；讷@得的不變像元樣本，以MERSI 探測(cè)結(jié)果為基準(zhǔn)，經(jīng)光譜匹配修正后，與VIRR 探測(cè)結(jié)果進(jìn)行正交回歸以得到交叉定標(biāo)系數(shù)。

圖2 基于IR-MAD不變像元的FY-3B/VIRR與MERSI交叉定標(biāo)流程圖Fig.2 Cross calibration flow chart of FY-3B/VIRR and MERSI based on IR-MAD No-change pixels

3.1 數(shù)據(jù)篩選和前處理

FY-3B/VIRR 與MERSI同平臺(tái)觀測(cè)，且可實(shí)現(xiàn)每天一次的全球掃描覆蓋，因此以每日獲取的研究區(qū)域同時(shí)相觀測(cè)場(chǎng)景為一組匹配圖像數(shù)據(jù)對(duì)。兩傳感器獲取的通道DN 值數(shù)據(jù)中，均存在失效數(shù)據(jù)。IR-MAD方法基于統(tǒng)計(jì)學(xué)原理展開(kāi)，輸入模型的兩組數(shù)據(jù)維度一致，故需要將兩圖像中的無(wú)效像元對(duì)應(yīng)剔除。

IR-MAD 算法用于檢測(cè)場(chǎng)景中的不變像元，變化像元會(huì)在算法迭代過(guò)程中逐步剔除，考慮到算法的運(yùn)行效率，對(duì)于場(chǎng)景中顯著的非穩(wěn)定像元，如云、沙塵等目標(biāo)應(yīng)在算法運(yùn)行前予以剔除。研究區(qū)域顯著變化場(chǎng)景主要是云，VIRR 和MERSI 的L1 數(shù)據(jù)沒(méi)有云掩膜數(shù)據(jù)，因此本文使用VIRR 的通道數(shù)據(jù)，采用閾值判別法進(jìn)行初步的云識(shí)別。首先利用VIRR 的1 通道（0.58—0.68 μm）和2 通道（0.84—0.89 μm）兩個(gè)波段的輻射比值區(qū)分云與晴空區(qū)域，再使用2 通道（0.84—0.89 μm）和6通道（1.55—1.64 μm）的輻射比值區(qū)分地上云層與地面高反射率物體（齊修東和岳繼博，2016）。由于去云處理的主要目的是提高算法的效率，因此云識(shí)別結(jié)果存在部分的云漏判不會(huì)影響到不變像元的檢測(cè)，這種簡(jiǎn)單的閾值法去云效果是可接受的。

水體及海洋等區(qū)域雖然有較為穩(wěn)定的反射率特性，但其反射特性光譜差異較大，且對(duì)于研究選取的部分通道波段信號(hào)微弱。因此對(duì)于海洋水體像元目標(biāo)，利用VIRR 與MERSI數(shù)據(jù)中的陸地海洋掩膜（landseamask）數(shù)據(jù)進(jìn)行剔除。

研究使用的西北地區(qū)數(shù)據(jù)單幅場(chǎng)景大小為1400×3200 像素，北非地區(qū)單幅場(chǎng)景大小為1600×2600 像素，經(jīng)過(guò)無(wú)效點(diǎn)、云和海洋水體像元去除后，仍有數(shù)百萬(wàn)數(shù)量級(jí)的像元樣本，其中存在大量衛(wèi)星觀測(cè)天頂角過(guò)大的像元樣本。衛(wèi)星觀測(cè)天頂角過(guò)大會(huì)導(dǎo)致像元空間分辨率和探測(cè)輻射精度下降，為了提高IR-MAD 對(duì)不變像元的識(shí)別精度，使用衛(wèi)星觀測(cè)天頂角小于30°的像元樣本用于后續(xù)的統(tǒng)計(jì)分析。

FY-3B/VIRR 和MERSI同時(shí)觀測(cè)同一地物目標(biāo)時(shí)，具有相同的觀測(cè)幾何和大氣路徑，避免了水汽氣溶膠、臭氧等大氣因素的影響，使用衛(wèi)星運(yùn)行初期的固定定標(biāo)系數(shù)，由VIRR 探測(cè)器通道DN值計(jì)算像元表觀反射率：

式中，slopei、biasi分別是探測(cè)器第i通道對(duì)應(yīng)的固定定標(biāo)斜率與截距，DNi為探測(cè)器獲取的DN 值，d2為日地距離修正因子，θs為太陽(yáng)天頂角。

FY-3B/MERSI在軌定標(biāo)分為兩個(gè)階段：2013年3 月6 日前采用線性定標(biāo)模型，基于敦煌輻射校正場(chǎng)替代定標(biāo)結(jié)果進(jìn)行在軌校正，定標(biāo)精度在波長(zhǎng)小于1 μm 的窗區(qū)波段優(yōu)于3%，波長(zhǎng)大于1 μm 的波段優(yōu)于5%（除2.1 μm 波段）（孫凌 等，2012）；2013年3月6日之后MERSI采用基于全球多目標(biāo)輻射定標(biāo)跟蹤建立的日更新模型定標(biāo)，其1—4、10、11波段與Aqua MODIS的TOA反射率的相對(duì)平均偏差在2%以內(nèi)，6通道在5%以內(nèi)（Sun和Li，2014）。為進(jìn)一步分析MERSI數(shù)據(jù)的穩(wěn)定性與同期VIRR 數(shù)據(jù)質(zhì)量，對(duì)兩儀器測(cè)得的敦煌輻射校正場(chǎng)TOA 反射率進(jìn)行時(shí)間序列追蹤，并計(jì)算各通道的年際相對(duì)平均偏差（RSD）（表2）。圖3（a）所示為MERSI采用日更新定標(biāo)模型后第一個(gè)自然年的時(shí)間序列結(jié)果，散點(diǎn)序列較為穩(wěn)定，其1，2，3，10，11通道的RSD 值小于4%，4，6 通道約為5%；圖3（b）為同時(shí)期VIRR 觀測(cè)結(jié)果，時(shí)間序列存在較為明顯的震蕩情況，研究選取的各通道RSD 值在6%—8%，數(shù)據(jù)穩(wěn)定性較差。MERSI 在定標(biāo)精度及數(shù)據(jù)穩(wěn)定性上較VIRR 有顯著優(yōu)勢(shì)，雖然整體精度較MODIS仍有一定差距，但考慮到MERSI與VIRR同平臺(tái)觀測(cè)的一致性優(yōu)勢(shì)以及應(yīng)用到IR-MAD 匹配場(chǎng)景的便利性，將其作為同期VIRR 的交叉定標(biāo)參考基準(zhǔn)是可行的。

圖3 2013年3月—2014年3月FY-3B/VIRR與MERSI敦煌場(chǎng)地TOA反射率時(shí)間序列Fig.3 TOA reflectance time series of FY-3B/VIRR and MERSI in Dunhuang site from March 2013 to March 2014

表2 2013年3月—2014年3月VIRR與MERSI敦煌場(chǎng)TOA反射率序列相對(duì)標(biāo)準(zhǔn)偏差Table 2 RSD of TOA series between VIRR and MERSI in Dunhuang Site from March 2013 to March 2014

3.2 IR-MAD不變像元檢測(cè)模型

將IR-MAD不變像元檢測(cè)應(yīng)用于同一時(shí)間VIRR和MERSI采集的同一場(chǎng)景的兩幅n通道多光譜圖像。

分別用向量F=(F1，…，F(xiàn)n)T、G=(G1，…，Gn)T表示VIRR、MERSI 多光譜圖像中DN 值轉(zhuǎn)換后的表觀反射率，對(duì)所有譜帶進(jìn)行線性組合以獲得典型變量U、V：

式中，m和n是常向量，可通過(guò)求解耦合廣義特征值方程得出：

式中，ρ為典型變量U、V的相關(guān)系數(shù)，∑ff、∑gg、∑fg、∑gf為圖像向量F、G的協(xié)方差矩陣。

MAD變量由典型變量U、V的差分得到，即：

由于MAD 變量具有線性尺度不變性，因此我們選擇像元樣本中滿足

單次的MAD 變換難以達(dá)到理想的不變像元檢測(cè)效果。通過(guò)單次迭代確定的不變概率對(duì)樣本數(shù)據(jù)進(jìn)行加權(quán)，賦予更高不變概率的像元更大的權(quán)重，使用典型相關(guān)分析確定下一次迭代的MAD 變量，多次迭代以獲得更好的不變像元檢測(cè)效果。對(duì)于每次迭代，不變概率權(quán)重值Pr可由觀測(cè)值Z的卡方分布檢驗(yàn)結(jié)果確定：

為尋求找到更大概率（90%以上）的不變像元，當(dāng)自由度n=6（6 組匹配通道）時(shí)，對(duì)90%、92.5%、95%、97.5%的不變概率置信度進(jìn)行測(cè)試分析，綜合考慮算法運(yùn)行效率和不變像元檢測(cè)效果，選擇t==1.635 作為不變概率決策閾值，當(dāng)觀測(cè)值Z＜t=1.635時(shí)，即可認(rèn)為該像元樣本為不變像元的置信度高于95%，可用于后續(xù)交叉定標(biāo)分析。

對(duì)于IR-MAD算法，需要設(shè)定典型變量相關(guān)系數(shù)變化閾值、最大迭代次數(shù)和最小NCPs 數(shù)量3 個(gè)迭代停止閾值（王俊偉 等，2019），使其在算法優(yōu)化收斂以及小概率檢測(cè)失效時(shí)停止對(duì)某一圖像對(duì)的迭代處理，以保證算法的自動(dòng)運(yùn)行。本研究設(shè)定兩次典型變量相關(guān)系數(shù)的變化差值小于0.001 時(shí)即可視為算法收斂并停止迭代，并設(shè)定算法最大迭代次數(shù)為30 次。同時(shí)，為保證有足夠的不變像元樣本用于后續(xù)的交叉定標(biāo)回歸分析，設(shè)定最小的NCPs 數(shù)量為400，當(dāng)某次迭代后NCPs 數(shù)量小于此數(shù)值，迭代停止。

IR-MAD 算法從像元樣本的輻射信息出發(fā)，構(gòu)造通道間線性組合，消除單個(gè)傳感器不同通道以及不同傳感器非匹配通道間的相關(guān)性，經(jīng)過(guò)多次不變概率加權(quán)迭代，獲得兩幅觀測(cè)場(chǎng)景中的不變像元，不變像元輻射量具有對(duì)VIRR 與MERSI匹配通道間的最大線性相關(guān)性。對(duì)不變像元的空間分布特征進(jìn)行分析（圖4），圖4（a）為2013 年6 月3 日獲取的不變像元檢測(cè)結(jié)果，寬幅場(chǎng)景圖像包含兩個(gè)FY-3B軌道，算法檢測(cè)范圍包含被分割的左右兩塊觀測(cè)天頂角＜30°的區(qū)域，檢測(cè)得到的不變像元（紅色標(biāo)記點(diǎn)）主要分布于新疆塔克拉瑪干沙漠、內(nèi)蒙巴丹吉林沙漠及騰格里沙漠區(qū)域。圖4（b）為研究區(qū)域內(nèi)的像元不變概率分布圖，通過(guò)對(duì)獲取的長(zhǎng)時(shí)間序列不變像元結(jié)果進(jìn)行概率頻次的疊加得到，紅色表示像元的不變概率更高。時(shí)間序列結(jié)果顯示，高概率不變像元主要聚集在塔克拉瑪干沙漠，敦煌沙漠，巴丹吉林沙漠，騰格里沙漠和一些小戈壁地區(qū)，這與使用VIRR 以及MODIS等單個(gè)傳感器不同時(shí)相場(chǎng)景獲取的分布情況相一致（Hu等，2020）。相較于場(chǎng)景內(nèi)其他目標(biāo)，不變像元具有更好的時(shí)間穩(wěn)定性，基于穩(wěn)定定標(biāo)場(chǎng)的交叉定標(biāo)也常選取沙漠目標(biāo)進(jìn)行，這體現(xiàn)了IRMAD方法在目標(biāo)選用上的合理性。

圖4 中國(guó)西北地區(qū)不變像元分布結(jié)果Fig.4 NCPs distribution in Northwest China region

3.3 光譜匹配

VIRR與MERSI對(duì)應(yīng)通道的光譜響應(yīng)（圖5）存在差異，即對(duì)于相同的入瞳輻射量，兩傳感器會(huì)得到不同的測(cè)量值。因此需要對(duì)傳感器對(duì)應(yīng)通道進(jìn)行光譜匹配，以修正由于光譜響應(yīng)差異造成的定標(biāo)誤差。同平臺(tái)觀測(cè)條件下（相同的地表、大氣和觀測(cè)時(shí)間幾何條件），兩傳感器匹配通道的入同輻射量之比稱為光譜匹配因子（SBAF）。使用SCIAMACHY儀器星下點(diǎn)觀測(cè)模式的數(shù)據(jù)作為高光譜樣本（圖5），分別與VIRR 和MERSI的光譜響應(yīng)函數(shù)卷積，得到相應(yīng)傳感器的入瞳輻亮度，關(guān)系式為

圖5 FY-3B/VIRR、MERSI 通道光譜響應(yīng)函數(shù)及SCIAMACHY高光譜樣本Fig.5 Channel SRF of FY-3B/VIRR and MERSI，IAMACHY hyperspectral samples

式中，Rscia為SCIAMACHY 光譜樣本輻亮度，fi_sensor為儀器i通道的光譜響應(yīng)函數(shù)，Ri_sensor為卷積后的儀器入瞳輻亮度。

建立VIRR 與MERSI匹配通道入瞳幅亮度間的關(guān)系式：

式中，R表示輻亮度，A、B為光譜匹配因子，i、j分別為VIRR與MERSI的匹配通道序號(hào)，對(duì)于VIRR與MERSI 的單對(duì)單匹配通道使用式（8），VIRR 的8 通道和MERSI 的10、11 通道使用式（9），通過(guò)最小二乘回歸擬合計(jì)算SBAF，結(jié)果如表3。

表3 FY-3B VIRR與MERSI光譜匹配因子（SBAF）Table 3 Spectral band adjustment factor of FY-3B VIRR and MERSI

3.4 回歸分析與定標(biāo)計(jì)算

對(duì)IR-MAD 檢測(cè)獲取的不變像元，以MERSI探測(cè)結(jié)果作為輻射基準(zhǔn)，利用光譜匹配因子計(jì)算VIRR 表觀反射率準(zhǔn)真值，將其與VIRR 探測(cè)結(jié)果建立線性擬合關(guān)系，從而獲得交叉定標(biāo)系數(shù)。

式中，a為定標(biāo)斜率，b為定標(biāo)截距。ρVIRR為基于發(fā)射初期固定定標(biāo)斜率計(jì)算得到的表觀反射率，因此定標(biāo)系數(shù)a，b是相較于發(fā)射初期固定定標(biāo)系數(shù)的相對(duì)定標(biāo)結(jié)果。

本研究區(qū)域廣闊，單日?qǐng)D像對(duì)結(jié)果僅包含場(chǎng)景中符合觀測(cè)天頂角限制條件的部分區(qū)域信息，為盡可能在單次回歸中包含更豐富的不變像元樣本，將連續(xù)多日檢測(cè)得到的不變像元合并處理?；趥鞲衅鬏椛漤憫?yīng)的短期穩(wěn)定性，這種合并處理是可行的，并將有助于擴(kuò)大像元樣本的反射率值范圍，有效提升回歸效果與定標(biāo)精度。考慮到風(fēng)云三號(hào)衛(wèi)星標(biāo)稱軌道回歸周期為5.5 天，同時(shí)對(duì)于識(shí)別出的不變像元樣本，Canty 等（2004）分析得出正交回歸是更好的回歸處理方法，對(duì)連續(xù)5天檢測(cè)獲取的不變像元合并進(jìn)行正交回歸分析。圖6為2011 年4 月8 日—12 日不變像元各組匹配通道正交回歸結(jié)果。

圖6 2011年4月8日—12日不變像元表觀反射率正交回歸結(jié)果Fig.6 Orthogonal regression results of TOA reflectance of NCPs from April 8 to 12，2011

4 交叉定標(biāo)結(jié)果與分析

4.1 交叉定標(biāo)質(zhì)量精度分析

基于上述的IR-MAD 不變像元檢測(cè)和交叉定標(biāo)回歸分析，對(duì)2011 年1 月21 日至2018 年11 月14 日數(shù)據(jù)進(jìn)行分析處理，共得到中國(guó)西北地區(qū)、北非地區(qū)的5 日合并正交回歸結(jié)果數(shù)量分別為523、551 組。為保證回歸結(jié)果的質(zhì)量，對(duì)每個(gè)回歸結(jié)果的NCPs 樣本表觀反射率相關(guān)系數(shù)進(jìn)行考察，圖7 給出了所有回歸結(jié)果8-10_11 匹配通道的NCPs 樣本反射率相關(guān)系數(shù)分布情況，相關(guān)系數(shù)值大于0.99，對(duì)其余匹配通道相關(guān)系數(shù)小于0.99 的結(jié)果進(jìn)行剔除。

圖7 8—10_11匹配通道的圖像對(duì)NCPs樣本相關(guān)系數(shù)Fig.7 NCPs sample correlation coefficients of image pairs with 8-10_11 matching channels

為驗(yàn)證基于IR-MAD 不變像元交叉定標(biāo)方法所得到的VIRR 交叉定標(biāo)結(jié)果精度，對(duì)2012 年1 月15 日—2 月15 日中國(guó)西北、北非地區(qū)IR-MAD 交叉定標(biāo)結(jié)果與2012 年2 月1 日VIRR 業(yè)務(wù)定標(biāo)值進(jìn)行對(duì)比，如圖8所示。需指出的是，本文后續(xù)表中所列斜率結(jié)果均為相對(duì)定標(biāo)系數(shù)與發(fā)射初期固定定標(biāo)系數(shù)的乘積結(jié)果。如表4所列具體數(shù)值，中國(guó)西北地區(qū)IR-MAD 自然月內(nèi)交叉定標(biāo)斜率的相對(duì)標(biāo)準(zhǔn)偏差值小于1.5%，表明了定標(biāo)結(jié)果的穩(wěn)定性。除VIRR 7 通道外，其余各通道的IR-MAD 中國(guó)西北與北非的定標(biāo)結(jié)果一致性優(yōu)于2%，說(shuō)明所述方法基本不受地理區(qū)域選擇限制，有較好的通用性。VIRR 7 通道為短波藍(lán)色光波段，受大氣和氣溶膠分子散射影響明顯，雖然用于IR-MAD 分析的圖像對(duì)大氣狀況幾乎一致，但同期中國(guó)西北和北非兩區(qū)域大氣狀況存在顯著差異，這造成了兩區(qū)域所得VIRR 7 通道IR-MAD 定標(biāo)結(jié)果相較于其余通道的較大偏差（3.73%）。

圖8 中國(guó)西北、北非地區(qū)IR-MAD定標(biāo)結(jié)果與2012年2月1日業(yè)務(wù)定標(biāo)結(jié)果對(duì)比Fig.8 Comparison of calibration results between IR-MAD in Northwest China and North Africa with business calibration results on February 1，2012

表4 中國(guó)西北、北非地區(qū)IR-MAD交叉定標(biāo)結(jié)果與VIRR業(yè)務(wù)定標(biāo)對(duì)比Table 4 Comparison of cross calibration results between IR-MAD in northwest China and north Africa with VIRR business calibration

同時(shí)，IR-MAD 西北地區(qū)各通道的交叉定標(biāo)結(jié)果與業(yè)務(wù)定標(biāo)結(jié)果的相對(duì)偏差均小于4%，1 通道（紅色波段）和2 通道（近紅外波段）的結(jié)果最為接近，符合這兩個(gè)通道的較為穩(wěn)定的特征。

VIRR 業(yè)務(wù)定標(biāo)頻次很低，為進(jìn)一步驗(yàn)證定標(biāo)結(jié)果精度，將IR-MAD 結(jié)果與基于敦煌輻射校正場(chǎng)的FY-3B/VIRR 定標(biāo)結(jié)果（呂佳彥 等，2017）進(jìn)行比較，結(jié)果如表5。兩方法定標(biāo)結(jié)果顯示，VIRR 1，2通道的相對(duì)偏差基本在2%以內(nèi)，7，8，9 通道的相對(duì)偏差低于5%。VIRR 7，8，9 這3 個(gè)通道是藍(lán)綠光波段，相較于紅色、近紅外和短波紅外波段，VIRR 和MERSI 在這幾個(gè)波段有較大的衰變情況，MERSI 的業(yè)務(wù)定標(biāo)結(jié)果對(duì)這種衰變的輻射校正精度在一定程度上影響了IR-MAD 方法的交叉定標(biāo)精度。此外，兩方法結(jié)果顯示的各通道定標(biāo)斜率的年際變化趨勢(shì)一致，其中6通道（短波紅外）表現(xiàn)出與其他5個(gè)通道不一致的趨勢(shì)，具體為先增大后回落，此情況發(fā)生的原因可能是短波紅外通道電子增益出現(xiàn)隨機(jī)躍動(dòng)，導(dǎo)致定標(biāo)系數(shù)存在波動(dòng)，胡秀清等（2013）對(duì)風(fēng)云三號(hào)A 星MERSI 的兩個(gè)短波紅外通道響應(yīng)追蹤中也發(fā)現(xiàn)了類似現(xiàn)象。

表5 IR-MAD不變像元交叉定標(biāo)與敦煌場(chǎng)地定標(biāo)結(jié)果對(duì)比Table 5 Comparison of cross calibration results between IR-MAD NCPs and Dunhuang site calibration

表6 給出了2012 年不變像元樣本數(shù)量及VIRR定標(biāo)后的不變像元反射率較MERSI 觀測(cè)的相對(duì)平均偏差?？梢悦黠@看出，冬季月份的不變像元樣本數(shù)量高于夏季月份，這是由于冬季晴空概率高且水汽含量低，場(chǎng)景中可用于統(tǒng)計(jì)分析的像元數(shù)量更多。定標(biāo)校正后，兩儀器表觀反射率的相對(duì)平均偏差小于2%。

表6 2012年不變像元樣本VIRR定標(biāo)反射率較MERSI觀測(cè)的相對(duì)平均偏差Table 6 Relative mean deviation of VIRR calibration reflectivity of NCPs compared with MERSI observation in 2012

4.2 長(zhǎng)時(shí)間序列監(jiān)測(cè)分析

圖9考察了2011年1月21日—2018年11月14日VIRR 各通道相對(duì)定標(biāo)斜率的長(zhǎng)時(shí)間序列趨勢(shì)，分別基于中國(guó)西北地區(qū)和北非地區(qū)得到的IR-MAD交叉定標(biāo)結(jié)果顯現(xiàn)出高度的一致性。從圖9中可以明顯發(fā)現(xiàn)，北非地區(qū)有更好的散點(diǎn)擬合線聚合情況，這與北非地區(qū)優(yōu)異的地表反射率穩(wěn)定性和空間均勻性相關(guān)。

圖9 2011年1月21日—2018年11月14日VIRR相對(duì)定標(biāo)斜率長(zhǎng)時(shí)間序列Fig.9 Long time series of VIRR relative calibration slope from January 21，2011 to November 14，2018

除VIRR 6 通道外，其余各通道自研究數(shù)據(jù)的初始時(shí)間（2011 年1 月21 日）開(kāi)始，距相對(duì)值1存在較大的差異。這是因?yàn)榛貧w分析所用的VIRR像元表觀反射率通過(guò)衛(wèi)星發(fā)射初期使用的固定定標(biāo)系數(shù)（發(fā)射前定標(biāo)結(jié)果）計(jì)算獲取，而在軌運(yùn)行一段時(shí)間后，儀器經(jīng)歷發(fā)射階段和在軌運(yùn)行期間的復(fù)雜環(huán)境，在軌輻射性能與發(fā)射前定標(biāo)結(jié)果存在較大差異。VIRR 的通道響應(yīng)衰減趨勢(shì)與波長(zhǎng)表現(xiàn)出一定的相關(guān)性，整體上呈現(xiàn)波長(zhǎng)越短，衰減越大的趨勢(shì)。VIRR 的7 和8 通道在軌運(yùn)行前兩年期間，響應(yīng)衰減幅度較小，第3—5 年出現(xiàn)大幅度的衰減，之后趨于相對(duì)平緩，這與呂佳彥等（2017）基于敦煌場(chǎng)地定標(biāo)的FY-3B/VIRR 通道衰減規(guī)律一致。

本文提出的方法可進(jìn)行近乎不間斷的自動(dòng)交叉定標(biāo)實(shí)施，在VIRR 2 和6 通道的中國(guó)西北地區(qū)定標(biāo)結(jié)果中發(fā)現(xiàn)了明顯的季節(jié)性波動(dòng)現(xiàn)象；2通道的北非地區(qū)結(jié)果也存在類似現(xiàn)象，但幅度較小；兩地區(qū)的1 通道結(jié)果也存在小幅度的季節(jié)性波動(dòng)。首先，本文IR-MAD 算法使用的同時(shí)相單日?qǐng)D像對(duì)大氣條件一致，因此未對(duì)TOA 反射率進(jìn)行大氣校正，大氣的季節(jié)性變化會(huì)帶來(lái)相應(yīng)的序列跟蹤不確定度；IR-MAD 算法從像素級(jí)展開(kāi)統(tǒng)計(jì)分析，不需要預(yù)先獲取地表的先驗(yàn)知識(shí)，不變像元的BRDF 季節(jié)性變化以及對(duì)地表均勻性產(chǎn)生影響的季節(jié)性降雪、沙塵等天氣活動(dòng)也是定標(biāo)序列波動(dòng)的重要影響因素。此外，孫凌等（2013）對(duì)FY-3A/MERSI 的多場(chǎng)地定標(biāo)跟蹤以及何靈莉等（2020）對(duì)FY-3C/VIRR 的自動(dòng)化定標(biāo)序列中均發(fā)現(xiàn)了季節(jié)性波動(dòng)的情況。除地表和大氣因素外，儀器本身的溫度響應(yīng)、通道光譜響應(yīng)特性也會(huì)造成不同季節(jié)的響應(yīng)差異，如VIRR 6 波段對(duì)云雪較為敏感，中國(guó)西北地區(qū)相較于北非地區(qū)，存在季節(jié)性的降雪活動(dòng)，由于在預(yù)處理中未進(jìn)行積雪識(shí)別和去除，這一定程度上導(dǎo)致了兩地區(qū)在6通道序列結(jié)果的不同表現(xiàn)。北非地區(qū)主要為低緯度區(qū)域，而中國(guó)西北地區(qū)為中緯度區(qū)域，其太陽(yáng)活動(dòng)和高度角的季節(jié)性變化更為明顯，因此近紅外波段的響應(yīng)衰減季節(jié)性變化也可能與近紅外波段太陽(yáng)輻射和地表反射率的分光譜特征季節(jié)相關(guān)性有關(guān)。

5 結(jié)論

本文針對(duì)風(fēng)云三號(hào)衛(wèi)星光學(xué)成像儀的定標(biāo)問(wèn)題和歷史數(shù)據(jù)再定標(biāo)需求，以FY-3B 衛(wèi)星的VIRR（待定標(biāo)）和MERSI（參考基準(zhǔn)）為研究對(duì)象，提出了基于IR-MAD 不變像元的交叉定標(biāo)方法。該方法利用同平臺(tái)傳感器獲取的同時(shí)相衛(wèi)星場(chǎng)景圖像，經(jīng)數(shù)據(jù)篩選與預(yù)處理后，通過(guò)IR-MAD 不變像元檢測(cè)模型自動(dòng)檢測(cè)場(chǎng)景中的不變像元目標(biāo)，計(jì)算光譜匹配因子以修正兩傳感器間的光譜響應(yīng)差異，建立不變像元樣本表觀反射率的線性擬合關(guān)系，得到接近8年的時(shí)間序列定標(biāo)結(jié)果。該方法有別于傳統(tǒng)的交叉定標(biāo)方法，不依賴衛(wèi)星過(guò)境人工選取的偽不變定標(biāo)場(chǎng)，也不受限于地面同步觀測(cè)的苛刻條件，能有效提高傳感器間的交叉定標(biāo)頻次。研究結(jié)果表明，IR-MAD 方法獲取的交叉定標(biāo)結(jié)果與業(yè)務(wù)定標(biāo)結(jié)果的差異小于4%；同時(shí)，定標(biāo)精度和通道波段有一定的相關(guān)性，與敦煌沙漠場(chǎng)定標(biāo)結(jié)果相比，VIRR 紅色和近紅外通道的定標(biāo)結(jié)果優(yōu)于3個(gè)藍(lán)綠短波通道。此外，中國(guó)西北和北非地區(qū)獲得的定標(biāo)結(jié)果長(zhǎng)時(shí)間序列表現(xiàn)出良好的一致性。

IR-MAD 方法對(duì)于不變像元目標(biāo)的檢測(cè)是自動(dòng)進(jìn)行的，對(duì)于處理長(zhǎng)時(shí)間序列數(shù)據(jù)和歷史衛(wèi)星數(shù)據(jù)再定標(biāo)有積極的作用。同時(shí)算法獲取的不變像元目標(biāo)在空間上是不連續(xù)的，雖然不能完整分析其地表的空間特性，但大量的不變像元樣本有效增加了反射率的覆蓋范圍，可適用于大場(chǎng)景數(shù)據(jù)，改善了傳統(tǒng)交叉定標(biāo)方法觀測(cè)目標(biāo)場(chǎng)地反射率范圍小、樣本單一的問(wèn)題。

目前基于全球定標(biāo)場(chǎng)網(wǎng)或自動(dòng)化輻射計(jì)的連續(xù)定標(biāo)方法，對(duì)于數(shù)量龐大的衛(wèi)星儀器，每個(gè)儀器均實(shí)施連續(xù)定標(biāo)的成本代價(jià)是巨大的。本文研究的IR-MAD 方法，在具有一個(gè)連續(xù)定標(biāo)精準(zhǔn)儀器的情況下，可實(shí)現(xiàn)對(duì)同平臺(tái)其他儀器的高效相對(duì)響應(yīng)監(jiān)測(cè)，解決寬動(dòng)態(tài)、全視場(chǎng)以及數(shù)據(jù)處理效率方面的問(wèn)題，有效降低定標(biāo)成本。IR-MAD方法已實(shí)現(xiàn)單個(gè)傳感器響應(yīng)衰變分析以及同平臺(tái)不同傳感器間相對(duì)響應(yīng)的交叉定標(biāo)，相比于以上兩種情況，異平臺(tái)傳感器間存在觀測(cè)幾何、時(shí)間、大氣以及光譜差異等誤差來(lái)源。在寬視角條件下，異平臺(tái)傳感器的觀測(cè)幾何條件匹配并不困難，光譜差異可通過(guò)光譜匹配因子加以校正，在輻射來(lái)源穩(wěn)定的情況下，時(shí)間因素的主要誤差來(lái)源于大氣變化，因此大氣差異的校正是實(shí)現(xiàn)異平臺(tái)儀器間IR-MAD 交叉定標(biāo)的關(guān)鍵，借助精細(xì)的大氣參量和輻射傳輸模型，我們對(duì)異平臺(tái)傳感器間的IRMAD算法表現(xiàn)有較好的期望。

本文方法在儀器的非線性響應(yīng)和場(chǎng)景的空間特征分析方面仍有不足。一些研究者們提出基于核典型相關(guān)分析（KCCA）的IR-MAD 方法和對(duì)MAD 變量進(jìn)行最大自相關(guān)因子（MAF）分析將有助于解決非線性和空間特性問(wèn)題，這將是我們優(yōu)化核心算法的工作方向。