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

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

0 引言

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

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

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

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

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

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

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

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

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

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

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

給定數(shù)據(jù)中的經(jīng)度和緯度坐標（X和Y）以及相應(yīng)的離散點高程，則高程數(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ù)中混合誤差的疊加形式，可將包含隨機噪聲和系統(tǒng)誤差的結(jié)構(gòu)模型描述為：

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1.3 評估指標

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

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

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

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

其中：

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

2 實驗與結(jié)果分析

2.1 航天飛機雷達數(shù)據(jù)與實驗設(shè)計

如圖2所示，實驗使用數(shù)據(jù)集為航天飛機雷達地形測繪任務(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 個區(qū)域文件，單文件包含3601×3601 個采樣點的高度數(shù)據(jù)［18］與對應(yīng)經(jīng)緯度坐標。

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

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

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

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

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

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

2.2.1 模擬實驗

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

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

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

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

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

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

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

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

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

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

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

2.2.2 真實實驗

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

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

該樣區(qū)地形表面存在大范圍平緩地形，進行坡度坡向計算可以發(fā)現(xiàn)（圖8），該區(qū)域坡度變化較小，原始數(shù)據(jù)中存在大范圍的隨機噪聲導(dǎo)致坡度結(jié)果混亂，經(jīng)過本文算法處理恢復(fù)了平緩區(qū)域的地表信息，同時保留了部分起伏的地形細節(jié)。同時，該樣區(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 算法時間復(fù)雜度分析

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)語

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