楊錚園， 申立勇

(中國科學院大學數(shù)學科學學院，北京101408)

數(shù)據(jù)點擬合在計算機輔助幾何設計、計算機圖形學和數(shù)控加工等領(lǐng)域都有廣泛的應用，擬合中較常使用B樣條曲線。B樣條曲線具有較好的造型能力[1-2]。在擬合算法及提高精確度和效率方面，前人已做了許多研究。通常把平方距離和作為擬合誤差，這樣擬合過程就可以描述成一個求最小擬合誤差的最優(yōu)化問題。在優(yōu)化迭代中，怎樣估計擬合誤差是一個關(guān)鍵的問題。目前已有很多估計擬合誤差的方法，包括PDM，TDM，SDM，GTDM和CDM[3-4]等，在擬合效果和速度上都各有優(yōu)缺點。Yang等[5]提出了隱式B樣條曲線擬合算法，適用于更一般的情形，比如幾條分離的閉曲線的擬合。Fl?ry[6]提出了在障礙存在下的B樣條曲線曲面擬合算法，該算法能夠進行點云邊界的重構(gòu)，并推廣到避開一般真實障礙物的曲線曲面重構(gòu)。本文針對實際實驗數(shù)據(jù)，采用多約束條件下的樣條擬合，處理油氣爆炸界限的估計問題，得到不同氧氣濃度下油氣濃度和爆炸強度曲線，應用于燃爆預判。

近年來，工業(yè)生產(chǎn)過程危險化學品事故頻繁發(fā)生，事故造成了生命和財產(chǎn)的重大損失。因此，危險化學品事故的預測是眾多學者積極研究的主題[7-9]。Chang和Lin[10]整理了在2006年之前的40年間有記載的工廠事故，其中約74%的事故發(fā)生在石油精煉廠、石油港口和油庫。近年來這類事故的數(shù)量還在上升，包括2010年大連港油庫爆炸事故。油罐爆炸可能會導致多米諾骨牌效應，加劇意外事故所造成的損害[11]。油罐爆炸主要是油氣燃爆引起的，因此對油氣爆炸范圍進行分析非常必要，良好的分析結(jié)果將有助于防止油氣爆炸事故的發(fā)生。之前人們主要應用不同分類器對某種油氣的實驗數(shù)據(jù)進行分析，通過訓練得到合適的分類器參數(shù)，從而可以對新數(shù)據(jù)進行是否屬于可以燃爆范圍的預判。常用的分離器包括支持向量機(SVM)方法[12]、邏輯回歸(LR)方法[13]等。這類基于統(tǒng)計分析的分類器一般適合給出邏輯變量分類，即是否可燃爆。例如目前常用的燃氣容積比作為因變量，而它的可燃性作為目標變量。但是實際生產(chǎn)中人們不僅關(guān)注燃氣是否爆炸，還關(guān)注其爆炸程度強弱，以便對后援工作進行有效安排。因此，本文基于此實際需求，對某油罐燃氣可燃性以及燃爆強度進行分析，其影響因素涉及其中的油氣濃度和氧氣濃度等。我們通過在特制的密閉管道內(nèi)進行油氣爆炸試驗，從而獲取在常溫常壓條件下不同油氣濃度、氧氣濃度預混氣體的爆炸極限及爆炸強度值(圖1)。然后分析不同氧氣濃度下的油氣燃爆界，以及可燃爆情況下的爆炸強度。

圖1 實驗裝置示意圖(單位：mm)

本文在關(guān)注油氣數(shù)據(jù)擬合中，首先根據(jù)問題的物理化學背景和實驗數(shù)據(jù)特征，對擬合曲線的預期性質(zhì)進行估計。然后把這些性質(zhì)轉(zhuǎn)化為曲線擬合中的約束條件，對擬合曲線進行帶約束的優(yōu)化。包括，不同氧氣濃度下的爆炸強度曲線應該互不相交；爆炸強度曲線是凸曲線；爆炸強度曲線應該為安全曲線等。這些要求在轉(zhuǎn)化為約束條件時又有不同的轉(zhuǎn)化和實現(xiàn)方式，需要不同的嘗試。本文首先回顧了B樣條曲線擬合的一般過程、SDM誤差項估計方法和Fl?ry的障礙約束方法，然后根據(jù)數(shù)據(jù)特征和問題背景導出曲線組擬合的約束條件，最后在實際實驗數(shù)據(jù)上實現(xiàn)算法，在分析、比較擬合效果后選擇合適的方案。

1 實驗數(shù)據(jù)和樣條擬合

1.1 實驗數(shù)據(jù)采集

在實驗裝置中(圖1)，首先測量實驗油氣(HC)濃度、氧氣(O2)濃度和惰性氣體(CO2)濃度，惰性氣體主要是用來調(diào)整前兩種實驗氣體濃度。然后通過點火裝置，再由傳感器測量并記錄引爆后的各自氣體濃度，包括油氣、氧氣、一氧化碳、二氧化碳。在3個位置安裝了濃度傳感器，將數(shù)據(jù)平均值作為最終實驗數(shù)據(jù)進行分析。另外，還安裝了5個壓強傳感器，測量爆炸瞬間的壓強，其平均值作為實驗壓強值，該瞬時壓強稱之為燃氣爆炸強度。而燃氣燃爆界是指固定氧氣濃度下油氣發(fā)生爆炸的濃度邊界，通常用上下界區(qū)間表示。表1和表2分別給出了引爆后發(fā)生爆炸和沒有爆炸的兩組實驗數(shù)據(jù)；表3為爆炸瞬間壓強表。

整個實驗過程中，共采集、整理了44組數(shù)據(jù)(圖2)。經(jīng)初步分析發(fā)現(xiàn)，隨著氧氣濃度的增大，油氣爆炸范圍增大、爆炸強度也增強。因此，將數(shù)據(jù)按氧氣濃度分類，然后對每一組數(shù)據(jù)進行擬合。但是每組數(shù)據(jù)擬合不是獨立的，不同氧氣濃度下爆炸曲線有相互約束關(guān)系。

表1 爆炸實驗數(shù)據(jù)(No.21)

表2 未爆炸實驗數(shù)據(jù)(No.19)

表3 爆炸瞬間壓強(MPa)

圖2 實驗數(shù)據(jù)

1.2 B樣條曲線擬合

設f:?→?2:t?f(t)是平面上的一條參數(shù)曲線。已知有n個控制點的B樣條曲線的一般表示為：

其中，Bk,i(t)∈?是k次B樣條基函數(shù)，是曲線的控制點。令：

則式(1)可以寫成矩陣表示：

給出平面上的一組無序數(shù)據(jù)點pk,k=1,…,m，用一條B樣條曲線來擬合這些點，使得這條曲線盡可能地描述這些點的軌跡形狀。這一過程可以描述為如下一個非線性優(yōu)化問題：已知點集求一組控制點c,i=1,…,n,∈?2，使i得目標函數(shù)：

的值最小。其中|d(pk,f)|是pk到曲線f(t)的距離，ωfs是平滑項。

式(3)的最優(yōu)化問題可以用迭代過程解決。我們不妨設f(tk)是已知曲線f(t)上距離pk最近的點，pk到曲線f(t)的距離定義為|d(pk,f(tk))|。對控制點向量C做位移為δC的擾動，則原始曲線f(t)變?yōu)椋?/p>

pk到曲線的距離未知，但是可以用合理的方法來估計。優(yōu)化問題的關(guān)鍵在于如何估計曲線的擬合誤差d2(pk,）。Wang等[3]提出的SDM方法給出的誤差項：

其中，Tk和Nk分別是f(tk)點的單位切向量和單位法向量，ρk是f(tk)點的曲率半徑。|dk|=||pk-f(tk)||，當ρk與曲率中心K在曲線f(t)同側(cè)時dk＞0，異側(cè)時dk＜0。顯然dk＜ρk。式(5)可展開計算并求和，

其中：

另一方面，平滑項部分在實際應用中常用最小化薄板能量，使用二階導數(shù)的平方[14]，因此fs也是C的二次函數(shù)：

于是式(3)的無約束最優(yōu)化問題等價于在現(xiàn)有值C的情況下，每次迭代中最小化如下函數(shù)：

根據(jù)以上記號，B樣條曲線擬合算法可描述如下：

輸入：需要擬合的一組二維無序數(shù)據(jù)點pk,k=1,…,m,∈?2。

輸出：一條B樣條擬合曲線f(t)。

(1) 根據(jù)所要擬合的數(shù)據(jù)點，合理設置初始控制點,i=1,…,n，根據(jù)式(2)計算出初始的B樣條曲線f0(t)。

(2) 對每個數(shù)據(jù)點pk，計算其在當前的曲線fr(t)上的最近點fr(tk)。

(3) 通過求解式(8)的最小值，計算擬合曲線控制點Cr的位移δCr。則新的擬合曲線為fr+1(t)=LT(t)(Cr+δCr)。如果新的擬合曲線滿足了事先設定的條件（比如迭代步數(shù)或者誤差閾值的控制），則轉(zhuǎn)到第4步；否則轉(zhuǎn)到第2步。

(4) 輸出擬合曲線fr+1(t)。

在樣條擬合的一般算法的基礎上，F(xiàn)l?ry[6]提出了在一般障礙物存在的情況下擬合點云的方法。

同上文類似的記號，并記f()處的向外單位法向量n)，則式(9)的約束等價于：

由于每個數(shù)據(jù)點對應的與當前的擬合曲線的控制點相關(guān)，故可寫成f(t)的控制點向量C的函數(shù)，并記式(10)左側(cè)的線性近似為′(C)，則：

因此，約束(10)可以近似寫成：

′(C)對每個控制點ci求偏導：

其中，I是二階單位矩陣，J是將切向量向曲線外旋轉(zhuǎn)π/2的矩陣。經(jīng)整理后可寫出梯度矩陣?Clk′(C)，式(12)的具體推導過程可參見文獻[6]。

2 多約束的B樣條曲線擬合

根據(jù)爆炸數(shù)據(jù)的物理化學特征，將不同的擬合約束轉(zhuǎn)化為合理的數(shù)學條件，并對條件做進一步簡化。

2.1 多條互不相交的曲線擬合

爆炸反應的原理顯示，在相同的油氣濃度下，爆炸強度的大小跟氧氣濃度有關(guān)。不同氧氣濃度下的爆炸強度曲線應該互不相交。因此在實際應用于爆炸數(shù)據(jù)的擬合中，需要針對曲線性質(zhì)的特殊要求，對擬合曲線進行帶約束的優(yōu)化。對于本文的實驗數(shù)據(jù)，不能直接應用上一節(jié)的方法。因為Fl?ry的算法處理的障礙物都是處于靜止狀態(tài)下，擬合的多條曲線分別作為其他曲線的邊界約束在每一次迭代中會發(fā)生變化，因此約束曲線上取的樣點作為障礙點也在變化。

上述問題的擬合可以用兩種方法來完成。首先很直接地想到第一種方法：先擬合其中一條曲線，比如f2(t2)，然后擬合另一條曲線f1(t1)時，在f2(t2)上取足夠多的樣點作為障礙點集，即可用上一節(jié)的帶約束的算法來擬合f1(t1)。這樣的做法能夠滿足曲線互不相交，但是曲線擬合順序的不同可能導致產(chǎn)生不同的擬合效果，尤其在兩組點云重疊部分會有明顯的偏差，較先擬合的曲線對擬合結(jié)果的影響較大。

第二種方法旨在解決第一種方法中可能出現(xiàn)的問題，避免因曲線擬合的順序不同產(chǎn)生較大的偏差，對{p1,k1}和{p2,k2}同時進行擬合。把兩個控制點向量C1和C2寫成一個列向量，第r步迭代中，在當前曲線f2(t2)上取足夠多的樣點f2(t2,k′),k′=1,…,m′。下一步迭代中由于曲線f2(t2)的變化，f2(t2,k′)也隨之移動，記為(t2,k′)，要求(t2,k′)都在(t1)的內(nèi)側(cè)。于是可轉(zhuǎn)化為解問題：

記號定義同上文。曲線f2(t2)上的樣點f2(t2,k′)可以寫成C2的函數(shù)，故是C1和C2的函數(shù)。對屬于C1的控制點的偏導可參照式(12)，對屬于C2的控制點c2,i的偏導：

B2,i(t2)是f2(t2)的B樣條基函數(shù)。于是式(13)的約束可寫成：

多組數(shù)據(jù)的相互約束擬合可基于上述兩組數(shù)據(jù)擬合的約束條件。lk對不與當前擬合曲線有邊界約束關(guān)系的曲線的控制點的偏導均為0，因此整理后也可寫成如(14)的線性不等式約束。

2.2 凸曲線的擬合

根據(jù)燃氣爆炸實驗的特點，在氧氣濃度固定的情況下，油氣濃度達到爆炸區(qū)間的某一值時，爆炸強度可達到最大，爆炸強度曲線應是凸曲線。在上一小節(jié)的的約束擬合中加入新的約束用以保證擬合曲線的上凸性質(zhì)。不妨設擬合函數(shù)f(t)需要保持向外凸的性質(zhì)，故在曲線的任一點處的曲率中心均在曲線的內(nèi)側(cè)。這一條件也可敘述為f′′(t)始終與f(t)的向內(nèi)法向量同向。沿用上文的記號，J是將f′(t)向曲線外旋轉(zhuǎn)π/2的矩陣，則JT是將f′(t)向曲線內(nèi)旋轉(zhuǎn)π/2的矩陣。曲線f(t)=L(t)TC應滿足：

寫成控制點向量C的約束：

上式的左側(cè)發(fā)生擾動δC后在C處展開并線性化：

式(16)就是保持擬合曲線外凸性質(zhì)的約束條件。在實際應用中，可以在每一次迭代中在f(t)上取足夠多的樣點，用式(16)的約束求解子問題的最優(yōu)化問題。

上述條件實現(xiàn)中，算法需要在f(t)上取足夠多的樣點，因此存在如何控制所取樣點數(shù)目的問題。下面我們根據(jù)樣條的變差縮減性質(zhì)直接給出控制點約束，從而導出線性的簡化條件。設擬合函數(shù)f(t)是外凸，由B樣條曲線變差縮減性，曲線的凸性受控于其控制多邊形凸性。如前文，J是將向曲線外旋轉(zhuǎn)π/2的矩陣，則JT=-J是向曲線內(nèi)旋轉(zhuǎn)π/2的矩陣。外凸曲線f(t)應滿足：

上式的左側(cè)發(fā)生擾動δC后在C處展開并線性近似，設ai=ci+1-ci,i=1,…,n-1，整理后得：

式(18)就是用變差縮減性質(zhì)保證曲線外凸的約束條件，只使用了控制點之間的關(guān)系，比直接用曲線上點的曲率構(gòu)成的式(16)簡單得多。

2.3 外邊界約束的曲線擬合

油氣爆炸曲線擬合的目的在于對有爆炸危險的濃度預警，在應用中可能出現(xiàn)兩種錯誤：一種是可燃爆數(shù)據(jù)被判定為安全；另一種則是不可燃爆數(shù)據(jù)被判定為可燃爆。在實際生產(chǎn)中，前一種錯誤顯然更加危險，一旦誘發(fā)燃爆，帶來的損失將遠遠大于后一種錯誤。因此，爆炸強度曲線應為安全曲線，即至少保證安全曲線之外沒有任何燃爆數(shù)據(jù)點。另一方面，還需考慮預警的準確性，即希望擬合曲線盡量靠近實驗數(shù)據(jù)。綜合兩個要求，需對爆炸數(shù)據(jù)作外邊界擬合。對于每一組數(shù)據(jù)的擬合曲線，把這組數(shù)據(jù)作為障礙點約束用于擬合數(shù)據(jù)點的外邊界。數(shù)據(jù)點的邊界擬合是Fl?ry算法的一種特殊情形，通過設置d(pk,f(t))＜0或d(pk,f(t))＞0即可擬合出的外邊界或內(nèi)邊界。約束條件應為：

3 數(shù)據(jù)擬合結(jié)果與分析

通過油氣爆炸實驗中獲取了四組數(shù)據(jù)點，分別記為{p20},{p18},{p16},{p15}，表示實驗反應在氧氣濃度為20%、18%、16%、15%采集到的爆炸壓強。用橫坐標表示爆炸反應前氣體中油氣的濃度，縱坐標表示爆炸時測得的平均壓強，需要用四條B樣條曲線f20(t),f18(t),f16(t),f15(t)分別擬合這四組爆炸數(shù)據(jù)。首先用B樣條曲線擬合的一般方法分別擬合四組數(shù)據(jù)點，可以粗略地看一下爆炸曲線的走勢，得到的擬合結(jié)果如圖3。

圖3 一般B樣條擬合

從圖3可以看出，曲線無凸性且兩條曲線有交點，不符合燃氣燃爆特性。依據(jù)曲線不相交約束的擬合算法，使用第一種方法，先對其中一組數(shù)據(jù)點做曲線擬合，然后對其他數(shù)據(jù)點分別做擬合時，在已擬合好的曲線上取若干點作為邊界約束。用第一種方法得到自上而下的約束擬合圖4(a)和自下而上的約束擬合圖4(b)，相應的曲線擬合優(yōu)先順序分別是f20(t)＞f18(t)＞f16(t)＞f15(t)和f15(t)＞f16(t)＞f18(t)＞f20(t)。從圖中可以看出，由于曲線擬合順序的不同，圖4(a)和圖4(b)有不同的擬合效果，這里較明顯的不同表現(xiàn)在曲線f18(t)和f16(t)上。與圖3沒有約束情況下的擬合結(jié)果相比，顯然在圖4(a)中，f18(t)和圖3中相應的曲線形狀大致相同，而f16(t)則因為f18(t)的限制而有明顯的變化，在圖4(b)中則相反。這一結(jié)果驗證了先擬合的曲線對整體擬合結(jié)果影響較大，而后擬合的曲線對先擬合的曲線沒有約束效果。

圖4 帶邊界約束的爆炸數(shù)據(jù)擬合

用第二種方法，式(14)，四條擬合曲線相互約束，互為上/下邊界，擬合結(jié)果如圖5。與圖4(a)和圖4(b)的按順序分別擬合的效果相比，圖5的相互擬合在直觀上有折中的效果，尤其在f18(t)和f16(t)上表現(xiàn)得比較明顯。從原理上來說，實驗測得的數(shù)據(jù)由于一些未知原因產(chǎn)生誤差，因此不能確定哪一組數(shù)據(jù)比較準確或者靠近理論值。如果較先擬合的數(shù)據(jù)不準確，那么整體的擬合效果就會受到較大的影響。因此，采用第二種方法將曲線彼此作為約束條件能夠較好地平衡這種關(guān)系，減小由于測量誤差導致的擬合誤差的增大。

圖5 互為邊界約束的爆炸數(shù)據(jù)擬合

從理論上看，固定氧氣的濃度，實驗爆炸時的壓強對于油氣的濃度作圖會呈現(xiàn)上凸的形狀。前面的擬合曲線都呈現(xiàn)出波動的狀態(tài)，要避免這種波動可以加上上凸曲線的約束。如果只用式(18)的約束來做擬合，得到圖6(a)，與圖3類似，曲線之間有交點。圖6(b)在式(18)的基礎上加上式(14)，更加符合理論上對爆炸壓強擬合曲線走向的預期。

為了得到爆炸強度的安全曲線，我們對四組數(shù)據(jù)點做相互約束的外邊界擬合，得到圖7(a)，直觀上是將爆炸點都包在曲線內(nèi)部。但如前討論，曲線呈現(xiàn)出不合理波動狀態(tài)。加上了上凸約束之后的圖7(b)與圖6(b)相比，圖7(b)的爆炸強度曲線形狀相似度更高一些。綜合考慮到擬合曲線的準確性、曲線形狀與真實情況的比較、以及危險預警的安全性，圖7(b)是本文推薦的爆炸強度曲線的擬合結(jié)果，當氧氣濃度為15%、16%、18%、20%時，其對應的油氣燃爆界限分別為[0.98,1.81]、[0.95,1.91]、[0.75,2.53]、[0.72,3.02]。值得注意的是，具有凸性的安全曲線之間沒有交點，因此這組實驗數(shù)據(jù)擬合中，不必添加互為邊界的約束。

圖6 上凸約束的爆炸數(shù)據(jù)擬合

圖7 爆炸強度的安全曲線擬合

4 結(jié)束語

本文針對某燃氣燃爆實驗數(shù)據(jù)，給出了多約束條件樣條擬合的一個應用實例。該實例本質(zhì)上是給定數(shù)據(jù)的多條件曲線擬合問題，本文主要在于約束條件數(shù)學化、條件合理簡化以及擬合算法的實現(xiàn)。經(jīng)過分析和比較，我們選擇具有上凸性質(zhì)的安全曲線擬合為最終分析方案。擬合得到安全曲線和爆炸界限都比較合理，并且得到合作消防公司的認可。更為重要的是，該方法可以推廣到多因素燃爆模型。實際工程生產(chǎn)中，燃爆影響因素除了本文涉及的各種混合氣體濃度和壓強，還可以有溫度、濕度等。添加新的影響因素，在我們的分析方法中轉(zhuǎn)化為增加樣條擬合維數(shù)，即多約束曲線擬合轉(zhuǎn)化為(超)曲面擬合。

[1]Pottmann H,Leopoldseder S,Hofer M.Approximation with active B-spline curves and surfaces [C]//10th Pacific Conference on Computer Graphics and Applications,Proceedings,2002: 8-25.

[2]Shen Liyong.Error Bounded Conic Spline Approximation for NC Code [C]//Proc.of International Conference on Machine Vision (ICMV2011),SPIE 8349,2011: 1-7.

[3]Wang Wenping,Pottman H,Liu Yang.Fitting B-spline curves to point clouds by curvature-based squared distance minimization [J].ACM Transactions on Graphic,2006,25(2): 214-238.

[4]Liu Yang,Wang Wenping.A revisit to least squares orthogonal distance fitting of parametric curves and surfaces [J].Lecture Notes Computer Science,2008,4975: 384-397.

[5]Yang Zhouwang,Deng Jiansong,Chen Falai.Fitting unorganized point clouds with active implicit B-spline curves [J].The Visual Computer,2005,21(8-10):831-839.

[6]Fl?ry S.Fitting curves and surfaces to point clouds in the presence of obstacles [J].Computer Aided Geometric Design,2009,26(2): 192-202.

[7]Rasbash D.Review of explosion and fire hazard of liquefied petroleum gas [J].Fire Safety Journal,1980,2:223-236.

[8]Bakke J R,Van Wingerden K,Hoorelbeke P,Brewerton B.A study on the effect of trees on gas explosions [J].Journal of Loss Prevention in the Process Industries,2010,23(6): 878-884.

[9]Cheng Jianwei,Yang Shengqiang.Improved coward explosive triangle for determining explosibility of mixture gas [J].Process Safety and Environmental Protection,2011,89(2): 89-94.

[10]Chang J I,Lin Chengchung.A study of storage tank accidents [J].Journal of Loss Prevention in the Process Industries,2006,19: 51-59.

[11]CCPS.Guidelines for Chemical Process Quantitative Risk Analysis [M].2nd ed.Wiley-AIChE,1999: 57-283.

[12]Vapnik V.The Nature of Statistical Learning Theory [M].NY,Springer,1999: 137-170.

[13]王濟川,郭志剛.Logistic回歸模型——方法與應用[M].北京: 高等教育出版社,2001: 146-181.

[14]Brunnet G,Hagen H,Santarelli P.Variational design of curves and surfaces [J].Surveys on Mathematics for Indestry,1993,3: 1-27.