李 冬，廖和琴，田 雨

(1.天津大學(xué) 電氣與自動(dòng)化工程學(xué)院，天津300072;2.國(guó)家海洋技術(shù)中心 海洋測(cè)量傳感器技術(shù)研究室，天津300112)

0 引 言

由于化石燃料燃燒、森林砍伐以及其他人為因素，大氣中的CO2濃度顯著提高［1］。海洋吸收了1980～1994 年之間排入大氣中CO2的28%～34%［2，3］。海洋表面pH 值較工業(yè)革命之前下降了0.1，預(yù)計(jì)到2100 年，海洋表層pH 還會(huì)下降0.3 ～0.5［4］。海洋酸化被認(rèn)為是對(duì)海洋生物系統(tǒng)尤其是貝殼類和鈣質(zhì)類生物的主要威脅［5］。海洋層pH 測(cè)量是研究海洋酸化最重要的手段，pH 傳感器的動(dòng)態(tài)特性是一個(gè)重要特性，對(duì)于測(cè)量方法和測(cè)量數(shù)據(jù)處理有重要意義。

SBE18 是由美國(guó)海鳥公司推出的海洋測(cè)量研究用pH傳感器。在經(jīng)常標(biāo)定的情況下，SBE18 pH 值測(cè)量精度能達(dá)到0.1，是海洋測(cè)量領(lǐng)域目前的精度最高，測(cè)量最穩(wěn)定的pH 傳感器之一。

CTD 傳感器中溫度傳感器和電導(dǎo)率傳感器動(dòng)態(tài)特性一直是研究的重點(diǎn)［6］，pH 傳感器動(dòng)態(tài)響應(yīng)也偶有涉及［7，8］。研究CTD 傳感器動(dòng)態(tài)特性時(shí)，通常將一階系統(tǒng)作為傳感器的動(dòng)態(tài)模型，將時(shí)間常數(shù)作為傳感器動(dòng)態(tài)特性的參數(shù)。實(shí)際測(cè)量發(fā)現(xiàn)SBE18動(dòng)態(tài)特性與一階系統(tǒng)相差甚遠(yuǎn)，時(shí)間常數(shù)無法完全描述SBE18 的動(dòng)態(tài)過程。本文提出了辨識(shí)SBE18 動(dòng)態(tài)模型的方法，辨識(shí)出了SBE18 的動(dòng)態(tài)模型，并分析了環(huán)境因素對(duì)傳感器動(dòng)態(tài)模型參數(shù)的影響。

1 傳感器動(dòng)態(tài)特性測(cè)試平臺(tái)與測(cè)試方法

傳感器動(dòng)態(tài)特性測(cè)試平臺(tái)槽體分為A，B 兩區(qū)，中間用密封隔板隔開。兩區(qū)均安裝有加熱器、熱交換器、感應(yīng)式溫度電導(dǎo)率傳感器、溶解氧傳感器、pH 傳感器，能分別實(shí)時(shí)監(jiān)控兩區(qū)水體的溫度、電導(dǎo)率、溶解氧和pH 等參數(shù)，并能將兩區(qū)水溫分別控制在設(shè)定值的±0.1 ℃偏差范圍內(nèi)。測(cè)量裝置的結(jié)構(gòu)框圖如圖1。

圖1 傳感器動(dòng)態(tài)響應(yīng)測(cè)量平臺(tái)結(jié)構(gòu)框圖Fig 1 Structure block diagram of sensor dynamic response measurement platform

測(cè)量前，將兩區(qū)水溫控制在20 ℃;分別加入適量鄰苯二甲酸氫鉀和硼砂緩沖劑;隔板打開后，A，B 兩區(qū)水體在隔板處形成厚度為2～5 cm，pH 梯度為80～200/m 的水層，該水層是測(cè)量傳感器階躍響應(yīng)的理想位置。測(cè)量時(shí)，打開密封隔板，同時(shí)傳感器在驅(qū)動(dòng)電機(jī)的牽引下以速度0.1，0.5，1.0，1.5 m/s 穿過pH 躍層，從A 區(qū)運(yùn)動(dòng)到B 區(qū);由于該躍層厚度不能完全忽略不計(jì)，這相當(dāng)于給傳感器一個(gè)非理想的階躍輸入。測(cè)量過程中傳感器輸出信號(hào)通過電纜傳遞到采集板卡，采集板卡將模擬信號(hào)轉(zhuǎn)換為數(shù)字信號(hào)，存儲(chǔ)在上位機(jī)上。上位機(jī)數(shù)據(jù)處理軟件計(jì)算出傳感器的動(dòng)態(tài)響應(yīng)參數(shù)。

2 SBE18 動(dòng)態(tài)特性參數(shù)測(cè)定原理

本文首先通過pH 傳感器的定標(biāo)公式，證明可以使用輸出電壓計(jì)算傳感器的動(dòng)態(tài)特性參數(shù)，然后給出了各階備選模型的時(shí)域微分方程，提出給定參數(shù)條件下求解各階模型時(shí)域微分方程的數(shù)值解法，最后闡述了使用遺傳算法求解最優(yōu)模型參數(shù)的方法。

2.1 SBE18 pH 傳感器定標(biāo)公式

SBE18 輸出電壓為0～5 V，其定標(biāo)公式為

其中，Vout為傳感器輸出電壓，pHoffset為pH 偏置電壓，pHslope為pH 擬合斜率，K 為以開氏溫標(biāo)計(jì)算的環(huán)境溫度，R氣體常數(shù)，F(xiàn) 為法拉第常數(shù)。由式(1)可計(jì)算出pH 值。pHoffset，pHslope，R，F(xiàn) 均為常數(shù)，在一次測(cè)量過程中，A，B 兩槽水體溫度設(shè)為同一值，pH 值與傳感器的輸出電壓呈線性關(guān)系，因此，本文中采用傳感器輸出電壓—時(shí)間曲線來測(cè)定傳感器的動(dòng)態(tài)響應(yīng)參數(shù)。

2.2 備選模型集時(shí)域微分方程

設(shè)傳感器的傳遞函數(shù)為H(s)，傳感器的輸入信號(hào)為X(s)、輸出信號(hào)為Y(s)，則Y(s)=H(s)·X(s)，在不考慮傳感器放大系數(shù)和輸入信號(hào)幅值的情況下，不妨設(shè)傳感器的開環(huán)放大系數(shù)為1，即H(0)=1，則

其中，H1，H2，H3分別為一階、二階、三階系統(tǒng)模型，τ1，τ2，τ3為其參數(shù);t0輸入信號(hào)躍變時(shí)刻，y(0)為傳感器輸出穩(wěn)定初值，y(∞)為傳感器輸出穩(wěn)定終值。H1，H2，H3對(duì)應(yīng)的時(shí)域微分方程分別為

2.3 給定參數(shù)備選模型時(shí)域解數(shù)值解法

本文采用Bogacki-Shampine 算法求解時(shí)域微分方程。Bogacki-Shampine 方法是一種具有前后一致特性的三階龍格庫塔法，內(nèi)嵌一個(gè)二階龍格庫塔法，可實(shí)現(xiàn)變步長(zhǎng)計(jì)算［9，10］，在容許誤差較大時(shí)，求解效果好于Dormand-Prince算法［11］。對(duì)于任意一組τ1τ2τ3t0，使用Bogacki-Shampine算法求解式(7)～式(9)，即可得到各階模型的時(shí)域數(shù)值解yc(tc(i))，其中，tc是時(shí)間序列，yc是幅值序列。該算法的基本問題是在時(shí)間區(qū)間［t0，tf］內(nèi)求解初值問題y'=F(t，y)，y(t0)=y0。Bogacki-Shampine 算法在［tn，tn+1］之間的遞推公式如下

其中，yn為初值問題在tn處的數(shù)值解，hn為步長(zhǎng)，即hn=tn+1－tn。這里，zn+1為精確解的二階近似，yn+1為精確解的三階近似，因此，zn+1和yn+1的差值可用作變步長(zhǎng)。

2.4 模型誤差判定

本文將模型誤差分為幅值誤差Ev和階躍時(shí)刻誤差Et。幅值誤差Ev為傳感器實(shí)際輸出序列y(t(i))與模型時(shí)域數(shù)值解yc(tc(i))在時(shí)間區(qū)間［t0.05，t0.95］內(nèi)的均方差和傳感器輸出變化幅值的比值，即

階躍時(shí)刻誤差Et為模型階躍開始時(shí)刻t0與實(shí)際階躍開始時(shí)刻tr0之差值的絕對(duì)值，即

其中，實(shí)際階躍開始時(shí)刻tr0為穩(wěn)定初值f(0)與實(shí)際輸出曲線的最后一個(gè)交點(diǎn)。

2.5 備選模型最優(yōu)參數(shù)優(yōu)化方法

本文使用Matlab 遺傳算法求解各階系統(tǒng)的最優(yōu)τ1τ2τ3t0［12，13］。遺傳算法是一種受進(jìn)化論啟發(fā)的優(yōu)化算法，不同于確定性求解方法，遺傳算法通過受控的隨機(jī)方式向最優(yōu)解趨近［14］，非常適合求解空間大，峰值單一，系統(tǒng)模型復(fù)雜，不要求絕對(duì)全局最優(yōu)的問題。

遺傳算法求解的基本問題是

3 實(shí)驗(yàn)數(shù)據(jù)分析處理

SBE18 pH 傳感器輸出信號(hào)為電壓信號(hào)，容易受到電磁輻射和工頻交流電的干擾，本文使用軟件對(duì)原始數(shù)據(jù)進(jìn)行50 Hz 理想低通濾波。

3.1 最佳傳感器動(dòng)態(tài)模型的確定

為了評(píng)價(jià)各階模型的適應(yīng)性，對(duì)同一組傳感器輸出的數(shù)據(jù)分別使用各階模型進(jìn)行優(yōu)化，各階模型優(yōu)化曲線與實(shí)際輸出曲線進(jìn)行對(duì)比如圖2。評(píng)價(jià)模型優(yōu)劣指標(biāo)有三個(gè)，1)優(yōu)化后模型的幅值誤差Ev;2)優(yōu)化后模型的階躍時(shí)刻誤差Et;3)模型參數(shù)數(shù)量N。Ev表征了實(shí)際輸出曲線與優(yōu)化曲線在取樣區(qū)間內(nèi)的重合度;Et表征了實(shí)際曲線與優(yōu)化曲線階躍開始時(shí)刻的差值，由于實(shí)際非理想階躍輸入的影響，該差值是必然存在的;模型參數(shù)數(shù)量N 決定了模型的復(fù)雜度，求取參數(shù)所需的時(shí)間與應(yīng)用參數(shù)的難度。各階系統(tǒng)適應(yīng)性評(píng)價(jià)指標(biāo)如表1 所示。

圖2 各階系統(tǒng)優(yōu)化曲線與實(shí)際輸出曲線對(duì)比Fig 2 Comparison of optimized curve and real output curve of each order system

表1 使用不同傳感器模型的優(yōu)化結(jié)果誤差Tab 1 Optimized result error of different kind of models

由以上數(shù)據(jù)分析可知，幅值誤差Ev和階躍時(shí)刻誤差Et都隨模型的階數(shù)增大而減小，但模型階數(shù)越高確定模型所需的參數(shù)N 也會(huì)越多;二，三階模型的幅值誤差Ev和階躍時(shí)刻誤差Et遠(yuǎn)小于一階系統(tǒng);三階系統(tǒng)模型的幅值誤差Ev略小于二階系統(tǒng)，兩系統(tǒng)的幅值誤差Ev都足夠小。求解復(fù)雜模型參數(shù)耗費(fèi)時(shí)間多，同時(shí)復(fù)雜模型也會(huì)給使用帶來不便，因此，本文將二階系統(tǒng)選定為SBE18 的動(dòng)態(tài)響應(yīng)模型。

3.2 環(huán)境因素對(duì)模型參數(shù)的影響

實(shí)際傳感器的動(dòng)態(tài)特性參數(shù)不是一組固定值，它會(huì)隨著傳感器的使用環(huán)境改變而改變。影響pH 傳感器動(dòng)態(tài)特性的因素主要有兩個(gè):1)傳感器的附近水流速度;2)水體pH 梯度變化方向。本文中分別在0.1，0.5，1.0，1.5 m/s 四個(gè)速度和pH 由高到低、由低到高兩個(gè)方向上求取了傳感器模型的參數(shù)，結(jié)果如表2 所示。

表2 不同流速和pH 變化方向下的模型參數(shù)Tab 2 Model parameters under different velocity of flow and pH change direction

隨著流速增大模型參數(shù)均減小，且流速大于1.0 m/s以后傳感器參數(shù)趨于穩(wěn)定，因此，維持1 m/s 左右的流速有助于傳感器輸出快速穩(wěn)定。傳感器的動(dòng)態(tài)參數(shù)受傳感器輸入信號(hào)變化方向影響十分顯著，pH 值由高到低變化時(shí)，傳感器的動(dòng)態(tài)特性參數(shù)明顯大于當(dāng)pH 值由低到高變化時(shí)的數(shù)值。

4 結(jié) 論

SBE18 的參數(shù)指標(biāo)里僅標(biāo)明了特定流速下時(shí)間常數(shù)，沒有標(biāo)明該值的具體測(cè)量環(huán)境，而且時(shí)間常數(shù)也無法描述傳感器的具體動(dòng)態(tài)過程，這給不少相關(guān)海洋科技工作者帶來很大的不便。本文深入研究了SBE18 動(dòng)態(tài)參數(shù)的測(cè)量方法，分別分析了一，二，三階系統(tǒng)模型對(duì)傳感器動(dòng)態(tài)過程的適應(yīng)性，得出二階系統(tǒng)模型與SBE18 動(dòng)態(tài)特性符合較好;同時(shí)，求取不同流速和pH 梯度變化方下傳感器二階系統(tǒng)模型參數(shù)，并分析了兩者對(duì)模型參數(shù)的影響。通過大量的實(shí)驗(yàn)數(shù)據(jù)分析得出，模型參數(shù)隨流速增大而減小，且流速大于1m/s 后傳感器動(dòng)態(tài)特性參數(shù)趨于穩(wěn)定;pH 梯度變化方向?qū)δＰ蛥?shù)影響大，pH 由高到低變化時(shí)傳感器動(dòng)態(tài)特性參數(shù)明顯大于pH 由低到高時(shí)的變化。

［1］ Wootton J T，Pfister C A，F(xiàn)orester J D.Dynamic patterns and ecological impacts of declining ocean pH in a high-resolution multiyear dataset［C］∥Proceedings of the National Academy of Sciences，2008:18848－18853.

［2］ Millero F J.The marine inorganic carbon cycle［J］.Chemical Reviews，2007，107(2):308－341.

［3］ Sabine C L，F(xiàn)eely R A，Gruber N，et al.The oceanic sink for anthropogenic CO2［J］.Science，2004，305(5682):367－371.

［4］ Caldeira K，Wickett M E.Oceanography:Anthropogenic carbon and ocean pH［J］.Nature，2003，425(6956):365－365.

［5］ Caldeira K，Wickett M E.Ocean model predictions of chemistry changes from carbon dioxide emissions to the atmosphere and ocean［J］.Journal of Geophysical Research:Oceans，2005，110(C9):1978－2012.

［6］ Gregg M C，Hess W C.Dynamic response calibration of sea-bird temperature and conductivity probes［J］.Journal of Atmospheric and Oceanic Technology，1985，2(3):304－313.

［7］ Oelbner W，Zosel J，Berthold F，et al.Investigation of the dynamic response behaviour of ISFET pH sensors by means of laser Doppler velocimetry(LDV)［J］.Sensors and Actuators B:Chemical，1995，27(1):345－348.

［8］ Hara H，Ohta T.Dynamic response of a Ta2O5-gate pH-sensitive field-effect transistor［J］.Sensors and Actuators B:Chemical，1996，32(2):115－119.

［9］ Bogacki Przemyslaw，Shampine Lawrence F.A 3(2)pair of Runge-Kutta formulas［J］.Applied Mathematics Letters，1989，2(4):321－325.

［10］Shampine L F，Reichelt M W.The Matlab ode suite［J］.SIAM Journal on Scientific Computing，1997，18(1):1－22.

［11］Dormand J R，Prince P J.A family of embedded Runge-Kutta formulae［J］.Journal of Computational and Applied Mathematics，1980，6(1):19－26.

［12］Conn A R，Gould N I M，Toint P.A globally convergent augmented Lagrangian algorithm for optimization with general constraints and simple bounds［J］.SIAM Journal on Numerical Analysis，1991，28(2):545－572.

［13］Conn A，Gould N，Toint P.A globally convergent Lagrangian barrier algorithm for optimization with general inequality constraints and simple bounds［J］.Mathematics of Computation of the American Mathematical Society，1997，66(217):261－288.

［14］Child B F M，Venugopal V.Optimal configurations of wave energy device arrays［J］.Ocean Engineering，2010，37(16):1402－1417.