熊楓 彭志敏 王振 丁艷軍 呂俊復 杜艷君
1) (華北電力大學控制與計算機工程學院,北京 102206)
2) (清華大學能源與動力工程系,電力系統與發電設備控制與仿真國家重點實驗室,北京 100084)
3) (清華大學能源與動力工程系,熱科學與動力工程教育部重點實驗室,北京 100084)
H2S 作為一種有毒且腐蝕性較強的氣體污染物,實現其濃度的準確測量意義重大.實際工業過程中,H2S 測量常受其他排放產物的干擾,本文基于腔衰蕩吸收光譜技術(CRDS)通過掃描6336—6339 cm—1 范圍內的吸收光譜,實現了H2S/CO2/CO 三組分物質濃度的同步測量,為實際工業過程中物質干擾下的H2S 濃度測量提供新思路.首先,對不同采樣長度下提取衰蕩時間的準確性進行分析,發現衰蕩信號的采樣長度約為衰蕩時間的8 倍時,衰蕩時間提取效果最好;通過不同壓力對比實驗確定最佳實驗壓力工況為50 kPa,并將最佳采樣長度與壓力工況應用于H2S 濃度測量.隨后,改變H2S 濃度對CO2/CO 干擾下系統對痕量H2S 濃度的測量效果進行檢驗,并對不同稀釋比例下濃度測量結果的線性度進行分析.最后,對本文CRDS 系統的檢測限進行分析,通過對4 組低濃度H2S 光譜的信噪比進行分析,得到H2S 的檢出限為6.9 ppb (1 ppb=10—9);通過對系統長期測量結果進行Allan 方差分析,得到系統對H2S 物質濃度的檢測下限約在2 ppb 左右.
硫化氫(H2S)是一種毒性較強的氣體污染物,濃度為25—50 ppm (1 ppm=10—6)的H2S 會對呼吸道產生強烈刺激,引發結膜炎,濃度>1000 ppm 時就會使人瞬間猝死[1-3].此外,H2S 還具有強腐蝕性,會對生產設備造成巨大損害[4-6].以火電行業為例,隨著低氮燃燒技術與循環流化床的廣泛應用,爐膛內部H2S,CO 等還原性氣體濃度大幅增大,這不僅會加劇鍋爐結焦結渣的問題;其中的H2S 還會對水冷壁造成強烈的高溫腐蝕,影響機組安全運行,嚴重時會引起水冷壁“爆管”,造成重大安全事故[7,8],故實現H2S 氣體濃度的準確測量具有重大意義.
然而,工業過程中H2S 實際排放濃度低、波動大且吸附性較強,實際測量還面臨許多問題.氣相色譜法[9-11]與電化學法[12]都可以實現H2S 氣體濃度的測量,但上述方法均具有無法回避的局限性.氣相色譜技術測量精度較高,但樣品需要進行預處理,導致單次測量耗時較長,因此該技術僅適用于實驗室條件下的氣體濃度測量.電化學法的測量系統簡單,但在實際應用過程中易受到SO2等組分的干擾,測量精度較低且使用壽命較短,難以滿足實際工業過程低濃度H2S 的診斷需求.
吸收光譜技術作為一種非接觸式的測量手段,具有系統簡單,響應速度快等優勢,近年來,基于該技術的H2S 測量也得到了較大的發展.Mathieu等[13]利用激波管對H2S-O2的反應機理進行了研究與分析.張楊等[14]對H2S 吸收信號的二次諧波進行測量.何岸等[15]將波長調制法與Herriott 池結合,設計開發了一種便攜式H2S 檢測儀,并在六盤水大灣煤礦進行了實際測量.Guo 等[16]利用Herriott 池對H2S 吸收光譜的二次諧波信號進行測量,實現了天然氣中H2S 氣體濃度的測量.王振等[17]將腔衰蕩吸收光譜技術與Herriott 池連用,建立了一種針對單一組分、寬量程、免標定的H2S 氣體濃度檢測方法.彭志敏等[18]選擇不同譜線,利用Herriott 池對低濃度CO 與H2S 的體積分數進行了動態測量.上述工作均針對H2S 濃度測量的相關問題進行了研究,但也存在相應局限性.首先,針對波長調制光譜而言,實際二次諧波幅值與濃度的關系需要通過繁瑣的標定過程得到,在實際應用中存在一定的局限性;更重要的是,上述工作大多關注H2S 的單組分測量,然而H2S 的吸收光譜測量常與其他排放產物,如CO2,CO,H2O 等的吸收譜線交疊,如何排除其他物質吸收譜線的干擾或實現多組分的同步測量有待進一步研究.
針對目前H2S 濃度測量研究的局限性,本文通過腔衰蕩吸收光譜技術實現了H2S/CO2/CO 三組分物質濃度的同步測量,為實際工業過程中,在其他排放產物干擾下的H2S 濃度測量提供了新思路.本文首先針對不同采樣長度下衰蕩時間提取的準確程度進行分析,利用修正的連續積分法(CSI 法)分別對實驗與仿真信號的衰蕩時間進行提取,對所提取衰蕩時間的平均值(準度),標準差(精度)與消耗時間(速度)進行分析,發現當采樣長度為衰蕩時間的8 倍時,提取效果最好并作為實驗信號的采樣長度.選擇6336—6339 cm—1范圍對三組分的吸收光譜測量,通過不同壓力下的探究實驗確定50 kPa 為最佳壓力工況并用于后續實驗.通過實驗探究CO2/CO 干擾下系統對不同濃度H2S 的測量效果;進一步對不同稀釋比例下三組分濃度測量結果的線性度進行比較,并對誤差產生原因進行分析.最后,對本文CRDS 系統的檢測限進行分析,通過對四組低濃度H2S 光譜的信噪比進行分析,得到H2S 的檢出限約為6.9 ppb;通過對系統長期運行結果進行Allan 方差分析,得到本文系統對H2S 的檢測下限約在2 ppb 左右.
腔衰蕩吸收光譜技術利用兩片反射率極高的鏡片(R>99.99%)實現激光在腔內的往復反射,極大提升有效光程.光束在腔內往復反射過程中,各次出射信號隨時間指數衰減,見圖1.往返n次后出射光強In為:In=I0·R2n·exp(—2nδd).光線一次往復用時t=2L/c,則第n次出射光線在腔內傳播時間為t=2nL/c,時刻t出射光線強度I(t)為:In=I0·exp[—t·c·(—lnR+δd)/L].由于鏡片反射率R→1,有lnR≈ — (1—R),則時刻t出射的光線強度I(t)可以表達為

圖1 腔衰蕩光譜基本原理Fig.1.Basic schematic of cavity ring-down spectroscopy.
其中,I0為入射光線強度,c(cm/s)為光速,L(cm)為光程,即兩球面反射鏡中心距離,R為鏡片反射率,δ(cm—1)為腔內總損耗系數,d(cm)為腔內介質有效吸收長度光程(d≤L),常認為與光程近似相等,即d≈L.
出射信號I(t)的衰蕩時間與腔內氣體的吸收k(v)有關.測量無吸收時的空腔衰蕩時間τ0(v)和存在氣體吸收時的衰蕩時間τ(v),通過兩者差值可計算得到吸收系數k(v)[19-22]:
其中,φ(v) (cm)為線型函數,由譜線位置v0(cm—1)與真實展寬Δνline(cm—1)決定,包括碰撞展寬、多普勒效應等;P(atm)是絕對壓力;X為體積濃度;S(cm—2·atm—1)為單位壓力線強度;在較窄的波長范圍內,τ0(v)可認為是常數τ0.若測量工況(P,T)和譜線參數(S)已知,則可通過衰蕩時間τ(v)確定待測氣體濃度X.
本文所用腔衰蕩吸收光譜測量系統分為兩個部分: CRDS 光譜測量系統[23-25]與氣體配置系統,如圖2 所示.CRDS 部分,激光控制器驅動DFB激光器(四川騰光)輸出激光,經光纖隔離器與聲光調制器后進入準直器,準直后出射光線進入衰蕩腔.衰蕩腔由一對鍍有多層介質膜的高反鏡(R>99.99%,1570—1655 nm,ATF)組成,鏡面中心間距約50 cm,腔內鍍有硅烷鈍化層以降低待測氣體的吸附與腐蝕作用.腔體一側安裝有環形壓電陶瓷用以驅動鏡片往復運動,實現腔長掃描,進而完成腔內縱模與輸入激光模式的匹配.由于激光持續輸入,此時并無衰蕩事件產生.輸出光線經探測器轉化為電信號,一路輸入采集卡,另一路輸入脈沖信號發生器.當縱模信號強度超過設定值時,脈沖發生器就會輸出方波信號;方波信號經射頻發生器轉化后輸入并觸發聲光調制器.此時,聲光調制器被觸發,關斷激光輸出,失去激光輸入后,衰蕩腔出射光強隨時間指數衰減,產生衰蕩事件.衰蕩信號經采集卡后輸入計算機,通過CSI 法提取衰蕩時間.一次衰蕩事件完成后,計算機驅動激光控制器步進電流,實現波長掃描.為保證光譜信噪比,本文實驗光譜為100 次測量的平均結果.

圖2 腔衰蕩光譜測量系統(PC: 計算機,LC: 激光器控制器,DFB: 激光器,ISO: 光纖隔離器,AOM: 聲光調制器,RDC: 衰蕩腔,PZT: 壓電陶瓷,APD: 雪崩式光電探測器,PG: 脈沖信號發生器,RF: 射頻發生器,DAQ: 數據采集系統,RD: 衰蕩信號,Trig: 觸發信號)Fig.2.Cavity ring-down spectroscopy measurement system (PC: personal computer,LC: laser controller,DFB: DFB laser,ISO:fiber isolator,AOM: acousto-optic modulator,RDC: ring-down cavity,PZT: piezoceramics,APD: avalanche photodiode,PG: pulse generator,RF: radio frequency generator,DAQ: digital acquisition,RD: ring-down signal,Trig: trigger).
配置不同濃度待測氣體時,利用混合氣罐對H2S,CO2,CO 標準氣體與純N2進行混合,通過調整四種組分的分壓力來配置組分濃度不同的待測氣體.配氣后靜置六個小時以確保各組分混合均勻.實驗時,待測氣體壓力由氣室連接的壓力傳感器測量得到,氣體溫度近似為室溫,由水銀壓力計測量得到.
衰蕩時間的提取是腔衰蕩吸收光譜技術的關鍵,直接關系到CRDS 測量的速度與準確性.衰蕩時間是通過對實際測量的指數信號進行數據處理而得到,本文使用修正的連續積分法(CSI 法)來提取衰蕩時間.Matheson[26]將連續積分(SI)用于對實際信號的快速處理,而未修正的連續積分法針對實際離散數據進行處理時存在系統誤差,Halmer 等[27]在其研究過程中發現該問題并對連續積分進行修正,提出了CSI 方法,廣泛應用于指數信號的快速擬合,其具體計算方法見附錄.
實際測量的衰蕩信號可以表示為:It(t)=A·exp(-t/τ)+B+ε,其中,A為信號幅值,B為時間無關的信號偏移基線,ε為噪聲.CSI 法為三參數方法,無需扣除基線B,對實際衰蕩信號進行計算可直接得到A,B與τ.為方便討論,定義無量綱常數Tm=(tN-t0)/τ為衰蕩信號的采樣長度,其中,tN表示衰蕩信號結束時刻,t0表示開始時刻,Tm則表示采樣時長與信號衰蕩時間的倍數關系.衰蕩信號與采樣長度Tm如圖3 所示.

圖3 不同衰蕩信號采樣長度Fig.3.Different sampling length of ring-down signal.
不同Tm下,衰蕩時間τ的提取效果并不相同,為確定最佳采樣長度,對仿真與實驗兩類信號進行分析.實驗采集100 組衰蕩信號原始數據,工況為純N2,50 kPa、室溫(~26 ℃),衰蕩時間τ約160 μs(真實值為159.8 μs),噪聲標準差約1%,仿真信號參考實驗信號的測量效果并添加1%白噪聲.分別對100 組實驗與仿真信號進行處理,將τ提取結果的平均值、標準差στ與計算消耗時間Tcal作為準度、精度與速度的評價標準,不同Tm下的τ提取效果如圖4 所示,紅(藍)實線表示對仿真信號提取結果,黑色實心點表示對實驗信號處理結果.

圖4 不同采樣長度下,衰蕩時間提取效果Fig.4.The result of ring-down time extraction in different sampling length.
從仿真信號的提取效果可以看出,Tm過短,信號包含數據點過少,衰蕩時間τ的標準差στ較大,說明τ提取結果的隨機波動較大,提取精度較差;同時由于采樣區間為信號前段,數據衰減快速,平均值小于真實值,τ提取準度較差.隨著數據長度即Tm增大,用于提取τ的數據量增加,στ持續減小,逐漸接近真實值,精度與準度逐漸提高;當Tm> 5·τ后,仿真信號的στ與趨于穩定值.由于數據處理的速度一定,計算耗時Tcal隨著Tm近似正比例增大.從仿真信號的τ提取效果可以看出,當Tm過短時(Tm< 3),τ提取精度與準度均較差,隨著Tm增大,τ提取效果逐步提升;當Tm足夠長時(Tm> 5),Tm增大對精度與準度的改善作用較小,而計算耗時Tcal隨Tm持續增大.
實驗信號衰蕩時間τ的提取效果與仿真信號的趨勢基本一致.對比仿真信號,實驗信號的平均值在Tm較小時的偏差更大且存在波動,需要更大的Tm來使得趨于穩定.由于實驗中的非理想因素與其他噪聲,實驗信號的標準差στ大于對應Tm處仿真信號的στ.此外,當Tm過大時,衰蕩信號尾端信號已經極弱,實驗信號中噪聲的影響已經大于實際信號對τ提取的正向作用,實驗信號的標準差στ反而隨Tm的增大而有所升高,精度下降.針對實驗信號進行分析,當Tm過短時(Tm< 8),τ偏差較大,準度不足;而當Tm過長時(Tm> 10),στ開始上升,精度下降,數據處理所耗時間Tcal較長.綜合來看,當Tm≈ 8τ時,趨近于穩定值,στ也近似最小值,兼顧準度與精度,并且所耗時間Tcal較短,τ提取效果最好.因此,在后續同步測量實驗中均選取8 倍衰蕩時間作為采樣長度.
考慮實際工業排放中H2S 濃度極低且易受到其他排放產物的影響,如CO2、CO,選擇6335—6339 cm—1范圍作為目標區域,實現H2S/CO2/CO三組分物質濃度的同步測量,待測區域譜線與線強度分布如圖5 所示.該范圍內H2S 的線強度比CO2/CO 大一個數量級,可以實現在CO2/CO 干擾下,較低濃度H2S 的測量.此外,3 種目標物質的主要吸收譜線位置分立,可以較清晰地分辨三者的吸收譜線.但該范圍包含譜線較多,無明顯的未吸收區,需要在實驗前對空腔衰蕩時間τ0進行測量.

圖5 待測區域譜線展示與不同壓力測量結果Fig.5.Display of spectral lines in measured region and the measured spectra in different pressures.
定義吸收系數峰值與殘差標準差的比值為信噪比SNR,為確定信噪比最高的實驗工況,對不同壓力下三組分的吸收光譜進行測量.H2S/CO2/CO濃度配置為250/250/20 ppm,測量吸收系數k(v)光譜與對應譜線強度S如圖5(a)和(b)所示.譜線參數來自Hitran 數據庫,使用Voigt 線型對吸收光譜進行擬合.可以看出,不同壓力下測量效果不同,隨著壓力的降低,光譜信噪比逐漸提升.這與吸收光譜的一般規律相吻合,由于Voigt 線型對實際吸收譜線進行擬合時存在結構性殘差,壓力越大,譜線展寬越大,結構性殘差覆蓋范圍越廣,殘差標準差越大,信噪比越低.但25 kPa 下測量信噪比最低,主要是由激光器電流調制系數過大(約0.17 cm—1/mA)導致的.控制器輸出電流的精度一定,激光器電流調制系數過大導致輸出激光的頻率出現較大波動,導致激光匹配到諧振腔目標縱模的相鄰模式,腔與激光的模式匹配出現問題.壓力較大(不小于50 kPa)時,譜線展寬較大,吸收譜線兩翼處斜率相對較小,頻率不穩定性所帶來的影響相對較小;而25 kPa 時譜線展寬已經較小,吸收峰兩翼處斜率最大,輸出激光頻率的不穩定性對其影響最為嚴重.通過對比不同壓力下光譜信噪比,發現測量系統在壓力為50 kPa 時信噪比最高,選擇該壓力為后續同步測量的壓力工況.
該波長范圍還包含H2O 的吸收譜線,但其線強度較弱,約比3 種目標物質的譜線強度低3 個數量級左右.不同壓力實驗未考慮H2O 吸收的影響,從殘差中可以看出,雖然該范圍內水吸收極弱,但仍然會對測量結果產生影響,綠框內的單峰形殘差就是水的吸收.由于H2O 濃度并不足夠高且線強度極低,H2O 吸收極弱,濃度測量的誤差極大,故無法實現H2S,CO2,CO 與H2O 的四組分同步測量,只能將水的吸收作為干擾項而消除.后續不同組分濃度實驗均考慮H2O 吸收的影響對測量結果進行修正.
對混有CO2/CO 干擾組分的不同濃度H2S 氣體進行測量.測量壓力為50 kPa,固定CO2/CO濃度為250 ppm,探究不同濃度H2S 的測量效果,設置H2S 濃度為20 ppm,10 ppm,5 ppm,測量結果如圖6 所示,相對誤差小于2.5%,偏差較小.此外,殘差具有“w”形結構,證明實際測量結果可靠,這是由于Voigt 線型為Gauss 線型與Lorentz 線型的簡單卷積,使用Voigt 線型來近似實際吸收譜線會有結構性殘差;實際吸收還存在Dicke 收斂、速度依賴等物理現象,使得實際譜線線型比Voigt 線型峰值增高,兩翼收窄,導致“w”形殘差出現[28-31].該組實驗SNR 均在大于200,表明實驗測量及光譜擬合具有較高精度,滿足濃度測量要求.3 組實驗H2S 濃度的測量結果分別為: 18.2 ppm,9.4 ppm,4.8 ppm,而CO2/ CO 測量結果為: 230/236 ppm,236/241 ppm,240/243 ppm.實驗表明,采用本文所選譜線及實驗工況,可在H2S 濃度低于干擾氣體CO2/CO 近2 個數量級的環境下,實現三組分同步測量.

圖6 僅改變H2S 濃度測量結果Fig.6.the spectra in only H2S concentration changed.
通過實驗結果可以看出,隨著H2S 濃度的減小,組分濃度測量的偏差減小,H2S 光譜殘差對測量光譜總體殘差貢獻最大.分析原因為H2S 譜線參數的不確定度較大.Hitran 數據庫[32]中,H2S 線強度誤差等級為4 (10% < 偏差 < 20%)、5 (5% <偏差 < 10%),而CO2/CO 線強度誤差等級為6(2% < 偏差 < 5%)、7 (1% < 偏差 < 2%),H2S 譜線參數不確定度最大.
考慮到所用標氣的不確定度與配氣過程的非理想因素,氣體各組分的實際濃度并不完全等于配氣濃度的理論值.為檢驗三組分濃度測量的準確性,對不同濃度下三組分濃度測量結果與理論值之間的線性度進行分析.保持三組分物質濃度比例不變,H2S/CO2/CO 濃度在250/250/20 ppm 與25/25/2 ppm 之間設置4 個濃度梯度.不同濃度下的光譜測量結果與線性度如圖7 所示,實驗結果的光譜信噪比在185—220 范圍內,滿足濃度測量要求.對不同稀釋比例下的濃度測量結果進行線性度分析,可以看出三組分濃度關系的線性回歸系數R2均大于0.999,線性度較好,證明在該范圍內測量系統對三組分濃度變化的響應具有較好的線性度.此外,系統對不同組分濃度變化的響應效果并不完全相同,其中,H2S 的響應曲線的R2與一次項系數均最小,濃度測量線性度較差,測量值偏小的系統誤差最明顯;CO2/CO 的測量效果相對較好,通過R2與一次項系數可以看出,CO2測量的線性稍好,而CO 測量濃度偏小的系統誤差更小.

圖7 變物質濃度測量光譜與線性度分析Fig.7.The spectra of changing concentration and analysis of linearity.
實驗中測量濃度均小于氣體配置時的理論濃度,針對系統誤差進行分析,認為可能存在如下3 個方面問題: 1)物質的吸附.H2S 本身的吸附效應最為明顯,待測組分的吸附導致混合氣體實際濃度略低于理論配氣濃度值.2)譜線參數存在誤差.Hitran 數據庫譜線本身不確定度會導致濃度測量的誤差.3)空腔衰蕩信號與吸收信號分時測量,測量空腔衰蕩時間與存在待測氣體但無吸收時的基底衰蕩時間并不完全相同,空腔衰蕩時間的偏差導致測量結果系統偏差.
為進一步分析本文CRDS 系統對H2S 的檢測限,對低濃度H2S 的光譜信噪比與測量系統的Allan 方差進行研究.光譜信噪比反映波長掃描范圍內吸收光譜測量的整體效果.對濃度為2 ppm,1 ppm,500 ppb,250 ppb 的H2S 進行測量,結果如圖8 所示.由于H2S 濃度極低,水的吸收較為明顯,在對實驗數據進行處理時,同時考慮水吸收帶來的影響.4 組實驗測量的信噪比分別為: 98,83,49,36,若定義信噪比為1 時,對應濃度為物質的檢出限,則4 組實驗對應的H2S 檢出限分別為:20.5 ppb,12.1 ppb,10.2 ppb,6.9 ppb.高濃度測量結果對應的檢出限更高,這是因為濃度越大,吸收越強,voigt 線型引入的“w”形結構殘差越大,使得信噪比并不隨濃度等比例增大.選擇“w”形殘差影響最小、濃度為250 ppb 的測量結果作為評價標準,H2S 的檢出限為6.9 ppb.

圖8 低濃度H2S 測量結果Fig.8.The spectra of H2S in low concentration.
最后,再通過Allan 方差[33]對系統的檢測下限進行分析.通過對固定波長處(H2S 吸收峰值附近)的衰蕩時間進行長期測量并對1/c·τ做Allan方差分析,可以得到測量系統吸收系數k(v)的檢測下限.在室溫(~26 ℃),50 kPa 壓力,純N2的條件下,對測量系統的空腔衰蕩時間進行測量,測量結果與Allan 方差如圖9 所示.系統測量時長約15 min,采樣時間間隔約為45 ms.可以看出,對1200 點進行平均,即積分時間約為54 s 時,系統可達到最低的吸收系數檢測下限,約3.5×10—11cm—1.結合線強度信息將吸收系數轉化為對應H2S 濃度,約為2 ppb.Allan 方差的分析結果與前由信噪比得到的H2S 的檢出限基本一致.進一步增加平均次數,測量效果反而下降,Allan 方差曲線上升,表明系統長期測量的衰蕩時間存在漂移,這是由于激光頻率飄移,以及腔長、溫度、壓力等因素變化所導致.可以看出,CRDS 具有極低的檢測下限,但由于其對光在腔內的諧振條件與腔的穩定性要求較高,長期運行會出現衰蕩時間的系統性漂移.

圖9 吸收系數Allan 方差Fig.9.Allan variance of absorption coefficient.
本文利用腔衰蕩吸收光譜技術實現了CO2/CO干擾下,H2S 物質濃度的準確測量.首先針對不同采樣長度下衰蕩時間提取的準確程度進行分析,利用CSI 算法對實驗與仿真信號的衰蕩時間進行提取,對所提取衰蕩時間的平均值(準度),標準差(精度)與消耗時間(速度)進行分析,發現采樣長度為衰蕩時間的8 倍時,提取效果最好并作為同步測量實驗信號的采樣長度.選擇6336—6339 cm—1范圍對三組分的吸收光譜測量,通過不同壓力實驗,發現50 kPa 下信噪比最高并作為后續實驗的壓力工況.通過僅改變H2S 濃度的實驗來對CO2/CO干擾下系統對痕量H2S 濃度的測量效果進行檢驗,并對不同稀釋比例下濃度測量結果的線性度進行分析.通過對濃度測量的誤差分析,發現物質吸附性,譜線參數的不確定度與空腔衰蕩時間分時測量為導致系統誤差主要因素.最后,對本文CRDS 系統的檢測限進行分析,通過對4 組低濃度H2S 光譜的信噪比進行分析,得到H2S 的檢出限約為6.9 ppb;通過對系統長期運行結果進行Allan 方差分析,得到本文系統對H2S 的檢測下限約在2 ppb左右.
附錄
CSI 法通過求解如下形式矩陣獲得系數a1,進而計算衰蕩時間τ:
矩陣中的各元素可由實際測量信號It(t)計算得到.It(t)為包含N個數據點,等時間間隔 Δt的離散數據,在ti時刻的數據值為It(ti),i=1,2,···,N,各元素具體計算公式如下: