呂安強 李永倩 李 靜 張 旭
(1.華北電力大學電氣與電子工程學院 保定 071003 2.華北電力大學信息與網絡管理中心 保定 071003)
我國擁有3.2 萬公里的海岸線和分布在沿海的6 500 多個大小島嶼,海底電纜應用廣泛,進行海底電纜的日常維護尤其是故障的早期發(fā)現與維修具有重要意義。近年來,船錨等漁具造成海纜損傷的機械事故頻發(fā),占海纜機械故障的80%以上[1]。嚴重錨害事故可導致斷路、短路或接地故障,一般能夠及時發(fā)現;而輕度損傷不能立即顯現,極具隱蔽性,往往會隨著時間的推移演變成滲水漏電、接地等故障,造成嚴重后果[2-4]。因此,對海纜進行實時監(jiān)測并判斷錨害程度勢在必行。目前的電纜電氣故障監(jiān)測與分析方法較多[5-7],但對海纜錨害程度的監(jiān)測與分析未見報道。國內有學者提出基于布里淵散射原理的分布式光纖傳感技術監(jiān)測海纜所受外力和內部溫度變化的方法,可初步實現海纜的溫度和應變監(jiān)測[8-10],為利用光纖應變判斷海纜錨害程度奠定了技術基礎。
利用光纖應變判斷海纜錨害程度,需要分析錨害過程中光纖應變的變化規(guī)律,使用有限元法能夠模擬海纜錨害過程,獲得光纖應變與海纜機械損傷程度之間的關系。目前已有學者對海底光纜進行了有限元仿真的探索。林開泉等人使用有限元仿真方法對海底光纜結構與抗錨害水平的關系進行了研究,為其結構的優(yōu)化設計提供了參考[11];張雁龍等人通過有限元法仿真了海底光纜在錨害過程中的變形情況,但研究未涉及海底光纜損傷程度與光纖應變之間的聯系[12,13];陸瑩等人針對分布式光纖傳感器在線監(jiān)測高壓海底電纜運行狀況的方法,利用有限元軟件建立了交聯聚乙烯(Cross Linked PolyEthylene,XLPE)高壓海底電纜模型,仿真分析海底電纜受外力損壞時內部物理量的變化[14],證明了利用分布式光纖傳感技術監(jiān)測海纜故障的可行性。
本文利用有限元分析法建立了110kV 光電復合海纜的有限元模型,詳細敘述了影響建模仿真精度和計算時間的關鍵要素處理方法,對海纜錨害過程進行了仿真和分析,獲得了光纖和鎧裝的應力、應變數據;建立了光纖應變和鎧裝應力的時間函數;對錨害程度進行了分級,并給出了相應的光纖應變閾值,為利用分布式光纖傳感技術進行海纜錨害監(jiān)測、預測和錨害程度判斷提供了理論參考。
福建省某海峽內敷設有YJQ41 型110kV 單芯XLPE 光電復合海纜,擬利用分布式光纖應變傳感設備進行海纜應變測量。目前的分布式光纖傳感技術能夠測量的物理量有溫度、應變和振動,其中,ROTDR、OFDR 可實現分布式溫度測量,BOTDR、BOTDA、COTDR 可實現分布式溫度和應變測量;POTDR、干涉式光纖傳感器可實現分布式振動測量。能測量應變的三種設備價格都比較高,BOTDR是高性價比的選擇,其傳感通道不需閉環(huán),即使監(jiān)測過程中光纖出現斷點,斷點前的光路仍能正常測量。本文依托該海纜監(jiān)測項目,對海纜錨害程度進行研究。海纜結構如圖1 所示。
圖1 海纜剖面圖Fig.1 Cross section of submarine cable
海纜結構的真實建模對錨害仿真結果的正確性起決定作用,但海纜結構的復雜性往往導致仿真時間過長。為了兼顧以上兩點,應在保證力學結構真實性的前提下,盡量簡化海纜有限元模型,以縮短仿真計算時間。海纜結構中,導體屏蔽、絕緣屏蔽、阻水帶和瀝青防腐層的厚度均不足海纜直徑的1%,且非剛性材料,在碰撞過程中起到的保護作用小,故對其做省略處理。簡化前后的海纜各層尺寸對比見下表。
表 海纜幾何尺寸參數Tab. Geometry parameters of submarine cable
錨害過程中,鎧裝層和滌綸樹脂(Poly Ethylene Tereph-thalate,PET)填充條起主要保護作用,它們以一定的節(jié)距絞合在高比密度聚乙烯(High-Density Poly Ethylene,HDPE)護套外面,呈復雜的空間結構,不是規(guī)則的幾何體,建模復雜。本文提出掃略螺旋線生成體的思路[15],通過點、線、面、體的步驟生成絞合鋼絲、光單元和PET 填充條,組成鎧裝層和PET 填充層。船錨的質量和錨纜接觸面積對錨害程度有影響,又因為船錨質地堅硬,錨害時不會發(fā)生形變,為減少計算時間,將船錨簡化為梯形體。
有限元模型參數的設置決定了仿真的準確性和計算時間。有限元仿真中的積分運算占據了大量的計算時間,其中單點積分的計算效率是全積分的8倍。因為海纜內銅導體、XLPE 絕緣、鉛合金護套、HDPE 護套、黃銅帶、光單元、PET 填充條、鋼絲鎧裝層和船錨都屬于三維實體,所以本文對以上實體選用支持所有非線性特性和單點積分的8 節(jié)點三維實體單元SOLID164,以兼顧仿真準確性和計算速度。海纜的繩被層和外被層厚度小,本文選用SHELL163 薄殼單元,以減少由于使用單點積分而產生的沙漏能。
海纜結構中的銅導體、鉛合金護套、黃銅帶和鎧裝層由于存在材料非線性問題,使用 ANSYS/LS-DYNA 中提供的雙線性隨動強化(BKIN)材料模型,使用彈性斜率和塑性斜率兩個直線段模擬彈塑性材料的本構方程。海纜的其他部分使用計算速度快的線彈性材料模型,其應力-應變關系滿足胡克定律。此外,船錨使用剛體材料模型,其在接觸和受力過程中不會產生變形,可提高計算速度30%左右。
本文依據拓撲正確性、幾何保持、特性一致、單元性狀優(yōu)良等原則進行網格劃分[16]。通過控制模型的輪廓線和輪廓面調整網格數量,兼顧計算的準確性和經濟性,最終生成的有限元模型單元數約9萬個,是通過軟件自動生成單元數的11%,在保證準確仿真的前提下大大提高了計算效率,網格劃分效果如圖2 所示。
圖2 有限元網格劃分效果Fig.2 Effect of finite element meshing
船錨對海纜的主要損傷方式為撞擊和鉤掛,撞擊時,海纜受力由上至下,經過接觸、拖動和擠壓過程,受力時間較短;鉤掛時,海纜受力由下至上,同樣經過以上過程,但受力時間相對較長。海纜埋設在海床下2m,海床淤泥質地稀軟,所以兩種情況下海纜和船錨所受阻力近似相等,兩種損傷方式可以合并為一種進行仿真。本文研究的110kV XLPE光電復合海纜敷設海域最大水深為25m,根據杭州灣水域案例和淺海船錨使用情況[17,18],實驗船錨選取660kg 規(guī)格的霍爾錨。
接觸定義的準確性決定了仿真的精度。錨害過程中,海纜的外表面會與內層及船錨外表面接觸,受力增大到一定程度時會產生穿透,因此,本文根據對稱罰函數法在接觸面之間引入界面接觸力,使用單面接觸類型。又因為錨害過程中同時存在剛體與變形體及變形體之間的接觸,接觸方向判斷困難,所以采用計算機自動識別方式,這是在人工無法準確判斷接觸集合時最有效的方法。
錨害發(fā)生位置附近一定長度的海纜會產生位移,而更遠距離處的海纜無位移。因此,對海纜兩端截面上的節(jié)點施加約束,使其位置固定,讓海纜的其他部分隨錨害發(fā)生變形與移動,模擬錨害的實際情況。另外,海纜敷設于海床下2m 的淤泥中,自身重力和洋流的影響平均分布于海纜全長上,由于錨害僅作用于海纜局部,海纜的應力、應變變化也僅產生于局部,所以仿真時忽略自重和洋流的影響。
利用建立的模型進行錨害仿真,得到海纜在錨害后的位移云圖如圖3 所示。
圖3 接觸過程中的位移變化Fig.3 Displacement change in contact process
由圖3 可見,在整個錨害過程中,海纜被船錨拖動并擠壓,船錨接觸位置海纜位移最大,海纜位移量與施加的船錨載荷位移量吻合。另外,沙漏能是衡量仿真正確性的重要指標,準確的仿真應保證沙漏能不超過內能的10%[19]。為了減少計算時間,本文采用了單點積分,將仿真時間縮短了80%,但同時也導致了沙漏能的產生。本文仿真結果的內能與沙漏能變化如圖4 所示。由圖4 可見,沙漏能被控制在內能的1.6%以內,保證了計算的準確性。
圖4 內能/沙漏能-時間曲線Fig.4 Internal energy/hourglass energy-time profile
海纜直徑大,兩根光單元在海纜內部成螺旋絞合狀,海纜軸向不同位置處,光單元處于海纜橫截面圓周的不同位置上,因此,錨害時船錨與海纜的接觸點與光單元在海纜橫截面圓周上的相對位置是隨機的。為了保證仿真的全面性,需考慮船錨從不同角度接觸海纜的情況。本文以30°為步進,進行了12 組仿真。
鎧裝層在海纜保護中起主要作用,鎧裝層的失效意味著其保護海纜內部結構功能的喪失,因此本文通過鎧裝層損傷判斷海纜的機械損傷。目前有4種強度理論可描述材料的失效情況,其中,形狀改變比能理論以Von Mises 應力為依據,最適合于鋼、銅、鋁等塑性金屬材料[20]。本文利用鎧裝層應力判斷鎧裝層損傷情況,進而判斷錨害程度。根據計算結果分別繪制12 組不同方向接觸過程中的鎧裝層應力-時間曲線和光纖應變-時間曲線,如圖5 所示。為清楚地展示鎧裝層應力與光纖應變的變化,使用K-means 算法[21]得到曲線的聚類中心,如圖5 中星狀曲線所示。該聚類中心與各曲線的誤差平方和最小,能夠反映各方向接觸時應力和應變隨時間變化的整體趨勢。
圖5 鎧裝層應力和光纖應變隨時間的變化Fig.5 Variation of armored layer stress and optical fiber strain over time
由圖5 可知,錨害過程中鎧裝層應力和光纖應變在不同接觸角度時幅度有差異,但整體變化趨勢基本一致,說明海纜錨害過程中以上兩個參量的變化是有規(guī)律的,聚類中心曲線可以描述錨害過程。鎧裝層應力經過了上升、保持、上升三個階段;第一階段,錨害剛剛開始,鎧裝層應力逐漸增加,鎧裝鋼絲處于彈性范圍內,應力隨時間近似呈線性變化;第二階段,鎧裝層應力超過了屈服應力,鎧裝產生了塑性應變,海纜內層結構產生位移,鎧裝鋼絲變得相對松散,應力下降;第三階段,錨害的加劇使鎧裝鋼絲再次拉緊受力,應力再次線性增加。仿真數據與材料的彈塑性理論相符。光纖是彈性材料,光纖應變隨錨害程度的加劇增加。相對于鎧裝的變化過程,第一階段,光纖余長導致應變基本不變;第二階段,光單元受力變長消耗了光纖余長,光纖應變開始增加,此后一直上升??梢岳霉饫w應變變化描述海纜錨害過程。
為了能通過光纖應變反映海纜的錨害過程,為基于光纖傳感的海纜狀態(tài)監(jiān)測提供理論支持,本文分別將光纖應變、鎧裝層應力和鎧裝層塑性應變的K-means 聚類結果按照各自最大值做歸一化處理,發(fā)現在鎧裝層應力變化過程中有四個關鍵轉折點,分別記為T1、T2、T3和T4,如圖6 所示。
圖6 歸一化的光纖應變/鎧裝層應力/鎧裝層塑性應變-時間曲線及擬合結果Fig.6 Normalized optical fiber strain/ armored layer stress/armored layer plastic strain-time profile and fitting results
0~T1過程中,船錨接觸并拖動海纜,鎧裝層承受的應力處于彈性范圍之內成線性增加,還沒有發(fā)生塑性應變;此過程中,鎧裝層對光單元具有一定的保護作用,加之光纖余長的存在,導致光纖應變增加很少;在T1時刻,鎧裝層應力達到屈服應力,鎧裝層開始產生不可恢復的塑性應變。T1~T2過程中,鎧裝層應力和塑性應變隨著錨害程度加劇而增加,達到峰值后進入保持狀態(tài);鎧裝層擠壓光單元,光纖應變在光纖余長消耗完后大幅上升;在 T2時刻,鎧裝層變形后進入強化階段。T2~T3過程中,船錨擠壓導致海纜橫截面變扁,鎧裝層結構變松散,鎧裝鋼絲長度相對增加,鎧裝層應力成下降趨勢;鎧裝層塑性應變保持不變;光纖應變增加速度變緩;在T3時刻鎧裝層應力下降至最低值。T3~T4過程中,錨害程度加劇使進入強化階段的鎧裝層再次產生應變,應力增加,但強化導致上升趨勢較0~T1時段變緩;塑性應變仍保持不變;光纖應變上升速度比上一階段有所增加;在T4時刻,鎧裝層應力再次增加到新的屈服應力數值,產生第二次塑性應變,即發(fā)生更嚴重的變形。T4時刻之后,鎧裝層應力和塑性應變、光纖應變都快速上升,達到頂點后鎧裝層鋼絲進入材料強度極限,發(fā)生斷裂,光纖也會斷裂。分別將不同時段的鎧裝層應力和光纖應變對時間步擬合得出擬合結果為
式中 s(t)——鎧裝層應力;
t——時間;
RMSE——鎧裝層應力的擬合均方根誤差。
式中 f(t)——光纖應變;
RMSE——光纖應變的擬合均方根誤差。
實際運行的海纜敷設于海床上或海床下的淤泥中,受洋流、潮汐、巖石摩擦、海床坡度變化等多種因素影響,導致鎧裝層沿海纜軸向的應力分布是變化的,但這些因素的影響相對于錨害而言是大時間尺度的,即在錨害發(fā)生的時間內,可以認為以上因素的影響是固定不變的,相當于海纜上某一點處的鎧裝應力在錨害前有一個初始值,初始值可能是變化的,但錨害過程中鎧裝應力和光纖應變的變化趨勢是有規(guī)律的。因此,公式中的常數項可根據實際情況修正。
受多種因素影響,錨害對海纜的損壞程度各有不同。為了能夠利用光纖應變反映海纜損壞程度,合理評估海纜狀態(tài),有必要給出表征海纜損壞程度的光纖應變閾值。圖7 所示是利用海纜截面變形情況直觀展示了海纜的錨害過程。
圖7 海纜截面變化Fig.7 Cross-section changes of submarine cable
圖7 中,T1時刻,海纜鎧裝層開始屈服并產生塑性應變,但海纜的整體結構基本不變,HDPE 護套以內形狀基本保持,海纜仍能正常運行,此時光纖應變處于0.12%~0.14%范圍,可以將光纖應變達到0.12%作為判斷海纜輕度損傷的閾值。T2時刻,鎧裝層發(fā)生了較大變形,內部結構在其保護下仍比較完整,但HDPE 護套以內形狀發(fā)生了變化,運行中可能發(fā)生局部放電,長期運行有可能導致漏電流的產生;此時,光纖應變增大至0.35%~0.42%范圍,可以將光纖應變達到0.35%作為判斷海纜中度損傷的閾值。T4時刻,海纜鎧裝和內部結構都產生了嚴重擠壓變形,鎧裝層的保護作用接近極限,海纜不能正常運行;此時光纖應變約為1.1%,可將此值作為判斷海纜重度損傷的閾值。此后,隨著錨害程度的加劇,海纜將產生更嚴重的扭曲變形,甚至有斷裂的危險。
海纜體積和重量大,進行實體錨害試驗困難,目前國內海纜廠家和研究院所都不具備試驗條件。由于本文研究的是海纜內復合光纖的應變隨海纜錨害的變化情況,進而反映海纜的錨害程度,所以可用海纜模型試驗驗證光纖應變仿真結果的正確性。
本文用塑膠圓管代替海纜光單元以內的海纜本體,根據實際海纜的直徑,按比例將裸光纖以一定節(jié)距絞合纏繞于塑膠圓管外壁,沿光纖軸向用膠帶固定光纖于圓管外壁,避免錨害時光纖產生不符合實際的軸向移位,同時起到保護光纖的作用,實物照片如圖8 所示。將海纜模型兩端固定,從中間施加外力模擬錨害故障,用BOTDR 監(jiān)測光纖上的應變變化,測量的應變曲線如圖 9 所示。圖 9a是BOTDR 測量的光纖應變分布曲線。由圖可見,由于往塑膠圓管上纏繞光纖時,對光纖施加了一定的預應力,所以光纖的初始應變不為零。隨著船錨拖動海纜位移的增加,海纜模型上光纖的應變成上升趨勢。將海纜模型上光纖的應變在長度上取平均,畫出它隨錨害位移增加的變化曲線如圖9b 所示。由圖可見,隨著位移的增加,光纖應變先增加,而后增幅變緩,之后增幅又上升,此變化趨勢與仿真結果的T1~T2、T2~T3和T3~T4過程一致,證明了仿真的正確性。
圖8 海纜模型照片Fig.8 Picture of submarine cable model
圖9 不同錨害位移下的光纖應變變化Fig.9 Optical fiber strain variation at different anchor damage displacement
利用有限元分析法建立海纜模型,可為分析錨害過程提供方便、有效的手段,并可獲得比實體實驗更豐富的數據。通過分析仿真數據,得出以下結論:
(1)海纜在錨害過程中,鎧裝層應力呈先增加、后減少、再增加的分段變化趨勢,在超過材料的屈服應力后會發(fā)生強化現象,產生不可逆的塑性應變,之后應力以相對較低斜率繼續(xù)增加。相同錨害過程中,鎧裝層應力和光纖應變在不同錨害角度下的數值存在差異,但它們具有相同的變化趨勢,可用各變化曲線的聚類中心描述這一趨勢。
(2)海纜鎧裝層應力和光纖應變的變化趨勢在多次錨害仿真中具有可重復性,可利用分段函數從理論上表示海纜錨害的過程,為利用分布式光纖傳感技術進行海纜錨害監(jiān)測、預測和錨害程度判斷提供理論參考。
(3)結合鎧裝層應力、塑性應變和海纜整體變形情況,可將海纜錨害分為輕度損傷、中度損傷和重度損傷三個等級,對應的光纖應變閾值為0.12%、0.35%和1.1%,作為判斷海纜錨害程度的依據。
[1]Kordahi M E,Shapiro S,Lucas G.Trends in submarine cable system faults[C].Submarine Optical Conference,2007,37.
[2]王新超,蘇秀蘋,張麗麗,等.短暫擾動注入法對于小電流接地系統線路狀態(tài)在線監(jiān)測的探討[J].電力系統保護與控制,2010,38(6):141-145.Wang Xinchao,Su Xiuping,Zhang Lili,et al.Discussion on the online monitoring for non-solidly grounded system power line based on a temporary disturbance injection[J].Power System Protection and Control,2010,38(6):141-145.
[3]周文,毛志強,毛志芳,等.一起線路接地故障引起機組主變保護誤動的事故分析[J].電力系統保護與控制,2010,38(12):149-152.Zhou Wen,Mao Zhiqiang,Mao Zhifang,et al.Analysis of the main transformer differential protection maloperation of generator units caused by single-phase ground fault[J].Power System Protection and Control,2010,38(12):149-152.
[4]朱韜析,彭武.天廣直流輸電系統線路高阻接地故障研究[J].電力系統保護與控制,2009,37(23):137-140.Zhu Taoxi,Peng Wu.Research on high impedance earth fault of Tian-Guang HVDC transmission project[J].Power System Protection and Control,2009,37(23):137-140.
[5]王雅群,尹毅,李旭光,等.等溫松弛電流用于10kV XLPE 電纜壽命評估的方法[J].電工技術學報,2009,24(9):33-37.Wang Yaqun,Yin Yi,Li Xuguang,et al.The method of lifetime evaluation on 10kV XLPE cables by isothermal relaxation current[J].Transactions of China Electrotechnical Society,2009,24(9):33-37.
[6]束洪春,彭仕欣.基于短窗數據 S 變換能量的纜-線混合配電網絡故障選線方法[J].電工技術學報,2009,24(10):152-159.Shu Hongchun,Peng Shixin.A fault line detection algorithm for distribution network of overhead line and underground cable mixed lines using s-transform energy from short window data[J].Transactions of China Electrotechnical Society,2009,24(10):152-159.
[7]吳曉文,舒乃秋,李洪濤,等.氣體絕緣輸電線路溫升數值計算及相關因素分析[J].電工技術學報,2013,28(1):65-72.Wu Xiaowen,Shu Naiqiu,Li Hongtao,et al.Temperature rise numerical calculation and correlative factors analysis of gas-insulated transmission lines[J].Transactions of China Electrotechnical Society,2013,28(1):65-72.
[8]蔣奇,徐于超,康彥森,等.基于分布式布里淵光纖散射傳感的海底動力電纜監(jiān)測技術研究[J].化工自動化及儀表,2009,36(4):41-43.Jiang Qi,Xu Yuchao,Kang Yansen,et al.Technological study on distributed brillouin fiber sensor monitoring of high voltage power cable in seabed[J].Control and Instruments In Chemical Industry,2009,36(4):41-43.
[9]楊黎鵬.基于光纖布里淵散射的分布式傳感海底電纜在線監(jiān)測技術研究[J].船海工程,2009,38(3):139-142.Yang Lipeng.Study on-line monitor technique based on brillouin scatting distributed optical fiber sensing[J].Ship &Ocean Engineering,2009,38(3):139-142.
[10]周學軍,王紅霞,趙四新.BOTDR 在海底光纜監(jiān)測中的應用[J].光纖與電纜及其應用技術,2006,40(4):42-44.Zhou Xuejun,Wang Hongxia,Zhao Sixin.BOTDR and its application in supervisionof submarine fiber optic cables[J].Optical Fiber &Electric Cable and Their Applications,2006,40(4):42-44.
[11]林開泉,王紅霞,劉紅亮,等.海底光纜錨害的有限元分析[J].電線電纜,2010,39(6):38-41.Lin Kaiquan,Wang Hongxia,Liu Hongliang,et al.Finite element analysis of the anchor-caused fault on the submarine optical fiber cables[J].Wire &Cable,2010,39(6):38-41.
[12]張雁龍,王紅霞,薛馳.海底光纜受壓過程的數值模擬[C].中國通信學會2009 年光纜電纜學術年會論文集,2009.
[13]張雁龍,王紅霞,王瑛劍.海底光纜錨害成因的實驗研究[J].光纖與電纜及其應用技術,2009,43(2):41-43.Zhang Yanlong,Wang Hongxia,Wang Yingjian.Experimental research on submarine fiber optic cable faults caused by anchor[J].Optical Fiber &Electric Cable and Their Applications,2009,43(2):41-43.
[14]陸瑩,黃輝.基于分布式光纖傳感技術的高壓海底電纜外力損壞仿真[J].電氣技術,2012,13(12):87-90.Lu Ying,Huang Hui.External force damage simulation of submarine high-voltage cable based on distributed fiber optic sensor technology[J].Electrical Engineering,2012,13(12):87-90.
[15]張旭,尹成群,呂安強,等.光電復合海纜有限元建模與仿真[J].計算機仿真,2013,30(2):80-85.Zhang Xu,Yin Chengqun,Lü Anqiang,et al.Method of constructing finite element model and simulating of photoelectric composite submarine cable[J].Computer Simulation,2013,30(2):80-85.
[16]胡于進,王璋奇.有限元分析及應用[M].北京:清華大學出版社,2009.
[17]顧玉娟.海底電纜的安全性探討[C].救撈專業(yè)委員會 2003 年學術交流會論文集,2003:126-127.
[18]劉小虎,易煒.船錨損壞海底敷埋電纜的風險評估方法[C].第七屆全國工程結構安全防護學術會議論文集,2009:268-272.
[19]趙海鷗.LS-DYNA 動力分析指南[M].北京:兵器工業(yè)出版社,2003.
[20]白鵬偉,史青錄,程結結,等.推土機傳動軸強度有限元分析[J].工程機械,2012,43(11):36-39.Bai Pengwei,Shi Qinglu,Cheng Jiejie,et al.Strength finite element analysis of bulldozer drive shaft[J].Construction Machinery and Equipment,2012,43(11):36-39.
[21]孫謙,姚建剛,李欣然,等.基于聚類趨勢分析與逐步回歸的電鐵牽引負載負序源模型研究[J].電機工程學報,2012,32(34):120-128.Sun Qian,Yao Jiangang,Li Xinran,et al.Study on negative sequence source model of electrified railway traction load based on clustering tendency analysis and stepwise regression[J].Proceedings of the CSEE,2012,32(34):120-128.