去除DMSP/OLS影像模糊的RTSVD算法

2019-08-07 02:04:46于志文鄧開元畢旋旋

測繪通報 2019年7期

于志文，任超,2，鄧開元，畢旋旋

(1. 桂林理工大學測繪地理信息學院，廣西桂林 541004; 2. 廣西空間信息與測繪重點實驗室，廣西桂林 541004)

DMSP/OLS(the defense meteorological satellite program operational linescan system)影像是一種應用廣泛的夜光遙感數據，在社會經濟數據挖掘方面具有重要價值。該影像容易獲取且可重復監測，但不當的數據采集與管理機制導致影像模糊，一般通過選取閾值或結合其他數據來抑制模糊。文獻[1]提出89%累計發光閾值濾波，有效緩解了燈光溢出問題,但僅針對城市面積提取，且丟失了很多細節信息。文獻[2—3]通過對DMSP/OLS影像設定3種不同閾值來獲取城市邊界，得到了比較精確的城市邊界，但是處理過程過于煩瑣。文獻[4]利用DMSP/OLS輻射定標夜間燈光強度數據校正,有效擴展了強度值的動態范圍，但并未改善模糊現象。文獻[5]將NDVI和地表溫度結合起來，在削弱城市中心DN值飽和的同時緩解了模糊現象，但輔助數據需經過復雜的計算。文獻[6]提出了Pct影像濾波方法。Pct影像容易獲取且結合Pct影像濾波的去模糊處理方法便于實現，但僅依靠Pct影像濾波不足以消除燈光邊緣部分的模糊及噪聲，還需結合其他方法。

綜上，本文提出一種結合Pct影像濾波處理的正則化截斷奇異值分解(RTSVD)去模糊算法。首先，根據光源像素發光頻率必定高于非光源像素的規律，利用Pct數據比較像素發光頻率，關閉圖像中所有非光源像素，并設定15%發光頻率閾值去除偶然影響。然后，利用高斯分布擬合點擴散函數(PSF)，并用自反邊界條件[7]構造模糊矩陣。最后，繪制L曲線計算截斷參數，采用RTSVD法對圖像進行去模糊處理。

1 DMSP/OLS數據介紹

DMSP/OLS數據始于20世紀70年代美國國防部開展的美國軍事氣象衛星計劃。DMSP衛星系統由6顆衛星組成，采用雙星運行機制，兩顆衛星同時運行。衛星運行在高度為833 km的太陽同步軌道上，運行周期為101 min，掃描寬幅為3000 km，每天收集14個軌道信息，可提供4次全球覆蓋。OLS傳感器是DMSP衛星的主要傳感器，擁有可見光近紅外(VNIR:0.4～1.1 μm)和熱紅外(TIR:10～13 μm)兩個通道[8]，近地點分辨率可達500 m，掃描帶邊緣分辨率大于1 km。

本文以平均燈光影像和Pct影像兩種DMSP/OLS數據為研究對象，數據來源于美國國家地球物理數據中心。DMSP/OLS非輻射定標夜間燈光強度數據包括無云觀測頻數影像、平均燈光影像、穩定燈光影像及Pct影像。其中，平均燈光影像是未經濾波處理的可見光強度平均值影像，是模糊問題的主要研究對象；穩定燈光影像是去除偶然噪聲且經過降噪濾波的影像，圖像質量優于平均燈光影像，可以用于影像去模糊效果評估；Pct影像是無云觀測強度與燈光探測百分比的乘積，可用于監測燈光的發光頻率。

2 平均燈光影像模糊原因

2.1 數據采集部分

DMSP衛星正下方的視場是正圓[8]，掃描運動使視場拉伸為橢圓。位于視場質心處的精細像素大小為0.56 km×0.56 km，遠小于視場面積。在視場移動過程中，同一固定光源會落在多個重疊的視場范圍內，形成以光源為中心的視場范圍大小的模糊橢圓，并且由于衛星誤算，模糊橢圓的質心存在隨機偏移，導致DMSP影像中光源的地理位置偏移。文獻[8]通過在全球不同地區移動便攜式光源，測算出每夜DMSP影像中光源位置與GPS觀測位置的距離與方位，所測的28個點全部向北偏移，偏移量均值為2.9 km,大約為一個平滑像素的寬度。文獻[5]進一步推測出偏移的標準差為1～1.12 km。

2.2 數據管理部分

DMSP衛星的技術有限，不能存儲浮點型小數，因此每個精細像素都是取整后的結果。同時，DMSP衛星也不能直接存儲精細像素的高分辨率信息，因此存儲時對精細像素進行平滑，把5×5像素的精細像素作為一個平滑像素。但DMSP衛星對精細像素平滑時每夜的起始計算位置并不統一，導致每夜的平滑像素圖像不會完全重疊，這也是年平均影像模糊的原因之一。

3 基于Pct影像濾波的RTSVD算法

3.1 Pct影像濾波

Pct影像是由無云探測的平均DN值乘以探測百分比頻率得到的，可以代表燈光的百分比頻率。依據燈光影像成像原理，光源像素的發光頻率必定高于非光源像素，可利用Pct影像關閉發光頻率小于周圍像素的像元，以去除非光源像素誤差。

首先篩選出像素xi,j周圍發光頻率最高的像素xmax，比較xi,j與xmax為

(1)

把經過式(1)處理后的Pct影像與模糊影像相乘，得到光源像素影像。但Pct影像濾波不能完全關閉燈光邊緣部分的模糊像素，而且無法去除影像中的噪聲。因此，本文在Pct影像濾波的基礎上提出RTSVD算法進一步處理平均燈光影像。

3.2 RTSVD算法

RTSVD算法是融合了正則化思想的截斷奇異值分解法，目的在于去除影像中的噪聲并恢復圖像質量。對于m行n列的模糊圖像，一般有線性模糊模型

Ax=b

(2)

式中，x、b分別為真實圖像與模糊圖像排列成的列向量，長度為N=mn；A為模糊矩陣，A∈RN×N。由此可以得出基礎去模糊模型

x=A-1b

(3)

模糊矩陣A由點擴散函數與邊界條件兩個成分決定。其中，點擴散函數可以利用高斯分布來擬合[7，9]。根據模糊原因，模糊矩陣A的邊界條件設為自反邊界條件[10]。然后，利用奇異值分解進一步完善去模糊模型。對于m×n階模糊矩陣A，存在一個分解使得

A=USVT

(4)

式中，U和V為正交矩陣，分別為左奇異向量陣和右奇異向量陣，滿足UTU=Im，VTV=In；S為奇異值矩陣，為對角矩陣，σ1、σ2、…、σm為對角線上的奇異值，滿足

σ1≥σ2≥…≥σm≥0

(5)

式(3)可寫成

(6)

式中，ui為左奇異向量；vi為右奇異向量。模糊矩陣A為病態矩陣，條件數cond(A)=σ1/σn非常大，奇異值衰減的速度非?？?，如圖1所示。

與小奇異值對應的誤差分量被放大，干擾圖像真實信息，公式描述如下[11]

(7)

式中，Δ為誤差。為了削弱誤差Δ的干擾，運用RTSVD算法，用截取的奇異值矩陣Sk來替代原奇異值矩陣Sm，公式如下

(8)

式中，k為截斷參數。RTSVD算法的誤差來源于正則化與噪聲兩個方面，在選擇截斷參數時，應考慮在抑制噪聲的同時盡可能縮小正則化誤差。

3.3 L曲線

截斷參數的選擇對于RTSVD算法至關重要，一般采用GCV或L曲線來確定截斷參數。但針對圖像模糊問題時，GCV難以計算出截斷誤差的最小值。因此，本文采用L曲線確定截斷參數，利用三次樣條曲線擬合離散點。L曲線是以log-log為尺度，在以正則化解范數為縱軸，對應的殘差范數為橫軸的直角坐標系中，對一系列截斷參數繪制的曲線[12]。L曲線上曲率最大的點對應的截斷參數為最優解。令ρ、η分別為殘差與解的對數，L曲線的曲率K公式表達為

(9)

在截斷奇異值分解中，L曲線是由若干離散點組成的，要計算出最優解，必須在保證曲線整體形狀不變的條件下定義一條與離散點相關聯的可微、光滑的曲線。因此，采用三次樣條曲線逼近L曲線離散點,計算連續曲線上曲率最大的點從而確定截斷參數。

4 算例分析

為了驗證RTSVD算法對DMSP圖像去模糊的有效性,裁剪大小為417×223像素的臺灣地區F142000期平均燈光影像進行去模糊試驗。為了評價本文方法效果，選取穩定燈光影像和89%累計發光閾值濾波影像與試驗結果進行比較,如圖2所示。

除了對試驗結果進行目視觀察評價外，引入峰值信噪比(PSNR)、圖像信息熵、灰度平均梯度(GMG)與邊緣強度作為客觀評價標準。由表1可以看出，原始的平均燈光影像(圖2(a))存在嚴重的模糊現象。依據L曲線選取截斷參數k=67，處理得到圖像(圖2(d))與原始影像相比，峰值信噪比增大2.505 3，信息熵增加0.197 4，灰度平均梯度增加1.023 9，邊緣強度增大9.855 3，與穩定燈光影像相比，各項指標均有所增長。對于直接頻率濾波處理的影像(圖2(c))，峰值信噪比提高3.443 5，信息熵降低0.482 8，表明該方法雖然有效降低了噪聲與模糊，但嚴重丟失了影像的細節信息。這也說明結合Pct影像濾波的RTSVD算法在消除模糊的同時保留了影像的細節信息。