999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

衛星分群的抗差Kalman濾波在GPS/BDS融合精密單點定位中的應用

2016-04-15 03:19:22張傳定胡小工宋葉志馬紹龍唐成盼
中國慣性技術學報 2016年6期
關鍵詞:融合

王 琰,張傳定,胡小工,宋葉志,馬紹龍,唐成盼,李 冉

(1. 解放軍信息工程大學 地理空間信息學院,鄭州 450052;2. 北京衛星導航中心,北京 100094;3. 中科院上海天文臺,上海 200030;4. 南京陸軍指揮學院 軍隊指揮系,南京 210045)

衛星分群的抗差Kalman濾波在GPS/BDS融合精密單點定位中的應用

王 琰1,2,張傳定2,胡小工3,宋葉志3,馬紹龍4,唐成盼3,李 冉3

(1. 解放軍信息工程大學 地理空間信息學院,鄭州 450052;2. 北京衛星導航中心,北京 100094;3. 中科院上海天文臺,上海 200030;4. 南京陸軍指揮學院 軍隊指揮系,南京 210045)

抗差Kalman濾波是控制GNSS動態導航定位中觀測異常的有效算法,當應用到GPS/BDS實時動態精密單點定位(Precise Point Positioning, PPP)時,會出現某些歷元定位精度甚至不如單一系統定位精度高,這主要是因為同一接收機接收的不同種類衛星觀測量的隨機特性不同,使得觀測量驗后殘差的分布特性不一致,抗差估計時隨機特性不同的觀測量驗后殘差互比,反而對某一系統優質數據也進行了降權,導致定位結果出現偏差,減弱了 GPS/BDS融合精密單點定位的優勢。針對這一問題,提出了衛星分群的抗差Kalman濾波算法,并應用到GPS/BDS融合精密單點定位中,算法的核心是在每一歷元觀測數據質量控制時根據衛星類型分類構建方差膨脹因子,給出了算法的實施步驟,最后通過MGEX實測數據進行了驗證,結果表明算法應用到GPS/BDS融合精密單點定位中,相較傳統的抗差Kalman濾波算法在東、北、天三個方向分別提高了34.6%、33.3%、31.0%,同時表明該算法提高了GPS/BDS融合精密單點定位的可靠性。

抗差估計;Kalman濾波;衛星分類定權;GPS/BDS;精密單點定位

精密單點定位(Precise Point Positioning,PPP)技術利用衛星精密星歷和鐘差產品獲得全球任一接收機在ITRF下的絕對坐標,該技術經過了15年的發展,應用前景廣闊[1]。北斗衛星導航系統(BeiDou Navigation Satellite System,BDS)是我國獨立自主發展的導航系統,目前已經為亞太地區提供導航定位服務[2],GPS/BDS融合PPP吸引了一系列專家學者的廣泛關注,相較單一導航系統的PPP解,GPS/BDS融合PPP增加了觀測冗余,改善了定位的DOP值,能夠提高PPP的精度和可靠性,縮短初始化時間,其優勢毋庸置疑[2-7]。文獻[2~7]研究了多系統融合PPP的函數模型,并采用實測數據評估了多系統融合PPP的定位精度與收斂時間,闡述了多系統融合 PPP的優勢;文獻[4][5]對GPS/BDS PPP的定位精度與收斂時間進行了比較。PPP技術還有許多關鍵算法需要改進,文獻[8]對影響PPP定位的關鍵因素進行了分析研究。

眾所周知,GNSS實時數據質量控制是保證動態導航定位精度的關鍵,對于GPS/BDS融合動態PPP,要盡量利用可以利用的觀測數據,但要保證參與PPP解算的觀測數據的質量,若質量較差的觀測數據參與平差,非但無法起到增加觀測冗余的作用,反而影響定位精度[9]。實時PPP隨著歷元向前處理,過程不可逆,因此無法像事后PPP根據前后歷元觀測量時間序列特性進行質量控制,故實時PPP質量控制更為重要[10]。

Kalman濾波在動態數據處理中應用較廣,實時PPP一般也采用Kalman濾波,但是標準Kalman濾波無法對數據中的粗差進行有效處理。為了控制Kalman濾波遞推過程中的觀測異常,有兩種思路:一種是文獻[11]提出的與Kalman濾波并行操作的誤差探測、診斷與修復的DIA(Detection, Identification, Adaptation)方法,該算法的理論基礎是將粗差歸于均值漂移的粗差探測技術;另一種是將抗差估計與 Kalman濾波相結合的抗差Kalman濾波,抗差Kalman濾波的理論基礎是基于將粗差歸于方差膨脹的穩健估計技術,不同于DIA算法對“值得懷疑”的異常數據進行硬性拒絕,抗差Kalman濾波通過對觀測數據“降權”、“保權”、“拒絕”三種方式盡量利用可以利用的觀測數據,從而保證定位結果的精度和可靠性,因此在GNSS動態數據處理中應用較廣,其中以楊元喜等提出的抗差自適應濾波理論最具代表,該算法能能夠更有效地控制異常影響,提高動態濾波精度,相關研究很多不再詳述[12-13]。

抗差Kalman濾波在PPP中廣泛應用,文獻[9][10]對實時PPP的質量控制進行了詳細研究,文獻[13]提出了一種改進的抗差Kalman濾波方法并應用到精密單點定位中。對抗差估計的原理進行研究發現,抗差估計實際是根據參與平差的觀測量驗后殘差進行互比,對驗后殘差離群的觀測量進行降權,這個前提是參與定位的觀測量驗后殘差向量同類、同分布。GPS/BDS融合PPP與基線解的模式不同,衛星星歷誤差無法通過差分消除或減弱,IGS提供的GPS衛星精密星歷和鐘差精度較高,而BDS采用了混合星座,不同種類的衛星星歷精度存在偏差,同時軌道高度不同,外加測量噪聲本身隨機特性的不同,都造成GPS/BDS融合PPP觀測量驗后殘差很難做到同類、同分布,因此本文提出了一種衛星分群的抗差Kalman濾波算法,并應用到GPS/BDS融合PPP中,介紹了算法的實施流程,最后采用MGEX實測數據驗證了算法的有效性。

1 多系統融合PPP

1.1 GPS/BDS融合PPP的觀測方程

GPS/BDS融合PPP示意圖如圖 1。PPP一般采用消電離層組合偽距和載波相位觀測值,消去電離層一階項誤差。衛星軌道和鐘差固定(一般采用IGS分析中心提供的精密星歷和鐘差產品),衛星鐘差產品中包含了衛星端的偽距硬件延遲,接收機端的偽距硬件延遲被接收機鐘差吸收。GPS/BDS融合PPP中GPS的觀測方程如式(1)[1,3-4]:

式中:偽距和相位硬件延遲與模糊度參數無法分離,這也是PPP中模糊度無法固定的原因。PPP中模糊度參數一般采用浮點解:

前述表明接收機鐘差參數會吸收接收機端的偽距硬件延遲,而此延遲與信號頻率和導航系統相關,因此兩個系統在接收機端會產生偽距硬件延遲之差,也即碼偏差(Differenced Code Bias,DCB),另外由于不同導航系統時間基準有差異,因而在 GPS/BDS融合PPP時,同一個接收機對不同的導航系統觀測方程中會采用不同的接收機鐘差參數。融合PPP的觀測方程一般以GPS系統的接收機鐘差作為基準,其他系統的觀測方程中增加ISB參數,該參數吸收了DCB和時間基準的系統偏差[3-5]。因此 GPS/BDS融合 PPP中BDS的觀測方程為

圖1 GPS/BDS融合PPP示意圖Fig.1 Schematic of combining GPS/BDS with PPP

1.2 擴展Kalman濾波(EKF)

GPS/BDS融合PPP相較單一導航系統PPP具有很大的優勢,但是作者在實際處理時發現,當某一系統的觀測數據出現異常,若處理不當可能導致融合PPP的定位精度甚至不如單一導航系統的定位精度高。下面采用實測數據的算例進行說明。

1.3 試驗分析

為了避免 EKF中由于載體運動導致狀態參數出現異常的情況,采用靜態仿動態的數據處理模式,以MGEX(Multi-GNSS Experiment)兩個靜態監測站GMSD和JFNG站2014年年積日070天的數據進行說明,GMSD和JFNG站分別位于日本和中國,兩站接收GPS/BDS衛星情況如圖 2所示,采用GFZ和WHU提供的星歷產品進行PPP解算[14-15],比較以下三種方案的定位結果:方案1:單GPS PPP,采用IGS精密星歷和30 s采樣率的精密鐘差產品;方案2:GPS/BDS融合 PPP,采用 GFZ提供的星歷產品;方案 3:GPS/BDS融合PPP,采用WHU提供的星歷產品。

圖2 可觀測衛星數Fig.2 Number of visible satellites

三種方案解算策略一致,僅僅是采用的星歷產品不同,對模糊度收斂后的定位結果進行統計,圖 3、圖 4分別為GMSD與JFNG站為三種方案每天50 min之后的PPP定位結果時間序列。

從圖3和圖4的結果可以看出:

1)GMSD站的定位結果說明GPS/BDS融合PPP能夠提高定位精度,方案2、方案3定位結果的時間序列比方案1更加穩定,不會出現方案1中紅色點跡的“毛刺”現象。同時圖4中JFNG站方案1的定位序列在19:00-21:00天方向的定位結果較差,而同時期方案3的定位結果無異常,這是因為該段時間GPS觀測數據存在異常值,BDS觀測數據的加入彌補了這一不足,充分說明了GPS/BDS融合PPP的優勢。

圖3 GMSD站三種方案定位結果Fig.3 Positioning results of three schemes (site GMSD)

2)圖 4中JFNG站方案1的定位序列在 12:00前后無明顯異常,說明GPS觀測量無明顯異常,但是方案3的定位結果序列在這個時間段明顯異常,分析是 BDS的觀測數據存在異常造成了 GPS/BDS融合PPP定位結果在該時間段變差。

圖4 JFNG站三種方案定位結果Fig.4 Positioning results of three schemes (site JFNG)

表1 GMSD與JFNG站三種方案的定位結果(RMS)Tab.1 Positioning results (RMS) of three schemes (site GMSD and JFNG)

3)圖 4中JFNG站方案2在15:00-21:00時段內東方向的定位精度比方案3的要差,分析是由于該段時間GFZ提供的BDS軌道和鐘差的精度不高,PPP數據處理模式與基線解不同,衛星端的星歷誤差無法通過差分消除,導致GPS/BDS融合PPP的定位精度甚至不如單 GPS系統的好。對于這一問題可以采用Helmert方差分量估計解決,擬在其他文章中另行介紹,本文不做討論。

本文主要對前2個問題進行討論。當兩大系統每個歷元的觀測數據都存在冗余且均無異常觀測量時,BDS/GPS融合PPP的定位精度明顯要優于單GPS的定位精度。但是當其中一個系統的觀測數據出現異常(星歷精度的問題或者觀測量本身的問題),導致另一系統的觀測數據殘差整體變大,使定位結果出現偏差,甚至不如單一導航系統的定位精度。

根據Helmert方差分量估計的思想,利用觀測量的驗后殘差計算不同種類觀測值的方差因子,從而對哪一類觀測量出現了異常進行定位。式(6)為近似Helmert方差分量估計的Forstner公式,其中為第i類觀測量的方差分量估計值,ni為當前歷元該類觀測量個數。按照文獻[12][13]的研究,歷史的觀測異常誤差會對Kalman濾波結果造成影響,從而影響對當前歷元方差分量估計值的計算,故本文選取當前歷元的觀測量驗后殘差構建方差因子。若某一類觀測在某一歷元的方差分量估計值變大,說明當前歷元此類觀測值中存在異常。將GPS/BDS融合PPP的衛星按照類型分為四類:GPS、BDS GEO、BDS IGSO、BDS MEO。圖5為JFNG站方案3中四類衛星的時間序列。

圖5 JFNG站第i類觀測量方差分量估計值的時間序列Fig.5 Variance-component estimation valueof type-i observables in site JFNG

圖5中10:00至14:00之間GPS、GEO、IGSO三類衛星觀測量的值均變大,前邊分析了是由于BDS的觀測量存在異常,GPS觀測量的驗后殘差也相應變大,分別對JFNG站GPS/BDS的觀測量驗后殘差做分析。圖 6為JFNG站10:00-15:00之間GPS/BDS觀測量驗后殘差序列:上圖為GPS的結果,在12:00與14:00前后存在觀測量殘差出現離群的現象,但是不明顯;下圖為BDS的結果,明顯地在11~12 h期間C01/C04星觀測量出現異常,基于前后歷元時間序列進行粗差探測,判定是 C01星觀測量出現了異常,導致圖 5中12:00左右BDS GEO衛星的方差分量估計值出現異常,相同歷元的觀測量驗后殘差均出現離群狀態(GPS以及C04星的結果可以說明)。

圖6 JFNG站10:00-15:00之間GPS/BDS觀測量殘差Fig.6 Post-fit phase residuals of GPS/BDS (time 10:00-15:00, site JFNG)

綜上,對于GPS/BDS融合PPP,要避免觀測數據中上述的異常才能發揮融合PPP的優勢,抗差Kalman濾波是一種有效手段,但是若按照常規抗差 Kalman濾波計算等價權時,圖 5中12:00前后GPS觀測量驗后殘差也變大,C01觀測量在降權的同時GPS衛星的觀測量也整體降權,這就更造成了定位的偏差。因此本文提出了衛星分群的抗差 Kalman濾波算法來解決該問題,下面對算法進行介紹。

2 衛星分群抗差Kalman濾波

2.1 抗差Kalman濾波

抗差M-LS濾波極值條件[12]:

比較式(4)(8)兩者的差異僅僅是在測量噪聲的方差協方差矩陣由Qk變為了等價方差協方差矩陣,從而引起濾波增益矩陣Kk發生變化。鑒于 GPS/BDS融合PPP所有觀測量是不相關的,將Qk對角線上第i個元素Qki乘以方差膨脹因子,就置換為。方差膨脹因子函數如式(9),該式利用IGG III權因子函數的倒數構建[12-13]。

式中:K0、K1根據經驗選取,在本文計算中,K0=1.5,K1=3.0;vi是當前歷元某一相位觀測量的驗后殘差,通過當前歷元觀測殘差向量序列v=(v1,v2,…,vk)將當前vi化成標準化觀測殘差,的構建有多種方式,本文選取式(10)[13],其中σv是v的中誤差。

2.2 衛星分群方差膨脹因子

前述分析知,GPS/BDS融合PPP中不同種類衛星的觀測量驗后殘差序列v的隨機特性不同,將當前歷元觀測殘差向量序列v按照衛星類型分為四類v=(vG,vCG,vCI,vCM),G代表GPS衛星,CG代表BDS中GEO衛星,CI代表BDS中IGSO衛星,CM代表BDS中MEO衛星,不同類型衛星的觀測量驗后殘差計算等價權因子如式(11):

式中:σvG、σvCG、σvCI、σvCM為各類衛星觀測量驗后殘差的中誤差。每次抗差迭代時僅對該群觀測量中驗后殘差最大的觀測量進行方差膨脹,降低其對參數估計的貢獻[13]。

3 試驗分析

為了檢驗本文算法應用到GPS/BDS融合PPP中是否有效,收集了GMSD與JFNG站2014年年積日070-076共7天的數據進行試驗分析,設計了三種方案:

方案1:標準Kalman濾波,采用最簡單的3σ準則剔除每一歷元的粗差數據。

方案2:常規的抗差Kalman濾波。

方案 3:采用本文提出的衛星分群抗差 Kalman濾波。

其中,GPS/BDS衛星精密星歷和精密鐘差選用WHU提供的產品。三種方案除每個歷元數據質量控制模塊不同外,其余數據處理策略均相同,從定位精度與可靠性兩個方面評估本文算法的有效性。

3.1 定位精度

統計每天模糊度收斂后(每天50 min后)定位序列的RMS(結果如圖 7),三種方案7天定位序列RMS的平均值如表 2所示。

圖 7 三種方案定位序列的RMS(上圖為GMSD站,下圖為JFNG站)Fig.7 Positioning results RMS of three schemes (The top one is the result of site GMSD, and the bottom one is the result of site JFNG)

表 2 三種方案定位序列RMS的平均值Tab.2 Mean values of positioning results by the three schemes

從圖 7兩個站多天的定位結果可以看出,方案3的定位結果最優,方案2次之,方案1的結果最差。表 2的統計結果顯示本文算法相較傳統的抗差Kalman濾波算法平均在東方向提高了34.6%,北方向提高了33.3%,天方向提高了31.0%。

3.2 可靠性

統計三種方案出現 GPS/BDS融合解比 GPS或BDS單獨定位精度差情況的概率,作為衡量三種算法可靠性的一個重要指標??煽啃员容^的方法:對每天模糊度收斂后的定位序列進行統計(每天50 min,7天共19 460個歷元),統計所有歷元中三種方案定位結果比GPS或者BDS單獨定位精度差的百分比,若方案有效,則出現定位下降的百分比較小,結果如圖 8所示。圖 9為方案3這7天定位序列與GPS、BDS單獨定位的時間序列。

采用本文提出的衛星分群的抗差 Kalman濾波算法,與傳統的Kalman濾波法相比,GPS/BDS融合PPP定位結果的可靠性增強,出現融合解的結果比單獨解的結果差的概率變小,但是仍存在某些歷元出現GPS/BDS融合PPP的結果比單GPS、BDS的結果差的情況。

圖8 三種方案可靠性比較Fig.8 Reliability comparison among three schemes

圖9 GMSD站三種模式的定位序列Fig.9 Positioning sequence of the three schemes (site GMSD)

4 結 論

傳統的抗差Kalman濾波算法在 GPS/BDS融合PPP數據處理中應用會存在些許不足:實時質量控制時因為不同隨機特性的觀測量互比,可能導致優質觀測量過分降權的問題,導致定位結果甚至不如單一導航系統的定位精度高。本文采用MGEX實測數據對該問題進行了分析說明,并提出了改進的衛星分群抗差Kalman濾波算法。

算法的基本原理是根據衛星種類分群構建方差膨脹因子,避免同一接收機不同種類衛星的觀測量互比,造成其中某一系統的正常觀測量過分降權。本文給出了算法的實施步驟,最后采用MGEX實測數據檢驗了算法的有效性。與傳統的抗差Kalman濾波算法相比,本文算法對GPS/BDS融合PPP的定位精度有一定提高,并且提高了GPS/BDS融合PPP的可靠性。

(References):

[1] Zumberge J, Heflin M, Jefferson D, Watkins M, et al. Precise point positioning for the efficient and robust analysis of GPS data from large networks[J]. Journal of Geophysical Research, 1997, 102 (B3): 5005-5017.

[2] 楊元喜. 北斗衛星導航系統的進展、貢獻與挑戰[J]. 測繪學報, 2010, 39(1): 1-6. Yang Yuan-xi. Progress, contribution and challenge of compass/BeiDou satellite navigation system[J]. Acta Geodaetica et Cartographica Sinica, 2010, 39(1): 1-6.

[3] Cai C, Gao Y. Modeling and assessment of combined GPS/GLONASS precise point positioning[J]. GPS Solut, 2013, 17(2): 223-236.

[4] Zhang X, Li X. Instantaneous re-initialization in real-time kinematic PPP with cycle slip fixing[J]. GPS Solutions, 2012, 16(3): 315-327.

[5] Li X, Ge M, Dai X, et al. Accuracy and reliability of multi-GNSS real-time precise positioning: GPS, GLONASS, BeiDou, and Galileo[J]. J Geod, 2015, 89(6): 607-635.

[6] 米洋, 陳家斌, 劉紅光, 等. 一種采用小波神經網絡的GPS精密單點定位方法[J]. 中國慣性技術學報, 2016, 24(3): 337-341. Mi Yang, Chen Jia-bin, Liu Hong-guang, et al. GPS precise point positioning method using wavelet neural network[J]. Journal of Chinese Inertial Technology, 2016, 24(3): 337-341.

[7] 劉帥, 孫付平, 李海峰, 等. 前后向平滑算法在精密單點定位/INS緊組合數據后處理中的應用[J]. 中國慣性技術學報, 2015, 23(1): 85-91. Liu Shuai, Sun Fu-ping, Li Hai-feng, et al. Forwardbackward-smoothing algorithm with application to tightly coupled PPP/INS data post-processing[J]. Journal of Chinese Inertial Technology, 2015, 23(1): 85-91.

[8] Seepersad G, Bisnath S. Challenges in assessing PPP performance[J]. Journal of Applied Geodesy, 2014; 8(3): 205-222.

[9] 蔡華, 趙齊樂, 孫漢榮, 等. GNSS實時數據質量控制[J]. 武漢大學學報信息科學版, 2011, 36(7): 820-824. Cai Hua, Zhao Qi-le, Sun Han-rong, et al. GNSS real-time quality control[J]. Geomatics and Information Science of Wuhan University, 2015, 36(7): 820-824.

[10] 張小紅, 郭斐, 李盼, 等.GNSS精密單點定位中的實時質量控制[J]. 武漢大學學報(信息科學版), 2012, 37(8): 940-944. Zhang Xiao-hong, Guo Fei, Li Pan, et al. Real-time quality control procedure for GNSS precise point positioning[J]. Geomatics and Information Science of Wuhan University, 2012, 37(8): 940-944.

[11] Teunissen P J G. Quality control in navigation systems[J]. IEEE Aerospace and Electronic System Magazine, 1990, 5(7): 35-41.

[12] Yang Yuan-xi. Robust Bayesian estimation [J]. Bulletine Geodesique, 1991, 65: 145-150.

[13] 張小紅, 潘宇明, 左翔, 等. 一種改進的抗差 Kalman濾波方法在精密單點定位中的應用[J]. 武漢大學學報信息科學版, 2015, 40(7): 858-864. Zhang Xiao-hong, Pan Yu-ming, Zuo Xiang, et al. An improved robust Kalman filtering and its application in PPP[J]. Geomatics and Information Science of Wuhan University, 2015, 40(7): 858-864.

[14] Uhlemann M, Gendt G, Ramatschi M, et al. GFZ global multi-GNSS network and data processing results[C]//International Association of Geodesy Symposia. 2016: 673-679.

[15] Guo Jing, Xu Xiao-long, Zhao Qi-le, et al. Precise orbit determination for quad-constellation satellites at Wuhan University: strategy, result validation, and comparison[J]. J Geod, 2016, 90(2): 143-159.

Robust Kalman filtering based on different satellite types and it’s application in GPS/BDS precise point positioning

WANG Yan1,2, ZHANG Chuan-ding2, HU Xiao-gong3, SONG Ye-zhi3, MA Shao-long4, TANG Cheng-pan3, LI Ran3
(1. Institute of Surveying and Mapping, Information Engineering University, Zhengzhou 450052, China; 2. Beijing Satellite Navigation Center, Beijing 100094; 3. Shanghai Astronomical Observatory, Chinese Academy of Sciences, Shanghai 200030, China; 4. Nanjing Army Command College, Nanjing 210045, China)

A robust Kalman filtering is used to control the abnormal errors of observations in GNSS dynamic navigation and positioning. However, when it is used in combined GPS/BDS dynamic precise point positioning (PPP), the positioning accuracy may be worse than that of the single system PPP. Since the random characteristics of different types of satellite’s observables are different, the post-fit residuals of these observables are different either. When the post-fit residuals are used to calculate the equivalent weight factor in robust estimation, the high quality data may be down weighted, leading to the deviation of the positioning result. So the advantage of combining GPS/BDS with PPP is weakened. To solve this problem, the equivalent weight factor calculation based on different satellite types in robust Kalman filtering is proposed and the execution process of the algorithm is given. Finally, the experiment based on MGEX data is carried out. Compared with the robust Kalman filtering, the proposed algorithm improves the positioning accuracy and reliability of GPS/BDS with PPP. The positioning accuracy in ENU directions are increased by 34.6%, 33.3%, and 31.0%, respectively.

robust estimation; Kalman filtering; satellite classification weighting; GPS/BDS; precise point positioning

P227

:A

2016-08-25;

:2016-11-28

國家自然科學基金(41374038,41204022,41504018)

王琰(1990—),男,博士研究生,從事測量數據處理理論與方法研究。E-mail: wang1yan.hi@163.com

1005-6734(2016)06-0769-06

10.13695/j.cnki.12-1222/o3.2016.06.013

猜你喜歡
融合
一次函數“四融合”
兩個壓縮體融合為一個壓縮體的充分必要條件
村企黨建聯建融合共贏
今日農業(2021年19期)2022-01-12 06:16:36
融合菜
寬窄融合便攜箱TPFS500
寬窄融合便攜箱IPFS500
從創新出發,與高考數列相遇、融合
寬窄融合便攜箱IPFS500
《融合》
現代出版(2020年3期)2020-06-20 07:10:34
“四心融合”架起頤養“幸福橋”
福利中國(2015年4期)2015-01-03 08:03:38
主站蜘蛛池模板: 欧美午夜小视频| 国产香蕉国产精品偷在线观看| 日本人妻一区二区三区不卡影院 | 一本色道久久88| 国产精品护士| 在线看片中文字幕| 久久免费视频播放| 白浆视频在线观看| 国产99精品久久| 99热这里只有精品国产99| 一级爆乳无码av| 日韩欧美网址| 日本黄色a视频| 国产在线精品香蕉麻豆| 999福利激情视频| 亚洲AV无码一区二区三区牲色| 国产成人AV综合久久| 国产精品3p视频| 亚洲国产精品国自产拍A| 福利在线一区| 91香蕉视频下载网站| 国产成人亚洲精品色欲AV| 欧美日韩一区二区三区在线视频| 久久香蕉国产线看精品| 日韩免费毛片| 国产亚洲精久久久久久无码AV| 操美女免费网站| 亚洲精品爱草草视频在线| 国产精品夜夜嗨视频免费视频| 亚洲娇小与黑人巨大交| 亚洲日韩高清在线亚洲专区| 91久久青青草原精品国产| 亚洲日韩在线满18点击进入| 久久精品日日躁夜夜躁欧美| 久久精品嫩草研究院| 国产一级妓女av网站| 欧美国产菊爆免费观看| 色天堂无毒不卡| 欧美不卡视频在线观看| 91视频首页| 亚洲va在线∨a天堂va欧美va| 亚洲精品国产综合99| 五月丁香在线视频| 99精品高清在线播放| 国产91在线|中文| 在线看片中文字幕| 午夜精品影院| 日本五区在线不卡精品| 亚洲黄色网站视频| 伊人久久影视| 国产一区成人| 中文字幕有乳无码| 亚洲日韩AV无码精品| 久久国产精品影院| 婷婷六月综合网| 狠狠色香婷婷久久亚洲精品| 成人精品免费视频| 久久久久久久97| 国产精品自在线拍国产电影| 欧美激情一区二区三区成人| 18禁高潮出水呻吟娇喘蜜芽| 中文字幕免费在线视频| 国禁国产you女视频网站| 亚洲日韩精品欧美中文字幕| 国产另类视频| 黄色在线网| 精品国产99久久| 日韩精品一区二区三区大桥未久| 26uuu国产精品视频| 亚洲天堂精品视频| 国产尤物在线播放| 亚洲va欧美ⅴa国产va影院| 亚洲高清中文字幕| 日韩人妻无码制服丝袜视频| 国产人成乱码视频免费观看| 国产成人一区免费观看| 欧美精品不卡| 在线欧美日韩国产| 1级黄色毛片| 无码AV日韩一二三区| 国产精品无码久久久久久| 国产在线观看一区二区三区|