[發(fā)明專利]一種基于多年時(shí)間序列數(shù)據(jù)的遙感圖像云污染去除方法

有效

• 申請(qǐng)?zhí)枺?/td>: 202010617922.2
• 申請(qǐng)日: 2020-06-30
• 公開(kāi)（公告）號(hào): CN111798394B
• 公開(kāi)（公告）日: 2022-10-14
• 發(fā)明（設(shè)計(jì)）人: 曹入尹;陳洋
• 申請(qǐng)（專利權(quán)）人: 電子科技大學(xué)
• 主分類號(hào): G06T5/00
• 分類號(hào): G06T5/00
• 代理公司: 電子科技大學(xué)專利中心 51203
• 代理人: 鄧?yán)?/td>

1.一種基于多年時(shí)間序列數(shù)據(jù)的遙感圖像云污染去除方法，其特征在于，該方法包括以下步驟：

S1.采集3年以上的原始Landsat圖像數(shù)據(jù)集作為參考圖像，確定需要消除云污染的Landsat圖像作為目標(biāo)圖像；

S2.目標(biāo)圖像云像元(x,y)第b波段的反射率值包含兩類估計(jì)值：長(zhǎng)時(shí)間估計(jì)值和短時(shí)間估計(jì)值，通過(guò)權(quán)重結(jié)合起來(lái)得到云像元(x,y)第b波段的反射率值R(x,y,b)，表示為：

R(x,y,b)＝R_l(x,y,b)×W_l(x,y,b)+R_s(x,y,b)×W_s(x,y,b) (1)

其中，R_l(x,y,b)表示基于多年時(shí)間序列反射率圖像的估計(jì)值，作為長(zhǎng)時(shí)間估計(jì)值，R_S(x,y,b)表示基于單幅參考圖像的估計(jì)值，作為短時(shí)間估計(jì)值，W_l(x,y,b)和W_s(x,y,b)分別是這兩種估計(jì)值的權(quán)重，并且它們的和為1；S3.計(jì)算長(zhǎng)時(shí)間估計(jì)值R_l(x,y,b)：

基于Landsat反射率時(shí)間序列圖像數(shù)據(jù)集的自相關(guān)，長(zhǎng)時(shí)間估計(jì)值R_l(x,y,b)表示為：

其中，m是在云像元(x,y)位置處，多年的Landsat圖像數(shù)據(jù)集中獲得的無(wú)云清晰Landsat圖像總數(shù)，對(duì)于不同的像元，對(duì)應(yīng)的m數(shù)值不同，R_i(x,y,b)表示云像元(x,y)位置處第i個(gè)無(wú)云清晰圖像第b波段的反射率值，a_i(x,y,b)和a₀(x,y,b)分別表示線性回歸方程的斜率和截距，其中i＝1，2...m；如果云像元(x,y)與其鄰域內(nèi)相鄰像元之間具有相似的反射率時(shí)間序列變化軌跡，那么它們之間有相同的回歸參數(shù)；通過(guò)以下步驟確定云像元(x,y)與其鄰域內(nèi)相鄰像元之間的相似性：

A.利用非監(jiān)督分類器ISODATA對(duì)云像元(x,y)及其鄰域內(nèi)的所有像元進(jìn)行自動(dòng)聚類；B.計(jì)算與云像元(x,y)屬于同一類的相鄰像元的反射率時(shí)間序列之間的絕對(duì)差(DIFF)和線性相關(guān)(COR)，分別表示為：

其中，time_series_(x,y,b)＝[R₁(x,y,b),R₂(x,y,b),…,R_m(x,y,b)]，表示云像元(x,y)位置處m個(gè)無(wú)云清晰圖像的順序集合；R_i(x_j,y_j,b)表示第j個(gè)相鄰像元(x_j,y_j)位置處第i個(gè)無(wú)云清晰圖像第b波段的反射率值，表示第j個(gè)相鄰像元(x_j,y_j)位置處m個(gè)無(wú)云清晰圖像的順序集合，其中j＝1,2,…n，n是與云像元(x,y)屬于同一個(gè)類別的相鄰像元(x_j,y_j)的總數(shù)；

由于中有些值可能會(huì)被云污染，用線性插值填補(bǔ)被云污染的值；通過(guò)將絕對(duì)差DIFF_{j_(x,y)}和線性相關(guān)COR_{j_(x,y,b)}通過(guò)公式(5)組合在一起，得到新定義的相似性值S_{j_(x,y,b)}：

考慮這n個(gè)相鄰像元之間的相似性，通過(guò)最小化公式(6)的目標(biāo)函數(shù)來(lái)估計(jì)公式(2)中的回歸參數(shù)：

其中，R(x_j,y_j,b)表示第j個(gè)相鄰像元(x_j,y_j)第b波段的反射率值；采用最小二乘法求解式(6)，求出a₀(x,y,b),a₁(x,y,b),…,a_m(x,y,b)后代入式(2)中便可以得到長(zhǎng)時(shí)間估計(jì)值R_l(x,y,b)；

S4.計(jì)算短時(shí)間估計(jì)值R_S(x,y,b)：

選取沒(méi)有云污染的單幅參考圖像R_r，其短時(shí)間估計(jì)值表示為：

R_s(x,y,b)＝α(x,y,b)×R_r(x,y,b)+β(x,y,b) (7)

其中，R_r(x,y,b)表示參考圖像上云像元(x,y)位置處第b波段的反射率值；α(x,y,b)和β(x,y,b)分別表示線性回歸的斜率和截距，分別通過(guò)最小化以下目標(biāo)函數(shù)得到:

其中，(x_is,y_is)為云像元(x,y)位置處其鄰域內(nèi)的第is個(gè)相似像元，is＝1,2,…h(huán)，R(x_is,y_is,b)和R_r(x_is,y_is,b)分別表示目標(biāo)圖像和參考圖像中相似像元(x_is,y_is)位置處第b波段的反射率值；w_is為每個(gè)相似像元的權(quán)重，相似值越高其權(quán)重越大，w_is由空間距離D_is和光譜差值S_is共同確定：

其中，g是波段數(shù)目，b＝1,2...g，x、y分別為云像元(x,y)位置處的橫坐標(biāo)和縱坐標(biāo)，x_is、y_is分別為相似像元(x_is，y_is)的橫坐標(biāo)和縱坐標(biāo)；將D_is和S_is歸一化后，由以下公式得到相似像元的權(quán)重值w_is：

因此，通過(guò)用最小二乘法求解式(8)得到α(x,y,b)和β(x,y,b)，并將其代入式(7)后便能得到短時(shí)間估計(jì)值R_S(x,y,b)；

S5.計(jì)算長(zhǎng)時(shí)間估計(jì)值R_l(x,y,b)和短時(shí)間估計(jì)值R_S(x,y,b)的權(quán)重：

長(zhǎng)時(shí)間估計(jì)值R_l(x,y,b)和短時(shí)間估計(jì)值R_S(x,y,b)的權(quán)重根據(jù)預(yù)測(cè)誤差確定；目標(biāo)圖像云覆蓋區(qū)域的真實(shí)反射率值是未知的；因此，采用短時(shí)間估計(jì)方法中所使用的單幅參考圖像R_r計(jì)算預(yù)測(cè)誤差；為了避免在以下的解釋中出現(xiàn)混淆，使用縮寫(xiě)CR_T和NCR_T分別表示在目標(biāo)圖像中云覆蓋的區(qū)域及其鄰域，用CR_R和NCR_R分別表示在參考圖像中對(duì)應(yīng)的區(qū)域；

采用步驟S3中的長(zhǎng)時(shí)間估計(jì)方法計(jì)算CR_R中云像元(x,y)位置處第b波段反射率的長(zhǎng)時(shí)間估計(jì)值R_{r_l}(x,y,b)和對(duì)應(yīng)的預(yù)測(cè)誤差ε_l(x,y,b)：

ε_l(x,y,b)＝|R_{r_l}(x,y,b)-R_r(x,y,b)| (11)

短時(shí)間估計(jì)值的預(yù)測(cè)誤差計(jì)算步驟為：

隨機(jī)選取NCR_T中所有無(wú)云像元的一半作為基礎(chǔ)輔助值，假設(shè)選擇的一半像元的數(shù)目為p，對(duì)NCR_R中對(duì)應(yīng)位置的像元(x_inei,y_inei)計(jì)算其短時(shí)間估計(jì)值R_{r_s}(x_inei,y_inei,b)，inei＝1,2...p，對(duì)這些像元的預(yù)測(cè)誤差進(jìn)行加權(quán)平均間接得到短時(shí)間估計(jì)的預(yù)測(cè)誤差值ε_s(x,y,b)：

其中，R_r(x_inei,y_inei,b)表示NCR_R中像元(x_inei,y_inei)位置處第b波段反射率的真實(shí)值，W_inei表示像元(x_inei,y_inei)的權(quán)重值，采用同公式(10)一樣的方式通過(guò)空間距離和光譜距離計(jì)算得到W_inei；

長(zhǎng)時(shí)間估計(jì)值和短時(shí)間估計(jì)值的權(quán)重W_l(x,y,b)和W_s(x,y,b)由以下公式得到：

W_l(x,y,b)＝[1/ε_l(x,y,b)]/[1/ε_l(x,y,b)+1/ε_s(x,y,b)]

S6.將長(zhǎng)時(shí)間估計(jì)值、短時(shí)間估計(jì)值及其對(duì)應(yīng)的權(quán)重值代入公式(1)中，最終得到目標(biāo)圖像中云像元(x,y)第b波段的反射率值R(x,y,b)。