專利名稱：航空高光譜遙感圖像光譜域噪聲自檢測與去除方法
技術領域：
本發(fā)明涉及一種航空高光譜遙感圖像光譜域噪聲自檢測與去除方法，該方法針對國產(chǎn)高光譜遙感器的高光譜圖像進行處理，可以檢測并去除高光譜圖像光譜域中噪聲，屬工程科學技術中的遙感技術應用領域。
背景技術：
在遙感圖像數(shù)字處理中，圖像噪聲的去除是一項較重要而基礎的工作。傳統(tǒng)的遙感圖像去噪與濾波總是從空間域對遙感圖像進行噪聲檢測與濾波以消除噪聲。高光譜圖像具有空間圖像和地物光譜兩個方面的信息，因而噪聲對高光譜圖像的影響最終也會表現(xiàn)為空間域噪聲和光譜域噪聲兩方面。傳統(tǒng)的多光譜遙感空間噪聲去除方法把高光譜每個波段圖像看作是孤立的單幅圖像，往往在空間域濾波后，光譜域中噪聲并沒有得到去除，有時反而會增強光譜域中噪聲。反射率轉換是高光譜圖像預處理的重要步驟之一。
目前高光譜圖像反射率轉換的途徑主要有三種基于地面同步實測光譜的經(jīng)驗線形法、基于大氣輻射傳輸模型的轉換方法和同時使用大氣輻射傳輸模型和地面實測光譜的混合方法，我國目前使用的一般是基于地面同步實測光譜的經(jīng)驗線性法。然而由于高光譜圖像成像時自然光照明條件影響、地面地形影響、混合像元問題等原因引入各種噪聲，表現(xiàn)為轉換后的反射率曲線呈鋸齒狀。這種噪聲存在于光譜域中，是需要檢測并予以去除的。
Boardman，J W.注意到ATREM得到的光譜表現(xiàn)出一種波段間的系統(tǒng)誤差。這種誤差表現(xiàn)為鋸齒狀噪聲，使得光譜看起來很不自然。他認為這些連貫誤差是增益誤差、儀器校正誤差和大氣與太陽模型影響的累積。在此基礎上他提出了所謂經(jīng)驗平面場優(yōu)化反射率變換EFFORT(the Empirical Flat Field Optimal ReflectanceTransformation)方法。經(jīng)驗平面場優(yōu)化反射率變換EFFORT方法通過高光譜圖像本身的統(tǒng)計找到一種輕微的改正(這種改正的增益系數(shù)接近1，偏移系數(shù)接近0)，這種改正可以使誤差減小進而提高反射率的精度，使得反射率曲線有了一個視覺上的改進。EFFORT也可以使用與圖像中對應地物的實測光譜參與改正，以保存一些光譜中的明顯變化(比如植被紅邊)。
John R.Jensen也提到過這個問題，他認為反射率轉換后光譜中會存在相當大的噪聲，并采用一種線性校正方法稱為單光譜增強SSE(Single Spectrum Enhancement)的方法(在美國高光譜“即時大氣校正”ACORN軟件中包含了此種方法)。該方法使用了同一地面目標(應為研究區(qū)域的主要地物類型)的轉換后反射率曲線與實測光譜反射率曲線，并將轉換后光譜擬合到實測地物光譜上。
常用的曲線平滑法包括均值濾波平滑法、加權平均法與最小二乘平滑法，這些方法也可以用來去除光譜域噪聲。但如果使用目前常用的一維空間曲線平滑方法進行光譜噪聲去除，往往會在去除噪聲的同時，也將有用的細微光譜吸收等光譜特征同時去除了。而經(jīng)驗平面場優(yōu)化反射率變換EFFORT和單光譜增強SSE方法需要一定的地面實測資料，且處理的高光譜圖像區(qū)域中地表覆蓋類型不能過于復雜，其應用受到一定限制。
中國于80年代中后期開始著手發(fā)展自己的高光譜圖像分析系統(tǒng)，90年代初，中科院上海技術物理所完成了模塊式航空成像光譜儀MAIS新型模塊化航空成像光譜儀實用化系統(tǒng)；在“九五”期間，中科院上海技術物理所又相繼研制成功128波段實用型模塊化成像光譜儀系統(tǒng)(OMIS)、244波段推帚式高光譜成像儀(PHI)。這些成果標志著我國在高光譜系統(tǒng)領域取得了長足進展。然而由于高光譜系統(tǒng)本身的限制，高光譜圖像在光譜域仍然會存在著較大噪聲。為促進高光譜技術的進一步應用，有必要針對高光譜數(shù)據(jù)的特點，研制出一種用于光譜域噪聲自檢測與去除的方法，用以提高數(shù)據(jù)的質量。

發(fā)明內(nèi)容
為了克服上述不足之處，本發(fā)明的主要目的旨在提供一種適用于光譜噪聲評定的變量、判定適用于光譜噪聲評定的條件及找到適用于高光譜圖像光譜域濾波器的方法，實現(xiàn)檢測區(qū)分高光譜圖像光譜域噪聲或細微光譜特征，對光譜域較大噪聲去除，對光譜域的細微特征或較小噪聲不予處理或進行更小程度處理的航空高光譜遙感圖像光譜域噪聲自檢測與去除方法。
本發(fā)明要解決的技術問題是要解決如何找到一種適用于光譜噪聲評定的變量、判定適用于光譜噪聲評定的條件問題；要解決如何找到適用于高光譜圖像光譜域濾波器的方法系統(tǒng)軟件及有關硬件等技術問題。
本發(fā)明解決其技術問題所采用的技術方案是航空高光譜遙感圖像光譜域噪聲自檢測與去除方法是通過一種適用于光譜噪聲評定的變量、判定適用于光譜噪聲評定的條件及找到適用于高光譜圖像光譜域濾波器的方法，實現(xiàn)檢測區(qū)分高光譜圖像光譜域噪聲或細微光譜特征，對光譜域較大噪聲去除，對光譜域的細微特征或較小噪聲不予處理或進行更小程度處理，其具體工作步驟是
步驟1高光譜遙感圖像高光譜遙感圖像為標準格式的機載航空高光譜遙感反射率圖像；步驟2圖像光譜反射率的提取a).高光譜遙感圖像模塊的輸出信號傳送到圖像光譜反射率曲線模塊；b).打開高光譜圖像，讀取高光譜圖像的基本信息其中設置波段之間寬度為bw納米，圖像存儲類型為整型數(shù)int、浮點數(shù)float或其他，像元值存儲字節(jié)長度nc，文件頭字節(jié)數(shù)no，圖像波段數(shù)nb，圖像的寬度上像元數(shù)wp和高度上像元數(shù)hp；c).計算反射率值，建立二維數(shù)組由圖像基本信息可以得到圖像第i個像元的第j個波段的反射率值，從文件頭偏移值ro由下式計算ro＝no+(wp*hp)*(j-1)*nc+i-1 (1)建立二維數(shù)組rb[][]，rb維數(shù)分別為wp*hp和nb，從p＝1到p＝wp*hp，b＝1到b＝nb依次讀取每個像元每個波段的反射率值，并存儲在rb[p][b]中；d).二維數(shù)組rb[][]為提取出的圖像光譜反射率曲線；步驟3圖像光譜反射率曲線圖像光譜反射率曲線為步驟2中提取出的數(shù)組rb[][]，圖像光譜反射率曲線模塊的輸出信號一路傳送到一次濾波后反射率曲線模塊的輸入端，另一路傳送到二次濾波后反射率曲線模塊的輸入端；步驟4反射率二階導數(shù)的計算a).圖像光譜反射率曲線模塊的輸出信號傳送到反射率曲線二階導數(shù)模塊；
b).用高光譜圖像光譜的二階導數(shù)作為光譜噪聲評定的變量其中光譜導數(shù)技術包括對反射光譜進行數(shù)學模擬和計算不同階數(shù)的導數(shù)值；確定光譜彎曲點及最大最小反射率的波長位置；光譜導數(shù)處理強調(diào)曲線的變化和壓縮均值影響；c).計算光譜的二階導數(shù)光譜的二階導數(shù)計算公式如下d2sdλ2|j≈s(λi)-2s(λj)+s(λk)(Δλ)2]]>(2)式中Δλ＝λk-λj＝λj-λi，λk＞λj＞λi；λi為高光譜圖像第i個波段的波長，Δλ為高光譜圖像的波段之間寬度bw，s(λi)為像元反射率光譜曲線第i個波段的反射率值；d).具體計算二階導數(shù)值對圖像光譜反射率曲線數(shù)組rb[][]，對第p個像元其光譜二階導數(shù)計算時，第b個波段的二階導數(shù)值為(rb[p][b-1]-2*rb[p][b]+rb[p][b+1])/bw2(3)從p＝1到p＝wp*hp，b＝2到b＝nb-1依次進行光譜二階導數(shù)計算，并將計算結果保存在數(shù)組sd[][]中，維數(shù)分別為wp*hp和nb；在數(shù)組sd[p][b]中，當b＝1和b＝nb時，數(shù)組元素值為0；e).數(shù)組sd[][]為反射率曲線二階導數(shù)；步驟5反射率曲線二階導數(shù)反射率曲線二階導數(shù)為步驟4中提取出的數(shù)組sd[][]；步驟6高光譜圖像光譜域噪聲自檢測及反射率曲線噪聲判定系數(shù)的計算a).反射率曲線二階導數(shù)模塊的輸出信號傳送到反射率曲線噪聲判定系數(shù)模塊；b).用高光譜圖像光譜二階導數(shù)的標準差作為光譜噪聲的判定條件其中光譜的低階導數(shù)處理對噪聲影響敏感性較低，高階微分對噪聲影響敏感度高；一階導數(shù)反映了反射率曲線的斜率，二階導數(shù)反映了噪聲的實際分布情況；c).用光譜二階導數(shù)對光譜曲線進行噪聲影響程度檢測1).進行光譜二階導數(shù)計算將二階導數(shù)值s″(λi)與其平均值μ之差和給定的某一閾值作比較；2).判定該波段噪聲的大小；3).確定該波段濾波平滑窗的大??；4).判斷不等式對某一波段i而言，如果不等式|s″(λi)-μ|＞σ (4)成立的話，那么該波段反射率值認為是有較大噪聲存在，反射率曲線噪聲判定系數(shù)為1，表明該波段噪聲較大，式中σ是光譜二階導數(shù)的標準差；如果不等式|s″(λi)-μ|≤σ (5)成立的話，那么該波段反射率值認為存在噪聲較小，反射率曲線噪聲判定系數(shù)為0；d).具體計算二階導數(shù)的平均值1).對反射率曲線二階導數(shù)數(shù)組sd[][]，第p個像元的反射率曲線二階導數(shù)的平均值msd為msd=(Σb=1nbsd[p][b])/nb]]>(6)
2).第p個像元的反射率曲線二階導數(shù)的標準差ssd為ssd=[Σb=1nb(sd[p][b]-msd)]/nb]]>(7)3).建立二維數(shù)組de[][]其維數(shù)分別為wp*hp和nb；4).判斷不等式對第p個像元的第b個波段而言，如果不等式|sd[p][b]-msd|＞ssd (8)成立的話，de[p][b]＝1；如果不等式|sd[p][b]-msd|≤ssd (9)成立的話，de[p][b]＝0；從p＝1到p＝wp*hp，b＝1到b＝nb依次進行反射率曲線噪聲判定系數(shù)計算，并將計算結果保存在數(shù)組de[][]中；e).數(shù)組de[][]為反射率曲線噪聲判定系數(shù)；步驟7反射率曲線噪聲判定系數(shù)反射率曲線噪聲判定系數(shù)為步驟6中計算得到的數(shù)組de[][]；步驟8噪聲的去除及圖像光譜反射率曲線噪聲的第一次濾波a).反射率曲線噪聲判定系數(shù)模塊的輸出信號傳送到一次濾波后反射率曲線模塊；b).用賽維特斯基-高勒Savitzky-Golay平滑濾波器；c).用簡化的最小二乘擬合卷積方法對曲線進行平滑處理；d).計算平滑后曲線各階導數(shù)1).簡化后的最小二乘卷積方程式如下Yj=Σi=-mmCiyj+iN]]>(10)
式中，y是原始光譜值，Y是平滑后光譜值；Ci是平滑窗口中第i個光譜值的系數(shù)，N是卷積中點值個數(shù)，j是沿原始數(shù)據(jù)縱坐標數(shù)據(jù)列的計算點下標；2).平滑濾波法計算的卷積點為25點，并計算到光譜的第六階導數(shù)；3).各階導數(shù)平滑系數(shù)的計算公式其零階導數(shù)二次或三次多項式擬合的任意點數(shù)平滑窗系數(shù)計算公式如下ps(0)=3(3m2+3m-1-5s2)(2m+3)(2m+1)(2m-1),s=-m,...,0,...m]]>(12)式中，m為平滑窗的寬度的一半；e).檢測在應用光譜學中光譜處理方法之前，對光譜中誤差嚴重程度進行檢測；1).對噪聲嚴重的地方進行高點數(shù)的平滑；2).對噪聲較小部分應用較小平滑窗進行平滑；f).具體計算第一次濾波后圖像反射率曲線1).建立一維數(shù)組sg1[11]；2).用來存儲11點賽維特斯基-高勒濾波器系數(shù)；3).計算sg1[]的值按如下公式sg1[6+m]=3(3*52+3*5-1-5*m2)(2*5+3)(2*5+1)(2*5-1),m=-5,...0,...5]]>(13)計算出sg1[]的值；4).建立二維數(shù)組rf1[][]其維數(shù)分別為wp*hp和nb；5)判斷如果第p個像元的第b個波段的反射率曲線噪聲判定系數(shù)de[p][b]＝1，那么對第p個像元的第b個波段反射率值進行光譜域濾波，按如下公式rf1[p][b]=Σi=-55sg1[6+i]*rb[p][b+i]11]]>(14)式中5＜b＜nb-5；如果第p個像元的第b個波段的反射率曲線噪聲判定系數(shù)de[p][b]＝0，按如下公式rf1[p][b]＝rb[p][b] (15)對圖像中每個像元逐波段進行計算；g).二維數(shù)組rf1[][]為第一次濾波后，去除了較大噪聲的圖像反射率曲線；步驟9一次濾波后反射率曲線一次濾波后反射率曲線為步驟8中計算出的二維數(shù)組rf[][]；步驟10噪聲的去除及圖像光譜反射率曲線噪聲的第二次濾波(10)a).一次濾波后反射率曲線模塊的輸出信號傳送到二次濾波后反射率曲線模塊；b).具體計算第二次濾波后圖像反射率曲線1).建立一維數(shù)組sg2[5]；2).用來存儲5點賽維特斯基-高勒濾波器系數(shù)；3).計算sg2[]的值按如下公式sg1[3+m]=3(3*22+3*2-1-2*m2)(2*2+3)(2*2+1)(2*2-1),m=-2,...0,...2]]>(16)計算出sg2[]的值；
4).建立二維數(shù)組rf2[][]其維數(shù)分別為wp*hp和nb；5)判斷對第p個像元的第b個波段反射率值進行光譜域濾波，按如下公式rf2[p][b]=Σi=-55sg1[3+i]*rf1[p][b+i]5]]>(17)式中2＜b＜nb-1；當b＝1，2，nb-1，nb時，按如下公式rf2[p][b]＝rb[p][b] (18)對圖像中每個像元逐波段進行計算；c).二維數(shù)組rf2[][]為第二次濾波后的圖像反射率曲線；步驟11二次濾波后反射率曲線二次濾波后反射率曲線為步驟10中得到的二維數(shù)組rf2[][]；步驟12由濾波后反射率曲線重構高光譜圖像a).二次濾波后反射率曲線模塊的輸出信號傳送到光譜域濾波后高光譜遙感圖像模塊；b).建立新文件按航空高光譜圖像的格式標準，高光譜圖像文件大小為(no+wp*hp*nb)*nc個字節(jié)，建立新文件；c).寫入新文件將步驟1中讀取出的高光譜圖像文件頭的no*nc個字節(jié)，寫入新文件中；d).寫入對p＝1，…，wp*hp，依次寫入rf2[p][1]，然后再從p＝1，…，wp*hp，依次寫入rf2[p][2]，直到從p＝1，…，wp*hp，依次寫入rf2[p][nb]；
e).重構完成高光譜圖像文件重構完成；步驟13光譜域濾波后高光譜遙感圖像光譜域濾波后高光譜遙感圖像為步驟12中重構的高光譜圖像文件。
本發(fā)明的有益效果是本方法能有效的去除光譜中存在噪聲，并能保留光譜原有的大部分光譜特征，是地物實測光譜和高光譜圖像預處理的一種有效手段；按本方法對高光譜圖像進行光譜域噪聲的自檢測與濾波，可以在去除高光譜圖像光譜域噪聲的同時去除因光譜噪聲引起的空間域噪聲。


下面結合

和實施例對本發(fā)明進一步說明。
附圖1為本發(fā)明總體流程示意圖；附圖2為本發(fā)明原始圖像中植被反射率光譜示意圖；附圖3為本發(fā)明用本方法對噪聲自檢測與去除后圖像中植被反射率光譜示意圖；附圖4為本發(fā)明含有噪聲的圖像反射率光譜示意圖；附圖5為本發(fā)明五點賽維特斯基-高勒Savitzky-Golay濾波后光譜曲線示意圖；附圖6為本發(fā)明十一點賽維特斯基-高勒Savitzky-Golay濾波后光譜曲線示意圖；附圖7為用本發(fā)明方法處理后光譜曲線示意圖；附圖8為本發(fā)明實測反射率光譜示意圖；附圖9為本發(fā)明經(jīng)驗線性法得到的圖像反射率光譜示意圖；附圖10為本發(fā)明最小噪聲變換MNF去噪后光譜示意圖；附圖11為本發(fā)明經(jīng)驗平面場優(yōu)化反射率變換EFFORT方法處理后光譜示意圖；
附圖12為本發(fā)明光譜均值濾波后光譜示意圖；附圖13為本發(fā)明用本方法處理后光譜示意圖；附圖中標號說明1-高光譜遙感圖像；2-圖像光譜反射率的提?。?-圖像光譜反射率曲線；4-反射率二階導數(shù)的計算；5-反射率曲線二階導數(shù)；6-噪聲判定系數(shù)的計算；7-反射率曲線噪聲判定系數(shù)；8-反射率曲線噪聲的第一次濾波；9-一次濾波后反射率曲線；10-噪聲的第二次濾波；11-二次濾波后反射率曲線；12-重構高光譜圖像；13-光譜域濾波后高光譜遙感圖像；
具體實施例方式請參閱附圖1所示，本發(fā)明解決的最重要問題是找到方法，能檢測區(qū)分高光譜圖像光譜域噪聲或細微光譜特征，對光譜域較大噪聲去除，對光譜域的細微特征或較小噪聲不予處理或進行更小程度處理。本技術分三個內(nèi)容1)找到一種適用于光譜噪聲評定的變量；2)判定適用于光譜噪聲評定的條件；3)找到適用于高光譜圖像光譜域的濾波器。
本發(fā)明通過一種適用于光譜噪聲評定的變量、判定適用于光譜噪聲評定的條件及找到適用于高光譜圖像光譜域濾波器的方法，實現(xiàn)檢測區(qū)分高光譜圖像光譜域噪聲或細微光譜特征，對光譜域較大噪聲去除，對光譜域的細微特征或較小噪聲不予處理或進行更小程度處理，其具體工作步驟是步驟1高光譜遙感圖像(1)高光譜遙感圖像(1)為標準格式的機載航空高光譜遙感反射率圖像(可參考中科院上海技術物理研究所的標準)；步驟2圖像光譜反射率的提取(2)高光譜遙感圖像(1)模塊的輸出信號傳送到圖像光譜反射率曲線(3)模塊；打開高光譜圖像，讀取高光譜圖像的基本信息波段之間寬度bw(單位為納米)，圖像存儲類型(整型數(shù)(int)、浮點數(shù)(float)或其他)，像元值存儲字節(jié)長度nc(整型數(shù)為2個字節(jié)，浮點數(shù)為4個字節(jié))，文件頭字節(jié)數(shù)no，圖像波段數(shù)nb，圖像的寬度上像元數(shù)wp和高度上像元數(shù)hp。由圖像基本信息可以得到圖像第i個像元的第j個波段的反射率值在從文件頭偏移值ro由下式計算ro＝no+(wp*hp)*(j-1)*nc+i-1 (1)建立二維數(shù)組rb[][]，rb維數(shù)分別為wp*hp和nb。從p＝1到p＝wp*hp，b＝1到b＝nb依次讀取每個像元每個波段的反射率值，并存儲在rb[p][b]中，二維數(shù)組rb[][]即為提取出的圖像光譜反射率曲線。
步驟3圖像光譜反射率曲線(3)圖像光譜反射率曲線(3)即為步驟2中提取出的數(shù)組rb[][]，圖像光譜反射率曲線(3)模塊的輸出信號一路傳送到一次濾波后反射率曲線(9)模塊的輸入端，另一路傳送到二次濾波后反射率曲線(11)模塊的輸入端。
步驟4反射率二階導數(shù)的計算(4)圖像光譜反射率曲線(3)模塊的輸出信號傳送到反射率曲線二階導數(shù)(5)模塊；
本技術采用高光譜圖像光譜的二階導數(shù)作為光譜噪聲評定的變量。光譜導數(shù)技術包括對反射光譜進行數(shù)學模擬和計算不同階數(shù)的導數(shù)值以迅速地確定光譜彎曲點及最大最小反射率的波長位置。光譜導數(shù)處理強調(diào)曲線的變化和壓縮均值影響。光譜導數(shù)技術對光譜信噪比非常敏感。光譜的二階導數(shù)計算公式如下d2sdλ2|j≈s(λi)-2s(λj)+s(λk)(Δλ)2]]>(2)其中Δλ＝λk-λj＝λj-λi，λk＞λj＞λi。λi為高光譜圖像第i個波段的波長，Δλ為高光譜圖像的波段之間寬度bw(見步驟2)，s(λi)為像元反射率光譜曲線第i個波段的反射率值。
具體計算步驟為對圖像光譜反射率曲線數(shù)組rb[][]，對第p個像元其光譜二階導數(shù)計算時，第b個波段的二階導數(shù)值為(rb[p][b-1]-2*rb[p][b]+rb[p][b+1])/bw2(3)從p＝1到p＝wp*hp，b＝2到b＝nb-1依次進行光譜二階導數(shù)計算，并將計算結果保存在數(shù)組sd[][]中，維數(shù)分別為wp*hp和nb。注意，由于公式(3)的限制，在數(shù)組sd[p][b]中，當b＝1和b＝nb時，數(shù)組元素值為0。數(shù)組sd[][]即為反射率曲線二階導數(shù)。
步驟5反射率曲線二階導數(shù)(5)反射率曲線二階導數(shù)即為步驟4中提取出的數(shù)組sd[][]。
步驟6高光譜圖像光譜域噪聲自檢測-反射率曲線噪聲判定系數(shù)的計算(6)反射率曲線二階導數(shù)(5)模塊的輸出信號傳送到反射率曲線噪聲判定系數(shù)(7)模塊；本技術采用高光譜圖像光譜二階導數(shù)的標準差作為光譜噪聲的判定條件。研究表明，光譜的低階導數(shù)處理對噪聲影響敏感性較低，而高階微分對噪聲影響敏感度高。但一階導數(shù)也同時反映了反射率曲線的斜率，二階導數(shù)則很好的反映了噪聲的實際分布情況。因此可以考慮用光譜二階導數(shù)對光譜曲線進行噪聲影響程度檢測。首先進行光譜二階導數(shù)計算，將二階導數(shù)值s″(λi)的與其平均值μ之差和給定的某一閾值作比較，以判定該波段噪聲的大小，從而確定該波段濾波平滑窗的大小。計算結果表明，對某一波段i而言，如果不等式|s″(λi)-μ|＞σ (4)成立的話，那么該波段反射率值認為是有較大噪聲存在，反射率曲線噪聲判定系數(shù)為1，表明該波段噪聲較大，這里的σ是光譜二階導數(shù)的標準差。如果不等式|s″(λi)-μ|≤σ (5)成立的話，那么該波段反射率值認為存在噪聲較小，反射率曲線噪聲判定系數(shù)為0。
具體計算步驟如下對反射率曲線二階導數(shù)數(shù)組sd[][]，第p個像元的反射率曲線二階導數(shù)的平均值msd為msd=(Σb=1nbsd[p][b])/nb]]>(6)第p個像元的反射率曲線二階導數(shù)的標準差ssd為ssd=[Σb=1nb(sd[p][b]-msd)]/nb]]>(7)建立二維數(shù)組de[][]，其維數(shù)分別為wp*hp和nb。對第p個像元的第b個波段而言，如果不等式|sd[p][b]-msd|＞ssd (8)成立的話，de[p][b]＝1。如果不等式|sd[p][b]-msd|≤ssd (9)成立的話，de[p][b]＝0。從p＝1到p＝wp*hp，b＝1到b＝nb依次進行反射率曲線噪聲判定系數(shù)計算，并將計算結果保存在數(shù)組de[][]中，數(shù)組de[][]即為反射率曲線噪聲判定系數(shù)。
步驟7反射率曲線噪聲判定系數(shù)(7)反射率曲線噪聲判定系數(shù)(7)即為步驟6中計算得到的數(shù)組de[][]。
步驟8噪聲的去除-圖像光譜反射率曲線噪聲的第一次濾波(8)反射率曲線噪聲判定系數(shù)(7)模塊的輸出信號傳送到一次濾波后反射率曲線(9)模塊；本技術采用賽維特斯基-高勒(Savitzky-Golay)平滑濾波器。使用簡化的最小二乘擬合卷積方法對曲線進行平滑處理并可計算平滑后曲線各階導數(shù)。簡化后的最小二乘卷積一般方程式如下Yj=Σi=-mmCiyj+iN]]>(10)式中，y是原始光譜值，Y是平滑后光譜值。Ci是平滑窗口中第i個光譜值的系數(shù)，N是卷積中點值個數(shù)，j是沿原始數(shù)據(jù)縱坐標數(shù)據(jù)列的計算點下標。賽維特斯基-高勒(Savitzky-Golay)平滑濾波法可以計算的卷積點最高為25點，并可計算到光譜的第六階導數(shù)。Madden改正了原賽維特斯基-高勒(Savitzky-Golay)系數(shù)的一些錯誤，給出了各階導數(shù)平滑系數(shù)的計算公式。其零階導數(shù)二次或三次多項式擬合的任意點數(shù)平滑窗系數(shù)計算公式如下ps(0)=3(3m2+3m-1-5s2)(2m+3)(2m+1)(2m-1),s=-m,...,0,...m]]>(12)式中，m為平滑窗的寬度的一半。使用Madden提供的方程可以計算零到六階導數(shù)的最小二乘擬合卷積系數(shù)值，從而計算光譜的零到六階平滑后的導數(shù)值。由于卷積公式本身限制，光譜兩端各m個點值沒有進行計算。因此，得到平滑后的光譜要比原來的光譜要短。
賽維特斯基-高勒(Savitzky-Golay)平滑濾波因為可保留光譜的一些細微特征(如光譜吸收峰)，因此在光譜學中應用非常廣泛，但由于高光譜遙感數(shù)據(jù)與光譜學中光譜的差別，光譜學中應用的方法不一定能直接用于高光譜遙感數(shù)據(jù)的分析中。光譜學中的數(shù)據(jù)是在嚴格的實驗室控制條件下進行采集的，可近似認為數(shù)據(jù)的誤差呈一定分布。而高光譜數(shù)據(jù)在成像時由于受自然光照明條件、地面地形起伏、大氣衰減和鏡頭的工作狀態(tài)等影響，數(shù)據(jù)的誤差呈現(xiàn)不確定性。在應用光譜學中光譜處理方法之前需要對光譜中誤差嚴重程度進行檢測，對噪聲嚴重的地方進行更高點數(shù)的平滑，而對噪聲較小部分應用較小平滑窗進行平滑。
具體計算步驟如下建立一維數(shù)組sg1[11]用來存儲11點賽維特斯基-高勒濾波器系數(shù)。按如下公式sg1[6+m]=3(3*52+3*5-1-5*m2)(2*5+3)(2*5+1)(2*5-1),m=-5,...0,...5]]>(13)計算出sg1[]的值。建立二維數(shù)組rf1[][]，其維數(shù)分別為wp*hp和nb。如果第p個像元的第b個波段的反射率曲線噪聲判定系數(shù)de[p][b]＝1，那么對第p個像元的第b個波段反射率值進行光譜域濾波，按如下公式rf1[p][b]=Σi=-55sg1[6+i]*rb[p][b+i]11]]>(14)注意，式中5＜b＜nb-5。如果第p個像元的第b個波段的反射率曲線噪聲判定系數(shù)de[p][b]＝0，按如下公式rf1[p][b]＝rb[p][b] (15)對圖像中每個像元逐波段進行計算。得到的二維數(shù)組rf1[][]即為第一次濾波后去除了較大噪聲的圖像反射率曲線。
步驟9一次濾波后反射率曲線(9)
一次濾波后反射率曲線即為步驟8中計算出的二維數(shù)組rf[][]。
步驟10噪聲的去除-圖像光譜反射率曲線噪聲的第二次濾波(10)一次濾波后反射率曲線(9)模塊的輸出信號傳送到二次濾波后反射率曲線(11)模塊；具體計算步驟如下建立一維數(shù)組sg2[5]用來存儲5點賽維特斯基-高勒濾波器系數(shù)。按如下公式sg1[3+m]=3(3*22+3*2-1-2*m2)(2*2+3)(2*2+1)(2*2-1),m=-2,...0,...2]]>(16)計算出sg2[]的值。建立二維數(shù)組rf2[][]，其維數(shù)分別為wp*hp和nb。對第p個像元的第b個波段反射率值進行光譜域濾波，按如下公式rf2[p][b]=Σi=-55sg1[3+i]*rf1[p][b+i]5]]>(17)注意，式中2＜b＜nb-1。當b＝1，2，nb-1，nb時，按如下公式rf2[p][b]＝rb[p][b] (18)對圖像中每個像元逐波段進行計算。得到的二維數(shù)組rf2[][]即為第二次濾波后的圖像反射率曲線。
步驟11二次濾波后反射率曲線(11)二次濾波后反射率曲線(11)即為步驟10中得到的二維數(shù)組rf2[][]。
步驟12由濾波后反射率曲線重構高光譜圖像(12)二次濾波后反射率曲線(11)模塊的輸出信號傳送到光譜域濾波后高光譜遙感圖像模塊；按航空高光譜圖像的格式標準，高光譜圖像文件大小為(no+wp*hp*nb)*nc個字節(jié)，建立新文件。將步驟1中讀取出的高光譜圖像文件頭的no*nc個字節(jié)，寫入新文件中。然后對p＝1，…，wp*hp，依次寫入rf2[p][1]，然后再從p＝1，…，wp*hp，依次寫入rf2[p][2]，直到從p＝1，…，wp*hp，依次寫入rf2[p][nb]。高光譜圖像文件重構完成。
步驟13光譜域濾波后高光譜遙感圖像(13)光譜域濾波后高光譜遙感圖像(13)即為步驟12中重構的高光譜圖像文件。
本發(fā)明的具體實施例按照本發(fā)明的技術流程圖各步驟，對上海2003年世博會區(qū)域PHI高光譜圖像進行光譜域噪聲的自檢測與去除，處理效果比較好。
請參閱附圖2所示，為原始圖像中植被反射率光譜；請參閱附圖3所示，為用本方法對噪聲自檢測與去除后圖像中植被反射率光譜。從圖3中可以看出，本專利方法DSGF對高光譜圖像最大的作用首先體現(xiàn)在對光譜平滑上。未經(jīng)光譜域噪聲去除的圖像光譜曲線由于受各種噪聲影響而呈現(xiàn)鋸齒狀(圖2)，在400nm-500nm與750nm-980nm之間波段尤其明顯。經(jīng)本方法處理后，圖像光譜在這兩個光譜區(qū)間明顯平滑許多(圖3)，同時從圖3還可以看到，光譜曲線仍保持了一些如吸收峰等細微的光譜特征。
請參閱附圖4所示，為含有噪聲的圖像反射率光譜；請參閱附圖5所示，為五點賽維特斯基-高勒(Savitzky-Golay)濾波后光譜曲線圖；請參閱附圖6所示，為十一點賽維特斯基-高勒(Savitzky-Golay)濾波后光譜曲線圖；
請參閱附圖7所示，為本發(fā)明方法處理后光譜曲線圖。賽維特斯基-高勒(Savitzky-Golay)濾波的效果由參與濾波的點數(shù)決定，點數(shù)越高，平滑程度也越高，但初始光譜的特征就越難以保留。從圖4中可以看出，初始反射率光譜曲線的兩個矩形框區(qū)域呈現(xiàn)不同程度的鋸齒狀，其中尤其以近950nm附近波段為據(jù)。這兩個區(qū)域的光譜值表明該波段區(qū)域含有噪聲而且噪聲的程度不一。從圖5中看到，5點濾波使得噪聲情況有所緩解，同時也保持了大部分的光譜細微特征，但在噪聲較嚴重區(qū)域仍然呈鋸齒狀；從圖6中看到，11點濾波后的光譜非常平滑，然而一些正常的光譜特征也被掩蓋；從圖7中看到，本發(fā)明方法處理后的光譜曲線在噪聲嚴重的區(qū)域能較徹底的清除噪聲，同時也能保持光譜原有的細微特征，這證明本發(fā)明方法是一種較好的高光譜圖像光譜域噪聲自檢測與去除方法。
為了與高光譜其他一些噪聲去除方法進行比較，我們還使用了多種方法對高光譜圖像進行處理。
請參閱附圖8所示，為實測反射率光譜；請參閱附圖9所示，為經(jīng)驗線性法得到的圖像反射率光譜；請參閱附圖10所示，為最小噪聲變換MNF去噪后光譜；請參閱附圖11所示，為經(jīng)驗平面場優(yōu)化反射率變換EFFORT方法處理后光譜；請參閱附圖12所示，為光譜均值濾波后光譜；請參閱附圖13所示，為本發(fā)明方法處理后光譜。
從圖9可以看出，與圖8的圖像光譜相比，利用經(jīng)驗線性法得到的圖像反射率包含有較大噪聲，如果在高光譜遙感應用中直接利用這種包含較大噪聲的反射率值的話，無疑會使得應用結果也呈現(xiàn)噪聲值。因此，在進一步高光譜應用處理步驟之前，對圖像進行光譜域的噪聲去除是有必要的。圖10為最小噪聲變換MNF后，對最小噪聲變換MNF變換結果中波段進行噪聲去除，然后進行MNF逆變換后得到的圖像光譜，最小噪聲變換MNF方法應用正交變化進行噪聲分離，總體上說來，對噪聲的去除效果還是比較不錯的，但位于950nm波段附近的反射率值仍然有鋸齒狀噪聲存在。
圖11中可以看到經(jīng)驗平面場優(yōu)化反射率變換EFFORT方法對光譜處理的效果不是很理想。圖12和圖13表明，光譜均值濾波和本專利方法DSGF方法濾波對光譜的處理效果相對其他方法而言要好。而光譜均值濾波后的光譜過于平滑，雖然在視覺上與實測光譜非常的相符，但由于高光譜的多種噪聲影響，加上高光譜的波段寬度比地物光譜儀的波段寬度要大，因而圖像光譜與實測光譜應該是有著較大的區(qū)別的，這些區(qū)別應該表現(xiàn)在光譜的一些細微特征上，光譜均值濾波雖然保留了光譜的整體趨勢，但是光譜的細微特征不能辨別。
為了進一步比較各種處理方法處理后圖像光譜與原圖像光譜和實測光譜的相似程度，對上述六種光譜兩兩之間底光譜歐幾里德距離和光譜角距離進行計算，得出了如表1和表2的結果。表1四種方法處理后的光譜與實測光譜和原圖像光譜的歐幾里德距離

(接下頁)
表2四種方法處理后的光譜與實測光譜和原圖像光譜的光譜角距離

由歐幾里德距離與光譜角距離的定義可知，距離值越小，兩光譜越相似。從表1、表2中可以看出，圖十三中所示用本發(fā)明方法處理后光譜與均值濾波光譜相比，圖十三所示光譜與實測光譜的距離和與圖像光譜的距離都要小于均值濾波光譜與實測光譜和圖像光譜的距離。這說明，均值濾波后光譜過于平滑，而本發(fā)明方法處理后光譜保留了光譜的大部分特征，在進行光譜匹配時最具有優(yōu)勢。
權利要求
1.一種航空高光譜遙感圖像光譜域噪聲自檢測與去除方法，其特征在于該方法通過一種適用于光譜噪聲評定的變量、判定適用于光譜噪聲評定的條件及找到適用于高光譜圖像光譜域濾波器的方法，實現(xiàn)檢測區(qū)分高光譜圖像光譜域噪聲或細微光譜特征，對光譜域較大噪聲去除，對光譜域的細微特征或較小噪聲不予處理或進行更小程度處理，其具體工作步驟是步驟1高光譜遙感圖像(1)高光譜遙感圖像(1)為標準格式的機載航空高光譜遙感反射率圖像；步驟2圖像光譜反射率的提取(2)a).高光譜遙感圖像(1)模塊的輸出信號傳送到圖像光譜反射率曲線(3)模塊；b).打開高光譜圖像，讀取高光譜圖像的基本信息其中設置波段之間寬度為bw納米，圖像存儲類型為整型數(shù)int、浮點數(shù)float或其他，像元值存儲字節(jié)長度nc，文件頭字節(jié)數(shù)no，圖像波段數(shù)nb，圖像的寬度上像元數(shù)wp和高度上像元數(shù)bp；c).計算反射率值，建立二維數(shù)組由圖像基本信息可以得到圖像第i個像元的第j個波段的反射率值，從文件頭偏移值ro由下式計算ro＝no+(wp*hp)*(j-1)*nc+i-1(1)建立二維數(shù)組rb[][]，rb維數(shù)分別為wp*hp和nb，從p＝1到p＝wp*hp，b＝1到b＝nb依次讀取每個像元每個波段的反射率值，并存儲在rb[p][b]中；d).二維數(shù)組rb[][]為提取出的圖像光譜反射率曲線；步驟3圖像光譜反射率曲線(3)圖像光譜反射率曲線(3)為步驟2中提取出的數(shù)組rb[][]，圖像光譜反射率曲線(3)模塊的輸出信號一路傳送到一次濾波后反射率曲線(9)模塊的輸入端，另一路傳送到二次濾波后反射率曲線(11)模塊的輸入端；步驟4反射率二階導數(shù)的計算(4)a).圖像光譜反射率曲線(3)模塊的輸出信號傳送到反射率曲線二階導數(shù)(5)模塊；b).用高光譜圖像光譜的二階導數(shù)作為光譜噪聲評定的變量其中光譜導數(shù)技術包括對反射光譜進行數(shù)學模擬和計算不同階數(shù)的導數(shù)值；確定光譜彎曲點及最大最小反射率的波長位置；光譜導數(shù)處理強調(diào)曲線的變化和壓縮均值影響；c).計算光譜的二階導數(shù)光譜的二階導數(shù)計算公式如下d2sdλ2|j≈s(λi)-2s(λj)+s(λk)(Δλ)2---(2)]]>式中Δλ＝λk-λj＝λj-λi，λk＞λj＞λi；λi為高光譜圖像第i個波段的波長，Δλ為高光譜圖像的波段之間寬度bw，s(λi)為像元反射率光譜曲線第i個波段的反射率值；d).具體計算二階導數(shù)值對圖像光譜反射率曲線數(shù)組rb[][]，對第p個像元其光譜二階導數(shù)計算時，第b個波段的二階導數(shù)值為(rb[p][b-1]-2*rb[p][b]+rb[p][b+1])/bw2(3)從p＝1到p＝wp*hp，b＝2到b＝nb-1依次進行光譜二階導數(shù)計算，并將計算結果保存在數(shù)組sd[][]中，維數(shù)分別為wp*hp和nb；在數(shù)組sd[p][b]中，當b＝1和b＝nb時，數(shù)組元素值為0；e).數(shù)組sd[][]為反射率曲線二階導數(shù)；步驟5反射率曲線二階導數(shù)(5)反射率曲線二階導數(shù)(5)為步驟4中提取出的數(shù)組sd[][]；步驟6高光譜圖像光譜域噪聲自檢測及反射率曲線噪聲判定系數(shù)的計算(6)a).反射率曲線二階導數(shù)(5)模塊的輸出信號傳送到反射率曲線噪聲判定系數(shù)(7)模塊；b).用高光譜圖像光譜二階導數(shù)的標準差作為光譜噪聲的判定條件其中光譜的低階導數(shù)處理對噪聲影響敏感性較低，高階微分對噪聲影響敏感度高；一階導數(shù)反映了反射率曲線的斜率，二階導數(shù)反映了噪聲的實際分布情況；c).用光譜二階導數(shù)對光譜曲線進行噪聲影響程度檢測1).進行光譜二階導數(shù)計算將二階導數(shù)值s″(λi)與其平均值μ之差和給定的某一閾值作比較；2).判定該波段噪聲的大小；3).確定該波段濾波平滑窗的大??；4).判斷不等式對某一波段i而言，如果不等式|s″(λi)-μ|＞σ(4)成立的話，那么該波段反射率值認為是有較大噪聲存在，反射率曲線噪聲判定系數(shù)為1，表明該波段噪聲較大，式中σ是光譜二階導數(shù)的標準差；如果不等式|s″(λi)-μ|≤σ (5)成立的話，那么該波段反射率值認為存在噪聲較小，反射率曲線噪聲判定系數(shù)為0；d).具體計算二階導數(shù)的平均值1).對反射率曲線二階導數(shù)數(shù)組sd[][]，第p個像元的反射率曲線二階導數(shù)的平均值msd為msd=(Σb=1nbsd[p][b])/nb---(6)]]>2).第p個像元的反射率曲線二階導數(shù)的標準差ssd為-3-ssd=[Σb=1nb(sd[p][b]-msd)]/nb---(7)]]>3).建立二維數(shù)組de[][]其維數(shù)分別為wp*hp和nb；4).判斷不等式對第p個像元的第b個波段而言，如果不等式|sd[p][b]-msd|＞ssd (8)成立的話，de[p][b]＝1；如果不等式|sd[p][b]-msd|≤ssd (9)成立的話，de[p][b]＝0；從p＝1到p＝wp*hp，b＝1到b＝nb依次進行反射率曲線噪聲判定系數(shù)計算，并將計算結果保存在數(shù)組de[][]中；e).數(shù)組de[][]為反射率曲線噪聲判定系數(shù)；步驟7反射率曲線噪聲判定系數(shù)(7)反射率曲線噪聲判定系數(shù)(7)為步驟6中計算得到的數(shù)組de[][]；步驟8噪聲的去除及圖像光譜反射率曲線噪聲的第一次濾波(8)a).反射率曲線噪聲判定系數(shù)(7)模塊的輸出信號傳送到一次濾波后反射率曲線(9)模塊；b).用賽維特斯基-高勒Savitzky-Golay平滑濾波器；c).用簡化的最小二乘擬合卷積方法對曲線進行平滑處理；d).計算平滑后曲線各階導數(shù)1).簡化后的最小二乘卷積方程式如下Yj=Σi=-mmCiyj+iN----(10)]]>式中，y是原始光譜值，Y是平滑后光譜值；Ci是平滑窗口中第i個光譜值的系數(shù)，N是卷積中點值個數(shù)，j是沿原始數(shù)據(jù)縱坐標數(shù)據(jù)列的計算點下標；2).平滑濾波法計算的卷積點為25點，并計算到光譜的第六階導數(shù)；3).各階導數(shù)平滑系數(shù)的計算公式其零階導數(shù)二次或三次多項式擬合的任意點數(shù)平滑窗系數(shù)計算公式如下ps(0)=3(3m2+3m-1-5s2)(2m+3)(2m+1)(2m-1),s=-m,...0,...m---(12)]]>式中，m為平滑窗的寬度的一半；e).檢測在應用光譜學中光譜處理方法之前，對光譜中誤差嚴重程度進行檢測；1).對噪聲嚴重的地方進行高點數(shù)的平滑；2).對噪聲較小部分應用較小平滑窗進行平滑；f).具體計算第一次濾波后圖像反射率曲線1).建立一維數(shù)組sg1[11]；2).用來存儲11點賽維特斯基-高勒濾波器系數(shù)；3).計算sg1[]的值按如下公式sg1[6+m]=3(3*52+3*5-1-5*m2)(2*5+3)(2*5+1)(2*5-1)m=-5,...0,...5---(13)]]>計算出sg1[]的值；4).建立二維數(shù)組rff1[][]其維數(shù)分別為wp*hp和nb；5)判斷如果第p個像元的第b個波段的反射率曲線噪聲判定系數(shù)de[p][b]＝1，那么對第p個像元的第b個波段反射率值進行光譜域濾波，按如下公式rf1[p][b]=Σi=-55sg1[6+i]*rb[p][b+i]11---(14)]]>式中5＜b＜nb-5；如果第p個像元的第b個波段的反射率曲線噪聲判定系數(shù)de[p][b]＝0，按如下公式rf1[p][b]＝rb[p][b](15)對圖像中每個像元逐波段進行計算；g).二維數(shù)組rf1[][]為第一次濾波后，去除了較大噪聲的圖像反射率曲線；步驟9一次濾波后反射率曲線(9)一次濾波后反射率曲線(9)為步驟8中計算出的二維數(shù)組rf[][]；步驟10噪聲的去除及圖像光譜反射率曲線噪聲的第二次濾波(10)a).一次濾波后反射率曲線(9)模塊的輸出信號傳送到二次濾波后反射率曲線(11)模塊；b).具體計算第二次濾波后圖像反射率曲線1).建立一維數(shù)組sg2[5]；2).用來存儲5點賽維特斯基-高勒濾波器系數(shù)；3).計算sg2[]的值按如下公式sg1[3+m]=3(3*22+3*2-1-2*m2)(2*2+3)(2*2+1)(2*2-1)m=-2,...0,...2---(16)]]>計算出sg2[]的值；4).建立二維數(shù)組rf2[][]其維數(shù)分別為wp*hp和nb；5)判斷對第p個像元的第b個波段反射率值進行光譜域濾波，按如下公式rf2[p][b]=Σi=-55sg1[3+i]*rf1[p][b+i]5---(17)]]>式中2＜b＜nb-1；當b＝1，2，nb-1，nb時，按如下公式rf2[p][b]＝rb[p][b](18)對圖像中每個像元逐波段進行計算；c).二維數(shù)組rf2[][]為第二次濾波后的圖像反射率曲線；步驟11二次濾波后反射率曲線(11)二次濾波后反射率曲線(11)為步驟10中得到的二維數(shù)組rf2[][]；步驟12由濾波后反射率曲線重構高光譜圖像(12)a).二次濾波后反射率曲線(11)模塊的輸出信號傳送到光譜域濾波后高光譜遙感圖像模塊；b).建立新文件按航空高光譜圖像的格式標準，高光譜圖像文件大小為(no+wp*hp*nb)*nc個字節(jié)，建立新文件；c).寫入新文件將步驟1中讀取出的高光譜圖像文件頭的no*nc個字節(jié)，寫入新文件中；d).寫入對p＝1，...，wp*hp，依次寫入rf2[p][1]，然后再從p＝1，...，wp*hp，依次寫入rf2[p][2]，直到從p＝1，...，wp*hp，依次寫入rf2[p][nb]；e).重構完成高光譜圖像文件重構完成；步驟13光譜域濾波后高光譜遙感圖像(13)光譜域濾波后高光譜遙感圖像(13)為步驟12中重構的高光譜圖像文件。
全文摘要
一種涉及航空高光譜遙感圖像光譜域噪聲自檢測與去除方法，尤指一種針對國產(chǎn)高光譜遙感器的高光譜圖像進行處理，可以檢測并去除高光譜圖像光譜域中噪聲，屬工程科學技術中的遙感技術應用領域。該方法是通過一種適用于光譜噪聲評定的變量、判定適用于光譜噪聲評定的條件及找到適用于高光譜圖像光譜域濾波器的方法，實現(xiàn)檢測區(qū)分高光譜圖像光譜域噪聲或細微光譜特征，對光譜域較大噪聲去除；主要解決如何找到適用于高光譜圖像光譜域濾波器的方法系統(tǒng)軟件及有關硬件等技術問題。本發(fā)明的優(yōu)點是能有效的去除光譜中存在噪聲，并能保留光譜原有的大部分光譜特征，是地物實測光譜和高光譜圖像預處理的一種有效手段。
文檔編號G01S17/00GK1710445SQ200510027528
公開日2005年12月21日 申請日期2005年7月5日 優(yōu)先權日2005年7月5日
發(fā)明者王強, 束炯 申請人:華東師范大學