?

基于正態變換的變樣本容量Shewhart比例控制圖設計

2024-03-13 13:08周思陽柯楚賢
計算機集成制造系統 2024年2期
關鍵詞:樣本容量警戒警報

周思陽,張 瑩 ,柯楚賢

(1.交通運輸部水運科學研究所,北京 100088;2.武漢理工大學 交通與物流工程學院,湖北 武漢 430063)

0 引言

控制圖作為一種有效的統計過程控制(Statistical Process Control, SPC)工具,被廣泛應用于科學研究和工業生產。在一些生產場景中,成分間的比例關系是一種關鍵質量特性,需要對其進行監控以保證產品質量的穩定性,這些場景大致分為兩類:①生產過程中混合物的成分之間存在特定的比例關系,如食品生產中原料的配比[1]和港口配煤作業中原煤的配比[2];②生產過程中涉及物理或化學反應等特定操作前后產品成分的比例,如在電池回收過程中可回收部分的重量比例[3]。因此近年來,如何利用控制圖監控兩個變量之間的比例關系逐漸成為研究熱點。

在比例控制圖研究中,CELANO等[1]首先在樣本容量n=1的情況下對兩個正態變量之比進行監控,提出Shewhart-RZ 控制圖,并對其統計性能進行了研究;隨后,這項研究由CELANO等[4]擴展到樣本容量n>1的情形;NGUYEN等[5]提出單邊Shewhart-RZ控制圖,并利用線性協變量模型研究了測量誤差對其性能的影響;ZAIDI等[6]和TRAN等[7]將對兩個變量間比例的監控擴展到對成分數據的監控,從而可以同時對混合物中的多種成分進行監控。當比例特性偏移較大時,Shewhart-RZ控制圖能較快發現過程失控,但由于其僅能利用當前樣本信息,當比例特性的偏移較小時,往往需要較長時間才能發現過程失控。研究中通常采用附加運行規則、與合格品鏈長(Conforming Run Length,CRL)控制圖相結合、變抽樣區間和變樣本容量等策略中的一種或多種對Shewhart控制圖進行優化,以改善其在偏移較小時的性能。其中針對Shewhart-RZ控制圖,CELANO等[8]和NGUYEN等[9]分別將雙邊和單邊Shewhart-RZ控制圖與合格品鏈長控制圖結合;TRAN等[10]和TRAN[11]采用不同附加運行規則對雙邊Shewhart-RZ控制圖進行優化;NGUYEN等[12]采用變抽樣區間策略對單邊Shewhart-RZ控制圖進行優化設計。目前采用變樣本容量策略對單邊Shewhart-RZ控制圖進行優化的研究比較缺乏,有些研究采用兩種記憶性控制圖,即指數移動加權平均(Exponentially Weighted Moving Average, EWMA)和累計和(CUmulative SUM, CUSUM)控制圖,對二元正態變量間的比例關系進行監控[13-17],計算結果表明其性能優于Shewhart-RZ控制圖。

分析發現,動態控制圖中的變樣本容量策略能有效提高Shewhart型控制圖的性能,并具有易于實施的良好特性。然而,目前比例控制圖研究多集中于靜態控制圖領域,對動態控制圖的研究較少,且均未考慮變樣本容量策略在比例控制圖中的應用。本文提出一種單邊變樣本容量Shewhart比例控制圖,即VSS Shewhart-RZ*控制圖,來監控二元正態變量間的比例關系,并比較其與單邊Shewhart-RZ控制圖的性能,結果表明相同情況下VSS Shewhart-RZ*控制圖具有更小的ARL。

1 比例Z分布

假設W=(X,Y)T為二元正態隨機向量,W~N(μW,∑W),μW為均值向量,∑W為協方差矩陣,具體計算公式如下:

(1)

(2)

式中ρ為變量X,Y之間的相關性系數。此時變量X,Y的變異系數可定義為γX=σX/μX,γY=σY/μY,同時記X,Y的標準差之比為ω=σX/σY。

令X,Y的比值為Z,Z=X/Y,由于無法得到比例Z的累計分布函數的準確表達式,CELANO等[4]對其近似估計為

(3)

式中:Φ(·)為標準正態分布的累積分布函數;A和B為參數z,γX,γY,ω,ρ的函數,

(4)

(5)

比例Z的概率分布函數可近似估計為

(6)

式中φ(·)為標準正態分布的概率分布函數。

比例Z累積分布函數的反函數可由式(3)推導得到

(7)

式中C1,C2,C3是參數γX,γY,ω,ρ的函數:

(8)

(9)

(10)

為了在生產過程中監控比例Z,在時刻i(i=1,2,…)抽取樣本容量為n的樣本{Wi,1,Wi,2,…,Wi,n},Wi,j=(Xi,j,Yi,j)T~N(μW,i,∑W,i)(j=1,…,n)為二元正態隨機向量,其均值向量和協方差矩陣分別為

(11)

(12)

與CELANO等[4]的研究相同,對樣本數據做以下假設:①不同抽樣時刻下產品的規格可以發生變化,即當i≠k時,有μW,i≠μW,k,∑W,i≠∑W,k;②對于不同的樣本,變量X,Y的變異系數γX,γY為已知常數,即對樣本{Wi,1,Wi,2,…,Wi,n}可得σX,i=γX×μX,i,σY,i=γY×μY,i;③當過程處于受控狀態時,變量X,Y的比例為已知常數z0=μX,i/μY,i(i=1,2,…)。此時觀測統計量為

(13)

(14)

(15)

2 VSS Shewhart-RZ*控制圖

2.1 VSS Shewhart-RZ*控制圖的實施流程

由于比例Z的分布不對稱,本文設計了兩個單邊VSS Shewhart-RZ*控制圖分別監控比例向上或向下的偏移,以避免性能損失。在抽樣時樣本容量n有兩種可能的取值,即小樣本容量nS和大樣本容量nL,nS

(1)在進行控制圖優化設計時,需要計算n(i)=nS時的警戒限和控制限{(LWLS,LCLS),(UWLS,UCLS)},以及n(i)=nL時的警戒限和控制限{(LWLL,LCLL),(UWLL,UCLL)},因此需要花費更多的計算資源和時間對控制圖的參數組合進行尋優。

(2)在實際使用過程中,使用者需要檢查每個樣本點所處的位置,以判斷是否需要更換控制限和警戒限組合,增加了使用者的操作難度。

(16)

此時上單邊VSS Shewhart-RZ*控制圖的控制限和警戒限可固定為

UCL=z0×KU,

(17)

UWL=z0×WU。

(18)

與之類似,下單邊VSS Shewhart-RZ*控制圖的控制限和警戒限為

LCL=z0×KD,

(19)

LWL=z0×WD。

(20)

式中:KU,KD為控制限系數,WU,WD為警戒限系數,且滿足WU≤KU和KD≤WD。

上單邊VSS Shewhart-RZ*控制圖的區域劃分如圖1a所示,工作機制如下:

下單邊VSS Shewhart-RZ*控制圖的區域劃分如圖1b所示,工作機制與上單邊控制圖的工作機制類似,不再贅述。

2.2 VSS Shewhart-RZ*控制圖的性能指標

在變樣本容量控制圖研究中,常采用平均運行長度ARL和平均樣本容量ASS作為控制圖的性能指標。ARL指從過程開始運行到控制圖發出過程失控信號時所經過的平均樣本數,受控狀態下的ARL0越大越好,此時控制圖虛發警報的概率降低,可避免多余的停機檢查;失控狀態下的ARL1越小越好,此時控制圖漏發警報的概率降低,可以快速檢測到異常因素。ASS為樣本容量的均值,ASS=E(n(i)),當控制圖的性能達到需求時,ASS越小越好,更小的平均樣本數表示抽樣檢測成本更小。比較VSS型控制圖與固定樣本容量的靜態控制圖時,使受控狀態下的ARL0和ASS0保持一致,比較失控狀態下的ARL1和ASS1,更小的ARL1或ASS1表示擁有更好的性能。

假設過程失控時,過程中存在異常因素使受控狀態下的比例z0偏移至z1=τ×z0,其中τ>0,τ的取值表示偏移的大小和方向,τ>1對應比例向上的偏移,τ<1對應比例向下的偏移,過程處于受控狀態時τ=1。同時考慮失控狀態下,變量(X,Y)間的相關系數也可能發生偏移,即從受控狀態下的ρ=ρ0偏移至ρ=ρ1。

采用馬爾可夫鏈推導單邊VSS Shewhart-RZ*控制圖的性能指標。將單邊VSS Shewhart-RZ*控制圖建模為一個含有3種狀態的馬爾可夫鏈,其中前兩種為受控狀態,分別對應使用小樣本容量n(i)=nS的情形和使用大樣本容量n(i)=nL的情形,第3種為失控狀態,在馬爾可夫鏈模型中為吸收態。馬爾可夫鏈的一步轉移概率矩陣

(21)

式中:Q為2×2的轉移概率矩陣,1=(1,1)T,r為2×1的向量,r=1-Q1,即矩陣P中每一行的和為1;pS(n(i))和pL(n(i))分別為從當前狀態轉移至小樣本容量和大樣本容量的概率,n(i)∈{nS,nL}表示當前狀態。對于上單邊VSS Shewhart-RZ*控制圖,

(22)

(23)

對于下單邊VSS Shewhart-RZ*控制圖,

(24)

(25)

此時控制圖的ARL為

ARL=qT(I-Q)-11。

(26)

式中:I為2×2的單位矩陣,q為初始概率向量,1=(1,1)T。當設定第1組樣本為小樣本容量,即n(1)=nS時,q=(1,0)T;當設定第1組樣本為大樣本容量,即n(1)=nL時,q=(0,1)T。

計算ASS時需要先將轉移矩陣P轉換成P*,當n(1)=nS時,

(27)

當n(1)=nL時

(28)

與矩陣P不同,P*沒有吸收態,當馬爾可夫鏈到達失控狀態時,消除過程中的異常因素后立刻重啟,使馬爾可夫鏈恢復至初始狀態。由P*定義的馬爾可夫鏈有穩態概率向量π=(πS,πL,πOOC)T,可由下式求得:

(29)

其中穩態概率(πS,πL,πOOC)與樣本容量(nS,nL,n(1))相聯系。當n(1)=nS時,矩陣R可由矩陣P*先轉置,再將對角線元素減1,并將第1行元素替換為1,得到

(30)

當n(1)=nL時,矩陣R可由矩陣P*先轉置,對角線元素減1后,將第2行元素替換為1,得到

(31)

此時ASS可由下式計算:

ASS=(nS,nL,n(1))π。

(32)

2.3 VSS Shewhart-RZ*控制圖的優化設計

VSS Shewhart-RZ*控制圖的優化目標是在特定的比例偏移τ和相關系數ρ1下,確定最優的控制限系數KU(或KD)、警戒限系數WU(或WD)和樣本容量組合(nS,nL),使受控狀態下控制圖的平均運行長度ARL0=L0,平均樣本容量ASS0=n0,同時使失控狀態下的ARL1最小。其中:L0與受控狀態下控制圖誤發警報的概率α密切相關,L0=1/α;失控狀態下的ARL1與漏發警報的概率β相關,ARL1=1/(1-β)??刂茍D的優化過程即在給定可接受的誤發警報概率α下,使失控狀態下控制圖漏發警報的概率β最小。由于在實際應用中樣本容量的最大值一般會有限制,與之前的研究[21-23]類似,本文設置樣本容量可取的最大值為31,則nS,nL需要滿足1≤nS

minARL1(nS,nL,WU,KU,γX,γY,ρ1,τ,z0)。

s.t.

ARL0(nS,nL,WU,KU,γX,γY,ρ0,τ=1,z0)=L0;

ASS0(nS,nL,WU,KU,γX,γY,ρ0,τ=1,z0)=n0;

WU≤KU;

1≤nS

nS,nL∈N+。

(33)

下單邊VSS Shewhart-RZ*控制圖的數學模型與上單邊類似,僅警戒限和控制限之間大小關系的約束略有不同,需將其改寫為KD≤WD。其中L0的取值一般需要事先給定,基于3σ原理常取L0=370.4,本文為了方便比較,與之前研究[4]保持一致,取ARL0=200。n0的取值為實際過程中受控狀態下樣本容量的期望值,為方便比較,本文取n0∈{5,15}。

由于非凸的混合整數非線性規劃問題屬于NP-hard問題,求解非常困難,目前缺少有效的求解算法,本文將整數變量分離,使問題化簡為對多個非線性方程組來求解,具體過程為(如圖2):

(1)確定γX,γY,ρ0,ρ1,τ,z0,n0,L0等參數的取值。

(2)令nS=1。

(3)令nL=n0+1。

(4)在當前參數下,采用MATLAB軟件中的fslove函數求解由兩個等式約束構成的非線性方程組,得到一組警戒限和控制限變量的取值(WU,KU)或(WD,KD)。

(5)根據當前(nS,nL,WU,KU)或(nS,nL,WD,KD)的取值計算ARL1指標。

(6)保持nS不變,令nL=nL+1。

(7)重復(4)~(6),直至nL=31。

(8)令nS=nS+1。

(9)重復(3)~(8),直至nS=n0-1。

3 數值分析及性能比較

在不同參數設置下求得單邊VSS Shewhart-RZ*控制圖的最優決策變量組合,計算ARL1性能指標,并與固定樣本容量(Fixed Sample Size,FSS)的單邊Shewhart-RZ控制圖[5]、單邊Synthetic-RZ控制圖[9]和單邊VSI Shewhart-RZ控制圖[12]進行比較。同時由2.2節可知,n(1)取值會對控制圖性能產生影響,目前變樣本容量控制圖文獻中大多只研究了n(1)=nS的情況,對n(1)=nL的情形研究較少,因此將在不同n(1)取值下研究單邊VSS Shewhart-RZ*控制圖的性能。為方便對比,采用與文獻[4]的參數組合,如n0∈{5,15},γX∈{0.01,0.2},γY∈{0.01,0.2},ρ∈{0.0,±0.4,±0.8},τ∈{0.9,0.95,0.98,0.99,1.01,1.02,1.05,1.1},考慮γX=γY,γX≠γY,ρ0=ρ1和ρ0≠ρ1等場景。

圖3和圖4所示為單邊Shewhart-RZ控制圖和單邊VSS Shewhart-RZ*控制圖(包括n(1)=nS和n(1)=nL兩種情形,沒有特別指明時默認n(1)=nS),在過程失控且X,Y的相關系數保持不變(ρ0=ρ1)時,ARL1的值隨比例偏移變化的情況。圖3主要展示變量X,Y的變異系數相等(γX=γY)時的情形,其中左邊兩列圖為γX=γY=0.01的情形,右邊兩列圖為γX=γY=0.2的情形。圖4主要展示γX≠γY時的情形,其中左邊兩列圖為(γX=γY)=(0.01,0.2)的情形,右邊兩列圖為(γX=γY)=(0.2,0.01)的情形。單邊VSS Shewhart-RZ*控制圖與單邊Synthetic-RZ控制圖、單邊VSI Shewahrt-RZ控制圖的性能比較在圖3中進行了展示,其中單邊Synthetic-RZ控制圖的性能指標為ARL,單邊VSI Shewahrt-RZ控制圖的性能指標為平均報警時間ATS。從圖3和圖4可得以下結論:

(1)單邊VSS Shewhart-RZ*控制圖性能受(γX,γY)和ρ影響。變量(X,Y)的變異系數(γX,γY)取值越小,控制圖的性能越好。例如當ρ=-0.8,n0=5,τ=0.9時,若(γX,γY)=(0.01,0.01)則ARL1=1,若(γX,γY)=(0.2,0.2)則ARL1=19.8。ρ>0時控制圖的性能優于ρ<0時的性能,即變量(X,Y)之間呈正相關時控制圖的性能更好。例如當(γX,γY)=(0.2,0.2),n0=5,τ=0.9時,若ρ=-0.8則ARL1=19.8,若ρ=0.8則ARL1=2.6。

(2)對于相同百分比的偏移ΔZ=100×|τ-1|,當γX=γY時單邊VSS Shewhart-RZ*控制圖對下偏移的檢測性能略好于對上偏移的檢測性能。例如當(γX,γY)=(0.2,0.2),n0=5,ρ=-0.8,τ=0.99(或τ=1.01),即ΔZ=1%時,可得下單邊控制圖的ARL1=170.3(上單邊控制圖ARL1=170.6),如圖3所示。當γX≠γY時,單邊VSS Shewhart-RZ*控制圖對上下偏移的性能趨勢取決于變異系數(γX,γY)之間的大小關系。當γX<γY時,控制圖對下偏移的檢測性能更好;當γX>γY時,控制圖對上偏移的檢測性能更好。例如當(γX,γY)=(0.01,0.2),n0=5,ρ=-0.8,τ=0.99(或τ=1.01)時,ARL1=137.8(或ARL1=157.4);當(γX,γY)=(0.2,0.01),n0=5,ρ=-0.8,τ=0.99(或τ=1.01)時,ARL1=157(或ARL1=138.3),如圖4所示。

圖5所示為單邊Shewhart-RZ控制圖和單邊VSS Shewhart-RZ*控制圖,當過程失控X,Y的相關系數發生變化(ρ0≠ρ1,n0=5)時,ARL1隨比例偏移變化的情況。受控狀態下ρ0=±0.4,受異常因素影響偏移至ρ1=0.5×ρ0或ρ1=2×ρ0。從圖5可以觀察到以下趨勢:增強負相關性會提高單邊VSS Shewhart-RZ*控制圖對比例偏移的敏感性,例如(γX,γY)=(0.01,0.2),ρ1=0.5×ρ0=-0.2,τ=0.99時,ARL1=148.1,如果ρ1=ρ0=-0.4則ARL1=136.6,如果ρ1=2×ρ0=-0.8則ARL1=117.5;增強正相關性會降低控制圖對比例偏移的敏感性,例如(γX,γY)=(0.01,0.2),ρ1=0.5×ρ0=0.2,τ=0.99時,ARL1=123.3,如果ρ1=ρ0=0.4則ARL1=134.1,如果ρ1=2×ρ0=-0.8則ARL1=160.5。

在圖3~圖5中,通過對比可知,對于單邊VSS Shewhart-RZ*控制圖,第1次抽樣采用大樣本容量(n(1)=nL)的性能優于第1次抽樣采用小樣本容量(n(1)=nS),在(γX,γY)=(0.01,0.01)、偏移較小時性能優勢尤為顯著。n(1)=nL時,雖然單邊VSS Shewhart-RZ*控制圖擁有更小的ARL1,但是同時也擁有更大的ASS1,例如(γX,γY)=(0.01,0.01),n0=5,ρ0=ρ1=-0.8,τ=0.99,如果n(1)=nS則ARL1=4.7,ASS1=8.7,如果n(1)=nL則ARL1=1.7,ASS1=28.9。

在圖3~圖5中,通過對比可知,單邊VSS Shewhart-RZ*控制圖的性能優于單邊FSS Shewhart-RZ控制圖,說明變樣本容量策略起到了很好的優化效果。當變量X,Y的變異系數均較小時,VSS策略對小偏移的優化效果更為明顯;當變量X,Y的變異系數有一者較大時,VSS策略對大偏移的優化效果更為明顯。在一些場景下,FSS Shewhart-RZ控制圖與VSS Shewhart-RZ*控制圖都能取得很好的效果,使ARL1=1,但是VSS Shewhart-RZ*控制圖可以使ASS1

在圖3中通過對比可知,在以下兩種情況下,單邊VSS Shewhart-RZ*控制圖的性能優于單邊Synthetic-RZ控制圖:

(1)變量X,Y的變異系數較小,即(γX,γY)=(0.01,0.01)時,n(1)=nL的單邊VSS Shewhart-RZ*控制圖性能優于單邊Synthetic-RZ控制圖,且偏移較小時優勢更為顯著。例如(γX,γY)=(0.01,0.01),ρ0=ρ1=-0.8,n0=5,τ=1.01時,單邊Synthetic-RZ控制圖ARL1=5.2,n(1)=nL的單邊VSS Shewhart-RZ*控制圖ARL1=1.8。

(2)當變量X,Y的變異系數較大,即(γX,γY)=(0.2,0.2),比例偏移較大(τ>1.1或τ<0.9)且受控狀態下的樣本容量較小(n0=5)時,單邊VSS Shewhart-RZ*控制圖同樣擁有更好的性能。例如(γX,γY)=(0.2,0.2),ρ0=ρ1=0,n0=5,τ=1.1時,單邊Synthetic-RZ控制圖ARL1=14.1,n(1)=nL的單邊VSS Shewhart-RZ*控制圖ARL1=8.4。

在VSS Shewhart-RZ*控制圖與VSI Shewhart-RZ控制圖的性能對比中,以平均報警時間ATS=E(h)×ARL作為性能指標。其中E(h)為平均抽樣間隔,對于VSS Shewhart-RZ*控制圖,其抽樣間隔h為固定值,取h=1,則ATS=ARL。因為VSI Shewhart-RZ控制圖在設計時滿足受控狀態下平均抽樣區間E(h)0=1的約束,所以可以通過比較失控狀態下VSI Shewhart-RZ控制圖ATS1指標與VSS Shewhart-RZ*控制圖ARL1指標的數值,來對比兩者的性能。圖3展示了抽樣區間組合(hS,hL)=(0.1,1.1),τ∈{0.95,0.98,0.99,1.01,1.02,1.05}時,VSI Shewhart-RZ控制圖失控狀態下的ATS1[12]??芍谝韵聝煞N情況下,單邊VSS Shewhart-RZ*控制圖的性能優于VSI Shewhart-RZ控制圖:

(1)當變量X,Y的變異系數較小,且相關系數ρ≤0、偏移和受控狀態下的樣本容量較小時,例如 (γX,γY)=(0.01,0.01),n0=5 ,ρ0=ρ1=-0.4 ,τ=1.01時,單邊VSI Shewhart-RZ控制圖ATS1=5.9,單邊VSS Shewhart-RZ*控制圖ARL1=4.0。

(2)當變量X,Y的變異系數較大,且偏移較大時,例如(γX,γY)=(0.2,0.2)n0=5,ρ0=ρ1=0,τ=0.95時,單邊VSI Shewhart-RZ控制圖ATS1=64.7,單邊VSS Shewhart-RZ*控制圖ARL1=53.5。

4 案例說明

表1 案例樣本數據

針對τ=1.01的偏移,在ARL0=200,n0=5的條件下對采用的4種控制圖進行優化設計,得到最優參數。對于單邊Shewhart-RZ控制圖計算可得UCL=1.015 4;對于單邊Synthetic-RZ控制圖,參考文獻[9]計算可得合格品鏈長子圖的控制限為H=4,單邊Shewhart-RZ子圖的控制限為UCL=1.010 7;對于單邊VSI Shewhart-RZ控制圖,參考文獻[12]計算可得抽樣間隔組合(hS,hf)=(0.1,1.1),警戒限UWL=1.007 5,控制限UCL=1.015 4;對于單邊VSS Shewhart-RZ*控制圖(n(1)=nS),根據第2章計算可得樣本容量組合(nS,nL)=(3,20),警戒限UWL=1.159 0,控制限UCL=2.575 8。

根據樣本數據及控制圖參數繪制4種比例控制圖如圖6所示,其中Synthetic-RZ控制圖僅繪制了單邊Shewhart-RZ子圖。由圖6可知,在前10個受控樣本下,所有控制圖均未誤發警報;當生產過程失控后,根據4種控制圖各自的運行規則,單邊Shewhart-RZ控制圖在第14個樣本點發出警報,報警時間tSh=4;單邊Synthetic-RZ控制圖在第12個樣本點發出警報,報警時間tSyn=2;單邊VSI Shewhart-RZ控制圖在第14個樣本點發出警報,報警時間tVSI=2.4;單邊VSS Shewhart-RZ*控制圖在第11~14樣本點發出警報,報警時間tVSS=1。當過程處于失控狀態時,4種比例控制圖均在較短時間內發出了警報,說明了其有效性。其中VSS Shewhart-RZ*控制圖在第11個點處及時發出警報,報警時間最短,且連續4點報警,具有很好的穩定性,效果最好。實際應用中,生產者應在VSS Shewhart-RZ*控制圖發出警報后及時停產檢查,在消除引起比例偏移的異常因素后再恢復生產,以免產生大量不合格品,保證產品質量的穩定性。

VSS Shewhart-RZ*控制圖性能優秀,且易理解、易實施,可廣泛應用于工程實踐。首先其原理簡單易理解,即檢測統計量距離控制限較遠,過程失控的可能性較小時,采用小樣本容量抽樣,以減少抽樣成本;檢測統計量距離控制限較近,過程失控的可能性較大時,采用大樣本容量抽樣,以更好地了解當前過程狀態。其次,VSS策略在應用中僅需改動每次抽樣的產品數量,便于操作,可廣泛應用于人工抽樣檢測、自動化抽樣檢測等作業場景。

5 結束語

猜你喜歡
樣本容量警戒警報
基于北斗三號的人防警報控制系統及應用
采用無核密度儀檢測壓實度的樣本容量確定方法
假期終結者
步兵班前進——警戒(XV)
步兵班前進——警戒(ⅩⅣ)
步兵班前進——警戒(XII)
步兵班前進——警戒(Ⅶ)
是誰的責任?
拉響夏日警報定格無痕跡美肌
廣義高斯分布參數估值與樣本容量關系
91香蕉高清国产线观看免费-97夜夜澡人人爽人人喊a-99久久久无码国产精品9-国产亚洲日韩欧美综合