基于總變分低秩組稀疏的全球雷達(dá)數(shù)據(jù)修復(fù)算法

葛晨宇，董良，許伊昆，常毅，張宏鳴

（西北農(nóng)林科技大學(xué)信息工程學(xué)院，陜西楊凌 712100）

0 引言

航天飛機(jī)雷達(dá)地形測(cè)繪任務(wù)（Shuttle Radar Terrain Mission，SRTM）通過(guò)離散高程點(diǎn)采集來(lái)數(shù)字化模擬地形表面［1］，在過(guò)去十年中被廣泛應(yīng)用于地形建模、地理信息系統(tǒng)和數(shù)據(jù)計(jì)算等領(lǐng)域中；然而，由于衛(wèi)星在數(shù)據(jù)采集過(guò)程中雷達(dá)干涉儀桅桿的殘余運(yùn)動(dòng)［2］、地形表面反射率的突然變化［3］等因素，生成后的數(shù)據(jù)中存在混合分布的隨機(jī)噪聲和系統(tǒng)誤差，具體表現(xiàn)為尖峰、斑點(diǎn)和多向條紋誤差。這類誤差的存在對(duì)數(shù)據(jù)的可靠性和后續(xù)應(yīng)用造成嚴(yán)重影響。

國(guó)內(nèi)外研究者先后基于該問(wèn)題進(jìn)行了數(shù)據(jù)修復(fù)方法研究，主要方法可分為4類：基于空間域的方法［4］、基于域的組合方法［5］、基于統(tǒng)計(jì)的方法［6-7］和基于傅里葉變換的方法［8-9］。基于空間域的方法通過(guò)計(jì)算空間變化的噪聲方差構(gòu)造數(shù)據(jù)窗口大小不同的濾波器去除數(shù)據(jù)中的隨機(jī)噪聲，但全局平滑操作將不可避免地影響局部非噪聲分量，使后續(xù)應(yīng)用與計(jì)算出現(xiàn)與實(shí)際不符的現(xiàn)象，同時(shí)這種類型的方法不適用于消除系統(tǒng)誤差；基于域的組合方法通過(guò)交叉濾波來(lái)實(shí)現(xiàn)去條紋，但該類方法在消除條帶誤差的同時(shí)也會(huì)導(dǎo)致條帶方向上的數(shù)據(jù)丟失，使其地形表面出現(xiàn)斷層現(xiàn)象，并且該方法針對(duì)條帶誤差的不同分量需要手動(dòng)修改不同的內(nèi)核函數(shù)；基于統(tǒng)計(jì)的方法通過(guò)將實(shí)地采集的真實(shí)高程數(shù)據(jù)與目標(biāo)數(shù)據(jù)集進(jìn)行校準(zhǔn)匹配，獲得了誤差的先驗(yàn)信息從而對(duì)目標(biāo)數(shù)據(jù)進(jìn)行修復(fù)，但顯然對(duì)于全球范圍，這種適用于局部區(qū)域的修復(fù)方法不再適用；基于傅里葉變換的方法分離數(shù)據(jù)中的低頻和高頻分量，目標(biāo)是去除錯(cuò)誤信息的頻率分量，但當(dāng)實(shí)際數(shù)據(jù)中的錯(cuò)誤信息具有復(fù)雜多樣的頻率成分時(shí)，數(shù)據(jù)的修復(fù)效果存在問(wèn)題。Gallant［10］對(duì)上述方法進(jìn)行綜合應(yīng)用，利用基于自適應(yīng)濾波與傅里葉變換的方法修復(fù)了澳大利亞區(qū)域數(shù)據(jù)，但由于上述提到的方法局限性，其產(chǎn)品手冊(cè)［11］中提到部分修復(fù)后數(shù)據(jù)存在過(guò)平滑與誤差殘留現(xiàn)象，這些問(wèn)題會(huì)對(duì)數(shù)據(jù)的進(jìn)一步使用造成影響。

近年來(lái)，隨著圖像修復(fù)技術(shù)的發(fā)展，基于優(yōu)化的方法［12-16］被提出，該類方法將誤差消除視為逆向修復(fù)問(wèn)題，通過(guò)求解正則化模型來(lái)進(jìn)行誤差消除，如改進(jìn)的單向總變差模型、低秩模型和全局稀疏模型。這些方法在醫(yī)學(xué)圖像、高光譜圖像和自然圖像等領(lǐng)域產(chǎn)生了出色的誤差消除效果，但是仍存在一定問(wèn)題，如由于過(guò)多地關(guān)注數(shù)據(jù)的非局部方面和單一噪聲的特性，在面對(duì)多類混合誤差時(shí)存在其他類噪聲殘留的現(xiàn)象。因此，數(shù)據(jù)恢復(fù)的效果更多地取決于任務(wù)中誤差的復(fù)雜程度，這導(dǎo)致在真實(shí)數(shù)據(jù)中的復(fù)雜情況下，如非正則方向條紋誤差與隨機(jī)噪聲混合時(shí)，誤差消除效果并不優(yōu)異。

對(duì)于數(shù)據(jù)修復(fù)或去噪算法的評(píng)價(jià)，先前研究更多關(guān)注于數(shù)據(jù)全局結(jié)構(gòu)相似度以及噪聲水平此類宏觀評(píng)價(jià)，而忽視了能夠反映地形數(shù)據(jù)細(xì)節(jié)和地形特征表達(dá)能力的相關(guān)指標(biāo)。在基于高程數(shù)據(jù)的后續(xù)評(píng)估［17-18］中發(fā)現(xiàn)，目前宏觀評(píng)價(jià)優(yōu)異的航天飛機(jī)雷達(dá)地形測(cè)繪任務(wù)（SRTM）、MERIT 數(shù)字高程模型（MERIT-Digital Elevation Model，MERIT-DEM）等數(shù)據(jù)集在坡度（Slope）、坡向（Aspect）的計(jì)算中表現(xiàn)不佳，而坡度、坡向等與數(shù)據(jù)一階導(dǎo)數(shù)相關(guān)的地形因子作為局部細(xì)節(jié)信息的體現(xiàn)，擁有對(duì)局部數(shù)據(jù)變化的強(qiáng)敏感性，是評(píng)價(jià)數(shù)據(jù)的合理指標(biāo)；因此本文考慮將坡度、坡向作為數(shù)據(jù)修復(fù)算法后續(xù)應(yīng)用方面的評(píng)價(jià)指標(biāo)。

考慮到數(shù)據(jù)中細(xì)節(jié)信息的恢復(fù)效果，本文研究從優(yōu)化角度出發(fā)，分析了局部混合誤差的固有特征并構(gòu)建誤差的低秩稀疏混合錯(cuò)誤模型，同時(shí)利用變分思想避免非誤差分量的消除，最后利用凸優(yōu)化算法對(duì)模型進(jìn)行求解保證收斂性從而實(shí)現(xiàn)消除數(shù)據(jù)中的混合誤差。同時(shí)，在此基礎(chǔ)上，結(jié)合一階地形因子的計(jì)算原理，考慮算法中單方向約束對(duì)數(shù)據(jù)修復(fù)的影響，對(duì)數(shù)據(jù)梯度方向進(jìn)行總變分（Total Variation，TV）正則約束，使數(shù)據(jù)局部范圍變化差值趨于正常，減少由單方向變分導(dǎo)致的局部數(shù)據(jù)值銳利過(guò)度，使得地貌邊緣與實(shí)際數(shù)據(jù)更加相符。實(shí)驗(yàn)結(jié)果表明，與現(xiàn)有處理方法對(duì)數(shù)據(jù)進(jìn)行處理的結(jié)果相比較，本文采用的基于TV 約束的低秩組稀疏（Low-Rank Group Sparsity_Total Variation，LRGS_TV）算法在宏觀表達(dá)和局部特征評(píng)估中效果更好。本文的主要工作如下：

1）從基于優(yōu)化思想的角度出發(fā)，分析SRTM 數(shù)據(jù)中混合誤差的固有特征，構(gòu)建混合誤差模型；

2）考慮到SRTM 數(shù)據(jù)后續(xù)地形因子計(jì)算，對(duì)混合誤差模型進(jìn)行改進(jìn)，改善局部數(shù)據(jù)銳利過(guò)度的問(wèn)題；

3）利用交替方向乘子優(yōu)化對(duì)低秩組稀疏模型進(jìn)行求解，保證了模型的收斂性。

1 總變分約束低秩組稀疏算法

高程模型是衛(wèi)星雷達(dá)收集數(shù)據(jù)生成的離散點(diǎn)云，并用空間網(wǎng)格結(jié)構(gòu)來(lái)表示實(shí)際地形的空間分布［19］，后續(xù)數(shù)據(jù)計(jì)算以及場(chǎng)景建模都基于此空間結(jié)構(gòu)進(jìn)行，因此，可以將SRTM 數(shù)據(jù)每個(gè)弧秒的局部數(shù)據(jù)網(wǎng)格近似認(rèn)為是以經(jīng)度和緯度為坐標(biāo)、以高度為坐標(biāo)值的矩陣，據(jù)此可構(gòu)建數(shù)據(jù)的局部數(shù)據(jù)結(jié)構(gòu)。

1.1 數(shù)據(jù)結(jié)構(gòu)

給定數(shù)據(jù)中的經(jīng)度和緯度坐標(biāo)（X和Y）以及相應(yīng)的離散點(diǎn)高程，則高程數(shù)據(jù)集可以表示為：

其中：R是數(shù)據(jù)集；i和j代表經(jīng)緯度網(wǎng)格中的行數(shù)和列數(shù)；Xi和Yj是相應(yīng)網(wǎng)格的經(jīng)度和緯度；EXi，Yj是相應(yīng)經(jīng)緯度處的混合高程值。

基于數(shù)據(jù)中混合誤差的疊加形式，可將包含隨機(jī)噪聲和系統(tǒng)誤差的結(jié)構(gòu)模型描述為：

其中：E是包含混合誤差的高程數(shù)據(jù)；T是真實(shí)高程值；e(T(xi，yj))是相應(yīng)經(jīng)緯度處的混合誤差分量。

而混合誤差中多向條帶誤差與T在條帶的切線方向密切相關(guān)，因此混合錯(cuò)誤可以表示為：

其中：S(T(xi，yj)) 為多向條紋誤差，N為隨機(jī)噪聲?；谶@一數(shù)據(jù)結(jié)構(gòu)，本文研究技術(shù)路線如圖1所示，主體分為3部分：低秩變換框架、混合誤差剔除算法與優(yōu)化器系統(tǒng)。

圖1 全球雷達(dá)數(shù)據(jù)修復(fù)算法技術(shù)路線Fig.1 Technical roadmap of global-scale radar data restoration algorithm

首先，利用低秩稀疏逼近提取局部混合誤差矩陣，搜索誤差矩陣最低秩時(shí)的方向角，從而使用低秩紋理映射算法獲取混合誤差中條帶結(jié)構(gòu)最低秩方向時(shí)的變換因子構(gòu)建方向變換框架，而不直接進(jìn)行全局處理。這是因?yàn)闂l帶結(jié)構(gòu)分量的低秩特性在正則方向時(shí)，其低秩和稀疏特性最為優(yōu)異［20］，易于約束優(yōu)化。

其次，利用數(shù)據(jù)在局部范圍低秩方向上的唯一性，正則化全局多方向條帶錯(cuò)誤結(jié)構(gòu)并使用變分思想進(jìn)行單向約束；同時(shí)使用加權(quán)核范數(shù)的非局部自相似性來(lái)消除隨機(jī)噪聲，結(jié)合TV 正則對(duì)梯度進(jìn)行約束，使數(shù)據(jù)局部范圍變化差值趨于正常，數(shù)據(jù)細(xì)節(jié)與實(shí)際更加相符，從而減少對(duì)后續(xù)數(shù)據(jù)計(jì)算的影響。

最后，由于模型中正則化項(xiàng)不可微且不可分割，很難直接從恢復(fù)模型中求解恢復(fù)數(shù)據(jù)［21］。為了解決這個(gè)問(wèn)題，設(shè)計(jì)基于交替方向乘子法優(yōu)化器，保證了算法的有效性與模型的收斂性［22］。

1.2 LRGS_TV數(shù)據(jù)修復(fù)算法

梯度對(duì)數(shù)據(jù)中高度值的突然變化非常敏感［23］，這表明基于高程數(shù)據(jù)計(jì)算得到的一階或二階導(dǎo)數(shù)相關(guān)指標(biāo)將受到梯度數(shù)據(jù)異常的影響，且由于數(shù)據(jù)中混合誤差具有空間自相關(guān)性，因此全局范圍的多方向條帶在局部區(qū)域表現(xiàn)為單一方向，主要影響與條帶方向垂直方向上的梯度數(shù)據(jù)。

根據(jù)這一特性，對(duì)局部條帶進(jìn)行低秩紋理映射后正則方向條帶將會(huì)對(duì)水平方向數(shù)據(jù)造成破壞，則數(shù)據(jù)的水平梯度以及正則方向的條帶誤差‖ ?y(S°τ)‖1的L1稀疏表示可作為數(shù)據(jù)與條帶誤差自身的約束項(xiàng)，使用迭代收縮閾值算法（Iterative Shrinkage Thresholding Algorithm，ISTA）［24］可對(duì)其進(jìn)行優(yōu)化求解，如式（4）～（5）所示：

其中：C是拉格朗日乘數(shù)；k是正向懲罰參數(shù)；λ是正則化參數(shù)；軟閾值操作函數(shù)為soft(x，y)=sign(x)max{x-y，0}。

由于衛(wèi)星雷達(dá)在飛行過(guò)程中的偶發(fā)異常振動(dòng)，導(dǎo)致數(shù)據(jù)局部區(qū)域通常會(huì)出現(xiàn)多種誤差混合現(xiàn)象，因此對(duì)后續(xù)計(jì)算影響情況也更為復(fù)雜。條帶誤差是由于衛(wèi)星雷達(dá)干涉儀桅桿的殘余運(yùn)動(dòng)引起的高度規(guī)則波動(dòng)，且條帶誤差主要集中在東北到西南和西北到東南方向。根據(jù)這一特性可發(fā)現(xiàn)，低秩紋理映射后正則方向的條帶具有規(guī)律的分布，因此在局部區(qū)域內(nèi)，在排除隨機(jī)噪聲的干擾情況下，條帶誤差具有組稀疏特征，可約束表示為，使用軟閾值算法進(jìn)行求解［24］。Sl是局部區(qū)域排除隨機(jī)噪聲干擾的正則方向條帶誤差，可表示為S°τ=Sl+N，其中‖·‖2，1是L2，1范 數(shù)，，τ是變換因子。

同時(shí)，在無(wú)隨機(jī)噪聲干擾情況下，干涉儀引起的高程變化在局部區(qū)域具有相似的高程值，因此條帶誤差可視為系統(tǒng)誤差，可以利用誤差的低秩特性與常規(guī)數(shù)據(jù)區(qū)分開(kāi)，可用核范數(shù)‖Sl‖*表示，然后使用奇異值閾值（Singular-Value Threshold，SVT）［25］算法進(jìn)行求解，如式（6）：

其中：UσΣV*是矩陣的奇異值分解；U和V為正交矩陣；σΣ是奇異值對(duì)角矩陣。尖峰誤差通常分布在地形平緩區(qū)域，主要是因?yàn)槠骄弲^(qū)域地表反射率的突然變化會(huì)導(dǎo)致隨機(jī)噪聲的產(chǎn)生，這一特性表明尖峰誤差可視為隨機(jī)噪聲從而使用非局部相似性進(jìn)行分離。隨機(jī)噪聲對(duì)局部區(qū)域數(shù)據(jù)的低秩性會(huì)造成干擾，這對(duì)該區(qū)域條帶誤差的低秩性檢測(cè)造成不利影響，同時(shí)隨機(jī)噪聲在兩個(gè)梯度方向上都干擾了數(shù)據(jù)，這同樣是導(dǎo)致基于梯度的后續(xù)數(shù)據(jù)計(jì)算出現(xiàn)異常的重要因素，因此將其稀疏表示約束為‖S°τ-Sl‖1。

根據(jù)以上特征約束得到的低秩組稀疏模型可對(duì)數(shù)據(jù)中存在的混合誤差實(shí)現(xiàn)有效消除，但在對(duì)于基本地形特征如坡向、坡度的表現(xiàn)中卻仍然不盡人意。探索坡度坡向公式［26］可得坡度坡向的計(jì)算定義如下。

坡向（Aspect）是地表上一點(diǎn)切平面的法線矢量在水平面的投影與過(guò)該點(diǎn)正北方向的夾角，在z=f(x，y)曲面上其計(jì)算式為式（7）：

坡度（Slope）是地表任意一點(diǎn)切平面與水平面的夾角，通常簡(jiǎn)化的差分公式如式（8）：

由坡度坡向公式可得，雖然上述模型中采用單梯度方向約束‖?xT‖1保證了數(shù)據(jù)更多細(xì)節(jié)保留，但另一梯度方向的忽略，導(dǎo)致了地貌邊緣梯度突兀，使出現(xiàn)算法模型引入的誤差，從而影響局部區(qū)域地形特征的計(jì)算。因此本文在單梯度方向約束消除混合噪聲之后，進(jìn)行TV 約束V(E)，使梯度不統(tǒng)一導(dǎo)致的局部高程變化躍變現(xiàn)象減輕，減少對(duì)后續(xù)計(jì)算的影響。其中：

綜合上述各特征約束項(xiàng)，將基于TV約束的低秩組稀疏模型表示為式（10）：

使用交替方向乘子算法進(jìn)行優(yōu)化求解，偽代碼如下。

算法1 基于ADMM優(yōu)化的修復(fù)算法。

1.3 評(píng)估指標(biāo)

定量評(píng)估指標(biāo)采用峰值信噪比（Peak Signal-to-Noise Ratio，PSNR）、結(jié)構(gòu)相似性（Structural SIMilarity，SSIM）指數(shù)以及通過(guò)數(shù)據(jù)計(jì)算所得坡度、坡向［27-28］。坡向信息能反映出地形表面局部的傾斜方向，不同的數(shù)值代表了對(duì)應(yīng)的方向，表示0°～360°的方向角范圍。坡度信息能反映出地形表面局部的傾斜角度，坡度值的范圍是0°～90°。

原始高程與修復(fù)后高程之差的F 范數(shù)的平方是原始數(shù)據(jù)與恢復(fù)的數(shù)據(jù)之間的均方誤差，因此適合本文研究數(shù)據(jù)的PSNR評(píng)估指標(biāo)為式（11）：

其中：MSE代表原始數(shù)據(jù)E和恢復(fù)后數(shù)據(jù)T之間高程的均方誤差；MAXE代表原始數(shù)據(jù)中的最大高程。PSNR值越大，數(shù)據(jù)失真越小。

本文實(shí)驗(yàn)的SSIM設(shè)定為式（12）：

其中：

對(duì)于不同算法的比較，除數(shù)據(jù)修復(fù)效果的PSNR、SSIM 以及坡度、坡向等定量評(píng)估，還需要考慮不同算法的時(shí)間復(fù)雜度，將在2.2節(jié)中進(jìn)行討論。

2 實(shí)驗(yàn)與結(jié)果分析

2.1 航天飛機(jī)雷達(dá)數(shù)據(jù)與實(shí)驗(yàn)設(shè)計(jì)

如圖2所示，實(shí)驗(yàn)使用數(shù)據(jù)集為航天飛機(jī)雷達(dá)地形測(cè)繪任務(wù)SRTM1（https：//e4ftl01.cr.usgs.gov/MEASURES/SRTMGL1.003/2000.02.11/），數(shù)據(jù)范圍在南緯56°到北緯61°，數(shù)據(jù)精度為1″（空間分辨率30 m×30 m）。數(shù)據(jù)集共分為14520 個(gè)區(qū)域文件，單文件包含3601×3601 個(gè)采樣點(diǎn)的高度數(shù)據(jù)［18］與對(duì)應(yīng)經(jīng)緯度坐標(biāo)。

圖2 SRTM數(shù)據(jù)與實(shí)驗(yàn)樣區(qū)Fig.2 SRTM data and experimental sample areas

現(xiàn)有研究表明，隨機(jī)噪聲和條紋誤差主要集中在平緩地形區(qū)域［29］，因平緩區(qū)域流向以及坡度表現(xiàn)不顯著，不利于對(duì)數(shù)據(jù)修復(fù)效果進(jìn)行定量評(píng)價(jià)，因此在地形細(xì)節(jié)豐富區(qū)域增加設(shè)置模擬實(shí)驗(yàn)進(jìn)行算法評(píng)估。

將LRGS_TV 與幾種主流降噪方法進(jìn)行比較，分別為：1）總變分（TV）［30］。利用含噪數(shù)據(jù)的總變分比無(wú)噪數(shù)據(jù)的總變分大這一特性，使數(shù)據(jù)的總變化在涉及噪聲的約束條件下被最小化。2）單方向總變分（Unidirectional TV，UTV）［31］。在誤差特征約束方面選擇比全局?jǐn)?shù)據(jù)更可靠的保真度項(xiàng)，將約束集中在條帶誤差的梯度方向，保留原始數(shù)據(jù)沒(méi)有被干擾的梯度方向不處理。3）低秩單圖像分解（Low-Rank-based Single-Image Decomposition，LRSID）［32］。對(duì)條帶的結(jié)構(gòu)特征進(jìn)行詳細(xì)分析，以期望將原始數(shù)據(jù)與條帶誤差完美分離。4）本文采用的低秩組稀疏（Low-Rank Group Sparsity，LRGS）模型［33］。

在模擬實(shí)驗(yàn)中，分別在各地區(qū)選取地形紋理細(xì)節(jié)豐富的N02E13（非洲）、N38W83（北美洲）、N57E123（亞歐大陸）、S36E148（澳洲）區(qū)域作為樣區(qū)1、2、3、4，將該四樣區(qū)數(shù)據(jù)作為模擬實(shí)驗(yàn)的相對(duì)真實(shí)值數(shù)據(jù)，在其中添加隨機(jī)寬度和隨機(jī)分布的傾斜條帶誤差并混合隨機(jī)強(qiáng)度全局隨機(jī)噪聲作為混合誤差干擾后數(shù)據(jù)，模擬實(shí)際中可能出現(xiàn)的極端情況。設(shè)置模擬實(shí)驗(yàn)有利于進(jìn)行不同算法對(duì)混合誤差干擾后數(shù)據(jù)處理對(duì)比，即向原始數(shù)據(jù)中添加混合誤差生成退化后數(shù)據(jù)進(jìn)行誤差去除實(shí)驗(yàn)時(shí)，原始數(shù)據(jù)就可以作為相對(duì)真實(shí)值數(shù)據(jù)與被各個(gè)降噪方法處理后的數(shù)據(jù)進(jìn)行PSNR、SSIM 以及坡向、坡度的定量對(duì)比評(píng)估以驗(yàn)證算法有效性。

在真實(shí)實(shí)驗(yàn)中，選取數(shù)據(jù)集中存在實(shí)際混合誤差的S35W62（南美洲）區(qū)域作為樣區(qū)5，其中包含局部小范圍傾斜條紋、全局大范圍傾斜條紋和隨機(jī)噪聲。

2.2 算法修復(fù)效果比較與分析

2.2.1 模擬實(shí)驗(yàn)

從局部放大中可以發(fā)現(xiàn)，樣區(qū)1 中的條紋誤差與隨機(jī)噪聲嚴(yán)重影響了數(shù)據(jù)信息的表達(dá)，如圖3（b）所示。對(duì)添加誤差后的數(shù)據(jù)使用幾種主流降噪方法進(jìn)行了處理，從細(xì)節(jié)方面進(jìn)行評(píng)估，圖3（c）中的TV 對(duì)于數(shù)據(jù)中的隨機(jī)噪聲消除優(yōu)異，但明顯出現(xiàn)地形細(xì)節(jié)丟失現(xiàn)象，這主要是因?yàn)門V無(wú)法對(duì)具體特征進(jìn)行約束，而全局統(tǒng)一平滑處理將不可避免產(chǎn)生過(guò)平滑的結(jié)果，同樣因?yàn)闊o(wú)法對(duì)特定特征進(jìn)行約束，其處理后數(shù)據(jù)中條紋誤差存在一定程度的殘留。UTV對(duì)于條帶誤差以及隨機(jī)噪聲都有一定程度的抑制作用但效果有限，這主要因?yàn)閮H進(jìn)行單梯度方向約束，對(duì)于特定方向的誤差消除有一定效果，但對(duì)于多方向混合誤差無(wú)法完全消除，且單向強(qiáng)約束對(duì)圖像的特定方向結(jié)構(gòu)產(chǎn)生了影響（圖3（d））。LRSID抑制了部分隨機(jī)噪聲的表達(dá)且對(duì)于數(shù)據(jù)結(jié)構(gòu)影響較小，但對(duì)于多方向混合誤差效果欠佳（圖3（e））。LRGS 在消除混合誤差方面顯示出更好的結(jié)果，且能夠保留較多的地形細(xì)節(jié)（圖3（f）），但與原數(shù)據(jù)相比，在地形中仍存在一定程度的扭曲現(xiàn)象，這主要是因?yàn)閱翁荻确较虻募s束雖然對(duì)數(shù)據(jù)細(xì)節(jié)的保留起到了一定程度的作用，但對(duì)另一梯度方向的忽略，導(dǎo)致了地貌邊緣梯度過(guò)于突兀，這一問(wèn)題對(duì)后續(xù)計(jì)算應(yīng)用將存在不利影響。圖3（g）中的LRGS_TV 在消除混合誤差和地形細(xì)節(jié)保留方面都是最佳的，并且對(duì)于局部高程變化躍變現(xiàn)象的約束將減少對(duì)后續(xù)計(jì)算的影響，這一點(diǎn)將在坡度坡向的評(píng)估中展示。

圖3 樣區(qū)1中不同算法誤差消除效果對(duì)比Fig.3 Error elimination effect comparison of different algorithms in sample area 1

從定量評(píng)估的角度來(lái)看，LRGS_TV 在PSNR 和SSIM 指標(biāo)上相較其他方法都具有更好的表現(xiàn)（表1），這與細(xì)節(jié)評(píng)估的結(jié)論相一致。

表1 不同算法在各樣區(qū)的峰值信噪比和結(jié)構(gòu)相似性評(píng)估結(jié)果Tab.1 PSNR and SSIM evaluation results of different algorithms in different sample areas

觀察上述各種方法對(duì)于數(shù)據(jù)的處理結(jié)果可以發(fā)現(xiàn)，LRGS_TV 相較于其他方法能保證X和Y雙方向同步處理，充分考慮到了地形數(shù)據(jù)中細(xì)節(jié)變化，保持了更多的地形細(xì)節(jié)特征。為了進(jìn)一步評(píng)估該方法處理后是否對(duì)數(shù)據(jù)的后續(xù)計(jì)算有影響，對(duì)修復(fù)數(shù)據(jù)前后的坡度坡向提取結(jié)果進(jìn)行評(píng)估。

由于TV方法對(duì)隨機(jī)噪聲消除徹底，處理后的坡度結(jié)果大致分布清晰（圖4（c）），但由于無(wú)法對(duì)具體特征進(jìn)行約束，且存在局部數(shù)據(jù)細(xì)節(jié)丟失與條帶誤差殘留現(xiàn)象，導(dǎo)致該處計(jì)算所得坡度數(shù)據(jù)不可避免地有丟失現(xiàn)象。UTV和LRSID由于無(wú)法有效恢復(fù)高程數(shù)據(jù)，因此其二者坡度結(jié)果不佳（圖4（d）～（e））。經(jīng)過(guò)LRGS_TV 算法對(duì)帶有噪聲的數(shù)據(jù)處理后提取到如圖4（g）所示的坡度結(jié)果，其總體分布的趨勢(shì)和原始坡度一致，保留了相應(yīng)的地形細(xì)節(jié)信息。

圖4 樣區(qū)1中不同算法消除數(shù)據(jù)誤差后計(jì)算所得坡度結(jié)果Fig.4 Slope results calculated by different algorithms after eliminating data errors in sample area 1

相較于原始數(shù)據(jù)的坡度結(jié)果，分別對(duì)比了添加噪聲和各方法去除噪聲后提取的坡度結(jié)果與原始數(shù)據(jù)的坡度作差的情況，坡度作差后的數(shù)據(jù)累計(jì)頻率統(tǒng)計(jì)情況如圖5 所示，其中LRGS_TV 修復(fù)數(shù)據(jù)后提取的坡度相較于原始地形，92.13%的數(shù)據(jù)的坡度差小于8°，在局部特征表現(xiàn)以及定量評(píng)估中都優(yōu)于現(xiàn)有基于TV、UTV 以及低秩分解（LRSID）的方法，能在一定程度上恢復(fù)復(fù)雜干擾下的坡度信息。

圖5 坡度差值累計(jì)頻率Fig.5 Cumulative frequency of slope difference

在樣區(qū)1未加入噪聲前，如圖6（a）所示，不同方向角度的坡向能反映地表的傾斜趨勢(shì)。添加隨機(jī)噪聲和條紋噪聲后，提取的坡向信息混亂，丟失了地形細(xì)節(jié)且出現(xiàn)明顯的條紋狀分布（圖6（b））。使用TV 方法處理后的坡向結(jié)果出現(xiàn)明顯的分布和方向混亂現(xiàn)象（圖6（c）），這表明TV處理數(shù)據(jù)在視覺(jué)效果上有所提升，全局范圍與源數(shù)據(jù)相似，但由于數(shù)據(jù)的過(guò)平滑使局部區(qū)域細(xì)節(jié)丟失，導(dǎo)致在后續(xù)計(jì)算中結(jié)果與實(shí)際不符。UTV 和LRSID 由于無(wú)法完全恢復(fù)高程數(shù)據(jù)，因此其二者坡向結(jié)果不佳（圖6（d）～（e））。經(jīng)過(guò)LRGS_TV 方法對(duì)帶有噪聲的數(shù)據(jù)處理后提取到如圖6（g）所示的坡向結(jié)果，其總體分布和原始結(jié)果一致，不僅能剔除兩種噪聲，還能保留相應(yīng)的地形細(xì)節(jié)信息。

圖6 樣區(qū)1中不同算法消除數(shù)據(jù)誤差后計(jì)算所得坡向結(jié)果Fig.6 Aspect results calculated by different algorithms after eliminating data errors in sample area 1

2.2.2 真實(shí)實(shí)驗(yàn)

樣區(qū)5中存在全局條帶誤差與局部條帶誤差相交叉，同時(shí)混合全局隨機(jī)噪聲（圖7）。從局部地形特征放大的結(jié)果來(lái)看，TV 方法對(duì)于數(shù)據(jù)中的隨機(jī)噪聲消除效果較好，但部分細(xì)節(jié)丟失。從高程值來(lái)看，全局范圍過(guò)平滑現(xiàn)象仍然存在。UTV 和LRSID方法由于無(wú)法約束隨機(jī)方向的條紋誤差，因此對(duì)于多方向混合誤差消除效果欠佳（圖7（c）～（d））。LRGS消除了混合誤差，但在地形中存在一定程度的邊緣扭曲現(xiàn)象，這是因?yàn)閱翁荻确较虻募s束導(dǎo)致了地貌邊緣梯度過(guò)于銳利。針對(duì)這一問(wèn)題，LRGS_TV 對(duì)局部高程值躍變的約束改善了地形邊緣處的效果，且在消除混合誤差和地形細(xì)節(jié)保留方面都是最佳的。

圖7 樣區(qū)5中不同算法誤差消除效果對(duì)比Fig.7 Error elimination effect comparison of different algorithm in sample area 5

該樣區(qū)地形表面存在大范圍平緩地形，進(jìn)行坡度坡向計(jì)算可以發(fā)現(xiàn)（圖8），該區(qū)域坡度變化較小，原始數(shù)據(jù)中存在大范圍的隨機(jī)噪聲導(dǎo)致坡度結(jié)果混亂，經(jīng)過(guò)本文算法處理恢復(fù)了平緩區(qū)域的地表信息，同時(shí)保留了部分起伏的地形細(xì)節(jié)。同時(shí)，該樣區(qū)中的斜條紋明顯，原始數(shù)據(jù)提取的坡度中也存在明顯的條紋狀分布，修復(fù)算法處理后消除了全局的條紋誤差。

圖8 樣區(qū)5中LRGS_TV得到的坡度坡向結(jié)果Fig.8 Results of slope and aspect obtained by LRGS_TV in sample area 5

2.3 算法時(shí)間復(fù)雜度分析

除了對(duì)修復(fù)效果的評(píng)估與分析，算法的時(shí)間復(fù)雜度也是需要討論的一部分。為了比較實(shí)驗(yàn)中各算法的時(shí)間復(fù)雜度，分別在50×50、100×100、250×250、512×512、1000×1000 大小的數(shù)據(jù)矩陣上進(jìn)行算法運(yùn)行，算法的具體運(yùn)行時(shí)間如表2所示。

從表2 中可以發(fā)現(xiàn)，隨著數(shù)據(jù)規(guī)模的增加，基于變分思想的TV 與UTV 算法運(yùn)行時(shí)間均能保持較低水平，而基于優(yōu)化分解角度的LRSID 與本文的LRGS、LRGS_TV 算法其運(yùn)行時(shí)間隨著數(shù)據(jù)規(guī)模的增大顯著增加。

表2 不同數(shù)據(jù)規(guī)模下不同算法修復(fù)數(shù)據(jù)所消耗的時(shí)間 單位：sTab.2 Time consumed by different algorithms to restore data under different data sizes unit：s

從算法原理角度分析，TV 與UTV 定義為梯度幅值的積分，是一個(gè)依靠梯度下降對(duì)數(shù)據(jù)進(jìn)行平滑的各向異性的模型，旨在使得相鄰像素的差值較小，是一種全局平滑思想。該類算法全局平滑，不考慮局部存在的誤差特性，模型簡(jiǎn)單且易于收斂，因此算法的運(yùn)行時(shí)間較少，但相應(yīng)數(shù)據(jù)恢復(fù)效果十分依賴誤差的種類，且從實(shí)驗(yàn)中可看出，數(shù)據(jù)細(xì)節(jié)容易與誤差一起被刪除。

LRSID 將修復(fù)問(wèn)題轉(zhuǎn)換為數(shù)據(jù)分解問(wèn)題，算法主要考慮全局范圍具有正則方向的條紋和大面積隨機(jī)噪聲，將數(shù)據(jù)層和噪聲層分解處理，是一種全局?jǐn)?shù)據(jù)分解、優(yōu)化求解的思想。相較于前兩種全局統(tǒng)一平滑的思想，LRSID 的特征約束模型更為復(fù)雜，考慮到了不同類型的噪聲分量特性，收斂速度隨數(shù)據(jù)規(guī)模的增大降低，因此運(yùn)行時(shí)間隨數(shù)據(jù)規(guī)模的增加而增加。但該算法對(duì)于正則方向條紋特性的約束使得在面對(duì)數(shù)據(jù)中垂直、水平方向條紋誤差時(shí)消除效果優(yōu)異。

本文LRGS 與LRGS_TV 同樣從優(yōu)化角度出發(fā)，考慮了真實(shí)情況下，數(shù)據(jù)中存在局部非正則方向條紋誤差與隨機(jī)噪聲混合的問(wèn)題，將全局?jǐn)?shù)據(jù)分解為多局部數(shù)據(jù)分別進(jìn)行優(yōu)化處理。相較于全局統(tǒng)一優(yōu)化，本文算法針對(duì)復(fù)雜情況下全局?jǐn)?shù)據(jù)中存在多局部非正則方向誤差與隨機(jī)噪聲混合的問(wèn)題處理效果顯著；但同樣的，模型復(fù)雜度進(jìn)一步提升，收斂速度隨數(shù)據(jù)規(guī)模的增加進(jìn)一步下降，會(huì)不可避免地導(dǎo)致算法運(yùn)行時(shí)間成倍增長(zhǎng)。

3 模型性能和應(yīng)用指標(biāo)討論

SRTM 可基本反映地球表面的地形起伏，其值在局部區(qū)域變化平穩(wěn)，在非局部區(qū)域是自相似的，這表示全球高程數(shù)據(jù)應(yīng)該具有低秩特征；但是，衛(wèi)星采集和傳輸過(guò)程中引入的混合誤差會(huì)對(duì)高程數(shù)據(jù)造成影響。以往的研究普遍介紹，高程數(shù)據(jù)的誤差是一個(gè)統(tǒng)一的整體或一個(gè)單頻分量。本文研究發(fā)現(xiàn)全球高程數(shù)據(jù)中存在的混合誤差表現(xiàn)形式多樣，但主要可以分為兩類：系統(tǒng)誤差和隨機(jī)噪聲，具體表現(xiàn)為多方向條帶誤差以及尖峰和斑點(diǎn)問(wèn)題，且系統(tǒng)錯(cuò)誤的低秩組稀疏特征在剔除隨機(jī)噪聲的局部區(qū)域表現(xiàn)更為明顯?，F(xiàn)有研究中的數(shù)據(jù)處理后發(fā)生過(guò)平滑現(xiàn)象，主要是因?yàn)殡p梯度方向處理導(dǎo)致的，因此將處理限制在與條帶相交的方向上的高程數(shù)據(jù)有利于細(xì)節(jié)的保留。

然而，梯度數(shù)據(jù)對(duì)海拔高度突然變化非常敏感。由坡度坡向公式入手分析可知，單方向梯度的約束雖然保證了數(shù)據(jù)圖像更多的細(xì)節(jié)保留，但另一梯度方向的忽略，將會(huì)導(dǎo)致地貌邊緣處的梯度值過(guò)度過(guò)于銳利，這將會(huì)對(duì)后續(xù)計(jì)算造成不利影響。

本文算法在低秩組稀疏約束下，優(yōu)化了梯度方向的處理，使梯度不統(tǒng)一導(dǎo)致的局部高程變化躍變現(xiàn)象減輕，減少了對(duì)后續(xù)計(jì)算的影響；并通過(guò)對(duì)修復(fù)后的高程數(shù)據(jù)進(jìn)行坡度坡向提取，來(lái)驗(yàn)證該算法的修復(fù)效果。結(jié)合坡度坡向的計(jì)算式可知，二者作為評(píng)判指標(biāo)對(duì)于梯度變化敏感，需要同時(shí)考慮水平和垂直兩個(gè)梯度方向的高程變化情況，對(duì)基于高程數(shù)據(jù)后續(xù)應(yīng)用的檢驗(yàn)是嚴(yán)格的。

對(duì)于坡度指標(biāo)，其數(shù)值含義主要是在某一確定方向下的高程差與柵格中心距的反正切角度，這就要求數(shù)據(jù)修復(fù)過(guò)程盡可能小地影響高程數(shù)據(jù)的相對(duì)值，而坡向的計(jì)算主要受中心柵格與其他鄰域柵格的帶權(quán)高程差影響。LRGS_TV 算法在處理過(guò)程中的梯度約束，能夠在一定程度上保證水平和垂直兩個(gè)方向高程同步變化。

4 結(jié)語(yǔ)

本文分析了局部混合誤差的固有特征并構(gòu)建了誤差的稀疏低秩混合錯(cuò)誤模型，使用凸優(yōu)化算法對(duì)模型進(jìn)行求解，提出了LRGS_TV 算法以解決數(shù)據(jù)中的混合誤差問(wèn)題。實(shí)驗(yàn)結(jié)果表明，基于總變分（TV）約束的低秩組稀疏（LRGS_TV）算法在特征細(xì)節(jié)和定量評(píng)估中更好地修復(fù)了原始數(shù)據(jù)，并在基于高程的后續(xù)應(yīng)用計(jì)算中有更好的表現(xiàn)。下一步工作將研究公開(kāi)發(fā)布修復(fù)后全球數(shù)據(jù)；由于復(fù)雜模型與局部特征的引入，算法收斂的時(shí)間復(fù)雜度方面仍需改善，能根據(jù)模型的特點(diǎn)開(kāi)發(fā)出更快收斂的算法也是未來(lái)工作的重點(diǎn)之一；同時(shí)，平坦地區(qū)數(shù)據(jù)的微小變化導(dǎo)致坡向提取結(jié)果存在一定偏差，后續(xù)的研究中也將進(jìn)一步優(yōu)化。