SBE18 pH 傳感器動態特性研究*

李 冬，廖和琴，田 雨

(1.天津大學 電氣與自動化工程學院，天津300072;2.國家海洋技術中心 海洋測量傳感器技術研究室，天津300112)

0 引 言

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

SBE18 是由美國海鳥公司推出的海洋測量研究用pH傳感器。在經常標定的情況下，SBE18 pH 值測量精度能達到0.1，是海洋測量領域目前的精度最高，測量最穩定的pH 傳感器之一。

CTD 傳感器中溫度傳感器和電導率傳感器動態特性一直是研究的重點［6］，pH 傳感器動態響應也偶有涉及［7，8］。研究CTD 傳感器動態特性時，通常將一階系統作為傳感器的動態模型，將時間常數作為傳感器動態特性的參數。實際測量發現SBE18動態特性與一階系統相差甚遠，時間常數無法完全描述SBE18 的動態過程。本文提出了辨識SBE18 動態模型的方法，辨識出了SBE18 的動態模型，并分析了環境因素對傳感器動態模型參數的影響。

1 傳感器動態特性測試平臺與測試方法

傳感器動態特性測試平臺槽體分為A，B 兩區，中間用密封隔板隔開。兩區均安裝有加熱器、熱交換器、感應式溫度電導率傳感器、溶解氧傳感器、pH 傳感器，能分別實時監控兩區水體的溫度、電導率、溶解氧和pH 等參數，并能將兩區水溫分別控制在設定值的±0.1 ℃偏差范圍內。測量裝置的結構框圖如圖1。

圖1 傳感器動態響應測量平臺結構框圖Fig 1 Structure block diagram of sensor dynamic response measurement platform

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

2 SBE18 動態特性參數測定原理

本文首先通過pH 傳感器的定標公式，證明可以使用輸出電壓計算傳感器的動態特性參數，然后給出了各階備選模型的時域微分方程，提出給定參數條件下求解各階模型時域微分方程的數值解法，最后闡述了使用遺傳算法求解最優模型參數的方法。

2.1 SBE18 pH 傳感器定標公式

SBE18 輸出電壓為0～5 V，其定標公式為

其中，Vout為傳感器輸出電壓，pHoffset為pH 偏置電壓，pHslope為pH 擬合斜率，K 為以開氏溫標計算的環境溫度，R氣體常數，F 為法拉第常數。由式(1)可計算出pH 值。pHoffset，pHslope，R，F 均為常數，在一次測量過程中，A，B 兩槽水體溫度設為同一值，pH 值與傳感器的輸出電壓呈線性關系，因此，本文中采用傳感器輸出電壓—時間曲線來測定傳感器的動態響應參數。

2.2 備選模型集時域微分方程

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

其中，H1，H2，H3分別為一階、二階、三階系統模型，τ1，τ2，τ3為其參數;t0輸入信號躍變時刻，y(0)為傳感器輸出穩定初值，y(∞)為傳感器輸出穩定終值。H1，H2，H3對應的時域微分方程分別為

2.3 給定參數備選模型時域解數值解法

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

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

2.4 模型誤差判定

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

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

其中，實際階躍開始時刻tr0為穩定初值f(0)與實際輸出曲線的最后一個交點。

2.5 備選模型最優參數優化方法

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

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

3 實驗數據分析處理

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

3.1 最佳傳感器動態模型的確定

為了評價各階模型的適應性，對同一組傳感器輸出的數據分別使用各階模型進行優化，各階模型優化曲線與實際輸出曲線進行對比如圖2。評價模型優劣指標有三個，1)優化后模型的幅值誤差Ev;2)優化后模型的階躍時刻誤差Et;3)模型參數數量N。Ev表征了實際輸出曲線與優化曲線在取樣區間內的重合度;Et表征了實際曲線與優化曲線階躍開始時刻的差值，由于實際非理想階躍輸入的影響，該差值是必然存在的;模型參數數量N 決定了模型的復雜度，求取參數所需的時間與應用參數的難度。各階系統適應性評價指標如表1 所示。

圖2 各階系統優化曲線與實際輸出曲線對比Fig 2 Comparison of optimized curve and real output curve of each order system

表1 使用不同傳感器模型的優化結果誤差Tab 1 Optimized result error of different kind of models

由以上數據分析可知，幅值誤差Ev和階躍時刻誤差Et都隨模型的階數增大而減小，但模型階數越高確定模型所需的參數N 也會越多;二，三階模型的幅值誤差Ev和階躍時刻誤差Et遠小于一階系統;三階系統模型的幅值誤差Ev略小于二階系統，兩系統的幅值誤差Ev都足夠小。求解復雜模型參數耗費時間多，同時復雜模型也會給使用帶來不便，因此，本文將二階系統選定為SBE18 的動態響應模型。

3.2 環境因素對模型參數的影響

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

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

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

4 結 論

SBE18 的參數指標里僅標明了特定流速下時間常數，沒有標明該值的具體測量環境，而且時間常數也無法描述傳感器的具體動態過程，這給不少相關海洋科技工作者帶來很大的不便。本文深入研究了SBE18 動態參數的測量方法，分別分析了一，二，三階系統模型對傳感器動態過程的適應性，得出二階系統模型與SBE18 動態特性符合較好;同時，求取不同流速和pH 梯度變化方下傳感器二階系統模型參數，并分析了兩者對模型參數的影響。通過大量的實驗數據分析得出，模型參數隨流速增大而減小，且流速大于1m/s 后傳感器動態特性參數趨于穩定;pH 梯度變化方向對模型參數影響大，pH 由高到低變化時傳感器動態特性參數明顯大于pH 由低到高時的變化。

［1］ Wootton J T，Pfister C A，Forester 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，Feely 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.