王 月，蔣 泉，2*

(1.南通大學(xué) 理學(xué)院，江蘇 南通 226019；2.南通大學(xué) 交通與土木工程學(xué)院，江蘇 南通 226019)

在工程中很多問題均歸結(jié)于偏微分方程的求解[1-3]，其正的問題研究在理論和算法上已經(jīng)較為成熟，但在反問題方面還需進(jìn)一步開展深入的研究，例如：在已知部分區(qū)域邊界條件的情況下，由觀測(cè)點(diǎn)的函數(shù)值反演出Poisson 方程的未知的邊界處函數(shù)值；或在函數(shù)中存在部分參數(shù)未知情況下，由觀測(cè)點(diǎn)測(cè)量值及邊界條件反演出函數(shù)的參數(shù)。Yang等[4]人利用簡化Tikhonov 正則化方法探討了半帶型區(qū)域中含有一個(gè)變量的二維Poisson 方程未知源識(shí)別反問題，得到問題的正則近似解，并給出了正則解和精確解之間的H?lder 型誤差估計(jì)；Hamad 等[5]人研究了Poisson 方程源函數(shù)反演迭代算法，其算法能夠根據(jù)邊界處的測(cè)量值反演未知源函數(shù)的近似估計(jì)，并得到了良好的結(jié)果；Koulouri 等[6]人通過貝葉斯近似方法研究了Poisson 方程的逆源問題，并將該方法用于腦電圖成像，得到相對(duì)準(zhǔn)確的結(jié)果；Rond 等[7]人將高斯去噪算法用于Poisson 噪聲反問題中，解決了Poisson 定向問題的信噪比范圍問題。

基本解方法(method of fundamental solutions，MFS)的實(shí)質(zhì)是基本解函數(shù)的疊加[8-9]，在具體計(jì)算中不需要積分，可大大減少計(jì)算量，節(jié)約計(jì)算時(shí)間；同時(shí)具備應(yīng)用性較強(qiáng)、易推廣至高維領(lǐng)域的優(yōu)點(diǎn)。Sun等[10]人利用MFS 對(duì)含邊界噪聲的熱傳導(dǎo)反問題進(jìn)行了研究，其數(shù)值算例證明了該方法的穩(wěn)定性和精確性；Karageorghis 等[11]人將MFS 引入到熱彈性反問題，對(duì)MFS 形成的高維病態(tài)方程的正則化參數(shù)選取進(jìn)行了分析，其數(shù)值結(jié)果具有較好的穩(wěn)定性；Johansson 等[12]人對(duì)非定常熱傳導(dǎo)反問題進(jìn)行了研究，進(jìn)一步說明了MFS 具有高度的穩(wěn)定性、耗費(fèi)計(jì)算資源少的優(yōu)點(diǎn)。

Kalman 濾波算法是處理反問題最為直接有效的方法之一，可以在系統(tǒng)狀態(tài)空間表示的基礎(chǔ)上，有效地處理測(cè)量數(shù)據(jù)噪聲帶來的不確定性，最終求得系統(tǒng)狀態(tài)的最優(yōu)估計(jì)，具有結(jié)構(gòu)簡單和易于實(shí)現(xiàn)的優(yōu)點(diǎn)[13]。目前，Kalman 濾波器在工程系統(tǒng)中得到了廣泛的應(yīng)用，包括導(dǎo)航定位系統(tǒng)[14-16]、無人機(jī)的實(shí)時(shí)追蹤[17]、雷達(dá)的目標(biāo)追蹤[18-19]等方面。同時(shí)在大氣觀測(cè)、儀器檢測(cè)[20]等方面也得到了廣泛的應(yīng)用，如：Rincón 等[21]人結(jié)合Kalman 濾波器對(duì)輸出進(jìn)行偏差校正，給出了天氣模型預(yù)測(cè)(WRF-ARW)的后處理分析，提高了整體預(yù)測(cè)值的正確率；Baptista 等[22]人研究了Kalman 濾波對(duì)剩余使用壽命估計(jì)值的適用性，結(jié)果表明該技術(shù)具有較好的精度和收斂性，大幅度提高了設(shè)備壽命回歸模型的精度；Kang 等[23]人給出了基于慣性測(cè)量單元估計(jì)姿態(tài)和陀螺偏差的無跡Kalman 濾波算法，并通過實(shí)驗(yàn)和仿真進(jìn)一步驗(yàn)證了該算法的有效性。

本文應(yīng)用Kalman 濾波技術(shù)，對(duì)Poisson 方程數(shù)值求解問題中的不完備邊界及未知參數(shù)反問題進(jìn)行分析和研究。利用基本解方法對(duì)調(diào)和函數(shù)的邊界值問題進(jìn)行計(jì)算，對(duì)獲得的高度不適定高維線性方程利用截?cái)嗥娈愔捣纸夥ㄟM(jìn)行數(shù)值求解；進(jìn)一步通過Kalman 濾波技術(shù)，利用有限個(gè)觀測(cè)點(diǎn)的函數(shù)測(cè)量值最終反演出未知邊界處函數(shù)值和未知參數(shù)值。

1 基本方程

1.1 基本解方法

二維平面內(nèi)區(qū)域Ω 上的泊松方程為

式中：φ 為區(qū)域中未知函數(shù)值；Δ2為Laplace 算子；q為已知函數(shù)；λ 為已知參數(shù)；Ω1為Dirichlet 邊界；Ω2為Neumann 邊界。

根據(jù)基本解方法原理，在二維問題中，區(qū)域Ω內(nèi)的任何一點(diǎn)P，式(1a)中函數(shù)φ(P)可寫成基本解函數(shù)的疊加[19-20]，

式中：Ci為待定系數(shù)；Q 為區(qū)域Ω 外的源點(diǎn)；n 為源點(diǎn)總數(shù)目；φ*(P)為式(1)的特解；G1(P，Qi)為齊次調(diào)和函數(shù)φ 解的基本函數(shù)，具體形式為

其中ri為點(diǎn)P 和Qi之間的距離。在直角坐標(biāo)系中，如P 和Qi的坐標(biāo)為(x，y)和(xi，yi)，則有如下關(guān)系：

對(duì)于特解φ*(P)，當(dāng)q 為常數(shù)時(shí)，可寫成如下形式[24-25]：

當(dāng)(xq，yq)處存在強(qiáng)度為q0的源點(diǎn)時(shí)，特解可為

根據(jù)邊界條件(1b)和(1c)，式(2)中待定常數(shù)Ci可由如下的線性方程組求出：

式中：n 為邊界點(diǎn)數(shù)；n1為Dirichlet 邊界點(diǎn)數(shù)；n -n1為Neumann 邊界點(diǎn)數(shù)。由上述方程組求出待定系數(shù)后，任意一點(diǎn)的函數(shù)值就可以按式(2)得到。

1.2 Kalman 濾波基本方程

根據(jù)Kalman 濾波理論，t 狀態(tài)下的系統(tǒng)可由t -1 狀態(tài)進(jìn)行描述，其狀態(tài)方程和觀測(cè)方程可寫成[13]

式中：Xt為狀態(tài)向量；c 為系統(tǒng)矩陣；Wt為系統(tǒng)噪聲；Zt為觀測(cè)向量；Ht為觀測(cè)矩陣，且當(dāng)觀測(cè)點(diǎn)的數(shù)目為m 時(shí)，Ht為m 維列向量；Vt為測(cè)量噪聲。

Kalman 濾波計(jì)算過程如下所示：

1)預(yù)測(cè)系統(tǒng)狀態(tài)

2)預(yù)測(cè)系統(tǒng)均方差

3)修正卡爾曼混合系數(shù)

4)修正系統(tǒng)狀態(tài)

5)修正系統(tǒng)均方差

式中：I 為單位矩陣；a 為系統(tǒng)參數(shù)，對(duì)于靜態(tài)問題一般為單位矩陣；為先驗(yàn)估計(jì)；為后驗(yàn)估計(jì)；和Pt為預(yù)測(cè)系統(tǒng)均方差和修正系統(tǒng)均方差；Kt為卡爾曼增益矩陣；Q 為過程噪聲方差；R 為測(cè)量噪聲方差。

對(duì)于非線性測(cè)量系統(tǒng)，觀測(cè)方程(7b)通常用下式來描述：

當(dāng)h(Xt)足夠平滑時(shí)，為了線性化，可在先驗(yàn)估計(jì)Xt=處泰勒展開成[26]

將觀測(cè)矩陣Ht寫成

略去高次項(xiàng)，有如下關(guān)系：

定義新的觀測(cè)向量ηt，滿足

可將觀測(cè)方程(7b)改寫成

則式(14)可應(yīng)用線性系統(tǒng)的Kalman 濾波，稱為擴(kuò)展Kalman 濾波。

2 數(shù)值算例

2.1 未知邊界問題反演

如圖1 所示，以直角坐標(biāo)系中2a × 2h(a = 10，h = 10)矩形區(qū)域Ω 為例，在區(qū)域中函數(shù)φ 滿足Δ2φ= 0，邊界條件定義如下：

圖1 邊界問題區(qū)域Fig.1 Domain of the problem

式中φ(1)、φ(2)和φ(3)分別為已知邊界函數(shù)值。根據(jù)式(15)，該區(qū)域內(nèi)函數(shù)存在如下形式的解析解φana：

根據(jù)式(16)，區(qū)域Ω 在x 軸以下的邊界并未給出，此問題的邊界條件不完整，僅利用基本解方法計(jì)算不能得出準(zhǔn)確的數(shù)值解。

不失一般性，令式(3)中λ = 1，取4 個(gè)觀測(cè)點(diǎn)(如圖1 所示)，以觀測(cè)點(diǎn)處的解析解為測(cè)量值，利用Kalman 濾波方程(8a)—(8e)進(jìn)行計(jì)算，此時(shí)c 為4 ×4 單位矩陣，Ht為4×1 維向量。在具體計(jì)算過程中，定義測(cè)量后驗(yàn)估計(jì)和先驗(yàn)估計(jì)的相對(duì)誤差δ1為

函數(shù)計(jì)算值與解析解相對(duì)誤差δ2為

式中n 為所取測(cè)試點(diǎn)總數(shù)。

根據(jù)式(8a)—(8e)，Kalman 濾波過程為

2)根據(jù)基本解方法可以算出待定系數(shù)Ci，并由式(16)計(jì)算觀測(cè)點(diǎn)的實(shí)際值；

4)第2 次計(jì)算時(shí)，重復(fù)以上步驟，直到相對(duì)誤差δ1小于0.01%時(shí)結(jié)束計(jì)算。

對(duì)4 個(gè)觀測(cè)點(diǎn)以不同的初始預(yù)測(cè)值(-30和-50)進(jìn)行Kalman 濾波，結(jié)果如圖2 所示。從圖2可以看出：在不同初始預(yù)測(cè)值情況下，通過濾波均能收斂于解析解，且相對(duì)誤差δ1可控制在0.01%以內(nèi)；濾波速度快，濾波3 次后即可收斂至實(shí)際值，并在收斂過程中，具有良好的魯棒性。

圖2 未知邊界問題Kalman 濾波結(jié)果Fig.2 Kalman filtering results of unknown boundary problems

圖3(a)和(b)分別給出了研究區(qū)域內(nèi)不完備邊界y=-h 處的函數(shù)及其導(dǎo)數(shù)?φ/?n 的數(shù)值解和解析解的對(duì)比情況。由圖3 可知，通過Kalman 濾波，在未知邊界上函數(shù)數(shù)值解與解析解的相對(duì)誤差δ2在5%以內(nèi)，充分證明了Kalman 濾波技術(shù)在未知邊界反演問題中的有效性。

圖3 邊界y=-10 處計(jì)算值與解析解對(duì)比Fig.3 Comparison of calculated value and analytical value at boundary y=-10

2.2 未知參數(shù)反演

如圖1 所示，同樣以2a × 2h(a=10，h=10)矩形區(qū)域內(nèi)的二維Piosson 方程為例，進(jìn)行參數(shù)反演研究，為擴(kuò)展Kalman 問題。考慮到λ 為常數(shù)，取目標(biāo)參數(shù)ε=1/λ 進(jìn)行Kalman 濾波。本問題的邊界條件定義如下：

式中φ(4)、φ(5)、φ(6)和φ(7)分別為已知邊界函數(shù)值。則該區(qū)域內(nèi)φ 函數(shù)同樣存在解析解(16)，且常數(shù)目標(biāo)值ε=8(λ=1/8)。

同樣在區(qū)域內(nèi)取4 個(gè)觀測(cè)點(diǎn)，分別利用不同初始預(yù)測(cè)值ε=20 和ε=50 進(jìn)行濾波，結(jié)果如圖4(a)和(b)所示。由圖4 可以看出，參數(shù)的初始值在經(jīng)過兩次濾波之后就趨于實(shí)際值ε=8，且相對(duì)誤差δ1在0.01%以內(nèi)，表明該方法對(duì)未知參數(shù)反演問題有效。

圖4 未知參數(shù)問題Kalman 濾波結(jié)果Fig.4 Kalman filtering results for unknown parameter problems

利用反演參數(shù)對(duì)y=0 處的函數(shù)φ 及其導(dǎo)數(shù)進(jìn)行對(duì)比，結(jié)果如圖5 所示。由圖5 可以看出，在邊界y=0 上得到的函數(shù)及其導(dǎo)數(shù)的計(jì)算值與解析解相對(duì)誤差在0.05%以內(nèi)，證明了本方法的有效性。

圖5 y=0 處計(jì)算值與解析解邊界條件對(duì)比Fig.5 Comparison of calculated value and analytical value at boundary y=0

3 結(jié)論

本文采用基本解方法和Kalman 濾波技術(shù)相結(jié)合，對(duì)存在不完備邊界和未知參數(shù)的Poisson 方程的反問題進(jìn)行研究和計(jì)算，通過觀測(cè)點(diǎn)的測(cè)量值反演出未知的邊界值或參數(shù)值。通過數(shù)值算例可以看出，利用該方法得到的數(shù)值解與解析解之間的誤差小，收斂速度快，魯棒性強(qiáng)，顯示了本方法的有效性。該方法可為電磁學(xué)、滲流、熱傳導(dǎo)等工程領(lǐng)域具體問題的數(shù)學(xué)建模、數(shù)值分析提供參考，為反問題的分析研究提供思路。