艾藝紅, 邱 東

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

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

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

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

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

1 問題描述

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

|yi-y(xi)|≤θ，

(1)

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

(2)

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

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

當k=3時，問題(1)即對應了二階導數(shù)的逼近問題.利用下式可以分段構造泛函的極小值：

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)可得相應的線性系統(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，

此時，(3)所示的k=3時問題可以轉化為：

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

(5)

(6)

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

(7)

3 仿真分析

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

3.1 二階導數(shù)近似實驗

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

3.2 不連續(xù)性檢測

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

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

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

4 結 論

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

參考文獻

[1] 鄭志靜, 王璐. 一種基于階次轉換的分數(shù)階微分方程近似解法[J]. 湘潭大學自然科學學報, 2017, 39(4): 10-13.

[2] 張強, 齊興斌. 利用同倫攝動法的四階微分方程數(shù)值解求解方法[J]. 湘潭大學自然科學學報, 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.