同倫攝動稀疏正則化方法及其在

1、同倫攝動稀疏正則化方法及其在初始地形地貌重構(gòu)中的應(yīng)用 摘要：稀疏正則化方法在參數(shù)重構(gòu)中起到了越來越重要的作用。與傳統(tǒng)的正則化方法相比，稀疏正則化方法能較好地重構(gòu)稀疏變量。由于稀疏正則化的不可微性，需要對已有的經(jīng)典算法進行改進。本文構(gòu)建同倫攝動稀疏正則化方法克服標準稀疏正則化的不可微性，并應(yīng)用該方法應(yīng)用到項目問題中，能夠有效地重構(gòu)山體初始表面。數(shù)值實驗表明，所提出的方法是收斂和穩(wěn)定的。 關(guān)鍵詞：稀疏正則化；同倫攝動稀疏正則化方法；參數(shù)重構(gòu)；重構(gòu)山體初始表面1. 引言 參數(shù)重構(gòu)在許多應(yīng)用中起到重要作用，如波動率和政策參數(shù)重構(gòu)1, 2，在其他工程實踐領(lǐng)域也有重要的應(yīng)用，例如：圖像去噪3, 4，地震信

2、號重構(gòu)5, 6以及心電信號重構(gòu) 7, 8。一般情況下，參數(shù)重構(gòu)問題是不適定的，也就是說，即使測量數(shù)據(jù)的噪聲水平很小，也可能導致重構(gòu)結(jié)果嚴重偏離真實參數(shù)9。正則化方法主要就是為了克服參數(shù)重構(gòu)的不適定性，通過選取合適的正則化方法能夠抑制觀測數(shù)據(jù)中誤差對于參數(shù)重構(gòu)的不良影響，獲得較為準確的重構(gòu)值。使用最為普遍的正則化方法是吉洪諾夫正規(guī)化方法，它的目標函數(shù)是由擬合項和懲罰項構(gòu)成。吉洪諾夫正則化方法已被用于許多參數(shù)辨識問題中，很多學者對其數(shù)值計算方法進行了研究，如Landweber方法10，高斯牛頓法11，Newton-Kaczmarz方法12和多尺度平滑方法13，這些方法能夠有效地重構(gòu)光滑參數(shù)。隨著經(jīng)

3、濟和金融理論的發(fā)展，波動率和經(jīng)濟政策參數(shù)的重構(gòu)已經(jīng)廣泛應(yīng)用于許多實際問題中。在實際應(yīng)用中，很多需要重構(gòu)的參數(shù)都是稀疏的，即參數(shù)的非零元素的數(shù)目非常有限，遠遠小于零元素的數(shù)目。即使重構(gòu)參數(shù)稀疏程度不夠，也可以利用小波和曲波變換使參數(shù)稀疏化。本文先對該方法在經(jīng)濟和金融領(lǐng)域中的應(yīng)用進行研究，因為這兩個領(lǐng)域中的參數(shù)通?？梢苑殖梢阎糠趾臀粗糠?。已知部分通常與已有的經(jīng)濟和金融政策相對應(yīng)，而未知部分通過參數(shù)重構(gòu)，再將結(jié)果進行經(jīng)濟學和金融學解釋，能夠為政策制定部門提供行之有效的對策建議。通過這兩類問題的研究表明該方法能夠應(yīng)用到參數(shù)分解成已知部分和未知部分的問題中。然后再將該方法應(yīng)用到較為復(fù)雜的山體表面重構(gòu)

4、問題中。傳統(tǒng)的吉洪諾夫正則化方法對于重構(gòu)稀疏參數(shù)效果很不好，而稀疏正則化方法卻能較好地重構(gòu)稀疏參數(shù)，但是稀疏正則化是不可微的，因此需要采用一些技巧來克服這一困難，典型的方法是Bregman迭代14, 15-17。本文構(gòu)造同倫攝動稀疏正則化方法，達到提高算法精度和提高計算效率的目的。數(shù)值實驗表明同倫攝動稀疏正則化方法是收斂和穩(wěn)定的。2. 稀疏正則化方法在實際應(yīng)用中，對于光滑參數(shù)重構(gòu)，吉洪諾夫正則化方法具有良好的收斂性和穩(wěn)定性。但是，在稀疏參數(shù)重構(gòu)時，吉洪諾夫正則化方法的重構(gòu)效果很差，很難滿足工程實踐的要求。因此，需要采用稀疏正則化方法進行參數(shù)重構(gòu)。2.1 吉洪諾夫正則化方法參數(shù)重構(gòu)問題可以歸結(jié)為

5、非線性算子方程(1)的形式：(1)其中， 分別代表非線性算子、需要重構(gòu)的參數(shù)和觀測數(shù)據(jù)。在實際問題中，往往還需要考慮觀測數(shù)據(jù)和理想數(shù)據(jù)之間的噪音水平，即：(2)其中，分別代表真實的觀測數(shù)據(jù)和噪音水平。 參數(shù)重構(gòu)的難點在于不適定性，很小的噪音水平 也會使得重構(gòu)結(jié)果嚴重背離真實的物理參數(shù)，從而造成結(jié)果的無意義。解決這一難點最重要的方法就是吉洪諾夫正則化方法，與之對應(yīng)的吉洪諾夫正則化目標函數(shù)定義為：(3)其中， 是在范數(shù)意義下的數(shù)據(jù)擬合項, 是起到穩(wěn)定作用的懲罰項, 參數(shù)為正則化參數(shù), 該參數(shù)主要是起到平衡數(shù)據(jù)擬合項和懲罰項的作用。吉洪諾夫正則化方法主要是求解下面的優(yōu)化問題.(4)最小值滿足如下的歐

6、拉方程：(5)其中， 是F-導數(shù)。在求解方程(5)的時候，采用的最普遍的方法是Landweber方法，該方法可以寫成下面的表達式：,(6)其中，表示迭代次數(shù)。方程(6)是著名的Landweber迭代格式，該數(shù)值格式的顯著優(yōu)點是穩(wěn)定性特別好，但是，收斂速度很慢，不適合應(yīng)用到大型實際問題中。另一個重要的方法就是高斯-牛頓方法，該方法收斂速度較快，但是，不如Landweber方法穩(wěn)定。本文主要是在Landweber方法基礎(chǔ)上研究同倫攝動稀疏正則化方法，整個過程可以推廣到高斯-牛頓方法上。2.2 稀疏正則化方法為了能夠有效地重構(gòu)稀疏變量，將標準的吉洪諾夫正則化方法進行改進，使得吉洪諾夫正則化目標函數(shù)轉(zhuǎn)

7、換為如下形式:,(7)其中，表示向量的非零元素的個數(shù)。近些年，由于稀疏正則化方法能夠有效地重構(gòu)稀疏變量，使得反問題領(lǐng)域?qū)W者越來越多地重視該方法在實際問題中的應(yīng)用。該方法的難點在于，泛函(7)的懲罰項是不可微的，而且該問題還是NP不可解問題。為了解決NP不可解難點，采用下面的形式進行替代：,(8)其中，范數(shù)表示. 在泛函(8)中, 利用范數(shù)代替了原有的范數(shù)，這樣的改進使得計算效率得到了顯著的提高。盡管進行了上述的改進，但是(8)的懲罰項依然是不可微的，因此，還要在此基礎(chǔ)上進一步改進，將(8)中的懲罰項改為下面的帶有阻尼系數(shù)的范數(shù)：,(9)其中，利用代替范數(shù), 輔助參數(shù)是一個正實數(shù)。與吉洪諾夫正則

8、化方法相仿，將Landweber方法應(yīng)用到(9)中，得到下面的迭代格式：.(10)3. 同倫攝動稀疏正則化方法 為了提高數(shù)值迭代格式(10)的收斂速度，將同倫攝動反演算法與稀疏正則化方法有機地融合到一起，達到提高重構(gòu)效果和收斂速度的目的。構(gòu)造如下的同倫映射 ,(11)其中，代表嵌入?yún)?shù), 迭代算法的初始值。因此，有下面的關(guān)系式：(12)將寫成關(guān)于嵌入?yún)?shù)的冪級數(shù)形式：(13)目標函數(shù)(9)的解可以表示為：.(14)將方程(11)中的在點寫成泰勒級數(shù)的形式：(15)可以得到下面的關(guān)系式：.(16)和.(17)根據(jù)(17)，對于含有噪音的觀測數(shù)據(jù)，采用前兩項構(gòu)造稀疏變量:,(18)其中，表示迭代次

9、數(shù)。 采用第一項構(gòu)造稀疏變量，則有下面的迭代格式：,(19)迭代格式(19)正是前面介紹的Landweber方法，而迭代格式(18)可以看做是在其上面進行的修正，將更多的信息應(yīng)用到迭代算法中，可以提高計算效率。4 數(shù)值算例4.1 重構(gòu)隱含波動率在金融領(lǐng)域中，基于布萊克一斯科爾斯期權(quán)定價模型(Black-Scholes (B-S) option pricing model)重構(gòu)隱含波動率一直是重要的研究方向。利用有限差分方法求解該模型的正演問題，在重構(gòu)隱含波動率的時候采用同倫攝動稀疏正則化方法，達到提高重構(gòu)精度和計算效率的目的。布萊克一斯科爾斯期權(quán)定價模型表明金融衍生品的價格變化滿足標準幾何布朗

10、運動。該模型的邊界條件隨著金融衍生品的變化而變化。當邊界條件給定的時候，金融衍生品的價格可以通過求解布萊克一斯科爾斯期權(quán)定價模型獲得。當金融衍生品選定為期權(quán)時，本文定義正演問題為求解期權(quán)的價格，反演則是重構(gòu)隱含波動率。根據(jù)期權(quán)的行權(quán)方式的不同，可以分為歐式期權(quán)和美式期權(quán)，本節(jié)以歐式期權(quán)為研究對象。 在時間上的歐式看漲期權(quán), 設(shè) 是歐式看漲期權(quán)在區(qū)域上的價格, 滿足下面的方程：,(20) 其中，: 股票價格, : 執(zhí)行價格, : 利率, : 紅利, : 到期日, : 時間, : 隱含波動率。.反問題可以表述為: 從觀測數(shù)據(jù), 重構(gòu)隱含波動率.定義非線性向量值函數(shù), 即：. 為了測試本文提出的方法

11、和標準的吉洪諾夫正則化方法之間的差別，設(shè)置如下參數(shù)：到期日, 股票價格，利率, 執(zhí)行價格. 為了測試算法的抑制噪音的能力，在觀測數(shù)據(jù)上添加1%的高斯噪音。真實的隱含波動率如下：(21)將真實的波動率分解為兩部分：,(22)其中，分別表示常值隱含波動率和稀疏隱含波動率。設(shè)是已知的, 需要重構(gòu)部分，該部分為：(23) 目標函數(shù)(9) 可以寫成：,(24)采用同倫攝動稀疏正則化方法重構(gòu)結(jié)果為： (25)采用標準吉洪諾夫正則化方法重構(gòu)結(jié)果為： (26)從上面的重構(gòu)結(jié)果可以看出，同倫攝動稀疏正則化方法能夠提高標準吉洪諾夫正則化方法對于重構(gòu)稀疏隱含波動率的精度。4.2 托達羅模型的政策參數(shù)重構(gòu)托達羅模型是

12、著名的描述人口流動的經(jīng)濟學模型，在勞動經(jīng)濟學領(lǐng)域具有重要的地位。托達羅模型將人口流動的數(shù)量和不同地區(qū)的收入差距聯(lián)系到一起：(27)其中，分別表示人口流動的數(shù)量和不同地區(qū)之間的收入差距。函數(shù)是單調(diào)遞增的，即. 隨著經(jīng)濟的發(fā)展，標準的托達羅模型(27)已經(jīng)不能準確地描述人口流動和收入差距之間的關(guān)系，還需要考慮政策參數(shù)的作用，將方程(27)改為下面的形式： (28)其中，表示政府部門的政策參數(shù)。 設(shè)有個農(nóng)村地區(qū)和個城市地區(qū)，第個農(nóng)村區(qū)域和第個城市區(qū)域之間的收入差距表示為(), 第個和 第個城市區(qū)域之間收入差距表示為(). 人口轉(zhuǎn)移到第個城市區(qū)域的數(shù)量用表示。 對于第城市區(qū)域, 政策參數(shù)可以分解為兩部

13、分, 其中分別表示政策參數(shù)對于農(nóng)村和城市的作用效果。為了簡化討論，假設(shè)從城市區(qū)域轉(zhuǎn)移到農(nóng)村區(qū)域的人口數(shù)量為零。在實際應(yīng)用中，函數(shù)有多種表達方式，由于本文主要是測試同倫攝動稀疏正則化方法的有效性，因此，函數(shù)選取為線性函數(shù)的形式。托達羅模型可以寫成下面的形式：(29)其中，、和 已知，政策參數(shù)向量是未知的。 令 ,(30)和,.(31) 方程 (30) 可以寫成下面的形式：.(32)值得注意的是方程(32)是欠定的，因此，需要正則化方法求解該方程。將政策參數(shù)進行分解 , 其中分別表示以前的政策參數(shù)(已知)和稀疏政策參數(shù)(未知)。目標函數(shù)改寫為：.(33)在數(shù)值算例中， 設(shè), 的每個元素都是. 為了

14、測試算法的抑制噪音的能力，在觀測數(shù)據(jù)上添加1%的高斯噪音。重構(gòu)5個稀疏政策參數(shù)，每一個 有且只有一個非零元素，大小為0.5. 同倫攝動稀疏正則化方法能夠精確地將非零元素的位置判斷出來，但是標準的吉洪諾夫正則化方法就不具備這樣的能力。表1給出了兩種不同正則化方法重構(gòu)的均方差。表1：兩種不同方法的均方差Table 1: Mean Square Errors非零元素的位置同倫攝動稀疏正則化方法吉洪諾夫正則化方法25.25%10.07%45.19%12.96%64.69%15.78%84.92%14.06%104.71%9.93%4.3重構(gòu)初始地形地貌 地形地貌演變模型將引起地球表面變化的內(nèi)因和外因有

15、機的結(jié)合到一起，影響地形地貌演變的內(nèi)因是地球內(nèi)部板塊速度場, 而外因則非常多. 主要有以下三個因素：(I) 地球表面擴散過程(diffusion of hillslope topography): 這里所說的擴散過程是一個包含很多復(fù)雜因素的統(tǒng)稱，例如：風化作用(weathering), 坡面沖刷(slope wash), 地面水流(overland flow), 土壤滑動(soil creep)和基巖滑坡引起的質(zhì)量損失(mass wasting by bedrock-involved landsliding); (II) 河流沖刷對于基床的切割(bedrock incision)；(III)

16、冰川的融化影響(melting effect).建立地形地貌演變數(shù)學模型是屬于地形地貌學的范疇，描述上面因素的數(shù)學公式一般都是高度非線性的，將這些具體模型的模型進行統(tǒng)一考慮，可以得到如下模型： (34)模型正演定義為： 給 定 地 球 表 面 的 速 度 場和 初 始 地 形 地 貌 分 布 函數(shù)求解地形地貌隨時間變化關(guān)系，寫成數(shù)學表達式為： (36)模型反演的主要目的是利用現(xiàn)在的地形地貌，重構(gòu)初始地形地貌，可以表示為： (37)值得注意的是我們這里假設(shè)山體的隨時間運動的速度場是已知的，將研究的重點放到了地形地貌運動模型上，并沒有放到地球內(nèi)部熱場分布模型上，該部分的研究已經(jīng)在上一年度完場。方程

17、 (37) 可以寫成下面的線性形式：.(38)值得注意的是方程(38)是不適定的，因此，需要正則化方法求解該方程。將初始地形地貌進行分解 , 其中分別表示通過其他探測手段已知的地形地貌成分(尺度較大)和細節(jié)地形地貌成分(尺度較小)。目標函數(shù)改寫為：.(39)在數(shù)值算例中，在小波基下表示為稀疏系數(shù)，重構(gòu)這些稀疏系數(shù)。 為了測試算法的抑制噪音的能力，在觀測數(shù)據(jù)上添加1%的高斯噪音。同倫攝動稀疏正則化方法能夠精確地將非零元素的位置判斷出來，但是標準的吉洪諾夫正則化方法就不具備這樣的能力。表2給出了兩種不同正則化方法重構(gòu)的均方差。從表2中還可以看出新方法更加穩(wěn)定。表2：兩種不同方法的均方差Table2

18、: Mean Square Errors實驗次數(shù)同倫攝動稀疏正則化方法吉洪諾夫正則化方法12.21%8.07%22.16%8.16%32.23%7.98%42.39%9.01%52.31%8.58%5 結(jié)論本文構(gòu)造了同倫攝動稀疏正則化方法，提高了標準吉洪諾夫正則化方法重構(gòu)稀疏變量的精度，同時也提高了計算效率，并將該方法應(yīng)用到重構(gòu)初始地形地貌問題中。對于實際情況而言，絕大部分參數(shù)變量都可以分解為已知部分和稀疏部分，因此只需要重構(gòu)對應(yīng)的稀疏部分即可，說明該方法具有較好的應(yīng)用前景。同時，數(shù)值算例表明本文構(gòu)造的算法具有較高的精度和收斂速度。參考文獻1 Egger H, Engl H W. Tikhon

19、ov regularization applied to the inverse problem of option pricing: convergence analysis and rateJ. Inverse Problems, 2005, 21: 10271045.2 Trong D D, Thanh D, Lan N N, Uyen P H. Calibration of the purely T-dependent Black-Scholes implied volatilityJ. Applicable Analysis, 2014, 93: 859-874.3 Wu C D.

20、A new denoising method based on multistage median filter and nonsubsampled contourlet transformJ. Journal of Computational Information Systems, 2012, 8(12): 4881-4888.4 Zhao Q, Ye B, Cao J. Image denoising based on improved non-local means and nonsubsampled contourlet transform wiener filteringJ. Jo

21、urnal of Computational Information Systems, 2010, 6 (2): 601-609.5 Baig A M, M. Campillo M, Brenguier F. Denoising seismic noise cross correlationsJ. Journal of Geophysical Research, 2009, 114: doi: 10.1029/2008JB006085. 6 Hennenfent G, Herrmann F J. Simply denoise: wavefield reconstruction via jitt

22、ered undersamplingJ. Geophysics, 2008, 73 (3): V19-V28.7 Istepanian R H, Petrocian A A. Optimal zonal wavelet-based ECG data compression for a mobile telecardiology systemJ. IEEE Transactions on Information Technology in Biomedicine, 2000, 4 (3): 200-211.8 Sameni R, Shamsollahi M B, Jutten C, Cliffo

23、rd G.D. A nonlinear Bayesian filtering framework for ECG denoisingJ. IEEE Transactions on Biomedical Engineering, 2007, 54 (12): 2172-2185.9 Engl H W, Hanke M, Neubauer A. Regularization of Inverse ProblemsM. Kluwer: Dordrecht, 1996.10 Hanke M, Neubauer A, Scherzer O. A convergence analysis of the Landweber iteration for nonlinear ill-posed problemsJ. Numer. Math., 1995, 72: 21-37.11 Blaschke B, Neubauer A, Scherzer O. On convergence rates for the iteratively regularized Gauss-Newton methodJ. IMA J. Numer. Anal., 1997, 17: 421-436.12 Burger M, Kaltenbacher B. R