專利名稱：遙感影像的處理方法及系統的制作方法
技術領域：
本發明涉及遙感影像處理的技術領域，特別是指通過使用多于一 幅遙感影像生成新的遙感影像的處理方法及系統。
背景技術：
綜合應用多傳感器遙感數據在空間、時間、光譜分辨率上的不同
優勢的技術，即多源遙感數據融合技術從20世紀80年代開始發展， 目的是通過對多種傳感器數據的處理，提高遙感數據的空間和光譜 分辨率。目前較為成熟的融合方法有IHS變換、主成分變換、小波 變換以及Brovey變換、高通濾波等，但是這些方法主要用于融合同 一或相近時刻的高空間分辨率全色波段影像和較低空間分辨率的多 光鐠波段影像，其結果圖像為同 一時刻的高空間分辨率多光譜影像。 由于高空間分辨率的影像并不能和所有低空間分辨率影像同時獲 取，這些傳統的多源數據融合方法并不能同時提高遙感數據的空間 和時間分辨率。即，無法生成高時空分辨率的反射率數據。為此， Gao等人在2006年提出了 STARFM方法，該方法能夠得到較為準確 的高時空分辨率反射率影像，但其生成的影像的細節保持的不好。
發明內容
有鑒于此，本發明的主要目的在于提供一種遙感影像的處理方 法，以實現能生成將細節保持的更好的影像。
為解決上述技術問題，本發明提供了一種遙感影像的處理方法， 設定某時刻的所求影像AO的一像元為目標像元，確定目標像元對應 的地表區域為目標區域，獲取同一地表區域的具有高空間分辨率的 第一組影像和具有低空間分辨率的第二組影像，從所述兩組影像中 對應的選擇至少兩對不同時刻的影像，按以下步驟分別計算各目標 像元的反射率以生成AO影像a，根據所選的影像計算所述目標區 域在兩組影像中的反射率變化量之間的轉換系數；b,從第二組影像
6中選擇與所述影像A 0時刻對應的影像B 0和至少 一 個其他時刻的影 像；根據所選擇的影像計算所述目標區域在各所述其他時刻的影像 中的反射率相對于其在影像BO中的反射率的變化量；c，根據所述 轉換系數將所述各反射率的變化量轉換為在第 一組影像中的各反射 率變化量；和d,從第一組影像中選擇與所述其它時刻對應的影像， 根據轉換得到的各反射率變化量和目標區域在所選影像中的反射率 計算所述目標區域的反射率。
由至少兩對時刻相對應的影像就可以得到時刻對應的兩組影像 中地表區域反射率的變化信息之間的關系，即轉換系數。根據所挑 的影像經由步驟a得到目標區域的轉換系數。因為變化信息之間的 關系在較短時間范圍內可以認為是穩定的。所以只要所選的各影像 的時刻和影像AO時刻相差不大，則可以根據轉換系數經由步驟c將 步驟b所得的目標區域在第二組影像中的反射率變化量(簡稱第二 變化量)轉換為目標區域在第一組影像中的反射率變化量(簡稱第 一變化量)。由此，由步驟d所得的反射率是基于更精細的第一組 影像中的反射率和轉換后的變化量生成的。這樣生成的AO影像就很 好的保持了第一組影像中的細節，避免所生成的影像在視覺上有被 平滑的效果。
優選的是，所述步驟a包括al，將所選的第一組影像中的與所 述目標區域對應的像元i殳定為中心像元；a2， 乂人所述中心像元各自 所在影像中分別篩選出至少 一 個與所述中心像元屬于同類地物的相 似^f象元；和a3,確定各相似4象元所對應的地表區域為相似區域，根 據所述相似區域在所選影像中的反射率計算所述轉換系數。
混合像元反射率的變化量為各端元反射率的變化量的綜合反映。 所以端元可以提供比混合像元更為準確的變化信息。又因為作為端 元的相似像元與中心像元屬于同類地物，所以，相似像元能夠提供 最為準確的變化信息，由此可以得到更為準確的轉換系數。其中， 中心像元也算是相似像元之一，在未能篩選出其它相似像元時，中 心像元是唯一的相似4象元。尤其是對于目標區域中存在石皮碎斑塊和小地物的地表區域，其影像在第二組影像中為混合像元。其提供的 變化信息(第二變化量)被由各相似像元得到的轉換系數轉換修正 成與目標像元屬于同類地物的端元的變化信息(第一變化量)。由 此，目標像元的反射率準確度更高了 ，尤其是對破碎斑塊和小地物 具有更好的融合效果。
優選的是，所述步驟a3之前還包括a21,根據所選的影像計算 各相似像元的權重；所述步驟b包括bl,根據BO影像和所述第二 組的其他時刻的影像計算各所述相似區域在第二組影像中的反射率 變化量；b2，根據得到的所述權重將步驟bl得到的各反射率變化量 加權平均計算所述目標區域在各所述其他時刻的影像中的反射率相 對于其在影像B0中的反射率的變化量。
目標區域的第二變化量不是直接由目標區域在B0影像和所述第 二組的其他時刻的影像的反射率求差得到，而是由各相似像元按照 各自權重共同來提供。由此修正了 BO影像和第二組的其他時刻的影 像中目標區域對應的混合像元所提供的變化信息。從而得到更準確 的第二變化量。
優選的是，所述步驟a3中所述轉換系數的計算是將各所述相似 區域在所選影像中的反射率以在第二組影像中的反射率的變化量為 自變量，在第一組影像中的反射率的變化量為因變量進行回歸分析。
利用回歸分析可以消除因噪聲等原因引起的誤差，從而得到更穩 健的變換系數，減少可能的誤差。
優選的是，所述步驟a3之前還包括a22，根據所選的影像計算 各相似像元的權重；所述回歸分析是根據各相似像元的權重進行的 力口權回歸分析。
各相似像元距離其中心像元的空間距離，以及其的光譜特征和中 心像元的光譜特征的一致性等原因都導致了各相似像元提供變化信 息時，存在貢獻大小的差異。在回歸分析時，利用各相似像元的權 重進行加權回歸分析，可以體現出這些差異。進而得到更準確的轉 換系數，以進一步提高目標像元的反射率準確度。
8優選的是，所述各相似像元的權重的計算包括以下步驟a221， 根據所選影像計算各相似像元所有波段光譜向量和第二組影像中各 所述相似區域的像元所有波段光譜向量的一致程度；a222,根據所 述一致程度計算各相似像元的權重。
所述第一組影像的像元為高空間分辨率像元，所述第二組影像的 像元為低空間分辨率像元，高空間分辨率像元對應地表區域所屬地 物在對應的低空間分辨率像元對應地表區域的各種地物中所占比例 即為高空間分辨率像元的純度。所以相似像元的純度越大，也就意 味著其所提供的變化信息越準確。從而通過步驟a221和步驟a222 可以提高能提供更準確變化信息的純度大的相似像元的權重。而純 度大小的計算則是由所有波段光譜向量的 一 致程度計算得到。
優選的是，所述步驟a2還包括將各篩選出的相似像元所對應 的地表區域取交集，以取交集后的地表區域所對應的像元為相似像 元。
有的地物隨時間會發生變化，其光譜特征也相應發生變化，只利 用 一組時刻對應的影像篩選出的相似像元可能會因此篩選出已和中 心像元不屬于同一地物的相似像元。所以取交集后可以篩選出的各 時刻的均和中心像元都是同 一地物的相似像元。從而得到更可靠的 目標像元。
本發明還相應的提供了 一種遙感影像的處理系統，根據所述兩組 影像求得高空間分辨率的A 0影像，所述A 0影像不屬于第 一 組影像， 并由至少一個待計算的像元組成，所述處理系統包括影像選擇模 塊101 ,用于獲取同一地表區域的具有高空間分辨率的第一組影像 和具有低空間分辨率的第二組影像，從所述兩組影像中對應的選擇 至少兩對不同時刻的影像；目標像元設定模塊102 ,用于設定某時 刻的所求影像AO的一像元為目標像元，確定目標像元對應的地表區 域為目標區域；目標像元反射率計算模塊104 ,用于計算目標像元 的反射率；和影像生成模塊108 ,用于根據目標像元反射率計算模 塊104計算得到的各目標像元的反射率生成AO影像；所述目標像元反射率計算模塊104包括以下子模塊轉換系數模塊103 ，根 據所選的影像計算所述目標區域在兩組影像中的反射率變化量之間 的轉換系數；第二變化量計算模塊105 ,用于從第二組影像中選擇 與所述影像AO時刻對應的影像BO和至少一個其他時刻的影像，根 據所選擇的影像計算所述目標區域在各所述其他時刻的影像中的反 射率相對于其在影像BO中的反射率的變化量；變化量轉換模塊 106 ,用于根據所述轉換系數將所述各第二變化量計算模塊105得 到的各變化量轉換為在第一組影像中的各反射率變化量；和反射率 計算模塊107 ，用于從第一組影像中選擇與所述其它時刻對應的影
像，根據轉換得到的各反射率變化量和目標區域在所選影像中的反 射率計算所述目標區域的反射率。
優選的是，所述轉換系數模塊103包括中心像元設定模塊 1031 ，用于將所選的第一組影像中的與所述目標區域對應的像元設 定為中心像元；相似像元篩選模塊1032 ,用于從所述中心像元各 自所在影像中分別篩選出至少 一個與所述中心像元屬于同類地物的 相似像元；和轉換系數計算模塊1033 ，用于確定各相似像元所對 應的地表區域為相似區域，根據所述相似區域在所選影像中的反射 率計算所述轉換系數。
優選的是，所述目標像元反射率計算模塊104還包括權重模 塊，用于根據所選的影像計算各相似像元的權重；所述第二變化量 計算模塊105包括相似像元第二變化量計算模塊1051 ,用于從 第二組影像中選擇與所述某時刻對應的影像BO和至少一個其他時 刻的影像，根據所選擇的影像計算各所述相似區域在第二組影像中 的反射率變化量；第二變化量加權平均模塊1052 ,根據得到的所 述權重將相似像元第二變化量計算模塊1051得到的各反射率變化 量加權平均計算所述目標區域在各所述其他時刻的影像中的反射率 相對于其在影像BO中的反射率的變化量。
圖1為遙感影像的處理方法的實施例流程圖；圖2為篩選相似像元的流程圖3為取交集篩選相似像元示意圖4為本實施例中計算相似像元的權重的流程圖5為得到轉換系數的流程圖6為轉換系數的計算示意圖7為計算預測時刻的目標像元的高空間反射率的流程圖8為處理系統的結構圖9為轉換系數模塊103結構圖10為第二變化量計算模塊105結構圖11為一個實施例中權重模塊109結構圖12為另一個實施例中權重模塊109結構圖13為相似像元篩選模塊1032結構圖14為調整模塊1034結構圖
圖15為模擬物候變化的情況；
圖16為模擬線狀地物的情況；
圖17為MDCM和STARFM對圖16(b)中線狀地物的反射率的 預測結果；
圖18為模擬小面積地物的情況；
圖19為MDCM和STARFM對圖18(b)中各圓形小地物的預測 誤差；
圖20為MODIS上排和Landsat下排假彩色合成影像；
圖21為分別由STARFM和MDCM方法生成的7月11日Landsat
分辨率的假彩色合成影像；
圖22為真實影像與重建影像的各波段散點圖。
具體實施例方式
在對實施方式進行說明之前，先對本發明的理論基礎進行說明。 對于同一區域來自不同傳感器的遙感數據經過輻射定標、幾何配 準和大氣糾正后，這些數據之間具有一定的可比性和相關性。但是， 由于不同傳感器在波段寬度、光譜響應函數、獲取時刻的大氣狀況等方面存在差異，使得這些多源數據間存在一定的系統偏差。在對 系統誤差進行糾正的同時，如何利用它們之間的相關性實現多源數 據的融合構成了本發明的核心。以下分兩種情況介紹本發明理論基 礎。
對于均一地表情況假設低空間分辨率、高時間分辨率(略寫為 低空間分辨率)和高空間分辨率、低時間分辨率(略寫為高空間分 辨率)傳感器有類似的光譜波段設置。當低空間分辨率影像重采樣 為與高空間分辨率影像相同的空間分辨率(即相同的像元大小)和坐 標系統時，對于同一且均一的地物而言，低空間分辨率像元內僅為 一種地物類型純像元，此時在波段S上高空間分辨率反射率與低空 間分辨率反射率之間的關系即為高、低空間分辨率影像上反射率之 間的相對輻射定標關系，該關系僅由兩種傳感器特性差異(波段寬度 及光譜響應函數)和成像時刻大氣狀況差異所決定， 一般可視為以下 線性關系
,力,用="x C(x,.,力,^， 5) + 6 i
其中，F、 C分別代表高、低空間分辨率的像元反射率，a,，力) 是像元位置，S為光譜波段，？k代表了影像獲取時刻(以日為單位)， "、6分別代表高空間分辨率和低空間分辨率影像相對輻射定標的增 益和偏離值。對于同一且均一的地物而言，l式所示關系在任何影像 獲耳又時刻都成立。
如果已經獲取了 Z。時刻波段B的高、低空間分辨率的影像和/p 時刻同一波段的低空間分辨率影像，那么/p時刻該波段高空間分辨 率的像元反射率可由式2來表示~'，乃^，5) = """，^ )" 2
而/。時刻高、低空間分辨率影像的像元反射率關系為
F(、 ， ^, ^ ， B) = a x C(x,,力，~, B) + 6 3
由式2和式3相減并移項整理可得到式4: ^",力,^，5) = 7^,，3;乂0,5) + "><((:"，力,^,￡)一(:",3;-。,5)) 4
上式可以理解為/p時刻高空間分辨率的像元反射率由？o時刻高 空間分辨率的反射率加上 相對？。時刻的反射率變化量。式中fl代表高空間分辨率和低空間分辨率影像相對輻射定標的增益，也可理 解為高空間分辨率和低空間分辨率像元反射率之間的轉換系數，該 系數由兩種傳感器特性差異和成像時刻大氣狀況差異所決定。如果 兩種傳感器成像時刻相近，可視為大氣狀況相同，則該系數僅由兩 種傳感器特性差異所決定，且不隨時間發生變化。
如果獲取任意兩個時刻Zm、 /n的兩對高、低空間分辨率影像，那 么通過對這兩個時刻的高、低空間分辨率圖象上的各波段純像元反
射率進行線性回歸即可求出各波段的轉換系數a,通過該系數和式4
即可由 一 系列時間連續的低空間分辨率影像生成一 系列時間連續的 高空間分辨率的反射率影像。
對于非均一地表情況由于傳感器空間分辨率的限制以及地物的
復雜多樣性，混合像元普遍存在于低空間分辨率影像上，對地表地
物類型多樣、分布雜亂的區域尤其如此。此時，低空間分辨率像元 與其對應的高空間分辨率像元的反射率并不是對同 一 種地物的反
應，二者之間的差異不僅僅是輻射定標的差異，此時基于均一地表
純像元的式4并不適合非均一地表混合像元的情況。
混合像元的光譜特征是其內部各類地物光譜端元光譜的綜合反 映，所以兩個時刻/m、 ^之間其反射率的變化為各類端元光譜變化的 綜合反映。目前混合像元光譜混合模型最常用的是線性光譜混合模 型，即利用 一 個線性關系來表達混合像元和各端元光譜之間的關系。
假設低空間分辨率影像上fm、 ？n時刻的某波段混合像元反射率分別為
4、 K,根據線性光譜混合模型，IV K可由下式表達
其中，M為該混合像元包含的端元數目，《M和；^分別表示^、
^時刻第i種端元的光譜，S為誤差項，且考慮S在Zm、 Zn時刻不變。
由式5可得到變化量其中下式中的力是第/個端元的比例；
13,=i 6
假設/m、 ^的時間間隔在一定范圍內時，各端元反射率的變化可
以近似為線性變化，則端元光譜《 可由式7表示其中下式中的n' 是第/個端元反射率的線性變化率；X, =。xA"X,.w 7
其中，Vm,將式7代入式6可得到下式^="￡〃 8 已知第A:端元/m、/n時刻的反射率，則由式7可求出△〖& = ^"_義" 將式9代入式8，可得%^ = + 10
由于假設Zm、 ^的時間間隔在一定范圍內，各端元比例/和各端 元反射率的線性變化率r在,m、 Zn時刻之間可視為穩定，則式10的 右邊為一穩定變量，用V表示，其意義為A類地物端元反射率變化 量與混合像元反射率變化量的比值變化比例系數。
對于非均一地表，低空間分辨率像元為混合像元，其對應的高空
間分辨率影像上的各像元視為端元。根據式10， Zm、 Z。時刻的高空間 分辨率像元反射率尸(JC,.，》，^， B)和F(X,.， 乂， ^， 5)與對應的低空
間分辨率像元反射率C(x,.， 乂.， l 5)和C(x,.，為，"，B)滿足以下關系..
通過式11可知，/m、 ^時刻的各波段高空間分辨率反射率與對應 的低空間分辨率反射率為線性關系，進行線性回歸即可求出各波段 的變化比例系數v。
如果在fm到rn內，已經獲取了 /。時刻波段B的高、低空間分辨 率的影像和fp時刻同一波段的低空間分辨率影像。根據以上的假設， ^相比"時刻該波段高空間分辨率的像元與低空間分辨率像元的反
射率變化量之間的比例系數仍為v則c^，乂，/p，旬-c"，乂，/。，5) 12 通過式12即可求解出^時刻的高空間分辨率反射率F"，力，,p,5卜F(x";^。,萬)+ vx(C(x"力，^5) —C"，力J。,5)) 13
雖然公式13與公式4的形式相同，但其意義和適用情況不同。 式4表示的是純像元不同分辨率反射率之間相對定標的關系，在任 何時刻都成立，其計算的結果最為準確。而式13表示的是混合像元 中不同分辨率反射率變化量之間的比例關系，根據假設只在^和& 之間或其附近時刻才成立，而且是基于短時間范圍內各地物反射率 的變化為線'f生的假設，所以其計算的結果是近似解。
基于上述的理論基礎，當低空間分辨率像元為純像元時，可以利 用式4求解出^時刻的高空間分辨率反射率，而當低空間分辨率像 元為混合像元時利用式13求解。但純像元計算的結果較混合像元準 確可靠。為盡量利用純像元的信息，以及減小影像受到云、大氣等 污染給計算結果帶來的不確定性，本實施例中采取窗口內求解的辦 法。即，以像元XW2， ， ^W2)為中心，取寬度為W的窗口，先篩選出
與中心像元屬同類地物的像元相似像元，然后利用式14計算中心像 元^時刻的高空間分辨率反射率
F(Uw2，^5",( 2,尺/"^) + iXxrX(C(H^S)-C(《，乂，^鄰
m 14 其中，iV為窗口內與中心4象元屬于同類地物的相似^象元數目，『 為各相似像元的權重，F為高、低空間分辨率的轉換系數。式14的 意義是/p時刻的反射率由已知時刻的反射率加上^時刻相對已知/0 時刻的變化量決定，而該變化量由窗口內的相似像元共同預測。
在本實施例之前，通常需要先對獲取的高、低空間分辨率影像進 行預處理，使得多源數據之間具有相似的波段設置，相同的像元大 小和圖幅尺寸以及相同的坐標系統。采用常用的軟件和方法即可實 現。下文將結合附圖對遙感影像的處理方法的實施方式進行詳細說 明。
本實施方式中遙感影像為對同一地表攝像而得到的兩組影像，其 中，第一組影像具有高空間分辨率，第二組影像具有低空間分辨率， 根據所述兩組影像求得tp時刻的高空間分辨率的A影像，所述A影 像由至少一個待計算的像元組成。如圖1所示，遙感影像的處理方法包括以下步驟步驟100，從所述遙感影像中挑選至少兩對時刻
對應的影像。本實施例中，選用的是已知的tm、 tn時刻的兩對影像。
從第二組影像中挑選其時刻與tp時刻所對應影像；步驟101，設定 一個待計算的像元為目標像元，確定目標像元的應的地表區域為目 標區域；步驟102，以目標區域對應的像元為中心像元，篩選出已知 時刻的高空間分辨率影像中所述中心像元的窗口內的同類地物像元 得到相似像元，其中，中心像元也算是相似像元之一，在未能篩選 出其它相似像元時，中心像元是唯一的相似像元；步驟104，計算所 述相似像元的權重。本實施例中，根據/m、 /w時刻的兩幅低空間分辨 率影像，計算出步驟102中篩選出的相似像元的權重；步驟106,根 據所述權重采用加權算法對已知時刻的高、低空間分辨率反射率進 行線性回歸計算得到轉換系數；步驟108,計算預測時刻的目標像元 的高空間分辨率反射率，本實施例中，根據轉換系數和預測時刻tp 的低空間分辨率影像計算已知時刻高空間分辨率影像的時間權重， 并根據已知時刻高空間分辨率影像按照所述時間權重加權計算得到 預測時刻的目標像元的高空間分辨率反射率；步驟IIO,逐次計算所 有目標像元并最終生成A影像。下文將結合附圖對上述各個步驟進
行詳細說明o
圖2為本實施例中篩選相似像元的流程圖，如圖所示包括 步驟1022,以目標區域對應的像元為中心4象元，篩選出各已知 時刻的高空間分辨率影像中所述中心像元的窗口內的與目標像元屬 于同類地物的像元得到相似像元。在目標像元( )為中心
的窗口內，與要計算的中心像元屬于同類地物的高空間分辨率像元 才能提供比較正確可靠的反射率變化信息。可以使用以下兩種方法
篩選出相似像元(a)對高空間分辨率影像進行非監督分類，和中 心像元類型相同的像元為相似像元；(b)4吏用閾值判斷，以窗口內 其它像元反射率與中心像元反射率之差來判斷是否屬于同類地物。 這兩種篩選相似像元的方法具有共同的意義，即是以光譜相似的像 元為同一類地物。但是又有不同之處閾值法是在窗口內尋找與中心像元的光譜相似的像元，中心像元作為了搜尋的中心，而且其篩
選條件是局部適用的，隨著中心像元的位置變化而變；而分類是在 全影像進行，如果分類有誤差，那么判斷出的相似像元可能具有很 大的錯誤。所以在本實施例中，使用閾值法來篩選相似像元。在本 實施例中，判斷窗口內某像元所有波段都滿足式15的條件，則該像 元被確 定為中 心像元的相似像元 。
其中，o(S)為&時刻高空間分辨率影像所有像元S波段的標準差， m為預估的地物類別數，比如影像覆蓋的區域主要有植被、棵地、 水體、巖石4類地物，則可將m值設為4。當然，也可以用所有波 段中的部分波段進行篩選或者采用方法(a)進行篩選。
步驟1024，將篩選出的兩組已知時刻的相似像元取交集。有的 地物隨時間會發生變化，其光譜特征也相應發生變化，只利用一個 時相的影像篩選出的相似像元可能會出現錯誤。比如窗口內有棵地 和農作物兩種地物類型，如果中心像元是農作物，^時刻農作物并 未生長，其光譜特征與棵地相同，此時篩選出來的相似像元是棵地， 如果"時刻農作物已生長，則在"時刻的影像上能篩選出正確的農 作物像元。鑒于此，本實施例中分別利用fM、 /"時刻兩幅高空間分辨 率影像篩選相似像元，并將篩選結果取交集，以此提高篩選結果的 準確性。如圖3所示為取交集篩選相似像元示意圖，"*"表示中心 像元，"O"表示中心像元的相似像元。
圖4為本實施例中計算相似像元權重的流程，如圖所示包括
步驟10402,計算相似像元的純度大小。權重『決定了窗口內篩 選出的各相似像元對計算目標像元反射率的貢獻。在本實施例中權 重『由各相似像元的純度大小和離中心像元的空間距離遠近共同決 定。純度越高、空間距離越近的相似像元權重越大，這是因為均一 地表的純像元能提供最準確的變化信息，而且認為空間上離中心像 元越近的相似像元的變化情況和中心像元越 一 致。
高空間分辨率影像上的某像元代表了某種地物類型，而其對應位
17置低空間分辨率影像的像元包含了該類地物，純度即是指低空間分 辨率 <象元中該種地物所占比重。如果在均 一 地表情況下低空間分辨
率像元全由該種地物覆蓋，則認為是純像元。比如Landsat影像上某
像元所覆蓋的30mx 30m的區域為小麥，其對應位置的MODIS影像
上像元所覆蓋的500mx 500m的區域也全為小麥，則認為該像元為
小麥純像元。由于不同的地物具有不同的光譜特征，所以可以通過
比較高空間分辨率像元和其對應位置的低空間分辨率像元的光譜曲
線的一致程度來判斷純度大小。如果高空間分辨率像元和其對應的
低空間分辨率像元的光譜一致程度越高，則認為純度越大。而高、
低空間分辨率像元的光譜曲線之間的 一致性可由相關系數來刻畫，
所以篩選出的各相似像元的純度i 可由式16計算 及=譜(F,，C,)
'一V^y,(C,) 16
F, ={F(x,,_y,，,m,JB1),...,F(x,.,;;,,,m,JBJ,F(x,j,,, ,JB1),...,JF(x,，j;,， ,5 )} 其中向量Fi， Ci分別表示第/個相似像元高、低空間分辨率的光
鐠向量，cov為求協方差，D為求方差。i 值的范圍為-l到l, iM直
越大表示純度越高。之所以將G、 4時刻的光譜曲線放在一起計算相 關系數，是考慮到地物隨物候的變化，光譜特征也會發生變化，只 使用一個時相的光譜可能會出現錯誤。比如小麥在返青之前其光譜 特征和棵地一致，如果某小麥高空間分辨率像元對應的低空間分辨 率像元除了小麥外還包含棵地，則用小麥返青前的獲取的影像計算 會得到較高的純度，顯然是不正確的。當加入小麥返青后獲取的影 像進行判斷時，就會得到更加正確的純度。由上可以看出，只要獲 得至少兩組已知時刻的高空間分辨率像元的所有波段光譜向量和在 其對應位置同時的低空間分辨率像元的所有波段光譜向量，就可以 根據其一致程度來計算純度大小。
利用純度大小作為計算權重的一個要素，較之目前采用同時刻 高、低空間分辨率反射率之差表示的光譜特性作為計算權重的要素將更加有效的利用相似像元并且得到的結果也更加準確。
步驟10404,判斷相似像元中是否有純像元，有則進入步驟 10406，沒有則進入步驟10408。
步驟10406，設定純像元的權重為1。在本實施例中，權重W的 取值范圍是0到1，所有相似像元的權重之和為1。 R=l的像元為 純像元。為了盡量利用純像元，將純像元的權重設為1。即，中心像 元的反射率變化信息完全由這些純像元提供。在另一個實施例中， 當相似像元中有P個純像元存在即R=l的像元時，規定這些純像元 的權重W為1/P，其他相似像元的權重為0。當然，純像元的判斷也 可以是在計算權重過程中的其它步驟中實現。例如在步驟10402之 前或者在步驟10408之后。甚至不判斷是否有純像元也可以。
步驟10408,計算相似像元距離中心像元的空間距離。相似像元 ，")與窗口中心像元;w2, ， JV2)的空間距離^由式17計算
《=1 + V"W2 - A )2 + - X )2 2) "
上式包含了對空間距離的歸一化，在窗口 w內各相似像元的空間距 離《的取值范圍為1到l+2G's,值越大，該相似像元離中心像元越遠。 步驟10410，計算相似像元的權重。接下來將純度和空間距離這 兩個不相關的量聯系在一起，來確定各相似像元的權重『。上文已 經提到純度越高、距離越近的像元能提供更準確的信息，即A越大， d越小，其權重『越大。本實施例中，通過公式18，將純度和空間 距離結合為總距離D: A=(l-^>《 18
Z)越大的相似像元為目標像元反射率的計算貢獻越小，再將D的倒
『'=(i/A)/i](i/A)
數歸一化得到各相似像元的權重『式19: '=1 19
當然，建立純度、空間和權重的關系可以采用多種形式，只要能
實現純度越大，距中心像元越近的相似像元的權重越大即可。 圖5為得到轉換系數的流程，如圖所示包括 步驟1062,計算轉換系數。由公式4和公式13可知，轉換系數
K包含了相對定標的增益a和高低空間分辨率變化量的比例系數v,
19二者都可由已經獲取的fw、"時刻高、低空間分辨率各波段反射率進 行線性回歸得到。為了減小噪聲污染等對轉換系數計算帶來的不確 定性，將所有相似像元^、 ^時刻高、低空間分辨率反射率進行線性
回歸，回歸直線的斜率即是轉換系數F。為了利用更可靠的相似像元 的信息計算轉換系數，采用加權最小二乘法進行線性回歸，權重為 各相似像元的『。如圖6所示為轉換系數的計算示意圖。某窗口內 有12個像元和中心像元屬于同類地物包括中心像元本身，虛線框分 別標注出"、4時刻的反射率，可見從^到,"時刻，高空間分辨率 反射率和其對應的低空間分辨率反射率都有所增加，但是高空間分
辨率反射率增加幅度較大，是低空間分辨率的6.448倍^=0.915, 尸<0.001。這說明該類地物該波段反射率在/w、"時刻之間發生了較 大的變化，比如植被生長導致的近紅外反射率急劇增加。根據理論 基礎的論述，各相似像元都反映了低、高空間分辨率反射率之間的 轉換系數，但由于噪聲等影響，通常不能只利用某一個相似像元來 計算V，這樣的結果帶有很大的不確定性，而線性回歸能利用所有 相似像元的信息得到 一個比較穩健的轉換系數。又因為各相似像元 存在貢獻大小的差異，為了體現這些差異，在線性回歸之中進一步 還考慮各樣本的權重，本實施例采用的是加權最小二乘法。當然， 也可以釆用其它的回歸方法。
步驟1064，根據不確定性調整轉換系數。由于幾何配準、大氣 修正等數據處理過程都給反射率帶來一 定的誤差，如果低空間分辨 率像元^、"時刻之間的反射率變化太小，在誤差范圍內，則進行加 權回歸得到的斜率J/具有很大的不確定性。這可能是以下兩種情況 低空間分辨率像元內反射率發生變化的地物所占面積比例很小，在 整個低空間分辨率像元尺度上體現不出這種變化；也可能是低空間 分辨率像元內有的地物反射率增加，有的地物反射率減少，而增加 和減少的量相等，使得在低空間分辨率像元尺度上體現不出變化來。 為解決這個問題，在本實施例中，在所有相似像元的低空間分辨率 反射率變化均值小于至少兩個已知時刻反射率的變化量的不確定性時將轉換系數設為1。假設反射率各波段的不確定性U為該波段最大 值的1%,而且這種不確定性在各影像各波段間是互相獨立的，那么
兩個已知時刻反射率的變化量的不確定性為u = AX 20
其中，Utm是Um分別是^、"時刻反射率的不確定性。如果所有 相似像元的低空間分辨率反射率變化均值小于式20計算出的不確定 性，則只能將第二變化量平均分配到其內部的每一個高空間分辨率 像元上，即變化量的轉換系數r為l。當然，已知時刻反射率的變 化量的不確定性也可以通過其它統計方法得到。對式20的不同表達 方式，以及其它的變換的計算方式都應包含在本發明的保護范圍之 內。這樣就避免了兩個時刻間反射率的變化太小時，因為反射率自 身帶有一定的不確定性，也就是誤差，而導致計算出不真實的轉換 系數V。
步驟1066，修正奇異的轉換系數。由于影像自身噪聲的存在， 導致絕對值較大的轉換系數出現。為修正個別較奇異的轉換系數， 對全幅影像所有像元的轉換系數進行統計，在本實施例中將均值上 下2倍標準差以外的轉換系數r設定為1。當然，也可以根據需要將 其它預設倍數的標準差以外的轉換系數設定為1。
圖7為計算預測時刻的目標像元的高空間反射率的流程，如圖所 示包括步驟10802，根據預測時刻的低空間分辨率影像計算各相似 像元相對已知時刻的第二變化量；步驟10804,根據所述轉換系數將 所述第二變化量轉換為第一變化量；步驟10806，將各相似像元的第 一變化量按所述權重加權平均后得到目標像元的第一變化量；和， 步驟10808 ，將所述目標像元的第 一變化量與已知時刻的高空間分辨 率反射率相加得到基于不同已知時刻的預測時刻目標像元的高空間 分辨率反射率。通過上述4個步驟完成了式14的計算。以兩組分別 為"、^的已知時刻所獲取的時刻之間或附近的任意^時刻的低空間 分辨率影像，計算出^時刻的高空間分辨率反射率為例。根據式14, 由已知的Zw、"時刻高、低空間分辨率影像即可計算出 時刻的高空間分辨率反射率，此時/w、 ^時刻的高、低空間分辨率反射率都已知， 所以可以分別由/w、 時刻的高空間分辨率反射率計算出Zp時刻的高
空間分辨率反射率Fjcw/2, — 5,將這兩個計算結果分別記
作Fwxw/2, ~， B、 i^W2， JV2， ~, 5。當然，只要已知轉
換系數、預測時刻的低空間分辨率影像和各相似單元的權重就可以 在已知時刻的高空間分辨率影像的基礎上計算得到預測時刻的目標
像元的高空間分辨率反射率。對式14的不同表達方式，以及其它的
變換的計算方式都應包含在本發明的保護范圍之內。
步驟10810,加權計算得到預測時刻的高空間分辨率反射率。為 充分利用G、 /"時刻的高空間分辨率反射率的信息，將兩個計算結果 加權求和得到最后的預測結果。假設/p時刻的低空間分辨率反射率 相對已知時刻變化越小，則其計算的結果越準確可靠。這是因為如 果低空間分辨率影像上沒有體現出反射率的變化，則認為高空間分 辨率影像也不會有反射率的變化。基于以上假設，兩個計算結果的 時間權重T由式21計算
<formula>formula see original document page 22</formula>得到時間權重T后，^時刻的高空間分辨率反射率F:cw/2，少w， ^， S可由下式計算
F(^2,3^2，,p,5^7;xFm(^2,;^2,/p，5) + 7;xF"(^2,:^2^,5) 22
當然，也可以采用其它方法來計算基于不同已知時刻計算得到的
高空間分辨率反射率的權重。例如，如式23，用預測時刻/。和已知
時刻4低空間分辨率反射率之差來計算權重大小。
<formula>formula see original document page 22</formula>
至此得到了最終的預測時刻的目標像元的高空間分辨率反射率。
僅是單獨應用轉換系數配合步驟108就能改善對小地物和線狀地物 的計算精度。而相似像元的篩選以及權重計算的改進則有助于整體
精度的提高。為了更好的對本實施方式更進一步說明和驗證，下面對本實施方 式分別進行模擬數據和真實數據的檢驗。
模擬數據檢驗。為了更好的檢驗方法的精確性和可靠性下面將本
實施例運用于簡單的模擬數據。并且為了便于和GAO的STARFM 方法進行比較，利用GAO文中的模擬數據進行檢驗。在下文中本實 施方式簡稱為MDCM。
物候變化的情況。物候是引起反射率變化的主要原因，尤其在植 被覆蓋率較大的區域。如圖15所示為模擬物候變化的情況。其中 (a)-(c)中的圓形區域的反射率固定為0.05,周圍區域的反射率從O.l 變為0.2再到0.4;(d)-(f)分別由(a)-(c)聚合而成的低空間分辨率影像； (g)和(h)分別是利用MDCM和STARFM對(b)的預測結果。筒單假設 只存在兩類地物，分別為水體和植被，而且形狀固定不變。假設圓 形的水體具有固定的反射率0.05，而水體四周植被的反射率從0.1 圖15(a)增加到0.2圖15(b)再到0.4圖15(c),代表隨著物候變化，植 被的生長導致的近紅外反射率增加。圖15(d)- (f)是將圖15(a)- (c)按 17*17聚合成低空間分辨率的影像(a)- (c)類比于30米分辨率的 landsat, (d)- (f)類比于500米分辨率的MODIS。利用圖15(a)和圖 15(c)-(f)預測圖15(b)。圖15(g)和(h)分別是MDCM方法和STARFM 的結果，圖15(g)和(h)與圖15(b)完全一致，說明兩種方法生成的高 空間分辨率反射率都能準確地反映由于物候引起的反射率變化。物 候變化所引起的反射率變化是真實影像中最普遍存在的現象，說明 MDCM方法和STARFM都能運用于真實影像，生成高時空分辨率反 射率影像。
線狀地物的情況。在土地覆蓋或土地利用類型中，線狀地物相當 普遍，比如道路、橋梁、河流等，而人眼對線狀地物的敏感程度較 高，生成的高分辨率影像中線狀地物的正確與否直接影響影像的視 覺效果。如圖16所示為模擬線狀地物的情況。圖16。 (a)-(c)中的圓 形區域的反射率固定為0.05,直線的反射率固定為0.5,周圍區域的 反射率從O.l變為0.2再到0.4; (d)-(f)分別由(a)-(c)聚合而成的低空
23間分辨率影^象；(g)、 (h)為MDCM和STARFM對(b)中線狀地物反射 率的預測結果。
圖16(a)- (c)在圖15(a)-(c)的基礎上在對角線上增加一線狀地物， 反射率固定為0.5,模擬一條道路陸地之上和橋梁水體之上，圖16(d)-(f)分別由圖16(a)- (c)聚合而成的低空間分辨率影像。利用圖16(b) 之外的5幅影像重新生成圖16(b),圖16(g)顯示了 MDCM和STARFM 重新生成的圖16(b),可見MDCM和STARFM都能定性地預測出該 線狀地物。如圖17所示為MDCM和STARFM對圖16(b)中線狀地 物的反射率的預測結果。從圖17所示兩種方法對圖16(b)的線狀地 物預測的定量結果看，STARFM并不能完全正確地得到線狀地物的 反射率，錯誤主要出現在線狀地物周圍區域發生了變化的一段圓形 地物之外，而MDCM方法能得到完全正確的線狀地物的反射率。
小面積地物的情況。小面積地物在實際影像中也是非常普遍存在 的，特別是在地物斑塊比較破碎的區域，比如小塊農田或小面積水 體等。如果這些小面積地物的尺寸小于了低空間分辨率像元的尺度， 則在低空間分辨率影像中不能清晰觀察到它們的形狀以及得到它們 的反射率信息，而重新生成的高空間分辨率影像希望能準確地反映 這些小地物的形狀和反射率信息。如圖18所示為模擬小面積地物的 情況。(a)-(c)中的圓形區域的反射率從O.l變為0.2再到0.4,各圖中 4個圓形地物的直徑依次為5、 10、 15、 20個高空間分辨率單位，周 圍區域的反射率固定為0.05; (d)-(f)分別由(a)-(c)聚合而成的低空間 分辨率影像；(g)和(h)分別為STARFM和MDCM方法對(b)的預測結果。
假設只存在兩類地物，圖18(a)-(c)中的圓形地物及周圍背景。圖 中4個圓形地物的直徑依次為5、 10、 15、 20的4個小面積地物， 圓形地物的反射率從O.l圖18(a)增加到0.2圖18(b)再到0.4圖18(c), 周圍區域的反射率固定為0.05。圖18(d)-(f)是將圖18(a)-(c)按17*17 聚合成低空間分辨率的影像。圖18(g)和(h)分別是STARFM和 MDCM方法對圖18(b)的重建結果，而如圖19所示為MDCM和STARFM對圖18(b)中各圓形小地物的預測誤差。圖19顯示了二者 對不同直徑小地物的預測平均誤差小地物各像元反射率預測值與真 實值之差的平均值，可見STARFM只有直徑為20的圓形地物的反 射率是正確的，其他小于低空間分辨率像元尺度17的地物的反射率 都比真實值偏小，而MDCM方法對所有尺寸的地物都能生成正確的 反射率。原因主要是STARFM需要純像元提供變化信息，而小于低 空間分辨率像元尺度的小面積地物在低空間分辨率影像上不存在純 像元，所以生成的反射率出現誤差，而MDCM方法如果在窗口內找 不到純像元時，通過這些小面積地物反射率變化量和其所在的低空 間分辨率像元反射率變化量之間的比例系數，將非純像元提供的變 化信息加以了修正，所以能準確地生成這些小面積地物的反射率。
真實數據檢驗。將MDCM方法運用于真實的Landsat-7和MODIS 遙感影像。同模擬數據一樣，為了與Gao等的STARFM進行比較， 使用 Gao 等人文中的遙感數據進行檢驗數據可在 http:〃ledaps.nascom.nasa.gov/ledaps/Tools/StarFM.htm下載。Landsat-7 ETM+數據的空間分辨率為30米，MODIS數據的空間分辨率為500 米，波段為綠光、紅光、近紅外3個波段MODIS的4、 1、 2波段， 分別對應ETM+的2、 3、 4波段。
數據所覆蓋區域的地理位置是54° N， 104° W,大小為12kmx 12km,主要地表覆蓋類型為植被、水體、棵地、道路等，植被的生 長季較短，物候變化較快。影像共有3組，獲取時刻分別為2001年 5月24日、2001年7月11日、2001年8月12日。如圖20為MODIS 上排和Landsat下排假彩色合成影像近紅-紅-綠波段對應紅-綠-藍彩 色合成。從左到右獲取時刻依次是2001年5月24日、2001年7月 11日、2001年8月12日。上面一排為MODIS影像，下面一排為 Landsat影像。利用5月24日和8月12日的兩組影像和7月11曰 的MODIS影像生成7月11日的Landsat分辨率的影像，并將結果與 真實的7月11日Landsat影像進行比較，來評價MDCM方法的精度。 根據該區域特點，在本實施例中將窗口 w設置為3 x 3個MODIS像元即50 x 50個Landsat像元，地物類型數量w設置為4。圖21顯示 了分別由STARFM和MDCM方法生成的7月11日Landsat分辨率 的假彩色合成影像。如圖21所示為真實的7月11日Landsat影像左 和由STARFM中和MDCM方法右生成的Landsat分辨率影像。
由圖21可見，MDCM方法生成的影像和真實影像視覺上很接近， 細節特征保持較好，而STARFM的結果有些細節特征沒有突出，在 視覺上有影像被平滑的效果。這主要是因為微小地物在對應的低空 間分辨率上沒有對應的純像元存在，對于這些小地物^象元，STARFM 不能生成正確高分辨率反射率，而MDCM方法能克服這個問題。
為定量比較兩種方法結果的精度，現計算出預測反射率和真實反 射率之間的絕對誤差，即預測反射率和真實反射率之差的絕對值。 如圖22所示為真實影像與重建影像的各波段散點圖。上一排為 STARFM計算的結果，下一排為MDCM計算的結果。R是相關系數， AD為平均絕對誤差。
圖22顯示了兩種方法的結果各波段與真實值之間的散點圖。對 于綠光波段green band, MDCM方法的結果相比真實值有所偏高， 使得整幅影像的絕對誤差均值AD為0.0039,略高于STARFM的 0.0033;對于紅光波段red band, MDCM方法的結果比STARFM略 好一點，預測值與真實值的相關系數R為0.90,絕對誤差為0.0034, 而STARFM的R為0.89,絕對誤差為0.0037,從散點圖來看，MDCM 方法的結果偏高，而STARFM的結果偏低；對于近紅外波段NIR band, MDCM方法的結果明顯優于STARFM,預測值與真實值的相 關系數R為0.95，絕對誤差為0.0103,而STARFM的R為0.85， 絕對誤差為0.0155,從散點圖來看，MDCM方法的結果和真實值基 本一致，而STARFM有很多像元的反射率值比真實值低，而這些低 估的像元正是影像中的小面積地物板塊。
本發明還相應的提供了 一種遙感影像的處理系統。遙感影像為對 同一地表攝像而得到的兩組影像，其中，第一組影像具有高空間分 辨率，第二組影像具有低空間分辨率，根據所述兩組影像求得tp時刻的高空間分辨率的A影像，所述A影像由至少一個待計算的像元 組成。下面結合圖8~圖14對遙感影像的處理系統的實施方式進行介 紹。
如圖8所示，遙感影像的處理系統包括以下模塊影像選擇模塊 101 ，用于獲取同 一地表區域的具有高空間分辨率的第 一組影像和具 有低空間分辨率的第二組影像，從所述兩組影像中對應的選擇至少 兩對不同時刻的影像；目標像元設定模塊102,用于設定某時刻的所 求影像AO的一像元為目標像元，確定目標像元對應的地表區域為目 標區域；目標像元反射率計算模塊104,用于計算目標像元的反射率； 和，影像生成模塊108,用于根據目標像元反射率計算模塊104計算 得到的各目標像元的反射率生成AO影像。
所述目標像元反射率計算模塊104包括以下子模塊
轉換系數模塊103，根據所選的影像計算所述目標區域的反射率 變化量在兩組影像中之間的轉換系數。如圖9所示，在本實施例中， 所述轉換系數模塊103包括中心像元設定模塊1031，用于將所選 的第一組影像中的與所述目標區域對應的像元設定為中心像元；相 似像元篩選模塊1032，用于從所述中心像元各自所在影像中分別篩 選出至少 一 個與所述中心像元屬于同類地物的相似像元；轉換系數 計算才莫塊1033,用于確定各相似像元所對應的地表區域為相似區域， 根據所述相似區域在所選影像中的反射率計算所述轉換系數。調整 模塊1034，用于調整誤差導致的可能不正確的轉換系數和修正奇異 的轉換系數。調整模塊1034將在下文結合圖14進行詳細介紹。在 本實施例中，所述轉換系數計算模塊1033中所述轉換系數的計算是 將所述所獲得的反射率以第二變化量為自變量，第一變化量為因變 量根據各相似像元的權重進行加權回歸分析。相似像元的篩選的具 體實現可以參見上文中步驟1022的說明，在此不再贅述。
權重模塊109，用于根據所選的影像計算各相似像元的權重。下 文將結合圖11和圖12進行詳細介紹。
第二變化量計算模塊105，用于從第二組影像中選擇與所述影像AO時刻對應的影像B0和至少一個其他時刻的影像，根據所選擇的 影像計算所述目標區域在各所述其他時刻的影像中的反射率相對于 其在影像BO中的反射率的變化量。如圖10所示，所述第二變化量 計算模塊105包括相似像元第二變化量計算模塊1051,用于從第 二組影像中選擇與所述某時刻對應的影像BO和至少一個其他時刻 的影像，根據所選擇的影像計算各所述相似區域在第二組影像中的 反射率變化量；第二變化量加權平均模塊1052,根據得到的所述權 重將相似像元第二變化量計算模塊1051得到的各反射率變化量加權
平均計算所述目標區域在各所述其他時刻的影像中的反射率相對于 其在影像BO中的反射率的變化量。
變化量轉換模塊106，用于根據所述轉換系數將所述各第二變化 量計算模塊105得到的各變化量轉換為在第一組影像中的各反射率 變化量。
反射率計算模塊107,用于根據變化量轉換模塊106得到的各反 射率變化量和目標區域在第 一組相應時刻的各影像中的反射率計算 所述目標區域的反射率。在本實施例中經由第二變化量加權平均模 塊1052、變化量轉換模塊106和反射率計算模塊107完成了式14 的計算，得到了目標像元的反射率。影像生成模塊108,用于逐次計 算所有待計算的像元的反射率生成AO影像。
在另一個實施例中，目標像元反射率計算模塊104還包括時間 權重計算模塊111，用于根據第二變化量計算模塊105所得到的變化 量計算所得到的所述目標區域的各反射率的時間權重。反射率計算
模塊107還包括反射率加權計算模塊，用于根據所述時間權重， 將各反射率加權平均得到所求目標像元的反射率。時間權重的計算 和最終目標像元反射率的計算的具體實現可以參見上文中步驟 10810中的說明，在此不再贅述。
如圖ll所示，所述權重模塊109包括純度計算模塊1091,用 于根據所選影像計算各相似像元所有波段光譜向量和第二組影像中 各所述相似區域的像元所有波段光譜向量的 一致程度；純像元判斷
28模塊1093,用于判斷相似像元中是否有純像元，并當有純像元時將 該純像元權重設定為最大值，即，所有變化信息均取自該純像元， 當沒有純像元時不調整各相似像元的權重；空間距離計算模塊1094, 用于計算各相似像元距離各自中心像元的空間距離；權重計算模塊 1095，用于根據各相似像元的純度和距各自中心像元的空間距離計 算各相似像元的權重。權重模塊109中各模塊的具體實現可以參見上 文中步驟104的說明，在此不再贅述。
如圖12所示，在另一個實施例中，所述權重模塊109包括純度 計算模塊1091，用于根據所述至少兩對時刻對應的影像計算各相似像 元所有波段光譜向量和與其時刻對應的第二組影像中地表區域對應的 像元所有波段光譜向量的一致程度；相似像元權重計算模塊1092，用 于根據所述一致程度計算各相似像元的權重。具體純度計算和權重計 算可以參考權重模塊109的上一個實施例，在此不再贅述。
如圖13所示，在另一個實施例中，所述相似4象元篩選才莫塊1032 還包括交集篩選模塊10321,用于將各篩選出的相似像元所對應的地 表區域取交集，以取交集后的地表區域所對應的像元為相似像元。交 集篩選模塊10321的具體實現可以參見上文中步驟1024的說明，在此 不再贅述。
如圖14所示，所述調整模塊1034包括均值判斷模塊10341,用 于判斷各相似像元對應的第二變化量的均值是否小于反射率自身誤 差；誤差調整模塊10342，用于當均值判斷模塊(10341 )判斷為小于 時，轉換系數設定為使第一變化量等于第二變化量，當均值判斷模塊 (10341 )判斷為不小于時，轉換系數保持不變；奇異轉換系數判斷模 塊10343，用于判斷轉換系數是否奇異；奇異轉換系數修正模塊10344, 用于將奇異的轉換系數修正成。調整模塊1034的具體實現可以參見上 文中步驟1064和步驟1066的說明，在此不再贅述。
以上所述僅為本發明的較佳實施例而已，并不用以限制本發明， 凡在本發明的精神和原則之內，所作的任何修改、等同替換、改進等， 均應包含在本發明的保護范圍之內。
29
權利要求
1. 一種遙感影像的處理方法，其特征在于，設定某時刻的所求影像A0的一像元為目標像元，確定目標像元對應的地表區域為目標區域，獲取同一地表區域的具有高空間分辨率的第一組影像和具有低空間分辨率的第二組影像，從所述兩組影像中對應的選擇至少兩對不同時刻的影像，按以下步驟分別計算各目標像元的反射率以生成A0影像a，根據所選的影像計算所述目標區域的反射率變化量在兩組影像中之間的轉換系數；b，從第二組影像中選擇與所述影像A0時刻對應的影像B0和至少一個其他時刻的影像；根據所選擇的影像計算所述目標區域在各所述其他時刻的影像中的反射率相對于其在影像B0中的反射率的變化量；c，根據所述轉換系數將所述各反射率的變化量轉換為在第一組影像中的各反射率變化量；和d，從第一組影像中選擇與所述其它時刻對應的影像，根據轉換得到的各反射率變化量和目標區域在所選影像中的反射率計算所述目標區域的反射率。
2. 根據權利要求l所述方法，其特征在于，所述步驟a包括 al,將所選的第一組影像中的與所述目標區域對應的像元設定為中心像元；a2,從所述中心像元各自所在影像中分別篩選出至少 一個與所述 中心像元屬于同類地物的相似像元；和a3，確定各相似像元所對應的地表區域為相似區域，根據所述相 似區域在所選影像中的反射率計算所述轉換系數。
3. 根據權利要求2所述方法，其特征在于，所述步驟a3之前還 包括a21,根據所選的影像計算各相似像元的權重；所述步驟b包括b 1 ，根據B0影像和所述第二組的其他時刻的影像計算各所述相 似區域在第二組影像中的反射率變化量；b2,根據得到的所述權重將步驟bl得到的各反射率變化量加權 平均計算所述目標區域在各所述其他時刻的影像中的反射率相對于其在影像BO中的反射率的變化量。
4. 根據權利要求2所述方法，其特征在于，所述步驟a3中所述 轉換系數的計算是將各所述相似區域在所選影像中的反射率以在第 二組影像中的反射率的變化量為自變量，在第一組影像中的反射率 的變化量為因變量進行回歸分析。
5. 根據權利要求4所述方法，其特征在于，所述步驟a3之前還 包括a22,根據所選的影像計算各相似像元的權重；所述回歸分析是才艮據各相似像元的權重進行的加權回歸分析。
6. 根據權利要求3或5所述方法，其特征在于，所述各相似像 元的權重的計算包括以下步驟a221 ，根據所選影像計算各相似像元所有波段光譜向量和第二組 影像中各所述相似區域的像元所有波段光譜向量的 一 致程度； a222,根據所述一致程度計算各相似像元的權重。
7. 根據權利要求2所述方法，其特征在于，所述步驟a2還包括 將各篩選出的相似像元所對應的地表區域取交集，以取交集后的地表區域所對應的^f象元為相似4象元。
8. 遙感影像的處理系統，其特征在于，根據所述兩組影像求得 高空間分辨率的AO影像，所述AO影像不屬于第一組影像，并由至 少一個待計算的像元組成，所述處理系統包括影像選擇模塊(101),用于獲取同一地表區域的具有高空間分 辨率的第 一 組影像和具有低空間分辨率的第二組影像，從所述兩組影像中對應的選擇至少兩對不同時刻的影像；目標像元設定模塊(102),用于設定某時刻的所求影像A0的一像元為目標-像元，確定目標像元對應的地表區域為目標區域；目標像元反射率計算模塊(104),用于計算目標像元的反射率；和影像生成模塊(108 ),用于根據目標像元反射率計算模塊(104 ) 計算得到的各目標像元的反射率生成AO影像；所述目標像元反射率計算模塊(104)包括以下子模塊轉換系數模塊(103)，根據所選的影像計算所述目標區域的反 射率變化量在兩組影像中之間的轉換系數；第二變化量計算模塊(105)，用于從第二組影像中選擇與所述 影像A 0時刻對應的影像B 0和至少 一 個其他時刻的影像，根據所選 擇的影像計算所述目標區域在各所述其他時刻的影像中的反射率相 對于其在影像BO中的反射率的變化量；變化量轉換模塊(106),用于根據所述轉換系數將所述各第二 變化量計算模塊(105)得到的各變化量轉換為在第一組影像中的各 反射率變化量；和反射率計算模塊(107)，用于從第一組影像中選擇與所述其它 時刻對應的影像，根據變化量轉換模塊(106)得到的各反射率變化 量和目標區域在所選影像中的反射率計算所述目標區域的反射率。
9. 根據權利要求8所述系統，其特征在于，所述轉換系數模塊 (103 )包括中心像元設定模塊(1031),用于將所選的第一組影像中的與所 述目標區域對應的像元設定為中心像元；相似像元篩選模塊(1032),用于從所述中心像元各自所在影像 中分別篩選出至少 一個與所述中心像元屬于同類地物的相似像元；和轉換系數計算模塊(1033 ),用于確定各相似像元所對應的地表 區域為相似區域，根據所述相似區域在所選影像中的反射率計算所述 轉換系數。
10. 根據權利要求9所述系統，其特征在于，所述目標像元反射 率計算模塊(104)還包括權重模塊，用于根據所選的影像計算各相似像元的權重； 所述第二變化量計算模塊(105)包括相似像元第二變化量計算模塊(1051 ),用于從第二組影像中選 擇與所述某時刻對應的影像B0和至少一個其他時刻的影像，根據所 選擇的影像計算各所述相似區域在第二組影像中的反射率變化量；第二變化量加權平均模塊(1052)，根據得到的所述權重將相似 像元第二變化量計算模塊(1051 )得到的各反射率變化量加權平均計 算所述目標區域在各所述其他時刻的影像中的反射率相對于其在影 像BO中的反射率的變化量。
全文摘要
本發明公開了一種遙感影像的處理方法，根據空間分辨率不同的兩組影像求得某時刻的所求影像A0。首先根據空間分辨率不同的兩組影像計算目標區域的反射率變化量在這兩組影像之間的轉換系數，其次根據所述轉換系數將低空間分辨率影像中的變化量轉換成高空間分辨率影像中的變化量，最后由此根據已知的高空間分辨率影像得到所求影像A0的目標像元的反射率。還相應的提供了遙感影像的處理系統。本發明能生成將細節保持的更好的影像。
文檔編號G06K9/62GK101482929SQ20091007937
公開日2009年7月15日 申請日期2009年3月9日 優先權日2009年3月9日
發明者朱孝林, 晉 陳, 陳仲新, 陳學弘 申請人:中國農業科學院農業資源與農業區劃研究所;北京師范大學