艾藝紅, 邱 東

(1.重慶工商大學(xué) 派斯學(xué)院，重慶 401520；2.重慶郵電大學(xué) 理學(xué)院，重慶 400065)

函數(shù)微分方程的求解與運(yùn)算一直以來(lái)都具有非常高的理論研究意義和應(yīng)用價(jià)值[1-2].其中，作為微分方程中重要的研究?jī)?nèi)容，數(shù)值微分是一種利用離散點(diǎn)的函數(shù)值逼近該函數(shù)導(dǎo)數(shù)的方法，它被廣泛應(yīng)用于圖像處理[3]、倒偵測(cè)系統(tǒng)[4]、管道建模與分析[5]等領(lǐng)域中.

現(xiàn)有數(shù)值微分方法主要包括兩類(lèi)：噪聲知識(shí)方法[6-7]，啟發(fā)式方法[8-9].與噪聲知識(shí)方法相比，啟發(fā)式方法并不使用那些可用性較差的先驗(yàn)信息，而是基于正則化方法從問(wèn)題和數(shù)據(jù)中提取這些信息，從而實(shí)現(xiàn)數(shù)值微分的準(zhǔn)確運(yùn)算[10].在現(xiàn)有研究中，蟻群優(yōu)化[11]、L曲線(xiàn)[12]等都是比較成熟的方法，但這些方法在數(shù)值微分中的應(yīng)用卻很少.

在現(xiàn)有方法中，平滑樣條曲線(xiàn)是一種應(yīng)用非常廣泛的數(shù)值微分方法.目前，已有很多學(xué)者研究了正則化參數(shù)選擇問(wèn)題，并給出了相關(guān)的研究成果.文獻(xiàn)[13]對(duì)廣義交叉驗(yàn)證法進(jìn)行了優(yōu)化，通過(guò)引入交替最小化算法提升了參數(shù)選取精度.文獻(xiàn)[14]在傳統(tǒng)廣義交叉驗(yàn)證的基礎(chǔ)上，提出了一種自適應(yīng)參數(shù)獲取方法，提高了算法的精度和魯棒性.文獻(xiàn)[15]提出了一種基于跟蹤平衡差異原理的啟發(fā)式方法，并借助樣條工具箱執(zhí)行相關(guān)仿真程序以驗(yàn)證所提方法在參數(shù)選取方面的優(yōu)越性.

本文針對(duì)樣條曲線(xiàn)中的正則化參數(shù)選擇問(wèn)題，提出一種參數(shù)選取方法：利用基于Hansen的正規(guī)化數(shù)據(jù)包的L曲線(xiàn)對(duì)函數(shù)進(jìn)行平滑處理，采用de Boor算法對(duì)平滑處理后的函數(shù)進(jìn)一步優(yōu)化，進(jìn)而借助正則化方法選擇最優(yōu)參數(shù).基于Matlab仿真軟件中的de Boor樣條工具箱和正則化工具箱，驗(yàn)證了所提啟發(fā)式算法在平滑曲線(xiàn)參數(shù)正則化選擇中的有效性.

1 問(wèn)題描述

在給出解決方案之前，需要先對(duì)研究的問(wèn)題進(jìn)行簡(jiǎn)要敘述.假設(shè)y=y(x) 為[0,1]之間的函數(shù)，n為[0,1]之間的自然數(shù)，x0

|yi-y(xi)|≤θ，

(1)

式中：i=1, …,n.所研究的問(wèn)題即為基于測(cè)得的樣本yi來(lái)構(gòu)造函數(shù)的近似值及其特定階次導(dǎo)數(shù).該問(wèn)題可描述為：

(2)

式中：γ為正則化參數(shù)，當(dāng)γ=0時(shí)，極小值為插值數(shù)據(jù)的曲線(xiàn)，當(dāng)γ→∞時(shí)，極小值趨近于數(shù)據(jù)的最適合直線(xiàn)；g(x)∈Hk(0,1)，g(0)=y0，g(1)=yn；Hk(0,1)的具體表達(dá)式可參見(jiàn)文獻(xiàn)[15].

2 新型正則化參數(shù)選擇方法

當(dāng)k=3時(shí)，問(wèn)題(1)即對(duì)應(yīng)了二階導(dǎo)數(shù)的逼近問(wèn)題.利用下式可以分段構(gòu)造泛函的極小值：

g(x)=αi+βi(x-xi)+χi(x-xi)2+δi(x-xi)3+εi(x-xi)4+φi(x-xi)5

.

(3)

式中：xi=i/n.由式(3)可得相應(yīng)的線(xiàn)性系統(tǒng)矩陣形式為：

(4)

式中：d=1/n，α=(α1,α2,…, αn-1)T,χ= (χ1,χ2,…, χn-1)T，ε=(ε1,ε2,…,εn-1)T,y=(y1,y2,…,yn-1)T,z=(z1,z2,…,zn-1)T，

此時(shí)，(3)所示的k=3時(shí)問(wèn)題可以轉(zhuǎn)化為：

min(qU1(gq)+(1-q)V1(gq)).

(5)

(6)

q=1/(1+d4R(W1)/R(W2)),qR(·)=6(1-q)ATB2A.

(7)

3 仿真分析

為了驗(yàn)證所提方法在數(shù)值微分參數(shù)選擇中的優(yōu)勢(shì)，在Matlab中進(jìn)行了仿真分析.仿真中，網(wǎng)格大小n=10，σ=δ2，U(g)=δ4.對(duì)于L曲線(xiàn)方法，使用文獻(xiàn)[15]中的正則化工具箱，并對(duì)L曲線(xiàn)、廣義交叉驗(yàn)證、de Boor等方法在選擇正則化參數(shù)時(shí)的性能進(jìn)行了對(duì)比分析.

3.1 二階導(dǎo)數(shù)近似實(shí)驗(yàn)

實(shí)驗(yàn)中，選取d=1/n=0.1.在Matlab中的函數(shù)rand=(size(x))中增加噪聲δ=3×10-2，得到rand=δ·(size(x))，用來(lái)模擬所提方法在提取數(shù)據(jù)時(shí)可以達(dá)到的精確程度，其中x為網(wǎng)格點(diǎn)向量.對(duì)比了de Boor方法和噪聲采集方法的性能.實(shí)驗(yàn)結(jié)果如圖1所示.由實(shí)驗(yàn)結(jié)果可知，與噪聲采集方法相比，所提方法雖然沒(méi)有利用與噪聲相關(guān)的任何信息，但是卻可以實(shí)現(xiàn)更好的跟蹤性能，導(dǎo)數(shù)逼近誤差幾乎為零，而噪聲采集方法則始終存在著較大的逼近誤差.

3.2 不連續(xù)性檢測(cè)

為了驗(yàn)證所提方法在分段函數(shù)中參數(shù)選擇的有效性，考慮以下函數(shù)：

對(duì)兩種情況下的分段函數(shù)進(jìn)行了仿真分析，包括：a=0.1，b=0.5，c=0.8和a=0.2，b=0.65，c=0.9.對(duì)分段函數(shù)來(lái)說(shuō)，由于其具備明顯的跳躍性，因此不僅具有一定的典型性，還可以充分說(shuō)明所提方法在數(shù)值微分參數(shù)選擇中的高精確性和強(qiáng)適應(yīng)性.在實(shí)驗(yàn)中，通過(guò)近似一階導(dǎo)數(shù)來(lái)檢測(cè)函數(shù)的不連續(xù)點(diǎn)，噪聲為δ=5×10-3的不連續(xù)點(diǎn)跳躍.

由圖2可知，由于L曲線(xiàn)選擇的q幾乎等于1，其結(jié)果是漸近正態(tài)計(jì)算解，因此該方法并不利于不連續(xù)點(diǎn)的檢測(cè)，在a=0.2，b=0.65，c=0.9的函數(shù)中，L曲線(xiàn)甚至?xí)霈F(xiàn)選擇失敗.相反，de Boor方法獲得了非常好的實(shí)驗(yàn)結(jié)果.

4 結(jié) 論

針對(duì)離散函數(shù)導(dǎo)數(shù)逼近，采用光滑樣條函數(shù)以及不連續(xù)性預(yù)測(cè)對(duì)參數(shù)進(jìn)行正則化選取，在不依賴(lài)噪聲測(cè)量的同時(shí)實(shí)現(xiàn)了參數(shù)的準(zhǔn)確選擇.通過(guò)將L曲線(xiàn)與de Boor方法進(jìn)行有機(jī)整合，可以提高算法的精確度.在二階導(dǎo)數(shù)逼近問(wèn)題中進(jìn)行了數(shù)值仿真實(shí)驗(yàn),實(shí)驗(yàn)結(jié)果表明所提方法具有非常好的選取結(jié)果，甚至優(yōu)于基于噪聲的參數(shù)選取方法.

參考文獻(xiàn)

[1] 鄭志靜, 王璐. 一種基于階次轉(zhuǎn)換的分?jǐn)?shù)階微分方程近似解法[J]. 湘潭大學(xué)自然科學(xué)學(xué)報(bào), 2017, 39(4): 10-13.

[2] 張強(qiáng), 齊興斌. 利用同倫攝動(dòng)法的四階微分方程數(shù)值解求解方法[J]. 湘潭大學(xué)自然科學(xué)學(xué)報(bào), 2017, 39(2): 5-8.

[3] ELLIOTT C M,SMITHEMAN S A. Numerical analysis of the TV regularization and H-1 fidelity model for decomposing an image into cartoon plus texture[J]. Communications on Pure & Applied Analysis, 2017, 6(4): 917-936.

[4] WAGNER J,MAZUREK P,MORAWSKI R Z, et al. Regularized numerical differentiation of depth-sensor data in a fall detection system[C]//IEEE International Conference on Computational Intelligence & Virtual Environments for Measurement Systems & Applications, 2017: 234-236.

[5] MAYSTRENKO A V,SVETLAKOV A A,GANDZHA T V, et al. Application of numerical signal differentiation methods to determine statio-narity of a process[J]. Petroleum & Coal, 2017, 59(3): 311-318.

[6] SZOSTOK T. Functional inequalities involving numerical differentiation formulas of order two[C]//Bulletin of the Malaysian Mathematical Sciences Society, 2017: 1-14.

[7] GAO W,ZHANG R. Multiquadric trigonometric spline quasi-interpolation for numerical differentiation of noisy data: a stochastic perspective[J]. Numerical Algorithms, 2017, 12(3): 1-17.

[8] ALBANI V,CEZARO A D,ZUBELLI J P. On the choice of the Tikhonov regularization parameter and the discretization level: a discrepancy-based strategy[J]. Inverse Problems & Imaging, 2017, 10(1): 1-25.

[9] ALJAMAL M F,ALOMARI A K,GOCKENBACH M S. Smoothing via elliptic operators with application to edge detection[J]. Inverse Problems in Science & Engineering, 2017,18(4): 1-20.

[10] DAVYDOV O,SCHABACK R. Error bounds for kernel-based numerical differentiation[J]. Numerische Mathematik, 2016, 132(2): 243-269.

[11] LAI M,TONG X. A metaheuristic method for vehicle routing problem based on improved ant colony optimization and Tabu search[J]. Journal of Industrial & Management Optimization, 2017, 8(2): 469-484.

[13] ZHANG X,JAVIDI B,NG M K. Automatic regularization parameter selection by generalized cross-validation for total variational Poisson noise removal[J]. Appl Opt, 2017, 56(9): 35-47.

[14] HUANG G, LI J, LUO C, et al. Regularization parameter adaptive acquisition based on improved GCV method and its application in pre-stack AVO inversion[C]//Seg Technical Program Expanded, 2017: 748-752.

[15] FANG H,CHANG Y,ZHOU G, et al. Iteratively reweighted blind deconvolution with adaptive regularization parameter estimation[C]//IEEE Access, 2017,99: 1-1.